!

! $Id: ppm3d.F 2197 2015-02-09 07:13:05Z emillour $

!

cFrom lin@explorer.gsfc.nasa.gov Wed Apr 15 17:44:44 1998

cDate: Wed, 15 Apr 1998 11:37:03 -0400

cFrom: lin@explorer.gsfc.nasa.gov

cTo: Frederic.Hourdin@lmd.jussieu.fr

cSubject: 3D transport module of the GSFC CTM and GEOS GCM

cThis code is sent to you by S-J Lin, DAO, NASA-GSFC

cNote: this version is intended for machines like CRAY

C-90. No multitasking directives implemented.

C ********************************************************************

C

C TransPort Core for Goddard Chemistry Transport Model (G-CTM), Goddard

C Earth Observing System General Circulation Model (GEOS-GCM), and Data

C Assimilation System (GEOS-DAS).

C

C ********************************************************************

C

C Purpose: given horizontal winds on  a hybrid sigma-p surfaces,

C          one call to tpcore updates the 3-D mixing ratio

C          fields one time step (NDT). [vertical mass flux is computed

C          internally consistent with the discretized hydrostatic mass

C          continuity equation of the C-Grid GEOS-GCM (for IGD=1)].

C

C Schemes: Multi-dimensional Flux Form Semi-Lagrangian (FFSL) scheme based

C          on the van Leer or PPM.

C          (see Lin and Rood 1996).

C Version 4.5

C Last modified: Dec. 5, 1996

C Major changes from version 4.0: a more general vertical hybrid sigma-

C pressure coordinate.

C Subroutines modified: xtp, ytp, fzppm, qckxyz

C Subroutines deleted: vanz

C

C Author: Shian-Jiann Lin

C mail address:

C                 Shian-Jiann Lin*

C                 Code 910.3, NASA/GSFC, Greenbelt, MD 20771

C                 Phone: 301-286-9540

C                 E-mail: lin@dao.gsfc.nasa.gov

C

C *affiliation:

C                 Joint Center for Earth Systems Technology

C                 The University of Maryland Baltimore County

C                 NASA - Goddard Space Flight Center

C References:

C

C 1. Lin, S.-J., and R. B. Rood, 1996: Multidimensional flux form semi-

C    Lagrangian transport schemes. Mon. Wea. Rev., 124, 2046-2070.

C

C 2. Lin, S.-J., W. C. Chao, Y. C. Sud, and G. K. Walker, 1994: A class of

C    the van Leer-type transport schemes and its applications to the moist-

C    ure transport in a General Circulation Model. Mon. Wea. Rev., 122,

C    1575-1593.

C

C ****6***0*********0*********0*********0*********0*********0**********72

C

      implicit none

c     rajout de d�clarations

c      integer Jmax,kmax,ndt0,nstep,k,j,i,ic,l,js,jn,imh,iad,jad,krd

c      integer iu,iiu,j2,jmr,js0,jt

c      real dtdy,dtdy5,rcap,iml,jn0,imjm,pi,dl,dp

c      real dt,cr1,maxdt,ztc,d5,sum1,sum2,ru

C

C ********************************************************************

C

C =============

C INPUT:

C =============

C

C Q(IMR,JNP,NLAY,NC): mixing ratios at current time (t)

C NC: total # of constituents

C IMR: first dimension (E-W); # of Grid intervals in E-W is IMR

C JNP: 2nd dimension (N-S); # of Grid intervals in N-S is JNP-1

