SUBROUTINE TRANSF implicit none C **** C **** CALCULATE QUANTITIES NEEDED FOR FIELD-LINE INTEGRATIONS C **** AND TRANSFORM TO GEOMAGNETIC COORDINATE SYSTEM. C **** include "params.h" include "consts.h" include "dynamo.h" include "fieldz.h" include "cterp.h" include "coefm.h" include "trig.h" include "dynphi.h" include "transmag.h" ! ! Local: integer :: i,k,j,n,kk,je,jj real :: rl1,rl2,z0 ! TINTx workspace must be in common to allow "backward" reference ! of TINT3 in DO 39 below: real :: TINT1(JMAXM),TINT2(JMAXM),TINT3(JMAXM) common/tint/ tint1,tint2,tint3 ! ! Externals: real,external :: sddot ! in util.F 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 **** C **** DATA RL1,RL2,Z0/5.E5,3.E5,90.E5/ 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 = 1,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,KMAXP) = Z(I,J,KMAXP) 2 CONTINUE C **** C **** EXTEND FIELDS DOWN TO 90KM INSERTING 3 EXTRA LEVELS C **** SET THREE EQUALLY SPACED LEVELS FOR Z C **** TAKE U, V AND W TO BE CONSTANT C **** EXTRAPOLATE SIGMAS EXPONENTIALLY C **** DO 3 K = -2,0 DO 3 J = 1,JMAXG DO 3 I = 1,IMAXGP ZZ(I,J,K) = Z0+FLOAT(K+2)*(ZZ(I,J,1)-Z0)/3. WW(I,J,K) = WW(I,J,1) 3 CONTINUE DO 4 K = -2,0 DO 4 J = 1,JMAXG DO 4 I = 1,IMAXGP SSIGMA1(I,J,K) = SSIGMA1(I,J,1)*EXP((ZZ(I,J,K)+ZZ(I,J,K+1)- 1 ZZ(I,J,1)-ZZ(I,J,2))/(2.*RL1)) SSIGMA2(I,J,K) = SSIGMA2(I,J,1)*EXP((ZZ(I,J,K)+ZZ(I,J,K+1)- 1 ZZ(I,J,1)-ZZ(I,J,2))/(2.*RL2)) 4 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 = -2,KMAX KK = K IF(KK.LT.1)KK = 1 DO 30 J = 1,JMAXG DO 30 I = 1,IMAXGP ADOTV(I,J,K,N) = AV(I,J,1,N)*U(I,J,KK)+ 1 AV(I,J,2,N)*V(I,J,KK) AXV(I,J,K,N) = AV(I,J,1,N)*V(I,J,KK)-AV(I,J,2,N)*U(I,J,KK) 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 = -2,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) = ZZB(I,J)/BMOD(I,J) 33 CONTINUE C **** C **** CALCULATE VALUES OF THESE FIELDS AT GEOGRAPHIC POLES C **** DO 5 K = -2,(KMAX+3)*9+4 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+3)+7)*(JMAXG+2)-1 SSIGMA1(IMAXGP,J,-2) = SSIGMA1(1,J,-2) 6 CONTINUE C **** C **** TRANSFORM NEEDED FIELDS TO GEOMAGNETIC COORDINATE SYSTEM C **** ONE MAGNETIC LATITUDE AT A TIME C **** !$OMP PARALLEL DO PRIVATE(j,k,n,i) DO 7 J = 2,JMAXM-1 ! IF(J.EQ.JMAXM/2+1)GO TO 7 if (j.ne.jmaxm/2+1) then ! write(6,"('transf: j=',i2)") j C **** C **** TGCM FIELDS C **** ! DO 8 K = -2,(KMAX+3)*9+4 ! CALL GRDINT(SIGMA1M(1,K),SSIGMA1(1,0,K),IG,JG,WT,IMAXGP, ! 1 IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J) ! 8 CONTINUE ! do k=-2,kmax call grdint(sigma1m(1,k),ssigma1(1,0,k),ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j) call grdint(sigma2m(1,k),ssigma2(1,0,k),ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j) call grdint(wm(1,k),ww(1,0,k),ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j) do n=1,2 call grdint(adotvm(1,k,n),adotv(1,0,k,n),ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j) call grdint(axvm(1,k,n),axv(1,0,k,n),ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j) enddo call grdint(vxbm(1,k),vxb(1,0,k),ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j) enddo ! ! Save sigma1m on magnetic grid: call addfsech('SIGMA1M','SIGMA1M',' ', | sigma1m,imaxmp,kmaxp+3,kmaxp+3,j) call addfsech('SIGMA2M','SIGMA2M',' ', | sigma2m,imaxmp,kmaxp+3,kmaxp+3,j) do k=1,2 call grdint(bxam(1,k),bxa(1,0,k),ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j) call grdint(adotam(1,k),adota(1,0,k),ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j) enddo do k=-2,kmaxp call grdint(zm(1,k),zz(1,0,k),ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j) enddo call grdint(a1a2m,a1dta2,ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j) call grdint(siniam,sini,ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j) C **** C **** STORE ZM IN 3-D ARRAY ZZM(IMAXMP,JMAXM,-2:KMAXP) C **** FOR USE IN THREED C **** DO 50 K = -2,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) endif 7 CONTINUE C **** C **** EQUATORIAL VALUES OF ZM C **** DO 51 K = -2,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 = -2,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 = -2,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 = -2,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 **** 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 **** PI = 4.*ATAN(1.) 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 **** ! jmaxm=97, jmaxmh=(jmaxm+1)/2=49, imaxm=80, tint3(jmaxm) JE = JMAXMH DO 45 J = 2,JE-1 ! 2,48 JJ = J+JE-1 ! 49,96 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 ! ! Save ZZM on magnetic grid to secondary histories: do j=1,jmaxm call addfsech('ZM','GEOPOTENTIAL HEIGHT (MAGNETIC)',' ', | zzm(:,j,:),imaxmp,kmaxp+3,kmaxp+3,j) enddo RETURN END C