C --- ------------------------------------------------------------------ SUBROUTINE T3STEP C --- time step executive routine, implement circular buffer shell game C --- Each of the routines called operates on a spherical grid shell, C --- i.e. for all (K,J) at the given I. C --- Results are stored into circular buffers of spherical grid shells. C --- The buffer pointers (indices, M?, L?) are rotated instead of data. C --- In the called routines, the following naming convention applies: C --- MB <--> I-1 (shell Below current) C --- MX <--> I (current shell, X marks the spot) C --- MA <--> I+1 (shell Above current) C --- LB <--> I-1/2 or I-1 ( for FI,SI and UP respectively) C --- LA <--> I+1/2 COMMON /BLOC4/ IM,JM,KM,IC1,IC2,JC1,JC2,KC1,KC2 I = 1 L1 = 1 M1 = 1 CALL FLUXES(I, M1) M2 = 2 CALL FLUXES(I+1, M2) CALL STEPFM(I,M1,M2, L1) I = 2 L2 = 2 M3 = 3 CALL FLUXES(I+1, M3) CALL STEPFM(I,M2,M3, L2) CALL STEPGM(I,M1,M2,M3 ) CALL STEPHM(I,M1,M2,M3 ) CALL STEPUP(I,L1,L2,M1,M2,M3 ) IF( IC2 .LT. 3 ) GO TO 11 DO 10 I= 3,IC2 LT = L1 L1 = L2 L2 = LT MT = M3 M3 = M1 M1 = M2 M2 = MT CALL FLUXES(I+1, M3) CALL STEPFM(I,M2,M3, L2) CALL STEPGM(I,M1,M2,M3 ) CALL STEPHM(I,M1,M2,M3 ) CALL STEPUP(I,L1,L2,M1,M2,M3 ) CALL STOREP(I-1, L2) 10 CONTINUE 11 CONTINUE CALL STOREP(IC2, L1) RETURN END C --- ------------------------------------------------------------------ SUBROUTINE FLUXES(I, M) C --- compute {UFGHS} from shell I primitives C --- into circular 2-d buffers {UFGHS}(M) INCLUDE 'igm_grid.ins' PARAMETER ( PI1 = 3.141593, PI4 = 4.*PI1, PI8 = 8.*PI1 ) COMMON /BLOC7/ GMS,RS,GAM,GAM1,GMM1,R,PI,XMU0,XMU1,DMU0,DMU1 COMMON /BUFFS/ & U(8,KX,JX,3), & F(8,KX,JX,3), & G(8,KX,JX,3), & H(8,KX,JX,3), & S(8,KX,JX,3), & FI(8,KX,JX,2), & SI(8,KX,JX,2), & GJ(8,KX,JX), & SJ(8,KX,JX), & HK(8,KX,JX), & SK(8,KX,JX), & UP(8,KX,JX,2) DO 20 J=JC1-1,JC2+1 STHM = SIN( TH(J) ) COTHM = 1./TAN( TH(J) ) DO 10 K=KC1-1,KC2+1 R1 = RAD(I) R2 = R1**2 ROP = RO(K,J,I) VRP = VR(K,J,I) VTP = VT(K,J,I) VPP = VP(K,J,I) BRP = BR(K,J,I) BTP = BT(K,J,I) BPP = BP(K,J,I) TPP = TP(K,J,I) PRP = 2.*ROP*R*TPP U(1,K,J,M) = R2*ROP U(2,K,J,M) = R2*ROP*VRP U(3,K,J,M) = R2*ROP*VTP U(4,K,J,M) = R2*ROP*VPP U(5,K,J,M) = R1*BRP U(6,K,J,M) = R1*BTP U(7,K,J,M) = R1*BPP TEMP1 = 2.*ROP*R*TPP/GAM1 TEMP2 = (VRP**2+VTP**2+VPP**2)/2. TEMP3 = (BRP**2+BTP**2+BPP**2)*DMU1 U(8,K,J,M) = R2*(TEMP1+ROP*TEMP2+TEMP3) F(1,K,J,M) = R2*ROP*VRP TEMP2 = ROP*VRP**2 TEMP3 = (-BRP**2+BTP**2+BPP**2)*DMU1 F(2,K,J,M) = R2*(PRP+TEMP2+TEMP3) F(3,K,J,M) = R2*(ROP*VRP*VTP-BRP*BTP*DMU0) F(4,K,J,M) = R2*(ROP*VRP*VPP-BRP*BPP*DMU0) F(5,K,J,M) = 0. F(6,K,J,M) = R1*(VRP*BTP-VTP*BRP) F(7,K,J,M) = R1*(VRP*BPP-VPP*BRP) VKE = 0.5*ROP*( VRP**2 + VTP**2 + VPP**2 ) GPK = GMM1*PRP + VKE TEMP2 = (BTP*DMU0)*(VRP*BTP-VTP*BRP) TEMP3 = -(BPP*DMU0)*(VPP*BRP-VRP*BPP) F(8,K,J,M) = R2*(VRP*GPK + TEMP2 + TEMP3) G(1,K,J,M) = R2*STHM*ROP*VTP G(2,K,J,M) = R2*STHM*(ROP*VRP*VTP-BRP*BTP*DMU0) TEMP2 = ROP*VTP**2 TEMP3 = (BRP**2-BTP**2+BPP**2)*DMU1 G(3,K,J,M) = R2*STHM*(PRP + TEMP2 + TEMP3) G(4,K,J,M) = R2*STHM*(ROP*VTP*VPP-BTP*BPP*DMU0) G(5,K,J,M) = R1*STHM*(VTP*BRP-VRP*BTP) G(6,K,J,M) = 0. G(7,K,J,M) = R1*STHM*(VTP*BPP-VPP*BTP) TEMP2 = -(BRP*DMU0)*(VRP*BTP-VTP*BRP) TEMP3 = (BPP*DMU0)*(VTP*BPP-VPP*BTP) G(8,K,J,M) = R2*STHM*(VTP*GPK + TEMP2 + TEMP3) H(1,K,J,M) = R2*ROP*VPP H(2,K,J,M) = R2*(ROP*VRP*VPP-BRP*BPP*DMU0) H(3,K,J,M) = R2*(ROP*VTP*VPP-BTP*BPP*DMU0) TEMP2 = ROP*VPP**2 TEMP3 = (BRP**2+BTP**2-BPP**2)*DMU1 H(4,K,J,M) = R2*(PRP + TEMP2 + TEMP3) H(5,K,J,M) = R1*(VPP*BRP-VRP*BPP) H(6,K,J,M) = R1*(VPP*BTP-VTP*BPP) H(7,K,J,M) = 0. TEMP2 = BRP*(VPP*BRP-VRP*BPP)*DMU0 TEMP3 = -BTP*(VTP*BPP-VPP*BTP)*DMU0 H(8,K,J,M) = R2*(VPP*GPK+TEMP2+TEMP3) S(1,K,J,M) = 0. TEMP1 = 2.*R1*PRP TEMP2 = -ROP*GMS + R1*BRP**2*DMU0 TEMP3 = ROP*R1*(VTP**2+VPP**2) S(2,K,J,M) = TEMP1+TEMP2+TEMP3 TEMP1 = R1*COTHM*PRP TEMP2 = (BRP**2+BTP**2-BPP**2)*R1*COTHM*DMU1 TEMP3 = R1*BRP*BTP*DMU0-ROP*R1*(VRP*VTP-VPP**2*COTHM) S(3,K,J,M) = TEMP1+TEMP2+TEMP3 TEMP1 = BTP*BPP*R1*COTHM*DMU0 TEMP2 = R1*BRP*BPP*DMU0 TEMP3 = -ROP*R1*(VRP*VPP+VTP*VPP*COTHM) S(4,K,J,M) = TEMP1+TEMP2+TEMP3 S(5,K,J,M) = 0. S(6,K,J,M) = 0. S(7,K,J,M) = (VTP*BPP-VPP*BTP)*COTHM S(8,K,J,M) = -ROP*VRP*GMS 10 CONTINUE 20 CONTINUE RETURN END C --- ------------------------------------------------------------------ SUBROUTINE STEPFM(I,MX,MA, L) C --- do Euler prediction of F & S at I+1/2 -> {FI,SI}(L) C --- using c-buffs {UFGHS}(MX & MA) INCLUDE 'igm_grid.ins' COMMON /BLOC7/ GMS,RS,GAM,GAM1,GMM1,R,PI,XMU0,XMU1,DMU0,DMU1 COMMON /BUFFS/ & U(8,KX,JX,3), & F(8,KX,JX,3), & G(8,KX,JX,3), & H(8,KX,JX,3), & S(8,KX,JX,3), & FI(8,KX,JX,2), & SI(8,KX,JX,2), & GJ(8,KX,JX), & SJ(8,KX,JX), & HK(8,KX,JX), & SK(8,KX,JX), & UP(8,KX,JX,2) REAL UT(8) RX = RAD(I) RA = RAD(I+1) R1 = ( RX + RA )/2. R2 = R1**2 FAC1 = DT/DR DO 20 J=JC1,JC2 FAC2 = DT/( SIN( TH(J) )*4.0*DTH ) FAC3 = DT/( SIN( TH(J) )*4.0*DPH ) COTHM = 1./TAN( TH(J) ) DO 10 K=KC1,KC2 DO 1 N=1,4 UT(N) = 0.5*( U(N,K,J,MA) + U(N,K,J,MX) ) & - FAC1*( F(N,K,J,MA) - F(N,K,J,MX) ) & - FAC2*( ( G(N,K,J+1,MX) - G(N,K,J-1,MX) )/RX & + ( G(N,K,J+1,MA) - G(N,K,J-1,MA) )/RA ) & - FAC3*( ( H(N,K+1,J,MX) - H(N,K-1,J,MX) )/RX & + ( H(N,K+1,J,MA) - H(N,K-1,J,MA) )/RA ) & + DT/4.*( S(N,K,J,MA) + S(N,K,J,MX) ) 1 CONTINUE DO 2 N=5,8 UT(N) = 0.5*( U(N,K,J,MA) + U(N,K,J,MX) ) & - FAC1*( F(N,K,J,MA) - F(N,K,J,MX) ) & - FAC2*( ( G(N,K,J+1,MX) - G(N,K,J-1,MX) )/RX & + ( G(N,K,J+1,MA) - G(N,K,J-1,MA) )/RA ) & - FAC3*( ( H(N,K+1,J,MX) - H(N,K-1,J,MX) )/RX & + ( H(N,K+1,J,MA) - H(N,K-1,J,MA) )/RA ) & + DT/4.*( S(N,K,J,MA) + S(N,K,J,MX) ) 2 CONTINUE ROP = UT(1)/R2 VRP = UT(2)/UT(1) VTP = UT(3)/UT(1) VPP = UT(4)/UT(1) BRP = UT(5)/R1 BTP = UT(6)/R1 BPP = UT(7)/R1 TEMP1 = VRP**2+VTP**2+VPP**2 TEMP2 = (BRP**2+BTP**2+BPP**2)*DMU1 TPP = GAM1/(2.*ROP*R)*(UT(8)/R2-0.5*ROP*TEMP1-TEMP2) PRP = 2.*ROP*R*TPP FI(1,K,J,L) = R2*ROP*VRP TEMP2 = ROP*VRP**2 TEMP3 = (-BRP**2+BTP**2+BPP**2)*DMU1 FI(2,K,J,L) = R2*(PRP+TEMP2+TEMP3) FI(3,K,J,L) = R2*(ROP*VRP*VTP-BRP*BTP*DMU0) FI(4,K,J,L) = R2*(ROP*VRP*VPP-BRP*BPP*DMU0) FI(5,K,J,L) = 0. FI(6,K,J,L) = R1*(VRP*BTP-VTP*BRP) FI(7,K,J,L) = R1*(VRP*BPP-VPP*BRP) VKE = 0.5*ROP*( VRP**2 + VTP**2 + VPP**2 ) GPK = GMM1*PRP + VKE TEMP2 = (BTP*DMU0)*(VRP*BTP-VTP*BRP) TEMP3 = -(BPP*DMU0)*(VPP*BRP-VRP*BPP) FI(8,K,J,L) = R2*(VRP*GPK + TEMP2 + TEMP3) SI(1,K,J,L) = 0. TEMP1 = 2.*R1*PRP TEMP2 = -ROP*GMS + R1*BRP**2*DMU0 TEMP3 = ROP*R1*(VTP**2+VPP**2) SI(2,K,J,L) = TEMP1+TEMP2+TEMP3 TEMP1 = R1*COTHM*PRP TEMP2 = (BRP**2+BTP**2-BPP**2)*R1*COTHM*DMU1 TEMP3 = R1*BRP*BTP*DMU0-ROP*R1*(VRP*VTP-VPP**2*COTHM) SI(3,K,J,L) = TEMP1+TEMP2+TEMP3 TEMP1 = BTP*BPP*R1*COTHM*DMU0 TEMP2 = R1*BRP*BPP*DMU0 TEMP3 = -ROP*R1*(VRP*VPP+VTP*VPP*COTHM) SI(4,K,J,L) = TEMP1+TEMP2+TEMP3 SI(5,K,J,L) = 0. SI(6,K,J,L) = 0. SI(7,K,J,L) = (VTP*BPP-VPP*BTP)*COTHM SI(8,K,J,L) = -ROP*VRP*GMS 10 CONTINUE 20 CONTINUE RETURN END C --- ------------------------------------------------------------------ SUBROUTINE STEPGM(I,MB,MX,MA ) C --- do Euler prediction of U for shell I into 2-d work array UJ C --- UJ is U at point (K, J+1/2, I) at time T + DT C --- UJ is computed from c-buffs {UFGHS}(K,J,{MB,MX,MA}) C --- then compute {GJ,SJ}(K,J), the flux and source vectors INCLUDE 'igm_grid.ins' PARAMETER ( PI1 = 3.141593, PI4 = 4.*PI1, PI8 = 8.*PI1 ) COMMON /BLOC7/ GMS,RS,GAM,GAM1,GMM1,R,PI,XMU0,XMU1,DMU0,DMU1 COMMON /BUFFS/ & U(8,KX,JX,3), & F(8,KX,JX,3), & G(8,KX,JX,3), & H(8,KX,JX,3), & S(8,KX,JX,3), & FI(8,KX,JX,2), & SI(8,KX,JX,2), & GJ(8,KX,JX), & SJ(8,KX,JX), & HK(8,KX,JX), & SK(8,KX,JX), & UP(8,KX,JX,2) REAL UT(8) R1 = RAD(I) R2 = R1**2 FAC1 = DT/(4.*DR) FAC3 = DT/(R1*4.*DPH) DO 20 J=JC1-1,JC2 STJ = SIN( TH(J) ) STJ1 = SIN( TH(J+1) ) THM = ( TH(J) + TH(J+1) )/2. FAC2 = DT/(R1*DTH*0.5*(STJ+STJ1)) STHM = SIN( THM ) COTHM = 1./TAN( THM ) DO 10 K=KC1,KC2 DO 1 N=1,4 UT(N) = 0.5*( U(N,K,J,MX) + U(N,K,J+1,MX) ) & -FAC1*( F(N,K,J,MA) - F(N,K,J,MB) & + F(N,K,J+1,MA) - F(N,K,J+1,MB) ) & -FAC2*( G(N,K,J+1,MX) - G(N,K,J,MX) ) & -FAC3*( ( H(N,K+1,J,MX) - H(N,K-1,J,MX) )/STJ & + ( H(N,K+1,J+1,MX) - H(N,K-1,J+1,MX) )/STJ1 ) & +DT/2.*( S(N,K,J,MX) + S(N,K,J+1,MX) ) 1 CONTINUE DO 2 N=5,8 UT(N) = 0.5*( U(N,K,J,MX) + U(N,K,J+1,MX) ) & -FAC1*( F(N,K,J,MA) - F(N,K,J,MB) & + F(N,K,J+1,MA) - F(N,K,J+1,MB) ) & -FAC2*( G(N,K,J+1,MX) - G(N,K,J,MX) ) & -FAC3*( ( H(N,K+1,J,MX) - H(N,K-1,J,MX) )/STJ & + ( H(N,K+1,J+1,MX) - H(N,K-1,J+1,MX) )/STJ1 ) & +DT/2.*( S(N,K,J,MX) + S(N,K,J+1,MX) ) 2 CONTINUE ROP = UT(1)/R2 VRP = UT(2)/UT(1) VTP = UT(3)/UT(1) VPP = UT(4)/UT(1) BRP = UT(5)/R1 BTP = UT(6)/R1 BPP = UT(7)/R1 TEMP1 = VRP**2+VTP**2+VPP**2 TEMP2 = (BRP**2+BTP**2+BPP**2)*DMU1 TPP = GAM1/(2.*ROP*R)*(UT(8)/R2-0.5*ROP*TEMP1-TEMP2) PRP = 2.*ROP*R*TPP GJ(1,K,J) = R2*STHM*ROP*VTP GJ(2,K,J) = R2*STHM*(ROP*VRP*VTP-BRP*BTP*DMU0) TEMP2 = ROP*VTP**2 TEMP3 = (BRP**2-BTP**2+BPP**2)*DMU1 GJ(3,K,J) = R2*STHM*(PRP + TEMP2 + TEMP3) GJ(4,K,J) = R2*STHM*(ROP*VTP*VPP-BTP*BPP*DMU0) GJ(5,K,J) = R1*STHM*(VTP*BRP-VRP*BTP) GJ(6,K,J) = 0. GJ(7,K,J) = R1*STHM*(VTP*BPP-VPP*BTP) VKE = 0.5*ROP*( VRP**2 + VTP**2 + VPP**2 ) GPK = GMM1*PRP + VKE TEMP2 = -(BRP*DMU0)*(VRP*BTP-VTP*BRP) TEMP3 = (BPP*DMU0)*(VTP*BPP-VPP*BTP) GJ(8,K,J) = R2*STHM*(VTP*GPK + TEMP2 + TEMP3) SJ(1,K,J) = 0. TEMP1 = 2.*R1*PRP TEMP2 = -ROP*GMS + R1*BRP**2*DMU0 TEMP3 = ROP*R1*(VTP**2+VPP**2) SJ(2,K,J) = TEMP1+TEMP2+TEMP3 TEMP1 = R1*COTHM*PRP TEMP2 = (BRP**2+BTP**2-BPP**2)*R1*COTHM*DMU1 TEMP3 = R1*BRP*BTP*DMU0-ROP*R1*(VRP*VTP-VPP**2*COTHM) SJ(3,K,J) = TEMP1+TEMP2+TEMP3 TEMP1 = BTP*BPP*R1*COTHM*DMU0 TEMP2 = R1*BRP*BPP*DMU0 TEMP3 = -ROP*R1*(VRP*VPP+VTP*VPP*COTHM) SJ(4,K,J) = TEMP1+TEMP2+TEMP3 SJ(5,K,J) = 0. SJ(6,K,J) = 0. SJ(7,K,J) = (VTP*BPP-VPP*BTP)*COTHM SJ(8,K,J) = -ROP*VRP*GMS 10 CONTINUE 20 CONTINUE RETURN END C --- ------------------------------------------------------------------ SUBROUTINE STEPHM(I,MB,MX,MA ) C --- do Euler prediction of U for shell I into 2-d work array UK C --- UK is U at point (K+1/2, J, I) at time T + DT C --- UK is computed from c-buffs {UFGHS}(K,J,{MB,MX,MA}) C --- then compute {GK,SK}(K,J), the flux and source vectors INCLUDE 'igm_grid.ins' PARAMETER ( PI1 = 3.141593, PI4 = 4.*PI1, PI8 = 8.*PI1 ) COMMON /BLOC7/ GMS,RS,GAM,GAM1,GMM1,R,PI,XMU0,XMU1,DMU0,DMU1 COMMON /BUFFS/ & U(8,KX,JX,3), & F(8,KX,JX,3), & G(8,KX,JX,3), & H(8,KX,JX,3), & S(8,KX,JX,3), & FI(8,KX,JX,2), & SI(8,KX,JX,2), & GJ(8,KX,JX), & SJ(8,KX,JX), & HK(8,KX,JX), & SK(8,KX,JX), & UP(8,KX,JX,2) REAL UT(8) R1 = RAD(I) R2 = R1**2 FAC1 = DT/(4.*DR) DO 20 J=JC1,JC2 STJ = SIN( TH(J) ) FAC2 = DT/(R1*STJ*4.*DTH) FAC3 = DT/(R1*STJ*DPH) COTHM = 1./TAN( TH(J) ) DO 10 K=KC1-1,KC2 DO 1 N=1,4 UT(N) = 0.5*( U(N,K,J,MX) + U(N,K+1,J,MX) ) & -FAC1*( F(N,K,J,MA) - F(N,K,J,MB) & + F(N,K+1,J,MA) - F(N,K+1,J,MB) ) & -FAC2*( G(N,K,J+1,MX) - G(N,K,J-1,MX) & + G(N,K+1,J+1,MX) - G(N,K+1,J-1,MX) ) & -FAC3*( H(N,K+1,J,MX) - H(N,K,J,MX) ) & +DT/2.*( S(N,K,J,MX) + S(N,K+1,J,MX) ) 1 CONTINUE DO 2 N=5,8 UT(N) = 0.5*( U(N,K,J,MX) + U(N,K+1,J,MX) ) & -FAC1*( F(N,K,J,MA) - F(N,K,J,MB) & + F(N,K+1,J,MA) - F(N,K+1,J,MB) ) & -FAC2*( G(N,K,J+1,MX) - G(N,K,J-1,MX) & + G(N,K+1,J+1,MX) - G(N,K+1,J-1,MX) ) & -FAC3*( H(N,K+1,J,MX) - H(N,K,J,MX) ) & +DT/2.*( S(N,K,J,MX) + S(N,K+1,J,MX) ) 2 CONTINUE ROP = UT(1)/R2 VRP = UT(2)/UT(1) VTP = UT(3)/UT(1) VPP = UT(4)/UT(1) BRP = UT(5)/R1 BTP = UT(6)/R1 BPP = UT(7)/R1 TEMP1 = VRP**2+VTP**2+VPP**2 TEMP2 = (BRP**2+BTP**2+BPP**2)*DMU1 TPP = GAM1/(2.*ROP*R)*(UT(8)/R2-0.5*ROP*TEMP1-TEMP2) PRP = 2.*ROP*R*TPP HK(1,K,J) = R2*ROP*VPP HK(2,K,J) = R2*(ROP*VRP*VPP-BRP*BPP*DMU0) HK(3,K,J) = R2*(ROP*VTP*VPP-BTP*BPP*DMU0) TEMP2 = ROP*VPP**2 TEMP3 = (BRP**2+BTP**2-BPP**2)*DMU1 HK(4,K,J) = R2*(PRP + TEMP2 + TEMP3) HK(5,K,J) = R1*(VPP*BRP-VRP*BPP) HK(6,K,J) = R1*(VPP*BTP-VTP*BPP) HK(7,K,J) = 0. VKE = 0.5*ROP*( VRP**2 + VTP**2 + VPP**2 ) GPK = GMM1*PRP + VKE TEMP2 = BRP*(VPP*BRP-VRP*BPP)*DMU0 TEMP3 = -BTP*(VTP*BPP-VPP*BTP)*DMU0 HK(8,K,J) = R2*(VPP*GPK+TEMP2+TEMP3) SK(1,K,J) = 0. TEMP1 = 2.*R1*PRP TEMP2 = -ROP*GMS + R1*BRP**2*DMU0 TEMP3 = ROP*R1*(VTP**2+VPP**2) SK(2,K,J) = TEMP1+TEMP2+TEMP3 TEMP1 = R1*COTHM*PRP TEMP2 = (BRP**2+BTP**2-BPP**2)*R1*COTHM*DMU1 TEMP3 = R1*BRP*BTP*DMU0-ROP*R1*(VRP*VTP-VPP**2*COTHM) SK(3,K,J) = TEMP1+TEMP2+TEMP3 TEMP1 = BTP*BPP*R1*COTHM*DMU0 TEMP2 = R1*BRP*BPP*DMU0 TEMP3 = -ROP*R1*(VRP*VPP+VTP*VPP*COTHM) SK(4,K,J) = TEMP1+TEMP2+TEMP3 SK(5,K,J) = 0. SK(6,K,J) = 0. SK(7,K,J) = (VTP*BPP-VPP*BTP)*COTHM SK(8,K,J) = -ROP*VRP*GMS 10 CONTINUE 20 CONTINUE RETURN END C --- ------------------------------------------------------------------ SUBROUTINE STEPUP(I,LB,LA,MB,MX,MA ) C --- do leap frog step of U(I) into UP(LB), using c-buffs C --- UP is U at next time step INCLUDE 'igm_grid.ins' COMMON /BUFFS/ & U(8,KX,JX,3), & F(8,KX,JX,3), & G(8,KX,JX,3), & H(8,KX,JX,3), & S(8,KX,JX,3), & FI(8,KX,JX,2), & SI(8,KX,JX,2), & GJ(8,KX,JX), & SJ(8,KX,JX), & HK(8,KX,JX), & SK(8,KX,JX), & UP(8,KX,JX,2) R1 = RAD(I) FAC1 = DT/(2.*DR) DO 20 J=JC1,JC2 FAC2 = DT/(2.*R1*SIN(TH(J))*DTH) FAC3 = DT/(2.*R1*SIN(TH(J))*DPH) DO 10 K=KC1,KC2 DO 1 N=1,8 UP(N,K,J,LB) = U(N,K,J,MX) & -FAC1*( .5*( F(N,K,J,MA) - F(N,K,J,MB) ) & + FI(N,K,J,LA) - FI(N,K,J,LB) ) & -FAC2*( .5*( G(N,K,J+1,MX) - G(N,K,J-1,MX) ) & + GJ(N,K,J) - GJ(N,K,J-1) ) & -FAC3*( .5*( H(N,K+1,J,MX) - H(N,K-1,J,MX) ) & + HK(N,K,J) - HK(N,K-1,J) ) & +DT/2.*( S(N,K,J,MX) & +( SI(N,K,J,LA) + SI(N,K,J,LB) & + SJ(N,K,J-1) + SJ(N,K,J) & + SK(N,K-1,J) + SK(N,K,J) )/6. ) 1 CONTINUE 10 CONTINUE 20 CONTINUE RETURN END C --- ------------------------------------------------------------------ SUBROUTINE STOREP(I,L) C --- store UP(L) back into grid at shell I INCLUDE 'igm_grid.ins' PARAMETER ( PI1 = 3.141593, PI4 = 4.*PI1, PI8 = 8.*PI1 ) COMMON /BLOC7/ GMS,RS,GAM,GAM1,GMM1,R,PI,XMU0,XMU1,DMU0,DMU1 COMMON /BUFFS/ & U(8,KX,JX,3), & F(8,KX,JX,3), & G(8,KX,JX,3), & H(8,KX,JX,3), & S(8,KX,JX,3), & FI(8,KX,JX,2), & SI(8,KX,JX,2), & GJ(8,KX,JX), & SJ(8,KX,JX), & HK(8,KX,JX), & SK(8,KX,JX), & UP(8,KX,JX,2) R1 = RAD(I) R2 = R1**2 DO 20 J=JC1,JC2 DO 10 K=KC1,KC2 ROP = UP(1,K,J,L)/R2 VRP = UP(2,K,J,L)/UP(1,K,J,L) VTP = UP(3,K,J,L)/UP(1,K,J,L) VPP = UP(4,K,J,L)/UP(1,K,J,L) BRP = UP(5,K,J,L)/R1 BTP = UP(6,K,J,L)/R1 BPP = UP(7,K,J,L)/R1 TEMP1 = VRP**2+VTP**2+VPP**2 TEMP2 = (BRP**2+BTP**2+BPP**2)*DMU1 TPP = GAM1/(2.*ROP*R)*(UP(8,K,J,L)/R2-0.5*ROP*TEMP1-TEMP2) RO(K,J,I) = ROP VR(K,J,I) = VRP VT(K,J,I) = VTP VP(K,J,I) = VPP BR(K,J,I) = BRP BT(K,J,I) = BTP BP(K,J,I) = BPP TP(K,J,I) = TPP 10 CONTINUE 20 CONTINUE RETURN END