c c SUBROUTINE TPCORE2D c c c ******************************************************************************* c Routine to do advection on 2D model grid, EVERYTHING IN DOUBLE PRECISION c -- E. Fleming (Aug. 96) cc ****************************************************************************** C C Purpose: advects constituents from given input v* and w* on a C-grid. c The current code assumes that the individual Courant numbers c are never greater than 1 which is well within the limits of the c 2D model circulation for a 12 hour time step. Inclusion of c the semi-lagrangian code from TPCORE.F is needed if Courant c numbers are greater than 1. This is relatively easy to do. c C Schemes: Piecewise parabolic method (PPM), with cross-derivative c terms to account for multidirectional advection c without the operator splitting that the Prather scheme imposes C (see Lin and Rood 1994, 1996; also see Colella and Woodward, 1984, c and Carpenter et al., 1990). C C This is adopted from the TPCORE.F routine (version 3.5) for 3D transport c written by Shian-Jiann Lin and Richard B. Rood for the GSFC 3D CTM, c the GEOS-GCM and the GSFC Data Assimilation System (GEOS-DAS). 1D and c 2D (x-y) advection codes utilizing the PPM written by S.J. Lin were also used. c c This code can accomodate the following options on the PPM: C C IORD= 3: monotonic PPM* (slightly improved PPM of Collela & Woodward 1984) - C uses 4th order reference slopes C 4: semi-monotonic PPM (same as 3, but overshoots are allowed) C 5: positive-definite PPM (constraint on the subgrid distribution is C only strong enough to prevent generation of negative values; C both overshoots & undershootes are possible). C 6: un-constrained PPM (nearly diffusion free; slightly faster but C positivity not quaranteed. Use this option only when the fields C and winds are very smooth). C C Very little difference in the 2D model simulations of CH4 was found between c these 4 options, since standard 2D model fields are rather smooth. Note, c however, that the un-constrained option did show a small area of negative c mixing ratios. Because of this, and since the numerical diffusion seems to be c quite small (actually smaller than Prather, probably because of the operator c splitting), Lin recommends just using the fully montonic c constraint (IORD=3), since this guarantees exact tracer correlations. C C Note that the implicit numerical diffusion decreases as IORD increases. C The last two options (IORD=5, 6) should only be used when there is C significant explicit diffusion (such as a turbulence parameterization). You C might get dispersive results otherwise. C C PPM is 4th order accurate when grid spacing is uniform (y & z) as C with the 2D model; 3rd order accurate for non-uniform grid. C Again, no space or time splitting is performed. C ******************************************************************** C c c XMR(L$,Z$58,T$) mixing ratio, defined in COMMON, THINGS ACT ON MIXING RATIO #include "com2d.h" c these are defined in VWMASS/STREAMF REAL*8 w1(18,59), WBAR(2,59), W(18,59), V(19,58) REAL*8 xpi, qp1, zk, dy0, dz0, cosc(18), cose(19), XMR(L$,Z$58,T$) real*8 dtdy1, dtdz1, p2de(59), ms(59), XMR0(L$,Z$58) real*8 dqv(18,58), dqw(18,58), dqt(18,58), ady(18,58), adz(18,58) real*8 px(L$), pcr(L$), ptt(18), rx(58), rcr(58), rtt(58) real*8 px1(0:L$+1), rx1(0:Z$58+1) real*8 vc(19), ymass(19), vav(L$), wc(59), zmass(59), wav(Z$58) real*8 fymr(L$+1), ffy(L$+1), fzmr(Z$58+1), ffz(Z$58+1) real*8 sum1, sum2, sumd, zzzz LOGICAL LMASS data IORD/3/ c define pressure and ms array at 1/2 grid points, number density in 1/cm3 at constant T = 239.27K do 103 ik=1,59 p2de(ik) = 1013.*DEXP(-.2844d0*(ik-1.)) 103 ms(ik) = p2de(ik)*1.d3/(28.97*2.87d6*239.27*1.66d-24) DTDY1 = DT/DY0 DTDZ1 = DT/DZ0 c Check Mass conservation (total number of molecules) before (and after) advection every 60 days, CO2 here C LMASS = .TRUE. LMASS = .FALSE. IF (LMASS) THEN if (mod(iday360-15,60) .eq. 0.0) then sum1 = 0.d0 do 990 ik=1,58 do 990 ij=1,18 990 sum1 = sum1 + xmr(ij,ik,31)*ms(ik)*DEXP(-zk)*cosc(ij) endif ENDIF c **************** loop through all transported species ********* DO 1000 INS=1,ITRANS c first load mixing ratios into temporary array do 150 ik=1,Z$58 do 150 ij=1,L$ 150 XMR0(ij,ik) = XMR(ij,ik,INS) C ***** Y SETUP -- COMPUTE X-TERMS **************** DO 160 ik=1,58 do 162 ij=1,18 162 px1(ij) = xmr0(ij,ik) c extend array one grid point each side, Lin says just use constant value px1(0) = px1(1) px1(L$+1) = px1(L$) do 164 ij=1,19 164 vc(ij) = v(ij,ik)*dtdy1 ! vc(19), Courant number do 166 ij=1,L$ ! vav(18), box-averged Courant number 166 vav(ij) = (vc(ij) + vc(ij+1))*.5 CALL CROSST(L$, VAV, PX1, PCR, PTT, IORD) !px1(0:L$+1), PCR, PTT(18)=total (PX1+PCR) x-terms do 168 ij=1,18 168 ady(ij,ik) = ptt(ij) 160 CONTINUE C ***** Z SETUP -- COMPUTE X-TERMS **************** DO 170 ij=1,18 do 172 ik=1,58 172 rx1(ik) = xmr0(ij,ik) c extend array one grid point each side, Lin says just use constant value rx1(0) = rx1(1) rx1(Z$58+1) = rx1(Z$58) do 174 ik=1,59 174 wc(ik) = w(ij,ik)*DTDZ1 ! wc(59), Courant number do 175 ik=1,Z$58 ! wav(58), box-averged Courant number 175 wav(ik) = (wc(ik) + wc(ik+1))*.5 CALL CROSST(Z$58, WAV, RX1, RCR, RTT, IORD) !rx1(0:59), RCR,RTT(58) = total(RX1+RCR) x-terms do 177 ik=1,58 177 adz(ij,ik) = rtt(ik) 170 CONTINUE C ******** Do Y ADVECTION [conservative (liberal) transport] ****** DO 210 ik=1,58 !! v(19,58), adz,dqv(18,58) do 212 ij=1,19 vc(ij) = v(ij,ik)*DTDY1 212 ymass(ij) = vc(ij)*ms(ik)*DEXP(-zk)*cose(ij) !! horizontal mass flux; ymass, vc, cose(19) do 215 ij=1,18 215 px(ij) = adz(ij,ik) ! px(18), use adz array here because of cross terms CALL FYPPM(VC, PX, fymr, IORD) do 600 ij=1,L$+1 ! constit. mass flux across box edge 600 ffy(ij) = fymr(ij)*ymass(ij) ! ymass=0 and so fx=0 at poles (90s, 90N) do 610 ij=1,L$ ! constit mass flux convergence at box center 610 dqv(ij,ik) = -(ffy(ij+1) - ffy(ij))/cosc(ij) 210 CONTINUE C ******** Z (vertical) ADVECTION **************** DO 220 ij=1,18 ! w(18,59), ady,dqw(18,58) do 223 ik=1,59 wc(ik) = w(ij,ik)*DTDZ1 !! ms, WC, zmass(59) = vertical mass flux 223 zmass(ik) = wc(ik)*ms(ik) do 225 ik=1,58 225 rx(ik) = ady(ij,ik) ! rx(58), use ady array here because of cross terms CALL FZPPM(WC, RX, fzmr, IORD) do 700 ik=1,Z$58+1 ! constit. mass flux across box edge - top/bottom 700 ffz(ik) = fzmr(ik)*zmass(ik) ! zmass=0 and so fx=0 at top and bottom boundary do 710 ik=1,Z$58 ! constit mass flux convergence at box center 710 dqw(ij,ik) = -(ffz(ik+1) - ffz(ik)) 220 CONTINUE c update new mixing ratio at box center, take old mixing ratio and add on change c in number density in both horiz and vertical, converted back to mixing ratio c (don't divide by cosine here, it's already been done) -- dqv, dqw, dqt (18,58) do 330 ik=1,58 do 330 ij=1,18 330 dqt(ij,ik) = (dqv(ij,ik) + dqw(ij,ik))/(ms(ik)*DEXP(-zk)) do 440 ik=1,Z$58 do 440 ij=1,L$ XMR(ij,ik,INS) = XMR0(ij,ik) + dqt(ij,ik) ! UPDATE w/ mixing ratio change c sometimes get negatives in solid HNO3, H2O c because of sharp density changes, readjust here zzzz = XMR(ij,ik,INS)*m(ij,ik) if (zzzz .le. 1.D-12) XMR(ij,ik,INS) = 1.D-12/m(ij,ik) 440 continue 1000 CONTINUE ! end species loop c Check Mass conservation (total number of molecules) (before) and after advection every 60 days, SPECIES 2 IF (LMASS) THEN if (mod(iday360-15,60) .eq. 0.0) then sum2 = 0.d0 do 992 ik=1,58 do 992 ij=1,18 992 sum2 = sum2 + xmr(ij,ik,31)*ms(ik)*DEXP(-zk)*cosc(ij) sumd = (sum2-sum1)/sum1 write(6,995) iday360, iyr, sum1, sum2, sumd 995 format(1x, 2I3, 3x, 1P3D24.15) endif ENDIF SAVE RETURN END CC ************************* SUBROUTINES ************************************************************** subroutine FYPPM(VC, PX, fluxt, IORD) integer L$ real*8 R3, R23, R24 parameter (R3 = 1.d0/3., R23 = 2.d0/3., R24 = 1./24.d0, L$=18) real*8 VC(L$+1), px(L$), DC(L$) real*8 fluxt(L$+1), AL(L$), AR(L$), A6(L$), tmp, pmax, pmin ILMT = IORD - 3 do 100 ij=2,L$-1 if (ij .ge. 3 .and. ij .le. L$-2) then tmp = (8.d0*(px(ij+1) - px(ij-1)) + px(ij-2) - px(ij+2))*R24 !4th order difference for 65S-65N else tmp = .25*(px(ij+1) - px(ij-1)) !2nd order difference for 75S,75N endif Pmax = DMAX1(px(ij-1), px(ij), px(ij+1)) - px(ij) Pmin = px(ij) - DMIN1(px(ij-1), px(ij), px(ij+1)) 100 DC(ij) = DSIGN(DMIN1(DABS(tmp), Pmax, Pmin), tmp) ! sign func: result = abs(1st)*sgn(2nd) DC(1) = 0.d0 DC(L$) = 0.d0 cc Compute values at box edges, for poles, extrapolate 1/2 grid point cc to get first guess for AL(1) and AR(L$) - these will be modified in LMTPPM cc just use constant gradient (ratio) for extrapolation, need sqrt since its 1/2 grid point c do 400 ij=2,L$ 400 AL(ij) = 0.5d0*(px(ij-1)+px(ij)) + (DC(ij-1) - DC(ij))*R3 AL(1) = DSQRT(px(1)/px(2))*px(1) do 402 ij=1,L$-1 402 AR(ij) = AL(ij+1) AR(L$) = DSQRT(px(L$)/px(L$-1))*px(L$) do 404 ij=1,L$ 404 A6(ij) = 3.d0*(px(ij)+px(ij) - (AL(ij)+AR(ij))) if (ILMT .le. 2) call LMTPPM(L$, PX, DC, A6, AR, AL, ILMT) ! adjusts AL, AR, A6 arrays do 450 ij = 2,L$ IF (VC(ij) .GT. 0.) then fluxt(ij) = AR(ij-1) + 0.5d0*VC(ij)*(AL(ij-1) - AR(ij-1) + c A6(ij-1)*(1.-R23*VC(ij)) ) else fluxt(ij) = AL(ij) - 0.5d0*VC(ij)*(AR(ij) - AL(ij) + c A6(ij)*(1.+R23*VC(ij))) ENDIF 450 continue c at poles, VC = 0 and YMASS = 0, so it doesn't matter c what fluxt(poles) is, just set to zero to define fluxt(1) = 0.d0 fluxt(L$+1) = 0.d0 RETURN END subroutine FZPPM(WC, PX, fluxt, IORD) integer Z$58 real*8 R3, R23, R24 parameter(R3= 1.d0/3., R23 = 2.d0/3., R24 = 1./24.d0, Z$58=58) real*8 WC(Z$58+1), px(Z$58), DC(Z$58), tmp, pmax, pmin real*8 fluxt(Z$58+1), AL(Z$58), AR(Z$58), A6(Z$58) ILMT = IORD - 3 do 100 ij=2,Z$58-1 if (ij .ge. 3 .and. ij .le. Z$58-2) then tmp = (8.d0*(px(ij+1) - px(ij-1)) + px(ij-2) - px(ij+2))*R24 !4th order difference for levs 3-56 else tmp = .25*(px(ij+1) - px(ij-1)) !2nd order difference for levs 2,57 endif Pmax = DMAX1(px(ij-1), px(ij), px(ij+1)) - px(ij) Pmin = px(ij) - DMIN1(px(ij-1), px(ij), px(ij+1)) 100 DC(ij) = DSIGN(DMIN1(DABS(tmp), Pmax, Pmin), tmp) ! sign func: result = abs(1st)*sgn(2nd) DC(1) = 0.d0 DC(Z$58) = 0.d0 cc Compute values at box edges, for top/bottom boundary, extrapolate 1/2 grid point cc to get first guess for AL(1) and AR(58) - these will be modified in LMTPPM cc just use constant gradient (ratio) for extrapolation, need sqrt since its 1/2 grid point c do 400 ij=2,Z$58 400 AL(ij) = 0.5d0*(px(ij-1)+px(ij)) + (DC(ij-1) - DC(ij))*R3 AL(1) = DSQRT(px(1)/px(2))*px(1) do 402 ij=1,Z$58-1 402 AR(ij) = AL(ij+1) AR(Z$58) = DSQRT(px(Z$58)/px(Z$58-1))*px(Z$58) do 404 ij=1,Z$58 404 A6(ij) = 3.d0*(px(ij)+px(ij) - (AL(ij)+AR(ij))) if (ILMT .le. 2) call LMTPPM(Z$58, PX, DC, A6, AR, AL, ILMT) ! adjusts AL, AR, A6 arrays do 450 ij = 2,Z$58 IF (WC(ij) .GT. 0.) then fluxt(ij) = AR(ij-1) + 0.5d0*WC(ij)*(AL(ij-1) - AR(ij-1) + c A6(ij-1)*(1.-R23*WC(ij)) ) else fluxt(ij) = AL(ij) - 0.5d0*WC(ij)*(AR(ij) - AL(ij) + c A6(ij)*(1.+R23*WC(ij))) ENDIF 450 continue c at top/bottom, WC = 0 and ZMASS = 0, so it doesn't matter c what fluxt there is, just set to zero to define fluxt(1) = 0.d0 fluxt(Z$58+1) = 0.d0 RETURN END subroutine LMTPPM(IM, PX, DCP, A6, AR, AL, ILMT) C C A6 = CURVATURE OF THE TEST PARABOLA C AR = RIGHT EDGE VALUE OF THE TEST PARABOLA C AL = LEFT EDGE VALUE OF THE TEST PARABOLA C DCP = 0.5 * MISMATCH C P = CELL-AVERAGED VALUE C IM = VECTOR LENGTH C C OPTIONS: C C ILMT = 0: FULL MONOTONICITY C ILMT = 1: SEMI-MONOTONIC CONSTRAINT (NO UNDERSHOOTS) C ILMT = 2: POSITIVE-DEFINITE CONSTRAINT C real*8 R12 parameter ( R12 = 1./12. ) real*8 A6(IM), AR(IM), AL(IM), px(IM), DCP(IM) real*8 da1, da2, A6DA, fmin C if (ILMT .eq. 0) then ! Full constraint do 100 i=1,IM if (DCP(i) .eq. 0.) then AR(i) = px(i) AL(i) = px(i) A6(i) = 0.d0 else da1 = AR(i) - AL(i) da2 = da1**2 A6DA = A6(i)*da1 if (A6DA .lt. -da2) then A6(i) = 3.*(AL(i)-px(i)) AR(i) = AL(i) - A6(i) elseif (A6DA .gt. da2) then A6(i) = 3.*(AR(i)-px(i)) AL(i) = AR(i) - A6(i) endif endif 100 continue elseif (ILMT .eq. 1) then ! Semi-monotonic constraint do 150 i=1,IM if (DABS(AR(i)-AL(i)) .GE. -A6(i)) go to 150 if (px(i) .lt. AR(i) .and. px(i) .lt. AL(i)) then AR(i) = px(i) AL(i) = px(i) A6(i) = 0.d0 elseif (AR(i) .gt. AL(i)) then A6(i) = 3.*(AL(i)-px(i)) AR(i) = AL(i) - A6(i) else A6(i) = 3.*(AR(i)-px(i)) AL(i) = AR(i) - A6(i) endif 150 continue elseif (ILMT .eq. 2) then ! Positive Definite constraint do 250 i=1,IM if (DABS(AR(i)-AL(i)) .GE. -A6(i)) go to 250 fmin = px(i) + 0.25*(AR(i)-AL(i))**2/A6(i) + A6(i)*R12 if (fmin .ge. 0.) go to 250 if (px(i) .lt. AR(i) .and. px(i) .lt. AL(i)) then AR(i) = px(i) AL(i) = px(i) A6(i) = 0. elseif (AR(i) .gt. AL(i)) then A6(i) = 3.*(AL(i)-px(i)) AR(i) = AL(i) - A6(i) else A6(i) = 3.*(AR(i)-px(i)) AL(i) = AR(i) - A6(i) endif 250 continue ENDIF RETURN END subroutine CROSST(IN, UA0, PX, PCR, PTT, IORD) ! px(0:IN+1), UA0(IN), pcr,ptt(in) c this computes the cross-derivative (SPLITTING) term which arise from the sequential c splitting process, and is critical to the stability of the scheme. It can be interpreted c as the convergence of a secondary oblique flux. It arises since the physical transport c can neither be perfectly modeled by the indepedent applications of the orthogonal w* c and v* fields, nor can it be perfectly modeled by a sequential process of w* advection c and then v* advection as was done by the Prather scheme. This term vanishes for a c uniform mixing ratio field. See Lin and Rood (1996) for a complete discussion of this. c Also, inclusion of the splitting terms are necessary for stability, ie, MAX(Cy, Cz) <1, c instead of Cy+Cz <1 (Cy, Cz are the Courant nums in y,z), the larger the Cy and Cz, the c larger the cross terms, UA0 = box averaged Courant number c real*8 ua0(in), px(0:IN+1), ptt(in), pcr(in), rut do 120 ij=1,IN iu = JIDINT(ua0(ij)) rut = ua0(ij) - iu iiu = ij - iu if (ua0(ij) .ge. 0.0) then pcr(ij) = px(iiu) + rut*(px(iiu-1) - px(iiu)) else pcr(ij) = px(iiu) + rut*(px(iiu) - px(iiu+1)) endif 120 continue do 300 ij=1,IN 300 ptt(ij) = 0.5*(pcr(ij) + px(ij)) RETURN END