--- trunk/Sources/dyn3d/ppm3d.f 2015/07/29 09:52:33 165 +++ trunk/Sources/dyn3d/PPM3d/ppm3d.f 2015/07/29 14:32:55 166 @@ -1,68 +1,3 @@ - -! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/ppm3d.F,v 1.1.1.1 2004/05/19 -! 12:53:07 lmdzadmin Exp $ - - -! From lin@explorer.gsfc.nasa.gov Wed Apr 15 17:44:44 1998 -! Date: Wed, 15 Apr 1998 11:37:03 -0400 -! From: lin@explorer.gsfc.nasa.gov -! To: Frederic.Hourdin@lmd.jussieu.fr -! Subject: 3D transport module of the GSFC CTM and GEOS GCM - - -! This code is sent to you by S-J Lin, DAO, NASA-GSFC - -! Note: this version is intended for machines like CRAY -! -90. No multitasking directives implemented. - - -! ******************************************************************** - -! TransPort Core for Goddard Chemistry Transport Model (G-CTM), Goddard -! Earth Observing System General Circulation Model (GEOS-GCM), and Data -! Assimilation System (GEOS-DAS). - -! ******************************************************************** - -! Purpose: given horizontal winds on a hybrid sigma-p surfaces, -! one call to tpcore updates the 3-D mixing ratio -! fields one time step (NDT). [vertical mass flux is computed -! internally consistent with the discretized hydrostatic mass -! continuity equation of the C-Grid GEOS-GCM (for IGD=1)]. - -! Schemes: Multi-dimensional Flux Form Semi-Lagrangian (FFSL) scheme based -! on the van Leer or PPM. -! (see Lin and Rood 1996). -! Version 4.5 -! Last modified: Dec. 5, 1996 -! Major changes from version 4.0: a more general vertical hybrid sigma- -! pressure coordinate. -! Subroutines modified: xtp, ytp, fzppm, qckxyz -! Subroutines deleted: vanz - -! Author: Shian-Jiann Lin -! mail address: -! Shian-Jiann Lin* -! Code 910.3, NASA/GSFC, Greenbelt, MD 20771 -! Phone: 301-286-9540 -! E-mail: lin@dao.gsfc.nasa.gov - -! *affiliation: -! Joint Center for Earth Systems Technology -! The University of Maryland Baltimore County -! NASA - Goddard Space Flight Center -! References: - -! 1. Lin, S.-J., and R. B. Rood, 1996: Multidimensional flux form semi- -! Lagrangian transport schemes. Mon. Wea. Rev., 124, 2046-2070. - -! 2. Lin, S.-J., W. C. Chao, Y. C. Sud, and G. K. Walker, 1994: A class of -! the van Leer-type transport schemes and its applications to the moist- -! ure transport in a General Circulation Model. Mon. Wea. Rev., 122, -! 1575-1593. - -! ****6***0*********0*********0*********0*********0*********0**********72 - SUBROUTINE ppm3d(igd, q, ps1, ps2, u, v, w, ndt, iord, jord, kord, nc, imr, & jnp, j1, nlay, ap, bp, pt, ae, fill, dum, umax) @@ -786,1243 +721,3 @@ END SUBROUTINE ppm3d ! ****6***0*********0*********0*********0*********0*********0**********72 -SUBROUTINE fzppm(imr, jnp, nlay, j1, dq, wz, p, dc, dqdt, ar, al, a6, flux, & - wk1, wk2, wz2, delp, kord) - PARAMETER (kmax=150) - PARAMETER (r23=2./3., r3=1./3.) - REAL wz(imr, jnp, nlay), p(imr, jnp, nlay), dc(imr, jnp, nlay), & - wk1(imr, *), delp(imr, jnp, nlay), dq(imr, jnp, nlay), & - dqdt(imr, jnp, nlay) - ! Assuming JNP >= NLAY - REAL ar(imr, *), al(imr, *), a6(imr, *), flux(imr, *), wk2(imr, *), & - wz2(imr, *) - - jmr = jnp - 1 - imjm = imr*jnp - nlaym1 = nlay - 1 - - lmt = kord - 3 - - ! ****6***0*********0*********0*********0*********0*********0**********72 - ! Compute DC for PPM - ! ****6***0*********0*********0*********0*********0*********0**********72 - - DO k = 1, nlaym1 - DO i = 1, imjm - dqdt(i, 1, k) = p(i, 1, k+1) - p(i, 1, k) - END DO - END DO - - DO k = 2, nlaym1 - DO i = 1, imjm - c0 = delp(i, 1, k)/(delp(i,1,k-1)+delp(i,1,k)+delp(i,1,k+1)) - c1 = (delp(i,1,k-1)+0.5*delp(i,1,k))/(delp(i,1,k+1)+delp(i,1,k)) - c2 = (delp(i,1,k+1)+0.5*delp(i,1,k))/(delp(i,1,k-1)+delp(i,1,k)) - tmp = c0*(c1*dqdt(i,1,k)+c2*dqdt(i,1,k-1)) - qmax = max(p(i,1,k-1), p(i,1,k), p(i,1,k+1)) - p(i, 1, k) - qmin = p(i, 1, k) - min(p(i,1,k-1), p(i,1,k), p(i,1,k+1)) - dc(i, 1, k) = sign(min(abs(tmp),qmax,qmin), tmp) - END DO - END DO - - - ! ****6***0*********0*********0*********0*********0*********0**********72 - ! Loop over latitudes (to save memory) - ! ****6***0*********0*********0*********0*********0*********0**********72 - - DO j = 1, jnp - IF ((j==2 .OR. j==jmr) .AND. j1/=2) GO TO 2000 - - DO k = 1, nlay - DO i = 1, imr - wz2(i, k) = wz(i, j, k) - wk1(i, k) = p(i, j, k) - wk2(i, k) = delp(i, j, k) - flux(i, k) = dc(i, j, k) !this flux is actually the monotone slope - END DO - END DO - - ! ****6***0*********0*********0*********0*********0*********0**********72 - ! Compute first guesses at cell interfaces - ! First guesses are required to be continuous. - ! ****6***0*********0*********0*********0*********0*********0**********72 - - ! three-cell parabolic subgrid distribution at model top - ! two-cell parabolic with zero gradient subgrid distribution - ! at the surface. - - ! First guess top edge value - DO i = 1, imr - ! three-cell PPM - ! Compute a,b, and c of q = aP**2 + bP + c using cell averages and delp - a = 3.*(dqdt(i,j,2)-dqdt(i,j,1)*(wk2(i,2)+wk2(i,3))/(wk2(i,1)+wk2(i, & - 2)))/((wk2(i,2)+wk2(i,3))*(wk2(i,1)+wk2(i,2)+wk2(i,3))) - b = 2.*dqdt(i, j, 1)/(wk2(i,1)+wk2(i,2)) - r23*a*(2.*wk2(i,1)+wk2(i,2)) - al(i, 1) = wk1(i, 1) - wk2(i, 1)*(r3*a*wk2(i,1)+0.5*b) - al(i, 2) = wk2(i, 1)*(a*wk2(i,1)+b) + al(i, 1) - - ! Check if change sign - IF (wk1(i,1)*al(i,1)<=0.) THEN - al(i, 1) = 0. - flux(i, 1) = 0. - ELSE - flux(i, 1) = wk1(i, 1) - al(i, 1) - END IF - END DO - - ! Bottom - DO i = 1, imr - ! 2-cell PPM with zero gradient right at the surface - - fct = dqdt(i, j, nlaym1)*wk2(i, nlay)**2/((wk2(i,nlay)+wk2(i, & - nlaym1))*(2.*wk2(i,nlay)+wk2(i,nlaym1))) - ar(i, nlay) = wk1(i, nlay) + fct - al(i, nlay) = wk1(i, nlay) - (fct+fct) - IF (wk1(i,nlay)*ar(i,nlay)<=0.) ar(i, nlay) = 0. - flux(i, nlay) = ar(i, nlay) - wk1(i, nlay) - END DO - - - ! ****6***0*********0*********0*********0*********0*********0**********72 - ! 4th order interpolation in the interior. - ! ****6***0*********0*********0*********0*********0*********0**********72 - - DO k = 3, nlaym1 - DO i = 1, imr - c1 = dqdt(i, j, k-1)*wk2(i, k-1)/(wk2(i,k-1)+wk2(i,k)) - c2 = 2./(wk2(i,k-2)+wk2(i,k-1)+wk2(i,k)+wk2(i,k+1)) - a1 = (wk2(i,k-2)+wk2(i,k-1))/(2.*wk2(i,k-1)+wk2(i,k)) - a2 = (wk2(i,k)+wk2(i,k+1))/(2.*wk2(i,k)+wk2(i,k-1)) - al(i, k) = wk1(i, k-1) + c1 + c2*(wk2(i,k)*(c1*(a1-a2)+a2*flux(i, & - k-1))-wk2(i,k-1)*a1*flux(i,k)) - END DO - END DO - - DO i = 1, imr*nlaym1 - ar(i, 1) = al(i, 2) - END DO - - DO i = 1, imr*nlay - a6(i, 1) = 3.*(wk1(i,1)+wk1(i,1)-(al(i,1)+ar(i,1))) - END DO - - ! ****6***0*********0*********0*********0*********0*********0**********72 - ! Top & Bot always monotonic - CALL lmtppm(flux(1,1), a6(1,1), ar(1,1), al(1,1), wk1(1,1), imr, 0) - CALL lmtppm(flux(1,nlay), a6(1,nlay), ar(1,nlay), al(1,nlay), & - wk1(1,nlay), imr, 0) - - ! Interior depending on KORD - IF (lmt<=2) CALL lmtppm(flux(1,2), a6(1,2), ar(1,2), al(1,2), wk1(1,2), & - imr*(nlay-2), lmt) - - ! ****6***0*********0*********0*********0*********0*********0**********72 - - DO i = 1, imr*nlaym1 - IF (wz2(i,1)>0.) THEN - cm = wz2(i, 1)/wk2(i, 1) - flux(i, 2) = ar(i, 1) + 0.5*cm*(al(i,1)-ar(i,1)+a6(i,1)*(1.-r23*cm)) - ELSE - cp = wz2(i, 1)/wk2(i, 2) - flux(i, 2) = al(i, 2) + 0.5*cp*(al(i,2)-ar(i,2)-a6(i,2)*(1.+r23*cp)) - END IF - END DO - - DO i = 1, imr*nlaym1 - flux(i, 2) = wz2(i, 1)*flux(i, 2) - END DO - - DO i = 1, imr - dq(i, j, 1) = dq(i, j, 1) - flux(i, 2) - dq(i, j, nlay) = dq(i, j, nlay) + flux(i, nlay) - END DO - - DO k = 2, nlaym1 - DO i = 1, imr - dq(i, j, k) = dq(i, j, k) + flux(i, k) - flux(i, k+1) - END DO - END DO -2000 END DO - RETURN -END SUBROUTINE fzppm - -SUBROUTINE xtp(imr, jnp, iml, j1, j2, jn, js, pu, dq, q, uc, fx1, xmass, & - iord) - DIMENSION uc(imr, *), dc(-iml:imr+iml+1), xmass(imr, jnp), fx1(imr+1), & - dq(imr, jnp), qtmp(-iml:imr+1+iml) - DIMENSION pu(imr, jnp), q(imr, jnp), isave(imr) - - imp = imr + 1 - - ! van Leer at high latitudes - jvan = max(1, jnp/18) - j1vl = j1 + jvan - j2vl = j2 - jvan - - DO j = j1, j2 - - DO i = 1, imr - qtmp(i) = q(i, j) - END DO - - IF (j>=jn .OR. j<=js) GO TO 2222 - ! ************* Eulerian ********** - - qtmp(0) = q(imr, j) - qtmp(-1) = q(imr-1, j) - qtmp(imp) = q(1, j) - qtmp(imp+1) = q(2, j) - - IF (iord==1 .OR. j==j1 .OR. j==j2) THEN - DO i = 1, imr - iu = float(i) - uc(i, j) - fx1(i) = qtmp(iu) - END DO - ELSE - CALL xmist(imr, iml, qtmp, dc) - dc(0) = dc(imr) - - IF (iord==2 .OR. j<=j1vl .OR. j>=j2vl) THEN - DO i = 1, imr - iu = float(i) - uc(i, j) - fx1(i) = qtmp(iu) + dc(iu)*(sign(1.,uc(i,j))-uc(i,j)) - END DO - ELSE - CALL fxppm(imr, iml, uc(1,j), qtmp, dc, fx1, iord) - END IF - - END IF - - DO i = 1, imr - fx1(i) = fx1(i)*xmass(i, j) - END DO - - GO TO 1309 - - ! ***** Conservative (flux-form) Semi-Lagrangian transport ***** - -2222 CONTINUE - - DO i = -iml, 0 - qtmp(i) = q(imr+i, j) - qtmp(imp-i) = q(1-i, j) - END DO - - IF (iord==1 .OR. j==j1 .OR. j==j2) THEN - DO i = 1, imr - itmp = int(uc(i,j)) - isave(i) = i - itmp - iu = i - uc(i, j) - fx1(i) = (uc(i,j)-itmp)*qtmp(iu) - END DO - ELSE - CALL xmist(imr, iml, qtmp, dc) - - DO i = -iml, 0 - dc(i) = dc(imr+i) - dc(imp-i) = dc(1-i) - END DO - - DO i = 1, imr - itmp = int(uc(i,j)) - rut = uc(i, j) - itmp - isave(i) = i - itmp - iu = i - uc(i, j) - fx1(i) = rut*(qtmp(iu)+dc(iu)*(sign(1.,rut)-rut)) - END DO - END IF - - DO i = 1, imr - IF (uc(i,j)>1.) THEN - ! DIR$ NOVECTOR - DO ist = isave(i), i - 1 - fx1(i) = fx1(i) + qtmp(ist) - END DO - ELSE IF (uc(i,j)<-1.) THEN - DO ist = i, isave(i) - 1 - fx1(i) = fx1(i) - qtmp(ist) - END DO - ! DIR$ VECTOR - END IF - END DO - DO i = 1, imr - fx1(i) = pu(i, j)*fx1(i) - END DO - - ! *************************************** - -1309 fx1(imp) = fx1(1) - DO i = 1, imr - dq(i, j) = dq(i, j) + fx1(i) - fx1(i+1) - END DO - - ! *************************************** - - END DO - RETURN -END SUBROUTINE xtp - -SUBROUTINE fxppm(imr, iml, ut, p, dc, flux, iord) - PARAMETER (r3=1./3., r23=2./3.) - DIMENSION ut(*), flux(*), p(-iml:imr+iml+1), dc(-iml:imr+iml+1) - DIMENSION ar(0:imr), al(0:imr), a6(0:imr) - INTEGER lmt - ! logical first - ! data first /.true./ - ! SAVE LMT - ! if(first) then - - ! correction calcul de LMT a chaque passage pour pouvoir choisir - ! plusieurs schemas PPM pour differents traceurs - ! IF (IORD.LE.0) then - ! if(IMR.GE.144) then - ! LMT = 0 - ! elseif(IMR.GE.72) then - ! LMT = 1 - ! else - ! LMT = 2 - ! endif - ! else - ! LMT = IORD - 3 - ! endif - - lmt = iord - 3 - - DO i = 1, imr - al(i) = 0.5*(p(i-1)+p(i)) + (dc(i-1)-dc(i))*r3 - END DO - - DO i = 1, imr - 1 - ar(i) = al(i+1) - END DO - ar(imr) = al(1) - - DO i = 1, imr - a6(i) = 3.*(p(i)+p(i)-(al(i)+ar(i))) - END DO - - IF (lmt<=2) CALL lmtppm(dc(1), a6(1), ar(1), al(1), p(1), imr, lmt) - - al(0) = al(imr) - ar(0) = ar(imr) - a6(0) = a6(imr) - - DO i = 1, imr - IF (ut(i)>0.) THEN - flux(i) = ar(i-1) + 0.5*ut(i)*(al(i-1)-ar(i-1)+a6(i-1)*(1.-r23*ut(i))) - ELSE - flux(i) = al(i) - 0.5*ut(i)*(ar(i)-al(i)+a6(i)*(1.+r23*ut(i))) - END IF - END DO - RETURN -END SUBROUTINE fxppm - -SUBROUTINE xmist(imr, iml, p, dc) - PARAMETER (r24=1./24.) - DIMENSION p(-iml:imr+1+iml), dc(-iml:imr+1+iml) - - DO i = 1, imr - tmp = r24*(8.*(p(i+1)-p(i-1))+p(i-2)-p(i+2)) - pmax = max(p(i-1), p(i), p(i+1)) - p(i) - pmin = p(i) - min(p(i-1), p(i), p(i+1)) - dc(i) = sign(min(abs(tmp),pmax,pmin), tmp) - END DO - RETURN -END SUBROUTINE xmist - -SUBROUTINE ytp(imr, jnp, j1, j2, acosp, rcap, dq, p, vc, dc2, ymass, fx, a6, & - ar, al, jord) - DIMENSION p(imr, jnp), vc(imr, jnp), ymass(imr, jnp), dc2(imr, jnp), & - dq(imr, jnp), acosp(jnp) - ! Work array - DIMENSION fx(imr, jnp), ar(imr, jnp), al(imr, jnp), a6(imr, jnp) - - jmr = jnp - 1 - len = imr*(j2-j1+2) - - IF (jord==1) THEN - DO i = 1, len - jt = float(j1) - vc(i, j1) - fx(i, j1) = p(i, jt) - END DO - ELSE - - CALL ymist(imr, jnp, j1, p, dc2, 4) - - IF (jord<=0 .OR. jord>=3) THEN - - CALL fyppm(vc, p, dc2, fx, imr, jnp, j1, j2, a6, ar, al, jord) - - ELSE - DO i = 1, len - jt = float(j1) - vc(i, j1) - fx(i, j1) = p(i, jt) + (sign(1.,vc(i,j1))-vc(i,j1))*dc2(i, jt) - END DO - END IF - END IF - - DO i = 1, len - fx(i, j1) = fx(i, j1)*ymass(i, j1) - END DO - - DO j = j1, j2 - DO i = 1, imr - dq(i, j) = dq(i, j) + (fx(i,j)-fx(i,j+1))*acosp(j) - END DO - END DO - - ! Poles - sum1 = fx(imr, j1) - sum2 = fx(imr, j2+1) - DO i = 1, imr - 1 - sum1 = sum1 + fx(i, j1) - sum2 = sum2 + fx(i, j2+1) - END DO - - sum1 = dq(1, 1) - sum1*rcap - sum2 = dq(1, jnp) + sum2*rcap - DO i = 1, imr - dq(i, 1) = sum1 - dq(i, jnp) = sum2 - END DO - - IF (j1/=2) THEN - DO i = 1, imr - dq(i, 2) = sum1 - dq(i, jmr) = sum2 - END DO - END IF - - RETURN -END SUBROUTINE ytp - -SUBROUTINE ymist(imr, jnp, j1, p, dc, id) - PARAMETER (r24=1./24.) - DIMENSION p(imr, jnp), dc(imr, jnp) - - imh = imr/2 - jmr = jnp - 1 - ijm3 = imr*(jmr-3) - - IF (id==2) THEN - DO i = 1, imr*(jmr-1) - tmp = 0.25*(p(i,3)-p(i,1)) - pmax = max(p(i,1), p(i,2), p(i,3)) - p(i, 2) - pmin = p(i, 2) - min(p(i,1), p(i,2), p(i,3)) - dc(i, 2) = sign(min(abs(tmp),pmin,pmax), tmp) - END DO - ELSE - DO i = 1, imh - ! J=2 - tmp = (8.*(p(i,3)-p(i,1))+p(i+imh,2)-p(i,4))*r24 - pmax = max(p(i,1), p(i,2), p(i,3)) - p(i, 2) - pmin = p(i, 2) - min(p(i,1), p(i,2), p(i,3)) - dc(i, 2) = sign(min(abs(tmp),pmin,pmax), tmp) - ! J=JMR - tmp = (8.*(p(i,jnp)-p(i,jmr-1))+p(i,jmr-2)-p(i+imh,jmr))*r24 - pmax = max(p(i,jmr-1), p(i,jmr), p(i,jnp)) - p(i, jmr) - pmin = p(i, jmr) - min(p(i,jmr-1), p(i,jmr), p(i,jnp)) - dc(i, jmr) = sign(min(abs(tmp),pmin,pmax), tmp) - END DO - DO i = imh + 1, imr - ! J=2 - tmp = (8.*(p(i,3)-p(i,1))+p(i-imh,2)-p(i,4))*r24 - pmax = max(p(i,1), p(i,2), p(i,3)) - p(i, 2) - pmin = p(i, 2) - min(p(i,1), p(i,2), p(i,3)) - dc(i, 2) = sign(min(abs(tmp),pmin,pmax), tmp) - ! J=JMR - tmp = (8.*(p(i,jnp)-p(i,jmr-1))+p(i,jmr-2)-p(i-imh,jmr))*r24 - pmax = max(p(i,jmr-1), p(i,jmr), p(i,jnp)) - p(i, jmr) - pmin = p(i, jmr) - min(p(i,jmr-1), p(i,jmr), p(i,jnp)) - dc(i, jmr) = sign(min(abs(tmp),pmin,pmax), tmp) - END DO - - DO i = 1, ijm3 - tmp = (8.*(p(i,4)-p(i,2))+p(i,1)-p(i,5))*r24 - pmax = max(p(i,2), p(i,3), p(i,4)) - p(i, 3) - pmin = p(i, 3) - min(p(i,2), p(i,3), p(i,4)) - dc(i, 3) = sign(min(abs(tmp),pmin,pmax), tmp) - END DO - END IF - - IF (j1/=2) THEN - DO i = 1, imr - dc(i, 1) = 0. - dc(i, jnp) = 0. - END DO - ELSE - ! Determine slopes in polar caps for scalars! - - DO i = 1, imh - ! South - tmp = 0.25*(p(i,2)-p(i+imh,2)) - pmax = max(p(i,2), p(i,1), p(i+imh,2)) - p(i, 1) - pmin = p(i, 1) - min(p(i,2), p(i,1), p(i+imh,2)) - dc(i, 1) = sign(min(abs(tmp),pmax,pmin), tmp) - ! North. - tmp = 0.25*(p(i+imh,jmr)-p(i,jmr)) - pmax = max(p(i+imh,jmr), p(i,jnp), p(i,jmr)) - p(i, jnp) - pmin = p(i, jnp) - min(p(i+imh,jmr), p(i,jnp), p(i,jmr)) - dc(i, jnp) = sign(min(abs(tmp),pmax,pmin), tmp) - END DO - - DO i = imh + 1, imr - dc(i, 1) = -dc(i-imh, 1) - dc(i, jnp) = -dc(i-imh, jnp) - END DO - END IF - RETURN -END SUBROUTINE ymist - -SUBROUTINE fyppm(vc, p, dc, flux, imr, jnp, j1, j2, a6, ar, al, jord) - PARAMETER (r3=1./3., r23=2./3.) - REAL vc(imr, *), flux(imr, *), p(imr, *), dc(imr, *) - ! Local work arrays. - REAL ar(imr, jnp), al(imr, jnp), a6(imr, jnp) - INTEGER lmt - ! logical first - ! data first /.true./ - ! SAVE LMT - - imh = imr/2 - jmr = jnp - 1 - j11 = j1 - 1 - imjm1 = imr*(j2-j1+2) - len = imr*(j2-j1+3) - ! if(first) then - ! IF(JORD.LE.0) then - ! if(JMR.GE.90) then - ! LMT = 0 - ! elseif(JMR.GE.45) then - ! LMT = 1 - ! else - ! LMT = 2 - ! endif - ! else - ! LMT = JORD - 3 - ! endif - - ! first = .false. - ! endif - - ! modifs pour pouvoir choisir plusieurs schemas PPM - lmt = jord - 3 - - DO i = 1, imr*jmr - al(i, 2) = 0.5*(p(i,1)+p(i,2)) + (dc(i,1)-dc(i,2))*r3 - ar(i, 1) = al(i, 2) - END DO - - ! Poles: - - DO i = 1, imh - al(i, 1) = al(i+imh, 2) - al(i+imh, 1) = al(i, 2) - - ar(i, jnp) = ar(i+imh, jmr) - ar(i+imh, jnp) = ar(i, jmr) - END DO - - ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - ! Rajout pour LMDZ.3.3 - ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - ar(imr, 1) = al(1, 1) - ar(imr, jnp) = al(1, jnp) - ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - - DO i = 1, len - a6(i, j11) = 3.*(p(i,j11)+p(i,j11)-(al(i,j11)+ar(i,j11))) - END DO - - IF (lmt<=2) CALL lmtppm(dc(1,j11), a6(1,j11), ar(1,j11), al(1,j11), & - p(1,j11), len, lmt) - - - DO i = 1, imjm1 - IF (vc(i,j1)>0.) THEN - flux(i, j1) = ar(i, j11) + 0.5*vc(i, j1)*(al(i,j11)-ar(i,j11)+a6(i,j11) & - *(1.-r23*vc(i,j1))) - ELSE - flux(i, j1) = al(i, j1) - 0.5*vc(i, j1)*(ar(i,j1)-al(i,j1)+a6(i,j1)*(1. & - +r23*vc(i,j1))) - END IF - END DO - RETURN -END SUBROUTINE fyppm - -SUBROUTINE yadv(imr, jnp, j1, j2, p, va, ady, wk, iad) - REAL p(imr, jnp), ady(imr, jnp), va(imr, jnp) - REAL wk(imr, -1:jnp+2) - - jmr = jnp - 1 - imh = imr/2 - DO j = 1, jnp - DO i = 1, imr - wk(i, j) = p(i, j) - END DO - END DO - ! Poles: - DO i = 1, imh - wk(i, -1) = p(i+imh, 3) - wk(i+imh, -1) = p(i, 3) - wk(i, 0) = p(i+imh, 2) - wk(i+imh, 0) = p(i, 2) - wk(i, jnp+1) = p(i+imh, jmr) - wk(i+imh, jnp+1) = p(i, jmr) - wk(i, jnp+2) = p(i+imh, jnp-2) - wk(i+imh, jnp+2) = p(i, jnp-2) - END DO - - IF (iad==2) THEN - DO j = j1 - 1, j2 + 1 - DO i = 1, imr - jp = nint(va(i,j)) - rv = jp - va(i, j) - jp = j - jp - a1 = 0.5*(wk(i,jp+1)+wk(i,jp-1)) - wk(i, jp) - b1 = 0.5*(wk(i,jp+1)-wk(i,jp-1)) - ady(i, j) = wk(i, jp) + rv*(a1*rv+b1) - wk(i, j) - END DO - END DO - - ELSE IF (iad==1) THEN - DO j = j1 - 1, j2 + 1 - DO i = 1, imr - jp = float(j) - va(i, j) - ady(i, j) = va(i, j)*(wk(i,jp)-wk(i,jp+1)) - END DO - END DO - END IF - - IF (j1/=2) THEN - sum1 = 0. - sum2 = 0. - DO i = 1, imr - sum1 = sum1 + ady(i, 2) - sum2 = sum2 + ady(i, jmr) - END DO - sum1 = sum1/imr - sum2 = sum2/imr - - DO i = 1, imr - ady(i, 2) = sum1 - ady(i, jmr) = sum2 - ady(i, 1) = sum1 - ady(i, jnp) = sum2 - END DO - ELSE - ! Poles: - sum1 = 0. - sum2 = 0. - DO i = 1, imr - sum1 = sum1 + ady(i, 1) - sum2 = sum2 + ady(i, jnp) - END DO - sum1 = sum1/imr - sum2 = sum2/imr - - DO i = 1, imr - ady(i, 1) = sum1 - ady(i, jnp) = sum2 - END DO - END IF - - RETURN -END SUBROUTINE yadv - -SUBROUTINE xadv(imr, jnp, j1, j2, p, ua, js, jn, iml, adx, iad) - REAL p(imr, jnp), adx(imr, jnp), qtmp(-imr:imr+imr), ua(imr, jnp) - - jmr = jnp - 1 - DO j = j1, j2 - IF (j>js .AND. j=0.) THEN - adx(i, j) = qtmp(iiu) + ru*(qtmp(iiu-1)-qtmp(iiu)) - ELSE - adx(i, j) = qtmp(iiu) + ru*(qtmp(iiu)-qtmp(iiu+1)) - END IF - END DO - END IF - - DO i = 1, imr - adx(i, j) = adx(i, j) - p(i, j) - END DO -1309 END DO - - ! Eulerian upwind - - DO j = js + 1, jn - 1 - - DO i = 1, imr - qtmp(i) = p(i, j) - END DO - - qtmp(0) = p(imr, j) - qtmp(imr+1) = p(1, j) - - IF (iad==2) THEN - qtmp(-1) = p(imr-1, j) - qtmp(imr+2) = p(2, j) - DO i = 1, imr - ip = nint(ua(i,j)) - ru = ip - ua(i, j) - ip = i - ip - a1 = 0.5*(qtmp(ip+1)+qtmp(ip-1)) - qtmp(ip) - b1 = 0.5*(qtmp(ip+1)-qtmp(ip-1)) - adx(i, j) = qtmp(ip) - p(i, j) + ru*(a1*ru+b1) - END DO - ELSE IF (iad==1) THEN - ! 1st order - DO i = 1, imr - ip = i - ua(i, j) - adx(i, j) = ua(i, j)*(qtmp(ip)-qtmp(ip+1)) - END DO - END IF - END DO - - IF (j1/=2) THEN - DO i = 1, imr - adx(i, 2) = 0. - adx(i, jmr) = 0. - END DO - END IF - ! set cross term due to x-adv at the poles to zero. - DO i = 1, imr - adx(i, 1) = 0. - adx(i, jnp) = 0. - END DO - RETURN -END SUBROUTINE xadv - -SUBROUTINE lmtppm(dc, a6, ar, al, p, im, lmt) - - ! A6 = CURVATURE OF THE TEST PARABOLA - ! AR = RIGHT EDGE VALUE OF THE TEST PARABOLA - ! AL = LEFT EDGE VALUE OF THE TEST PARABOLA - ! DC = 0.5 * MISMATCH - ! P = CELL-AVERAGED VALUE - ! IM = VECTOR LENGTH - - ! OPTIONS: - - ! LMT = 0: FULL MONOTONICITY - ! LMT = 1: SEMI-MONOTONIC CONSTRAINT (NO UNDERSHOOTS) - ! LMT = 2: POSITIVE-DEFINITE CONSTRAINT - - PARAMETER (r12=1./12.) - DIMENSION a6(im), ar(im), al(im), p(im), dc(im) - - IF (lmt==0) THEN - ! Full constraint - DO i = 1, im - IF (dc(i)==0.) THEN - ar(i) = p(i) - al(i) = p(i) - a6(i) = 0. - ELSE - da1 = ar(i) - al(i) - da2 = da1**2 - a6da = a6(i)*da1 - IF (a6da<-da2) THEN - a6(i) = 3.*(al(i)-p(i)) - ar(i) = al(i) - a6(i) - ELSE IF (a6da>da2) THEN - a6(i) = 3.*(ar(i)-p(i)) - al(i) = ar(i) - a6(i) - END IF - END IF - END DO - ELSE IF (lmt==1) THEN - ! Semi-monotonic constraint - DO i = 1, im - IF (abs(ar(i)-al(i))>=-a6(i)) GO TO 150 - IF (p(i)al(i)) THEN - a6(i) = 3.*(al(i)-p(i)) - ar(i) = al(i) - a6(i) - ELSE - a6(i) = 3.*(ar(i)-p(i)) - al(i) = ar(i) - a6(i) - END IF -150 END DO - ELSE IF (lmt==2) THEN - DO i = 1, im - IF (abs(ar(i)-al(i))>=-a6(i)) GO TO 250 - fmin = p(i) + 0.25*(ar(i)-al(i))**2/a6(i) + a6(i)*r12 - IF (fmin>=0.) GO TO 250 - IF (p(i)al(i)) THEN - a6(i) = 3.*(al(i)-p(i)) - ar(i) = al(i) - a6(i) - ELSE - a6(i) = 3.*(ar(i)-p(i)) - al(i) = ar(i) - a6(i) - END IF -250 END DO - END IF - RETURN -END SUBROUTINE lmtppm - -SUBROUTINE a2c(u, v, imr, jmr, j1, j2, crx, cry, dtdx5, dtdy5) - DIMENSION u(imr, *), v(imr, *), crx(imr, *), cry(imr, *), dtdx5(*) - - DO j = j1, j2 - DO i = 2, imr - crx(i, j) = dtdx5(j)*(u(i,j)+u(i-1,j)) - END DO - END DO - - DO j = j1, j2 - crx(1, j) = dtdx5(j)*(u(1,j)+u(imr,j)) - END DO - - DO i = 1, imr*jmr - cry(i, 2) = dtdy5*(v(i,2)+v(i,1)) - END DO - RETURN -END SUBROUTINE a2c - -SUBROUTINE cosa(cosp, cose, jnp, pi, dp) - DIMENSION cosp(*), cose(*) - - jmr = jnp - 1 - DO j = 2, jnp - ph5 = -0.5*pi + (float(j-1)-0.5)*dp - cose(j) = cos(ph5) - END DO - - jeq = (jnp+1)/2 - IF (jmr==2*(jmr/2)) THEN - DO j = jnp, jeq + 1, -1 - cose(j) = cose(jnp+2-j) - END DO - ELSE - ! cell edge at equator. - cose(jeq+1) = 1. - DO j = jnp, jeq + 2, -1 - cose(j) = cose(jnp+2-j) - END DO - END IF - - DO j = 2, jmr - cosp(j) = 0.5*(cose(j)+cose(j+1)) - END DO - cosp(1) = 0. - cosp(jnp) = 0. - RETURN -END SUBROUTINE cosa - -SUBROUTINE cosc(cosp, cose, jnp, pi, dp) - DIMENSION cosp(*), cose(*) - - phi = -0.5*pi - DO j = 2, jnp - 1 - phi = phi + dp - cosp(j) = cos(phi) - END DO - cosp(1) = 0. - cosp(jnp) = 0. - - DO j = 2, jnp - cose(j) = 0.5*(cosp(j)+cosp(j-1)) - END DO - - DO j = 2, jnp - 1 - cosp(j) = 0.5*(cose(j)+cose(j+1)) - END DO - RETURN -END SUBROUTINE cosc - -SUBROUTINE qckxyz(q, qtmp, imr, jnp, nlay, j1, j2, cosp, acosp, cross, ic, & - nstep) - - PARAMETER (tiny=1.E-60) - DIMENSION q(imr, jnp, nlay), qtmp(imr, jnp), cosp(*), acosp(*) - LOGICAL cross - - nlaym1 = nlay - 1 - len = imr*(j2-j1+1) - ip = 0 - - ! Top layer - l = 1 - icr = 1 - CALL filns(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, ipy, tiny) - IF (ipy==0) GO TO 50 - CALL filew(q(1,1,l), qtmp, imr, jnp, j1, j2, ipx, tiny) - IF (ipx==0) GO TO 50 - - IF (cross) THEN - CALL filcr(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, icr, tiny) - END IF - IF (icr==0) GO TO 50 - - ! Vertical filling... - DO i = 1, len - IF (q(i,j1,1)<0.) THEN - ip = ip + 1 - q(i, j1, 2) = q(i, j1, 2) + q(i, j1, 1) - q(i, j1, 1) = 0. - END IF - END DO - -50 CONTINUE - DO l = 2, nlaym1 - icr = 1 - - CALL filns(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, ipy, tiny) - IF (ipy==0) GO TO 225 - CALL filew(q(1,1,l), qtmp, imr, jnp, j1, j2, ipx, tiny) - IF (ipx==0) GO TO 225 - IF (cross) THEN - CALL filcr(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, icr, tiny) - END IF - IF (icr==0) GO TO 225 - - DO i = 1, len - IF (q(i,j1,l)<0.) THEN - - ip = ip + 1 - ! From above - qup = q(i, j1, l-1) - qly = -q(i, j1, l) - dup = min(qly, qup) - q(i, j1, l-1) = qup - dup - q(i, j1, l) = dup - qly - ! Below - q(i, j1, l+1) = q(i, j1, l+1) + q(i, j1, l) - q(i, j1, l) = 0. - END IF - END DO -225 END DO - - ! BOTTOM LAYER - sum = 0. - l = nlay - - CALL filns(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, ipy, tiny) - IF (ipy==0) GO TO 911 - CALL filew(q(1,1,l), qtmp, imr, jnp, j1, j2, ipx, tiny) - IF (ipx==0) GO TO 911 - - CALL filcr(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, icr, tiny) - IF (icr==0) GO TO 911 - - DO i = 1, len - IF (q(i,j1,l)<0.) THEN - ip = ip + 1 - - ! From above - - qup = q(i, j1, nlaym1) - qly = -q(i, j1, l) - dup = min(qly, qup) - q(i, j1, nlaym1) = qup - dup - ! From "below" the surface. - sum = sum + qly - dup - q(i, j1, l) = 0. - END IF - END DO - -911 CONTINUE - - IF (ip>imr) THEN - WRITE (6, *) 'IC=', ic, ' STEP=', nstep, ' Vertical filling pts=', ip - END IF - - IF (sum>1.E-25) THEN - WRITE (6, *) ic, nstep, ' Mass source from the ground=', sum - END IF - RETURN -END SUBROUTINE qckxyz - -SUBROUTINE filcr(q, imr, jnp, j1, j2, cosp, acosp, icr, tiny) - DIMENSION q(imr, *), cosp(*), acosp(*) - - icr = 0 - DO j = j1 + 1, j2 - 1 - DO i = 1, imr - 1 - IF (q(i,j)<0.) THEN - icr = 1 - dq = -q(i, j)*cosp(j) - ! N-E - dn = q(i+1, j+1)*cosp(j+1) - d0 = max(0., dn) - d1 = min(dq, d0) - q(i+1, j+1) = (dn-d1)*acosp(j+1) - dq = dq - d1 - ! S-E - ds = q(i+1, j-1)*cosp(j-1) - d0 = max(0., ds) - d2 = min(dq, d0) - q(i+1, j-1) = (ds-d2)*acosp(j-1) - q(i, j) = (d2-dq)*acosp(j) + tiny - END IF - END DO - IF (icr==0 .AND. q(imr,j)>=0.) GO TO 65 - DO i = 2, imr - IF (q(i,j)<0.) THEN - icr = 1 - dq = -q(i, j)*cosp(j) - ! N-W - dn = q(i-1, j+1)*cosp(j+1) - d0 = max(0., dn) - d1 = min(dq, d0) - q(i-1, j+1) = (dn-d1)*acosp(j+1) - dq = dq - d1 - ! S-W - ds = q(i-1, j-1)*cosp(j-1) - d0 = max(0., ds) - d2 = min(dq, d0) - q(i-1, j-1) = (ds-d2)*acosp(j-1) - q(i, j) = (d2-dq)*acosp(j) + tiny - END IF - END DO - ! ***************************************** - ! i=1 - i = 1 - IF (q(i,j)<0.) THEN - icr = 1 - dq = -q(i, j)*cosp(j) - ! N-W - dn = q(imr, j+1)*cosp(j+1) - d0 = max(0., dn) - d1 = min(dq, d0) - q(imr, j+1) = (dn-d1)*acosp(j+1) - dq = dq - d1 - ! S-W - ds = q(imr, j-1)*cosp(j-1) - d0 = max(0., ds) - d2 = min(dq, d0) - q(imr, j-1) = (ds-d2)*acosp(j-1) - q(i, j) = (d2-dq)*acosp(j) + tiny - END IF - ! ***************************************** - ! i=IMR - i = imr - IF (q(i,j)<0.) THEN - icr = 1 - dq = -q(i, j)*cosp(j) - ! N-E - dn = q(1, j+1)*cosp(j+1) - d0 = max(0., dn) - d1 = min(dq, d0) - q(1, j+1) = (dn-d1)*acosp(j+1) - dq = dq - d1 - ! S-E - ds = q(1, j-1)*cosp(j-1) - d0 = max(0., ds) - d2 = min(dq, d0) - q(1, j-1) = (ds-d2)*acosp(j-1) - q(i, j) = (d2-dq)*acosp(j) + tiny - END IF - ! ***************************************** -65 END DO - - DO i = 1, imr - IF (q(i,j1)<0. .OR. q(i,j2)<0.) THEN - icr = 1 - GO TO 80 - END IF - END DO - -80 CONTINUE - - IF (q(1,1)<0. .OR. q(1,jnp)<0.) THEN - icr = 1 - END IF - - RETURN -END SUBROUTINE filcr - -SUBROUTINE filns(q, imr, jnp, j1, j2, cosp, acosp, ipy, tiny) - DIMENSION q(imr, *), cosp(*), acosp(*) - ! logical first - ! data first /.true./ - ! save cap1 - - ! if(first) then - dp = 4.*atan(1.)/float(jnp-1) - cap1 = imr*(1.-cos((j1-1.5)*dp))/dp - ! first = .false. - ! endif - - ipy = 0 - DO j = j1 + 1, j2 - 1 - DO i = 1, imr - IF (q(i,j)<0.) THEN - ipy = 1 - dq = -q(i, j)*cosp(j) - ! North - dn = q(i, j+1)*cosp(j+1) - d0 = max(0., dn) - d1 = min(dq, d0) - q(i, j+1) = (dn-d1)*acosp(j+1) - dq = dq - d1 - ! South - ds = q(i, j-1)*cosp(j-1) - d0 = max(0., ds) - d2 = min(dq, d0) - q(i, j-1) = (ds-d2)*acosp(j-1) - q(i, j) = (d2-dq)*acosp(j) + tiny - END IF - END DO - END DO - - DO i = 1, imr - IF (q(i,j1)<0.) THEN - ipy = 1 - dq = -q(i, j1)*cosp(j1) - ! North - dn = q(i, j1+1)*cosp(j1+1) - d0 = max(0., dn) - d1 = min(dq, d0) - q(i, j1+1) = (dn-d1)*acosp(j1+1) - q(i, j1) = (d1-dq)*acosp(j1) + tiny - END IF - END DO - - j = j2 - DO i = 1, imr - IF (q(i,j)<0.) THEN - ipy = 1 - dq = -q(i, j)*cosp(j) - ! South - ds = q(i, j-1)*cosp(j-1) - d0 = max(0., ds) - d2 = min(dq, d0) - q(i, j-1) = (ds-d2)*acosp(j-1) - q(i, j) = (d2-dq)*acosp(j) + tiny - END IF - END DO - - ! Check Poles. - IF (q(1,1)<0.) THEN - dq = q(1, 1)*cap1/float(imr)*acosp(j1) - DO i = 1, imr - q(i, 1) = 0. - q(i, j1) = q(i, j1) + dq - IF (q(i,j1)<0.) ipy = 1 - END DO - END IF - - IF (q(1,jnp)<0.) THEN - dq = q(1, jnp)*cap1/float(imr)*acosp(j2) - DO i = 1, imr - q(i, jnp) = 0. - q(i, j2) = q(i, j2) + dq - IF (q(i,j2)<0.) ipy = 1 - END DO - END IF - - RETURN -END SUBROUTINE filns - -SUBROUTINE filew(q, qtmp, imr, jnp, j1, j2, ipx, tiny) - DIMENSION q(imr, *), qtmp(jnp, imr) - - ipx = 0 - ! Copy & swap direction for vectorization. - DO i = 1, imr - DO j = j1, j2 - qtmp(j, i) = q(i, j) - END DO - END DO - - DO i = 2, imr - 1 - DO j = j1, j2 - IF (qtmp(j,i)<0.) THEN - ipx = 1 - ! west - d0 = max(0., qtmp(j,i-1)) - d1 = min(-qtmp(j,i), d0) - qtmp(j, i-1) = qtmp(j, i-1) - d1 - qtmp(j, i) = qtmp(j, i) + d1 - ! east - d0 = max(0., qtmp(j,i+1)) - d2 = min(-qtmp(j,i), d0) - qtmp(j, i+1) = qtmp(j, i+1) - d2 - qtmp(j, i) = qtmp(j, i) + d2 + tiny - END IF - END DO - END DO - - i = 1 - DO j = j1, j2 - IF (qtmp(j,i)<0.) THEN - ipx = 1 - ! west - d0 = max(0., qtmp(j,imr)) - d1 = min(-qtmp(j,i), d0) - qtmp(j, imr) = qtmp(j, imr) - d1 - qtmp(j, i) = qtmp(j, i) + d1 - ! east - d0 = max(0., qtmp(j,i+1)) - d2 = min(-qtmp(j,i), d0) - qtmp(j, i+1) = qtmp(j, i+1) - d2 - - qtmp(j, i) = qtmp(j, i) + d2 + tiny - END IF - END DO - i = imr - DO j = j1, j2 - IF (qtmp(j,i)<0.) THEN - ipx = 1 - ! west - d0 = max(0., qtmp(j,i-1)) - d1 = min(-qtmp(j,i), d0) - qtmp(j, i-1) = qtmp(j, i-1) - d1 - qtmp(j, i) = qtmp(j, i) + d1 - ! east - d0 = max(0., qtmp(j,1)) - d2 = min(-qtmp(j,i), d0) - qtmp(j, 1) = qtmp(j, 1) - d2 - - qtmp(j, i) = qtmp(j, i) + d2 + tiny - END IF - END DO - - IF (ipx/=0) THEN - DO j = j1, j2 - DO i = 1, imr - q(i, j) = qtmp(j, i) - END DO - END DO - ELSE - - ! Poles. - IF (q(1,1)<0 .OR. q(1,jnp)<0.) ipx = 1 - END IF - RETURN -END SUBROUTINE filew