C NLAY: 3rd dimension (# of layers); vertical index increases from 1 at

C       the model top to NLAY near the surface (see fig. below).

C       It is assumed that 6 <= NLAY <= JNP (for dynamic memory allocation)

C

C PS1(IMR,JNP): surface pressure at current time (t)

C PS2(IMR,JNP): surface pressure at mid-time-level (t+NDT/2)

C PS2 is replaced by the predicted PS (at t+NDT) on output.

C Note: surface pressure can have any unit or can be multiplied by any

C       const.

C

C The pressure at layer edges are defined as follows:

C

C        p(i,j,k) = AP(k)*PT  +  BP(k)*PS(i,j)          (1)

C

C Where PT is a constant having the same unit as PS.

C AP and BP are unitless constants given at layer edges

C defining the vertical coordinate.

C BP(1) = 0., BP(NLAY+1) = 1.

C The pressure at the model top is PTOP = AP(1)*PT

C

C For pure sigma system set AP(k) = 1 for all k, PT = PTOP,

C BP(k) = sige(k) (sigma at edges), PS = Psfc - PTOP.

C

C Note: the sigma-P coordinate is a subset of Eq. 1, which in turn

C is a subset of the following even more general sigma-P-thelta coord.

C currently under development.

C  p(i,j,k) = (AP(k)*PT + BP(k)*PS(i,j))/(D(k)-C(k)*TE**(-1/kapa))

C

C                  /////////////////////////////////

C              / \ ------------- PTOP --------------  AP(1), BP(1)

C               |

C    delp(1)    |  ........... Q(i,j,1) ............

C               |

C      W(1)    \ / ---------------------------------  AP(2), BP(2)

C

C

C

C     W(k-1)   / \ ---------------------------------  AP(k), BP(k)

C               |

C    delp(K)    |  ........... Q(i,j,k) ............

C               |

C      W(k)    \ / ---------------------------------  AP(k+1), BP(k+1)

C

C

C

C              / \ ---------------------------------  AP(NLAY), BP(NLAY)

C               |

C  delp(NLAY)   |  ........... Q(i,j,NLAY) .........

C               |

C   W(NLAY)=0  \ / ------------- surface ----------- AP(NLAY+1), BP(NLAY+1)

C                 //////////////////////////////////

C

C U(IMR,JNP,NLAY) & V(IMR,JNP,NLAY):winds (m/s) at mid-time-level (t+NDT/2)

C U and V may need to be polar filtered in advance in some cases.

C

C IGD:      grid type on which winds are defined.

C IGD = 0:  A-Grid  [all variables defined at the same point from south

C                   pole (j=1) to north pole (j=JNP) ]

C

C IGD = 1  GEOS-GCM C-Grid

C                                      [North]

C

C                                       V(i,j)

C                                          |

C                                          |

C                                          |

C                             U(i-1,j)---Q(i,j)---U(i,j) [EAST]

C                                          |

C                                          |

C                                          |

C                                       V(i,j-1)

C

C         U(i,  1) is defined at South Pole.

C         V(i,  1) is half grid north of the South Pole.

C         V(i,JMR) is half grid south of the North Pole.

C

C         V must be defined at j=1 and j=JMR if IGD=1

C         V at JNP need not be given.

C

C NDT: time step in seconds (need not be constant during the course of

C      the integration). Suggested value: 30 min. for 4x5, 15 min. for 2x2.5

C      (Lat-Lon) resolution. Smaller values are recommanded if the model

C      has a well-resolved stratosphere.

C

C J1 defines the size of the polar cap:

C South polar cap edge is located at -90 + (j1-1.5)*180/(JNP-1) deg.

C North polar cap edge is located at  90 - (j1-1.5)*180/(JNP-1) deg.

C There are currently only two choices (j1=2 or 3).

C IMR must be an even integer if j1 = 2. Recommended value: J1=3.

C

C IORD, JORD, and KORD are integers controlling various options in E-W, N-S,

C and vertical transport, respectively. Recommended values for positive

C definite scalars: IORD=JORD=3, KORD=5. Use KORD=3 for non-

C positive definite scalars or when linear correlation between constituents

C is to be maintained.

C

C  _ORD=

C        1: 1st order upstream scheme (too diffusive, not a useful option; it

C           can be used for debugging purposes; this is THE only known "linear"

C           monotonic advection scheme.).

C        2: 2nd order van Leer (full monotonicity constraint;

C           see Lin et al 1994, MWR)

C        3: monotonic PPM* (slightly improved PPM of Collela & Woodward 1984)

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 & undershoots 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 *PPM: Piece-wise Parabolic Method

C

C Note that KORD <=2 options are no longer supported. DO not use option 4 or 5.

C for non-positive definite scalars (such as Ertel Potential Vorticity).

C

C The implicit numerical diffusion decreases as _ORD increases.

C The last two options (ORDER=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 No filter of any kind is applied to the constituent fields here.

C

C AE: Radius of the sphere (meters).

C     Recommended value for the planet earth: 6.371E6

C

C fill(logical):   flag to do filling for negatives (see note below).

C

C Umax: Estimate (upper limit) of the maximum U-wind speed (m/s).

C (220 m/s is a good value for troposphere model; 280 m/s otherwise)

C

C =============

C Output

C =============

C

C Q: mixing ratios at future time (t+NDT) (original values are over-written)

C W(NLAY): large-scale vertical mass flux as diagnosed from the hydrostatic

C          relationship. W will have the same unit as PS1 and PS2 (eg, mb).

C          W must be divided by NDT to get the correct mass-flux unit.

C          The vertical Courant number C = W/delp_UPWIND, where delp_UPWIND

C          is the pressure thickness in the "upwind" direction. For example,

C          C(k) = W(k)/delp(k)   if W(k) > 0;

C          C(k) = W(k)/delp(k+1) if W(k) < 0.

C              ( W > 0 is downward, ie, toward surface)

C PS2: predicted PS at t+NDT (original values are over-written)

C

C ********************************************************************

C NOTES:

C This forward-in-time upstream-biased transport scheme reduces to

C the 2nd order center-in-time center-in-space mass continuity eqn.

C if Q = 1 (constant fields will remain constant). This also ensures

C that the computed vertical velocity to be identical to GEOS-1 GCM

C for on-line transport.

C

C A larger polar cap is used if j1=3 (recommended for C-Grid winds or when

C winds are noisy near poles).

C

C Flux-Form Semi-Lagrangian transport in the East-West direction is used

C when and where Courant # is greater than one.

C

C The user needs to change the parameter Jmax or Kmax if the resolution

C is greater than 0.5 deg in N-S or 150 layers in the vertical direction.

C (this TransPort Core is otherwise resolution independent and can be used

C as a library routine).

C

C PPM is 4th order accurate when grid spacing is uniform (x & y); 3rd

C order accurate for non-uniform grid (vertical sigma coord.).

C

C Time step is limitted only by transport in the meridional direction.

C (the FFSL scheme is not implemented in the meridional direction).

C

C Since only 1-D limiters are applied, negative values could

C potentially be generated when large time step is used and when the

C initial fields contain discontinuities.

C This does not necessarily imply the integration is unstable.

C These negatives are typically very small. A filling algorithm is

C activated if the user set "fill" to be true.

C

C The van Leer scheme used here is nearly as accurate as the original PPM

C due to the use of a 4th order accurate reference slope. The PPM imple-

C mented here is an improvement over the original and is also based on

C the 4th order reference slope.

C

C ****6***0*********0*********0*********0*********0*********0**********72

C

C     User modifiable parameters

C

      integer,parameter :: Jmax = 361, kmax = 150

C

C ****6***0*********0*********0*********0*********0*********0**********72

C

C Input-Output arrays

C

      real Q(IMR,JNP,NLAY,NC),PS1(IMR,JNP),PS2(IMR,JNP),

     &     U(IMR,JNP,NLAY),V(IMR,JNP,NLAY),AP(NLAY+1),

     &     BP(NLAY+1),W(IMR,JNP,NLAY),NDT,val(NLAY),Umax

      integer IGD,IORD,JORD,KORD,NC,IMR,JNP,j1,NLAY,AE

      integer IMRD2

      real    PT

      logical  cross, fill, dum

C

C Local dynamic arrays

C

C

C Local static  arrays

C

      real DTDX(Jmax), DTDX5(Jmax), acosp(Jmax),

     &     cosp(Jmax), cose(Jmax), DAP(kmax),DBK(Kmax)

      data NDT0, NSTEP /0, 0/

      data cross /.true./

      REAL DTDY, DTDY5, RCAP

      INTEGER JS0, JN0, IML, JMR, IMJM

      SAVE DTDY, DTDY5, RCAP, JS0, JN0, IML,

     &     DTDX, DTDX5, ACOSP, COSP, COSE, DAP,DBK

C

      INTEGER NDT0, NSTEP, j2, k,j,i,ic,l,JS,JN,IMH

      INTEGER IU,IIU,JT,iad,jad,krd

      REAL r23,r3,PI,DL,DP,DT,CR1,MAXDT,ZTC,D5

      REAL sum1,sum2,ru

C

C *********** Initialization **********************

c

c

C

C controles sur les parametres

      endif

      endif

      endif

C

      endif

C

      ENDDO

C

C

C Compute analytic cosine at cell edges

      else

C Define cosine consistent with GEOS-GCM (using dycore2.0 or later)

      endif

C

C

C Inverse of the Scaled polar cap area.

C

      endif

C

      DT   = NDT

	endif

      endif

C

      ZTC = 0.

      else

C

      endif

C

C

c     print*,'dtdx=',dtdx(j)

      enddo

C

c      print*,'dtdy=',dtdy

C

c      write(6,*) 'J1=',J1,' J2=', J2

      endif

C

C *********** End Initialization **********************

C

C delp = pressure thickness: the psudo-density in a hydrostatic system.

            enddo

         enddo

      enddo

C

      endif

C

C Compute "tracer density"

C

C

C Convert winds on A-Grid to Courant # on C-Grid.

      else

C Convert winds on C-grid to Courant #

C

C

      endif

C

C Determine JS and JN

C

      endif

      enddo

      enddo

C

2222  continue

      endif

      enddo

      enddo

2233  continue

C

      enddo

      endif

C

C ******* Compute horizontal mass fluxes ************

C

C N-S component

      enddo

      enddo

C

C

C Poles

      enddo

C

      enddo

C

C E-W component

C

      enddo

      enddo

C

      enddo

C

C

C

C

      enddo

      enddo

C

      enddo

ccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c Rajouts pour LMDZ.3.3

ccccccccccccccccccccccccccccccccccccccccccccccccccccccc

         enddo

      enddo

      enddo

C

      enddo

      endif

C

C ****6***0*********0*********0*********0*********0*********0**********72

C

      enddo

C

C E-W advective cross term

C

      enddo

C

      enddo

C

      else

      endif

C

C

      enddo

C

C

      enddo

      enddo

      endif

C ****6***0*********0*********0*********0*********0*********0**********72

C Contribution from the N-S advection

      enddo

C

      enddo

C

C Add cross terms in the vertical direction.

	else

	endif

C

	else

	endif

      enddo

      enddo

      endif

C

      call xtp(IMR,JNP,IML,j1,j2,JN,JS,PU,DQ(1,1,k,IC),wk1(1,1,2)

      call ytp(IMR,JNP,j1,j2,acosp,RCAP,DQ(1,1,k,IC),wk1(1,1,1),CRY,

C

C

C ******* Compute vertical mass flux (same unit as PS) ***********

C

C 1st step: compute total column mass CONVERGENCE.

C

C

C

C

C 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption.

C Changes (increases) to surface pressure = total column mass convergence

C

C

C 3rd step: compute vertical mass flux from mass conservation principle.

C

C

C

C

C

C****6***0*********0*********0*********0*********0*********0**********72

      call FZPPM(IMR,JNP,NLAY,j1,DQ(1,1,1,IC),W,Q(1,1,1,IC),WK1,DPI,

C

C

C Recover tracer mixing ratio from "density" using predicted

C "air density" (pressure thickness) at time-level n+1

C

c            print*,'i=',i,'j=',j,'k=',k,'Q(i,j,k,IC)=',Q(i,j,k,IC)

      enddo

      enddo

      enddo

C

c     j=1 c'est le p�le Sud, j=JNP c'est le p�le Nord

      endif

C

      endif

C

      END

C

C****6***0*********0*********0*********0*********0*********0**********72

     &                 flux,wk1,wk2,wz2,delp,KORD)

      implicit none

      integer,parameter :: kmax = 150

      real,parameter :: R23 = 2./3., R3 = 1./3.

      integer IMR,JNP,NLAY,J1,KORD

      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)

C Assuming JNP >= NLAY

      real AR(IMR,*),AL(IMR,*),A6(IMR,*),flux(IMR,*),wk2(IMR,*),

     &     wz2(IMR,*)

      integer JMR,IMJM,NLAYM1,LMT,K,I,J

      real c0,c1,c2,tmp,qmax,qmin,a,b,fct,a1,a2,cm,cp

C

C

C

C ****6***0*********0*********0*********0*********0*********0**********72

C Compute DC for PPM

C ****6***0*********0*********0*********0*********0*********0**********72

C

C

C

C ****6***0*********0*********0*********0*********0*********0**********72

C Loop over latitudes  (to save memory)

C ****6***0*********0*********0*********0*********0*********0**********72

C

C

      enddo

      enddo

C

C****6***0*********0*********0*********0*********0*********0**********72

C Compute first guesses at cell interfaces

C First guesses are required to be continuous.

C ****6***0*********0*********0*********0*********0*********0**********72

C

C three-cell parabolic subgrid distribution at model top

C two-cell parabolic with zero gradient subgrid distribution

C at the surface.

C

C First guess top edge value

C three-cell PPM

C 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)) ) /

      b = 2.*DQDT(i,j,1)/(wk2(i,1)+wk2(i,2)) -

