/[lmdze]/trunk/Sources/dyn3d/PPM3d/ppm3d.f
ViewVC logotype

Diff of /trunk/Sources/dyn3d/PPM3d/ppm3d.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 166 by guez, Wed Jul 29 14:32:55 2015 UTC revision 169 by guez, Mon Sep 14 17:13:16 2015 UTC
# Line 1  Line 1 
1  SUBROUTINE ppm3d(igd, q, ps1, ps2, u, v, w, ndt, iord, jord, kord, nc, imr, &  module ppm3d_m
     jnp, j1, nlay, ap, bp, pt, ae, fill, dum, umax)  
2    
3    ! implicit none    implicit none
4    
5    ! rajout de déclarations  contains
   ! integer Jmax,kmax,ndt0,nstep,k,j,i,ic,l,js,jn,imh,iad,jad,krd  
   ! integer iu,iiu,j2,jmr,js0,jt  
   ! real dtdy,dtdy5,rcap,iml,jn0,imjm,pi,dl,dp  
   ! real dt,cr1,maxdt,ztc,d5,sum1,sum2,ru  
   
   ! ********************************************************************  
   
   ! =============  
   ! INPUT:  
   ! =============  
   
   ! Q(IMR,JNP,NLAY,NC): mixing ratios at current time (t)  
   ! NC: total number of constituents  
   ! IMR: first dimension (E-W); number of Grid intervals in E-W is IMR  
   ! JNP: 2nd dimension (N-S); number of Grid intervals in N-S is JNP-1  
   ! NLAY: 3rd dimension (number of layers); vertical index increases from 1  
   ! at  
   ! the model top to NLAY near the surface (see fig. below).  
   ! It is assumed that 6 <= NLAY <= JNP (for dynamic memory allocation)  
   
   ! PS1(IMR,JNP): surface pressure at current time (t)  
   ! PS2(IMR,JNP): surface pressure at mid-time-level (t+NDT/2)  
   ! PS2 is replaced by the predicted PS (at t+NDT) on output.  
   ! Note: surface pressure can have any unit or can be multiplied by any  
   ! const.  
   
   ! The pressure at layer edges are defined as follows:  
   
   ! p(i,j,k) = AP(k)*PT  +  BP(k)*PS(i,j)          (1)  
   
   ! Where PT is a constant having the same unit as PS.  
   ! AP and BP are unitless constants given at layer edges  
   ! defining the vertical coordinate.  
   ! BP(1) = 0., BP(NLAY+1) = 1.  
   ! The pressure at the model top is PTOP = AP(1)*PT  
   
   ! For pure sigma system set AP(k) = 1 for all k, PT = PTOP,  
   ! BP(k) = sige(k) (sigma at edges), PS = Psfc - PTOP.  
   
   ! Note: the sigma-P coordinate is a subset of Eq. 1, which in turn  
   ! is a subset of the following even more general sigma-P-thelta coord.  
   ! currently under development.  
   ! p(i,j,k) = (AP(k)*PT + BP(k)*PS(i,j))/(D(k)-C(k)*TE**(-1/kapa))  
   
   ! /////////////////////////////////  
   ! / \ ------------- PTOP --------------  AP(1), BP(1)  
   ! |  
   ! delp(1)    |  ........... Q(i,j,1) ............  
   ! |  
   ! W(1)    \ / ---------------------------------  AP(2), BP(2)  
   
   
   
   ! W(k-1)   / \ ---------------------------------  AP(k), BP(k)  
   ! |  
   ! delp(K)    |  ........... Q(i,j,k) ............  
   ! |  
   ! W(k)    \ / ---------------------------------  AP(k+1), BP(k+1)  
   
   
   
   ! / \ ---------------------------------  AP(NLAY), BP(NLAY)  
   ! |  
   ! delp(NLAY)   |  ........... Q(i,j,NLAY) .........  
   ! |  
   ! W(NLAY)=0  \ / ------------- surface ----------- AP(NLAY+1), BP(NLAY+1)  
   ! //////////////////////////////////  
   
   ! U(IMR,JNP,NLAY) & V(IMR,JNP,NLAY):winds (m/s) at mid-time-level (t+NDT/2)  
   ! U and V may need to be polar filtered in advance in some cases.  
   
   ! IGD:      grid type on which winds are defined.  
   ! IGD = 0:  A-Grid  [all variables defined at the same point from south  
   ! pole (j=1) to north pole (j=JNP) ]  
   
   ! IGD = 1  GEOS-GCM C-Grid  
   ! [North]  
   
   ! V(i,j)  
   ! |  
   ! |  
   ! |  
   ! U(i-1,j)---Q(i,j)---U(i,j) [EAST]  
   ! |  
   ! |  
   ! |  
   ! V(i,j-1)  
   
   ! U(i,  1) is defined at South Pole.  
   ! V(i,  1) is half grid north of the South Pole.  
   ! V(i,JMR) is half grid south of the North Pole.  
   
   ! V must be defined at j=1 and j=JMR if IGD=1  
   ! V at JNP need not be given.  
   
   ! NDT: time step in seconds (need not be constant during the course of  
   ! the integration). Suggested value: 30 min. for 4x5, 15 min. for 2x2.5  
   ! (Lat-Lon) resolution. Smaller values are recommanded if the model  
   ! has a well-resolved stratosphere.  
   
   ! J1 defines the size of the polar cap:  
   ! South polar cap edge is located at -90 + (j1-1.5)*180/(JNP-1) deg.  
   ! North polar cap edge is located at  90 - (j1-1.5)*180/(JNP-1) deg.  
   ! There are currently only two choices (j1=2 or 3).  
   ! IMR must be an even integer if j1 = 2. Recommended value: J1=3.  
   
   ! IORD, JORD, and KORD are integers controlling various options in E-W,  
   ! N-S,  
   ! and vertical transport, respectively. Recommended values for positive  
   ! definite scalars: IORD=JORD=3, KORD=5. Use KORD=3 for non-  
   ! positive definite scalars or when linear correlation between constituents  
   ! is to be maintained.  
   
   ! _ORD=  
   ! 1: 1st order upstream scheme (too diffusive, not a useful option; it  
   ! can be used for debugging purposes; this is THE only known "linear"  
   ! monotonic advection scheme.).  
   ! 2: 2nd order van Leer (full monotonicity constraint;  
   ! see Lin et al 1994, MWR)  
   ! 3: monotonic PPM* (slightly improved PPM of Collela & Woodward 1984)  
   ! 4: semi-monotonic PPM (same as 3, but overshoots are allowed)  
   ! 5: positive-definite PPM (constraint on the subgrid distribution is  
   ! only strong enough to prevent generation of negative values;  
   ! both overshoots & undershoots are possible).  
   ! 6: un-constrained PPM (nearly diffusion free; slightly faster but  
   ! positivity not quaranteed. Use this option only when the fields  
   ! and winds are very smooth).  
   
   ! *PPM: Piece-wise Parabolic Method  
   
   ! Note that KORD <=2 options are no longer supported. DO not use option 4  
   ! or 5.  
   ! for non-positive definite scalars (such as Ertel Potential Vorticity).  
   
   ! The implicit numerical diffusion decreases as _ORD increases.  
   ! The last two options (ORDER=5, 6) should only be used when there is  
   ! significant explicit diffusion (such as a turbulence parameterization).  
   ! You  
   ! might get dispersive results otherwise.  
   ! No filter of any kind is applied to the constituent fields here.  
   
   ! AE: Radius of the sphere (meters).  
   ! Recommended value for the planet earth: 6.371E6  
   
   ! fill(logical):   flag to do filling for negatives (see note below).  
   
   ! Umax: Estimate (upper limit) of the maximum U-wind speed (m/s).  
   ! (220 m/s is a good value for troposphere model; 280 m/s otherwise)  
   
   ! =============  
   ! Output  
   ! =============  
   
   ! Q: mixing ratios at future time (t+NDT) (original values are  
   ! over-written)  
   ! W(NLAY): large-scale vertical mass flux as diagnosed from the hydrostatic  
   ! relationship. W will have the same unit as PS1 and PS2 (eg, mb).  
   ! W must be divided by NDT to get the correct mass-flux unit.  
   ! The vertical Courant number C = W/delp_UPWIND, where delp_UPWIND  
   ! is the pressure thickness in the "upwind" direction. For example,  
   ! C(k) = W(k)/delp(k)   if W(k) > 0;  
   ! C(k) = W(k)/delp(k+1) if W(k) < 0.  
   ! ( W > 0 is downward, ie, toward surface)  
   ! PS2: predicted PS at t+NDT (original values are over-written)  
   
   ! ********************************************************************  
   ! NOTES:  
   ! This forward-in-time upstream-biased transport scheme reduces to  
   ! the 2nd order center-in-time center-in-space mass continuity eqn.  
   ! if Q = 1 (constant fields will remain constant). This also ensures  
   ! that the computed vertical velocity to be identical to GEOS-1 GCM  
   ! for on-line transport.  
   
   ! A larger polar cap is used if j1=3 (recommended for C-Grid winds or when  
   ! winds are noisy near poles).  
   
   ! Flux-Form Semi-Lagrangian transport in the East-West direction is used  
   ! when and where Courant number is greater than one.  
   
   ! The user needs to change the parameter Jmax or Kmax if the resolution  
   ! is greater than 0.5 deg in N-S or 150 layers in the vertical direction.  
   ! (this TransPort Core is otherwise resolution independent and can be used  
   ! as a library routine).  
   
   ! PPM is 4th order accurate when grid spacing is uniform (x & y); 3rd  
   ! order accurate for non-uniform grid (vertical sigma coord.).  
   
   ! Time step is limitted only by transport in the meridional direction.  
   ! (the FFSL scheme is not implemented in the meridional direction).  
   
   ! Since only 1-D limiters are applied, negative values could  
   ! potentially be generated when large time step is used and when the  
   ! initial fields contain discontinuities.  
   ! This does not necessarily imply the integration is unstable.  
   ! These negatives are typically very small. A filling algorithm is  
   ! activated if the user set "fill" to be true.  
   
   ! The van Leer scheme used here is nearly as accurate as the original PPM  
   ! due to the use of a 4th order accurate reference slope. The PPM imple-  
   ! mented here is an improvement over the original and is also based on  
   ! the 4th order reference slope.  
   
   ! ****6***0*********0*********0*********0*********0*********0**********72  
   
   ! User modifiable parameters  
   
   PARAMETER (jmax=361, kmax=150)  
   
   ! ****6***0*********0*********0*********0*********0*********0**********72  
   
   ! Input-Output arrays  
   
   
   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  
   
   ! Local dynamic arrays  
   
   REAL crx(imr, jnp), cry(imr, jnp), xmass(imr, jnp), ymass(imr, jnp), &  
     fx1(imr+1), dpi(imr, jnp, nlay), delp1(imr, jnp, nlay), &  
     wk1(imr, jnp, nlay), pu(imr, jnp), pv(imr, jnp), dc2(imr, jnp), &  
     delp2(imr, jnp, nlay), dq(imr, jnp, nlay, nc), va(imr, jnp), &  
     ua(imr, jnp), qtmp(-imr:2*imr)  
   
   ! Local static  arrays  
   
   REAL dtdx(jmax), dtdx5(jmax), acosp(jmax), cosp(jmax), cose(jmax), &  
     dap(kmax), dbk(kmax)  
   DATA ndt0, nstep/0, 0/  
   DATA cross/.TRUE./  
   SAVE dtdy, dtdy5, rcap, js0, jn0, iml, dtdx, dtdx5, acosp, cosp, cose, dap, &  
     dbk  
   
   
   jmr = jnp - 1  
   imjm = imr*jnp  
   j2 = jnp - j1 + 1  
   nstep = nstep + 1  
   
   ! *********** Initialization **********************  
   IF (nstep==1) THEN  
   
     WRITE (6, *) '------------------------------------ '  
     WRITE (6, *) 'NASA/GSFC Transport Core Version 4.5'  
     WRITE (6, *) '------------------------------------ '  
   
     WRITE (6, *) 'IMR=', imr, ' JNP=', jnp, ' NLAY=', nlay, ' j1=', j1  
     WRITE (6, *) 'NC=', nc, iord, jord, kord, ndt  
   
     ! controles sur les parametres  
     IF (nlay<6) THEN  
       WRITE (6, *) 'NLAY must be >= 6'  
       STOP  
     END IF  
     IF (jnp<nlay) THEN  
       WRITE (6, *) 'JNP must be >= NLAY'  
       STOP  
     END IF  
     imrd2 = mod(imr, 2)  
     IF (j1==2 .AND. imrd2/=0) THEN  
       WRITE (6, *) 'if j1=2 IMR must be an even integer'  
       STOP  
     END IF  
   
   
     IF (jmax<jnp .OR. kmax<nlay) THEN  
       WRITE (6, *) 'Jmax or Kmax is too small'  
       STOP  
     END IF  
   
     DO k = 1, nlay  
       dap(k) = (ap(k+1)-ap(k))*pt  
       dbk(k) = bp(k+1) - bp(k)  
     END DO  
