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