| 1 |
module ppm3d_m |
| 2 |
|
| 3 |
implicit none |
| 4 |
|
| 5 |
contains |
| 6 |
|
| 7 |
SUBROUTINE ppm3d(igd, q, ps1, ps2, u, v, w, ndt, iord, jord, kord, nc, imr, & |
| 8 |
jnp, j1, nlay, ap, bp, pt, ae, fill, umax) |
| 9 |
|
| 10 |
! INPUT: |
| 11 |
! ============= |
| 12 |
|
| 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 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
| 173 |
|
| 174 |
! User modifiable parameters |
| 175 |
|
| 176 |
integer, PARAMETER:: jmax=361, kmax=150 |
| 177 |
|
| 178 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
| 179 |
|
| 180 |
! Input-Output arrays |
| 181 |
|
| 182 |
integer imr |
| 183 |
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 |
|
| 310 |
|
| 311 |
dtdy = dt/(ae*dp) |
| 312 |
dtdy5 = 0.5*dtdy |
| 313 |
|
| 314 |
END IF |
| 315 |
|
| 316 |
! *********** End Initialization ********************** |
| 317 |
|
| 318 |
! delp = pressure thickness: the psudo-density in a hydrostatic system. |
| 319 |
DO k = 1, nlay |
| 320 |
DO j = 1, jnp |
| 321 |
DO i = 1, imr |
| 322 |
delp1(i, j, k) = dap(k) + dbk(k)*ps1(i, j) |
| 323 |
delp2(i, j, k) = dap(k) + dbk(k)*ps2(i, j) |
| 324 |
END DO |
| 325 |
END DO |
| 326 |
END DO |
| 327 |
|
| 328 |
|
| 329 |
IF (j1/=2) THEN |
| 330 |
DO ic = 1, nc |
| 331 |
DO l = 1, nlay |
| 332 |
DO i = 1, imr |
| 333 |
q(i, 2, l, ic) = q(i, 1, l, ic) |
| 334 |
q(i, jmr, l, ic) = q(i, jnp, l, ic) |
| 335 |
END DO |
| 336 |
END DO |
| 337 |
END DO |
| 338 |
END IF |
| 339 |
|
| 340 |
! Compute "tracer density" |
| 341 |
DO ic = 1, nc |
| 342 |
DO k = 1, nlay |
| 343 |
DO j = 1, jnp |
| 344 |
DO i = 1, imr |
| 345 |
dq(i, j, k, ic) = q(i, j, k, ic)*delp1(i, j, k) |
| 346 |
END DO |
| 347 |
END DO |
| 348 |
END DO |
| 349 |
END DO |
| 350 |
|
| 351 |
DO k = 1, nlay |
| 352 |
IF (igd==0) THEN |
| 353 |
! Convert winds on A-Grid to Courant number on C-Grid. |
| 354 |
CALL a2c(u(1,1,k), v(1,1,k), imr, jmr, j1, j2, crx, cry, dtdx5, dtdy5) |
| 355 |
ELSE |
| 356 |
! Convert winds on C-grid to Courant number |
| 357 |
DO j = j1, j2 |
| 358 |
DO i = 2, imr |
| 359 |
crx(i, j) = dtdx(j)*u(i-1, j, k) |
| 360 |
END DO |
| 361 |
END DO |
| 362 |
|
| 363 |
|
| 364 |
DO j = j1, j2 |
| 365 |
crx(1, j) = dtdx(j)*u(imr, j, k) |
| 366 |
END DO |
| 367 |
|
| 368 |
DO i = 1, imr*jmr |
| 369 |
cry(i, 2) = dtdy*v(i, 1, k) |
| 370 |
END DO |
| 371 |
END IF |
| 372 |
|
| 373 |
! Determine JS and JN |
| 374 |
js = j1 |
| 375 |
jn = j2 |
| 376 |
|
| 377 |
DO j = js0, j1 + 1, -1 |
| 378 |
DO i = 1, imr |
| 379 |
IF (abs(crx(i,j))>1.) THEN |
| 380 |
js = j |
| 381 |
GO TO 2222 |
| 382 |
END IF |
| 383 |
END DO |
| 384 |
END DO |
| 385 |
|
| 386 |
2222 CONTINUE |
| 387 |
DO j = jn0, j2 - 1 |
| 388 |
DO i = 1, imr |
| 389 |
IF (abs(crx(i,j))>1.) THEN |
| 390 |
jn = j |
| 391 |
GO TO 2233 |
| 392 |
END IF |
| 393 |
END DO |
| 394 |
END DO |
| 395 |
2233 CONTINUE |
| 396 |
|
| 397 |
IF (j1/=2) THEN ! Enlarged polar cap. |
| 398 |
DO i = 1, imr |
| 399 |
dpi(i, 2, k) = 0. |
| 400 |
dpi(i, jmr, k) = 0. |
| 401 |
END DO |
| 402 |
END IF |
| 403 |
|
| 404 |
! ******* Compute horizontal mass fluxes ************ |
| 405 |
|
| 406 |
! N-S component |
| 407 |
DO j = j1, j2 + 1 |
| 408 |
d5 = 0.5*cose(j) |
| 409 |
DO i = 1, imr |
| 410 |
ymass(i, j) = cry(i, j)*d5*(delp2(i,j,k)+delp2(i,j-1,k)) |
| 411 |
END DO |
| 412 |
END DO |
| 413 |
|
| 414 |
DO j = j1, j2 |
| 415 |
DO i = 1, imr |
| 416 |
dpi(i, j, k) = (ymass(i,j)-ymass(i,j+1))*acosp(j) |
| 417 |
END DO |
| 418 |
END DO |
| 419 |
|
| 420 |
! Poles |
| 421 |
sum1 = ymass(imr, j1) |
| 422 |
sum2 = ymass(imr, j2+1) |
| 423 |
DO i = 1, imr - 1 |
| 424 |
sum1 = sum1 + ymass(i, j1) |
| 425 |
sum2 = sum2 + ymass(i, j2+1) |
| 426 |
END DO |
| 427 |
|
| 428 |
sum1 = -sum1*rcap |
| 429 |
sum2 = sum2*rcap |
| 430 |
DO i = 1, imr |
| 431 |
dpi(i, 1, k) = sum1 |
| 432 |
dpi(i, jnp, k) = sum2 |
| 433 |
END DO |
| 434 |
|
| 435 |
! E-W component |
| 436 |
|
| 437 |
DO j = j1, j2 |
| 438 |
DO i = 2, imr |
| 439 |
pu(i, j) = 0.5*(delp2(i,j,k)+delp2(i-1,j,k)) |
| 440 |
END DO |
| 441 |
END DO |
| 442 |
|
| 443 |
DO j = j1, j2 |
| 444 |
pu(1, j) = 0.5*(delp2(1,j,k)+delp2(imr,j,k)) |
| 445 |
END DO |
| 446 |
|
| 447 |
DO j = j1, j2 |
| 448 |
DO i = 1, imr |
| 449 |
xmass(i, j) = pu(i, j)*crx(i, j) |
| 450 |
END DO |
| 451 |
END DO |
| 452 |
|
| 453 |
DO j = j1, j2 |
| 454 |
DO i = 1, imr - 1 |
| 455 |
dpi(i, j, k) = dpi(i, j, k) + xmass(i, j) - xmass(i+1, j) |
| 456 |
END DO |
| 457 |
END DO |
| 458 |
|
| 459 |
DO j = j1, j2 |
| 460 |
dpi(imr, j, k) = dpi(imr, j, k) + xmass(imr, j) - xmass(1, j) |
| 461 |
END DO |
| 462 |
|
| 463 |
DO j = j1, j2 |
| 464 |
DO i = 1, imr - 1 |
| 465 |
ua(i, j) = 0.5*(crx(i,j)+crx(i+1,j)) |
| 466 |
END DO |
| 467 |
END DO |
| 468 |
|
| 469 |
DO j = j1, j2 |
| 470 |
ua(imr, j) = 0.5*(crx(imr,j)+crx(1,j)) |
| 471 |
END DO |
| 472 |
! cccccccccccccccccccccccccccccccccccccccccccccccccccccc |
| 473 |
! Rajouts pour LMDZ.3.3 |
| 474 |
! cccccccccccccccccccccccccccccccccccccccccccccccccccccc |
| 475 |
DO i = 1, imr |
| 476 |
DO j = 1, jnp |
| 477 |
va(i, j) = 0. |
| 478 |
END DO |
| 479 |
END DO |
| 480 |
|
| 481 |
DO i = 1, imr*(jmr-1) |
| 482 |
va(i, 2) = 0.5*(cry(i,2)+cry(i,3)) |
| 483 |
END DO |
| 484 |
|
| 485 |
IF (j1==2) THEN |
| 486 |
imh = imr/2 |
| 487 |
DO i = 1, imh |
| 488 |
va(i, 1) = 0.5*(cry(i,2)-cry(i+imh,2)) |
| 489 |
va(i+imh, 1) = -va(i, 1) |
| 490 |
va(i, jnp) = 0.5*(cry(i,jnp)-cry(i+imh,jmr)) |
| 491 |
va(i+imh, jnp) = -va(i, jnp) |
| 492 |
END DO |
| 493 |
va(imr, 1) = va(1, 1) |
| 494 |
va(imr, jnp) = va(1, jnp) |
| 495 |
END IF |
| 496 |
|
| 497 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
| 498 |
DO ic = 1, nc |
| 499 |
|
| 500 |
DO i = 1, imjm |
| 501 |
wk1(i, 1, 1) = 0. |
| 502 |
wk1(i, 1, 2) = 0. |
| 503 |
END DO |
| 504 |
|
| 505 |
! E-W advective cross term |
| 506 |
DO j = j1, j2 |
| 507 |
IF (j>js .AND. j<jn) cycle |
| 508 |
|
| 509 |
DO i = 1, imr |
| 510 |
qtmp(i) = q(i, j, k, ic) |
| 511 |
END DO |
| 512 |
|
| 513 |
DO i = -iml, 0 |
| 514 |
qtmp(i) = q(imr+i, j, k, ic) |
| 515 |
qtmp(imr+1-i) = q(1-i, j, k, ic) |
| 516 |
END DO |
| 517 |
|
| 518 |
DO i = 1, imr |
| 519 |
iu = ua(i, j) |
| 520 |
ru = ua(i, j) - iu |
| 521 |
iiu = i - iu |
| 522 |
IF (ua(i,j)>=0.) THEN |
| 523 |
wk1(i, j, 1) = qtmp(iiu) + ru*(qtmp(iiu-1)-qtmp(iiu)) |
| 524 |
ELSE |
| 525 |
wk1(i, j, 1) = qtmp(iiu) + ru*(qtmp(iiu)-qtmp(iiu+1)) |
| 526 |
END IF |
| 527 |
wk1(i, j, 1) = wk1(i, j, 1) - qtmp(i) |
| 528 |
END DO |
| 529 |
END DO |
| 530 |
|
| 531 |
IF (jn/=0) THEN |
| 532 |
DO j = js + 1, jn - 1 |
| 533 |
|
| 534 |
DO i = 1, imr |
| 535 |
qtmp(i) = q(i, j, k, ic) |
| 536 |
END DO |
| 537 |
|
| 538 |
qtmp(0) = q(imr, j, k, ic) |
| 539 |
qtmp(imr+1) = q(1, j, k, ic) |
| 540 |
|
| 541 |
DO i = 1, imr |
| 542 |
iu = i - ua(i, j) |
| 543 |
wk1(i, j, 1) = ua(i, j)*(qtmp(iu)-qtmp(iu+1)) |
| 544 |
END DO |
| 545 |
END DO |
| 546 |
END IF |
| 547 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
| 548 |
! Contribution from the N-S advection |
| 549 |
DO i = 1, imr*(j2-j1+1) |
| 550 |
jt = float(j1) - va(i, j1) |
| 551 |
wk1(i, j1, 2) = va(i, j1)*(q(i,jt,k,ic)-q(i,jt+1,k,ic)) |
| 552 |
END DO |
| 553 |
|
| 554 |
DO i = 1, imjm |
| 555 |
wk1(i, 1, 1) = q(i, 1, k, ic) + 0.5*wk1(i, 1, 1) |
| 556 |
wk1(i, 1, 2) = q(i, 1, k, ic) + 0.5*wk1(i, 1, 2) |
| 557 |
END DO |
| 558 |
|
| 559 |
IF (cross) THEN |
| 560 |
! Add cross terms in the vertical direction. |
| 561 |
IF (iord>=2) THEN |
| 562 |
iad = 2 |
| 563 |
ELSE |
| 564 |
iad = 1 |
| 565 |
END IF |
| 566 |
|
| 567 |
IF (jord>=2) THEN |
| 568 |
jad = 2 |
| 569 |
ELSE |
| 570 |
jad = 1 |
| 571 |
END IF |
| 572 |
CALL xadv(imr, jnp, j1, j2, wk1(1,1,2), ua, js, jn, iml, dc2, iad) |
| 573 |
CALL yadv(imr, jnp, j1, j2, wk1(1,1,1), va, pv, w, jad) |
| 574 |
DO j = 1, jnp |
| 575 |
DO i = 1, imr |
| 576 |
q(i, j, k, ic) = q(i, j, k, ic) + dc2(i, j) + pv(i, j) |
| 577 |
END DO |
| 578 |
END DO |
| 579 |
END IF |
| 580 |
|
| 581 |
CALL xtp(imr, jnp, iml, j1, j2, jn, js, pu, dq(1,1,k,ic), wk1(1,1,2), & |
| 582 |
crx, fx1, xmass, iord) |
| 583 |
|
| 584 |
CALL ytp(imr, jnp, j1, j2, acosp, rcap, dq(1,1,k,ic), wk1(1,1,1), cry, & |
| 585 |
dc2, ymass, wk1(1,1,3), wk1(1,1,4), wk1(1,1,5), wk1(1,1,6), jord) |
| 586 |
|
| 587 |
END DO |
| 588 |
END DO |
| 589 |
|
| 590 |
! ******* Compute vertical mass flux (same unit as PS) *********** |
| 591 |
|
| 592 |
! 1st step: compute total column mass CONVERGENCE. |
| 593 |
|
| 594 |
DO j = 1, jnp |
| 595 |
DO i = 1, imr |
| 596 |
cry(i, j) = dpi(i, j, 1) |
| 597 |
END DO |
| 598 |
END DO |
| 599 |
|
| 600 |
DO k = 2, nlay |
| 601 |
DO j = 1, jnp |
| 602 |
DO i = 1, imr |
| 603 |
cry(i, j) = cry(i, j) + dpi(i, j, k) |
| 604 |
END DO |
| 605 |
END DO |
| 606 |
END DO |
| 607 |
|
| 608 |
DO j = 1, jnp |
| 609 |
DO i = 1, imr |
| 610 |
|
| 611 |
! 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption. |
| 612 |
! Changes (increases) to surface pressure = total column mass |
| 613 |
! convergence |
| 614 |
|
| 615 |
ps2(i, j) = ps1(i, j) + cry(i, j) |
| 616 |
|
| 617 |
! 3rd step: compute vertical mass flux from mass conservation |
| 618 |
! principle. |
| 619 |
|
| 620 |
w(i, j, 1) = dpi(i, j, 1) - dbk(1)*cry(i, j) |
| 621 |
w(i, j, nlay) = 0. |
| 622 |
END DO |
| 623 |
END DO |
| 624 |
|
| 625 |
DO k = 2, nlay - 1 |
| 626 |
DO j = 1, jnp |
| 627 |
DO i = 1, imr |
| 628 |
w(i, j, k) = w(i, j, k-1) + dpi(i, j, k) - dbk(k)*cry(i, j) |
| 629 |
END DO |
| 630 |
END DO |
| 631 |
END DO |
| 632 |
|
| 633 |
DO k = 1, nlay |
| 634 |
DO j = 1, jnp |
| 635 |
DO i = 1, imr |
| 636 |
delp2(i, j, k) = dap(k) + dbk(k)*ps2(i, j) |
| 637 |
END DO |
| 638 |
END DO |
| 639 |
END DO |
| 640 |
|
| 641 |
krd = max(3, kord) |
| 642 |
DO ic = 1, nc |
| 643 |
|
| 644 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
| 645 |
|
| 646 |
CALL fzppm(imr, jnp, nlay, j1, dq(1,1,1,ic), w, q(1,1,1,ic), wk1, dpi, & |
| 647 |
dc2, crx, cry, pu, pv, xmass, ymass, delp1, krd) |
| 648 |
|
| 649 |
|
| 650 |
IF (fill) CALL qckxyz(dq(1,1,1,ic), dc2, imr, jnp, nlay, j1, j2, cosp, & |
| 651 |
acosp, .FALSE., ic, nstep) |
| 652 |
|
| 653 |
! Recover tracer mixing ratio from "density" using predicted |
| 654 |
! "air density" (pressure thickness) at time-level n+1 |
| 655 |
|
| 656 |
DO k = 1, nlay |
| 657 |
DO j = 1, jnp |
| 658 |
DO i = 1, imr |
| 659 |
q(i, j, k, ic) = dq(i, j, k, ic)/delp2(i, j, k) |
| 660 |
END DO |
| 661 |
END DO |
| 662 |
END DO |
| 663 |
|
| 664 |
IF (j1/=2) THEN |
| 665 |
DO k = 1, nlay |
| 666 |
DO i = 1, imr |
| 667 |
! j=1 c'est le p\^ole Sud, j=JNP c'est le p\^ole Nord |
| 668 |
q(i, 2, k, ic) = q(i, 1, k, ic) |
| 669 |
q(i, jmr, k, ic) = q(i, jmp, k, ic) |
| 670 |
END DO |
| 671 |
END DO |
| 672 |
END IF |
| 673 |
END DO |
| 674 |
|
| 675 |
IF (j1/=2) THEN |
| 676 |
DO k = 1, nlay |
| 677 |
DO i = 1, imr |
| 678 |
w(i, 2, k) = w(i, 1, k) |
| 679 |
w(i, jmr, k) = w(i, jnp, k) |
| 680 |
END DO |
| 681 |
END DO |
| 682 |
END IF |
| 683 |
|
| 684 |
END SUBROUTINE ppm3d |
| 685 |
|
| 686 |
end module ppm3d_m |
Parent Directory
| 
Revision Log