! ! file from Ben Foster tgcm14 and include changes ! to the dynamo from Art Richmond ! #include "dims.h" SUBROUTINE TRANSF use cons_module,only: pi implicit none ! calculate quantities needed for field-line integrations ! and transform to geomagnetic coordinate system #include "params.h" #include "consdyn.h" #include "dynamo.h" #include "fieldz.h" #include "cterp.h" #include "coefm.h" #include "trig.h" #include "dynphi.h" #include "transmag.h" ! integer,parameter :: MSGUN=6, MLAT=37,MLON=74,MALT=2, NGRF=8, | LWK= MLAT*MLON*MALT*5 + MLAT+MLON+MALT real :: sinalat,r0or,rat,ddd real :: clm2(imaxgp,jmaxg) ! real :: dvec,dddarr,be3arr,gplat,gplon,gpalt,wk,corfac common/dvecddd/dvec(IMAXGP,JMAXG,3,2),dddarr(IMAXGP,JMAXG), |be3arr(IMAXGP,JMAXG),GPLAT(MLAT),GPLON(MLON),GPALT(MALT),WK(LWK) ! ! equatorial quantities real :: SIGMA1E,SIGMA2E,ADOTVE common/equat/SIGMA1E(IMAXMP,-2:KMAX),SIGMA2E(IMAXMP,-2:KMAX), 5 ADOTVE(IMAXMP,-2:KMAX,2) ! 3d arrays for calculation of height-integrated sheet-current density K_(q,lam) ! used only if icalkqlam = 1 common/ns3d/ sigma1m3d(imaxmp,jmaxm,-2:zkmxp), | sigma2m3d(imaxmp,jmaxm,-2:zkmxp), | bmodm3d(imaxmp,jmaxm), | adotv1m3d(imaxmp,jmaxm,-2:zkmxp), | adotv2m3d(imaxmp,jmaxm,-2:zkmxp), | a1a2m3d(imaxmp,jmaxm), | adotam3d(imaxmp,jmaxm), | sinim3d(imaxmp,jmaxm), | ed13d(imaxmp,jmaxm,-2:zkmxp),ed23d(imaxmp,jmaxm,-2:zkmxp) real :: sigma1m3d,sigma2m3d,bmodm3d,adotv1m3d,adotv2m3d, | a1a2m3d,adotam3d,sinim3d,ed13d,ed23d ! flag for calculation of K_(q,lam) common/kqlam/icalkqlam integer :: icalkqlam ! ! Local: integer :: i,k,j,n,kk,je,jj,ii,isolve,ijmr real :: rl1,rl2,z0 ! ! Diagnostics for secondary histories: real,dimension(zimxp,zkmxp) :: sigma1_plt,sigma2_plt,ww_plt, | zz_plt, ped_plt, hall_plt, z_plt, h_plt, u_plt, v_plt, w_plt, | adotv1_plt,adotv2_plt real,dimension(imaxmp,kmaxp+3) :: | zigm11_plt, zigmc_plt, zigm2_plt, zigm22_plt, rim1_plt, | rim2_plt, rhs_plt ! ! Externals: real,external :: sddot ! in util.F ! flag if height-integrated sheet-current density should be calculated ! icalkqlam = 1 ! will be calculated ! icalkqlam = 0 ! won't be calculated icalkqlam = 1 ! will be calculated ! set constants ! RL1,RL2 are rates at which SIGMA1 and SIGMA2 decay with height ! below bottom of model ! Z0 is lowest level for start of field line integration set in H0 consdyn.h DATA RL1,RL2/5.E5,3.E5/ z0 = h0 ! ! Plot input fields (dynamo.h): ! do j=1,zjmx ! do k=1,zkmx ! do i=1,zimx+1 ! ped_plt(i+2,k) = sigma1(i,j,k) ! 3-75 <= 1-73 ! hall_plt(i+2,k) = sigma2(i,j,k) ! z_plt(i+2,k) = z(i,j,k) ! h_plt(i+2,k) = h(i,j,k) ! u_plt(i+2,k) = u(i,j,k) ! v_plt(i+2,k) = v(i,j,k) ! w_plt(i+2,k) = w(i,j,k) ! enddo ! enddo ! do i=1,zimx+1 ! z_plt(i+2,zkmxp) = z(i,j,zkmxp) ! enddo ! call addfsech('SIGMAPED',' ',' ',ped_plt ,zimxp,zkmxp,zkmx,j) ! call addfsech('SIGMAHAL',' ',' ',hall_plt,zimxp,zkmxp,zkmx,j) ! call addfsech('ZPOTEN' ,' ',' ',z_plt ,zimxp,zkmxp,zkmxp,j) ! call addfsech('SCHEIGHT',' ',' ',h_plt ,zimxp,zkmxp,zkmx,j) ! call addfsech('UNVEL' ,' ',' ',u_plt ,zimxp,zkmxp,zkmx,j) ! call addfsech('VNVEL' ,' ',' ',v_plt ,zimxp,zkmxp,zkmx,j) ! call addfsech('WNVEL' ,' ',' ',w_plt ,zimxp,zkmxp,zkmx,j) ! enddo ! DO 35 J = 1,JMAXG DO 35 I = 1,IMAXGP sinalat = sin(ALATM(I,J)) ! sin(lam) clm2(I,J) = 1. - sinalat*sinalat ! cos^2(lam) SINI(I,J) = ZZB(I,J)/BMOD(I,J) ! sin(I_m) zzb defined in fieldz.h tgcm14 be3(I,J) = 1.e-9*be3arr(I,J) ! be3 is in T (be3arr in nT) 35 CONTINUE ! calculate quantities to be transformed to geomagnetic space ! 3D fields 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) 1 CONTINUE ! DO 2 J = 1,JMAXG DO 2 I = 1,IMAXGP ZZ(I,J,KMAXP) = Z(I,J,KMAXP) 2 CONTINUE ! extend fields down to 90km inserting 3 extra levels ! set three equally spaced levels for Z ! take U,V and W to be constant ! extrapolate sigmas exponentially DO 4 K = 0,-2,-1 DO 4 J = 1,JMAXG DO 4 I = 1,IMAXGP ZZ(I,J,K) = Z0+FLOAT(K+2)*(ZZ(I,J,1)-Z0)/3. 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 ! calculation of ADOTV = ue1,ue2 (m/s) DO 30 K = -2,KMAX KK = K IF(KK.LT.1)KK = 1 DO 30 J = 1,JMAXG DO 30 I = 1,IMAXGP r0or = R0/(R0 + .5*(ZZ(I,J,K)+ZZ(I,J,K+1)) - Z0) ! d_1 = (R_0/R)^1.5 rat = 1.e-2*r0or**1.5 ! 1/100 convertion in cm ADOTV(I,J,K,1) = rat*( ! A_1 dot V = fac( d_1(1) u 1 dvec(I,J,1,1)*U(I,J,KK) + ! + d_1(2) v + d_1(3) w 2 dvec(I,J,2,1)*V(I,J,KK) + 3 dvec(I,J,3,1)*W(I,J,KK)) ! rat = rat*sqrt((4.-3.*clm2(I,J))/(4.-3.*r0or*clm2(I,J))) ! Note: clm2 is being used here to represent the squared cosine of the ! quasi-dipole latitude, not of the M(90) latitude, since the wind ! values are aligned vertically, not along the field line. ! ADOTV(I,J,K,2) = rat*( ! A_2 dot V = fac( d_2(1) u 1 dvec(I,J,1,2)*U(I,J,KK) + ! + d_2(2) v + d_2(3) w 2 dvec(I,J,2,2)*V(I,J,KK) + 3 dvec(I,J,3,2)*W(I,J,KK)) 30 CONTINUE ! calculation of ADOTA(n) = d(n)**2/D ! A1DTA2 = (d(1) dot d(2)) /D DO 33 J = 1,JMAXG DO 33 I = 1,IMAXGP ddd = dddarr(I,J) ADOTA(I,J,1) = (dvec(I,J,1,1)**2 + dvec(I,J,2,1)**2 1 + dvec(I,J,3,1)**2)/ddd ADOTA(I,J,2) = (dvec(I,J,1,2)**2 + dvec(I,J,2,2)**2 1 + dvec(I,J,3,2)**2)/ddd A1DTA2(I,J) = (dvec(I,J,1,1)*dvec(I,J,1,2)+ 1 dvec(I,J,2,1)*dvec(I,J,2,2)+ 2 dvec(I,J,3,1)*dvec(I,J,3,2))/ddd 33 CONTINUE ! ! values at the poles ! loop over different arrays splitted ! original code from Art Richmond: ZZ(IMAXGP,0:JMAXGP,-2:KMAXP), ! SSIGMA1(IMAXGP,0:JMAXGP,-2:KMAX),SSIGMA2(IMAXGP,0:JMAXGP,-2:KMAX), ! ADOTV(IMAXGP,0:JMAXGP,-2:KMAX,2),ADOTA(IMAXGP,0:JMAXGP,2), ! A1DTA2(IMAXGP,0:JMAXGP),SINI(IMAXGP,0:JMAXGP), ! be3(IMAXGP,0:JMAXGP), ! DO 5 K = -2,KMAX ZZ(1,0,K) = (9.*sddot(IMAXG,UNIT,ZZ(1,1,K))- 1 sddot(IMAXG,UNIT,ZZ(1,2,K)))/ 2 (8.*FLOAT(IMAXG)) ZZ(1,JMAXGP,K) = 1 (9.*sddot(IMAXG,UNIT,ZZ(1,JMAXG,K))- 1 sddot(IMAXG,UNIT,ZZ(1,JMAXG-1,K)))/ 2 (8.*FLOAT(IMAXG)) SSIGMA1(1,0,K) = (9.*sddot(IMAXG,UNIT,SSIGMA1(1,1,K))- 1 sddot(IMAXG,UNIT,SSIGMA1(1,2,K)))/ 2 (8.*FLOAT(IMAXG)) SSIGMA1(1,JMAXGP,K) = (9.* 1 sddot(IMAXG,UNIT,SSIGMA1(1,JMAXG,K))- 1 sddot(IMAXG,UNIT,SSIGMA1(1,JMAXG-1,K)))/ 2 (8.*FLOAT(IMAXG)) SSIGMA2(1,0,K) = (9.*sddot(IMAXG,UNIT,SSIGMA2(1,1,K))- 1 sddot(IMAXG,UNIT,SSIGMA2(1,2,K)))/ 2 (8.*FLOAT(IMAXG)) SSIGMA2(1,JMAXGP,K) = (9.* 1 sddot(IMAXG,UNIT,SSIGMA2(1,JMAXG,K))- 1 sddot(IMAXG,UNIT,SSIGMA2(1,JMAXG-1,K)))/ 2 (8.*FLOAT(IMAXG)) ! calculate velocity distribution also at the poles ! used later for calculating K_{q lam} ! ADOTV(1,0,K,1) = (9.*ADOTV(1,1,K,1) - ADOTV(1,2,K,1))/8. ADOTV(1,JMAXGP,K,1) = (9.*ADOTV(1,JMAXG,K,1) - 1 ADOTV(1,JMAXG-1,K,1))/8. ADOTV(1,0,K,2) = (9.*ADOTV(1,1,K,2) - ADOTV(1,2,K,2))/8. ADOTV(1,JMAXGP,K,2) = (9.*ADOTV(1,JMAXG,K,2) - 1 ADOTV(1,JMAXG-1,K,2))/8. ! extend in longitude DO 5 I = 2,IMAXG ZZ(I,0,K) = ZZ(1,0,K) ZZ(I,JMAXGP,K) = ZZ(1,JMAXGP,K) SSIGMA1(I,0,K) = SSIGMA1(1,0,K) SSIGMA1(I,JMAXGP,K) = SSIGMA1(1,JMAXGP,K) SSIGMA2(I,0,K) = SSIGMA2(1,0,K) SSIGMA2(I,JMAXGP,K) = SSIGMA2(1,JMAXGP,K) ! ADOTV(I,0,K,1) = (9.*ADOTV(I,1,K,1) - 1 ADOTV(I,2,K,1))/8. ADOTV(I,JMAXGP,K,1) = (9.*ADOTV(I,JMAXG,K,1)- 1 ADOTV(I,JMAXG-1,K,1))/8. ADOTV(I,0,K,2) = (9.*ADOTV(I,1,K,2) - 1 ADOTV(I,2,K,2))/8. ADOTV(I,JMAXGP,K,2) = (9.*ADOTV(I,JMAXG,K,2) - 1 ADOTV(I,JMAXG-1,K,2))/8. 5 CONTINUE ! values at the poles ZZ(1,0,KMAXP)= (9.*sddot(IMAXG,UNIT,ZZ(1,1,KMAXP))- 1 sddot(IMAXG,UNIT,ZZ(1,2,KMAXP)))/ 2 (8.*FLOAT(IMAXG)) ZZ(1,JMAXGP,KMAXP) = (9.* 1 sddot(IMAXG,UNIT,ZZ(1,JMAXG,KMAXP))- 1 sddot(IMAXG,UNIT,ZZ(1,JMAXG-1,KMAXP)))/ 2 (8.*FLOAT(IMAXG)) ADOTA(1,0,1) = (9.*sddot(IMAXG,UNIT,ADOTA(1,1,1))- 1 sddot(IMAXG,UNIT,ADOTA(1,2,1)))/ 2 (8.*FLOAT(IMAXG)) ADOTA(1,JMAXGP,1) = (9.* 1 sddot(IMAXG,UNIT,ADOTA(1,JMAXG,1))- 1 sddot(IMAXG,UNIT,ADOTA(1,JMAXG-1,1)))/ 2 (8.*FLOAT(IMAXG)) ADOTA(1,0,2) = (9.*sddot(IMAXG,UNIT,ADOTA(1,1,2))- 1 sddot(IMAXG,UNIT,ADOTA(1,2,2)))/ 2 (8.*FLOAT(IMAXG)) ADOTA(1,JMAXGP,2) = (9.* 1 sddot(IMAXG,UNIT,ADOTA(1,JMAXG,2))- 1 sddot(IMAXG,UNIT,ADOTA(1,JMAXG-1,2)))/ 2 (8.*FLOAT(IMAXG)) A1DTA2(1,0) = (9.*sddot(IMAXG,UNIT,A1DTA2(1,1))- 1 sddot(IMAXG,UNIT,A1DTA2(1,2)))/ 2 (8.*FLOAT(IMAXG)) A1DTA2(1,JMAXGP) = (9.* 1 sddot(IMAXG,UNIT,A1DTA2(1,JMAXG))- 1 sddot(IMAXG,UNIT,A1DTA2(1,JMAXG-1)))/ 2 (8.*FLOAT(IMAXG)) SINI(1,0) = (9.*sddot(IMAXG,UNIT,SINI(1,1))- 1 sddot(IMAXG,UNIT,SINI(1,2)))/ 2 (8.*FLOAT(IMAXG)) SINI(1,JMAXGP) = (9.* 1 sddot(IMAXG,UNIT,SINI(1,JMAXG))- 1 sddot(IMAXG,UNIT,SINI(1,JMAXG-1)))/ 2 (8.*FLOAT(IMAXG)) BE3(1,0) = (9.*sddot(IMAXG,UNIT,BE3(1,1))- 1 sddot(IMAXG,UNIT,BE3(1,2)))/ 2 (8.*FLOAT(IMAXG)) BE3(1,JMAXGP) = (9.* 1 sddot(IMAXG,UNIT,BE3(1,JMAXG))- 1 sddot(IMAXG,UNIT,BE3(1,JMAXG-1)))/ 2 (8.*FLOAT(IMAXG)) ! extend in longitude DO 55 I = 2,IMAXG ZZ(I,0,KMAXP) = ZZ(1,0,KMAXP) ZZ(I,JMAXGP,KMAXP)= ZZ(1,JMAXGP,KMAXP) ADOTA(I,0,1) = ADOTA(1,0,1) ADOTA(I,JMAXGP,1) = ADOTA(1,JMAXGP,1) ADOTA(I,0,2) = ADOTA(1,0,2) ADOTA(I,JMAXGP,2) = ADOTA(1,JMAXGP,2) A1DTA2(I,0) = A1DTA2(1,0) A1DTA2(I,JMAXGP) = A1DTA2(1,JMAXGP) SINI(I,0) = SINI(1,0) SINI(I,JMAXGP) = SINI(1,JMAXGP) BE3(I,0) = BE3(1,0) BE3(I,JMAXGP) = BE3(1,JMAXGP) 55 CONTINUE ! periodic points ! original code from Art Richmond splitted in serveral arrays: ! DO 6 J = 0,JMAXG+1 DO 66 k = -2,KMAX ZZ(IMAXGP,J,k) = ZZ(1,J,k) SSIGMA1(IMAXGP,J,k) = SSIGMA1(1,J,k) SSIGMA2(IMAXGP,J,k) = SSIGMA2(1,J,k) ADOTV(IMAXGP,J,k,1) = ADOTV(1,J,k,1) ADOTV(IMAXGP,J,k,2) = ADOTV(1,J,k,2) 66 CONTINUE ZZ(IMAXGP,J,kmaxp)= ZZ(1,J,kmaxp) ADOTA(IMAXGP,J,1) = ADOTA(1,J,1) ADOTA(IMAXGP,J,2) = ADOTA(1,J,2) A1DTA2(IMAXGP,J) = A1DTA2(1,J) SINI(IMAXGP,J) = SINI(1,J) BE3(IMAXGP,J) = BE3(1,J) 6 CONTINUE ! ! Save geographic grid fields to secondary history: ! real,dimension(zimxp,zkmxp) :: sigma1_plt,sigma2_plt ! COMMON /transmag/ SSIGMA1(IMAXGP,0:JMAXGP,-2:KMAX), ! ! sigma1_plt(:,:) = 0. ! sigma2_plt(:,:) = 0. ! do j=1,zjmx ! do k=1,zkmx ! do i=1,zimx ! sigma1_plt(i,k) = ssigma1(i,j,k) ! sigma2_plt(i,k) = ssigma2(i,j,k) ! zz_plt(i,k) = zz(i,j,k) ! adotv1_plt(i,k) = adotv(i,j,k,1) ! adotv2_plt(i,k) = adotv(i,j,k,2) ! enddo ! i=1,zimx ! enddo ! k=1,zkmx ! call addfsech('SIGMA1',' ',' ',sigma1_plt,zimxp,zkmxp,zkmx,j) ! call addfsech('SIGMA2',' ',' ',sigma2_plt,zimxp,zkmxp,zkmx,j) ! call addfsech('ZZ' ,' ',' ',zz_plt ,zimxp,zkmxp,zkmx,j) ! call addfsech('ADOTV1',' ',' ',adotv1_plt,zimxp,zkmxp,zkmx,j) ! call addfsech('ADOTV2',' ',' ',adotv2_plt,zimxp,zkmxp,zkmx,j) ! enddo ! j=1,zjmx ! transform needed fields to geomagnetic coordinate system ! one magnetic latitude at a time ! original code from Art Richmond splitted in serveral arrays: ! equatorial values: used in INTGRLS for interpolating conductivities and winds ! along the magnetic field line J = JMAXM/2+1 DO 9 K = -2,KMAX CALL GRDINT(SIGMA1E(1,K),SSIGMA1(1,0,K),IG,JG,WT,IMAXGP, 1 IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J,0) SIGMA1E(IMAXMP,K) = SIGMA1E(1,K) CALL GRDINT(SIGMA2E(1,K),SSIGMA2(1,0,K),IG,JG,WT,IMAXGP, 1 IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J,0) SIGMA2E(IMAXMP,K) = SIGMA2E(1,K) CALL GRDINT(ADOTVE(1,K,1),ADOTV(1,0,K,1),IG,JG,WT,IMAXGP, 1 IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J,0) ADOTVE(IMAXMP,K,1) = ADOTVE(1,K,1) CALL GRDINT(ADOTVE(1,K,2),ADOTV(1,0,K,2),IG,JG,WT,IMAXGP, 1 IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM,J,0) ADOTVE(IMAXMP,K,2) = ADOTVE(1,K,2) ! 9 CONTINUE !!$OMP PARALLEL DO PRIVATE(j,k,n,i) ! for test not parallel ! DO 7 J = 2,JMAXM-1 ! field line integration for each latitude ! loop 8 splitted in different arrays do k=-2,kmax call grdint(sigma1m(1,k),ssigma1(1,0,k),ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j,0) call grdint(sigma2m(1,k),ssigma2(1,0,k),ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j,0) 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,0) enddo sigma1m(imaxmp,k) = sigma1m(1,k) ! periodic point lon=180 sigma2m(imaxmp,k) = sigma2m(1,k) ! periodic point lon=180 enddo ! call addfsech('SIGMA1M',' ',' ',sigma1m,imaxmp, ! | kmaxp+3,kmaxp+3,j) ! call addfsech('SIGMA2M',' ',' ',sigma2m,imaxmp, ! | kmaxp+3,kmaxp+3,j) ! cga ! Save sigma1m on magnetic grid (other latitudes): ! ! problem with runs **** commented out begin ! call addfsech('SIGMA1M','SIGMA1M',' ', ! | sigma1m,imaxmp,kmaxp+3,kmax+3,j) ! call addfsech('SIGMA2M','SIGMA2M',' ', ! | sigma2m,imaxmp,kmaxp+3,kmax+3,j) ! (poles): ! if(j.eq.jmaxm-1) then ! call addfsech('SIGMA1M','SIGMA1M',' ', ! | sigma1m,imaxmp,kmaxp+3,kmax+3,jmaxm) ! call addfsech('SIGMA2M','SIGMA2M',' ', ! | sigma2m,imaxmp,kmaxp+3,kmax+3,jmaxm) ! elseif(j.eq.2) then ! call addfsech('SIGMA1M','SIGMA1M',' ', ! | sigma1m,imaxmp,kmaxp+3,kmax+3,1) ! call addfsech('SIGMA2M','SIGMA2M',' ', ! | sigma2m,imaxmp,kmaxp+3,kmax+3,1) ! endif ! problem with runs **** commented out end ! cge do k=1,2 call grdint(adotam(1,k),adota(1,0,k),ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j,0) 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,0) enddo call grdint(a1a2m,a1dta2,ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j,0) call grdint(siniam,sini,ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j,0) call grdint(bmodm,be3,ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j,0) ! store ZM in 3-D array ZMM(IMAXMP,JMAXM,-2:KMAXP) for ! use in THREED if(icalkqlam.eq.1) then ! for calculation of k_(q,lam) DO I = 1,IMAXM ! bmodm3d(i,j) = bmodm(i) a1a2m3d(i,j) = a1a2m(i) adotam3d(i,j) = adotam(i,1) sinim3d(i,j) = siniam(i) DO K = -2,KMAX if(ZM(I,K).lt.z0) ZM(I,K) = Z0 ZZM(I,J,K) = ZM(I,K) ! ! store 3d arrays in magnetic coordinates ! sigma1m,sigma2m,bmodm,adotvm,a1a2m ! for postprocessing of the current density k_(q,lam) ! sigma1m3d(i,j,k) = sigma1m(i,k) sigma2m3d(i,j,k) = sigma2m(i,k) adotv1m3d(i,j,k) = adotvm(i,k,1) adotv2m3d(i,j,k) = adotvm(i,k,2) enddo k = kmax+1 if(ZM(I,K).lt.z0) ZM(I,K) = Z0 ZZM(I,J,K) = ZM(I,K) sigma1m3d(i,j,k) = sigma1m(i,k-1) sigma2m3d(i,j,k) = sigma2m(i,k-1) adotv1m3d(i,j,k) = adotvm(i,k-1,1) adotv2m3d(i,j,k) = adotvm(i,k-1,2) enddo ! periodic point bmodm3d(imaxmp,j) = bmodm3d(1,j) a1a2m3d(imaxmp,j) = a1a2m3d(1,j) adotam3d(imaxmp,j) = adotam3d(1,j) sinim3d(imaxmp,j) = sinim3d(1,j) sigma1m3d(imaxmp,j,:) = sigma1m3d(1,j,:) sigma2m3d(imaxmp,j,:) = sigma2m3d(1,j,:) adotv1m3d(imaxmp,j,:) = adotv1m3d(1,j,:) adotv2m3d(imaxmp,j,:) = adotv2m3d(1,j,:) else ! without calculation of k_(q,lam) DO I = 1,IMAXM DO K = -2,KMAX if(ZM(I,K).lt.z0) ZM(I,K) = Z0 ZZM(I,J,K) = ZM(I,K) enddo k = kmax+1 if(ZM(I,K).lt.z0) ZM(I,K) = Z0 ZZM(I,J,K) = ZM(I,K) enddo endif IF(J.EQ.JMAXM/2+1)GO TO 7 ! skip the equator (numerical problems) ! ! call intgrls to perform field line integrations and evaluate ! PDE coefficients and RHS for latitude J CALL INTGRLS(J) 7 CONTINUE ! next latitude of integration loop ! values for the poles j=1 and j=jmaxm for the 3d arrays if(icalkqlam.eq.1) then ! for calculation of k_(q,lam) do jj =1,2 j = jj if(jj.eq.2) j=jmaxm call grdint(a1a2m,a1dta2,ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j,0) call grdint(siniam,sini,ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j,0) call grdint(bmodm,be3,ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j,0) call grdint(adotam(1,1),adota(1,0,1),ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j,0) do k=-2,kmax call grdint(sigma1m(1,k),ssigma1(1,0,k),ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j,0) call grdint(sigma2m(1,k),ssigma2(1,0,k),ig,jg,wt,imaxgp, | imaxmp,imaxg,jmaxg+2,imaxm,jmaxm,j,0) 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,0) enddo sigma1m(imaxmp,k) = sigma1m(1,k) ! periodic point lon=180 sigma2m(imaxmp,k) = sigma2m(1,k) ! periodic point lon=180 enddo ! enddo k do i = 1,imaxm bmodm3d(i,j) = bmodm(i) a1a2m3d(i,j) = a1a2m(i) adotam3d(i,j) = adotam(i,1) sinim3d(i,j) = siniam(i) do k = -2,kmax sigma1m3d(i,j,k) = sigma1m(i,k) sigma2m3d(i,j,k) = sigma2m(i,k) adotv1m3d(i,j,k) = adotvm(i,k,1) adotv2m3d(i,j,k) = adotvm(i,k,2) enddo k = kmax+1 sigma1m3d(i,j,k) = sigma1m(i,k-1) sigma2m3d(i,j,k) = sigma2m(i,k-1) adotv1m3d(i,j,k) = adotvm(i,k-1,1) adotv2m3d(i,j,k) = adotvm(i,k-1,2) enddo ! enddo i ! bmodm3d(imaxmp,j) = bmodm3d(1,j) a1a2m3d(imaxmp,j) = a1a2m3d(1,j) adotam3d(imaxmp,j) = adotam3d(1,j) sinim3d(imaxmp,j) = sinim3d(1,j) sigma1m3d(imaxmp,j,:) = sigma1m3d(1,j,:) sigma2m3d(imaxmp,j,:) = sigma2m3d(1,j,:) adotv1m3d(imaxmp,j,:) = adotv1m3d(1,j,:) adotv2m3d(imaxmp,j,:) = adotv2m3d(1,j,:) enddo ! enddo jj endif ! ! in intgrls the following values are calculated: ! At this point, ZIGM11 is int[sig_p*d_1^2/D] ds, i.e. Sigma_(phi phi)/abs(sin Im) ! ZIGM22 is int[sig_p*d_2^2/D] ds, i.e. Sigma_(lam lam)*abs(sin Im) ! ZIGMC is int[sig_p*d_1*d_2/D] ds, i.e. Sigma_c ! ZIGM2 is int[sigma_h] ds, i.e. Sigma_h ! ! RIM1 is int[(sigma_h-sigma_p*d_1*d_2/D)u_e1 + sigma_p*d_1^2/D u_e2] *A(h_r)*B_e3 ds, i.e. K_(m phi)^D/abs(sin Im) ! RIM2 is int[(sigma_h+sigma_p*d_1*d_2/D)u_e2 - sigma_p*d_2^2/D u_e1] *A(h_r)*B_e3 ds, K_(m lam)^D ( minus in northern hemisphere ! Change sign of RIM(2) in S. hemisphere to be compatible with transf ! At this point, RIM2 is +-K_(m lam)^D ! Equatorial values: ! Assume that quantities primarily dependent on Pedersen conductivity ! have field-line integrals 1/4 as large as the averages for next-higher ! field lines; quantities primarily dependent on Hall conductivity ! have field-line integrals 0.12 as large as the averages for next-higher ! field lines. Exact values chosen should not be important for potential ! calculation, as long as they are physically reasonable and not too ! different from adjacent values. J = JMAXM/2+1 DO 5055 I = 1,IMAXM ZIGM11(I,J)= .125*(ZIGM11(I,J-1)+ ZIGM11(I,J+1)) ZIGM22(I,J)= .125*(ZIGM22(I,J-1)+ ZIGM22(I,J+1)) ZIGMC(I,J) = .125*(ZIGMC(I,J-1) + ZIGMC(I,J+1)) ZIGM2(I,J) = .06 *(ZIGM2(I,J-1) + ZIGM2(I,J+1)) RIM(I,J,1) = .06 *(RIM(I,J-1,1) + RIM(I,J+1,1)) RIM(I,J,2) = .06 *(RIM(I,J-1,2) + RIM(I,J+1,2)) 5055 CONTINUE ! include the boundary condition at the equator eq.(5.30) in ! Richmond (1995) Ionospheric Electrodynamics use. Mag. Apex Coord. J.Goemag.Geoelectr. 47,191-212 ! Sig_phiphi/abs(sin Im) = 0.5*Sig_cowling/abs(sin Im) ! = 0.5/abs(sin Im)*(Sig_phiphi - Sig_philam*sig_lamphi/Sig_lamlam) ! = 0.5/abs(sin Im)*(Sig_phiphi + (Sig_h-sig_c)*(Sig_h+sig_c)/Sig_lamlam) ! rim(1) / |sin I_m| = I_1 = R/2*(K_mphi - Sig_philam/Sig_lamlam*K_mlam) j = jmaxmh ! equator do i = 1,imaxm zigm11(i,j) = zigm11(i,j)+ (zigm2(i,j)-zigmc(i,j))* | (zigm2(i,j)+zigmc(i,j))/zigm22(i,j) rim(i,j,1) = rim(i,j,1) - (zigm2(i,j)-zigmc(i,j))/ | zigm22(i,j)*rim(i,j,2) zigm11(i,j) = 0.5*zigm11(i,j) rim(i,j,1) = 0.5*rim(i,j,1) enddo ! ! Using notation of Richmond (1995) on right-hand side below: ! Sigma_(phi phi) = ZIGM11*abs(sin I_m) ! Sigma_(lam lam) = ZIGM22/abs(sin I_m) ! Sigma_(phi lam) = +-(ZIGM2-ZIGMC) ! Sigma_(lam phi) = -+(ZIGM2+ZIGMC) ! K_(m phi)^D = RIM(1)*abs(sin I_m) ! K_(m lam)^D = +-RIM(2) ! ! ! transforming PDE from original apex (theta_a) to new apex grid (theta_0) ! which is equally spaced in magnetic latitude ! SCALE quantities to modified (0) magnetic latitude system, ! multiplying or dividing by abs(sin I_m) [inverse contained in DT1DTS] as ! necessary. Sign of K_(m lam)^D in southern hemisphere remains reversed. ! for the mixed terms the transformation from the integration and differentiation ! canceld out (ZIGMC, ZIGM2) ! DT1DTS : d theta_0/ d theta_a / abs(sin I_m) ! RCOS0S : cos(theta_0)/ cos(theta_a) DO 24 J = 2,JMAXM-1 corfac = RCOS0S(J)/DT1DTS(J) ! abs(I_m)*d theta_a/d theta_0 * cos(theta_0)/ cos(theta_a) DO 24 I = 1,IMAXM ZIGM11(I,J) = ZIGM11(I,J)*corfac ! abs(I_m)*d theta_a/d theta_0 * cos(theta_0)/ cos(theta_a) ZIGM22(I,J) = ZIGM22(I,J)/corfac ! 1/abs(I_m)*d theta_0/d theta_a * cos(theta_a)/ cos(theta_0) RIM(I,J,1) = RIM(I,J,1)/DT1DTS(J) ! abs(I_m)*d theta_a/d theta_0 RIM(I,J,2) = RIM(I,J,2)/RCOS0S(J) ! cos(theta_a)/ cos(theta_0) 24 CONTINUE DO 1012 J = 1,6*JMAXM ZIGM11(IMAXMP,J) = ZIGM11(1,J) ! periodic points 1012 CONTINUE ! the PDE is divided by 1/ DT0DTS ! Sigma_(phi phi) = ZIGM11/ RCOS0S * DT0DTS ! Sigma_(lam lam) = ZIGM22 * RCOS0S / DT0DTS ! Sigma_(phi lam) = +-(ZIGM2-ZIGMC) ! Sigma_(lam phi) = -+(ZIGM2+ZIGMC) ! K_(m phi)^D = RIM(1) * DT0DTS ! K_(m lam)^D = +-RIM(2) * RCOS0S ! compute polar values of ZZM (3-D geopotential height in magnetic coordinates) ! 4th order interpolation 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 ! extend over longitude 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 ! periodic points DO 54 K = -2,KMAX+1 DO 54 J = 1,JMAXM ZZM(IMAXMP,J,K) = ZZM(1,J,K) 54 CONTINUE ! compute polar values for the conductances ! 4th order interpolation ZIGM11(1, 1) = (4.*sddot(IMAXM,UNIT,ZIGM11(1, 2))- 1 sddot(IMAXM,UNIT,ZIGM11(1, 3)))/(3.*FLOAT(IMAXM)) ZIGM11(1,JMAXM) = (4.*sddot(IMAXM,UNIT,ZIGM11(1,JMAXM-1))- 1 sddot(IMAXM,UNIT,ZIGM11(1,JMAXM-2)))/(3.*FLOAT(IMAXM)) ZIGMC(1, 1) = (4.*sddot(IMAXM,UNIT,ZIGMC(1, 2))- 1 sddot(IMAXM,UNIT,ZIGMC(1, 3)))/(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, 1) = (4.*sddot(IMAXM,UNIT,ZIGM2(1, 2))- 1 sddot(IMAXM,UNIT,ZIGM2(1, 3)))/(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, 1) = (4.*sddot(IMAXM,UNIT,ZIGM22(1, 2))- 1 sddot(IMAXM,UNIT,ZIGM22(1, 3)))/(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)) ! extend over longitude DO 10 I = 2,IMAXM ZIGM11(I, 1) = ZIGM11(1, 1) ZIGM11(I,JMAXM) = ZIGM11(1,JMAXM) ZIGMC(I, 1) = ZIGMC(1, 1) ZIGMC(I, JMAXM) = ZIGMC(1, JMAXM) ZIGM2(I, 1) = ZIGM2(1, 1) ZIGM2(I, JMAXM) = ZIGM2(1, JMAXM) ZIGM22(I, 1) = ZIGM22(1, 1) ZIGM22(I, JMAXM) = ZIGM22(1,JMAXM) 10 CONTINUE ! RHS vector (I_1,I_2): average over poles DO 47 I = 1,IMAXM RIM(I, 1,1) = .5*(RIM(I, 2,1)-RIM(1+MOD(I-1+IMAXM/2, 1 IMAXM), 2,1)) RIM(I,JMAXM,1) = .5*(RIM(I,JMAXM-1,1)-RIM(1+MOD(I-1+IMAXM/2, 1 IMAXM),JMAXM-1,1)) RIM(I, 1,2) = .5*(RIM(I, 2,2)-RIM(1+MOD(I-1+IMAXM/2, 1 IMAXM), 2,2)) RIM(I,JMAXM,2) = .5*(RIM(I,JMAXM-1,2)-RIM(1+MOD(I-1+IMAXM/2, 1 IMAXM),JMAXM-1,2)) 47 CONTINUE ! periodic points DO 12 J = 1,6*JMAXM ZIGM11(IMAXMP,J) = ZIGM11(1,J) 12 CONTINUE ! ! Save to secondary history: ! do j=1,jmaxm ! do i=1,imaxmp ! zigm11_plt(i,:) = zigm11(i,j) ! zigmc_plt (i,:) = zigmc (i,j) ! zigm2_plt (i,:) = zigm2 (i,j) ! zigm22_plt(i,:) = zigm22(i,j) ! rim1_plt (i,:) = rim (i,j,1) ! rim2_plt (i,:) = rim (i,j,2) ! enddo ! i=1,imaxm ! call addfsech('ZIGM11',' ',' ',zigm11_plt,imaxmp,kmaxp+3,kmaxp,j) ! call addfsech('ZIGMC' ,' ',' ',zigmc_plt ,imaxmp,kmaxp+3,kmaxp,j) ! call addfsech('ZIGM2' ,' ',' ',zigm2_plt ,imaxmp,kmaxp+3,kmaxp,j) ! call addfsech('ZIGM22',' ',' ',zigm22_plt,imaxmp,kmaxp+3,kmaxp,j) ! call addfsech('RIM1' ,' ',' ',rim1_plt ,imaxmp,kmaxp+3,kmaxp,j) ! call addfsech('RIM2' ,' ',' ',rim2_plt ,imaxmp,kmaxp+3,kmaxp,j) ! enddo ! folding south on to northern hemisphere and calculation of RHS ! moved to subroutine dyncal RETURN END C