6    
7      pi = 4.*atan(1.)    SUBROUTINE ppm3d(igd, q, ps1, ps2, u, v, w, ndt, iord, jord, kord, nc, imr, &
8      dl = 2.*pi/float(imr)         jnp, j1, nlay, ap, bp, pt, ae, fill, umax)
     dp = pi/float(jmr)  
   
     IF (igd==0) THEN  
       ! Compute analytic cosine at cell edges  
       CALL cosa(cosp, cose, jnp, pi, dp)  
     ELSE  
       ! Define cosine consistent with GEOS-GCM (using dycore2.0 or later)  
       CALL cosc(cosp, cose, jnp, pi, dp)  
     END IF  
9    
10      DO j = 2, jmr      ! INPUT:
11        acosp(j) = 1./cosp(j)      ! =============
12      END DO  
13        ! Q(IMR,JNP,NLAY,NC): mixing ratios at current time (t)
14        ! NC: total number of constituents
15        ! IMR: first dimension (E-W); number of Grid intervals in E-W is IMR
16        ! JNP: 2nd dimension (N-S); number of Grid intervals in N-S is JNP-1
17        ! NLAY: 3rd dimension (number of layers); vertical index increases from 1
18        ! at
19        ! the model top to NLAY near the surface (see fig. below).
20        ! It is assumed that 6 <= NLAY <= JNP (for dynamic memory allocation)
21    
22        ! PS1(IMR,JNP): surface pressure at current time (t)
23        ! PS2(IMR,JNP): surface pressure at mid-time-level (t+NDT/2)
24        ! PS2 is replaced by the predicted PS (at t+NDT) on output.
25        ! Note: surface pressure can have any unit or can be multiplied by any
26        ! const.
27    
28        ! The pressure at layer edges are defined as follows:
29    
30        ! p(i,j,k) = AP(k)*PT  +  BP(k)*PS(i,j)          (1)
31    
32        ! Where PT is a constant having the same unit as PS.
33        ! AP and BP are unitless constants given at layer edges
34        ! defining the vertical coordinate.
35        ! BP(1) = 0., BP(NLAY+1) = 1.
36        ! The pressure at the model top is PTOP = AP(1)*PT
37    
38        ! For pure sigma system set AP(k) = 1 for all k, PT = PTOP,
39        ! BP(k) = sige(k) (sigma at edges), PS = Psfc - PTOP.
40    
41        ! Note: the sigma-P coordinate is a subset of Eq. 1, which in turn
42        ! is a subset of the following even more general sigma-P-thelta coord.
43        ! currently under development.
44        ! p(i,j,k) = (AP(k)*PT + BP(k)*PS(i,j))/(D(k)-C(k)*TE**(-1/kapa))
45    
46        ! Cf. ppm3d.txt.
47    
48        ! U(IMR,JNP,NLAY) & V(IMR,JNP,NLAY):winds (m/s) at mid-time-level (t+NDT/2)
49        ! U and V may need to be polar filtered in advance in some cases.
50    
51        ! IGD:      grid type on which winds are defined.
52        ! IGD = 0:  A-Grid  [all variables defined at the same point from south
53        ! pole (j=1) to north pole (j=JNP) ]
54    
55        ! IGD = 1  GEOS-GCM C-Grid
56        ! Cf. ppm3d.txt.
57    
58        ! U(i,  1) is defined at South Pole.
59        ! V(i,  1) is half grid north of the South Pole.
60        ! V(i,JMR) is half grid south of the North Pole.
61    
62        ! V must be defined at j=1 and j=JMR if IGD=1
63        ! V at JNP need not be given.
64    
65        ! NDT: time step in seconds (need not be constant during the course of
66        ! the integration). Suggested value: 30 min. for 4x5, 15 min. for 2x2.5
67        ! (Lat-Lon) resolution. Smaller values are recommanded if the model
68        ! has a well-resolved stratosphere.
69    
70        ! J1 defines the size of the polar cap:
71        ! South polar cap edge is located at -90 + (j1-1.5)*180/(JNP-1) deg.
72        ! North polar cap edge is located at  90 - (j1-1.5)*180/(JNP-1) deg.
73        ! There are currently only two choices (j1=2 or 3).
74        ! IMR must be an even integer if j1 = 2. Recommended value: J1=3.
75    
76        ! IORD, JORD, and KORD are integers controlling various options in E-W,
77        ! N-S,
78        ! and vertical transport, respectively. Recommended values for positive
79        ! definite scalars: IORD=JORD=3, KORD=5. Use KORD=3 for non-
80        ! positive definite scalars or when linear correlation between constituents
81        ! is to be maintained.
82    
83        ! _ORD=
84        ! 1: 1st order upstream scheme (too diffusive, not a useful option; it
85        ! can be used for debugging purposes; this is THE only known "linear"
86        ! monotonic advection scheme.).
87        ! 2: 2nd order van Leer (full monotonicity constraint;
88        ! see Lin et al 1994, MWR)
89        ! 3: monotonic PPM* (slightly improved PPM of Collela & Woodward 1984)
90        ! 4: semi-monotonic PPM (same as 3, but overshoots are allowed)
91        ! 5: positive-definite PPM (constraint on the subgrid distribution is
92        ! only strong enough to prevent generation of negative values;
93        ! both overshoots & undershoots are possible).
94        ! 6: un-constrained PPM (nearly diffusion free; slightly faster but
95        ! positivity not quaranteed. Use this option only when the fields
96        ! and winds are very smooth).
97    
98        ! *PPM: Piece-wise Parabolic Method
99    
100        ! Note that KORD <=2 options are no longer supported. DO not use option 4
101        ! or 5.
102        ! for non-positive definite scalars (such as Ertel Potential Vorticity).
103    
104        ! The implicit numerical diffusion decreases as _ORD increases.
105        ! The last two options (ORDER=5, 6) should only be used when there is
106        ! significant explicit diffusion (such as a turbulence parameterization).
107        ! You
108        ! might get dispersive results otherwise.
109        ! No filter of any kind is applied to the constituent fields here.
110    
111        ! AE: Radius of the sphere (meters).
112        ! Recommended value for the planet earth: 6.371E6
113    
114        ! fill(logical):   flag to do filling for negatives (see note below).
115    
116        ! Umax: Estimate (upper limit) of the maximum U-wind speed (m/s).
117        ! (220 m/s is a good value for troposphere model; 280 m/s otherwise)
118    
119        ! =============
120        ! Output
121        ! =============
122    
123        ! Q: mixing ratios at future time (t+NDT) (original values are
124        ! over-written)
125        ! W(NLAY): large-scale vertical mass flux as diagnosed from the hydrostatic
126        ! relationship. W will have the same unit as PS1 and PS2 (eg, mb).
127        ! W must be divided by NDT to get the correct mass-flux unit.
128        ! The vertical Courant number C = W/delp_UPWIND, where delp_UPWIND
129        ! is the pressure thickness in the "upwind" direction. For example,
130        ! C(k) = W(k)/delp(k)   if W(k) > 0;
131        ! C(k) = W(k)/delp(k+1) if W(k) < 0.
132        ! ( W > 0 is downward, ie, toward surface)
133        ! PS2: predicted PS at t+NDT (original values are over-written)
134    
135        ! ********************************************************************
136        ! NOTES:
137        ! This forward-in-time upstream-biased transport scheme reduces to
138        ! the 2nd order center-in-time center-in-space mass continuity eqn.
139        ! if Q = 1 (constant fields will remain constant). This also ensures
140        ! that the computed vertical velocity to be identical to GEOS-1 GCM
141        ! for on-line transport.
142    
143        ! A larger polar cap is used if j1=3 (recommended for C-Grid winds or when
144        ! winds are noisy near poles).
145    
146        ! Flux-Form Semi-Lagrangian transport in the East-West direction is used
147        ! when and where Courant number is greater than one.
148    
149        ! The user needs to change the parameter Jmax or Kmax if the resolution
150        ! is greater than 0.5 deg in N-S or 150 layers in the vertical direction.
151        ! (this TransPort Core is otherwise resolution independent and can be used
152        ! as a library routine).
153    
154        ! PPM is 4th order accurate when grid spacing is uniform (x & y); 3rd
155        ! order accurate for non-uniform grid (vertical sigma coord.).
156    
157        ! Time step is limitted only by transport in the meridional direction.
158        ! (the FFSL scheme is not implemented in the meridional direction).
159    
160        ! Since only 1-D limiters are applied, negative values could
161        ! potentially be generated when large time step is used and when the
162        ! initial fields contain discontinuities.
163        ! This does not necessarily imply the integration is unstable.
164        ! These negatives are typically very small. A filling algorithm is
165        ! activated if the user set "fill" to be true.
166    
167        ! The van Leer scheme used here is nearly as accurate as the original PPM
168        ! due to the use of a 4th order accurate reference slope. The PPM imple-
169        ! mented here is an improvement over the original and is also based on
170        ! the 4th order reference slope.
171    
172      ! Inverse of the Scaled polar cap area.      ! ****6***0*********0*********0*********0*********0*********0**********72
   
     rcap = dp/(imr*(1.-cos((j1-1.5)*dp)))  
     acosp(1) = rcap  
     acosp(jnp) = rcap  
   END IF  
   
   IF (ndt0/=ndt) THEN  
     dt = ndt  
     ndt0 = ndt  
   
     IF (umax<180.) THEN  
       WRITE (6, *) 'Umax may be too small!'  
     END IF  
     cr1 = abs(umax*dt)/(dl*ae)  
     maxdt = dp*ae/abs(umax) + 0.5  
     WRITE (6, *) 'Largest time step for max(V)=', umax, ' is ', maxdt  
     IF (maxdt<abs(ndt)) THEN  
       WRITE (6, *) 'Warning!!! NDT maybe too large!'  
     END IF  
   
     IF (cr1>=0.95) THEN  
       js0 = 0  
       jn0 = 0  
       iml = imr - 2  
       ztc = 0.  
     ELSE  
       ztc = acos(cr1)*(180./pi)  
   
       js0 = float(jmr)*(90.-ztc)/180. + 2  
       js0 = max(js0, j1+1)  
       iml = min(6*js0/(j1-1)+2, 4*imr/5)  
       jn0 = jnp - js0 + 1  
     END IF  
