! ---------------------------------------------------------------------- ! 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 [ [ &]]' 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