C **** START.3 C **** STARTS THE MODEL FROM HISTORY FILE START.4 C **** START.5 INCLUDE 'param.h' START.6 INCLUDE 'buffe.h' START.7 INCLUDE 'const.h' START.8 INCLUDE 'index.h' START.9 INCLUDE 'start.h' START.10 INCLUDE 'epara.h' START.11 INCLUDE 'eterp.h' EINPUT.9 INCLUDE 'edphi.h' EINPUT.10 DIMENSION POT(IMAXGP,0:JMAXGP,29),HHT(IMAXGP,JMAXG,29) DIMENSION DYNPOTN(73,0:37,29),HT(73,36,29),PHIMN(IMAXMP,JMAXM) C **** START.19 C **** DEFINE MODEL CONSTANTS START.20 C **** START.21 CALL CON START.22 C **** START.23 C **** CALCULATE DYNAMO PARAMETERS START.24 C **** START.25 CALL CNSTNT START.26 C **** START.120 C **** READ IN MAGNETIC COORDINATE TRANSFORMATION METRIX START.121 C **** START.122 CALL EINPUT START.123 C **** START.27 C **** ZERO SCM BUFFER AREA START.28 C **** START.29 IHIST=3 NHIST=24 OPEN(UNIT=IHIST,FILE='ehist',FORM='UNFORMATTED') START.37 c ***** START.98 c ***** READ IN THE DATA START.99 c ***** START.100 DO J=1,JMAX START.101 READ(IHIST)((F(I,K),I=1,IMAXP4),K=1,KMAXP1*NHIST+1) START.102 C **** START.103 C **** STORE F ARRAY INTO G ARRRAY START.104 C **** START.105 DO K=1,KMAXP1*NHIST+1 START.106 DO I=1,IMAXP4 START.107 G(I,K,J)=F(I,K) START.108 ENDDO START.109 ENDDO START.110 NPHIK=NPHI+1 START.111 NZK=NZ+1 DO K=1,KMAXP1 START.112 DO I=1,IMAX+1 START.113 DYNPOTN(I,J,K)=G(I+2,NPHIK,J) START.114 HT(I,J,K)=G(I+2,NZK,J) ENDDO START.115 NPHIK=NPHIK+1 START.116 NZK=NZK+1 ENDDO START.117 ENDDO START.118 C **** EINPUT.60 C **** INSERT POLAR POINTS IN DYNPOT EINPUT.61 C **** EINPUT.62 DO K=1,KMAXP1 EINPUT.63 TEMPA=0. EINPUT.64 TEMPB=0. EINPUT.65 DO I=1,73 EINPUT.66 TEMPA=TEMPA+DYNPOTN(I,1,K) EINPUT.67 TEMPB=TEMPB+DYNPOTN(I,2,K) EINPUT.68 ENDDO EINPUT.69 DYNPOTN(1,0,K)=(9.*TEMPA-TEMPB)/(8.*DBLE(73)) EINPUT.70 TEMPA=0. EINPUT.71 TEMPB=0. EINPUT.72 DO I=1,73 EINPUT.73 TEMPA=TEMPA+DYNPOTN(I,36,K) EINPUT.74 TEMPB=TEMPB+DYNPOTN(I,35,K) EINPUT.75 ENDDO EINPUT.76 DYNPOTN(1,37,K)=(9.*TEMPA-TEMPB)/(8.*DBLE(73)) EINPUT.77 DO I=2,73 EINPUT.78 DYNPOTN(I,0,K)=DYNPOTN(1,0,K) EINPUT.79 DYNPOTN(I,37,K)=DYNPOTN(1,37,K) EINPUT.80 ENDDO EINPUT.81 ENDDO EINPUT.82 C **** EINPUT.60 C **** LINEAR INTERPOLATE POTENTIAL TO HIGH RESOLUTION EINPUT.61 C **** EINPUT.62 OPEN(18,FILE='Pot-Height.dat',FORM='UNFORMATTED') OPEN(19,FILE='Pre-Height.dat',FORM='UNFORMATTED') DO K=1,KMAXP1 DO J=1,36 C **** C **** FIRST LINE IN LATITUDE C **** IF(J.EQ.1)THEN JJ=5+9*(J-1) DO I=1,73 II=9*(I-1)+1 IF(I.LT.73)THEN POT(II,JJ,K)=DYNPOTN(I,J,K) HHT(II,JJ,K)=HT(I,J,K) DO L=1,8 IL=II+L POT(IL,JJ,K)=DYNPOTN(I,J,K)+DBLE(L)*(DYNPOTN(I+1,J,K)- 1 DYNPOTN(I,J,K))/9. HHT(IL,JJ,K)=HT(I,J,K)+DBLE(L)*(HT(I+1,J,K)- 1 HT(I,J,K))/9. ENDDO ELSE POT(II,JJ,K)=DYNPOTN(I,J,K) HHT(II,JJ,K)=HT(I,J,K) ENDIF ENDDO ELSE C **** C **** NEXT LINE IN LATITUDE C **** J2=5+9*(J-1) J1=5+9*(J-2) DO I=1,73 II=9*(I-1)+1 IF(I.LT.73)THEN POT(II,J2,K)=DYNPOTN(I,J,K) c write(6,*)POT(II,J2,K),ii,j2,k HHT(II,J2,K)=HT(I,J,K) DO L=1,8 IL=II+L POT(IL,J2,K)=DYNPOTN(I,J,K)+ 1 DBLE(L)*(DYNPOTN(I+1,J,K)-DYNPOTN(I,J,K))/9. c write(6,*)POT(IL,J2,K),il,j2,k HHT(IL,J2,K)=HT(I,J,K)+ 1 DBLE(L)*(HT(I+1,J,K)-HT(I,J,K))/9. ENDDO ELSE POT(II,J2,K)=DYNPOTN(I,J,K) c write(6,*)POT(II,J2,K),ii,j2,k HHT(II,J2,K)=HT(I,J,K) ENDIF c pause ENDDO C **** C **** CALCULATE LATITUDINAL LINES IN BETWEEN C **** DO I=1,IMAXGP DO L=1,8 POT(I,J2-L,K)=POT(I,J2,K)-DBLE(L)*(POT(I,J2,K)-POT(I,J1,K))/9. HHT(I,J2-L,K)=HHT(I,J2,K)-DBLE(L)*(HHT(I,J2,K)-HHT(I,J1,K))/9. ENDDO ENDDO ENDIF ENDDO C **** C **** POLES C **** JJ=0 DO I=1,73 II=9*(I-1)+1 IF(I.LT.73)THEN POT(II,JJ,K)=DYNPOTN(I,0,K) DO L=1,8 IL=II+L POT(IL,JJ,K)=DYNPOTN(I,0,K)+DBLE(L)*(DYNPOTN(I+1,0,K)- 1 DYNPOTN(I,0,K))/9. ENDDO ELSE POT(II,JJ,K)=DYNPOTN(I,0,K) ENDIF ENDDO DO I=1,IMAXGP DO J=1,4 POT(I,J,K)=POT(I,JJ,K)+ 1 DBLE(2*(J-1)+1)*(POT(I,JJ+5,K)-POT(I,JJ,K))/5. ENDDO ENDDO JJ=JMAXGP DO I=1,73 II=9*(I-1)+1 IF(I.LT.73)THEN POT(II,JJ,K)=DYNPOTN(I,37,K) DO L=1,8 IL=II+L POT(IL,JJ,K)=DYNPOTN(I,37,K)+DBLE(L)*(DYNPOTN(I+1,37,K)- 1 DYNPOTN(I,37,K))/9. ENDDO ELSE POT(II,JJ,K)=DYNPOTN(I,37,K) ENDIF ENDDO DO I=1,IMAXGP DO J=1,4 POT(I,JMAXGP-J,K)=POT(I,JJ,K)- 1 DBLE(2*(J-1)+1)*(POT(I,JJ,K)-POT(I,JJ-5,K))/9. ENDDO ENDDO c do j=0,jmaxgp c do i=1,imaxgp c write(6,*)pot(i,j,k),i,j,k c enddo c pause c enddo C **** EINPUT.87 C **** TRANSFORM DYNPOT TO GEOMAGNETIC COORDINATES IN EINPUT.88 C **** PHIM3D(IMAXMP,JMAXM,-2:ZKMXP) EINPUT.89 C **** EINPUT.90 write(6,*)'ggg',k DO J=1,JMAXM EINPUT.92 CALL GRDINT(PHIMN(1,J),POT(1,0,K),IG,JG,WT,IMAXGP, EINPUT.93 1 IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J) EINPUT.94 ENDDO EINPUT.95 C **** EINPUT.98 C **** PERIODIC POINTS EINPUT.99 C **** EINPUT.100 DO J=1,JMAXM EINPUT.102 PHIMN(IMAXMP,J)=PHIMN(1,J) EINPUT.103 ENDDO EINPUT.105 J1=320 DO J=1,4 JJ=J1+J DO I=1,IMAXGP HHT(I,J,K)=HHT(I,5,K)-FLOAT(5-J)*(HHT(I,5,K)- 1 HHT(1+MOD(I-1+IMAXG/2,IMAXG),5,K))/DBLE(9) HHT(I,JJ,K)=HHT(I,J1,K)-FLOAT(5-J)*(HHT(I,J1,K)- 1 HHT(1+MOD(I-1+IMAXG/2,IMAXG),J1,K))/DBLE(9) ENDDO ENDDO c do j=1,jmaxm c do i=1,imaxmp c write(6,*)PHIMN(I,J),i,j c enddo c pause c enddo write(18)((PHIMN(I,J),I=1,IMAXMP),J=1,JMAXM) c do j=1,jmaxg c do i=1,imaxgp c write(19)HHT(I,J,K) c enddo c enddo c do j=1,jmaxg c do i=1,imaxgp c write(6,*)HHT(I,J,K),i,j,k c enddo c pause c enddo c write(19)((HHT(I,J,K),I=1,IMAXGP),J=1,JMAXG) ENDDO do kk=1,3 write(19)((HHT(I,J,1),I=1,IMAXGP),J=1,JMAXG) enddo do k=1,kmaxp1 write(19)((HHT(I,J,k),I=1,IMAXGP),J=1,JMAXG) enddo CLOSE(IHIST) START.120 STOP START.134 END START.135 C START.136 *DK BLOCKD.1 BLOCK DATA CONSTS BLOCKD.2 INCLUDE 'param.h' BLOCKD.3 INCLUDE 'const.h' BLOCKD.4 INCLUDE 'start.h' BLOCKD.5 DATA NFLDS/32/,NHIST/24/,ITAP/1/ BLOCKD.6 DATA KPOLE/1,-1,-1,1,1,1,1,1,0,1,1,0,0,1,1,1,1,1,-1,-1,1,1,1,1,1, BLOCKD.7 11,1,0,0,0,0,0/ BLOCKD.8 DATA KFLDS/16*3,1,10*3,5*2/ BLOCKD.9 DATA RMASS/32.,16.,28./ BLOCKD.10 DATA KUT/1,2,3,5,6,7,9,10,11,13,14,15,17,10*17,17,15,14,13,11, BLOCKD.11 110,9,7,6,5,3,2,1/ BLOCKD.12 DATA NESTIN/11,12,13,14/,NESTOT/21,22,23,24/ DATA NESTRD/'n1in','n2in','n3in','n4in'/ DATA NESTWT/'n1ot','n2ot','n3ot','n4ot'/ END BLOCKD.13 C BLOCKD.14 *DK CNSTNT CNSTNT.1 SUBROUTINE CNSTNT CNSTNT.2 C **** CNSTNT.3 C **** DEFINES AND CALCULATES CONSTANTS NEEDED BY PROGRAM CNSTNT.4 C **** CNSTNT.5 INCLUDE 'param.h' CNSTNT.6 INCLUDE 'epara.h' CNSTNT.7 INCLUDE 'econs.h' CNSTNT.8 DATA ALFA/1.668/ CNSTNT.9 DIMENSION TABL2(181,3:5),THETA0(JMAXM) CNSTNT.11 C **** CNSTNT.13 C **** NOW SET CONSTANTS CNSTNT.14 C **** CNSTNT.15 RE=6.378165D8 CNSTNT.16 H0=9.0D6 CNSTNT.17 H00=9.7D6 CNSTNT.18 R1=1.06D7 CNSTNT.24 R0=RE+H0 CNSTNT.25 R00=RE+H00 CNSTNT.26 PI=4*DATAN(1.D0) CNSTNT.28 DTR=PI/180. CNSTNT.29 RTD=180./PI CNSTNT.30 DLATM=PI/DBLE(JMAXM-1) CNSTNT.31 DLONM=2.*PI/DBLE(IMAXM) CNSTNT.32 C **** CNSTNT.33 C **** FILL ARRAY YLATM (EQUALLY SPACED IN THETAO BUT HOLDS CNSTNT.34 C **** CORRESPONDING VALUE OF THETAS) CNSTNT.35 C **** CNSTNT.36 DO J=1,JMAXM CNSTNT.37 C **** CNSTNT.38 C **** THETA0 = EQUALLY SPACED GRID VALUES CNSTNT.39 C **** CNSTNT.40 THETA0(J)=-PI/2.+DBLE(J-1)*DLATM CNSTNT.41 ENDDO CNSTNT.42 DO J=2,JMAXM-1 CNSTNT.43 C **** CNSTNT.44 C **** TANTH0 = ABS(TAN(THETA0)) CNSTNT.46 C **** CNSTNT.47 TANTH0=DABS(DTAN(THETA0(J))) CNSTNT.48 C **** CNSTNT.49 C **** HAMH0 = HA-H0 CNSTNT.50 C **** CNSTNT.51 HAMH0=R1*TANTH0+R0*TANTH0**(2+2*ALFA)/ CNSTNT.52 1 (1+TANTH0**2)**ALFA CNSTNT.53 C **** CNSTNT.54 C **** TANTHS = ABS(TAN(THETAS)) CNSTNT.55 C **** CNSTNT.56 TANTHS=DSQRT(HAMH0/R0) CNSTNT.57 C **** CNSTNT.58 C **** YLATM = TANTHS CNSTNT.59 C **** CNSTNT.60 YLATM(J)=DSIGN(DATAN(TANTHS),THETA0(J)) CNSTNT.61 C **** CNSTNT.62 C **** RCOS0S = COS(THETA0)/COS(THETAS) CNSTNT.63 C **** CNSTNT.64 RCOS0S(J)=DSQRT((1+TANTHS**2)/(1+TANTH0**2)) CNSTNT.65 C **** CNSTNT.66 C **** DT0DTS = D(THETA0)/D(THETAS) CNSTNT.67 C **** CNSTNT.68 C **** DT1DTS = DT0DTS DIVIDED BY ABS[SIN(Im)]. REMIANS FINITE CNSTNT.69 C **** AND NONZERO AY MAGNETIC EQUATOR. CNSTNT.70 C **** CNSTNT.71 TANTHS2=TANTHS**2 CNSTNT.72 DT1DTS(J)=(R0*DSQRT(1+4*TANTHS2)*(1+TANTHS2))/ CNSTNT.73 1 (R1*(1+TANTH0**2)+2*R0*TANTH0**(2*ALFA+1)* CNSTNT.74 2 (1+ALFA+TANTH0**2)/(1+TANTH0**2)**ALFA) CNSTNT.75 DT0DTS(J)=DT1DTS(J)*2*TANTHS/DSQRT(1+4*TANTHS2) CNSTNT.76 ENDDO CNSTNT.77 C **** CNSTNT.78 C **** NOW DO POLES CNSTNT.79 C **** CNSTNT.80 YLATM(1)=THETA0(1) CNSTNT.81 YLATM(JMAXM)=THETA0(JMAXM) CNSTNT.82 RCOS0S(1)=1. CNSTNT.83 RCOS0S(JMAXM)=1. CNSTNT.84 DT0DTS(1)=1. CNSTNT.85 DT0DTS(JMAXM)=1. CNSTNT.86 DT1DTS(1)=1. CNSTNT.87 DT1DTS(JMAXM)=1. CNSTNT.88 C **** CNSTNT.89 C **** FILL YLONM CNSTNT.90 C **** CNSTNT.91 DO I=1,IMAXMP CNSTNT.92 YLONM(I)=-PI+DBLE(I-1)*DLONM CNSTNT.93 ENDDO CNSTNT.94 DTHETA=PI/(2.*180.) CNSTNT.95 TABLE(1,1)=0. CNSTNT.96 TABLE(1,2)=0. CNSTNT.97 DO I=2,181 CNSTNT.98 TABLE(I,1)=TABLE(I-1,1)+DTHETA CNSTNT.99 ENDDO CNSTNT.100 DO I=2,180 CNSTNT.101 TABL2(I,4)=DTAN(TABLE(I,1)) CNSTNT.102 TABLE(I,2)=TABLE(I,1) CNSTNT.103 ENDDO CNSTNT.104 DO N=1,7 CNSTNT.105 DO I=2,180 CNSTNT.106 TABL2(I,3)=TABLE(I,2) CNSTNT.107 TABLE(I,2)=DTAN(TABL2(I,3)) CNSTNT.108 TABL2(I,5)=DSQRT(R1/R0*TABLE(I,2)+TABLE(I,2)** CNSTNT.109 1 (2*(1+ALFA))/(1+TABLE(I,2)**2)**ALFA) CNSTNT.110 TABLE(I,2)=TABL2(I,3)-(TABL2(I,5)-TABL2(I,4))*2*TABL2(I,5)/ CNSTNT.111 2 (R1/R0*(1+TABLE(I,2)**2)+2*TABLE(I,2)**(2*ALFA+1)* CNSTNT.112 3 (1+ALFA+TABLE(I,2)**2)/(1+TABLE(I,2)**2)**ALFA) CNSTNT.113 ENDDO CNSTNT.114 ENDDO CNSTNT.115 RETURN CNSTNT.116 END CNSTNT.117 C CNSTNT.118 *DK CON CON.1 SUBROUTINE CON CON.2 C **** CON.3 C **** DEFINES THE MODEL CONSTANTS CON.4 C **** CON.5 INCLUDE 'param.h' CON.6 INCLUDE 'buffe.h' CON.7 INCLUDE 'crfft.h' CON.8 INCLUDE 'const.h' CON.9 INCLUDE 'index.h' CON.10 INCLUDE 'start.h' CON.11 INCLUDE 'field.h' CON.12 DATA NPHYS/39/ CON.13 C **** CON.14 C **** RESOLUTION PARAMETERS CON.15 C **** CON.16 IMAX=ZIMX CON.17 JMAX=ZJMX CON.18 KMAX=ZKMX CON.19 C **** CON.20 C **** DERIVED CONSTANTS CON.21 C **** CON.22 PI=3.14159265358979 CON.23 DLAMDA=2*PI/IMAX CON.24 DPHI=PI/JMAX CON.25 DS=(ZST-ZSB)/KMAX CON.26 IMAXH=IMAX/2 CON.27 IMAXP2=IMAX+2 CON.28 IMAXP4=IMAX+4 CON.29 JMAXM1=JMAX-1 CON.31 JMAXM2=JMAX-2 CON.32 KMAXM1=KMAX-1 CON.33 KMAXP1=KMAX+1 CON.34 C **** CON.35 C **** BUFFER INFORMATION AND INITIALIZE POINTERS CON.36 C **** CON.37 LEN1=IMAXP4 CON.38 LEN2=IMAXP4*KMAX CON.39 LEN3=IMAXP4*KMAXP1 CON.40 NCBUF=0 CON.41 DO N=1,NFLDS CON.42 IF(N.EQ.17)THEN CON.43 K=1 CON.44 ELSE CON.45 K=KMAXP1 CON.46 ENDIF CON.47 NDEXA(N+1)=NCBUF CON.48 NCBUF=NCBUF+K CON.49 ENDDO CON.50 NCPHYS=0 CON.51 DO N=1,NPHYS CON.52 K=KMAXP1 CON.53 NDEXB(N+1)=NCPHYS CON.54 NCPHYS=NCPHYS+K CON.55 ENDDO CON.56 C **** CON.57 C **** COMPUTE COLUMN POINTERS FOR SCM BUFFERS CON.58 NBADDR(1)=1 CON.59 DO I=2,8 CON.60 NBADDR(I)=NBADDR(I-1)+NCBUF CON.61 ENDDO CON.62 NBAPHY=NBADDR(8)+NCBUF CON.63 DO N=1,NPHYS CON.64 NDEXB(N+1)=NDEXB(N+1)+NBAPHY CON.65 ENDDO CON.66 C **** CON.67 C **** OUTPUT POINTER INFORMATION CON.68 C **** CON.69 LIM=NFLDS+NPHYS+2 CON.70 WRITE(6,1005) (NDEXA(N),N=1,LIM) CON.71 1005 FORMAT(//' POINTER VALUES FOR EACH FIELD IN COMMON BLOCK NDEX', CON.72 A //(10X,5I5)) CON.73 C **** CON.78 C **** COMPUTE PHI VARIATION ARRAYS CON.79 C **** CON.80 DELT=86400/DBLE(IMAX) CON.81 JMAXH=JMAX/2 CON.82 OMEGA=7.292D-5 CON.83 DO J=1,JMAX CON.84 PHI=(J-JMAXH-0.5)*DPHI CON.85 CS(J)=DCOS(PHI) CON.86 SN(J)=DSIN(PHI) CON.87 TN(J)=DTAN(PHI) CON.88 COR(J)=2.*OMEGA*SN(J) CON.89 ENDDO CON.90 CSSP(1)=-CS(2) CON.91 CSSP(2)=-CS(1) CON.92 CSNP(1)=-CS(JMAX) CON.93 CSNP(2)=-CS(JMAXM1) CON.94 KUTSP(1)=KUT(2) CON.95 KUTSP(2)=KUT(1) CON.96 KUTNP(1)=KUT(JMAX) CON.97 KUTNP(2)=KUT(JMAXM1) CON.98 C **** CON.99 C **** CALCULATE SUN'S DECLINATION CON.100 C **** (WILL HAVE TO RECALCULATE THIS IF DAY IS ADVANCED) CON.101 C **** CON.102 DELTA=DATAN(DTAN(23.5*PI/180.)*DSIN(2.*PI*DBLE(IIDAY-80)/365.)) CON.103 C **** CON.104 C **** GRID CONSTANTS CON.105 C **** CON.106 C(1)=DLAMDA CON.107 C(2)=DPHI CON.108 C(3)=DS CON.109 C(6)=2.*C(4) CON.110 C(7)=1./C(6) CON.111 C **** CON.112 C **** FINITE DIFFERENCE OPERATOR CONSTANTS CON.113 C **** CON.114 C(10)=2./(3.*DLAMDA) CON.115 C(11)=1./(12.*DLAMDA) CON.116 C(12)=2./(3.*DPHI) CON.117 C(13)=1./(12.*DPHI) CON.118 C(14)=1./(2.*DLAMDA) CON.119 C(15)=1./(2.*DPHI) CON.120 C **** CON.121 C **** INTEGRATION CONSTANTS CON.122 C **** CON.123 C(18)=C(3)/(1.+.5*C(3)) CON.124 C(19)=(1.-.5*C(3))/(.1+.5*C(3)) CON.125 C(20)=1./DS CON.126 C(22)=DCOS(PI/3.) CON.127 C(23)=4.*PI/(24.*60.*60.) CON.128 C(24)=28.9 CON.129 C **** CON.130 C **** ANNUAL FREQUENCY CON.131 C **** CON.132 C(25)=C(23)/(2.*365.25) CON.133 C **** CON.134 C **** SHAPIRO SMOOTHER CONSTANT CON.135 C **** CON.136 C(26)=3.0D-2 CON.137 C **** CON.138 C **** EDDY DIFFUSION CON.139 C **** CON.140 C(27)=5.D-6 CON.141 C **** CON.142 C **** EDDY THERMAL CONDUCTION CON.143 C **** CON.144 C(28)=C(27) CON.145 C **** CON.146 C **** TIME SMOOTHING CONSTANT CON.147 C **** CON.148 ALPHA=.95 CON.149 C(30)=ALPHA CON.150 C(31)=.5*(1.-ALPHA) CON.151 C **** CON.152 C ***** SURPLUS HEAT PER EVENT CONVERTED FROM E.V. TO ERGS CON.153 C **** CON.154 C(45)=5.11*1.602D-12 CON.155 C **** CON.156 C **** PHYSICAL CONSTANTS CON.157 C **** CON.158 C **** EARTHS RADIUS CON.159 C **** CON.160 C(51)=6.37122D8 CON.161 C(52)=1./C(51) CON.162 C(53)=C(51)*C(51) CON.163 C **** CON.164 C **** GRAVITATIONAL CONSTANT CON.165 C **** CON.166 C(54)=870. CON.167 C(55)=1./C(54) CON.168 C **** CON.169 C **** EARTHS ANGULAR VELOCITY CON.170 C **** CON.171 C(56)=OMEGA CON.172 C **** CON.173 C **** GAS CONSTANT CON.174 C **** CON.175 C(57)=8.314D7 CON.176 C **** CON.177 C **** H*C CON.178 C **** CON.179 C(60)=1.9845D-16 CON.180 C **** CON.181 C **** G/R CON.182 C **** CON.183 C(65)=C(54)/C(57) CON.184 C **** CON.185 C **** STANDARD PRESSURE CON.186 C **** CON.187 C(81)=5.0D-4 CON.188 C **** CON.189 C **** BOLTZMANN CONST CON.190 C **** CON.191 C(84)=1.38D-16 CON.192 C **** CON.193 C **** AVOGRADO NUMBER CON.194 C **** CON.195 C(85)=6.023D23 CON.196 C **** CON.197 C **** EXP(-.5*DS) CON.198 C **** CON.199 C(86)=DEXP(-.5*DS) CON.200 C **** CON.201 C **** EXP(.5*DS) CON.202 C **** CON.203 C(87)=1./C(86) CON.204 C **** CON.205 C **** SUN'S DECLINATION: CON.206 C **** (WILL HAVE TO RECALCULATE IF DAY IS ADVANCED) CON.207 C **** CON.208 C(95)=DSIN(DELTA) CON.209 C(96)=DCOS(DELTA) CON.210 C **** CON.211 C **** KO COEFFICIENT OF EDDY VISCOSITY CON.212 C **** CON.213 C(109)=0.2 CON.214 C(110)=PI CON.215 C **** CON.216 C **** SAVE ORIGINAL DT VALUE CON.217 C **** CON.218 C(112)=DSQRT(PI) CON.219 C **** CON.220 C **** CALCULATE EXPS AND A ARRAYS CON.221 C **** CON.222 S=ZSB+.5*DS CON.223 EXPS(1)=DEXP(-S) CON.224 EXPDS=DEXP(-DS) CON.225 A(1)=1.D-6*DEXP(-7.-S) CON.226 DIFK(1)=C(27) CON.227 DIFT(1)=C(28) CON.228 XMUE(1)=C(27) CON.229 DO 1 K=2,KMAX CON.230 EXPS(K)=EXPS(K-1)*EXPDS CON.231 DIFK(K)=DIFK(K-1)*EXPDS CON.232 DIFT(K)=DIFT(K-1)*EXPDS CON.233 XMUE(K)=DIFK(K) CON.234 1 CONTINUE CON.235 DIFK(KMAXP1)=DIFK(KMAX)*EXPDS CON.236 DIFT(KMAXP1)=DIFT(KMAX)*EXPDS CON.237 XMUE(KMAXP1)=DIFK(KMAXP1) CON.238 RETURN CON.239 END CON.240 C CON.241 *DK EINPUT EINPUT.1 SUBROUTINE EINPUT EINPUT.2 INCLUDE 'param.h' EINPUT.3 INCLUDE 'const.h' EINPUT.4 INCLUDE 'epara.h' EINPUT.5 INCLUDE 'econs.h' EINPUT.6 INCLUDE 'efild.h' EINPUT.8 INCLUDE 'eterp.h' EINPUT.9 INCLUDE 'edphi.h' EINPUT.10 C **** EINPUT.11 C **** READ IN GEOMAGNETIC COORDINATE TRANSFORM METRIX EINPUT.12 C **** EINPUT.13 IMAG=23 EINPUT.14 OPEN(IMAG,FILE='magnet.dat',FORM='UNFORMATTED') EINPUT.15 DO I=1,IMAXGP*(JMAXGP+1) EINPUT.16 READ(IMAG)ALATM(I,0),ALONM(I,0),XB(I,0),YB(I,0),ZZB(I,0), EINPUT.17 1 BMOD(I,0) EINPUT.18 ENDDO EINPUT.19 DO I=1,IMAXGP*(JMAXGP+1)*4 EINPUT.20 READ(IMAG)RJAC(I,0,1,1) EINPUT.21 ENDDO EINPUT.22 DO I=1,IMAXGP*(JMAXGP+1)*6 EINPUT.23 READ(IMAG)AV(I,0,1,1) EINPUT.24 ENDDO EINPUT.25 DO J=1,JMAXM EINPUT.29 DO I=1,IMAXMP EINPUT.30 READ(IMAG)IG(I,J),JG(I,J) EINPUT.31 ENDDO EINPUT.32 ENDDO EINPUT.33 DO J=1,JMAXM EINPUT.34 DO I=1,IMAXMP EINPUT.35 READ(IMAG)WT(1,I,J),WT(2,I,J),WT(3,I,J),WT(4,I,J) EINPUT.36 ENDDO EINPUT.37 ENDDO EINPUT.38 DO I=1,IMAXGP*JMAXG*6 EINPUT.51 READ(IMAG)DVEC(I,1,1,1) EINPUT.52 ENDDO EINPUT.53 DO I=1,IMAXGP*JMAXG EINPUT.54 READ(IMAG)DDDARR(I,1) EINPUT.55 ENDDO EINPUT.56 DO I=1,IMAXGP*JMAXG EINPUT.57 READ(IMAG)BE3ARR(I,1) EINPUT.58 ENDDO EINPUT.59 write(6,*)IMAXGP,JMAXG,'done reading' CLOSE(IMAG) EINPUT.83 RETURN EINPUT.106 END EINPUT.107 C EINPUT.108 *DK GRIDINT GRDINT.1 SUBROUTINE GRDINT(AM,AG,IG,JG,WT,LG,LM,IMAXG,JMAXG,IMAXM,JMAXM, GRDINT.2 1JM) GRDINT.3 IMPLICIT REAL*8 (A-H,O-Z) GRDINT.4 C **** GRDINT.5 C **** TRANSFORM SCALAR FIELD GIVEN ON GEOGRAPHIC GRID TO GRDINT.6 C **** GEOMAGNETIC GRID USING INDICES AND WEIGHTS GENERATED GRDINT.7 C **** BY GRDSTM. GRDINT.8 C **** GRDINT.9 C **** INPUT: GRDINT.10 C **** AG(LG,1) IS 2-DIM SCALAR FIELD TO BE TRANSFORMED TO GRDINT.11 C **** GEOMAGNETIC GRID. GRDINT.12 C **** NB: PERIODIC POINT IS REPEATED SO THAT GRDINT.13 C **** AG(IMAXG+1,JG) = AG(1,JG) GRDINT.14 C **** (LM.GE.IMAXG+1) GRDINT.15 C **** IMAXG, JMAXG DIMENSIONS OF GEOGRAPHIC GRID GRDINT.16 C **** IMAXM, JMAXM DIMENSIONS OF GEOMAGNETIC GRID GRDINT.17 C **** LG IS FIRST DIMENSION OF AG IN CALLING PROGRAM GRDINT.18 C **** (LG.GE.IMAXG+1) GRDINT.19 C **** IG(LM,1) GIVES ROUNDED DOWN GEOGRAPHIC LONGITUDE INDEX GRDINT.20 C **** FOR EACH GEOMAGNETIC GRID POINT GRDINT.21 C **** JG(LM,1) GIVES ROUNDED DOWN GEOGRAPHIC LATITUDE INDEX GRDINT.22 C **** FOR EACH GEOMAGNETIC GRID POINT GRDINT.23 C **** WT(4,LM,1) WEIGHTS FOR 4 CORNERS OF EACH GEOMAGNETIC GRDINT.24 C **** GRID ELEMENT GRDINT.25 C **** LM IS FIRST DIMENSION OF ARRAYS AM, IG, JG, WT IN GRDINT.26 C **** CALLING PROGRAM GRDINT.27 C **** (LM.GE.IMAXM) GRDINT.28 C **** IEM IS FLAG INDICATING TYPE OF TRANSFORMATION GRDINT.29 C **** IEM = 1 FOR OUPUT FIELDS AT DIP EQUATOR GRDINT.30 C **** IEM = 0 FOR GLOBAL OUTPUT FIELDS GRDINT.31 C **** NOTE: IG, JG, WT ARE PRODUCED BY PREVIOUS CALL TO GRDINT.32 C **** GRDSTM. GRDINT.33 C **** GRDINT.34 C **** OUTPUT: GRDINT.35 C **** AM(LM,1) 2-DIM SCALAR FIELD TRANSFORMED TO GEOMAGNETIC GRDINT.36 C **** LATITUDE/LONGITUDE GRID GRDINT.37 C **** GRDINT.38 DIMENSION AM(LM,1),AG(LG,1),IG(LM,1),JG(LM,1),WT(4,LM,1) GRDINT.39 C **** GRDINT.40 C **** CARRY OUT INTERPOLATION GRDINT.41 C **** GRDINT.42 DO IM=1,IMAXM GRDINT.43 AM(IM,1)= GRDINT.44 1 AG(IG(IM,JM),JG(IM,JM))*WT(1,IM,JM)+ GRDINT.45 2 AG(IG(IM,JM)+1,JG(IM,JM))*WT(2,IM,JM)+ GRDINT.46 3 AG(IG(IM,JM)+1,JG(IM,JM)+1)*WT(3,IM,JM)+ GRDINT.47 4 AG(IG(IM,JM),JG(IM,JM)+1)*WT(4,IM,JM) GRDINT.48 ENDDO GRDINT.49 RETURN GRDINT.50 END GRDINT.51 C GRDINT.52