173    
174        ! User modifiable parameters
175    
176      DO j = 2, jmr      integer, PARAMETER:: jmax=361, kmax=150
       dtdx(j) = dt/(dl*ae*cosp(j))  
177    
178        dtdx5(j) = 0.5*dtdx(j)      ! ****6***0*********0*********0*********0*********0*********0**********72
     END DO  
179    
180        ! Input-Output arrays
181    
182      dtdy = dt/(ae*dp)      integer imr
183      dtdy5 = 0.5*dtdy      INTEGER igd, iord, jord, kord, nc, jnp, j1, nlay, ae
184        REAL q(imr, jnp, nlay, nc), ps1(imr, jnp), ps2(imr, jnp), &
185             u(imr, jnp, nlay), v(imr, jnp, nlay), ap(nlay+1), bp(nlay+1), &
186             w(imr, jnp, nlay), ndt, umax
187        INTEGER imrd2
188        REAL pt
189        LOGICAL cross, fill
190    
191        ! Local dynamic arrays
192    
193        REAL crx(imr, jnp), cry(imr, jnp), xmass(imr, jnp), ymass(imr, jnp), &
194             fx1(imr+1), dpi(imr, jnp, nlay), delp1(imr, jnp, nlay), &
195             wk1(imr, jnp, nlay), pu(imr, jnp), pv(imr, jnp), dc2(imr, jnp), &
196             delp2(imr, jnp, nlay), dq(imr, jnp, nlay, nc), va(imr, jnp), &
197             ua(imr, jnp), qtmp(-imr:2*imr)
198    
199        ! Local static  arrays
200    
201        REAL dtdx(jmax), dtdx5(jmax), acosp(jmax), cosp(jmax), cose(jmax), &
202             dap(kmax), dbk(kmax)
203        integer ndt0, nstep
204        DATA ndt0, nstep/0, 0/
205        DATA cross/.TRUE./
206        SAVE dtdy, dtdy5, rcap, js0, jn0, iml, dtdx, dtdx5, acosp, cosp, cose, dap, dbk
207        real cr1, d5, dl, dp, dt, dtdy, dtdy5, pi, rcap, ru, sum1, sum2, ztc
208        integer i, iad, ic, iiu, imh, imjm, iml, iu, j, j2, jad, jmp, jmr, jn, jn0
209        integer js, js0, jt, k, krd, l, maxdt
210    
211        jmr = jnp - 1
212        imjm = imr*jnp
213        j2 = jnp - j1 + 1
214        nstep = nstep + 1
215    
216        ! *********** Initialization **********************
217        IF (nstep==1) THEN
218    
219           WRITE (6, *) '------------------------------------ '
220           WRITE (6, *) 'NASA/GSFC Transport Core Version 4.5'
221           WRITE (6, *) '------------------------------------ '
222    
223           WRITE (6, *) 'IMR=', imr, ' JNP=', jnp, ' NLAY=', nlay, ' j1=', j1
224           WRITE (6, *) 'NC=', nc, iord, jord, kord, ndt
225    
226           ! controles sur les parametres
227           IF (nlay<6) THEN
228              WRITE (6, *) 'NLAY must be >= 6'
229              STOP
230           END IF
231           IF (jnp<nlay) THEN
232              WRITE (6, *) 'JNP must be >= NLAY'
233              STOP
234           END IF
235           imrd2 = mod(imr, 2)
236           IF (j1==2 .AND. imrd2/=0) THEN
237              WRITE (6, *) 'if j1=2 IMR must be an even integer'
238              STOP
239           END IF
240    
241    
242           IF (jmax<jnp .OR. kmax<nlay) THEN
243              WRITE (6, *) 'Jmax or Kmax is too small'
244              STOP
245           END IF
246    
247           DO k = 1, nlay
248              dap(k) = (ap(k+1)-ap(k))*pt
249              dbk(k) = bp(k+1) - bp(k)
250           END DO
251    
252           pi = 4.*atan(1.)
253           dl = 2.*pi/float(imr)
254           dp = pi/float(jmr)
255    
256           IF (igd==0) THEN
257              ! Compute analytic cosine at cell edges
258              CALL cosa(cosp, cose, jnp, pi, dp)
259           ELSE
260              ! Define cosine consistent with GEOS-GCM (using dycore2.0 or later)
261              CALL cosc(cosp, cose, jnp, pi, dp)
262           END IF
263    
264           DO j = 2, jmr
265              acosp(j) = 1./cosp(j)
266           END DO
267    
268           ! Inverse of the Scaled polar cap area.
269    
270           rcap = dp/(imr*(1.-cos((j1-1.5)*dp)))
271           acosp(1) = rcap
272           acosp(jnp) = rcap
273        END IF
274    
275        IF (ndt0/=ndt) THEN
276           dt = ndt
277           ndt0 = ndt
278    
279           IF (umax<180.) THEN
280              WRITE (6, *) 'Umax may be too small!'
281           END IF
282           cr1 = abs(umax*dt)/(dl*ae)
283           maxdt = dp*ae/abs(umax) + 0.5
284           WRITE (6, *) 'Largest time step for max(V)=', umax, ' is ', maxdt
285           IF (maxdt<abs(ndt)) THEN
286              WRITE (6, *) 'Warning!!! NDT maybe too large!'
287           END IF
288    
289           IF (cr1>=0.95) THEN
290              js0 = 0
291              jn0 = 0
292              iml = imr - 2
293              ztc = 0.
294           ELSE
295              ztc = acos(cr1)*(180./pi)
296    
297              js0 = float(jmr)*(90.-ztc)/180. + 2
298              js0 = max(js0, j1+1)
299              iml = min(6*js0/(j1-1)+2, 4*imr/5)
300              jn0 = jnp - js0 + 1
301           END IF
302    
303    
304           DO j = 2, jmr
305              dtdx(j) = dt/(dl*ae*cosp(j))
306    
307              dtdx5(j) = 0.5*dtdx(j)
308           END DO
309    
   END IF  
