IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) c10================================================================ c10 RAY TRACING PROGRAM FOR VLF WAVE (X-MODE) IN DE-MODEL c10 Dipole model is selected in REF by Call dipole c10 IGRF model is selected in REF by Call IGRF c10 kizami H is shown in Adam_update c10 Result of ray tracing (OUT1) is written in fort.21 c10 Instantaneous H is written in fort.25 c10======== initial conditions ================================== IYEAR=1990 IMNTH=1 IDATE=1 CALL IGRF90(IYEAR,IMNTH,IDATE) c fq00: wave frequency in kHz c rr00: starting radius distance from the earth center in km c th00: geomagnetic colatitude in deg c ph00: geomagnetic longitude in deg c es00: angle between the initial wave normal and geomagnetic c meiridian plane, in deg (Epsilon) c dl00: angle between the projected direction of the initial c wave normal onto the geomagnetic meridial plane and c the radial direction, in deg (Delta) rr00 = 6370.0+91.0d0 ph00 = 0.0d0 ameg = 1000.0d0 fq00 = 0.003d0*ameg dl00 = 0.0d0 es00 = 0.0d0 print*,'VLF_DE_IGRF_F.f',', rr00=', rr00,', fq00=',fq00, &', ph00=', ph00,', es00=', es00,', dl00=',dl00 DO 100 I=1,5 WRITE(21,*) write(25,*) th00 = dble(i-1)*1.0d0+40.0d0 print*, 'th00=', th00 CALL INIT(fq00,rr00,th00,ph00,es00,dl00) CALL OUT0 CALL ADAMS 100 CONTINUE STOP END ************************************************************************ SUBROUTINE INIT(fq00,rr00,th00,ph00,es00,dl00) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QINIT /R0,T0,P0,E0,D0,F0 save inited data inited / 0 / PAI = dASIN(1.0D0)*2.0D0 RAD = PAI/180.0D0 DEG = 1.0D0/RAD RE = 6370.0D0 RREF = 7370.0D0 RREFP = RREF RBTM = 9370.0D0 CVLCTY = 3.00D5 XME = 9.109534D-31 XMP = 1.6726485D-27 QE = 1.6021892D-19 EP = 8.854185D-12 XK = 1.380662D-23 GM = 3.98603D14 PF0(1) = QE**2/4.0D0/PAI**2/EP/ XME PF0(2) = QE**2/4.0D0/PAI**2/EP/ XMP PF0(3) = QE**2/4.0D0/PAI**2/EP/(XMP*4.0D0) PF0(4) = QE**2/4.0D0/PAI**2/EP/(XMP*16.0D0) GF0(1) =-QE/2.0D0/PAI/ XME *1.0D-3 GF0(2) = QE/2.0D0/PAI/ XMP *1.0D-3 GF0(3) = QE/2.0D0/PAI/(XMP* 4.0D0)*1.0D-3 GF0(4) = QE/2.0D0/PAI/(XMP*16.0D0)*1.0D-3 c XN0 = 32000.0D0 XN0 = 50000.0D0 c XN0 = 4000.0D0 THERM = 1000.0D0 SHINV(2) = XMP* (GM/RREF**2)/XK/THERM*1.0D-3 SHINV(3) = XMP* 4.0D0*(GM/RREF**2)/XK/THERM*1.0D-3 SHINV(4) = XMP*16.0D0*(GM/RREF**2)/XK/THERM*1.0D-3 SHINVP = SHINV(2) ETA(2) = 0.152D0 ETA(3) = 0.82D0 ETA(4) = 0.025D0 c Following RRLI: location of the lower edge of the ionosphere RRLI = 6370.0d0 + 90.0d0 c The following parameters areused in subroutine DENSLI HWLI = 140.0D0 XLPP = 4.0D0 ERPP = 0.2D0 DRR1 = 5.0D0 DTH1 = 0.05D0 DPH1 = 0.2D0 DRR2 = 2.5D0 DTH2 = 0.01D0 DPH2 = 0.1D0 c Following RRMAX & RRMIN: Upper & lower limit of ray tracing(Kmj c Following ERRR and ERRA; Initial values for the limiting error c and Hmin: Limit of integration increment RRMAX = 4.0*RE RRMIN = 6370.0d0 ERRR = 1.0D-5 ERRA = 1.0D-5 HMIN = 1.0D-10 c Folowing DISOUT: Every this count, the output data to be printed DISOUT = 200 c Following NLPMAX: Maximum number of integration loops NLPMAX = 30000 INDOUT = 21 INDMSH = 51 IEOF = -1 C...FREQ.(KHZ),GC-DIST.(KM),GM-COLAT.,GM-LONG.,EPS,DEL(DEG) F0 = fq00 R0 = rr00 T0 = th00 P0 = ph00 E0 = es00 D0 = dl00 if(inited.eq.0) then open(indmsh,file='MESH',form='unformatted') CALL MESH (INDMSH,IYEAR,IMNTH,IDATE) close(indmsh) CALL IGRF90(IYEAR,IMNTH,IDATE) CALL GTOM0 inited = 1 endif RETURN END ************************************************************************ SUBROUTINE ADAMS ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) *** CONSTANT COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY *** INPUT PARAMETER (NSPEC=4) COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QINIT /R0,T0,P0,E0,D0,F0 COMMON /QPATH /FKC,DEL,EPS,GMU,ALF,RAY *** OUTPUT COMMON /QCOND /L,LL,NN,NLOOP,A *** Variables below are only used in the ADAMS method. COMMON /QADAM /Y(8),P(8),F(8), & H,HOLD,WT(8),PHI(8,14),PSIOLD(12), & PSI(12),BETA(12),SIGMA(13),G(13), & ERKM1,ERK,K,KOLD, & KNEW,KS,IPHASE,IFAIL REAL*8 SUML(8) CALL OUT0 *** Initialization of Y CALL MTOG1(T0,P0) Y(1) = R0 Y(2) = T0 Y(3) = P0 EPS = E0 DEL = D0 FKC = F0 Y(4) = COS(EPS*RAD)*COS(DEL*RAD) Y(5) = COS(EPS*RAD)*SIN(DEL*RAD) Y(6) = SIN(EPS*RAD) Y(7) = 0.0D0 Y(8) = 0.0D0 L = 0 LL = 0 NN = 0 NLOOP = 0 ERR0 = MAX(ERRR,ERRA) ERRR = ERRR/ERR0 ERRA = ERRA/ERR0 * Derivative for Y CALL FUNCT(Y,F,1) yp20 = F(2) do L=1,8 WT(L)=ERRR*dABS(Y(L)+ERRA) SUML(L) = (F(L)/WT(L))**2 enddo SUM = 0.0D0 do L=1,8 SUM = SUM+SUML(L) enddo SUM = DSQRT(SUM) H = 0.00001D0*dSQRT(ERR0/SUM) HOLD = 0.0D0 K = 1 KOLD = 0 IPHASE = 1 IFAIL = 0 do L=1,8 PHI(L,1) = F(L) PHI(L,2) = 0.0D0 enddo CALL OUT1(Y) DO NLOOP=0,NLPMAX-1 * Adams-Bashforth Predictor CALL ADAM_B(NLOOP) * Adams-Moulton Corrector CALL ADAM_M(ERR0,NLOOP) *VOCL LOOP,SCALAR IF(LL.EQ.0) THEN IF(IFAIL.NE.0) THEN L=11 ELSE IF(L.GE.21) THEN LL=L ELSE IF(H.LT.HMIN) THEN L =23 LL=23 ELSE IF(Y(1).LT.RRMIN) THEN L =24 LL=24 ELSE IF(Y(1).GT.RRMAX) THEN L =25 LL=25 C ELSE IF((T0-90.0)*F(2).GT.0.0D0) THEN C L =26 C LL=26 C ELSE IF(yp20*yp(2).lt.0.0D0) THEN C L =27 C LL=27 ELSE L =0 NN=NN+1 ENDIF ENDIF ENDIF CALL OUT1(Y) IF (LL.NE.0) THEN write(*,*) "Stop Condition: ", LL GOTO 999 ENDIF do I=1,7 WT(I)=ERRR*ABS(Y(I)+ERRA) enddo ENDDO 999 CONTINUE write(*,*) 'Number of LOOP=',NLOOP,', Number of calculation=',A write(*,*) 'Last interval=',H,', IFAIL=',IFAIL write(21,*) write(25,*) RETURN END ************************************************************************ SUBROUTINE ADAM_B(NLOOP) * Calculation of Adams-Bashforth Predictor ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) COMMON /QADAM /Y(8),P(8),F(8), & H,HOLD,WT(8),PHI(8,14),PSIOLD(12), & PSI(12),BETA(12),SIGMA(13),G(13), & ERKM1,ERK,K,KOLD, & KNEW,KS,IPHASE,IFAIL real*8 ALPHA(12),V(12),W(12) save ALPHA,V,W G(1) = 1.0D0 G(2) = 0.5D0 SIGMA(1) = 1.0D0 * KS is the maximum order * In case step size is changed, KS->0. IF (H .NE.HOLD) KS = 0 IF (KS.LE.KOLD) KS = KS + 1 IF (K.GE.KS) THEN *** PSI do I=KS,K-1 PSIOLD(I) = PSI(I) enddo PSI(KS) = H*KS do I=KS+1,K PSI(I) = PSIOLD(I-1)+H enddo *** ALPHA ALPHA(KS) = 1.0D0/DBLE(KS) do I=KS+1,K ALPHA(I) = H/(PSIOLD(I-1)+H) enddo *** BETA BETA(KS) = 1.0D0 do I=KS+1,K BETA(I) = BETA(I-1)*PSI(I-1)/PSIOLD(I-1) enddo *** SIGMA SIGMA(KS+1) = 1.0D0 do I=KS+1,K SIGMA(I+1) = I*ALPHA(I)*SIGMA(I) enddo *** G (Integral of C) if (KS.EQ.0) write(0,*) 'Error! KS=0' if (KS.EQ.1) then do I=1,K V(I) = 1.0D0/DBLE(I*(I+1)) W(I) = V(I) enddo else if (K.gt.KOLD) then V(K) = 1.0D0/DBLE(K*(K+1)) do J=1,KS-2 V(K-J) = V(K-J)-ALPHA(J+1)*V(K-J+1) enddo endif do IQ=1,K+1-KS V(IQ) = V(IQ)-ALPHA(KS)*V(IQ+1) W(IQ) = V(IQ) enddo G(KS+1) = W(1) endif do I=KS+1,K do IQ=1,K+1-I W(IQ) = W(IQ)-ALPHA(I)*W(IQ+1) enddo G(I+1) = W(1) enddo ENDIF *** PREDICTOR do L=1,8 do I=KS+1,K PHI(L,I) = BETA(I)*PHI(L,I) enddo PHI(L,K+2) = PHI(L,K+1) PHI(L,K+1) = 0.0D0 SUM = 0.0D0 do I=K,1, -1 SUM = SUM + G(I)*PHI(L,I) PHI(L,I) = PHI(L,I) + PHI(L,I+1) enddo P(L) = Y(L) + H * SUM enddo *** Derivative CALL FUNCT(P,F,0) RETURN END ************************************************************************ SUBROUTINE ADAM_M(ERR0,NLOOP) * Calculation of Adams-Moulton Corrector ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) COMMON /QADAM /Y(8),P(8),F(8), & H,HOLD,WT(8),PHI(8,14),PSIOLD(12), & PSI(12),BETA(12),SIGMA(13),G(13), & ERKM1,ERK,K,KOLD, & KNEW,KS,IPHASE,IFAIL DIMENSION GSTR(13) DATA GSTR/0.500D0 ,0.0833D0 ,0.0417D0 ,0.0264D0 ,0.0188D0 , & 0.0143D0 ,0.0114D0 ,0.00936D0,0.00789D0,0.00679D0, & 0.00592D0,0.00524D0,0.00468D0/ KNEW = K *** ERK, ERKM1, ERKM2 ERK = 0.0D0 ERKM1 = 0.0D0 ERKM2 = 0.0D0 do L=1,8 TEMP1 = 1.0D0/WT(L) TEMP2 = F(L)-PHI(L,1) ERK = ERK+(TEMP2 *TEMP1)**2 IF (K.EQ.2) THEN ERKM1 = ERKM1+((PHI(L,K )+TEMP2)*TEMP1)**2 ELSE IF (K.GT.2) THEN ERKM1 = ERKM1+((PHI(L,K )+TEMP2)*TEMP1)**2 ERKM2 = ERKM2+((PHI(L,K-1)+TEMP2)*TEMP1)**2 ENDIF enddo TEMP = DABS(H)*DSQRT(ERK) ERK = TEMP*SIGMA(K+1)*GSTR(K) IF (k.EQ.2) THEN ERKM1 = DABS(H)*SIGMA(K )*GSTR(K-1)*DSQRT(ERKM1) IF(ERKM1.LE.0.5D0*ERK) KNEW = K-1 ELSE IF (K.GT.2) THEN ERKM1 = DABS(H)*SIGMA(K )*GSTR(K-1)*DSQRT(ERKM1) ERKM2 = DABS(H)*SIGMA(K-1)*GSTR(K-2)*DSQRT(ERKM2) IF(DMAX1(ERKM1,ERKM2).LE.ERK) KNEW = K-1 ENDIF *** Validity check ERR = TEMP*(G(K)-G(K+1)) IF(ERR.LE.ERR0) THEN CALL ADAM_UPDATE(ERR0) ELSE CALL ADAM_RETRY(ERR0) ENDIF RETURN END ************************************************************************ SUBROUTINE ADAM_UPDATE(ERR0) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) COMMON /QADAM /Y(8),P(8),F(8), & H,HOLD,WT(8),PHI(8,14),PSIOLD(12), & PSI(12),BETA(12),SIGMA(13),G(13), & ERKM1,ERK,K,KOLD, & KNEW,KS,IPHASE,IFAIL DIMENSION GSTR(13),TWO(13) DATA GSTR/0.500D0 ,0.0833D0 ,0.0417D0 ,0.0264D0 ,0.0188D0 , & 0.0143D0 ,0.0114D0 ,0.00936D0,0.00789D0,0.00679D0, & 0.00592D0,0.00524D0,0.00468D0/ DATA TWO/ 2.0D0, 4.0D0, 8.0D0, 16.0D0, 32.0D0, & 64.0D0, 128.0D0, 256.0D0, 512.0D0,1024.0D0, & 2048.0D0,4096.0D0,8192.0D0/ do L=1,8 Y(L) = P(L) + H * G(K+1) * (F(L)-PHI(L,1)) enddo IFAIL = 0 KOLD = K HOLD = H *** Derivative CALL FUNCT(Y,F,1) do L=1,8 PHI(L,K+1) = F(L)-PHI(L,1) PHI(L,K+2) = PHI(L,K+1)-PHI(L,K+2) enddo do I=1,K do L=1,8 PHI(L,I) = PHI(L,I)+PHI(L,K+1) enddo enddo ERKP1 = 0.0D0 IF(KNEW.EQ.K-1.OR.K.EQ.12) IPHASE=0 IF(IPHASE.EQ.1) THEN K = K+1 ERK = ERKP1 GOTO 500 ENDIF IF(KNEW.EQ.K-1) THEN K = K-1 ERK = ERKM1 GOTO 500 ENDIF IF(K+1.LE.KS) THEN do L=1,8 ERKP1 = ERKP1+(PHI(L,K+2)/WT(L))**2 enddo ERKP1 = DABS(H)*GSTR(K+1)*DSQRT(ERKP1) IF(K.LE.1) THEN IF(ERKP1.LT.0.5D0*ERK) THEN K = K+1 ERK = ERKP1 ENDIF GOTO 500 ENDIF IF(ERKM1.LE.DMIN1(ERK,ERKP1)) THEN K = K-1 ERK = ERKM1 ELSE IF(ERKP1.LT.ERK .AND. K.NE.12) THEN K = K+1 ERK = ERKP1 ENDIF ENDIF 500 CONTINUE *** H ->> IF (ERK*TWO(K+1).LT.0.5D0*ERR0 .OR. IPHASE.EQ.1) THEN HNEW = 1.01D0*H GOTO 1000 ENDIF *** <<- H IF (ERK.GT.0.5D0*ERR0) THEN R = (0.5D0*ERR0/ERK)**(1.0D0/DBLE(K+1)) HNEW = DMIN1(0.9D0,DMAX1(0.5D0,R)) * H GOTO 1000 ENDIF *** -> H <- HNEW = H 1000 CONTINUE H = HNEW WRITE(25,800) Y(1)-6370.0, HNEW, ERK 800 FORMAT(F18.2, F18.6, F18.8) RETURN END ************************************************************************ SUBROUTINE ADAM_RETRY(ERR0) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) COMMON /QADAM /Y(8),P(8),F(8), & H,HOLD,WT(8),PHI(8,14),PSIOLD(12), & PSI(12),BETA(12),SIGMA(13),G(13), & ERKM1,ERK,K,KOLD, & KNEW,KS,IPHASE,IFAIL IPHASE = 0 do I=1,K do L=1,8 PHI(L,I) = (PHI(L,I)-PHI(L,I+1))/BETA(I) enddo enddo do I = 2,K PSI(I-1) = PSI(I)-H enddo IFAIL = IFAIL+1 * Euler method IF (IFAIL.GE.3) KNEW = 1 K = KNEW IF (IFAIL.GT.3 .AND. 0.5D0*ERR0.LT.0.25D0*ERK) THEN H = DSQRT(0.5D0*ERR0/ERK) * H ELSE H = 0.5D0 * H ENDIF *** Derivative CALL FUNCT(Y,F,1) RETURN END ************************************************************************ SUBROUTINE OUT0 ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QINIT /R0,T0,P0,E0,D0,F0 COMMON /QOUT /XX,YY,ZZ,OO,PP COMMON /QTHPH0/COST0,SINT0,COSP0,SINP0 XX = R0*dSIN(T0*RAD)*dCOS(P0*RAD) YY = R0*dSIN(T0*RAD)*dSIN(P0*RAD) ZZ = R0*dCOS(T0*RAD) OO = 0.0D0 PP = 0.0D0 RETURN END ************************************************************************ SUBROUTINE OUT1(Y) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) DIMENSION Y(8) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QPATH /FKC,DEL,EPS,GMU,ALF,RAY COMMON /QCOND /L,LL,NN,NLOOP,A COMMON /QOUT /XX,YY,ZZ,OO,PP COMMON /QTHPH0/COST0,SINT0,COSP0,SINP0 COMMON /QREF /PMU,PMUS,PMU0, & PSI , & COPSI ,SIPSI , & COPSIS,SIPSIS, & PX(NSPEC),PY(NSPEC),PZ(NSPEC), & PA1,PA2,PA3, & PK1,PK2,PK4,PK5, & PA ,PB ,PC ,PD COMMON /QDENS /XNS (NSPEC), & DNSDRR(NSPEC), & DNSDTH(NSPEC), & DNSDPH(NSPEC) real*8 RR,TH,PH,DL RR = Y(1) TH = Y(2) PH = Y(3) CALL GTOM1(TH,PH) DL = DEL XXX = RR*dSIN(TH*RAD)*dCOS(PH*RAD) YYY = RR*dSIN(TH*RAD)*dSIN(PH*RAD) ZZZ = RR*dCOS(TH*RAD) DDD = dSQRT((XXX-XX)**2+(YYY-YY)**2+(ZZZ-ZZ)**2) OO = OO+DDD XX = XXX YY = YYY ZZ = ZZZ IF(L.NE.0 .OR. LL.NE.0) GOTO 1020 IF(OO.GE.PP) THEN write(21,1015) & REAL(90.-TH), & REAL(RR-6370.), & REAL(PH), & REAL(DL), & PSI, & 8.98*dsqrt(XNS(1)), & REAL(PMU), & REAL(Y(7)) c & Above time(ms)is the propagation delay c & ((dcos((90.0-TH)*RAD)**2+0.5)/1.5)* c & (XNS(1)+XNS(2)+XNS(3)+XNS(4)) 1015 FORMAT(F8.2,F9.2,F9.5,F7.1,F10.1,F12.3,F11.2,F12.4) PP=PP+DISOUT A=A+1 END IF 1020 CONTINUE RETURN END ************************************************************************ SUBROUTINE FUNCT(PATH,DPATH,IMDFY) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) real*8 path(8),dpath(8) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QPATH /FKC,DEL,EPS,GMU,ALF,RAY COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QCOND /L,LL,NN,NLOOP,a COMMON /QREF /PMU,PMUS,PMU0, & PSI , & COPSI ,SIPSI , & COPSIS,SIPSIS, & PX(NSPEC),PY(NSPEC),PZ(NSPEC), & PA1,PA2,PA3, & PK1,PK2,PK4,PK5, & PA ,PB ,PC ,PD COMMON /QIGRF /YR ,YT ,YP , & DYRDRR,DYTDRR,DYPDRR, & DYRDTH,DYTDTH,DYPDTH, & DYRDPH,DYTDPH,DYPDPH, & BB , & DBBDRR,DBBDTH,DBBDPH COMMON /QDENS /XNS (NSPEC), & DNSDRR(NSPEC), & DNSDTH(NSPEC), & DNSDPH(NSPEC) real*8 DMUDRR,DMUDTH,DMUDPH, & DK2DF ,DA3DF ,DENOM CALL REF(PATH) IF(IMDFY.EQ.1) THEN PATH(4) = PATH(4)*PMU/PMU0 PATH(5) = PATH(5)*PMU/PMU0 PATH(6) = PATH(6)*PMU/PMU0 DEL = dATAN2(PATH(5),PATH(4)) EPS = dATAN (PATH(6)/dSQRT(PATH(4)**2+PATH(5)**2)) DEL = DEL*DEG EPS = EPS*DEG END IF DMUDRR = 0.0D0 DMUDTH = 0.0D0 DMUDPH = 0.0D0 DK2DF = 0.0D0 DA3DF = 0.0D0 DENOM = (4.0D0*PA*PMUS+2.0D0*PB)*PMU DO 20 I=1,NSPEC DADX = - COPSIS DBDX = -(PK1/PZ(I)-PK2)*(1.0D0+COPSIS) DCDX = - PK5 PK2A3Y = 2.0D0*(PK2-PA3*PY(I)) DADX = DADX+ SIPSIS/PZ(I) DBDX = DBDX- PK2A3Y*SIPSIS/PZ(I) DCDX = DCDX+PK1*PK2A3Y /PZ(I) DMUDX = -(DADX*PMUS**2+DBDX*PMUS+DCDX)/DENOM DK2DY = -2.0D0 *PX(I)* PY(I) /PZ(I)**2 DK5DY = 2.0D0*PA3*PX(I)*(PZ(I)+2.0D0)/PZ(I)**2 DK5DY = DK5DY+2.0D0*DK2DY*PK2 DADY = DK2DY*SIPSIS DBDY = -PK1*DK2DY*(1.0D0+COPSIS)-DK5DY*SIPSIS DCDY = PK1*DK5DY DMUDY = -(DADY*PMUS**2+DBDY*PMUS+DCDY)/DENOM DXDRR = PX(I)*DNSDRR(I) DXDTH = PX(I)*DNSDTH(I) DXDPH = PX(I)*DNSDPH(I) DYDRR = PY(I)*DBBDRR DYDTH = PY(I)*DBBDTH DYDPH = PY(I)*DBBDPH DMUDRR = DMUDRR+DMUDX*DXDRR+DMUDY*DYDRR DMUDTH = DMUDTH+DMUDX*DXDTH+DMUDY*DYDTH DMUDPH = DMUDPH+DMUDX*DXDPH+DMUDY*DYDPH DK2DFI = 2.0D0 *PX(I) /PZ(I)**2 DA3DFI = (2.0D0-PZ(I))*PX(I)*PY(I)/PZ(I)**2 DK2DF = DK2DF+DK2DFI/FKC DA3DF = DA3DF+DA3DFI/FKC 20 CONTINUE DADPS = PA1+PA2 DBDPS = PK4-PK5 DMUDPS = (DADPS*PMUS+DBDPS)*PMUS DMUDPS = -2.0D0*COPSI*DMUDPS/DENOM DPSDRR = DYRDRR*PATH(4) & +DYTDRR*PATH(5) & +DYPDRR*PATH(6) DPSDTH = DYRDTH*PATH(4) & +DYTDTH*PATH(5) & +DYPDTH*PATH(6) DPSDPH = DYRDPH*PATH(4) & +DYTDPH*PATH(5) & +DYPDPH*PATH(6) DPSDRR = -DPSDRR/PMU DPSDTH = -DPSDTH/PMU DPSDPH = -DPSDPH/PMU if (XNS(1).NE.0.0d0) then DMUDRR = (DMUDRR+DMUDPS*DPSDRR)/PMU DMUDTH = (DMUDTH+DMUDPS*DPSDTH)/PMU DMUDPH = (DMUDPH+DMUDPS*DPSDPH)/PMU DMUDVR = DMUDPS*(PATH(4)*COPSI/PMU-YR) DMUDVT = DMUDPS*(PATH(5)*COPSI/PMU-YT) DMUDVP = DMUDPS*(PATH(6)*COPSI/PMU-YP) else DMUDRR = 0.0d0 DMUDTH = 0.0d0 DMUDPH = 0.0d0 DMUDVR = 0.0d0 DMUDVT = 0.0d0 DMUDVP = 0.0d0 endif DK1DF = 2.0D0* PA1/FKC DK4DF = PK1*DK2DF+PK2*DK1DF DK5DF = 2.0D0*(PK2*DK2DF-PA3*DA3DF) DCDF = PK1*DK5DF +PK5*DK1DF DADF = DK1DF* COPSIS +DK2DF*SIPSIS DBDF = -DK4DF*(1.0D0+COPSIS)-DK5DF *SIPSIS DMUDF = -((DADF*PMUS+DBDF)*PMUS+DCDF)/DENOM if (XNS(1).EQ.0.0d0) dmudf=0.0d0 GMU = PMU+FKC*DMUDF IF (PATH(1).LT.6460.0) THEN DMUDVR = 0.0d0 DMUDVR = 0.0d0 DMUDVP = 0.0d0 DMUDRR = 0.0d0 DMUDTH = 0.0d0 DMUDPH = 0.0d0 END IF DPATH(1) = PATH(4)-DMUDVR DPATH(2) = PATH(5)-DMUDVT DPATH(3) = PATH(6)-DMUDVP DPATH(1) = DPATH(1)/PMUS DPATH(2) = DPATH(2)/PMUS/PATH(1) DPATH(3) = DPATH(3)/PMUS/PATH(1)/dSIN(PATH(2)*RAD) DPATH(4) = DMUDRR & +PATH(5)*DPATH(2) & +PATH(6)*DPATH(3)*dSIN (PATH(2)*RAD) DPATH(5) = DMUDTH/PATH(1) & -PATH(5)*DPATH(1)/PATH(1) & +PATH(6)*DPATH(3)*dCOS (PATH(2)*RAD) DPATH(6) = DMUDPH/PATH(1)/dSIN(PATH(2)*RAD) & -PATH(6)*DPATH(1)/PATH(1) & -PATH(6)*DPATH(2)/dTAN(PATH(2)*RAD) DPATH(2) = DPATH(2)*DEG DPATH(3) = DPATH(3)*DEG DPATH(7) = 1.0D0/PMU/CVLCTY DPATH(8) = 1.0D0/CVLCTY if (XNS(1).eq.0.0d0) alf=0.0d0 ALF = -ATAN(DMUDPS*SIPSI/PMU)*DEG RAY = PSI+ALF RETURN END ************************************************************************ FUNCTION RIPARA(PSQ) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QREF /PMU,PMUS,PMU0, & PSI , & COPSI ,SIPSI , & COPSIS,SIPSIS, & PX(NSPEC),PY(NSPEC),PZ(NSPEC), & PA1,PA2,PA3, & PK1,PK2,PK4,PK5, & PA ,PB ,PC ,PD PAQ = PK1* dCOS(PSQ*RAD)**2 & +PK2* dSIN(PSQ*RAD)**2 PBQ =-PK4*(dCOS(PSQ*RAD)**2+1.0D0) & -PK5* dSIN(PSQ*RAD)**2 PCQ = PK1*PK5 PDQ = PBQ**2-4.0D0*PAQ*PCQ IF(PDQ.GE.0.0D0) THEN IF(PBQ.GE.0.0D0) THEN PMUSQ = (-PBQ-dSQRT(PDQ))/(2.0D0*PAQ) ELSE PMUSQ = (2.0D0*PCQ)/(-PBQ+dSQRT(PDQ)) END IF IF(PMUSQ.GT.0.0D0) THEN RIPARA = dSQRT(PMUSQ)*dCOS(PSQ*RAD) ELSE RIPARA = 1.0D50 END IF ELSE RIPARA = 1.0D50 END IF RETURN END ************************************************************************ SUBROUTINE REF(PATH) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) real*8 path(8) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QPATH /FKC,DEL,EPS,GMU,ALF,RAY COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QCOND /L,LL,NN,NLOOP,a COMMON /QREF /PMU,PMUS,PMU0, & PSI , & COPSI ,SIPSI , & COPSIS,SIPSIS, & PX(NSPEC),PY(NSPEC),PZ(NSPEC), & PA1,PA2,PA3, & PK1,PK2,PK4,PK5, & PA ,PB ,PC ,PD COMMON /QIGRF /YR ,YT ,YP , & DYRDRR,DYTDRR,DYPDRR, & DYRDTH,DYTDTH,DYPDTH, & DYRDPH,DYTDPH,DYPDPH, & FH , & DBBDRR,DBBDTH,DBBDPH COMMON /QDENS /XNS (NSPEC), & DNSDRR(NSPEC), & DNSDTH(NSPEC), & DNSDPH(NSPEC) CALL IGRF(PATH) c CALL DIPOLE(PATH) CALL DENS(PATH) DO 10 I=1,NSPEC PX(I) = PF0(I)*XNS(I)/FKC**2 PY(I) = GF0(I)*FH /FKC PZ(I) = PY(I)**2-1.0D0 10 CONTINUE PMU0 = dSQRT(PATH(4)**2+PATH(5)**2+PATH(6)**2) COPSI = (YR*PATH(4) & +YT*PATH(5) & +YP*PATH(6))/PMU0 COPSIS = COPSI**2 SIPSIS = dABS(1.0D0-COPSIS) SIPSI = dSQRT(SIPSIS) if (dabs(copsi).gt.1.0) copsi=1.0d0 PSI = dACOS(COPSI)*DEG PA1 = 0.0D0 PA2 = 0.0D0 PA3 = 0.0D0 DO 20 I=1,NSPEC PA1 = PA1+PX(I) PA2 = PA2+PX(I) /PZ(I) PA3 = PA3+PX(I)*PY(I)/PZ(I) 20 CONTINUE PK1 = 1.0D0-PA1 PK2 = 1.0D0+PA2 PK4 = PK1*PK2 PK5 = PK2**2-PA3**2 PA = PK1*COPSIS+PK2*SIPSIS PB =-PK4*(COPSIS+1.0D0)-PK5*SIPSIS PC = PK1*PK5 PD = DABS(PB**2-4.0D0*PA*PC) if (pd.lt.0.0d0) go to 1040 PMUS= (-PB-SQRT(PD))/2.0/PA c Above is Extraordinary mode(whistler mode) if (pmus.lt.0.0d0) go to 1050 PMU = dSQRT(PMUS) 1040 CONTINUE IF (PATH(1).LT.6370.0d0+90.0d0) THEN PMU=1.0d0 gmu=1.0d0 END IF 1050 CONTINUE RETURN END ************************************************************************ SUBROUTINE DENS(PATH) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) real*8 path(8) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QPATH /FKC,DEL,EPS,GMU,ALF,RAY COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QDENS /XNS (NSPEC), & DNSDRR(NSPEC), & DNSDTH(NSPEC), & DNSDPH(NSPEC) COMMON /QDNSDE/XNDE (NSPEC), & DNDDRR(NSPEC) COMMON /QDNSLI/XNLI ,DNLDRR COMMON /QDNSAO/XNAO ,DNADRR,DNADTH,DNADPH COMMON /QDNSOP/XNOP ,DNODRR,DNODTH,DNODPH call FLPRM(PATH) CALL DENSDE(PATH) CALL DENSLI(PATH) CALL DENSAO(PATH) CALL DENSOP(PATH) DO 10 IS=1,NSPEC XNS(IS) = XNDE(IS)*XNLI*XNAO*XNOP *... 1/N * dn/dr DNSDRR(IS) = DNDDRR(IS)+DNLDRR+DNADRR+DNODRR DNSDTH(IS) = DNADTH+DNODTH DNSDPH(IS) = DNADPH+DNODPH 10 CONTINUE RETURN END ************************************************************************ SUBROUTINE DENSAO(PATH) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) real*8 path(8) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QIGRF /YR ,YT ,YP , & DYRDRR,DYTDRR,DYPDRR, & DYRDTH,DYTDTH,DYPDTH, & DYRDPH,DYTDPH,DYPDPH, & BB , & DBBDRR,DBBDTH,DBBDPH COMMON /QDNSAO/XNAO ,DNADRR,DNADTH,DNADPH COMMON /QFLPRM/TH0 , & DT0DRR,DT0DTH,DT0DPH, & PH0 , & DP0DRR,DP0DTH,DP0DPH, & BB0 , & DB0DRR,DB0DTH,DB0DPH, & XLV , & DLVDRR,DLVDTH,DLVDPH XNAO = 1.0D0 DNADRR = 0.0D0 DNADTH = 0.0D0 DNADPH = 0.0D0 RETURN END ************************************************************************ SUBROUTINE DENSDE(PATH) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) real*8 path(8) COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QDNSDE/XNDE (NSPEC), & DNDDRR(NSPEC) GPH = RREF*(PATH(1)-RREF)/PATH(1) DGPHDR = (RREF/PATH(1))**2 ETEXP2 = ETA(2)*EXP(-GPH*SHINV(2)) ETEXP3 = ETA(3)*EXP(-GPH*SHINV(3)) ETEXP4 = ETA(4)*EXP(-GPH*SHINV(4)) SUM1 = ETEXP2 & + ETEXP3 & + ETEXP4 SUM2 = ETEXP2*SHINV(2) & + ETEXP3*SHINV(3) & + ETEXP4*SHINV(4) XNDE (1) = dSQRT(SUM1) XNDE (2) = ETEXP2/XNDE(1) XNDE (3) = ETEXP3/XNDE(1) XNDE (4) = ETEXP4/XNDE(1) DNDDRR(1) =-DGPHDR*SUM2/SUM1/2.0D0 DNDDRR(2) =-DGPHDR*SHINV(2)-DNDDRR(1) DNDDRR(3) =-DGPHDR*SHINV(3)-DNDDRR(1) DNDDRR(4) =-DGPHDR*SHINV(4)-DNDDRR(1) RETURN END ************************************************************************ SUBROUTINE DENSLI(PATH) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) real*8 path(8) COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QDNSLI/XNLI, DNLDRR IF ((PATH(1)-RRLI).GE.0) THEN ALTNRM = (PATH(1)-RRLI)/HWLI XNLI = 1.0D0-dEXP(-ALTNRM**2) DNLDRR = dEXP(-ALTNRM**2)*ALTNRM*2.0D0/HWLI/XNLI ELSE XNLI=0.0 DNLDRR=0.0 END IF RETURN END ************************************************************************ SUBROUTINE DENSOP(PATH) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) real*8 path(8) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QPATH /FKC,DEL,EPS,GMU,ALF,RAY COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QDNSOP/XNOP ,DNODRR,DNODTH,DNODPH XNOP = XN0 dNODRR = 0.0D0 DNODTH = 0.0d0 dnodph = 0.0d0 RETURN END ********************************************************************** SUBROUTINE DIPOLE(PATH) ********************************************************************** IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) DIMENSION PATH(7) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QIGRF /YR ,YT ,YP , & DYRDRR,DYTDRR,DYPDRR, & DYRDTH,DYTDTH,DYPDTH, & DYRDPH,DYTDPH,DYPDPH, & BB , & DBBDRR,DBBDTH,DBBDPH DATA GMDM/8.07D6/ RR = PATH(1) TH = PATH(2) PH = PATH(3) COSTH = COS(TH*RAD) SINTH = SIN(TH*RAD) COSTHS = 1.0D0+3.0D0*COSTH**2 YR =-2.0D0*COSTH/SQRT(COSTHS) YT =- SINTH/SQRT(COSTHS) YP = 0.0D0 DYRDRR = 0.0D0 DYTDRR = 0.0D0 DYPDRR = 0.0D0 DYRDTH =-2.0D0/COSTHS*YT DYTDTH = 2.0D0/COSTHS*YR DYPDTH = 0.0D0 DYRDPH = 0.0D0 DYTDPH = 0.0D0 DYPDPH = 0.0D0 BB = GMDM*SQRT(COSTHS)/RR**3 DBBDRR =-3.0D0/RR DBBDTH =-3.0D0*SINTH*COSTH/COSTHS DBBDPH = 0.0D0 RETURN END ********************************************************************** SUBROUTINE IGRF(PATH) ********************************************************************** IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) DIMENSION PATH(7) real*8 yr,yt,yp,dyrdrr,dytdrr,dypdrr,dyrdth,dytdth,dypdth, & dyrdph,dytdph,dypdph,bb,dbbdrr,dbbdth,dbbdph common /QIGRF/yr ,yt ,yp , & dyrdrr,dytdrr,dypdrr, & dyrdth,dytdth,dypdth, & dyrdph,dytdph,dypdph, & bb , & dbbdrr,dbbdth,dbbdph COMMON /QIGCOE/GG(1:8,0:8),HH(1:8,0:8),SS(0:8,0:8) PARAMETER (NIGRF=8) REAL*8 AARR(3:NIGRF+3), & COSTH,SINTH,COTTH, & COSMPH(0:NIGRF),SINMPH(0:NIGRF), & P (0:NIGRF,0:NIGRF), & DP (0:NIGRF,0:NIGRF), & DDP(0:NIGRF,0:NIGRF), & BR ,BT ,BP , & DBRDRR,DBTDRR,DBPDRR, & DBRDTH,DBTDTH,DBPDTH, & DBRDPH,DBTDPH,DBPDPH, & XX ,YY ,ZZ , & DXX,DYX,DZX, & DXY,DYY,DZY, & DXZ,DYZ,DZZ DATA AA/6371.2D0/ RAD = 3.1415926/180.0 P (0,0) = 1.0D0 DP (0,0) = 0.0D0 DDP(0,0) = 0.0D0 BR = 0.0D0 BT = 0.0D0 BP = 0.0D0 DBRDRR = 0.0D0 DBTDRR = 0.0D0 DBPDRR = 0.0D0 DBRDTH = 0.0D0 DBTDTH = 0.0D0 DBPDTH = 0.0D0 DBRDPH = 0.0D0 DBTDPH = 0.0D0 DBPDPH = 0.0D0 AARR (3) = (AA/PATH(1))**3 COSTH = COS(PATH(2)*RAD) SINTH = SIN(PATH(2)*RAD) COTTH = COSTH/SINTH DO 10 N=4,NIGRF+3 AARR(N) = (AA/PATH(1))*AARR(N-1) 10 CONTINUE DO 20 M=0,NIGRF COSMPH(M) = COS(M*PATH(3)*RAD) SINMPH(M) = SIN(M*PATH(3)*RAD) 20 CONTINUE DO 30 N=1,NIGRF XX = 0.0D0 YY = 0.0D0 ZZ = 0.0D0 DXX = 0.0D0 DYX = 0.0D0 DZX = 0.0D0 DXY = 0.0D0 DYY = 0.0D0 DZY = 0.0D0 DXZ = 0.0D0 DYZ = 0.0D0 DZZ = 0.0D0 P (N,N ) = SINTH*P (N-1,N-1) DP (N,N ) = SINTH*DP (N-1,N-1) & +COSTH*P (N-1,N-1) DDP(N,N ) = SINTH*DDP(N-1,N-1) & +COSTH*DP (N-1,N-1)*2.0D0 & -SINTH*P (N-1,N-1) P (N,N-1) = COSTH*P (N-1,N-1) DP (N,N-1) = COSTH*DP (N-1,N-1) & -SINTH*P (N-1,N-1) DDP(N,N-1) = COSTH*DDP(N-1,N-1) & -SINTH*DP (N-1,N-1)*2.0D0 & -COSTH*P (N-1,N-1) DO 40 M=N-2,0,-1 REALK = DBLE((N-1)**2-M**2)/DBLE((2*N-3)*(2*N-1)) P (N,M) = COSTH*P (N-1,M) & -REALK *P (N-2,M) DP (N,M) = COSTH*DP (N-1,M) & -SINTH*P (N-1,M) & -REALK *DP (N-2,M) DDP(N,M) = COSTH*DDP(N-1,M) & -SINTH*DP (N-1,M)*2.0D0 & -COSTH*P (N-1,M) & -REALK *DDP(N-2,M) 40 CONTINUE DO 50 M=0,N PP = SS(N,M)*P (N,M) DPP = SS(N,M)*DP (N,M) DDPP = SS(N,M)*DDP(N,M) GCHS = GG(N,M)*COSMPH(M) + HH(N,M)*SINMPH(M) GSHC = GG(N,M)*SINMPH(M) - HH(N,M)*COSMPH(M) XX = XX + GCHS*PP YY = YY + GCHS*DPP ZZ = ZZ + M* GSHC*PP DXX = DXX + GCHS*PP DYX = DYX + GCHS*DPP DZX = DZX + M* GSHC*PP DXY = DXY + GCHS*DPP DYY = DYY + GCHS*DDPP DZY = DZY + M* GSHC*(PP*COTTH-DPP) DXZ = DXZ + M* GSHC*PP DYZ = DYZ + M* GSHC*DPP DZZ = DZZ + M**2*GCHS*PP 50 CONTINUE BR = BR + (N+1)* AARR(N+2)*XX BT = BT + AARR(N+2)*YY BP = BP + AARR(N+2)*ZZ DBRDRR = DBRDRR + (N+1)*(N+2)*AARR(N+3)*DXX DBTDRR = DBTDRR + (N+2)*AARR(N+3)*DYX DBPDRR = DBPDRR + (N+2)*AARR(N+3)*DZX DBRDTH = DBRDTH + (N+1)* AARR(N+2)*DXY DBTDTH = DBTDTH + AARR(N+2)*DYY DBPDTH = DBPDTH + AARR(N+2)*DZY DBRDPH = DBRDPH + (N+1)* AARR(N+2)*DXZ DBTDPH = DBTDPH + AARR(N+2)*DYZ DBPDPH = DBPDPH + AARR(N+2)*DZZ 30 CONTINUE BR = BR BT =-BT BP = BP /SINTH DBRDRR =-DBRDRR/AA DBTDRR = DBTDRR/AA DBPDRR =-DBPDRR/AA/SINTH DBRDTH = DBRDTH DBTDTH =-DBTDTH DBPDTH =-DBPDTH /SINTH DBRDPH =-DBRDPH DBTDPH = DBTDPH DBPDPH = DBPDPH /SINTH BBS = BR**2+BT**2+BP**2 BB = SQRT(BBS) YR = BR/BB YT = BT/BB YP = BP/BB DBBDRR = (BR*DBRDRR & +BT*DBTDRR & +BP*DBPDRR)/BBS DBBDTH = (BR*DBRDTH & +BT*DBTDTH & +BP*DBPDTH)/BBS DBBDPH = (BR*DBRDPH & +BT*DBTDPH & +BP*DBPDPH)/BBS DYRDRR = (DBRDRR-BR*DBBDRR)/BB DYTDRR = (DBTDRR-BT*DBBDRR)/BB DYPDRR = (DBPDRR-BP*DBBDRR)/BB DYRDTH = (DBRDTH-BR*DBBDTH)/BB DYTDTH = (DBTDTH-BT*DBBDTH)/BB DYPDTH = (DBPDTH-BP*DBBDTH)/BB DYRDPH = (DBRDPH-BR*DBBDPH)/BB DYTDPH = (DBTDPH-BT*DBBDPH)/BB DYPDPH = (DBPDPH-BP*DBBDPH)/BB RETURN END ********************************************************************** SUBROUTINE IGRF90(IYEAR,IMNTH,IDATE) ********************************************************************** IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) COMMON /QIGCOE/GG(1:8,0:8),HH(1:8,0:8),SS(0:8,0:8) DIMENSION G(7,1:8,0:8),H(7,1:8,0:8) IF((IYEAR.LT.1965).OR.(IYEAR.GT.1995)) GOTO 999 DATE=IYEAR+(30.4176D0*(IMNTH-1)+(IDATE-1))/365.0D0 IF(IYEAR.GE.1990) THEN L=6 P=1.0D0 Q=DATE-1990.0D0 ELSE IF(IYEAR.GE.1985) THEN L=5 Q=(DATE-1985.0D0)/5.0D0 P=1.0D0-Q ELSE IF(IYEAR.GE.1980) THEN L=4 Q=(DATE-1980.0D0)/5.0D0 P=1.0D0-Q ELSE IF(IYEAR.GE.1975) THEN L=3 Q=(DATE-1975.0D0)/5.0D0 P=1.0D0-Q ELSE IF(IYEAR.GE.1970) THEN L=2 Q=(DATE-1970.0D0)/5.0D0 P=1.0D0-Q ELSE IF(IYEAR.GE.1965) THEN L=1 Q=(DATE-1965.0D0)/5.0D0 P=1.0D0-Q END IF C C THE SPHERICAL HARMONIC COEFFICIENTS IN UNITS OF TESLA (T) C DO 20 N=1,8 DO 30 M=0,N GG(N,M)=(G(L,N,M)*P+G(L+1,N,M)*Q)*1.0D-9 HH(N,M)=(H(L,N,M)*P+H(L+1,N,M)*Q)*1.0D-9 30 CONTINUE 20 CONTINUE C C THE SCHMIDT COEFFICIENT C SS(0,0)=1.0D0 DO 40 N=1,8 SS(N,0)=SS(N-1,0)*(2*N-1)/DBLE(N) SS(N,1)=SS(N,0)*SQRT(2*N/DBLE(N+1)) DO 50 M=2,N SS(N,M)=SS(N,M-1)*SQRT((N-M+1)/DBLE(N+M)) 50 CONTINUE 40 CONTINUE 999 RETURN C C THE SPHERICAL HARMONIC COEFFICIENTS OF IGRF ( 1991 REVISION ) C C *** G65 & H65 *** DATA (G(1,1,M),M=0,1)/-30334.,-2119./ & (G(1,2,M),M=0,2)/-1662.,2997.,1594./ & (G(1,3,M),M=0,3)/1297.,-2038.,1292.,856./ & (G(1,4,M),M=0,4)/957.,804.,479.,-390.,252./ & (G(1,5,M),M=0,5)/-219.,358.,254.,-31.,-157.,-62./ & (G(1,6,M),M=0,6)/45.,61.,8.,-228.,4.,1.,-111./ & (G(1,7,M),M=0,7)/75.,-57.,4.,13.,-26.,-6.,13.,1./ & (G(1,8,M),M=0,8)/13.,5.,-4.,-14.,0.,8.,-1.,11.,4./ DATA (H(1,1,M),M=0,1)/0.,5776./ & (H(1,2,M),M=0,2)/0.,-2016.,114./ & (H(1,3,M),M=0,3)/0.,-404.,240.,-165./ & (H(1,4,M),M=0,4)/0.,148.,-269.,13.,-269./ & (H(1,5,M),M=0,5)/0.,19.,128.,-126.,-97.,81./ & (H(1,6,M),M=0,6)/0.,-11.,100.,68.,-32.,-8.,-7./ & (H(1,7,M),M=0,7)/0.,-61.,-27.,-2.,6.,26.,-23.,-12./ & (H(1,8,M),M=0,8)/0.,7.,-12.,9.,-16.,4.,24.,-3.,-17./ C *** G70 & H70 *** DATA (G(2,1,M),M=0,1)/-30220.,-2068./ & (G(2,2,M),M=0,2)/-1781.,3000.,1611./ & (G(2,3,M),M=0,3)/1287.,-2091.,1278.,838./ & (G(2,4,M),M=0,4)/952.,800.,461.,-395.,234./ & (G(2,5,M),M=0,5)/-216.,359.,262.,-42.,-160.,-56./ & (G(2,6,M),M=0,6)/43.,64.,15.,-212.,2.,3.,-112./ & (G(2,7,M),M=0,7)/72.,-57.,1.,14.,-22.,-2.,13.,-2./ & (G(2,8,M),M=0,8)/14.,6.,-2.,-13.,-3.,5.,0.,11.,3./ DATA (H(2,1,M),M=0,1)/0.,5737./ & (H(2,2,M),M=0,2)/0.,-2047.,25./ & (H(2,3,M),M=0,3)/0.,-366.,251.,-196./ & (H(2,4,M),M=0,4)/0.,167.,-266.,26.,-279./ & (H(2,5,M),M=0,5)/0.,26.,139.,-139.,-91.,83./ & (H(2,6,M),M=0,6)/0.,-12.,100.,72.,-37.,-6.,1./ & (H(2,7,M),M=0,7)/0.,-70.,-27.,-4.,8.,23.,-23.,-11./ & (H(2,8,M),M=0,8)/0.,7.,-15.,6.,-17.,6.,21.,-6.,-16./ C *** G75 & H75 *** DATA (G(3,1,M),M=0,1)/-30100.,-2013./ & (G(3,2,M),M=0,2)/-1902.,3010.,1632./ & (G(3,3,M),M=0,3)/1276.,-2144.,1260.,830./ & (G(3,4,M),M=0,4)/946.,791.,438.,-405.,216./ & (G(3,5,M),M=0,5)/-218.,356.,264.,-59.,-159.,-49./ & (G(3,6,M),M=0,6)/45.,66.,28.,-198.,1.,6.,-111./ & (G(3,7,M),M=0,7)/71.,-56.,1.,16.,-14.,0.,12.,-5./ & (G(3,8,M),M=0,8)/14.,6.,-1.,-12.,-8.,4.,0.,10.,1./ DATA (H(3,1,M),M=0,1)/0.,5675./ & (H(3,2,M),M=0,2)/0.,-2067.,68./ & (H(3,3,M),M=0,3)/0.,-333.,262.,-223./ & (H(3,4,M),M=0,4)/0.,191.,-265.,39.,-288./ & (H(3,5,M),M=0,5)/0.,31.,148.,-152.,-83.,88./ & (H(3,6,M),M=0,6)/0.,-13.,99.,75.,-41.,-4.,11./ & (H(3,7,M),M=0,7)/0.,-77.,-26.,-5.,10.,22.,-23.,-12./ & (H(3,8,M),M=0,8)/0.,6.,-16.,4.,-19.,6.,18.,-10.,-17./ C *** G80 & H80 *** DATA (G(4,1,M),M=0,1)/-29992.,-1956./ & (G(4,2,M),M=0,2)/-1997.,3027.,1663./ & (G(4,3,M),M=0,3)/1281.,-2180.,1251.,833./ & (G(4,4,M),M=0,4)/938.,782.,398.,-419.,199./ & (G(4,5,M),M=0,5)/-218.,357.,261.,-74.,-162.,-48./ & (G(4,6,M),M=0,6)/48.,66.,42.,-192.,4.,14.,-108./ & (G(4,7,M),M=0,7)/72.,-59.,2.,21.,-12.,1.,11.,-2./ & (G(4,8,M),M=0,8)/18.,6.,0.,-11.,-7.,4.,3.,6.,-1./ DATA (H(4,1,M),M=0,1)/0.,5604./ & (H(4,2,M),M=0,2)/0.,-2129.,-200./ & (H(4,3,M),M=0,3)/0.,-336.,271.,-252./ & (H(4,4,M),M=0,4)/0.,212.,-257.,53.,-297./ & (H(4,5,M),M=0,5)/0.,46.,150.,-151.,-78.,92./ & (H(4,6,M),M=0,6)/0.,-15.,93.,71.,-43.,-2.,17./ & (H(4,7,M),M=0,7)/0.,-82.,-27.,-5.,16.,18.,-23.,-10./ & (H(4,8,M),M=0,8)/0.,7.,-18.,4.,-22.,9.,16.,-13.,-15./ C *** G85 & H85 *** DATA (G(5,1,M),M=0,1)/-29873.,-1905./ & (G(5,2,M),M=0,2)/-2072.,3044.,1687./ & (G(5,3,M),M=0,3)/1296.,-2208.,1247.,829./ & (G(5,4,M),M=0,4)/936.,780.,361.,-424.,170./ & (G(5,5,M),M=0,5)/-214.,355.,253.,-93.,-164.,-46./ & (G(5,6,M),M=0,6)/53.,65.,51.,-185.,4.,16.,-102./ & (G(5,7,M),M=0,7)/74.,-62.,3.,24.,-6.,4.,10.,0./ & (G(5,8,M),M=0,8)/21.,6.,0.,-11.,-9.,4.,4.,4.,-4./ DATA (H(5,1,M),M=0,1)/0.,5500./ & (H(5,2,M),M=0,2)/0.,-2197.,-306./ & (H(5,3,M),M=0,3)/0.,-310.,284.,-297./ & (H(5,4,M),M=0,4)/0.,232.,-249.,69.,-297./ & (H(5,5,M),M=0,5)/0.,47.,150.,-154.,-75.,95./ & (H(5,6,M),M=0,6)/0.,-16.,88.,69.,-48.,-1.,21./ & (H(5,7,M),M=0,7)/0.,-83.,-27.,-2.,20.,17.,-23.,-7./ & (H(5,8,M),M=0,8)/0.,8.,-19.,5.,-23.,11.,14.,-15.,-11./ C *** G90 & H90 *** DATA (G(6,1,M),M=0,1)/-29775.,-1851./ & (G(6,2,M),M=0,2)/-2136.,3058.,1693./ & (G(6,3,M),M=0,3)/1315.,-2240.,1246.,807./ & (G(6,4,M),M=0,4)/939.,782.,324.,-423.,142./ & (G(6,5,M),M=0,5)/-211.,353.,244.,-111.,-166.,-37./ & (G(6,6,M),M=0,6)/61.,64.,60.,-178.,2.,17.,-96./ & (G(6,7,M),M=0,7)/77.,-64.,4.,28.,1.,6.,10.,0./ & (G(6,8,M),M=0,8)/22.,5.,-1.,-11.,-12.,4.,4.,3.,-6./ DATA (H(6,1,M),M=0,1)/0.,5411./ & (H(6,2,M),M=0,2)/0.,-2278.,-380./ & (H(6,3,M),M=0,3)/0.,-287.,293.,-348./ & (H(6,4,M),M=0,4)/0.,248.,-240.,87.,-299./ & (H(6,5,M),M=0,5)/0.,47.,153.,-154.,-69.,98./ & (H(6,6,M),M=0,6)/0.,-16.,83.,68.,-52.,2.,27./ & (H(6,7,M),M=0,7)/0.,-81.,-27.,1.,20.,16.,-23.,-5./ & (H(6,8,M),M=0,8)/0.,10.,-20.,7.,-22.,12.,11.,-16.,-11./ C *** DG90 & DH90 *** DATA (G(7,1,M),M=0,1)/18.0,10.6/ & (G(7,2,M),M=0,2)/-12.9,2.4,0./ & (G(7,3,M),M=0,3)/3.3,-6.7,0.1,-5.9/ & (G(7,4,M),M=0,4)/0.5,0.6,-7.,0.5,-5.5/ & (G(7,5,M),M=0,5)/0.6,-0.1,-1.6,-3.1,-0.1,2.3/ & (G(7,6,M),M=0,6)/1.3,-0.2,1.8,1.3,-0.2,0.1,1.2/ & (G(7,7,M),M=0,7)/0.6,-0.5,-0.3,0.6,1.6,0.2,0.2,0.3/ & (G(7,8,M),M=0,8)/0.2,-0.7,-0.2,0.1,-1.1,0.,-0.1,-0.5,-0.6/ DATA (H(7,1,M),M=0,1)/0.,-16.1/ & (H(7,2,M),M=0,2)/0.,-15.8,-13.8/ & (H(7,3,M),M=0,3)/0.,4.4,1.6,-10.6/ & (H(7,4,M),M=0,4)/0.,2.6,1.8,3.1,-1.4/ & (H(7,5,M),M=0,5)/0.,-0.1,0.5,0.4,1.7,0.4/ & (H(7,6,M),M=0,6)/0.,0.2,-1.3,0.,-0.9,0.5,1.2/ & (H(7,7,M),M=0,7)/0.,0.6,0.2,0.8,-0.5,-0.2,0.,0./ & (H(7,8,M),M=0,8)/0.,0.5,-0.2,0.3,0.3,0.4,-0.5,-0.3,0.6/ END ************************************************************************ SUBROUTINE FLPRM(PATH) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) real*8 path(8) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QFLPRM/TH0 , & DT0DRR,DT0DTH,DT0DPH, & PH0 , & DP0DRR,DP0DTH,DP0DPH, & BB0 , & DB0DRR,DB0DTH,DB0DPH, & XLV , & DLVDRR,DLVDTH,DLVDPH REAL*8 RR ,TH ,PH , & DRR ,DTH ,DPH , & RRDRR ,THDTH ,PHDTH , & THDPH ,PHDPH ,TH0DRR, & TH0DTH,TH0DPH,PH0DRR, & PH0DTH,PH0DPH,BB0DRR, & BB0DTH,BB0DPH,XLVDRR, & XLVDTH,XLVDPH RR = PATH(1) TH = PATH(2) PH = PATH(3) CALL GTOM1(TH,PH) CALL HOKAN(RR,TH,PH,TH0,PH0,BB0,XLV) IF(ABS(XLV-XLPP).LT.ERPP) THEN DRR = DRR2 DTH = DTH2 DPH = DPH2 ELSE DRR = DRR1 DTH = DTH1 DPH = DPH1 END IF RRDRR = PATH(1)+DRR THDTH = PATH(2)+DTH PHDTH = PATH(3) CALL GTOM1(THDTH,PHDTH) THDPH = PATH(2) PHDPH = PATH(3)+DPH CALL GTOM1(THDPH,PHDPH) DTH = DTH*RAD DPH = DPH*RAD CALL HOKAN(RRDRR,TH ,PH ,TH0DRR,PH0DRR,BB0DRR,XLVDRR) CALL HOKAN(RR ,THDTH,PHDTH,TH0DTH,PH0DTH,BB0DTH,XLVDTH) CALL HOKAN(RR ,THDPH,PHDPH,TH0DPH,PH0DPH,BB0DPH,XLVDPH) DT0DRR = (TH0DRR-TH0)/DRR DT0DTH = (TH0DTH-TH0)/DTH DT0DPH = (TH0DPH-TH0)/DPH DP0DRR = (PH0DRR-PH0)/DRR DP0DTH = (PH0DTH-PH0)/DTH DP0DPH = (PH0DPH-PH0)/DPH DB0DRR = (BB0DRR-BB0)/DRR DB0DTH = (BB0DTH-BB0)/DTH DB0DPH = (BB0DPH-BB0)/DPH DLVDRR = (XLVDRR-XLV)/DRR DLVDTH = (XLVDTH-XLV)/DTH DLVDPH = (XLVDPH-XLV)/DPH RETURN END ************************************************************************ SUBROUTINE HOKAN(RR,TH,PH,TH0,PH0,BB0,XLV) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) REAL*8 RR ,TH ,PH ,TH0,PH0,BB0,XLV COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QCOND /L,LL,NN,NLOOP,A COMMON /QHOKAN/II (2),JJ (2),KK (2),THJ(2) &,PHJ(2),BBJ(2),XLJ(2),THK(2),PHK(2), &BBK(2),XLK(2) C PARAMETER (INMAX=75,JNMIN=15,JNMAX=75,KNMIN=4,KNMAX=32) PARAMETER (INMAX=75,JNMIN=5,JNMAX=75,KNMIN=2,KNMAX=34) COMMON /QMESH /NF,IN ,DPHMSH,JMIN,JMAX ,DTH0 ,KMIN,KMAX ,DTHMSH,PH &MSH(0:INMAX),RRMSH( JNMIN:JNMAX,KNMIN:KNMAX),THMSH( KNMIN:KNMAX),T &HREF(INMAX,JNMIN:JNMAX,KNMIN:KNMAX),PHREF(INMAX,JNMIN:JNMAX,KNMIN: &KNMAX),BBBTM(INMAX,JNMIN:JNMAX,KNMIN:KNMAX),XLVAL(INMAX,JNMIN:JNMA &X,KNMIN:KNMAX) II(1)=(PH-PHMSH(0))/DPHMSH II(2)=II(1)+1 KK(1)=TH/DTHMSH IF(KK(1).LT.KMIN) THEN L =33 KK(1)=KMIN TH=THMSH(KK(1)) ELSE IF(KK(1).GE.KMAX) THEN L =34 KK(1)=KMAX-1 TH=THMSH(KK(1)) END IF KK(2)=KK(1)+1 K=1 CALL HOKAN2(PH,RR,K) K=2 CALL HOKAN2(PH,RR,K) TH1 = THMSH(KK(1)) TH2 = THMSH(KK(2)) DTH1 =(TH-TH1)/DTHMSH DTH2 =(TH2-TH)/DTHMSH TH0 = DTH1*THK(2)+DTH2*THK(1) PH0 = DTH1*PHK(2)+DTH2*PHK(1) BB0 = DTH1*BBK(2)+DTH2*BBK(1) DTH1 = DTH1*(SIN(TH2*RAD)/SIN(TH*RAD))**2 DTH2 = DTH2*(SIN(TH1*RAD)/SIN(TH*RAD))**2 XLV = DTH1*XLK(2)+DTH2*XLK(1) RETURN END ************************************************************************ SUBROUTINE HOKAN2(PH,RR,K) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) REAL*8 RR,PH COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QCOND /L,LL,NN,NLOOP,A COMMON /QHOKAN/II (2),JJ (2),KK (2),THJ(2) &,PHJ(2),BBJ(2),XLJ(2),THK(2),PHK(2), &BBK(2),XLK(2) C PARAMETER (INMAX=75,JNMIN=15,JNMAX=75,KNMIN=4,KNMAX=32) PARAMETER (INMAX=75,JNMIN=5,JNMAX=75,KNMIN=2,KNMAX=34) COMMON /QMESH /NF,IN ,DPHMSH,JMIN,JMAX ,DTH0 ,KMIN,KMAX ,DTHMSH,PH &MSH(0:INMAX),RRMSH( JNMIN:JNMAX,KNMIN:KNMAX),THMSH( KNMIN:KNMAX),T &HREF(INMAX,JNMIN:JNMAX,KNMIN:KNMAX),PHREF(INMAX,JNMIN:JNMAX,KNMIN: &KNMAX),BBBTM(INMAX,JNMIN:JNMAX,KNMIN:KNMAX),XLVAL(INMAX,JNMIN:JNMA &X,KNMIN:KNMAX) JJ(1)=ASIN(SQRT(RE/RR)*SIN(THMSH(KK(K))*RAD))*DEG/DTH0 IF(JJ(1).LT.JMIN) THEN L =35 JJ(1)=JMIN RR=RRMSH(JJ(1),KK(K)) ELSE IF(JJ(1).GE.JMAX) THEN L =36 JJ(1)=JMAX-1 RR=RRMSH(JJ(1),KK(K)) END IF JJ(2)=JJ(1)+1 J=1 CALL HOKAN1(PH,J,K) J=2 CALL HOKAN1(PH,J,K) RR1 = RRMSH(JJ(1),KK(K)) RR2 = RRMSH(JJ(2),KK(K)) DRR1 =(RR1-RR)/(RR1-RR2) DRR2 =(RR-RR2)/(RR1-RR2) THK(K)= DRR1*THJ(2)+DRR2*THJ(1) PHK(K)= DRR1*PHJ(2)+DRR2*PHJ(1) BBK(K)= DRR1*BBJ(2)+DRR2*BBJ(1) XLK(K)= DRR1*XLJ(2)+DRR2*XLJ(1) RETURN END ************************************************************************ SUBROUTINE HOKAN1(PH,J,K) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) REAL*8 PH COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QCOND /L,LL,NN,NLOOP,A COMMON /QHOKAN/II (2),JJ (2),KK (2),THJ(2) &,PHJ(2),BBJ(2),XLJ(2),THK(2),PHK(2), &BBK(2),XLK(2) C PARAMETER (INMAX=75,JNMIN=15,JNMAX=75,KNMIN=4,KNMAX=32) PARAMETER (INMAX=75,JNMIN=5,JNMAX=75,KNMIN=2,KNMAX=34) COMMON /QMESH /NF,IN ,DPHMSH,JMIN,JMAX ,DTH0 ,KMIN,KMAX ,DTHMSH,PH &MSH(0:INMAX),RRMSH( JNMIN:JNMAX,KNMIN:KNMAX),THMSH( KNMIN:KNMAX),T &HREF(INMAX,JNMIN:JNMAX,KNMIN:KNMAX),PHREF(INMAX,JNMIN:JNMAX,KNMIN: &KNMAX),BBBTM(INMAX,JNMIN:JNMAX,KNMIN:KNMAX),XLVAL(INMAX,JNMIN:JNMA &X,KNMIN:KNMAX) IF((THREF(II(1),JJ(J),KK(K)).EQ.999.0D0).OR. & (THREF(II(2),JJ(J),KK(K)).EQ.999.0D0)) L=37 PH1 = PHMSH(II(1)) PH2 = PHMSH(II(2)) DPH1 =(PH-PH1)/DPHMSH DPH2 =(PH2-PH)/DPHMSH THJ(J) = DPH1*THREF(II(2),JJ(J),KK(K)) & +DPH2*THREF(II(1),JJ(J),KK(K)) PHJ(J) = DPH1*PHREF(II(2),JJ(J),KK(K)) & +DPH2*PHREF(II(1),JJ(J),KK(K)) BBJ(J) = DPH1*BBBTM(II(2),JJ(J),KK(K)) & +DPH2*BBBTM(II(1),JJ(J),KK(K)) XLJ(J) = DPH1*XLVAL(II(2),JJ(J),KK(K)) & +DPH2*XLVAL(II(1),JJ(J),KK(K)) RETURN END ************************************************************************ SUBROUTINE MESH(INDMSH,IYEAR,IMNTH,IDATE) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) C PARAMETER (INMAX=75,JNMIN=15,JNMAX=75,KNMIN=4,KNMAX=32) PARAMETER (INMAX=75,JNMIN=5,JNMAX=75,KNMIN=2,KNMAX=34) COMMON /QMESH /NF,IN ,DPHMSH,JMIN,JMAX ,DTH0 ,KMIN,KMAX ,DTHMSH,PH &MSH(0:INMAX),RRMSH( JNMIN:JNMAX,KNMIN:KNMAX),THMSH( KNMIN:KNMAX),T &HREF(INMAX,JNMIN:JNMAX,KNMIN:KNMAX),PHREF(INMAX,JNMIN:JNMAX,KNMIN: &KNMAX),BBBTM(INMAX,JNMIN:JNMAX,KNMIN:KNMAX),XLVAL(INMAX,JNMIN:JNMA &X,KNMIN:KNMAX) C PARAMETER (NFMAX=1109) PARAMETER (NFMAX=1258) DIMENSION JJJQ(NFMAX),KKKQ(NFMAX),RRMSHQ(NFMAX),THREFQ(NFMAX),PHRE &FQ(NFMAX),BBBTMQ(NFMAX),XLVALQ(NFMAX) DO 10 K=KNMIN,KNMAX DO 11 J=JNMIN,JNMAX DO 12 I=1,INMAX THREF(I,J,K) = 999.0D0 PHREF(I,J,K) = 999.0D0 BBBTM(I,J,K) = 0.0D0 XLVAL(I,J,K) = 0.0D0 12 CONTINUE 11 CONTINUE 10 CONTINUE READ(INDMSH) RREFQ,RBTMQ,RREFEQ,RBTMEQ,STEP0Q,DTH0,JMIN,JMAX, & DTHMSH,KMIN,KMAX,IHEMQ,IYEAR,IMNTH,IDATE READ(INDMSH) IN,NFMAXQ,NF READ(INDMSH) JJJQ,KKKQ READ(INDMSH) RRMSHQ DO 1010 IFF=1,NF RRMSH(JJJQ(IFF),KKKQ(IFF))=RRMSHQ(IFF) 1010 CONTINUE DO 20 K=KMIN,KMAX THMSH(K)=DTHMSH*K 20 CONTINUE DO 30 II=1,IN READ(INDMSH) THREFQ READ(INDMSH) PHREFQ READ(INDMSH) XLVALQ READ(INDMSH) BBBTMQ READ(INDMSH) PHMSH(II) DO 1020 IFF=1,NF THREF(II,JJJQ(IFF),KKKQ(IFF)) = THREFQ(IFF) PHREF(II,JJJQ(IFF),KKKQ(IFF)) = PHREFQ(IFF) BBBTM(II,JJJQ(IFF),KKKQ(IFF)) = BBBTMQ(IFF) XLVAL(II,JJJQ(IFF),KKKQ(IFF)) = XLVALQ(IFF) 1020 CONTINUE 30 CONTINUE DPHMSH=PHMSH(2)-PHMSH(1) DO 40 II=3,IN IF((PHMSH(II)-PHMSH(II-1)).NE.DPHMSH)WRITE(6,'('' MESH DATA ERRO &R !!! : II ='',I3)') II 40 CONTINUE PHMSH(0)=PHMSH(1)-DPHMSH RETURN END ************************************************************************ SUBROUTINE GTOM0 ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QTHPH0/COST0,SINT0,COSP0,SINP0 entry mtog0 TH0 = 11.2 D0*RAD PH0 = 70.75D0*RAD COST0 = COS(TH0) SINT0 = SIN(TH0) COSP0 = COS(PH0) SINP0 = SIN(PH0) RETURN END ************************************************************************ SUBROUTINE GTOM1(TH,PH) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QTHPH0/COST0,SINT0,COSP0,SINP0 TH = TH*RAD PH = PH*RAD COSTH = COS(TH) SINTH = SIN(TH) COSPH = COS(PH) SINPH = SIN(PH) QUI = SINTH*(COSPH*COSP0-SINPH*SINP0) TH = ACOS(COSTH*COST0+SINT0*QUI) AAA = SINTH*(SINPH*COSP0+COSPH*SINP0) BBB = COST0*QUI-SINT0*COSTH PH = ATAN2(AAA,BBB) TH = TH*DEG PH = PH*DEG RETURN END ************************************************************************ SUBROUTINE GTOM2(TH,PH,EPS,DEL) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QTHPH0/COST0,SINT0,COSP0,SINP0 TH = TH *RAD PH = PH *RAD EPS = EPS*RAD DEL = DEL*RAD COSTH = COS(TH) SINTH = SIN(TH) COSPH = COS(PH) SINPH = SIN(PH) QUI = SINTH*(COSPH*COSP0-SINPH*SINP0) TH = ACOS(COSTH*COST0+SINT0*QUI) AAA = SINTH*(SINPH*COSP0+COSPH*SINP0) BBB = COST0*QUI-SINT0*COSTH PH = ATAN2(AAA,BBB) COSTHN = COS(TH) SINTHN = SIN(TH) COSPHN = COS(PH) SINPHN = SIN(PH) ER = COS(EPS)*COS(DEL) ET = COS(EPS)*SIN(DEL) EP = SIN(EPS) ERN = SINTH *COSPH*ER + COSTH *COSPH *ET - SINPH *EP ETN = SINTH *SINPH*ER + COSTH *SINPH *ET + COSPH *EP EPN = COSTH *ER - SINTH *ET ER = COST0 *COSP0*ERN - COST0 *SINP0 *ETN - SINT0 *EPN ET = SINP0*ERN + COSP0 *ETN EP = SINT0 *COSP0*ERN - SINT0 *SINP0 *ETN + COST0 *EPN ERN = SINTHN*COSPHN*ER + SINTHN*SINPHN*ET + COSTHN*EP ETN = COSTHN*COSPHN*ER + COSTHN*SINPHN*ET - SINTHN*EP EPN =- SINPHN*ER + COSPHN*ET EPS = ATAN(EPN/SQRT(ETN**2+ERN**2)) DEL = ATAN2(ETN,ERN) TH = TH *DEG PH = PH *DEG EPS = EPS*DEG DEL = DEL*DEG RETURN END ************************************************************************ SUBROUTINE MTOG1(TH,PH) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QTHPH0/COST0,SINT0,COSP0,SINP0 TH = TH*RAD PH = PH*RAD COSTH = COS(TH) SINTH = SIN(TH) COSPH = COS(PH) SINPH = SIN(PH) TH = ACOS(COST0*COSTH-SINT0*SINTH*COSPH) AAA = -COST0*SINP0*SINTH*COSPH+COSP0*SINTH*SINPH-SINT0*SINP0*COSTH & BBB = COST0*COSP0*SINTH*COSPH+SINP0*SINTH*SINPH+SINT0*COSP0*COSTH PH = ATAN2(AAA,BBB) TH = TH*DEG PH = PH*DEG RETURN END ************************************************************************ SUBROUTINE MTOG2(TH,PH,EPS,DEL) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QTHPH0/COST0,SINT0,COSP0,SINP0 TH = TH *RAD PH = PH *RAD EPS = EPS*RAD DEL = DEL*RAD COSTH = COS(TH) SINTH = SIN(TH) COSPH = COS(PH) SINPH = SIN(PH) TH = ACOS(COST0*COSTH-SINT0*SINTH*COSPH) AAA =-COST0*SINP0*SINTH*COSPH+COSP0*SINTH*SINPH-SINT0*SINP0*COSTH BBB = COST0*COSP0*SINTH*COSPH+SINP0*SINTH*SINPH+SINT0*COSP0*COSTH PH = ATAN2(AAA,BBB) COSTHN = COS(TH) SINTHN = SIN(TH) COSPHN = COS(PH) SINPHN = SIN(PH) ER = COS(EPS)*COS(DEL) ET = COS(EPS)*SIN(DEL) EP = SIN(EPS) ERN = SINTH *COSPH *ER + COSTH *COSPH *ET - SINPH *EP ETN = SINTH *SINPH *ER + COSTH *SINPH *ET + COSPH *EP EPN = COSTH *ER - SINTH *ET ER = COST0 *COSP0 *ERN + SINP0 *ETN + SINT0 *COSP0*EPN ET =-COST0 *SINP0 *ERN + COSP0 *ETN - SINT0 *SINP0*EPN EP =-SINT0 *ERN + COST0 *EPN ERN = SINTHN*COSPHN*ER + SINTHN*SINPHN*ET + COSTHN *EP ETN = COSTHN*COSPHN*ER + COSTHN*SINPHN*ET - SINTHN *EP EPN =- SINPHN*ER + COSPHN*ET EPS = ATAN(EPN/SQRT(ETN**2+ERN**2)) DEL = ATAN2(ETN,ERN) TH = TH *DEG PH = PH *DEG EPS = EPS*DEG DEL = DEL*DEG RETURN END