#include "dims.h" SUBROUTINE TRANSF use cons_module,only: kmax,kmaxp1,pi implicit none SAVE C **** C **** CALCULATE QUANTITIES NEEDED FOR FIELD-LINE INTEGRATIONS C **** AND TRANSFORM TO GEOMAGNETIC COORDINATE SYSTEM. C **** #include "params.h" C **** C **** COMMON BLOCKS NEEDED: C **** /CONSTS/, /DYNAMO/, /FIELD/ C **** #include "consdyn.h" #include "dynamo.h" #include "field.h" #include "cterp.h" #include "coefm.h" #include "trig.h" #include "dynphi.h" C **** /WORK/ C **** C **** FIELDS TO BE TRANSFORMED TO GEOMAGNETIC SPACE C **** EXTRAPOLATED DOWN TO 90 KM C **** C **** C************************************************************* real :: ssigma1,ssigma2,zz,ww,adotv,axv,vxb,bxa,adota,a1dta2, | sini,sigma1m,sigma2m,zm,wm,adotvm,axvm,vxbm,bxam,adotam, | a1a2m,siniam,bmodm,azm,pm,sindm,cosdm,ram,aam,cosiam,csthdam, | rtadram,rrm,sinidm,cosidm,costhdm,rtramrm,tint1,tint2,tint3, | rr0 COMMON SSIGMA1(IMAXGP,0:JMAXGP,ZKBM:KMAX), C COMMON/NONAME/ SSIGMA1(IMAXGP,0:JMAXGP,ZKBM:KMAX), C************************************************************* 1 SSIGMA2(IMAXGP,0:JMAXGP,ZKBM:KMAX), 2 ZZ(IMAXGP,0:JMAXGP,ZKBM:kmaxp1),WW(IMAXGP,0:JMAXGP,ZKBM:KMAX), 3 ADOTV(IMAXGP,0:JMAXGP,ZKBM:KMAX,2), 4 AXV(IMAXGP,0:JMAXGP,ZKBM:KMAX,2),VXB(IMAXGP,0:JMAXGP,ZKBM:KMAX), 5 BXA(IMAXGP,0:JMAXGP,2),ADOTA(IMAXGP,0:JMAXGP,2), 6 A1DTA2(IMAXGP,0:JMAXGP),SINI(IMAXGP,0:JMAXGP), C **** C **** THESE SAME QUANTITIES IN GEOMAGNETIC SPACE C **** (ONE LATITUDE AT A TIME) C **** 1 SIGMA1M(IMAXMP,ZKBM:KMAX),SIGMA2M(IMAXMP,ZKBM:KMAX), 2 ZM(IMAXMP,ZKBM:kmaxp1),WM(IMAXMP,ZKBM:KMAX), 3 ADOTVM(IMAXMP,ZKBM:KMAX,2),AXVM(IMAXMP,ZKBM:KMAX,2), 4 VXBM(IMAXMP,ZKBM:KMAX),BXAM(IMAXMP,2),ADOTAM(IMAXMP,2), 5 A1A2M(IMAXMP),SINIAM(IMAXMP),BMODM(IMAXMP),AZM(IMAXMP,2), 6 PM(IMAXMP), C **** C **** QUANTITIES EVALUATED IN GEOMAGNETIC SPACE C **** C **** 2-D FIELDS C **** 1 SINDM(IMAXMP),COSDM(IMAXMP),RAM(IMAXMP),AAM(IMAXMP), 2 COSIAM(IMAXMP),CSTHDAM(IMAXMP),RTADRAM(IMAXMP), C **** C **** 3-D FIELDS C **** 1 RRM(IMAXMP,ZKBM:kmaxp1),SINIDM(IMAXMP,ZKBM:kmaxp1), 2 COSIDM(IMAXMP,ZKBM:kmaxp1),COSTHDM(IMAXMP,ZKBM:kmaxp1), 3 RTRAMRM(IMAXMP,ZKBM:kmaxp1), C **** C **** WORK SPACE C **** 1 TINT1(JMAXM),TINT2(JMAXM), 2 TINT3(JMAXM) C **** C **** SET CONSTANTS C **** RL1,RL2 ARE RATES AT WHICH SIGMA1 AND SIGMA2 DECAY C **** WITH HEIGHT BELOW BOTTOM OF MODEL C **** Z0 IS LOWEST Z LEVEL FOR START OF FIELD LINE C **** INTEGRATIONS C **** ! ! Local: real :: rl1=5.E5, rl2=3.E5, z0=30.E5 integer :: i,k,j,kk,n,jj,je real,external :: sddot C **** C **** POLAR POINTS OF MAGNETIC QUANTITIES TO BE TRAMSFORMED C **** TO GEOMAGNETIC SPACE C **** BMOD(1,0) = (9.*sddot(IMAXG,UNIT,BMOD(1,1))- 1 sddot(IMAXG,UNIT,BMOD(1,2)))/(8.*FLOAT(IMAXG)) BMOD(1,JMAXGP) = (9.*sddot(IMAXG,UNIT,BMOD(1,JMAXG))- 1 sddot(IMAXG,UNIT,BMOD(1,JMAXG-1)))/(8.*FLOAT(IMAXG)) P(1,0) = (9.*sddot(IMAXG,UNIT,P(1,1))- 1 sddot(IMAXG,UNIT,P(1,2)))/(8.*FLOAT(IMAXG)) P(1,JMAXGP) = (9.*sddot(IMAXG,UNIT,P(1,JMAXG))- 1 sddot(IMAXG,UNIT,P(1,JMAXG-1)))/(8.*FLOAT(IMAXG)) AV(1,0,3,1) = (9.*sddot(IMAXG,UNIT,AV(1,1,3,1))- 1 sddot(IMAXG,UNIT,AV(1,2,3,1)))/(8.*FLOAT(IMAXG)) AV(1,JMAXGP,3,1) = (9.*sddot(IMAXG,UNIT,AV(1,JMAXG,3,1))- 1 sddot(IMAXG,UNIT,AV(1,JMAXG-1,3,1)))/(8.*FLOAT(IMAXG)) AV(1,0,3,2) = (9.*sddot(IMAXG,UNIT,AV(1,1,3,2))- 1 sddot(IMAXG,UNIT,AV(1,2,3,2)))/(8.*FLOAT(IMAXG)) AV(1,JMAXGP,3,2) = (9.*sddot(IMAXG,UNIT,AV(1,JMAXG,3,2))- 1 sddot(IMAXG,UNIT,AV(1,JMAXG-1,3,2)))/(8.*FLOAT(IMAXG)) DO 35 I = 2,IMAXGP BMOD(I,0) = BMOD(1,0) BMOD(I,JMAXGP) = BMOD(1,JMAXGP) P(I,0) = P(1,0) P(I,JMAXGP) = P(1,JMAXGP) AV(I,0,3,1) = AV(1,0,3,1) AV(I,JMAXGP,3,1) = AV(1,JMAXGP,3,1) AV(I,0,3,2) = AV(1,0,3,2) AV(I,JMAXGP,3,2) = AV(1,JMAXGP,3,2) 35 CONTINUE C **** C **** CALCULATE QUANTITIES TO BE TRANSFORMED TO GEOMAGNETIC C **** SPACE C **** C **** 3-D FIELDS C **** DO 1 K = ZKBM,KMAX DO 1 J = 1,JMAXG DO 1 I = 1,IMAXGP SSIGMA1(I,J,K) = SIGMA1(I,J,K) SSIGMA2(I,J,K) = SIGMA2(I,J,K) ZZ(I,J,K) = Z(I,J,K) WW(I,J,K) = W(I,J,K) 1 CONTINUE DO 2 J = 1,JMAXG DO 2 I = 1,IMAXGP ZZ(I,J,kmaxp1) = Z(I,J,kmaxp1) 2 CONTINUE C **** C **** ADOTV = AV(X)*U + AV(Y)*V C **** V X A = U*AV(Y) - V*AV(X) C **** DO 30 N = 1,2 DO 30 K = ZKBM,KMAX DO 30 J = 1,JMAXG DO 30 I = 1,IMAXGP ADOTV(I,J,K,N) = AV(I,J,1,N)*U(I,J,K)+ 1 AV(I,J,2,N)*V(I,J,K) AXV(I,J,K,N) = AV(I,J,1,N)*V(I,J,K)-AV(I,J,2,N)*U(I,J,K) 30 CONTINUE C **** C **** 2-D FIELDS C **** C **** TINT1 = SIN(D), TINT2 = COS(D) C **** VXB = U*COS(D)-V*SIN(D) C **** BXA = B(X)*AV(Y) - B(Y)*AV(X) C **** ADOTA = AV(X)**2 + AV(Y)**2 C **** DO 31 J = 1,JMAXG DO 32 I = 1,IMAXGP TINT1(I) = SQRT(XB(I,J)**2+YB(I,J)**2) TINT2(I) = XB(I,J)/TINT1(I) TINT1(I) = YB(I,J)/TINT1(I) 32 CONTINUE DO 34 K = ZKBM,KMAX KK = K IF(KK.LT.1)KK = 1 DO 34 I = 1,IMAXGP VXB(I,J,K) = U(I,J,KK)*TINT2(I)-V(I,J,KK)*TINT1(I) 34 CONTINUE DO 31 N = 1,2 DO 31 I = 1,IMAXGP BXA(I,J,N) = TINT1(I)*AV(I,J,2,N)-TINT2(I)*AV(I,J,1,N) ADOTA(I,J,N) = AV(I,J,1,N)**2+AV(I,J,2,N)**2 31 CONTINUE C **** C **** A1DTA2 = AV(X,1)*AV(X,2) + AV(Y,1)*AV(Y,2) C **** SINI = ZB/BMOD C **** DO 33 J = 1,JMAXG DO 33 I = 1,IMAXGP A1DTA2(I,J) = AV(I,J,1,1)*AV(I,J,1,2)+ 1 AV(I,J,2,1)*AV(I,J,2,2) SINI(I,J) = ZB(I,J)/BMOD(I,J) 33 CONTINUE C **** C **** CALCULATE VALUES OF THESE FIELDS AT GEOGRAPHIC POLES C **** DO 5 K = ZKBM,(KMAX-ZKBM+1)*9+ZKBM+6 SSIGMA1(1,0,K) = (9.*sddot(IMAXG,UNIT,SSIGMA1(1,1,K))- 1 sddot(IMAXG,UNIT,SSIGMA1(1,2,K)))/(8.*FLOAT(IMAXG)) SSIGMA1(1,JMAXGP,K) = (9.*sddot(IMAXG,UNIT, 1 SSIGMA1(1,JMAXG,K))-sddot(IMAXG,UNIT, 2 SSIGMA1(1,JMAXG-1,K)))/(8.*FLOAT(IMAXG)) DO 5 I = 2,IMAXG SSIGMA1(I,0,K) = SSIGMA1(1,0,K) SSIGMA1(I,JMAXGP,K) = SSIGMA1(1,JMAXGP,K) 5 CONTINUE C **** C **** PERIODIC POINTS C **** DO 6 J = 0,(9*(KMAX-ZKBM+1)+7)*(JMAXG+2)-1 SSIGMA1(IMAXGP,J,ZKBM) = SSIGMA1(1,J,ZKBM) 6 CONTINUE C **** C **** TRANSFORM NEEDED FIELDS TO GEOMAGNETIC COORDINATE SYSTEM C **** ONE MAGNETIC LATITUDE AT A TIME C **** DO 7 J = 2,JMAXM-1 IF(J.EQ.JMAXM/2+1)GO TO 7 C **** C **** TGCM FIELDS C **** DO 8 K = ZKBM,(KMAX-ZKBM+1)*9+ZKBM+6 CALL GRDINT(SIGMA1M(1,K),SSIGMA1(1,0,K),IG,JG,WT,IMAXGP, 1 IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J) 8 CONTINUE C **** C **** STORE ZM IN 3-D ARRAY ZZM(IMAXMP,JMAXM,kmaxp1) C **** FOR USE IN THREED C **** DO 50 K = ZKBM,KMAX+1 DO 50 I = 1,IMAXM ZZM(I,J,K) = ZM(I,K) 50 CONTINUE C **** C **** MAGNETIC QUANTITIES C **** CALL GRDINT(BMODM,BMOD(1,0),IG,JG,WT,IMAXGP, 1 IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J) CALL GRDINT(PM,P(1,0),IG,JG,WT,IMAXGP, 1 IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J) DO 9 N = 1,2 CALL GRDINT(AZM(1,N),AV(1,0,3,N),IG,JG,WT,IMAXGP, 1 IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J) 9 CONTINUE C **** C **** CALL INTGRLS TO PERFORM FIELD LINE INTEGRATIONS C **** AND EVALUATE PDE COEFFICIENTS AND RHS FOR LATITUDE C **** JM C **** CALL INTGRLS(J) 7 CONTINUE C **** C **** EQUATORIAL VALUES OF ZM C **** DO 51 K = ZKBM,KMAX+1 DO 51 I = 1,IMAXM ZZM(I,JMAXMH,K) = .5*(ZZM(I,JMAXMH-1,K)+ZZM(I,JMAXMH+1,K)) 51 CONTINUE C **** C **** COMPUTE POLAR VALUES OF ZM C **** DO 52 K = ZKBM,KMAX+1 ZZM(1,1,K) = (4.*sddot(IMAXM,UNIT,ZZM(1,2,K))- 1 sddot(IMAXM,UNIT,ZZM(1,3,K)))/(3.*FLOAT(IMAXM)) ZZM(1,JMAXM,K) = (4.*sddot(IMAXM,UNIT,ZZM(1,JMAXM-1,K))- 1 sddot(IMAXM,UNIT,ZZM(1,JMAXM-2,K)))/(3.*FLOAT(IMAXM)) 52 CONTINUE C **** C **** EXTEND OVER LONGITUDE C **** DO 53 K = ZKBM,KMAX+1 DO 53 I = 2,IMAXM ZZM(I,1,K) = ZZM(1,1,K) ZZM(I,JMAXM,K) = ZZM(1,JMAXM,K) 53 CONTINUE C **** C **** PERIODIC POINTS C **** DO 54 K = ZKBM,KMAX+1 DO 54 J = 1,JMAXM ZZM(IMAXMP,J,K) = ZZM(1,J,K) 54 CONTINUE C **** C **** FOLD SOUTHERN HEMISPHERE OVER ON TO NORTHERN C **** DO 36 J = 2,JMAXMH-1 DO 36 I = 1,IMAXM ZIGM11(I,JMAXM+1-J) = ZIGM11(I,JMAXM+1-J)+ZIGM11(I,J) ZIGMC(I,JMAXM+1-J) = ZIGMC(I,JMAXM+1-J)+ZIGMC(I,J) ZIGM2(I,JMAXM+1-J) = ZIGM2(I,JMAXM+1-J)+ZIGM2(I,J) ZIGM22(I,JMAXM+1-J) = ZIGM22(I,JMAXM+1-J)+ZIGM22(I,J) RIM(I,JMAXM+1-J,1) = RIM(I,JMAXM+1-J,1)+RIM(I,J,1) RIM(I,JMAXM+1-J,2) = RIM(I,JMAXM+1-J,2)+RIM(I,J,2) 36 CONTINUE C **** C **** POLE C **** ZIGM11(1,JMAXM) = (4.*sddot(IMAXM,UNIT,ZIGM11(1,JMAXM-1))- 1 sddot(IMAXM,UNIT,ZIGM11(1,JMAXM-2)))/(3.*FLOAT(IMAXM)) ZIGMC(1,JMAXM) = (4.*sddot(IMAXM,UNIT,ZIGMC(1,JMAXM-1))- 1 sddot(IMAXM,UNIT,ZIGMC(1,JMAXM-2)))/(3.*FLOAT(IMAXM)) ZIGM2(1,JMAXM) = (4.*sddot(IMAXM,UNIT,ZIGM2(1,JMAXM-1))- 1 sddot(IMAXM,UNIT,ZIGM2(1,JMAXM-2)))/(3.*FLOAT(IMAXM)) ZIGM22(1,JMAXM) = (4.*sddot(IMAXM,UNIT,ZIGM22(1,JMAXM-1))- 1 sddot(IMAXM,UNIT,ZIGM22(1,JMAXM-2)))/(3.*FLOAT(IMAXM)) C RIM(1,JMAXM,1) = (4.*sddot(IMAXM,UNIT,RIM(1,JMAXM-1,1))- C 1 sddot(IMAXM,UNIT,RIM(1,JMAXM-2,1)))/(3.*FLOAT(IMAXM)) C RIM(1,JMAXM,2) = (4.*sddot(IMAXM,UNIT,RIM(1,JMAXM-1,2))- C 1 sddot(IMAXM,UNIT,RIM(1,JMAXM-2,2)))/(3.*FLOAT(IMAXM)) C **** C **** VECTOR (I1,I2) C **** CDIR$ IVDEP DO 47 I = 1,IMAXM RIM(I,JMAXM,1) = .5*(RIM(I,JMAXM-1,1)-RIM(1+MOD(I-1+IMAXM/2, 1 IMAXM),JMAXM-1,1)) RIM(I,JMAXM,2) = .5*(RIM(I,JMAXM-1,2)-RIM(1+MOD(I-1+IMAXM/2, 1 IMAXM),JMAXM-1,2)) 47 CONTINUE C **** C **** EXTEND OVER LONGITUDE C **** DO 10 I = 2,IMAXM ZIGM11(I,JMAXM) = ZIGM11(1,JMAXM) ZIGMC(I,JMAXM) = ZIGMC(1,JMAXM) ZIGM2(I,JMAXM) = ZIGM2(1,JMAXM) ZIGM22(I,JMAXM) = ZIGM22(1,JMAXM) C RIM(I,JMAXM,1) = RIM(1,JMAXM,1) C RIM(I,JMAXM,2) = RIM(1,JMAXM,2) 10 CONTINUE C **** C **** EQUATOR C **** DO 11 I =1,IMAXM ZIGM11(I,JMAXMH) = 0. ZIGMC(I,JMAXMH) = 0. ZIGM2(I,JMAXMH) = 0. ZIGM22(I,JMAXMH) = 0. RIM(I,JMAXMH,1) = 0. RIM(I,JMAXMH,2) = 0. 11 CONTINUE C **** C **** PERIODIC POINT C **** DO 12 J = 1,6*JMAXM ZIGM11(IMAXMP,J) = ZIGM11(1,J) 12 CONTINUE C **** C **** CALCULATE RHS OF PDE FROM RIM(1) AND RIM(2) C **** C **** FILL ARRAY TINT1 WITH COS(THETA0) C **** DO 37 J = 1,JMAXM TINT1(J) = COS(-PI/2.+(J-1)*DLATM) 37 CONTINUE C **** C **** TRANSFER RIM(1) TO TINT3 IN /WORK/ C **** JE = JMAXMH DO 45 J = 2,JE-1 JJ = J+JE-1 DO 38 I = 1,IMAXM TINT3(I) =RIM(I,JJ,1) 38 CONTINUE C **** C **** PERIODIC POINTS FOR TINT2 C **** DO 39 I =1,2 TINT3(I-2) = TINT3(I-2+IMAXM) TINT3(I+IMAXM) = TINT3(I) 39 CONTINUE C **** C **** PERFORM DIFFERENTIATION OF RIM(1) W.R.T. LAMDA C **** DO 40 I = 1,IMAXM RHS(I,J) = 1./(DLONM*TINT1(JE+J-1))*.5*(TINT3(I+1)- 1 TINT3(I-1)) 40 CONTINUE 45 CONTINUE C **** C **** PERFORM DIFFERENTIATION OF RIM(2) W.R.T. THETA0 C **** DO 41 J = JE+1,JMAXM-1 JJ = J-JE+1 DO 41 I = 1,IMAXM RHS(I,JJ) = RHS(I,JJ)+1./(DLATM*TINT1(J))* 1 .5*(RIM(I,J+1,2)*TINT1(J+1)- 2 RIM(I,J-1,2)*TINT1(J-1)) 41 CONTINUE C **** C **** NOW DEAL WITH JMAXM C **** RHS(1,JMAXMH) = -2./FLOAT(IMAXM)*sddot(IMAXM,UNIT, 1 RIM(1,JMAXM-1,2))/TINT1(JMAXM-1) C **** C **** AND EQUATOR C **** DO 48 I = 1,IMAXM RHS(I,1) = 1./DLATM*.5*TINT1(JMAXMH+1)*RIM(I,JMAXMH+1,2) 48 CONTINUE C **** C **** EXTEND OVER LONGITUDE C **** DO 43 I = 2,IMAXM RHS(I,JMAXMH) = RHS(1,JMAXMH) 43 CONTINUE C **** C **** PERIODIC POINTS C **** DO 44 J = 1,JMAXMH RHS(IMAXMP,J) = RHS(1,J) 44 CONTINUE C **** C **** SCALE (MULTIPLY BY R0*1.E-8 FOR RESULTS IN VOLTS) C **** DO 46 J = 1,JMAXMH DO 46 I = 1,IMAXMP RHS(I,J) = RHS(I,J)*R0*1.E-8 46 CONTINUE RETURN END C