310    
311    ! *********** End Initialization **********************         dtdy = dt/(ae*dp)
312           dtdy5 = 0.5*dtdy
   ! delp = pressure thickness: the psudo-density in a hydrostatic system.  
   DO k = 1, nlay  
     DO j = 1, jnp  
       DO i = 1, imr  
         delp1(i, j, k) = dap(k) + dbk(k)*ps1(i, j)  
         delp2(i, j, k) = dap(k) + dbk(k)*ps2(i, j)  
       END DO  
     END DO  
   END DO  
313    
314        END IF
315    
316    IF (j1/=2) THEN      ! *********** End Initialization **********************
     DO ic = 1, nc  
       DO l = 1, nlay  
         DO i = 1, imr  
           q(i, 2, l, ic) = q(i, 1, l, ic)  
           q(i, jmr, l, ic) = q(i, jnp, l, ic)  
         END DO  
       END DO  
     END DO  
   END IF  
317    
318    ! Compute "tracer density"      ! delp = pressure thickness: the psudo-density in a hydrostatic system.
   DO ic = 1, nc  
319      DO k = 1, nlay      DO k = 1, nlay
320        DO j = 1, jnp         DO j = 1, jnp
321          DO i = 1, imr            DO i = 1, imr
322            dq(i, j, k, ic) = q(i, j, k, ic)*delp1(i, j, k)               delp1(i, j, k) = dap(k) + dbk(k)*ps1(i, j)
323          END DO               delp2(i, j, k) = dap(k) + dbk(k)*ps2(i, j)
324        END DO            END DO
325           END DO
326      END DO      END DO
   END DO  