C

C Check if change sign

	else

	endif

C

C Bottom

C 2-cell PPM with zero gradient right at the surface

C

      fct = DQDT(i,j,NLAYM1)*wk2(i,NLAY)**2 /

C

C****6***0*********0*********0*********0*********0*********0**********72

C 4th order interpolation in the interior.

C****6***0*********0*********0*********0*********0*********0**********72

C

      AL(i,k) = wk1(i,k-1) + c1 + c2 *

     &        ( wk2(i,k  )*(c1*(A1 - A2)+A2*flux(i,k-1)) -

C      print *,'AL1',i,k, AL(i,k)

C

C      print *,'AR1',i,AR(i,1)

C

C      print *,'A61',i,A6(i,1)

C

C****6***0*********0*********0*********0*********0*********0**********72

C Top & Bot always monotonic

      call lmtppm(flux(1,NLAY),A6(1,NLAY),AR(1,NLAY),AL(1,NLAY),

C

C Interior depending on KORD

     &  call lmtppm(flux(1,2),A6(1,2),AR(1,2),AL(1,2),wk1(1,2),

C

C****6***0*********0*********0*********0*********0*********0**********72

C

      else

C        print *,'test2-0',i,j,wz2(i,1),wk2(i,2)

C        print *,'testCP',CP

C        print *,'test2',i, AL(i,2),AR(i,2),A6(i,2),R23

      endif

C

C

C

      end

C

      implicit none

      integer IMR,JNP,IML,j1,j2,JN,JS,IORD

      real PU,DQ,Q,UC,fx1,xmass

      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)

      integer jvan,j1vl,j2vl,j,i,iu,itmp,ist,imp

      real rut

C

      IMP = IMR + 1

C

C van Leer at high latitudes

C

C

      enddo

C

C ************* Eulerian **********

C

C

      ELSE

C

      else

      endif

C

      ENDIF

C

C

C

C ***** Conservative (flux-form) Semi-Lagrangian transport *****

C

2222  continue

C