#include "dims.h" SUBROUTINE DYN use cons_module,only: pi implicit none ! set up coefficients for dynamo PDE #include "params.h" #include "coefm.h" #include "dynamo.h" #include "fieldz.h" #include "cterp.h" #include "trig.h" #include "fcom.h" #include "ceee.h" #include "dynphi.h" ! global: unmodified coefficient for using modified mudpack (isolve=2) real :: cofum common/mudmd/cofum(imx0,jmx0,9) ! flag for calculation of K_(q,lam) set in subroutine transf common/kqlam/icalkqlam integer :: icalkqlam ! local: real :: ARRAY(-15:IMX0+16,JMX0),CS(JMX0) EQUIVALENCE (F(1,1),ARRAY(-15,1)) real :: C0(IMX0,JMX0,10),C1(IMX1,JMX1,9),C2(IMX2,JMX2,9), 1 C3(IMX3,JMX3,9),C4(IMX4,JMX4,9) EQUIVALENCE (CEE,C0),(CEE(NC1),C1),(CEE(NC2),C2),(CEE(NC3),C3), 1 (CEE(NC4),C4) real :: WK(IMX0,3),dlonm,dlatm,sym integer :: j,j0,jj,jjj,i,jntl,k,isolve,n real,dimension(imaxmp,zkmxp+3) :: | zigm11_plt, zigmc_plt, zigm2_plt, zigm22_plt, rim1_plt, | rim2_plt, rhs_plt ! specify solver typ ! isolve = 0 ! org. mud version 5. ! isolve = 1 ! muh hybrid solver (only as direct solver- slow) isolve = 2 ! modified mudpack solver ! transform needed fields to geomagnetic coordinates ! perform field-line integrations ! evaluate PDE coefficients and RHS ! the PDE is divided by 1/ DT0DTS (in dyncal divided by 1/cos(theta_0) ! 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 call transf ! ! calculate coefficient of PDE for each hemisphere if K_(q,lam) is calculated ! if(icalkqlam.eq.1) call nosocoef ! fold southern hemisphere over on to northern (was in transf.F version tgcm15) ! value at the equator is also folded therefore at the equatorial boundary ! condition factor 1/2 introduced ! size of arrays 49,97 from equator to poles (added) ! southern hemisphere 1,49 original values DO 36 J = 1,JMAXMH DO 36 I = 1,IMAXMP 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) ! Reverse sign of ZIGMC to be compatible with Cicely's ZIGMC(I,JMAXM+1-J) = -ZIGMC(I,JMAXM+1-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)) ! sign of K_(m lam)^D in southern hemisphere ! remains reversed; therefore + sign 36 CONTINUE ! calculate RHS of PDE from RIM(1) and RIM(2) (was in transf.F version tgcm15) ! [( d K_(m phi)^D / d phi /(cos(theta_m)?) + ! (d [ K_(m lam)^D * cos(lam_m)]/ d lam_m ) /cos ( lam_m) ] * R / (RCOS0S*DT0DTS) ! ~ J_(Mr)*r^2*cos(theta_m)/cos(theta_0)/DT0DTS call rhspde ! set array NC and magnetic latitude cosine array ! NC pointes to the start of the coefficient array for each level NC(1) = NC0 NC(2) = NC1 NC(3) = NC2 NC(4) = NC3 NC(5) = NC4 NC(6) = NCEE ! total length of coefficient array DLONM = 2.*PI/IMAXM ! Delta phi_m DLATM = PI/(JMAXM-1) ! Delta lam_m DO 10 J = 1,JMX0 CS(J) = COS(PI/2.-(JMX0-J)*DLATM) ! equator to poles (1,49): cos(0), cos(pi/2) 10 CONTINUE ! set up difference coefficients ! replace ZIGM11 by A, ZIGM22 by B, ZIGMC by C and ZIGM2 by D J0 = JMX0-JMAXMH ! jmx0 = jmaxmh (49) s. params.F DO 30 J = 1,JMAXMH ! latitudinal loop 1,49 JJ = JMAXMH+J-1 ! 49,97 JJJ = JMAXMH -J+1 ! 49,1 DO 31 I = 1,IMAXMP ! factor 4 from 5-point diff. stencil ZIGMC(I,JJ) = (ZIGMC(I,JJ)+ZIGM2(I,JJ))/(4.*DLATM*DLONM) ! Sigma_(phi lam)/( 4*Delta lam* Delta lon ) ZIGM2(I,JJ) = ZIGMC(I,JJ)-2.*ZIGM2(I,JJ)/(4.*DLATM*DLONM) ! Sigma_(phi lam)/( 4*Delta lam* Delta lon ) ZIGM22(I,JJ) = ZIGM22(I,JJ)*CS(J0+J)/DLATM**2 ! Sigma_(lam lam)*cos(lam_m)*DT0DTS/(Delta lam)^2 ZIGMC(I,JJJ) = -ZIGMC(I,JJ) ! -ZIGMC_north = southern hemis. 49,1 equator-pole ZIGM2(I,JJJ) = -ZIGM2(I,JJ) ! -ZIGM2_north = southern hemis. 49,1 equator-pole ZIGM22(I,JJJ)= ZIGM22(I,JJ) ! ZIGM22 = southern hemis. 49,1 equator-pole 31 CONTINUE ! ZIGM11 calculated for equator to include equatorial boundary condition DO 32 I = 1,IMAXMP ZIGM11(I,JJ) = ZIGM11(I,JJ)/(CS(J0+J)*DLONM**2) ! Sigma_(phi phi)/( cos(lam_m)*DT0DTS*(Delta lon)^2 ) ZIGM11(I,JJJ) = ZIGM11(I,JJ) ! ZIGM11 = southern hemis. 49,1 equator-pole 32 CONTINUE 30 CONTINUE ! end fo latitudinal loop 1,49 ! set ZIGM11 to zero at magnetic poles to avoid floating point ! exception (values at the poles are not used anyway) DO 34 I = 1,IMAXMP ZIGM11(I,1) = 0.0 ZIGM11(I,JMAXM) = 0.0 34 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,zkmxp+3,zkmxp,j) ! call addfsech('ZIGMC' ,' ',' ',zigmc_plt ,imaxmp,zkmxp+3,zkmxp,j) ! call addfsech('ZIGM2' ,' ',' ',zigm2_plt ,imaxmp,zkmxp+3,zkmxp,j) ! call addfsech('ZIGM22',' ',' ',zigm22_plt,imaxmp,zkmxp+3,zkmxp,j) ! call addfsech('RIM1' ,' ',' ',rim1_plt ,imaxmp,zkmxp+3,zkmxp,j) ! call addfsech('RIM2' ,' ',' ',rim2_plt ,imaxmp,zkmxp+3,zkmxp,j) ! enddo ! clear array for difference stencils at all levels CALL CLEAR(CEE,IMX0,JMX0) if(isolve.eq.2) then ! modified solver (modified and unmodified coefficients) cofum(:,:,:) = 0.0 ! initialize array ! calculate contribution to stencils from each PDE coefficient ARRAY(1,1) = 1. ! Sigma_(phi phi)/( cos(lam_m)*DT0DTS*(Delta lon)^2 ) SYM = 1. CALL STENCMD(ZIGM11,CS(1),IMX0,JMX0,SYM,CEE,cofum,ARRAY(1,1), | 1,WK) ARRAY(1,1) = 4. ! Sigma_(lam lam)*cos(lam_m)*DT0DTS/(Delta lam)^2 SYM = 1. CALL STENCMD(ZIGM22,CS(1),IMX0,JMX0,SYM,CEE,cofum,ARRAY(1,1), | 4,WK) ARRAY(1,1) = 2. ! Sigma_(phi lam)/( 4*Delta lam* Delta lon ) SYM = -1. CALL STENCMD(ZIGMC,CS(1),IMX0,JMX0,SYM,CEE,cofum,ARRAY(1,1), | 2,WK) ARRAY(1,1) = 3. ! Sigma_(lam phi)/( 4*Delta lam* Delta lon ) SYM = -1. CALL STENCMD(ZIGM2,CS(1),IMX0,JMX0,SYM,CEE,cofum,ARRAY(1,1), | 3,WK) else ! original solver (mud) or direct solver ! calculate contribution to stencils from each PDE coefficient ARRAY(1,1) = 1. ! Sigma_(phi phi)/( cos(lam_m)*DT0DTS*(Delta lon)^2 ) SYM = 1. CALL STENCIL(ZIGM11,CS(1),IMX0,JMX0,SYM,CEE,ARRAY(1,1),1,WK) ARRAY(1,1) = 4. ! Sigma_(lam lam)*cos(lam_m)*DT0DTS/(Delta lam)^2 SYM = 1. CALL STENCIL(ZIGM22,CS(1),IMX0,JMX0,SYM,CEE,ARRAY(1,1),4,WK) ARRAY(1,1) = 2. ! Sigma_(phi lam)/( 4*Delta lam* Delta lon ) SYM = -1. CALL STENCIL(ZIGMC,CS(1),IMX0,JMX0,SYM,CEE,ARRAY(1,1),2,WK) ARRAY(1,1) = 3. ! Sigma_(lam phi)/( 4*Delta lam* Delta lon ) SYM = -1. CALL STENCIL(ZIGM2,CS(1),IMX0,JMX0,SYM,CEE,ARRAY(1,1),3,WK) endif ! call RTHS to insert RHS in finest stencil CALL RTHS(RHS,IMX0,JMX0,JMAXMH,C0) ! set boundary condition at the pole CALL EDGES(C0,IMX0,JMX0) CALL EDGES(C1,IMX1,JMX1) CALL EDGES(C2,IMX2,JMX2) CALL EDGES(C3,IMX3,JMX3) CALL EDGES(C4,IMX4,JMX4) if(isolve.eq.2) CALL EDGES(cofum,IMX0,JMX0) ! divide stencil by cos(lam_0) (not RHS) CALL DIVIDE(C0,IMX0,JMX0,IMX0,JMX0,CS(1)) CALL DIVIDE(C1,IMX1,JMX1,IMX0,JMX0,CS(1)) CALL DIVIDE(C2,IMX2,JMX2,IMX0,JMX0,CS(1)) CALL DIVIDE(C3,IMX3,JMX3,IMX0,JMX0,CS(1)) CALL DIVIDE(C4,IMX4,JMX4,IMX0,JMX0,CS(1)) if(isolve.eq.2) call divimd(cofum,IMX0,JMX0,IMX0,JMX0,CS(1)) ! set value of solution to 1. at pole DO 20 I = 1,IMX0 C0(I,JMX0,10) = 1. 20 CONTINUE ! modify stencils and RHS so that Heelis potential ! ID inserted at high latitude if(isolve.eq.2) then CALL ADHEELMD else CALL ADHEEL endif ! do j=1,jmaxm ! write(6,"('dynamo: j=',i3,' rim(1)=',/,(6e12.4))") j,rim(:,j,1) ! write(6,"('dynamo: j=',i3,' rim(2)=',/,(6e12.4))") j,rim(:,j,2) ! enddo ! j=1,nmlat JNTL = 0 ! solving PDE ! solution returned in RIM(1) ! mudpack solver ! isolve = 0 org. mud v5. (modified stencils neq direct solution) ! isolve = 1 muh hybrid solver (no convergence => only as direct solver) ! isolve = 2 modified mud (residual calculated with unmodified stencils ! same solution as with direct solver, if same ! coefficient matrix is used) ! common/mudmd/cofum(imx0,jmx0,9) ! do j=1,jmx0 ! do n=1,9 ! write(6,"('dynamo: j=',i3,' n=',i2,' cofum(:,j,n)=',/, ! | (6e12.4))") j,n,cofum(:,j,n) ! enddo ! enddo ! call print_cee if(isolve.eq.0) then CALL MUD(RIM,JNTL) else if(isolve.eq.1) then call muh(RIM,JNTL) else if(isolve.eq.2) then call mudmod(rim,jntl) else write(6,*) 'dyncal: solver typ',isolve,' not implemented' stop endif ! copy output potential from RIM(1) to PHIM(IMAXMP,JMAXM) DO J = 1,JMAXM DO I = 1,IMAXMP PHIM(I,J) = RIM(I,J,1) enddo enddo ! am_02/02: calculate current for both hemisphere: ! R**2 * J_mr / dt0dts /rcos0s if(icalkqlam.eq.1) call nosocrrt ! call THREED to generate 3-D potential array in geomagnetic coordinates ! result returned in PHIM3D(IMAXMP,JMAXM,-2:KMAXP). ! also transform this field to geographic space in array ! DYNPOT(IMAXGP,0:JMAXGP,KMAXP) ! calculate 3D drift velocities ED13D and ED23D CALL THREED c ! am_02/02: calculate Je_1, K_(q,phi), K_(q,lam) if(icalkqlam.eq.1) call nosocrdens c RETURN END C !----------------------------------------------------------------------------------- subroutine rhspde use cons_module,only: pi implicit none ! ! am_02/02 calculate the RHS of the PDE from RIM(1) and RIM(2) ! this part was in transf.F before- semms to be more logical here #include "params.h" #include "consdyn.h" #include "coefm.h" ! 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 ! local integer :: i,j,jj,je real,dimension(imaxmp,zkmxp+3) :: rhs_plt ! externals: real,external :: sddot ! in util.F ! calculate RHS of PDE from RIM(1) and RIM(2) ! fill array TINT1 with cos(theta_0) DO 37 J = 1,JMAXM ! south to north pole TINT1(J) = COS(-PI/2.+(J-1)*DLATM) ! cos(theta_0) equally spaced 37 CONTINUE ! transfer RIM(1) to TINT3 in /WORK/ ! jmaxm=97, jmaxmh=(jmaxm+1)/2=49, imaxm=80, tint3(jmaxm) JE = JMAXMH DO 45 J = 2,JE-1 ! latitudinal loop: 2,48 from equator to north pole JJ = J+JE-1 ! 50,96 from equator to north pole (added hemispheres) DO 38 I = 1,IMAXM TINT3(I) = RIM(I,JJ,1) ! K_(m phi)^D / DT0DTS 38 CONTINUE ! periodic points for TINT3 DO 39 I =1,2 TINT3(I-2) = TINT3(I-2+IMAXM) TINT3(I+IMAXM) = TINT3(I) 39 CONTINUE ! perform differentiation of RIM(1) w.r.t.lambda ! ( d K_(m phi)^D / d phi )/ (cos (lam_0) / (d lam_0/ d lam_m) = ! ( d K_(m phi)^D / d phi )/ cos(lam_m) / DT0DTS / RCOS0S = ! ( d K_(m phi)^D(0) / d phi )/ cos(theta_0) 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 ! end of latitude loop ! perform differentiation of RIM(2) w.r.t. lam_0 ! +/- (d [ K_(m lam)^D * cos(lam_m)]/ d lam_0 ) /cos ( lam_0) = ! + (d [ K_(m lam)^D * cos(lam_m)]/ d lam_m ) /cos ( lam_m) / (RCOS0S*DT0DTS) = ! +/- (d [ K_(m lam)^D(0) * cos(lam_0)]/ d lam_0 ) /cos ( lam_0) = ! DO 41 J = JE+1,JMAXM-1 ! 50,96 from equator to poles JJ = J-JE+1 ! 2,48 (equator to poles) 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 ! polar value RHS(1,JMAXMH) = -2./FLOAT(IMAXM)*sddot(IMAXM,UNIT, 1 RIM(1,JMAXM-1,2))/TINT1(JMAXM-1) ! equator: ! include the boundary condition at the equator ! rhs(equator)/R = d (K_mphi^DT(0) - sig_philam/sig_lamlam*K_mlam^DT(0)) / d phi_m ! + d (cos lam_0 * K_mlam^DT(0))/ d lam_0 ! from Cicely's notes: ! I_1 = 0.5*(K_(m phi)^DT(0) - Sig_(phi lam)/Sig_(lam lam)*K_(ml am)^DT(0)) ! I_2 = K_(m lam)^DT(0) ! differentiate ! rhs = (I_1(i+1/2,j)-I_1(i-1/2,j))/dlonm + ! (2*cos(lam_0)_(j+1/2)*I_2(i,j+1/2))/dlat_0 ! ! i = 1 rhs(i,1) = 0.5/dlonm*(rim(i+1,je,1)-rim(imaxm,je,1)) rhs(i,1) = rhs(i,1) + 1./dlatm*(tint1(je)*rim(i,je,2)+ | tint1(je+1)*rim(i,je+1,2)) do 48 i = 2,imaxm-1 rhs(i,1) = 0.5/dlonm*(rim(i+1,je,1)-rim(i-1,je,1)) rhs(i,1) = rhs(i,1) + 1./dlatm*(tint1(je)*rim(i,je,2)+ | tint1(je+1)*rim(i,je+1,2)) 48 continue i = imaxm rhs(i,1) = 0.5/dlonm*(rim(1,je,1)-rim(i-1,je,1)) rhs(i,1) = rhs(i,1) + 1./dlatm*(tint1(je)*rim(i,je,2)+ | tint1(je+1)*rim(i,je+1,2)) ! extend over longitude DO 43 I = 2,IMAXM RHS(I,JMAXMH) = RHS(1,JMAXMH) 43 CONTINUE ! periodic points DO 44 J = 1,JMAXMH RHS(IMAXMP,J) = RHS(1,J) 44 CONTINUE ! scale (multiply by earth radius in meter = R0*1.E-2) ![( d K_(m phi)^D / d phi /(cos(theta_m)?) + ! (d [ K_(m lam)^D * cos(lam_m)]/ d lam_m ) /cos ( lam_m) ] * R / (RCOS0S*DT0DTS) ! ~ J_(Mr)*r^2*cos(theta_m)/cos(theta_0)/DT0DTS ! theta_m = theta_s DO 46 J = 1,JMAXMH DO 46 I = 1,IMAXMP RHS(I,J) = RHS(I,J)*R0*1.E-2 ! scale factor Earth radius in meters 46 CONTINUE ! rhs_plt = 0. ! do j=1,jmaxm ! if (j <= jmaxmh) then ! do i=1,imaxmp ! rhs_plt(i,:) = rhs(i,j) ! enddo ! endif ! call addfsech('RHS' ,' ',' ',rhs_plt,imaxmp,zkmxp+3,zkmxp,j) ! enddo ! return end !----------------------------------------------------------------------- subroutine print_cee implicit none #include "params.h" #include "ceee.h" real :: C0(IMX0,JMX0,10),C1(IMX1,JMX1,9),C2(IMX2,JMX2,9), 1 C3(IMX3,JMX3,9),C4(IMX4,JMX4,9) EQUIVALENCE (CEE,C0),(CEE(NC1),C1),(CEE(NC2),C2),(CEE(NC3),C3), 1 (CEE(NC4),C4) integer :: n,j ! do n=1,10 do j=1,jmx0 write(6,"('print_cee: j=',i3,' nmlat0=',i3,' n=',i2, | ' c0(:,j,n)=',/,(6e12.4))") j,jmx0,n,c0(:,j,n) enddo ! j=1,jmx0 enddo ! n=1,10 do n=1,9 do j=1,jmx1 write(6,"('print_cee: j=',i3,' nmlat1=',i3,' n=',i2, | ' c1(:,j,n)=',/,(6e12.4))") j,jmx1,n,c1(:,j,n) enddo ! j=1,jmx1 enddo ! n=1,9 do n=1,9 do j=1,jmx2 write(6,"('print_cee: j=',i3,' nmlat2=',i3,' n=',i2, | ' c2(:,j,n)=',/,(6e12.4))") j,jmx2,n,c2(:,j,n) enddo ! j=1,jmx2 enddo ! n=1,9 do n=1,9 do j=1,jmx3 write(6,"('print_cee: j=',i3,' nmlat3=',i3,' n=',i2, | ' c3(:,j,n)=',/,(6e12.4))") j,jmx3,n,c3(:,j,n) enddo ! j=1,jmx3 enddo ! n=1,9 do n=1,9 do j=1,jmx4 write(6,"('print_cee: j=',i3,' nmlat4=',i3,' n=',i2, | ' c4(:,j,n)=',/,(6e12.4))") j,jmx4,n,c4(:,j,n) enddo ! j=1,jmx4 enddo ! n=1,9 end subroutine print_cee