327    
   DO k = 1, nlay  
328    
329      IF (igd==0) THEN      IF (j1/=2) THEN
330        ! Convert winds on A-Grid to Courant number on C-Grid.         DO ic = 1, nc
331        CALL a2c(u(1,1,k), v(1,1,k), imr, jmr, j1, j2, crx, cry, dtdx5, dtdy5)            DO l = 1, nlay
332      ELSE               DO i = 1, imr
333        ! Convert winds on C-grid to Courant number                  q(i, 2, l, ic) = q(i, 1, l, ic)
334        DO j = j1, j2                  q(i, jmr, l, ic) = q(i, jnp, l, ic)
335          DO i = 2, imr               END DO
336            crx(i, j) = dtdx(j)*u(i-1, j, k)            END DO
337          END DO         END DO
       END DO  
   
   
       DO j = j1, j2  
         crx(1, j) = dtdx(j)*u(imr, j, k)  
       END DO  
   
       DO i = 1, imr*jmr  
         cry(i, 2) = dtdy*v(i, 1, k)  
       END DO  
338      END IF      END IF
339    
340      ! Determine JS and JN      ! Compute "tracer density"
341      js = j1      DO ic = 1, nc
342      jn = j2         DO k = 1, nlay
343              DO j = 1, jnp
344      DO j = js0, j1 + 1, -1               DO i = 1, imr
345        DO i = 1, imr                  dq(i, j, k, ic) = q(i, j, k, ic)*delp1(i, j, k)
346          IF (abs(crx(i,j))>1.) THEN               END DO
347            js = j            END DO
348            GO TO 2222         END DO
         END IF  
       END DO  
     END DO  
   
 2222 CONTINUE  
     DO j = jn0, j2 - 1  
       DO i = 1, imr  
         IF (abs(crx(i,j))>1.) THEN  
           jn = j  
           GO TO 2233  
         END IF  
       END DO  
349      END DO      END DO
 2233 CONTINUE  
350    
351      IF (j1/=2) THEN ! Enlarged polar cap.      DO k = 1, nlay
       DO i = 1, imr  
         dpi(i, 2, k) = 0.  
         dpi(i, jmr, k) = 0.  
       END DO  
     END IF  
   
     ! ******* Compute horizontal mass fluxes ************  
   
     ! N-S component  
     DO j = j1, j2 + 1  
       d5 = 0.5*cose(j)  
       DO i = 1, imr  
         ymass(i, j) = cry(i, j)*d5*(delp2(i,j,k)+delp2(i,j-1,k))  
       END DO  
     END DO  
352    
353      DO j = j1, j2         IF (igd==0) THEN
354        DO i = 1, imr            ! Convert winds on A-Grid to Courant number on C-Grid.
355          dpi(i, j, k) = (ymass(i,j)-ymass(i,j+1))*acosp(j)            CALL a2c(u(1,1,k), v(1,1,k), imr, jmr, j1, j2, crx, cry, dtdx5, dtdy5)
356        END DO         ELSE
357      END DO            ! Convert winds on C-grid to Courant number
358              DO j = j1, j2
359                 DO i = 2, imr
360                    crx(i, j) = dtdx(j)*u(i-1, j, k)
361                 END DO
362              END DO
363    
     ! Poles  
     sum1 = ymass(imr, j1)  
     sum2 = ymass(imr, j2+1)  
     DO i = 1, imr - 1  
       sum1 = sum1 + ymass(i, j1)  
       sum2 = sum2 + ymass(i, j2+1)  
     END DO  
