Alexander_KS/SAPFOR
4
2
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
Files
9d0fee77c73989977d65941478cb16a0311045a9
SAPFOR/tests/sapfor/parameter/magnit_3d.for

1514 lines
38 KiB
Plaintext
Raw Normal View History

added project
2023-09-14 19:43:13 +03:00
! ----------------------------------------------------------------------
! MAGNIT_3D. 14/06/2020
! ----------------------------------------------------------------------
IMPLICIT INTEGER*4 (I-N)
IMPLICIT REAL*8 (A-H,O-Z)
! ---
! include 'omp_lib.h'
!
! INTEGER(KIND=OMP_LOCK_KIND) :: LOCK
! ---
INTEGER*4 ,ALLOCATABLE:: MAP(:,:)
!
REAL*8 ,ALLOCATABLE:: RR(:),PP(:),ZZ(:)
REAL*8 ,ALLOCATABLE:: RS(:),PS(:),ZS(:)
!
REAL*8 ,ALLOCATABLE:: PHI1(:,:,:),PHI2(:,:,:),PHI3(:,:,:)
REAL*8 ,ALLOCATABLE:: ARX_RZ(:,:),ARY_RZ(:,:),ARZ_RZ(:,:)
REAL*8 ,ALLOCATABLE:: AR_RZ(:,:),AR_RP(:,:)
!
REAL*8 ,ALLOCATABLE:: HX(:,:,:),HY(:,:,:),HZ(:,:,:),HR(:,:,:),HP(:
&,:,:)
!
REAL*8 ,ALLOCATABLE:: BX(:,:,:),BY(:,:,:),BZ(:,:,:),BR(:,:,:),BP(:
&,:,:)
!
CHARACTER*8 :: SNUM
CHARACTER*256 :: VNAME
CHARACTER*1024 :: STR
! ---
COMMON /BMAIN/DR,DP,DZ,CFR,CFP,CFZ
COMMON /BMATR/VMUMAT(4,4),VJRM(4),VJZM(4)
COMMON /BITER/EPS,TAU,RKC,PHI0,ITM,ITP,ITS,K0
! ---
CPU_T1 = OMP_GET_WTIME ()
!
! CALL OMP_INIT_LOCK(LOCK)
! ---
WRITE (UNIT = *,FMT = *) 'command: magnit_3d <vname> [<nthreads> [
&<lprint>]]'
WRITE (UNIT = *,FMT = *)
! ---
VNAME = 'magnit_3d'
MTHR = OMP_GET_MAX_THREADS ()
NTHR = 1
LPP = 0
! ->
MXARG = IARGC ()
! ->
IF (MXARG .GT. 0) THEN
I = 1
STR = ' '
CALL GETARG(I,STR)
L = MIN0 (256,LENSTR (STR))
VNAME = STR(1:L)
ENDIF
!
LVNAME = LENSTR (VNAME)
! ->
IF (MXARG .GT. 1) THEN
I = 2
STR = ' '
CALL GETARG(I,STR)
READ (UNIT = STR,FMT = *) NTHR
ENDIF
! ->
IF (MXARG .GT. 2) THEN
I = 3
STR = ' '
CALL GETARG(I,STR)
READ (UNIT = STR,FMT = *) LPP
ENDIF
! ---
OPEN (UNIT = 2,FILE = VNAME(1:LVNAME)//'.log')
! ---
WRITE (UNIT = *,FMT = *) '---> Start'
WRITE (UNIT = 2,FMT = *) '---> Start'
! ---
WRITE (UNIT = *,FMT = *) 'MAGNIT_3D'
WRITE (UNIT = *,FMT = *) 'vname=',VNAME(1:LVNAME)
WRITE (UNIT = *,FMT = *) ' nthr=',NTHR
WRITE (UNIT = *,FMT = *) ' mthr=',MTHR
WRITE (UNIT = *,FMT = *) ' lpp=',LPP
!
WRITE (UNIT = 2,FMT = *) 'MAGNIT_3D'
WRITE (UNIT = 2,FMT = *) 'vname=',VNAME(1:LVNAME)
WRITE (UNIT = 2,FMT = *) ' nthr=',NTHR
WRITE (UNIT = 2,FMT = *) ' mthr=',MTHR
WRITE (UNIT = 2,FMT = *) ' lpp=',LPP
WRITE (UNIT = 2,FMT = *)
! ---
WRITE (UNIT = *,FMT = *) '---> Read'
WRITE (UNIT = 2,FMT = *) '---> Read'
! ---
OPEN (UNIT = 1,FILE = VNAME(1:LVNAME)//'.d')
! ->
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) R11
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) R12
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) Z11
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) Z12
! ->
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) R21
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) R22
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) Z21
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) Z22
! ->
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) R31
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) R32
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) Z31
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) Z32
! ->
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) R41
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) R42
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) Z41
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) Z42
! ->
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) RN
! ->
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) VJ0
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) VJX
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) VJY
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) VJZ
! ->
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) VMU0
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) VMU1
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) VMU2
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) VMU3
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) VMU4
! ->
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) NR0
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) NP0
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) NZ0
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) NRF
! ->
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) EPS
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) TAU
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) PHI0
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) RKC
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) ITM
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) ITP
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) ITS
!
CALL READSTR(1,STR,IER)
IF (IER .NE. 0) GOTO 999
READ (UNIT = STR,FMT = *) LP
!
LP = MAX0 (LP,LPP)
! ->
CLOSE (UNIT = 1)
! ->
PHI = DANGLE (VJX,VJY)
RHO = DSQRT (VJX * VJX + VJY * VJY)
!
CPH = DCOSM (PHI)
SPH = DSINM (PHI)
!
VJR = (+(CPH)) * VJX + SPH * VJY
VJP = (-(SPH)) * VJX + CPH * VJY
! ->
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' R11=',R11
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' R12=',R12
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' Z11=',Z11
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' Z12=',Z12
WRITE (UNIT = 2,FMT = *)
! ->
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' R21=',R21
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' R22=',R22
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' Z21=',Z21
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' Z22=',Z22
WRITE (UNIT = 2,FMT = *)
! ->
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' R31=',R31
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' R32=',R32
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' Z31=',Z31
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' Z32=',Z32
WRITE (UNIT = 2,FMT = *)
! ->
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' R41=',R41
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' R42=',R42
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' Z41=',Z41
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' Z42=',Z42
WRITE (UNIT = 2,FMT = *)
! ->
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' Rn=',RN
WRITE (UNIT = 2,FMT = *)
!
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' vJ0=',VJ0
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' vjx=',VJX
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' vjy=',VJY
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' vjz=',VJZ
WRITE (UNIT = 2,FMT = *)
!
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' phi=',PHI
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' rho=',RHO
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' vjr=',VJR
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' vjp=',VJP
WRITE (UNIT = 2,FMT = *)
!
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'vMu0=',VMU0
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'vMu1=',VMU1
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'vMu2=',VMU2
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'vMu3=',VMU3
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'vMu4=',VMU4
WRITE (UNIT = 2,FMT = *)
!
WRITE (UNIT = 2,FMT = '(1x,a5,1x,i11)') ' nr0=',NR0
WRITE (UNIT = 2,FMT = '(1x,a5,1x,i11)') ' np0=',NP0
WRITE (UNIT = 2,FMT = '(1x,a5,1x,i11)') ' nz0=',NZ0
WRITE (UNIT = 2,FMT = '(1x,a5,1x,i11)') ' nrf=',NRF
WRITE (UNIT = 2,FMT = *)
!
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' eps=',EPS
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' tau=',TAU
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'phi0=',PHI0
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' rkc=',RKC
WRITE (UNIT = 2,FMT = '(1x,a5,1x,i11)') ' itm=',ITM
WRITE (UNIT = 2,FMT = '(1x,a5,1x,i11)') ' itp=',ITP
WRITE (UNIT = 2,FMT = '(1x,a5,1x,i11)') ' its=',ITS
WRITE (UNIT = 2,FMT = '(1x,a5,1x,i11)') ' lp=',LP
WRITE (UNIT = 2,FMT = *)
! ---
WRITE (UNIT = *,FMT = *) '---> Normalization'
WRITE (UNIT = 2,FMT = *) '---> Normalization'
!
VJN = VJ0
VBN = VMU0 * VJ0
VPN = VJ0 * RN
!
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' vJn=',VJN
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' vBn=',VBN
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' vPn=',VPN
WRITE (UNIT = 2,FMT = *)
! ->
R11N = R11 / RN
R12N = R12 / RN
Z11N = Z11 / RN
Z12N = Z12 / RN
! ->
R21N = R21 / RN
R22N = R22 / RN
Z21N = Z21 / RN
Z22N = Z22 / RN
! ->
R31N = R31 / RN
R32N = R32 / RN
Z31N = Z31 / RN
Z32N = Z32 / RN
! ->
R41N = R41 / RN
R42N = R42 / RN
Z41N = Z41 / RN
Z42N = Z42 / RN
! ->
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'R11n=',R11N
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'R12n=',R12N
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'Z11n=',Z11N
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'Z12n=',Z12N
WRITE (UNIT = 2,FMT = *)
! ->
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'R21n=',R21N
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'R22n=',R22N
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'Z21n=',Z21N
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'Z22n=',Z22N
WRITE (UNIT = 2,FMT = *)
! ->
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'R31n=',R31N
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'R32n=',R32N
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'Z31n=',Z31N
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'Z32n=',Z32N
WRITE (UNIT = 2,FMT = *)
! ->
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'R41n=',R41N
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'R42n=',R42N
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'Z41n=',Z41N
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') 'Z42n=',Z42N
WRITE (UNIT = 2,FMT = *)
! ->
RA = DMIN1 (R11N,R21N,R31N,R41N)
RB = DMAX1 (R12N,R22N,R32N,R42N)
!
ZA = DMIN1 (Z11N,Z21N,Z31N,Z41N)
ZB = DMAX1 (Z12N,Z22N,Z32N,Z42N)
! ->
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' RA=',RA
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' RB=',RB
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' ZA=',ZA
WRITE (UNIT = 2,FMT = '(1x,a5,1x,1pe15.8)') ' ZB=',ZB
WRITE (UNIT = 2,FMT = *)
! ->
VMUMAT(1,1) = VMU1
VMUMAT(1,2) = 2D0 * VMU1 * VMU2 / (VMU1 + VMU2)
VMUMAT(1,3) = 2D0 * VMU1 * VMU3 / (VMU1 + VMU3)
VMUMAT(1,4) = 2D0 * VMU1 * VMU4 / (VMU1 + VMU4)
!
VMUMAT(2,1) = VMUMAT(1,2)
VMUMAT(2,2) = VMU2
VMUMAT(2,3) = 2D0 * VMU2 * VMU3 / (VMU2 + VMU3)
VMUMAT(2,4) = 2D0 * VMU2 * VMU4 / (VMU2 + VMU4)
!
VMUMAT(3,1) = VMUMAT(1,3)
VMUMAT(3,2) = VMUMAT(2,3)
VMUMAT(3,3) = VMU3
VMUMAT(3,4) = 2D0 * VMU3 * VMU4 / (VMU3 + VMU4)
!
VMUMAT(4,1) = VMUMAT(1,4)
VMUMAT(4,2) = VMUMAT(2,4)
VMUMAT(4,3) = VMUMAT(3,4)
VMUMAT(4,4) = VMU4
! ->
DO K = 1,4
DO I = 1,4
WRITE (UNIT = 2,FMT = '(1x,a4,2(i2),1x,1pe15.8)') 'vMu:',I,K
&,VMUMAT(I,K)
ENDDO
ENDDO
WRITE (UNIT = 2,FMT = *)
! ->
VJRM(1) = 0D0
VJRM(2) = (+(VJR))
VJRM(3) = (-(VJR))
VJRM(4) = 0D0
! ->
VJZM(1) = 0D0
VJZM(2) = (+(VJZ))
VJZM(3) = (-(VJZ))
VJZM(4) = 0D0
! ---
WRITE (UNIT = *,FMT = *) '---> Memory'
WRITE (UNIT = 2,FMT = *) '---> Memory'
! ---
!$SPF ANALYSIS(PARAMETER(NR0=130,NRF=2))
NR = NR0 * NRF
!$SPF ANALYSIS(PARAMETER(NP0=60))
NP = NP0 * NRF
!$SPF ANALYSIS(PARAMETER(NZ0=300))
NZ = NZ0 * NRF
!
WRITE (UNIT = *,FMT = '(1x,a4,3(1x,i8))') 'dim:',NR,NP,NZ
WRITE (UNIT = 2,FMT = '(1x,a4,3(1x,i8))') 'dim:',NR,NP,NZ
WRITE (UNIT = 2,FMT = *)
! ---
ALLOCATE(MAP(NR,NZ), STAT = IER)
!
ALLOCATE(RR(0:NR),PP(0:NP),ZZ(0:NZ), STAT = IER)
!
ALLOCATE(RS(NR),PS(NP),ZS(NZ), STAT = IER)
!
ALLOCATE(PHI1(NR,NP,NZ),PHI2(NR,NP,NZ),PHI3(NR,NP,NZ), STAT = IER)
!
ALLOCATE(ARX_RZ(NR,NZ),ARY_RZ(NR,NZ),ARZ_RZ(NR,NZ), STAT = IER)
!
ALLOCATE(AR_RZ(NR,NZ),AR_RP(NR,NP), STAT = IER)
!
ALLOCATE(HX(NR,NP,NZ),HY(NR,NP,NZ),HZ(NR,NP,NZ), STAT = IER)
ALLOCATE(HR(NR,NP,NZ),HP(NR,NP,NZ), STAT = IER)
!
ALLOCATE(BX(NR,NP,NZ),BY(NR,NP,NZ),BZ(NR,NP,NZ), STAT = IER)
ALLOCATE(BR(NR,NP,NZ),BP(NR,NP,NZ), STAT = IER)
! ---
WRITE (UNIT = *,FMT = *) '---> Grid & map'
WRITE (UNIT = 2,FMT = *) '---> Grid & map'
! ->
PI = 4D0 * DATAN (1D0)
!
DR = (RB - RA) / NR
DP = 2D0 * PI / NP
DZ = (ZB - ZA) / NZ
DM = DMIN1 (DR,DP,DZ)
! ->
WRITE (UNIT = *,FMT = '(1x,a4,4(1x,1pe15.8))') 'stp:',DR,DP,DZ,DM
WRITE (UNIT = 2,FMT = '(1x,a4,4(1x,1pe15.8))') 'stp:',DR,DP,DZ,DM
WRITE (UNIT = 2,FMT = *)
! ->
DM2 = DM * DM
DR2 = DR * DR
DP2 = DP * DP
DZ2 = DZ * DZ
!
CFR = DM2 / DR2
CFP = DM2 / DP2
CFZ = DM2 / DZ2
!
WRITE (UNIT = *,FMT = '(1x,a4,3(1x,1pe15.8))') 'cfs:',CFR,CFP,CFZ
WRITE (UNIT = 2,FMT = '(1x,a4,3(1x,1pe15.8))') 'cfs:',CFR,CFP,CFZ
! ->
DO I = 0,NR
RR(I) = RA + DR * I
ENDDO
!
DO J = 0,NP
PP(J) = DP * J
ENDDO
!
DO K = 0,NZ
ZZ(K) = ZA + DZ * K
ENDDO
! ->
DO I = 1,NR
RS(I) = 0.5D0 * (RR(I - 1) + RR(I))
ENDDO
!
DO J = 1,NP
PS(J) = 0.5D0 * (PP(J - 1) + PP(J))
ENDDO
!
DO K = 1,NZ
ZS(K) = 0.5D0 * (ZZ(K - 1) + ZZ(K))
ENDDO
! ->
K0 = NZ / 2
! ->
IF (LP .GT. 0) THEN
!
WRITE (UNIT = 2,FMT = '(5x,a1,8x,a2)') 'i','rr'
DO I = 0,NR
WRITE (UNIT = 2,FMT = '(1x,i8,1x,1pe15.8)') I,RR(I)
ENDDO
WRITE (UNIT = 2,FMT = *)
!
WRITE (UNIT = 2,FMT = '(5x,a1,8x,a2)') 'j','pp'
DO J = 0,NP
WRITE (UNIT = 2,FMT = '(1x,i8,1x,1pe15.8)') J,PP(J)
ENDDO
WRITE (UNIT = 2,FMT = *)
!
WRITE (UNIT = 2,FMT = '(5x,a1,8x,a2)') 'k','zz'
DO K = 0,NZ
WRITE (UNIT = 2,FMT = '(1x,i8,1x,1pe15.8)') K,ZZ(K)
ENDDO
WRITE (UNIT = 2,FMT = *)
!
WRITE (UNIT = 2,FMT = '(5x,a1,8x,a2)') 'i','rs'
DO I = 1,NR
WRITE (UNIT = 2,FMT = '(1x,i8,1x,1pe15.8)') I,RS(I)
ENDDO
WRITE (UNIT = 2,FMT = *)
!
WRITE (UNIT = 2,FMT = '(5x,a1,8x,a2)') 'j','ps'
DO J = 1,NP
WRITE (UNIT = 2,FMT = '(1x,i8,1x,1pe15.8)') J,PS(J)
ENDDO
WRITE (UNIT = 2,FMT = *)
!
WRITE (UNIT = 2,FMT = '(5x,a1,8x,a2)') 'k','zs'
DO K = 1,NZ
WRITE (UNIT = 2,FMT = '(1x,i8,1x,1pe15.8)') K,ZS(K)
ENDDO
WRITE (UNIT = 2,FMT = *)
!
ENDIF
! ->
DO K = 1,NZ
Z = ZS(K)
!
DO I = 1,NR
R = RS(I)
!
M = 4
!
IF (IN_RANGE_D (Z11N,Z12N,Z) .EQ. 0) THEN
IF (IN_RANGE_D (R11N,R12N,R) .EQ. 0) THEN
M = 1
GOTO 5
ENDIF
ENDIF
!
IF (IN_RANGE_D (Z21N,Z22N,Z) .EQ. 0) THEN
IF (IN_RANGE_D (R21N,R22N,R) .EQ. 0) THEN
M = 2
GOTO 5
ENDIF
ENDIF
!
IF (IN_RANGE_D (Z31N,Z32N,Z) .EQ. 0) THEN
IF (IN_RANGE_D (R31N,R32N,R) .EQ. 0) THEN
M = 3
GOTO 5
ENDIF
ENDIF
!
5 MAP(I,K) = M
!
ENDDO
ENDDO
! ->
CALL TPL2D_REC_I(10,'magint_3d_map.plt','TITLE = "MAP"','VARIABLES
& = "R" "Z" "MAP"',NR,NZ,RS,ZS,MAP)
! ->
CALL TPL2D_REC_IM(10,'magint_3d_map_xz.plt','TITLE = "MAP"','VARIA
&BLES = "X" "Z" "MAP"',NR,NZ,RS,ZS,MAP)
! ->
DO K = 1,NZ
DO I = 1,NR
M = MAP(I,K)
AR_RZ(I,K) = VMUMAT(M,M)
ENDDO
ENDDO
! ->
CALL TPL2D_REC_D(10,'magint_3d_vMu.plt','TITLE = "vMu"','VARIABLES
& = "R" "Z" "vMu"',NR,NZ,RS,ZS,AR_RZ)
! ---
WRITE (UNIT = *,FMT = *) '---> Iterative loop'
WRITE (UNIT = 2,FMT = *) '---> Iterative loop'
! ---
CPU_T2 = OMP_GET_WTIME ()
! ---
ZERO = 1D-15
! ---
IT = 0
!
IT1 = 0
RK1 = 1D0
! ->
DO K = 1,NZ
DO J = 1,NP
DO I = 1,NR
PHI1(I,J,K) = PHI0
PHI2(I,J,K) = PHI0
PHI3(I,J,K) = PHI0
ENDDO
ENDDO
ENDDO
! ->
!_omp parallel num_threads(nthr) private (nth,mth)
!OMP_GET_NUM_THREADS ()
NTH = 1
MTH = 0
!OMP_GET_THREAD_NUM ()
CALL MAIN_CALC(NTH,MTH,NR,NP,NZ,MAP,RR,RS,PHI1,PHI2,PHI3,RK,RK1,IT
&,IT1,LOCK)
!
!_omp end parallel
!
! --- results ---
!
20 CONTINUE
!
WRITE (UNIT = *,FMT = *) 'final result:'
WRITE (UNIT = 2,FMT = *) 'final result:'
!
WRITE (UNIT = *,FMT = '(4h it=,i8,4h rk=,1x,1pe12.5)') IT1,RK1
WRITE (UNIT = 2,FMT = '(4h it=,i8,4h rk=,1x,1pe12.5)') IT1,RK1
! ->
DO K = 1,NZ
DO I = 1,NR
AR_RZ(I,K) = PHI3(I,1,K)
ENDDO
ENDDO
!
DO J = 1,NP
DO I = 1,NR
AR_RP(I,J) = PHI3(I,J,K0)
ENDDO
ENDDO
!
CALL TPL2D_REC_D(10,'magint_3d_phi_rz.plt','TITLE = "phi*(r,z)"','
&VARIABLE = "R" "Z" "PHI"',NR,NZ,RS,ZS,AR_RZ)
!
CALL TPL2D_REC_DM(10,'magint_3d_phi_xz.plt','TITLE = "phi*(x,0,z)"
&','VARIABLE = "X" "Z" "PHI"',NR,NZ,RS,ZS,AR_RZ)
!
CALL TPL2D_REC_D(10,'magint_3d_phi_rp.plt','TITLE = "phi*(r,p)"','
&VARIABLE = "R" "P" "PHI"',NR,NP,RS,PS,AR_RP)
! ->
DO K = 1,NZ
DO J = 1,NP
DO I = 1,NR
! -
IF (I .EQ. 1) THEN
A1 = RS(I)
A2 = RS(I + 1)
A3 = RS(I + 2)
V1 = PHI3(I,J,K)
V2 = PHI3(I + 1,J,K)
V3 = PHI3(I + 2,J,K)
VP = FP2P1 (A1,A2,A3,V1,V2,V3,A1)
ELSE
IF (I .EQ. NR) THEN
A1 = RS(I - 2)
A2 = RS(I - 1)
A3 = RS(I)
V1 = PHI3(I - 2,J,K)
V2 = PHI3(I - 1,J,K)
V3 = PHI3(I,J,K)
VP = FP2P1 (A1,A2,A3,V1,V2,V3,A3)
ELSE
A1 = RS(I - 1)
A2 = RS(I)
A3 = RS(I + 1)
V1 = PHI3(I - 1,J,K)
V2 = PHI3(I,J,K)
V3 = PHI3(I + 1,J,K)
VP = FP2P1 (A1,A2,A3,V1,V2,V3,A2)
ENDIF
ENDIF
!
VHR = (-(VP))
! -
IF (J .EQ. 1) THEN
A1 = PS(NP)
V1 = PHI3(I,NP,K)
ELSE
A1 = PS(J - 1)
V1 = PHI3(I,J - 1,K)
ENDIF
!
A2 = PS(J)
V2 = PHI3(I,J,K)
!
IF (J .EQ. NP) THEN
A3 = PS(1)
V3 = PHI3(I,1,K)
ELSE
A3 = PS(J + 1)
V3 = PHI3(I,J + 1,K)
ENDIF
!
VP = FP2P1 (A1,A2,A3,V1,V2,V3,A2)
!
VHP = (-(VP)) / RS(J)
! -
IF (K .EQ. 1) THEN
A1 = ZS(K)
A2 = ZS(K + 1)
A3 = ZS(K + 2)
V1 = PHI3(I,J,K)
V2 = PHI3(I,J,K + 1)
V3 = PHI3(I,J,K + 2)
VP = FP2P1 (A1,A2,A3,V1,V2,V3,A1)
ELSE
IF (K .EQ. NZ) THEN
A1 = ZS(K - 2)
A2 = ZS(K - 1)
A3 = ZS(K)
V1 = PHI3(I,J,K - 2)
V2 = PHI3(I,J,K - 1)
V3 = PHI3(I,J,K)
VP = FP2P1 (A1,A2,A3,V1,V2,V3,A3)
ELSE
A1 = ZS(K - 1)
A2 = ZS(K)
A3 = ZS(K + 1)
V1 = PHI3(I,J,K - 1)
V2 = PHI3(I,J,K)
V3 = PHI3(I,J,K + 1)
VP = FP2P1 (A1,A2,A3,V1,V2,V3,A2)
ENDIF
ENDIF
!
VHZ = (-(VP))
! -
CPH = DCOSM (PS(J))
SPH = DSINM (PS(J))
!
VHX = VHR * CPH - VHP * SPH
VHY = VHR * SPH + VHP * CPH
! -
M = MAP(I,K)
! -
IF (M .EQ. 2 .OR. M .EQ. 3) THEN
VBX = VHX + VJX
VBY = VHY + VJY
VBZ = VHZ + VJZM(M)
VBR = VHR + VJR
VBP = VHP + VJP
ELSE
VMU = VMUMAT(M,M)
VBX = VMU * VHX
VBY = VMU * VHY
VBZ = VMU * VHZ
VBR = VMU * VHR
VBP = VMU * VHP
ENDIF
! -
HX(I,J,K) = VHX
HY(I,J,K) = VHY
HZ(I,J,K) = VHZ
HR(I,J,K) = VHR
HP(I,J,K) = VHP
! -
BX(I,J,K) = VBX
BY(I,J,K) = VBY
BZ(I,J,K) = VBZ
BR(I,J,K) = VBR
BP(I,J,K) = VBP
! -
IF (LP .GT. 0) THEN
!
WRITE (UNIT = 2,FMT = '(1x,a4,4(1x,i8))') 'ijkm',I,J,K
&,M
!
WRITE (UNIT = 2,FMT = '(1x,a4,5(1x,1pe13.6))') ' vH:',
&VHX,VHY,VHZ,VHR,VHP
!
WRITE (UNIT = 2,FMT = '(1x,a4,5(1x,1pe13.6))') ' vB:',
&VBX,VBY,VBZ,VBR,VBP
!
ENDIF
! -
ENDDO
ENDDO
ENDDO
! ->
DO K = 1,NZ
DO I = 1,NR
ARX_RZ(I,K) = HX(I,1,K)
ARY_RZ(I,K) = HY(I,1,K)
ARZ_RZ(I,K) = HZ(I,1,K)
ENDDO
ENDDO
!
CALL TPL2D_REC_D(10,'magint_3d_Hz_rz.plt','TITLE = "Hz(r,z)"','VAR
&IABLE = "R" "Z" "Hz"',NR,NZ,RS,ZS,ARZ_RZ)
!
CALL TPL2D_REC_DM3(10,'magint_3d_Hv_xz.plt','TITLE = "Hv(x,0,z)"',
&'VARIABLE = "X" "Z" "Hx" "Hy" "Hz" "Hm"',NR,NZ,RS,ZS,ARX_RZ,ARY_RZ
&,ARZ_RZ)
! ->
DO K = 1,NZ
DO I = 1,NR
ARX_RZ(I,K) = BX(I,1,K)
ARY_RZ(I,K) = BY(I,1,K)
ARZ_RZ(I,K) = BZ(I,1,K)
ENDDO
ENDDO
!
CALL TPL2D_REC_D(10,'magint_3d_Bz_rz.plt','TITLE = "Bz(r,z)"','VAR
&IABLE = "R" "Z" "Bz"',NR,NZ,RS,ZS,ARZ_RZ)
!
CALL TPL2D_REC_DM3(10,'magint_3d_Bv_xz.plt','TITLE = "Bv(x,0,z)"',
&'VARIABLE = "X" "Z" "Bx" "By" "Bz" "Bm"',NR,NZ,RS,ZS,ARX_RZ,ARY_RZ
&,ARZ_RZ)
! ---
CPU_T3 = OMP_GET_WTIME ()
!
CPU_DT1 = CPU_T2 - CPU_T1
CPU_DT2 = CPU_T3 - CPU_T2
!
WRITE (UNIT = 2,FMT = '(11h cpu_times:,2(1x,f10.3))') CPU_DT1,CPU_
&DT2
WRITE (UNIT = *,FMT = '(11h cpu_times:,2(1x,f10.3))') CPU_DT1,CPU_
&DT2
! ---
999 CONTINUE
!
CLOSE (UNIT = 2)
! ---
IF (ALLOCATED (MAP)) DEALLOCATE(MAP)
IF (ALLOCATED (RR)) DEALLOCATE(RR)
IF (ALLOCATED (PP)) DEALLOCATE(PP)
IF (ALLOCATED (ZZ)) DEALLOCATE(ZZ)
IF (ALLOCATED (RR)) DEALLOCATE(RS)
IF (ALLOCATED (PP)) DEALLOCATE(PS)
IF (ALLOCATED (ZZ)) DEALLOCATE(ZS)
IF (ALLOCATED (PHI1)) DEALLOCATE(PHI1)
IF (ALLOCATED (PHI2)) DEALLOCATE(PHI2)
IF (ALLOCATED (PHI3)) DEALLOCATE(PHI3)
IF (ALLOCATED (AR_RZ)) DEALLOCATE(AR_RZ)
IF (ALLOCATED (AR_RP)) DEALLOCATE(AR_RP)
IF (ALLOCATED (ARX_RZ)) DEALLOCATE(ARX_RZ)
IF (ALLOCATED (ARY_RZ)) DEALLOCATE(ARY_RZ)
IF (ALLOCATED (ARZ_RZ)) DEALLOCATE(ARZ_RZ)
IF (ALLOCATED (HX)) DEALLOCATE(HX)
IF (ALLOCATED (HY)) DEALLOCATE(HY)
IF (ALLOCATED (HZ)) DEALLOCATE(HZ)
IF (ALLOCATED (HR)) DEALLOCATE(HR)
IF (ALLOCATED (HP)) DEALLOCATE(HP)
IF (ALLOCATED (BX)) DEALLOCATE(BX)
IF (ALLOCATED (BY)) DEALLOCATE(BY)
IF (ALLOCATED (BZ)) DEALLOCATE(BZ)
IF (ALLOCATED (BR)) DEALLOCATE(BR)
IF (ALLOCATED (BP)) DEALLOCATE(BP)
! ---
STOP
END
! ----------------------------------------------------------------------
SUBROUTINE MAIN_CALC (NTH, MTH, NR, NP, NZ, MAP, RR, RS, PHI1, PHI
&2, PHI3, RK, RK1, IT, IT1, LOCK)
!
IMPLICIT INTEGER*4 (I-N)
IMPLICIT REAL*8 (A-H,O-Z)
!
! include 'omp_lib.h'
!
! INTEGER(KIND=OMP_LOCK_KIND) :: LOCK
!
INTEGER*4 :: MAP(NR,NZ)
!
REAL*8 :: RR(0:NR),RS(NR)
REAL*8 :: PHI1(NR,NP,NZ)
REAL*8 :: PHI2(NR,NP,NZ)
REAL*8 :: PHI3(NR,NP,NZ)
!
! ---
COMMON /BMAIN/DR,DP,DZ,CFR,CFP,CFZ
COMMON /BMATR/VMUMAT(4,4),VJRM(4),VJZM(4)
COMMON /BITER/EPS,TAU,RKC,PHI0,ITM,ITP,ITS,K0
! ---
CALL SECT1D(NTH,MTH,1,NZ,K1,K2)
!
! write(*,*) 'CALC_1:', nth, mth, k1, k2
! ---
10 CONTINUE
! ->
IF (MTH .EQ. 0) THEN
RK = 0D0
ENDIF
! ->
RK0 = 0D0
! ---
!$SPF ANALYSIS(PRIVATE(VMU,VIJMK,VIJPK))
DO K = K1,K2
DO J = 1,NP
DO I = 1,NR
! ->
MA_0 = MAP(I,K)
!
MA_L = 0
MA_R = 0
MA_B = 0
MA_T = 0
!
!- left
IF (I .GT. 1) MA_L = MAP(I - 1,K)
!- right
IF (I .LT. NR) MA_R = MAP(I + 1,K)
!- bottom
IF (K .GT. 1) MA_B = MAP(I,K - 1)
!- top
IF (K .LT. NZ) MA_T = MAP(I,K + 1)
! ->
VIJK = PHI1(I,J,K)
! ->
DIAG = 0D0
! ->
WRM = 0D0
WRP = 0D0
!
IF (I .GT. 1) THEN
!- Jacobi
VIMJK = PHI1(I - 1,J,K)
S1 = CFR * RR(I - 1) / RS(I)
IF (MA_0 .EQ. MA_L) THEN
VMU = VMUMAT(MA_0,MA_0)
WRM = S1 * VMU * (VIMJK - VIJK)
ELSE
VMU = VMUMAT(MA_L,MA_0)
WW1 = 0.5D0 * DR * VJRM(MA_L) / VMUMAT(MA_L,MA_L)
WW2 = 0.5D0 * DR * VJRM(MA_0) / VMUMAT(MA_0,MA_0)
WWW = WW1 + WW2
WRM = S1 * VMU * (VIMJK - VIJK + WWW)
ENDIF
DIAG = DIAG + S1 * VMU
ENDIF
!
IF (I .LT. NR) THEN
VIPJK = PHI1(I + 1,J,K)
S2 = CFR * RR(I) / RS(I)
IF (MA_0 .EQ. MA_R) THEN
VMU = VMUMAT(MA_0,MA_0)
WRP = S2 * VMU * (VIPJK - VIJK)
ELSE
VMU = VMUMAT(MA_0,MA_R)
WW1 = 0.5D0 * DR * VJRM(MA_0) / VMUMAT(MA_0,MA_0)
WW2 = 0.5D0 * DR * VJRM(MA_R) / VMUMAT(MA_R,MA_R)
WWW = WW1 + WW2
WRP = S2 * VMU * (VIPJK - VIJK - WWW)
ENDIF
DIAG = DIAG + S2 * VMU
ENDIF
! ->
IF (J .GT. 1) THEN
!- Jacobi
VIJMK = PHI1(I,J - 1,K)
ELSE
VIJMK = PHI1(I,NP,K)
ENDIF
!
IF (J .LT. NP) THEN
VIJPK = PHI1(I,J + 1,K)
ELSE
VIJPK = PHI1(I,1,K)
ENDIF
!
RS2 = RS(I) * RS(I)
VMU = VMUMAT(MA_0,MA_0)
S3 = CFP / RS2
WPM = S3 * VMU * (VIJMK - VIJK)
WPP = S3 * VMU * (VIJPK - VIJK)
DIAG = DIAG + S3 * VMU
! ->
WZM = 0D0
WZP = 0D0
!
IF (K .GT. 1) THEN
!- Jacobi
VIJKM = PHI1(I,J,K - 1)
IF (MA_0 .EQ. MA_B) THEN
VMU = VMUMAT(MA_0,MA_0)
WZM = CFZ * VMU * (VIJKM - VIJK)
ELSE
VMU = VMUMAT(MA_B,MA_0)
WW1 = 0.5D0 * DZ * VJZM(MA_B) / VMUMAT(MA_B,MA_B)
WW2 = 0.5D0 * DZ * VJZM(MA_0) / VMUMAT(MA_0,MA_0)
WWW = WW1 + WW2
WZM = CFZ * VMU * (VIJKM - VIJK + WWW)
ENDIF
DIAG = DIAG + CFZ * VMU
ENDIF
!
IF (K .LT. NZ) THEN
VIJKP = PHI1(I,J,K + 1)
IF (MA_0 .EQ. MA_T) THEN
VMU = VMUMAT(MA_0,MA_0)
WZP = CFZ * VMU * (VIJKP - VIJK)
ELSE
VMU = VMUMAT(MA_0,MA_T)
WW1 = 0.5D0 * DZ * VJZM(MA_0) / VMUMAT(MA_0,MA_0)
WW2 = 0.5D0 * DZ * VJZM(MA_T) / VMUMAT(MA_T,MA_T)
WWW = WW1 + WW2
WZP = CFZ * VMU * (VIJKP - VIJK - WWW)
DIAG = DIAG + CFZ * VMU
ENDIF
ENDIF
! ->
S0 = WRP + WRM + WPP + WPM + WZP + WZM
RK0 = DMAX1 (RK0,S0)
!
PHI2(I,J,K) = VIJK + TAU * (S0 / DIAG)
!
ENDDO
ENDDO
ENDDO
! ->
! write(*,*) 'CALC_2:', mth, rk0
! ---
!_omp barrier
! CALL OMP_SET_LOCK(LOCK)
RK = DMAX1 (RK,RK0)
! CALL OMP_UNSET_LOCK(LOCK)
!_omp barrier
! ---
IF (MTH .EQ. 0) THEN
IF (MOD (IT,ITP) .EQ. 0 .OR. RK .LE. EPS .OR. IT .EQ. ITM) THEN
WRITE (UNIT = *,FMT = '(4h it=,i8,4h rk=,1x,1pe12.5,1x,1pe12
&.5)') IT,RK,RK1
WRITE (UNIT = 2,FMT = '(4h it=,i8,4h rk=,1x,1pe12.5,1x,1pe12
&.5)') IT,RK,RK1
ENDIF
ENDIF
! ->
IF (RK .LT. RK1) THEN
IF (MTH .EQ. 0) THEN
RK1 = RK
IT1 = IT
ENDIF
!
DO K = K1,K2
DO J = 1,NP
DO I = 1,NR
S1 = PHI1(I,J,K)
PHI3(I,J,K) = S1
ENDDO
ENDDO
ENDDO
ENDIF
! ->
IF (RK .LE. EPS) GOTO 20
IF (IT .GE. ITM) GOTO 20
IF (RK .GT. RKC * RK1) GOTO 20
! ->
IF (MTH .EQ. 0) THEN
IT = IT + 1
ENDIF
! ->
DO K = K1,K2
DO J = 1,NP
DO I = 1,NR
S2 = PHI2(I,J,K)
PHI1(I,J,K) = S2
ENDDO
ENDDO
ENDDO
! ->
IF (MTH .EQ. 0) THEN
PHI1(1,1,K0) = PHI0
ENDIF
! ->
!_omp barrier
! ->
GOTO 10
! ->
20 CONTINUE
! ->
!_omp barrier
! ---
RETURN
END
! ----------------------------------------------------------------------
! Tecplot 2D, integer*4, (x,z) co-ordinates
! ----------------------------------------------------------------------
SUBROUTINE TPL2D_REC_IM (ICH, FNAME, TNAME, VNAME, N1, N2, X1, X2,
& X3)
!
IMPLICIT INTEGER*4 (I-N)
IMPLICIT REAL*8 (A-H,O-Z)
!
CHARACTER*(*) :: FNAME,TNAME,VNAME
!
REAL*8 :: X1(N1),X2(N2)
INTEGER*4 :: X3(N1,N2)
!
OPEN (UNIT = ICH,FILE = FNAME)
!
WRITE (UNIT = ICH,FMT = *) TNAME
WRITE (UNIT = ICH,FMT = *) VNAME
WRITE (UNIT = ICH,FMT = *) 'ZONE T="BIG" I=',2 * N1,' J=',N2,' F=P
&OINT'
!
DO J = 1,N2
V2 = X2(J)
DO I = N1,1,(-(1))
V1 = (-(X1(I)))
V3 = 1D0 * X3(I,J)
WRITE (UNIT = ICH,FMT = '(3(1x,1pe15.8))') V1,V2,V3
ENDDO
DO I = 1,N1
V1 = X1(I)
V3 = 1D0 * X3(I,J)
WRITE (UNIT = ICH,FMT = '(3(1x,1pe15.8))') V1,V2,V3
ENDDO
ENDDO
!
CLOSE (UNIT = ICH)
!
RETURN
END
! ----------------------------------------------------------------------
! Tecplot 2D, real*8, (x,z) co-ordinates, output scalar
! ----------------------------------------------------------------------
SUBROUTINE TPL2D_REC_DM (ICH, FNAME, TNAME, VNAME, N1, N2, X1, X2,
& X3)
!
IMPLICIT INTEGER*4 (I-N)
IMPLICIT REAL*8 (A-H,O-Z)
!
CHARACTER*(*) :: FNAME,TNAME,VNAME
!
REAL*8 :: X1(N1),X2(N2)
REAL*8 :: X3(N1,N2)
!
OPEN (UNIT = ICH,FILE = FNAME)
!
WRITE (UNIT = ICH,FMT = *) TNAME
WRITE (UNIT = ICH,FMT = *) VNAME
WRITE (UNIT = ICH,FMT = *) 'ZONE T="BIG" I=',2 * N1,' J=',N2,' F=P
&OINT'
!
DO J = 1,N2
V2 = X2(J)
DO I = N1,1,(-(1))
V1 = (-(X1(I)))
V3 = X3(I,J)
WRITE (UNIT = ICH,FMT = '(3(1x,1pe15.8))') V1,V2,V3
ENDDO
DO I = 1,N1
V1 = X1(I)
V3 = X3(I,J)
WRITE (UNIT = ICH,FMT = '(3(1x,1pe15.8))') V1,V2,V3
ENDDO
ENDDO
!
CLOSE (UNIT = ICH)
!
RETURN
END
! ----------------------------------------------------------------------
! Tecplot 2D, real*8, (x,z) co-ordinates, output of 3D-vector
! ----------------------------------------------------------------------
SUBROUTINE TPL2D_REC_DM3 (ICH, FNAME, TNAME, VNAME, N1, N2, X1, X2
&, X3, X4, X5)
!
IMPLICIT INTEGER*4 (I-N)
IMPLICIT REAL*8 (A-H,O-Z)
!
CHARACTER*(*) :: FNAME,TNAME,VNAME
!
REAL*8 :: X1(N1),X2(N2)
REAL*8 :: X3(N1,N2),X4(N1,N2),X5(N1,N2)
!
OPEN (UNIT = ICH,FILE = FNAME)
!
WRITE (UNIT = ICH,FMT = *) TNAME
WRITE (UNIT = ICH,FMT = *) VNAME
WRITE (UNIT = ICH,FMT = *) 'ZONE T="BIG" I=',2 * N1,' J=',N2,' F=P
&OINT'
!
DO J = 1,N2
V2 = X2(J)
DO I = N1,1,(-(1))
V1 = (-(X1(I)))
V3 = (-(X3(I,J)))
V4 = X4(I,J)
V5 = X5(I,J)
V6 = DSQRTM (V3 * V3 + V4 * V4 + V5 * V5)
IF (V6 .GE. 1E-15) THEN
V3 = V3 / V6
V4 = V4 / V6
V5 = V5 / V6
ENDIF
WRITE (UNIT = ICH,FMT = '(6(1x,1pe15.8))') V1,V2,V3,V4,V5,V6
ENDDO
DO I = 1,N1
V1 = X1(I)
V3 = X3(I,J)
V4 = X4(I,J)
V5 = X5(I,J)
V6 = DSQRTM (V3 * V3 + V4 * V4 + V5 * V5)
IF (V6 .GE. 1E-15) THEN
V3 = V3 / V6
V4 = V4 / V6
V5 = V5 / V6
ENDIF
WRITE (UNIT = ICH,FMT = '(6(1x,1pe15.8))') V1,V2,V3,V4,V5,V6
ENDDO
ENDDO
!
CLOSE (UNIT = ICH)
!
RETURN
END
Reference in New Issue Copy Permalink