364    
365      sum1 = -sum1*rcap            DO j = j1, j2
366      sum2 = sum2*rcap               crx(1, j) = dtdx(j)*u(imr, j, k)
367      DO i = 1, imr            END DO
       dpi(i, 1, k) = sum1  
       dpi(i, jnp, k) = sum2  
     END DO  
368    
369      ! E-W component            DO i = 1, imr*jmr
370                 cry(i, 2) = dtdy*v(i, 1, k)
371              END DO
372           END IF
373    
374      DO j = j1, j2         ! Determine JS and JN
375        DO i = 2, imr         js = j1
376          pu(i, j) = 0.5*(delp2(i,j,k)+delp2(i-1,j,k))         jn = j2
       END DO  
     END DO  
377    
378      DO j = j1, j2         DO j = js0, j1 + 1, -1
379        pu(1, j) = 0.5*(delp2(1,j,k)+delp2(imr,j,k))            DO i = 1, imr
380      END DO               IF (abs(crx(i,j))>1.) THEN
381                    js = j
382                    GO TO 2222
383                 END IF
384              END DO
385           END DO
386    
387      DO j = j1, j2  2222   CONTINUE
388        DO i = 1, imr         DO j = jn0, j2 - 1
389          xmass(i, j) = pu(i, j)*crx(i, j)            DO i = 1, imr
390        END DO               IF (abs(crx(i,j))>1.) THEN
391      END DO                  jn = j
392                    GO TO 2233
393                 END IF
394              END DO
395           END DO
396    2233   CONTINUE
397    
398      DO j = j1, j2         IF (j1/=2) THEN ! Enlarged polar cap.
399        DO i = 1, imr - 1            DO i = 1, imr
400          dpi(i, j, k) = dpi(i, j, k) + xmass(i, j) - xmass(i+1, j)               dpi(i, 2, k) = 0.
401        END DO               dpi(i, jmr, k) = 0.
402      END DO            END DO
403           END IF
404    
405      DO j = j1, j2         ! ******* Compute horizontal mass fluxes ************
       dpi(imr, j, k) = dpi(imr, j, k) + xmass(imr, j) - xmass(1, j)  
     END DO  
406    
407      DO j = j1, j2         ! N-S component
408        DO i = 1, imr - 1         DO j = j1, j2 + 1
409          ua(i, j) = 0.5*(crx(i,j)+crx(i+1,j))            d5 = 0.5*cose(j)
410        END DO            DO i = 1, imr
411      END DO               ymass(i, j) = cry(i, j)*d5*(delp2(i,j,k)+delp2(i,j-1,k))
412              END DO
413           END DO
414    
415      DO j = j1, j2         DO j = j1, j2
416        ua(imr, j) = 0.5*(crx(imr,j)+crx(1,j))            DO i = 1, imr
417      END DO               dpi(i, j, k) = (ymass(i,j)-ymass(i,j+1))*acosp(j)
418      ! cccccccccccccccccccccccccccccccccccccccccccccccccccccc            END DO
419      ! Rajouts pour LMDZ.3.3         END DO
     ! cccccccccccccccccccccccccccccccccccccccccccccccccccccc  
     DO i = 1, imr  
       DO j = 1, jnp  
         va(i, j) = 0.  
       END DO  
     END DO  
420    
421      DO i = 1, imr*(jmr-1)         ! Poles
422        va(i, 2) = 0.5*(cry(i,2)+cry(i,3))         sum1 = ymass(imr, j1)
423      END DO         sum2 = ymass(imr, j2+1)
424           DO i = 1, imr - 1
425              sum1 = sum1 + ymass(i, j1)
426              sum2 = sum2 + ymass(i, j2+1)
427           END DO
428    
429           sum1 = -sum1*rcap
430           sum2 = sum2*rcap
431           DO i = 1, imr
432              dpi(i, 1, k) = sum1
433              dpi(i, jnp, k) = sum2
434           END DO
435    
436           ! E-W component
437    
438           DO j = j1, j2
439              DO i = 2, imr
440                 pu(i, j) = 0.5*(delp2(i,j,k)+delp2(i-1,j,k))
441              END DO
442           END DO
443    
444      IF (j1==2) THEN         DO j = j1, j2
445        imh = imr/2            pu(1, j) = 0.5*(delp2(1,j,k)+delp2(imr,j,k))
446        DO i = 1, imh         END DO
         va(i, 1) = 0.5*(cry(i,2)-cry(i+imh,2))  
         va(i+imh, 1) = -va(i, 1)  
         va(i, jnp) = 0.5*(cry(i,jnp)-cry(i+imh,jmr))  
         va(i+imh, jnp) = -va(i, jnp)  
       END DO  
       va(imr, 1) = va(1, 1)  
       va(imr, jnp) = va(1, jnp)  
     END IF  
447    
448      ! ****6***0*********0*********0*********0*********0*********0**********72         DO j = j1, j2
449      DO ic = 1, nc            DO i = 1, imr
450                 xmass(i, j) = pu(i, j)*crx(i, j)
451              END DO
452           END DO
453    
454        DO i = 1, imjm         DO j = j1, j2
455          wk1(i, 1, 1) = 0.            DO i = 1, imr - 1
456          wk1(i, 1, 2) = 0.               dpi(i, j, k) = dpi(i, j, k) + xmass(i, j) - xmass(i+1, j)
457        END DO            END DO
458           END DO
       ! E-W advective cross term  
       DO j = j1, j2  
         IF (j>js .AND. j<jn) GO TO 250  
   
         DO i = 1, imr  
           qtmp(i) = q(i, j, k, ic)  
         END DO  
   
         DO i = -iml, 0  
           qtmp(i) = q(imr+i, j, k, ic)  
           qtmp(imr+1-i) = q(1-i, j, k, ic)  
         END DO  
   
         DO i = 1, imr  
           iu = ua(i, j)  
           ru = ua(i, j) - iu  
           iiu = i - iu  
           IF (ua(i,j)>=0.) THEN  
             wk1(i, j, 1) = qtmp(iiu) + ru*(qtmp(iiu-1)-qtmp(iiu))  
           ELSE  
             wk1(i, j, 1) = qtmp(iiu) + ru*(qtmp(iiu)-qtmp(iiu+1))  
           END IF  
           wk1(i, j, 1) = wk1(i, j, 1) - qtmp(i)  
         END DO  
 250   END DO  
459    
460        IF (jn/=0) THEN         DO j = j1, j2
461          DO j = js + 1, jn - 1            dpi(imr, j, k) = dpi(imr, j, k) + xmass(imr, j) - xmass(1, j)
462           END DO
463    
464           DO j = j1, j2
465              DO i = 1, imr - 1
466                 ua(i, j) = 0.5*(crx(i,j)+crx(i+1,j))
467              END DO
468           END DO
469    
470            DO i = 1, imr         DO j = j1, j2
471              qtmp(i) = q(i, j, k, ic)            ua(imr, j) = 0.5*(crx(imr,j)+crx(1,j))
472           END DO
473           ! cccccccccccccccccccccccccccccccccccccccccccccccccccccc
474           ! Rajouts pour LMDZ.3.3
475           ! cccccccccccccccccccccccccccccccccccccccccccccccccccccc
476           DO i = 1, imr
477              DO j = 1, jnp
478                 va(i, j) = 0.
479            END DO            END DO
480           END DO
481    
482            qtmp(0) = q(imr, j, k, ic)         DO i = 1, imr*(jmr-1)
483            qtmp(imr+1) = q(1, j, k, ic)            va(i, 2) = 0.5*(cry(i,2)+cry(i,3))
484           END DO
485    
486           IF (j1==2) THEN
487              imh = imr/2
488              DO i = 1, imh
489                 va(i, 1) = 0.5*(cry(i,2)-cry(i+imh,2))
490                 va(i+imh, 1) = -va(i, 1)
491                 va(i, jnp) = 0.5*(cry(i,jnp)-cry(i+imh,jmr))
492                 va(i+imh, jnp) = -va(i, jnp)
493              END DO
494              va(imr, 1) = va(1, 1)
495              va(imr, jnp) = va(1, jnp)
496           END IF
497    
498           ! ****6***0*********0*********0*********0*********0*********0**********72
499           DO ic = 1, nc
500    
501              DO i = 1, imjm
502                 wk1(i, 1, 1) = 0.
503                 wk1(i, 1, 2) = 0.
504              END DO
505    
506            DO i = 1, imr            ! E-W advective cross term
507              iu = i - ua(i, j)            DO j = j1, j2
508              wk1(i, j, 1) = ua(i, j)*(qtmp(iu)-qtmp(iu+1))               IF (j>js .AND. j<jn) cycle
509    
510                 DO i = 1, imr
511                    qtmp(i) = q(i, j, k, ic)
512                 END DO
513    
514                 DO i = -iml, 0
515                    qtmp(i) = q(imr+i, j, k, ic)
516                    qtmp(imr+1-i) = q(1-i, j, k, ic)
517                 END DO
518    
519                 DO i = 1, imr
520                    iu = ua(i, j)
521                    ru = ua(i, j) - iu
522                    iiu = i - iu
523                    IF (ua(i,j)>=0.) THEN
524                       wk1(i, j, 1) = qtmp(iiu) + ru*(qtmp(iiu-1)-qtmp(iiu))
525                    ELSE
526                       wk1(i, j, 1) = qtmp(iiu) + ru*(qtmp(iiu)-qtmp(iiu+1))
527                    END IF
528                    wk1(i, j, 1) = wk1(i, j, 1) - qtmp(i)
529                 END DO
530            END DO            END DO
         END DO  
       END IF  
       ! ****6***0*********0*********0*********0*********0*********0**********72  
       ! Contribution from the N-S advection  
       DO i = 1, imr*(j2-j1+1)  
         jt = float(j1) - va(i, j1)  
         wk1(i, j1, 2) = va(i, j1)*(q(i,jt,k,ic)-q(i,jt+1,k,ic))  
       END DO  
531    
532        DO i = 1, imjm            IF (jn/=0) THEN
533          wk1(i, 1, 1) = q(i, 1, k, ic) + 0.5*wk1(i, 1, 1)               DO j = js + 1, jn - 1
         wk1(i, 1, 2) = q(i, 1, k, ic) + 0.5*wk1(i, 1, 2)  
       END DO  
534    
535        IF (cross) THEN                  DO i = 1, imr
536          ! Add cross terms in the vertical direction.                     qtmp(i) = q(i, j, k, ic)
537          IF (iord>=2) THEN                  END DO
538            iad = 2  
539          ELSE                  qtmp(0) = q(imr, j, k, ic)
540            iad = 1                  qtmp(imr+1) = q(1, j, k, ic)
541          END IF  
542                    DO i = 1, imr
543                       iu = i - ua(i, j)
544                       wk1(i, j, 1) = ua(i, j)*(qtmp(iu)-qtmp(iu+1))
545                    END DO
546                 END DO
547              END IF
548              ! ****6***0*********0*********0*********0*********0*********0**********72
549              ! Contribution from the N-S advection
550              DO i = 1, imr*(j2-j1+1)
551                 jt = float(j1) - va(i, j1)
552                 wk1(i, j1, 2) = va(i, j1)*(q(i,jt,k,ic)-q(i,jt+1,k,ic))
553              END DO
554    
555          IF (jord>=2) THEN            DO i = 1, imjm
556            jad = 2               wk1(i, 1, 1) = q(i, 1, k, ic) + 0.5*wk1(i, 1, 1)
557          ELSE               wk1(i, 1, 2) = q(i, 1, k, ic) + 0.5*wk1(i, 1, 2)
           jad = 1  
         END IF  
         CALL xadv(imr, jnp, j1, j2, wk1(1,1,2), ua, js, jn, iml, dc2, iad)  
         CALL yadv(imr, jnp, j1, j2, wk1(1,1,1), va, pv, w, jad)  
         DO j = 1, jnp  
           DO i = 1, imr  
             q(i, j, k, ic) = q(i, j, k, ic) + dc2(i, j) + pv(i, j)  
558            END DO            END DO
         END DO  
       END IF  
559    
560        CALL xtp(imr, jnp, iml, j1, j2, jn, js, pu, dq(1,1,k,ic), wk1(1,1,2), &            IF (cross) THEN
561          crx, fx1, xmass, iord)               ! Add cross terms in the vertical direction.
562                 IF (iord>=2) THEN
563                    iad = 2
564                 ELSE
565                    iad = 1
566                 END IF
567    
568                 IF (jord>=2) THEN
569                    jad = 2
570                 ELSE
571                    jad = 1
572                 END IF
573                 CALL xadv(imr, jnp, j1, j2, wk1(1,1,2), ua, js, jn, iml, dc2, iad)
574                 CALL yadv(imr, jnp, j1, j2, wk1(1,1,1), va, pv, w, jad)
575                 DO j = 1, jnp
576                    DO i = 1, imr
577                       q(i, j, k, ic) = q(i, j, k, ic) + dc2(i, j) + pv(i, j)
578                    END DO
579                 END DO
580              END IF
581    
582              CALL xtp(imr, jnp, iml, j1, j2, jn, js, pu, dq(1,1,k,ic), wk1(1,1,2), &
583                   crx, fx1, xmass, iord)
584    
585        CALL ytp(imr, jnp, j1, j2, acosp, rcap, dq(1,1,k,ic), wk1(1,1,1), cry, &            CALL ytp(imr, jnp, j1, j2, acosp, rcap, dq(1,1,k,ic), wk1(1,1,1), cry, &
586          dc2, ymass, wk1(1,1,3), wk1(1,1,4), wk1(1,1,5), wk1(1,1,6), jord)                 dc2, ymass, wk1(1,1,3), wk1(1,1,4), wk1(1,1,5), wk1(1,1,6), jord)
587    
588           END DO
589      END DO      END DO
   END DO  
590    
591    ! ******* Compute vertical mass flux (same unit as PS) ***********      ! ******* Compute vertical mass flux (same unit as PS) ***********
592    
593    ! 1st step: compute total column mass CONVERGENCE.      ! 1st step: compute total column mass CONVERGENCE.
594    
595    DO j = 1, jnp      DO j = 1, jnp
596      DO i = 1, imr         DO i = 1, imr
597        cry(i, j) = dpi(i, j, 1)            cry(i, j) = dpi(i, j, 1)
598           END DO
599      END DO      END DO
   END DO  
600    
601    DO k = 2, nlay      DO k = 2, nlay
602      DO j = 1, jnp         DO j = 1, jnp
603        DO i = 1, imr            DO i = 1, imr
604          cry(i, j) = cry(i, j) + dpi(i, j, k)               cry(i, j) = cry(i, j) + dpi(i, j, k)
605        END DO            END DO
606           END DO
607      END DO      END DO
   END DO  
608    
609    DO j = 1, jnp      DO j = 1, jnp
610      DO i = 1, imr         DO i = 1, imr
611    
612        ! 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption.            ! 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption.
613        ! Changes (increases) to surface pressure = total column mass            ! Changes (increases) to surface pressure = total column mass
614        ! convergence            ! convergence
615    
616        ps2(i, j) = ps1(i, j) + cry(i, j)            ps2(i, j) = ps1(i, j) + cry(i, j)
617    
618        ! 3rd step: compute vertical mass flux from mass conservation            ! 3rd step: compute vertical mass flux from mass conservation
619        ! principle.            ! principle.
620    
621        w(i, j, 1) = dpi(i, j, 1) - dbk(1)*cry(i, j)            w(i, j, 1) = dpi(i, j, 1) - dbk(1)*cry(i, j)
622        w(i, j, nlay) = 0.            w(i, j, nlay) = 0.
623           END DO
624      END DO      END DO
   END DO  
625    
626    DO k = 2, nlay - 1      DO k = 2, nlay - 1
627      DO j = 1, jnp         DO j = 1, jnp
628        DO i = 1, imr            DO i = 1, imr
629          w(i, j, k) = w(i, j, k-1) + dpi(i, j, k) - dbk(k)*cry(i, j)               w(i, j, k) = w(i, j, k-1) + dpi(i, j, k) - dbk(k)*cry(i, j)
630        END DO            END DO
631           END DO
632      END DO      END DO
   END DO  
633    
634    DO k = 1, nlay      DO k = 1, nlay
635      DO j = 1, jnp         DO j = 1, jnp
636        DO i = 1, imr            DO i = 1, imr
637          delp2(i, j, k) = dap(k) + dbk(k)*ps2(i, j)               delp2(i, j, k) = dap(k) + dbk(k)*ps2(i, j)
638        END DO            END DO
639           END DO
640      END DO      END DO
   END DO  
641    
642    krd = max(3, kord)      krd = max(3, kord)
643    DO ic = 1, nc      DO ic = 1, nc
644    
645      ! ****6***0*********0*********0*********0*********0*********0**********72         ! ****6***0*********0*********0*********0*********0*********0**********72
646    
647      CALL fzppm(imr, jnp, nlay, j1, dq(1,1,1,ic), w, q(1,1,1,ic), wk1, dpi, &         CALL fzppm(imr, jnp, nlay, j1, dq(1,1,1,ic), w, q(1,1,1,ic), wk1, dpi, &
648        dc2, crx, cry, pu, pv, xmass, ymass, delp1, krd)              dc2, crx, cry, pu, pv, xmass, ymass, delp1, krd)
649    
650    
651      IF (fill) CALL qckxyz(dq(1,1,1,ic), dc2, imr, jnp, nlay, j1, j2, cosp, &         IF (fill) CALL qckxyz(dq(1,1,1,ic), dc2, imr, jnp, nlay, j1, j2, cosp, &
652        acosp, .FALSE., ic, nstep)              acosp, .FALSE., ic, nstep)
653    
654      ! Recover tracer mixing ratio from "density" using predicted         ! Recover tracer mixing ratio from "density" using predicted
655      ! "air density" (pressure thickness) at time-level n+1         ! "air density" (pressure thickness) at time-level n+1
656    
657      DO k = 1, nlay         DO k = 1, nlay
658        DO j = 1, jnp            DO j = 1, jnp
659          DO i = 1, imr               DO i = 1, imr
660            q(i, j, k, ic) = dq(i, j, k, ic)/delp2(i, j, k)                  q(i, j, k, ic) = dq(i, j, k, ic)/delp2(i, j, k)
661          END DO               END DO
662        END DO            END DO
663           END DO
664    
665           IF (j1/=2) THEN
666              DO k = 1, nlay
667                 DO i = 1, imr
668                    ! j=1 c'est le p\^ole Sud, j=JNP c'est le p\^ole Nord
669                    q(i, 2, k, ic) = q(i, 1, k, ic)
670                    q(i, jmr, k, ic) = q(i, jmp, k, ic)
671                 END DO
672              END DO
673           END IF
674      END DO      END DO
675    
676      IF (j1/=2) THEN      IF (j1/=2) THEN
677        DO k = 1, nlay         DO k = 1, nlay
678          DO i = 1, imr            DO i = 1, imr
679            ! j=1 c'est le pôle Sud, j=JNP c'est le pôle Nord               w(i, 2, k) = w(i, 1, k)
680            q(i, 2, k, ic) = q(i, 1, k, ic)               w(i, jmr, k) = w(i, jnp, k)
681            q(i, jmr, k, ic) = q(i, jmp, k, ic)            END DO
682          END DO         END DO
       END DO  
683      END IF      END IF
   END DO  
   
   IF (j1/=2) THEN  
     DO k = 1, nlay  
       DO i = 1, imr  
         w(i, 2, k) = w(i, 1, k)  
         w(i, jmr, k) = w(i, jnp, k)  
       END DO  
     END DO  
   END IF  
684    
685    RETURN    END SUBROUTINE ppm3d
 END SUBROUTINE ppm3d  
686    
687  ! ****6***0*********0*********0*********0*********0*********0**********72  end module ppm3d_m
Legend:
Removed from v.166  
changed lines
  Added in v.169
  ViewVC Help
Powered by ViewVC 1.1.21