| 1 | ! |
|---|
| 2 | ! $Id: ppm3d.F 5093 2024-07-21 11:07:18Z snguyen $ |
|---|
| 3 | ! |
|---|
| 4 | |
|---|
| 5 | cFrom lin@explorer.gsfc.nasa.gov Wed Apr 15 17:44:44 1998 |
|---|
| 6 | cDate: Wed, 15 Apr 1998 11:37:03 -0400 |
|---|
| 7 | cFrom: lin@explorer.gsfc.nasa.gov |
|---|
| 8 | cTo: Frederic.Hourdin@lmd.jussieu.fr |
|---|
| 9 | cSubject: 3D transport module of the GSFC CTM and GEOS GCM |
|---|
| 10 | |
|---|
| 11 | |
|---|
| 12 | cThis code is sent to you by S-J Lin, DAO, NASA-GSFC |
|---|
| 13 | |
|---|
| 14 | cNote: this version is intended for machines like CRAY |
|---|
| 15 | C-90. No multitasking directives implemented. |
|---|
| 16 | |
|---|
| 17 | |
|---|
| 18 | C ******************************************************************** |
|---|
| 19 | C |
|---|
| 20 | C TransPort Core for Goddard Chemistry Transport Model (G-CTM), Goddard |
|---|
| 21 | C Earth Observing System General Circulation Model (GEOS-GCM), and Data |
|---|
| 22 | C Assimilation System (GEOS-DAS). |
|---|
| 23 | C |
|---|
| 24 | C ******************************************************************** |
|---|
| 25 | C |
|---|
| 26 | C Purpose: given horizontal winds on a hybrid sigma-p surfaces, |
|---|
| 27 | C one call to tpcore updates the 3-D mixing ratio |
|---|
| 28 | C fields one time step (NDT). [vertical mass flux is computed |
|---|
| 29 | C internally consistent with the discretized hydrostatic mass |
|---|
| 30 | C continuity equation of the C-Grid GEOS-GCM (for IGD=1)]. |
|---|
| 31 | C |
|---|
| 32 | C Schemes: Multi-dimensional Flux Form Semi-Lagrangian (FFSL) scheme based |
|---|
| 33 | C on the van Leer or PPM. |
|---|
| 34 | C (see Lin and Rood 1996). |
|---|
| 35 | C Version 4.5 |
|---|
| 36 | C Last modified: Dec. 5, 1996 |
|---|
| 37 | C Major changes from version 4.0: a more general vertical hybrid sigma- |
|---|
| 38 | C pressure coordinate. |
|---|
| 39 | C Subroutines modified: xtp, ytp, fzppm, qckxyz |
|---|
| 40 | C Subroutines deleted: vanz |
|---|
| 41 | C |
|---|
| 42 | C Author: Shian-Jiann Lin |
|---|
| 43 | C mail address: |
|---|
| 44 | C Shian-Jiann Lin* |
|---|
| 45 | C Code 910.3, NASA/GSFC, Greenbelt, MD 20771 |
|---|
| 46 | C Phone: 301-286-9540 |
|---|
| 47 | C E-mail: lin@dao.gsfc.nasa.gov |
|---|
| 48 | C |
|---|
| 49 | C *affiliation: |
|---|
| 50 | C Joint Center for Earth Systems Technology |
|---|
| 51 | C The University of Maryland Baltimore County |
|---|
| 52 | C NASA - Goddard Space Flight Center |
|---|
| 53 | C References: |
|---|
| 54 | C |
|---|
| 55 | C 1. Lin, S.-J., and R. B. Rood, 1996: Multidimensional flux form semi- |
|---|
| 56 | C Lagrangian transport schemes. Mon. Wea. Rev., 124, 2046-2070. |
|---|
| 57 | C |
|---|
| 58 | C 2. Lin, S.-J., W. C. Chao, Y. C. Sud, and G. K. Walker, 1994: A class of |
|---|
| 59 | C the van Leer-type transport schemes and its applications to the moist- |
|---|
| 60 | C ure transport in a General Circulation Model. Mon. Wea. Rev., 122, |
|---|
| 61 | C 1575-1593. |
|---|
| 62 | C |
|---|
| 63 | C ****6***0*********0*********0*********0*********0*********0**********72 |
|---|
| 64 | C |
|---|
| 65 | subroutine ppm3d(IGD,Q,PS1,PS2,U,V,W,NDT,IORD,JORD,KORD,NC,IMR, |
|---|
| 66 | & JNP,j1,NLAY,AP,BP,PT,AE,fill,dum,Umax) |
|---|
| 67 | |
|---|
| 68 | implicit none |
|---|
| 69 | |
|---|
| 70 | c rajout de déclarations |
|---|
| 71 | c integer Jmax,kmax,ndt0,nstep,k,j,i,ic,l,js,jn,imh,iad,jad,krd |
|---|
| 72 | c integer iu,iiu,j2,jmr,js0,jt |
|---|
| 73 | c real dtdy,dtdy5,rcap,iml,jn0,imjm,pi,dl,dp |
|---|
| 74 | c real dt,cr1,maxdt,ztc,d5,sum1,sum2,ru |
|---|
| 75 | C |
|---|
| 76 | C ******************************************************************** |
|---|
| 77 | C |
|---|
| 78 | C ============= |
|---|
| 79 | C INPUT: |
|---|
| 80 | C ============= |
|---|
| 81 | C |
|---|
| 82 | C Q(IMR,JNP,NLAY,NC): mixing ratios at current time (t) |
|---|
| 83 | C NC: total # of constituents |
|---|
| 84 | C IMR: first dimension (E-W); # of Grid intervals in E-W is IMR |
|---|
| 85 | C JNP: 2nd dimension (N-S); # of Grid intervals in N-S is JNP-1 |
|---|
| 86 | C NLAY: 3rd dimension (# of layers); vertical index increases from 1 at |
|---|
| 87 | C the model top to NLAY near the surface (see fig. below). |
|---|
| 88 | C It is assumed that 6 <= NLAY <= JNP (for dynamic memory allocation) |
|---|
| 89 | C |
|---|
| 90 | C PS1(IMR,JNP): surface pressure at current time (t) |
|---|
| 91 | C PS2(IMR,JNP): surface pressure at mid-time-level (t+NDT/2) |
|---|
| 92 | C PS2 is replaced by the predicted PS (at t+NDT) on output. |
|---|
| 93 | C Note: surface pressure can have any unit or can be multiplied by any |
|---|
| 94 | C const. |
|---|
| 95 | C |
|---|
| 96 | C The pressure at layer edges are defined as follows: |
|---|
| 97 | C |
|---|
| 98 | C p(i,j,k) = AP(k)*PT + BP(k)*PS(i,j) (1) |
|---|
| 99 | C |
|---|
| 100 | C Where PT is a constant having the same unit as PS. |
|---|
| 101 | C AP and BP are unitless constants given at layer edges |
|---|
| 102 | C defining the vertical coordinate. |
|---|
| 103 | C BP(1) = 0., BP(NLAY+1) = 1. |
|---|
| 104 | C The pressure at the model top is PTOP = AP(1)*PT |
|---|
| 105 | C |
|---|
| 106 | C For pure sigma system set AP(k) = 1 for all k, PT = PTOP, |
|---|
| 107 | C BP(k) = sige(k) (sigma at edges), PS = Psfc - PTOP. |
|---|
| 108 | C |
|---|
| 109 | C Note: the sigma-P coordinate is a subset of Eq. 1, which in turn |
|---|
| 110 | C is a subset of the following even more general sigma-P-thelta coord. |
|---|
| 111 | C currently under development. |
|---|
| 112 | C p(i,j,k) = (AP(k)*PT + BP(k)*PS(i,j))/(D(k)-C(k)*TE**(-1/kapa)) |
|---|
| 113 | C |
|---|
| 114 | C ///////////////////////////////// |
|---|
| 115 | C / \ ------------- PTOP -------------- AP(1), BP(1) |
|---|
| 116 | C | |
|---|
| 117 | C delp(1) | ........... Q(i,j,1) ............ |
|---|
| 118 | C | |
|---|
| 119 | C W(1) \ / --------------------------------- AP(2), BP(2) |
|---|
| 120 | C |
|---|
| 121 | C |
|---|
| 122 | C |
|---|
| 123 | C W(k-1) / \ --------------------------------- AP(k), BP(k) |
|---|
| 124 | C | |
|---|
| 125 | C delp(K) | ........... Q(i,j,k) ............ |
|---|
| 126 | C | |
|---|
| 127 | C W(k) \ / --------------------------------- AP(k+1), BP(k+1) |
|---|
| 128 | C |
|---|
| 129 | C |
|---|
| 130 | C |
|---|
| 131 | C / \ --------------------------------- AP(NLAY), BP(NLAY) |
|---|
| 132 | C | |
|---|
| 133 | C delp(NLAY) | ........... Q(i,j,NLAY) ......... |
|---|
| 134 | C | |
|---|
| 135 | C W(NLAY)=0 \ / ------------- surface ----------- AP(NLAY+1), BP(NLAY+1) |
|---|
| 136 | C ////////////////////////////////// |
|---|
| 137 | C |
|---|
| 138 | C U(IMR,JNP,NLAY) & V(IMR,JNP,NLAY):winds (m/s) at mid-time-level (t+NDT/2) |
|---|
| 139 | C U and V may need to be polar filtered in advance in some cases. |
|---|
| 140 | C |
|---|
| 141 | C IGD: grid type on which winds are defined. |
|---|
| 142 | C IGD = 0: A-Grid [all variables defined at the same point from south |
|---|
| 143 | C pole (j=1) to north pole (j=JNP) ] |
|---|
| 144 | C |
|---|
| 145 | C IGD = 1 GEOS-GCM C-Grid |
|---|
| 146 | C [North] |
|---|
| 147 | C |
|---|
| 148 | C V(i,j) |
|---|
| 149 | C | |
|---|
| 150 | C | |
|---|
| 151 | C | |
|---|
| 152 | C U(i-1,j)---Q(i,j)---U(i,j) [EAST] |
|---|
| 153 | C | |
|---|
| 154 | C | |
|---|
| 155 | C | |
|---|
| 156 | C V(i,j-1) |
|---|
| 157 | C |
|---|
| 158 | C U(i, 1) is defined at South Pole. |
|---|
| 159 | C V(i, 1) is half grid north of the South Pole. |
|---|
| 160 | C V(i,JMR) is half grid south of the North Pole. |
|---|
| 161 | C |
|---|
| 162 | C V must be defined at j=1 and j=JMR if IGD=1 |
|---|
| 163 | C V at JNP need not be given. |
|---|
| 164 | C |
|---|
| 165 | C NDT: time step in seconds (need not be constant during the course of |
|---|
| 166 | C the integration). Suggested value: 30 min. for 4x5, 15 min. for 2x2.5 |
|---|
| 167 | C (Lat-Lon) resolution. Smaller values are recommanded if the model |
|---|
| 168 | C has a well-resolved stratosphere. |
|---|
| 169 | C |
|---|
| 170 | C J1 defines the size of the polar cap: |
|---|
| 171 | C South polar cap edge is located at -90 + (j1-1.5)*180/(JNP-1) deg. |
|---|
| 172 | C North polar cap edge is located at 90 - (j1-1.5)*180/(JNP-1) deg. |
|---|
| 173 | C There are currently only two choices (j1=2 or 3). |
|---|
| 174 | C IMR must be an even integer if j1 = 2. Recommended value: J1=3. |
|---|
| 175 | C |
|---|
| 176 | C IORD, JORD, and KORD are integers controlling various options in E-W, N-S, |
|---|
| 177 | C and vertical transport, respectively. Recommended values for positive |
|---|
| 178 | C definite scalars: IORD=JORD=3, KORD=5. Use KORD=3 for non- |
|---|
| 179 | C positive definite scalars or when linear correlation between constituents |
|---|
| 180 | C is to be maintained. |
|---|
| 181 | C |
|---|
| 182 | C _ORD= |
|---|
| 183 | C 1: 1st order upstream scheme (too diffusive, not a useful option; it |
|---|
| 184 | C can be used for debugging purposes; this is THE only known "linear" |
|---|
| 185 | C monotonic advection scheme.). |
|---|
| 186 | C 2: 2nd order van Leer (full monotonicity constraint; |
|---|
| 187 | C see Lin et al 1994, MWR) |
|---|
| 188 | C 3: monotonic PPM* (slightly improved PPM of Collela & Woodward 1984) |
|---|
| 189 | C 4: semi-monotonic PPM (same as 3, but overshoots are allowed) |
|---|
| 190 | C 5: positive-definite PPM (constraint on the subgrid distribution is |
|---|
| 191 | C only strong enough to prevent generation of negative values; |
|---|
| 192 | C both overshoots & undershoots are possible). |
|---|
| 193 | C 6: un-constrained PPM (nearly diffusion free; slightly faster but |
|---|
| 194 | C positivity not quaranteed. Use this option only when the fields |
|---|
| 195 | C and winds are very smooth). |
|---|
| 196 | C |
|---|
| 197 | C *PPM: Piece-wise Parabolic Method |
|---|
| 198 | C |
|---|
| 199 | C Note that KORD <=2 options are no longer supported. DO not use option 4 or 5. |
|---|
| 200 | C for non-positive definite scalars (such as Ertel Potential Vorticity). |
|---|
| 201 | C |
|---|
| 202 | C The implicit numerical diffusion decreases as _ORD increases. |
|---|
| 203 | C The last two options (ORDER=5, 6) should only be used when there is |
|---|
| 204 | C significant explicit diffusion (such as a turbulence parameterization). You |
|---|
| 205 | C might get dispersive results otherwise. |
|---|
| 206 | C No filter of any kind is applied to the constituent fields here. |
|---|
| 207 | C |
|---|
| 208 | C AE: Radius of the sphere (meters). |
|---|
| 209 | C Recommended value for the planet earth: 6.371E6 |
|---|
| 210 | C |
|---|
| 211 | C fill(logical): flag to do filling for negatives (see note below). |
|---|
| 212 | C |
|---|
| 213 | C Umax: Estimate (upper limit) of the maximum U-wind speed (m/s). |
|---|
| 214 | C (220 m/s is a good value for troposphere model; 280 m/s otherwise) |
|---|
| 215 | C |
|---|
| 216 | C ============= |
|---|
| 217 | C Output |
|---|
| 218 | C ============= |
|---|
| 219 | C |
|---|
| 220 | C Q: mixing ratios at future time (t+NDT) (original values are over-written) |
|---|
| 221 | C W(NLAY): large-scale vertical mass flux as diagnosed from the hydrostatic |
|---|
| 222 | C relationship. W will have the same unit as PS1 and PS2 (eg, mb). |
|---|
| 223 | C W must be divided by NDT to get the correct mass-flux unit. |
|---|
| 224 | C The vertical Courant number C = W/delp_UPWIND, where delp_UPWIND |
|---|
| 225 | C is the pressure thickness in the "upwind" direction. For example, |
|---|
| 226 | C C(k) = W(k)/delp(k) if W(k) > 0; |
|---|
| 227 | C C(k) = W(k)/delp(k+1) if W(k) < 0. |
|---|
| 228 | C ( W > 0 is downward, ie, toward surface) |
|---|
| 229 | C PS2: predicted PS at t+NDT (original values are over-written) |
|---|
| 230 | C |
|---|
| 231 | C ******************************************************************** |
|---|
| 232 | C NOTES: |
|---|
| 233 | C This forward-in-time upstream-biased transport scheme reduces to |
|---|
| 234 | C the 2nd order center-in-time center-in-space mass continuity eqn. |
|---|
| 235 | C if Q = 1 (constant fields will remain constant). This also ensures |
|---|
| 236 | C that the computed vertical velocity to be identical to GEOS-1 GCM |
|---|
| 237 | C for on-line transport. |
|---|
| 238 | C |
|---|
| 239 | C A larger polar cap is used if j1=3 (recommended for C-Grid winds or when |
|---|
| 240 | C winds are noisy near poles). |
|---|
| 241 | C |
|---|
| 242 | C Flux-Form Semi-Lagrangian transport in the East-West direction is used |
|---|
| 243 | C when and where Courant # is greater than one. |
|---|
| 244 | C |
|---|
| 245 | C The user needs to change the parameter Jmax or Kmax if the resolution |
|---|
| 246 | C is greater than 0.5 deg in N-S or 150 layers in the vertical direction. |
|---|
| 247 | C (this TransPort Core is otherwise resolution independent and can be used |
|---|
| 248 | C as a library routine). |
|---|
| 249 | C |
|---|
| 250 | C PPM is 4th order accurate when grid spacing is uniform (x & y); 3rd |
|---|
| 251 | C order accurate for non-uniform grid (vertical sigma coord.). |
|---|
| 252 | C |
|---|
| 253 | C Time step is limitted only by transport in the meridional direction. |
|---|
| 254 | C (the FFSL scheme is not implemented in the meridional direction). |
|---|
| 255 | C |
|---|
| 256 | C Since only 1-D limiters are applied, negative values could |
|---|
| 257 | C potentially be generated when large time step is used and when the |
|---|
| 258 | C initial fields contain discontinuities. |
|---|
| 259 | C This does not necessarily imply the integration is unstable. |
|---|
| 260 | C These negatives are typically very small. A filling algorithm is |
|---|
| 261 | C activated if the user set "fill" to be true. |
|---|
| 262 | C |
|---|
| 263 | C The van Leer scheme used here is nearly as accurate as the original PPM |
|---|
| 264 | C due to the use of a 4th order accurate reference slope. The PPM imple- |
|---|
| 265 | C mented here is an improvement over the original and is also based on |
|---|
| 266 | C the 4th order reference slope. |
|---|
| 267 | C |
|---|
| 268 | C ****6***0*********0*********0*********0*********0*********0**********72 |
|---|
| 269 | C |
|---|
| 270 | C User modifiable parameters |
|---|
| 271 | C |
|---|
| 272 | integer,parameter :: Jmax = 361, kmax = 150 |
|---|
| 273 | C |
|---|
| 274 | C ****6***0*********0*********0*********0*********0*********0**********72 |
|---|
| 275 | C |
|---|
| 276 | C Input-Output arrays |
|---|
| 277 | C |
|---|
| 278 | |
|---|
| 279 | real Q(IMR,JNP,NLAY,NC),PS1(IMR,JNP),PS2(IMR,JNP), |
|---|
| 280 | & U(IMR,JNP,NLAY),V(IMR,JNP,NLAY),AP(NLAY+1), |
|---|
| 281 | & BP(NLAY+1),W(IMR,JNP,NLAY),NDT,val(NLAY),Umax |
|---|
| 282 | integer IGD,IORD,JORD,KORD,NC,IMR,JNP,j1,NLAY,AE |
|---|
| 283 | integer IMRD2 |
|---|
| 284 | real PT |
|---|
| 285 | logical cross, fill, dum |
|---|
| 286 | C |
|---|
| 287 | C Local dynamic arrays |
|---|
| 288 | C |
|---|
| 289 | real CRX(IMR,JNP),CRY(IMR,JNP),xmass(IMR,JNP),ymass(IMR,JNP), |
|---|
| 290 | & fx1(IMR+1),DPI(IMR,JNP,NLAY),delp1(IMR,JNP,NLAY), |
|---|
| 291 | & WK1(IMR,JNP,NLAY),PU(IMR,JNP),PV(IMR,JNP),DC2(IMR,JNP), |
|---|
| 292 | & delp2(IMR,JNP,NLAY),DQ(IMR,JNP,NLAY,NC),VA(IMR,JNP), |
|---|
| 293 | & UA(IMR,JNP),qtmp(-IMR:2*IMR) |
|---|
| 294 | C |
|---|
| 295 | C Local static arrays |
|---|
| 296 | C |
|---|
| 297 | real DTDX(Jmax), DTDX5(Jmax), acosp(Jmax), |
|---|
| 298 | & cosp(Jmax), cose(Jmax), DAP(kmax),DBK(Kmax) |
|---|
| 299 | data NDT0, NSTEP /0, 0/ |
|---|
| 300 | data cross /.true./ |
|---|
| 301 | REAL DTDY, DTDY5, RCAP |
|---|
| 302 | INTEGER JS0, JN0, IML, JMR, IMJM |
|---|
| 303 | SAVE DTDY, DTDY5, RCAP, JS0, JN0, IML, |
|---|
| 304 | & DTDX, DTDX5, ACOSP, COSP, COSE, DAP,DBK |
|---|
| 305 | C |
|---|
| 306 | INTEGER NDT0, NSTEP, j2, k,j,i,ic,l,JS,JN,IMH |
|---|
| 307 | INTEGER IU,IIU,JT,iad,jad,krd |
|---|
| 308 | REAL r23,r3,PI,DL,DP,DT,CR1,MAXDT,ZTC,D5 |
|---|
| 309 | REAL sum1,sum2,ru |
|---|
| 310 | |
|---|
| 311 | JMR = JNP -1 |
|---|
| 312 | IMJM = IMR*JNP |
|---|
| 313 | j2 = JNP - j1 + 1 |
|---|
| 314 | NSTEP = NSTEP + 1 |
|---|
| 315 | C |
|---|
| 316 | C *********** Initialization ********************** |
|---|
| 317 | if(NSTEP==1) then |
|---|
| 318 | c |
|---|
| 319 | write(6,*) '------------------------------------ ' |
|---|
| 320 | write(6,*) 'NASA/GSFC Transport Core Version 4.5' |
|---|
| 321 | write(6,*) '------------------------------------ ' |
|---|
| 322 | c |
|---|
| 323 | WRITE(6,*) 'IMR=',IMR,' JNP=',JNP,' NLAY=',NLAY,' j1=',j1 |
|---|
| 324 | WRITE(6,*) 'NC=',NC,IORD,JORD,KORD,NDT |
|---|
| 325 | C |
|---|
| 326 | C controles sur les parametres |
|---|
| 327 | if(NLAY<6) then |
|---|
| 328 | write(6,*) 'NLAY must be >= 6' |
|---|
| 329 | stop |
|---|
| 330 | endif |
|---|
| 331 | if (JNP<NLAY) then |
|---|
| 332 | write(6,*) 'JNP must be >= NLAY' |
|---|
| 333 | stop |
|---|
| 334 | endif |
|---|
| 335 | IMRD2=mod(IMR,2) |
|---|
| 336 | if (j1==2.and.IMRD2/=0) then |
|---|
| 337 | write(6,*) 'if j1=2 IMR must be an even integer' |
|---|
| 338 | stop |
|---|
| 339 | endif |
|---|
| 340 | |
|---|
| 341 | C |
|---|
| 342 | if(Jmax<JNP .or. Kmax<NLAY) then |
|---|
| 343 | write(6,*) 'Jmax or Kmax is too small' |
|---|
| 344 | stop |
|---|
| 345 | endif |
|---|
| 346 | C |
|---|
| 347 | DO k=1,NLAY |
|---|
| 348 | DAP(k) = (AP(k+1) - AP(k))*PT |
|---|
| 349 | DBK(k) = BP(k+1) - BP(k) |
|---|
| 350 | ENDDO |
|---|
| 351 | C |
|---|
| 352 | PI = 4. * ATAN(1.) |
|---|
| 353 | DL = 2.*PI / REAL(IMR) |
|---|
| 354 | DP = PI / REAL(JMR) |
|---|
| 355 | C |
|---|
| 356 | if(IGD==0) then |
|---|
| 357 | C Compute analytic cosine at cell edges |
|---|
| 358 | call cosa(cosp,cose,JNP,PI,DP) |
|---|
| 359 | else |
|---|
| 360 | C Define cosine consistent with GEOS-GCM (using dycore2.0 or later) |
|---|
| 361 | call cosc(cosp,cose,JNP,PI,DP) |
|---|
| 362 | endif |
|---|
| 363 | C |
|---|
| 364 | do J=2,JMR |
|---|
| 365 | acosp(j) = 1. / cosp(j) |
|---|
| 366 | END DO |
|---|
| 367 | C |
|---|
| 368 | C Inverse of the Scaled polar cap area. |
|---|
| 369 | C |
|---|
| 370 | RCAP = DP / (IMR*(1.-COS((j1-1.5)*DP))) |
|---|
| 371 | acosp(1) = RCAP |
|---|
| 372 | acosp(JNP) = RCAP |
|---|
| 373 | endif |
|---|
| 374 | C |
|---|
| 375 | if(NDT0 /= NDT) then |
|---|
| 376 | DT = NDT |
|---|
| 377 | NDT0 = NDT |
|---|
| 378 | |
|---|
| 379 | if(Umax < 180.) then |
|---|
| 380 | write(6,*) 'Umax may be too small!' |
|---|
| 381 | endif |
|---|
| 382 | CR1 = abs(Umax*DT)/(DL*AE) |
|---|
| 383 | MaxDT = DP*AE / abs(Umax) + 0.5 |
|---|
| 384 | write(6,*)'Largest time step for max(V)=',Umax,' is ',MaxDT |
|---|
| 385 | if(MaxDT < abs(NDT)) then |
|---|
| 386 | write(6,*) 'Warning!!! NDT maybe too large!' |
|---|
| 387 | endif |
|---|
| 388 | C |
|---|
| 389 | if(CR1>=0.95) then |
|---|
| 390 | JS0 = 0 |
|---|
| 391 | JN0 = 0 |
|---|
| 392 | IML = IMR-2 |
|---|
| 393 | ZTC = 0. |
|---|
| 394 | else |
|---|
| 395 | ZTC = acos(CR1) * (180./PI) |
|---|
| 396 | C |
|---|
| 397 | JS0 = REAL(JMR)*(90.-ZTC)/180. + 2 |
|---|
| 398 | JS0 = max(JS0, J1+1) |
|---|
| 399 | IML = min(6*JS0/(J1-1)+2, 4*IMR/5) |
|---|
| 400 | JN0 = JNP-JS0+1 |
|---|
| 401 | endif |
|---|
| 402 | C |
|---|
| 403 | C |
|---|
| 404 | do J=2,JMR |
|---|
| 405 | DTDX(j) = DT / ( DL*AE*COSP(J) ) |
|---|
| 406 | |
|---|
| 407 | c print*,'dtdx=',dtdx(j) |
|---|
| 408 | DTDX5(j) = 0.5*DTDX(j) |
|---|
| 409 | enddo |
|---|
| 410 | C |
|---|
| 411 | |
|---|
| 412 | DTDY = DT /(AE*DP) |
|---|
| 413 | c print*,'dtdy=',dtdy |
|---|
| 414 | DTDY5 = 0.5*DTDY |
|---|
| 415 | C |
|---|
| 416 | c write(6,*) 'J1=',J1,' J2=', J2 |
|---|
| 417 | endif |
|---|
| 418 | C |
|---|
| 419 | C *********** End Initialization ********************** |
|---|
| 420 | C |
|---|
| 421 | C delp = pressure thickness: the psudo-density in a hydrostatic system. |
|---|
| 422 | do k=1,NLAY |
|---|
| 423 | do j=1,JNP |
|---|
| 424 | do i=1,IMR |
|---|
| 425 | delp1(i,j,k)=DAP(k)+DBK(k)*PS1(i,j) |
|---|
| 426 | delp2(i,j,k)=DAP(k)+DBK(k)*PS2(i,j) |
|---|
| 427 | enddo |
|---|
| 428 | enddo |
|---|
| 429 | enddo |
|---|
| 430 | |
|---|
| 431 | C |
|---|
| 432 | if(j1/=2) then |
|---|
| 433 | DO IC=1,NC |
|---|
| 434 | DO L=1,NLAY |
|---|
| 435 | DO I=1,IMR |
|---|
| 436 | Q(I, 2,L,IC) = Q(I, 1,L,IC) |
|---|
| 437 | Q(I,JMR,L,IC) = Q(I,JNP,L,IC) |
|---|
| 438 | END DO |
|---|
| 439 | END DO |
|---|
| 440 | END DO |
|---|
| 441 | endif |
|---|
| 442 | C |
|---|
| 443 | C Compute "tracer density" |
|---|
| 444 | DO IC=1,NC |
|---|
| 445 | DO k=1,NLAY |
|---|
| 446 | DO j=1,JNP |
|---|
| 447 | DO i=1,IMR |
|---|
| 448 | DQ(i,j,k,IC) = Q(i,j,k,IC)*delp1(i,j,k) |
|---|
| 449 | END DO |
|---|
| 450 | END DO |
|---|
| 451 | END DO |
|---|
| 452 | END DO |
|---|
| 453 | C |
|---|
| 454 | do k=1,NLAY |
|---|
| 455 | C |
|---|
| 456 | if(IGD==0) then |
|---|
| 457 | C Convert winds on A-Grid to Courant # on C-Grid. |
|---|
| 458 | call A2C(U(1,1,k),V(1,1,k),IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5) |
|---|
| 459 | else |
|---|
| 460 | C Convert winds on C-grid to Courant # |
|---|
| 461 | do j=j1,j2 |
|---|
| 462 | do i=2,IMR |
|---|
| 463 | CRX(i,J) = dtdx(j)*U(i-1,j,k) |
|---|
| 464 | END DO |
|---|
| 465 | END DO |
|---|
| 466 | |
|---|
| 467 | C |
|---|
| 468 | do j=j1,j2 |
|---|
| 469 | CRX(1,J) = dtdx(j)*U(IMR,j,k) |
|---|
| 470 | END DO |
|---|
| 471 | C |
|---|
| 472 | do i=1,IMR*JMR |
|---|
| 473 | CRY(i,2) = DTDY*V(i,1,k) |
|---|
| 474 | END DO |
|---|
| 475 | endif |
|---|
| 476 | C |
|---|
| 477 | C Determine JS and JN |
|---|
| 478 | JS = j1 |
|---|
| 479 | JN = j2 |
|---|
| 480 | C |
|---|
| 481 | do j=JS0,j1+1,-1 |
|---|
| 482 | do i=1,IMR |
|---|
| 483 | if(abs(CRX(i,j))>1.) then |
|---|
| 484 | JS = j |
|---|
| 485 | go to 2222 |
|---|
| 486 | endif |
|---|
| 487 | enddo |
|---|
| 488 | enddo |
|---|
| 489 | C |
|---|
| 490 | 2222 continue |
|---|
| 491 | do j=JN0,j2-1 |
|---|
| 492 | do i=1,IMR |
|---|
| 493 | if(abs(CRX(i,j))>1.) then |
|---|
| 494 | JN = j |
|---|
| 495 | go to 2233 |
|---|
| 496 | endif |
|---|
| 497 | enddo |
|---|
| 498 | enddo |
|---|
| 499 | 2233 continue |
|---|
| 500 | C |
|---|
| 501 | if(j1/=2) then ! Enlarged polar cap. |
|---|
| 502 | do i=1,IMR |
|---|
| 503 | DPI(i, 2,k) = 0. |
|---|
| 504 | DPI(i,JMR,k) = 0. |
|---|
| 505 | enddo |
|---|
| 506 | endif |
|---|
| 507 | C |
|---|
| 508 | C ******* Compute horizontal mass fluxes ************ |
|---|
| 509 | C |
|---|
| 510 | C N-S component |
|---|
| 511 | do j=j1,j2+1 |
|---|
| 512 | D5 = 0.5 * COSE(j) |
|---|
| 513 | do i=1,IMR |
|---|
| 514 | ymass(i,j) = CRY(i,j)*D5*(delp2(i,j,k) + delp2(i,j-1,k)) |
|---|
| 515 | enddo |
|---|
| 516 | enddo |
|---|
| 517 | C |
|---|
| 518 | do j=j1,j2 |
|---|
| 519 | DO i=1,IMR |
|---|
| 520 | DPI(i,j,k) = (ymass(i,j) - ymass(i,j+1)) * acosp(j) |
|---|
| 521 | END DO |
|---|
| 522 | END DO |
|---|
| 523 | C |
|---|
| 524 | C Poles |
|---|
| 525 | sum1 = ymass(IMR,j1 ) |
|---|
| 526 | sum2 = ymass(IMR,J2+1) |
|---|
| 527 | do i=1,IMR-1 |
|---|
| 528 | sum1 = sum1 + ymass(i,j1 ) |
|---|
| 529 | sum2 = sum2 + ymass(i,J2+1) |
|---|
| 530 | enddo |
|---|
| 531 | C |
|---|
| 532 | sum1 = - sum1 * RCAP |
|---|
| 533 | sum2 = sum2 * RCAP |
|---|
| 534 | do i=1,IMR |
|---|
| 535 | DPI(i, 1,k) = sum1 |
|---|
| 536 | DPI(i,JNP,k) = sum2 |
|---|
| 537 | enddo |
|---|
| 538 | C |
|---|
| 539 | C E-W component |
|---|
| 540 | C |
|---|
| 541 | do j=j1,j2 |
|---|
| 542 | do i=2,IMR |
|---|
| 543 | PU(i,j) = 0.5 * (delp2(i,j,k) + delp2(i-1,j,k)) |
|---|
| 544 | enddo |
|---|
| 545 | enddo |
|---|
| 546 | C |
|---|
| 547 | do j=j1,j2 |
|---|
| 548 | PU(1,j) = 0.5 * (delp2(1,j,k) + delp2(IMR,j,k)) |
|---|
| 549 | enddo |
|---|
| 550 | C |
|---|
| 551 | do j=j1,j2 |
|---|
| 552 | DO i=1,IMR |
|---|
| 553 | xmass(i,j) = PU(i,j)*CRX(i,j) |
|---|
| 554 | END DO |
|---|
| 555 | END DO |
|---|
| 556 | C |
|---|
| 557 | DO j=j1,j2 |
|---|
| 558 | DO i=1,IMR-1 |
|---|
| 559 | DPI(i,j,k) = DPI(i,j,k) + xmass(i,j) - xmass(i+1,j) |
|---|
| 560 | END DO |
|---|
| 561 | END DO |
|---|
| 562 | C |
|---|
| 563 | DO j=j1,j2 |
|---|
| 564 | DPI(IMR,j,k) = DPI(IMR,j,k) + xmass(IMR,j) - xmass(1,j) |
|---|
| 565 | END DO |
|---|
| 566 | C |
|---|
| 567 | DO j=j1,j2 |
|---|
| 568 | do i=1,IMR-1 |
|---|
| 569 | UA(i,j) = 0.5 * (CRX(i,j)+CRX(i+1,j)) |
|---|
| 570 | enddo |
|---|
| 571 | enddo |
|---|
| 572 | C |
|---|
| 573 | DO j=j1,j2 |
|---|
| 574 | UA(imr,j) = 0.5 * (CRX(imr,j)+CRX(1,j)) |
|---|
| 575 | enddo |
|---|
| 576 | ccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 577 | c Rajouts pour LMDZ.3.3 |
|---|
| 578 | ccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 579 | do i=1,IMR |
|---|
| 580 | do j=1,JNP |
|---|
| 581 | VA(i,j)=0. |
|---|
| 582 | enddo |
|---|
| 583 | enddo |
|---|
| 584 | |
|---|
| 585 | do i=1,imr*(JMR-1) |
|---|
| 586 | VA(i,2) = 0.5*(CRY(i,2)+CRY(i,3)) |
|---|
| 587 | enddo |
|---|
| 588 | C |
|---|
| 589 | if(j1==2) then |
|---|
| 590 | IMH = IMR/2 |
|---|
| 591 | do i=1,IMH |
|---|
| 592 | VA(i, 1) = 0.5*(CRY(i,2)-CRY(i+IMH,2)) |
|---|
| 593 | VA(i+IMH, 1) = -VA(i,1) |
|---|
| 594 | VA(i, JNP) = 0.5*(CRY(i,JNP)-CRY(i+IMH,JMR)) |
|---|
| 595 | VA(i+IMH,JNP) = -VA(i,JNP) |
|---|
| 596 | enddo |
|---|
| 597 | VA(IMR,1)=VA(1,1) |
|---|
| 598 | VA(IMR,JNP)=VA(1,JNP) |
|---|
| 599 | endif |
|---|
| 600 | C |
|---|
| 601 | C ****6***0*********0*********0*********0*********0*********0**********72 |
|---|
| 602 | do IC=1,NC |
|---|
| 603 | C |
|---|
| 604 | do i=1,IMJM |
|---|
| 605 | wk1(i,1,1) = 0. |
|---|
| 606 | wk1(i,1,2) = 0. |
|---|
| 607 | enddo |
|---|
| 608 | C |
|---|
| 609 | C E-W advective cross term |
|---|
| 610 | do j=J1,J2 |
|---|
| 611 | if(J>JS .and. J<JN) GO TO 250 |
|---|
| 612 | C |
|---|
| 613 | do i=1,IMR |
|---|
| 614 | qtmp(i) = q(i,j,k,IC) |
|---|
| 615 | enddo |
|---|
| 616 | C |
|---|
| 617 | do i=-IML,0 |
|---|
| 618 | qtmp(i) = q(IMR+i,j,k,IC) |
|---|
| 619 | qtmp(IMR+1-i) = q(1-i,j,k,IC) |
|---|
| 620 | enddo |
|---|
| 621 | C |
|---|
| 622 | DO i=1,IMR |
|---|
| 623 | iu = UA(i,j) |
|---|
| 624 | ru = UA(i,j) - iu |
|---|
| 625 | iiu = i-iu |
|---|
| 626 | if(UA(i,j)>=0.) then |
|---|
| 627 | wk1(i,j,1) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu)) |
|---|
| 628 | else |
|---|
| 629 | wk1(i,j,1) = qtmp(iiu)+ru*(qtmp(iiu)-qtmp(iiu+1)) |
|---|
| 630 | endif |
|---|
| 631 | wk1(i,j,1) = wk1(i,j,1) - qtmp(i) |
|---|
| 632 | END DO |
|---|
| 633 | 250 continue |
|---|
| 634 | END DO |
|---|
| 635 | C |
|---|
| 636 | if(JN/=0) then |
|---|
| 637 | do j=JS+1,JN-1 |
|---|
| 638 | C |
|---|
| 639 | do i=1,IMR |
|---|
| 640 | qtmp(i) = q(i,j,k,IC) |
|---|
| 641 | enddo |
|---|
| 642 | C |
|---|
| 643 | qtmp(0) = q(IMR,J,k,IC) |
|---|
| 644 | qtmp(IMR+1) = q( 1,J,k,IC) |
|---|
| 645 | C |
|---|
| 646 | do i=1,imr |
|---|
| 647 | iu = i - UA(i,j) |
|---|
| 648 | wk1(i,j,1) = UA(i,j)*(qtmp(iu) - qtmp(iu+1)) |
|---|
| 649 | enddo |
|---|
| 650 | enddo |
|---|
| 651 | endif |
|---|
| 652 | C ****6***0*********0*********0*********0*********0*********0**********72 |
|---|
| 653 | C Contribution from the N-S advection |
|---|
| 654 | do i=1,imr*(j2-j1+1) |
|---|
| 655 | JT = REAL(J1) - VA(i,j1) |
|---|
| 656 | wk1(i,j1,2) = VA(i,j1) * (q(i,jt,k,IC) - q(i,jt+1,k,IC)) |
|---|
| 657 | enddo |
|---|
| 658 | C |
|---|
| 659 | do i=1,IMJM |
|---|
| 660 | wk1(i,1,1) = q(i,1,k,IC) + 0.5*wk1(i,1,1) |
|---|
| 661 | wk1(i,1,2) = q(i,1,k,IC) + 0.5*wk1(i,1,2) |
|---|
| 662 | enddo |
|---|
| 663 | C |
|---|
| 664 | if(cross) then |
|---|
| 665 | C Add cross terms in the vertical direction. |
|---|
| 666 | if(IORD >= 2) then |
|---|
| 667 | iad = 2 |
|---|
| 668 | else |
|---|
| 669 | iad = 1 |
|---|
| 670 | endif |
|---|
| 671 | C |
|---|
| 672 | if(JORD >= 2) then |
|---|
| 673 | jad = 2 |
|---|
| 674 | else |
|---|
| 675 | jad = 1 |
|---|
| 676 | endif |
|---|
| 677 | call xadv(IMR,JNP,j1,j2,wk1(1,1,2),UA,JS,JN,IML,DC2,iad) |
|---|
| 678 | call yadv(IMR,JNP,j1,j2,wk1(1,1,1),VA,PV,W,jad) |
|---|
| 679 | do j=1,JNP |
|---|
| 680 | do i=1,IMR |
|---|
| 681 | q(i,j,k,IC) = q(i,j,k,IC) + DC2(i,j) + PV(i,j) |
|---|
| 682 | enddo |
|---|
| 683 | enddo |
|---|
| 684 | endif |
|---|
| 685 | C |
|---|
| 686 | call xtp(IMR,JNP,IML,j1,j2,JN,JS,PU,DQ(1,1,k,IC),wk1(1,1,2) |
|---|
| 687 | & ,CRX,fx1,xmass,IORD) |
|---|
| 688 | |
|---|
| 689 | call ytp(IMR,JNP,j1,j2,acosp,RCAP,DQ(1,1,k,IC),wk1(1,1,1),CRY, |
|---|
| 690 | & DC2,ymass,WK1(1,1,3),wk1(1,1,4),WK1(1,1,5),WK1(1,1,6),JORD) |
|---|
| 691 | C |
|---|
| 692 | END DO |
|---|
| 693 | END DO |
|---|
| 694 | C |
|---|
| 695 | C ******* Compute vertical mass flux (same unit as PS) *********** |
|---|
| 696 | C |
|---|
| 697 | C 1st step: compute total column mass CONVERGENCE. |
|---|
| 698 | C |
|---|
| 699 | do j=1,JNP |
|---|
| 700 | do i=1,IMR |
|---|
| 701 | CRY(i,j) = DPI(i,j,1) |
|---|
| 702 | END DO |
|---|
| 703 | END DO |
|---|
| 704 | C |
|---|
| 705 | do k=2,NLAY |
|---|
| 706 | do j=1,JNP |
|---|
| 707 | do i=1,IMR |
|---|
| 708 | CRY(i,j) = CRY(i,j) + DPI(i,j,k) |
|---|
| 709 | END DO |
|---|
| 710 | END DO |
|---|
| 711 | END DO |
|---|
| 712 | C |
|---|
| 713 | do j=1,JNP |
|---|
| 714 | do i=1,IMR |
|---|
| 715 | C |
|---|
| 716 | C 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption. |
|---|
| 717 | C Changes (increases) to surface pressure = total column mass convergence |
|---|
| 718 | C |
|---|
| 719 | PS2(i,j) = PS1(i,j) + CRY(i,j) |
|---|
| 720 | C |
|---|
| 721 | C 3rd step: compute vertical mass flux from mass conservation principle. |
|---|
| 722 | C |
|---|
| 723 | W(i,j,1) = DPI(i,j,1) - DBK(1)*CRY(i,j) |
|---|
| 724 | W(i,j,NLAY) = 0. |
|---|
| 725 | END DO |
|---|
| 726 | END DO |
|---|
| 727 | C |
|---|
| 728 | do k=2,NLAY-1 |
|---|
| 729 | do j=1,JNP |
|---|
| 730 | do i=1,IMR |
|---|
| 731 | W(i,j,k) = W(i,j,k-1) + DPI(i,j,k) - DBK(k)*CRY(i,j) |
|---|
| 732 | END DO |
|---|
| 733 | END DO |
|---|
| 734 | END DO |
|---|
| 735 | C |
|---|
| 736 | DO k=1,NLAY |
|---|
| 737 | DO j=1,JNP |
|---|
| 738 | DO i=1,IMR |
|---|
| 739 | delp2(i,j,k) = DAP(k) + DBK(k)*PS2(i,j) |
|---|
| 740 | END DO |
|---|
| 741 | END DO |
|---|
| 742 | END DO |
|---|
| 743 | C |
|---|
| 744 | KRD = max(3, KORD) |
|---|
| 745 | do IC=1,NC |
|---|
| 746 | C |
|---|
| 747 | C****6***0*********0*********0*********0*********0*********0**********72 |
|---|
| 748 | |
|---|
| 749 | call FZPPM(IMR,JNP,NLAY,j1,DQ(1,1,1,IC),W,Q(1,1,1,IC),WK1,DPI, |
|---|
| 750 | & DC2,CRX,CRY,PU,PV,xmass,ymass,delp1,KRD) |
|---|
| 751 | C |
|---|
| 752 | |
|---|
| 753 | if(fill) call qckxyz(DQ(1,1,1,IC),DC2,IMR,JNP,NLAY,j1,j2, |
|---|
| 754 | & cosp,acosp,.false.,IC,NSTEP) |
|---|
| 755 | C |
|---|
| 756 | C Recover tracer mixing ratio from "density" using predicted |
|---|
| 757 | C "air density" (pressure thickness) at time-level n+1 |
|---|
| 758 | C |
|---|
| 759 | DO k=1,NLAY |
|---|
| 760 | DO j=1,JNP |
|---|
| 761 | DO i=1,IMR |
|---|
| 762 | Q(i,j,k,IC) = DQ(i,j,k,IC) / delp2(i,j,k) |
|---|
| 763 | c print*,'i=',i,'j=',j,'k=',k,'Q(i,j,k,IC)=',Q(i,j,k,IC) |
|---|
| 764 | enddo |
|---|
| 765 | enddo |
|---|
| 766 | enddo |
|---|
| 767 | C |
|---|
| 768 | if(j1/=2) then |
|---|
| 769 | DO k=1,NLAY |
|---|
| 770 | DO I=1,IMR |
|---|
| 771 | c j=1 c'est le pôle Sud, j=JNP c'est le pôle Nord |
|---|
| 772 | Q(I, 2,k,IC) = Q(I, 1,k,IC) |
|---|
| 773 | Q(I,JMR,k,IC) = Q(I,JNP,k,IC) |
|---|
| 774 | END DO |
|---|
| 775 | END DO |
|---|
| 776 | endif |
|---|
| 777 | END DO |
|---|
| 778 | C |
|---|
| 779 | if(j1/=2) then |
|---|
| 780 | DO k=1,NLAY |
|---|
| 781 | DO i=1,IMR |
|---|
| 782 | W(i, 2,k) = W(i, 1,k) |
|---|
| 783 | W(i,JMR,k) = W(i,JNP,k) |
|---|
| 784 | END DO |
|---|
| 785 | END DO |
|---|
| 786 | endif |
|---|
| 787 | C |
|---|
| 788 | RETURN |
|---|
| 789 | END |
|---|
| 790 | C |
|---|
| 791 | C****6***0*********0*********0*********0*********0*********0**********72 |
|---|
| 792 | subroutine FZPPM(IMR,JNP,NLAY,j1,DQ,WZ,P,DC,DQDT,AR,AL,A6, |
|---|
| 793 | & flux,wk1,wk2,wz2,delp,KORD) |
|---|
| 794 | implicit none |
|---|
| 795 | integer,parameter :: kmax = 150 |
|---|
| 796 | real,parameter :: R23 = 2./3., R3 = 1./3. |
|---|
| 797 | integer IMR,JNP,NLAY,J1,KORD |
|---|
| 798 | real WZ(IMR,JNP,NLAY),P(IMR,JNP,NLAY),DC(IMR,JNP,NLAY), |
|---|
| 799 | & wk1(IMR,*),delp(IMR,JNP,NLAY),DQ(IMR,JNP,NLAY), |
|---|
| 800 | & DQDT(IMR,JNP,NLAY) |
|---|
| 801 | C Assuming JNP >= NLAY |
|---|
| 802 | real AR(IMR,*),AL(IMR,*),A6(IMR,*),flux(IMR,*),wk2(IMR,*), |
|---|
| 803 | & wz2(IMR,*) |
|---|
| 804 | integer JMR,IMJM,NLAYM1,LMT,K,I,J |
|---|
| 805 | real c0,c1,c2,tmp,qmax,qmin,a,b,fct,a1,a2,cm,cp |
|---|
| 806 | C |
|---|
| 807 | JMR = JNP - 1 |
|---|
| 808 | IMJM = IMR*JNP |
|---|
| 809 | NLAYM1 = NLAY - 1 |
|---|
| 810 | C |
|---|
| 811 | LMT = KORD - 3 |
|---|
| 812 | C |
|---|
| 813 | C ****6***0*********0*********0*********0*********0*********0**********72 |
|---|
| 814 | C Compute DC for PPM |
|---|
| 815 | C ****6***0*********0*********0*********0*********0*********0**********72 |
|---|
| 816 | C |
|---|
| 817 | do k=1,NLAYM1 |
|---|
| 818 | do i=1,IMJM |
|---|
| 819 | DQDT(i,1,k) = P(i,1,k+1) - P(i,1,k) |
|---|
| 820 | END DO |
|---|
| 821 | END DO |
|---|
| 822 | C |
|---|
| 823 | DO k=2,NLAYM1 |
|---|
| 824 | DO I=1,IMJM |
|---|
| 825 | c0 = delp(i,1,k) / (delp(i,1,k-1)+delp(i,1,k)+delp(i,1,k+1)) |
|---|
| 826 | c1 = (delp(i,1,k-1)+0.5*delp(i,1,k))/(delp(i,1,k+1)+delp(i,1,k)) |
|---|
| 827 | c2 = (delp(i,1,k+1)+0.5*delp(i,1,k))/(delp(i,1,k-1)+delp(i,1,k)) |
|---|
| 828 | tmp = c0*(c1*DQDT(i,1,k) + c2*DQDT(i,1,k-1)) |
|---|
| 829 | Qmax = max(P(i,1,k-1),P(i,1,k),P(i,1,k+1)) - P(i,1,k) |
|---|
| 830 | Qmin = P(i,1,k) - min(P(i,1,k-1),P(i,1,k),P(i,1,k+1)) |
|---|
| 831 | DC(i,1,k) = sign(min(abs(tmp),Qmax,Qmin), tmp) |
|---|
| 832 | END DO |
|---|
| 833 | END DO |
|---|
| 834 | |
|---|
| 835 | C |
|---|
| 836 | C ****6***0*********0*********0*********0*********0*********0**********72 |
|---|
| 837 | C Loop over latitudes (to save memory) |
|---|
| 838 | C ****6***0*********0*********0*********0*********0*********0**********72 |
|---|
| 839 | C |
|---|
| 840 | DO j=1,JNP |
|---|
| 841 | if((j==2 .or. j==JMR) .and. j1/=2) goto 2000 |
|---|
| 842 | C |
|---|
| 843 | DO k=1,NLAY |
|---|
| 844 | DO i=1,IMR |
|---|
| 845 | wz2(i,k) = WZ(i,j,k) |
|---|
| 846 | wk1(i,k) = P(i,j,k) |
|---|
| 847 | wk2(i,k) = delp(i,j,k) |
|---|
| 848 | flux(i,k) = DC(i,j,k) !this flux is actually the monotone slope |
|---|
| 849 | enddo |
|---|
| 850 | enddo |
|---|
| 851 | C |
|---|
| 852 | C****6***0*********0*********0*********0*********0*********0**********72 |
|---|
| 853 | C Compute first guesses at cell interfaces |
|---|
| 854 | C First guesses are required to be continuous. |
|---|
| 855 | C ****6***0*********0*********0*********0*********0*********0**********72 |
|---|
| 856 | C |
|---|
| 857 | C three-cell parabolic subgrid distribution at model top |
|---|
| 858 | C two-cell parabolic with zero gradient subgrid distribution |
|---|
| 859 | C at the surface. |
|---|
| 860 | C |
|---|
| 861 | C First guess top edge value |
|---|
| 862 | DO i=1,IMR |
|---|
| 863 | C three-cell PPM |
|---|
| 864 | C Compute a,b, and c of q = aP**2 + bP + c using cell averages and delp |
|---|
| 865 | a = 3.*( DQDT(i,j,2) - DQDT(i,j,1)*(wk2(i,2)+wk2(i,3))/ |
|---|
| 866 | & (wk2(i,1)+wk2(i,2)) ) / |
|---|
| 867 | & ( (wk2(i,2)+wk2(i,3))*(wk2(i,1)+wk2(i,2)+wk2(i,3)) ) |
|---|
| 868 | b = 2.*DQDT(i,j,1)/(wk2(i,1)+wk2(i,2)) - |
|---|
| 869 | & R23*a*(2.*wk2(i,1)+wk2(i,2)) |
|---|
| 870 | AL(i,1) = wk1(i,1) - wk2(i,1)*(R3*a*wk2(i,1) + 0.5*b) |
|---|
| 871 | AL(i,2) = wk2(i,1)*(a*wk2(i,1) + b) + AL(i,1) |
|---|
| 872 | C |
|---|
| 873 | C Check if change sign |
|---|
| 874 | if(wk1(i,1)*AL(i,1)<=0.) then |
|---|
| 875 | AL(i,1) = 0. |
|---|
| 876 | flux(i,1) = 0. |
|---|
| 877 | else |
|---|
| 878 | flux(i,1) = wk1(i,1) - AL(i,1) |
|---|
| 879 | endif |
|---|
| 880 | END DO |
|---|
| 881 | C |
|---|
| 882 | C Bottom |
|---|
| 883 | DO i=1,IMR |
|---|
| 884 | C 2-cell PPM with zero gradient right at the surface |
|---|
| 885 | C |
|---|
| 886 | fct = DQDT(i,j,NLAYM1)*wk2(i,NLAY)**2 / |
|---|
| 887 | & ( (wk2(i,NLAY)+wk2(i,NLAYM1))*(2.*wk2(i,NLAY)+wk2(i,NLAYM1))) |
|---|
| 888 | AR(i,NLAY) = wk1(i,NLAY) + fct |
|---|
| 889 | AL(i,NLAY) = wk1(i,NLAY) - (fct+fct) |
|---|
| 890 | if(wk1(i,NLAY)*AR(i,NLAY)<=0.) AR(i,NLAY) = 0. |
|---|
| 891 | flux(i,NLAY) = AR(i,NLAY) - wk1(i,NLAY) |
|---|
| 892 | END DO |
|---|
| 893 | |
|---|
| 894 | C |
|---|
| 895 | C****6***0*********0*********0*********0*********0*********0**********72 |
|---|
| 896 | C 4th order interpolation in the interior. |
|---|
| 897 | C****6***0*********0*********0*********0*********0*********0**********72 |
|---|
| 898 | C |
|---|
| 899 | DO k=3,NLAYM1 |
|---|
| 900 | DO i=1,IMR |
|---|
| 901 | c1 = DQDT(i,j,k-1)*wk2(i,k-1) / (wk2(i,k-1)+wk2(i,k)) |
|---|
| 902 | c2 = 2. / (wk2(i,k-2)+wk2(i,k-1)+wk2(i,k)+wk2(i,k+1)) |
|---|
| 903 | A1 = (wk2(i,k-2)+wk2(i,k-1)) / (2.*wk2(i,k-1)+wk2(i,k)) |
|---|
| 904 | A2 = (wk2(i,k )+wk2(i,k+1)) / (2.*wk2(i,k)+wk2(i,k-1)) |
|---|
| 905 | AL(i,k) = wk1(i,k-1) + c1 + c2 * |
|---|
| 906 | & ( wk2(i,k )*(c1*(A1 - A2)+A2*flux(i,k-1)) - |
|---|
| 907 | & wk2(i,k-1)*A1*flux(i,k) ) |
|---|
| 908 | C print *,'AL1',i,k, AL(i,k) |
|---|
| 909 | END DO |
|---|
| 910 | END DO |
|---|
| 911 | C |
|---|
| 912 | do i=1,IMR*NLAYM1 |
|---|
| 913 | AR(i,1) = AL(i,2) |
|---|
| 914 | C print *,'AR1',i,AR(i,1) |
|---|
| 915 | END DO |
|---|
| 916 | C |
|---|
| 917 | do i=1,IMR*NLAY |
|---|
| 918 | A6(i,1) = 3.*(wk1(i,1)+wk1(i,1) - (AL(i,1)+AR(i,1))) |
|---|
| 919 | C print *,'A61',i,A6(i,1) |
|---|
| 920 | END DO |
|---|
| 921 | C |
|---|
| 922 | C****6***0*********0*********0*********0*********0*********0**********72 |
|---|
| 923 | C Top & Bot always monotonic |
|---|
| 924 | call lmtppm(flux(1,1),A6(1,1),AR(1,1),AL(1,1),wk1(1,1),IMR,0) |
|---|
| 925 | call lmtppm(flux(1,NLAY),A6(1,NLAY),AR(1,NLAY),AL(1,NLAY), |
|---|
| 926 | & wk1(1,NLAY),IMR,0) |
|---|
| 927 | C |
|---|
| 928 | C Interior depending on KORD |
|---|
| 929 | if(LMT<=2) |
|---|
| 930 | & call lmtppm(flux(1,2),A6(1,2),AR(1,2),AL(1,2),wk1(1,2), |
|---|
| 931 | & IMR*(NLAY-2),LMT) |
|---|
| 932 | C |
|---|
| 933 | C****6***0*********0*********0*********0*********0*********0**********72 |
|---|
| 934 | C |
|---|
| 935 | DO i=1,IMR*NLAYM1 |
|---|
| 936 | IF(wz2(i,1)>0.) then |
|---|
| 937 | CM = wz2(i,1) / wk2(i,1) |
|---|
| 938 | flux(i,2) = AR(i,1)+0.5*CM*(AL(i,1)-AR(i,1)+A6(i,1)*(1.-R23*CM)) |
|---|
| 939 | else |
|---|
| 940 | C print *,'test2-0',i,j,wz2(i,1),wk2(i,2) |
|---|
| 941 | CP= wz2(i,1) / wk2(i,2) |
|---|
| 942 | C print *,'testCP',CP |
|---|
| 943 | flux(i,2) = AL(i,2)+0.5*CP*(AL(i,2)-AR(i,2)-A6(i,2)*(1.+R23*CP)) |
|---|
| 944 | C print *,'test2',i, AL(i,2),AR(i,2),A6(i,2),R23 |
|---|
| 945 | endif |
|---|
| 946 | END DO |
|---|
| 947 | C |
|---|
| 948 | DO i=1,IMR*NLAYM1 |
|---|
| 949 | flux(i,2) = wz2(i,1) * flux(i,2) |
|---|
| 950 | END DO |
|---|
| 951 | C |
|---|
| 952 | do i=1,IMR |
|---|
| 953 | DQ(i,j, 1) = DQ(i,j, 1) - flux(i, 2) |
|---|
| 954 | DQ(i,j,NLAY) = DQ(i,j,NLAY) + flux(i,NLAY) |
|---|
| 955 | END DO |
|---|
| 956 | C |
|---|
| 957 | do k=2,NLAYM1 |
|---|
| 958 | do i=1,IMR |
|---|
| 959 | DQ(i,j,k) = DQ(i,j,k) + flux(i,k) - flux(i,k+1) |
|---|
| 960 | END DO |
|---|
| 961 | END DO |
|---|
| 962 | 2000 continue |
|---|
| 963 | END DO |
|---|
| 964 | return |
|---|
| 965 | end |
|---|
| 966 | C |
|---|
| 967 | subroutine xtp(IMR,JNP,IML,j1,j2,JN,JS,PU,DQ,Q,UC, |
|---|
| 968 | & fx1,xmass,IORD) |
|---|
| 969 | implicit none |
|---|
| 970 | integer IMR,JNP,IML,j1,j2,JN,JS,IORD |
|---|
| 971 | real PU,DQ,Q,UC,fx1,xmass |
|---|
| 972 | real dc,qtmp |
|---|
| 973 | integer ISAVE(IMR) |
|---|
| 974 | dimension UC(IMR,*),DC(-IML:IMR+IML+1),xmass(IMR,JNP) |
|---|
| 975 | & ,fx1(IMR+1),DQ(IMR,JNP),qtmp(-IML:IMR+1+IML) |
|---|
| 976 | dimension PU(IMR,JNP),Q(IMR,JNP) |
|---|
| 977 | integer jvan,j1vl,j2vl,j,i,iu,itmp,ist,imp |
|---|
| 978 | real rut |
|---|
| 979 | C |
|---|
| 980 | IMP = IMR + 1 |
|---|
| 981 | C |
|---|
| 982 | C van Leer at high latitudes |
|---|
| 983 | jvan = max(1,JNP/18) |
|---|
| 984 | j1vl = j1+jvan |
|---|
| 985 | j2vl = j2-jvan |
|---|
| 986 | C |
|---|
| 987 | do j=j1,j2 |
|---|
| 988 | C |
|---|
| 989 | do i=1,IMR |
|---|
| 990 | qtmp(i) = q(i,j) |
|---|
| 991 | enddo |
|---|
| 992 | C |
|---|
| 993 | if(j>=JN .or. j<=JS) goto 2222 |
|---|
| 994 | C ************* Eulerian ********** |
|---|
| 995 | C |
|---|
| 996 | qtmp(0) = q(IMR,J) |
|---|
| 997 | qtmp(-1) = q(IMR-1,J) |
|---|
| 998 | qtmp(IMP) = q(1,J) |
|---|
| 999 | qtmp(IMP+1) = q(2,J) |
|---|
| 1000 | C |
|---|
| 1001 | IF(IORD==1 .or. j==j1 .or. j==j2) THEN |
|---|
| 1002 | DO i=1,IMR |
|---|
| 1003 | iu = REAL(i) - uc(i,j) |
|---|
| 1004 | fx1(i) = qtmp(iu) |
|---|
| 1005 | END DO |
|---|
| 1006 | ELSE |
|---|
| 1007 | call xmist(IMR,IML,Qtmp,DC) |
|---|
| 1008 | DC(0) = DC(IMR) |
|---|
| 1009 | C |
|---|
| 1010 | if(IORD==2 .or. j<=j1vl .or. j>=j2vl) then |
|---|
| 1011 | DO i=1,IMR |
|---|
| 1012 | iu = REAL(i) - uc(i,j) |
|---|
| 1013 | fx1(i) = qtmp(iu) + DC(iu)*(sign(1.,uc(i,j))-uc(i,j)) |
|---|
| 1014 | END DO |
|---|
| 1015 | else |
|---|
| 1016 | call fxppm(IMR,IML,UC(1,j),Qtmp,DC,fx1,IORD) |
|---|
| 1017 | endif |
|---|
| 1018 | C |
|---|
| 1019 | ENDIF |
|---|
| 1020 | C |
|---|
| 1021 | DO i=1,IMR |
|---|
| 1022 | fx1(i) = fx1(i)*xmass(i,j) |
|---|
| 1023 | END DO |
|---|
| 1024 | C |
|---|
| 1025 | goto 1309 |
|---|
| 1026 | C |
|---|
| 1027 | C ***** Conservative (flux-form) Semi-Lagrangian transport ***** |
|---|
| 1028 | C |
|---|
| 1029 | 2222 continue |
|---|
| 1030 | C |
|---|
| 1031 | do i=-IML,0 |
|---|
| 1032 | qtmp(i) = q(IMR+i,j) |
|---|
| 1033 | qtmp(IMP-i) = q(1-i,j) |
|---|
| 1034 | enddo |
|---|
| 1035 | C |
|---|
| 1036 | IF(IORD==1 .or. j==j1 .or. j==j2) THEN |
|---|
| 1037 | DO i=1,IMR |
|---|
| 1038 | itmp = INT(uc(i,j)) |
|---|
| 1039 | ISAVE(i) = i - itmp |
|---|
| 1040 | iu = i - uc(i,j) |
|---|
| 1041 | fx1(i) = (uc(i,j) - itmp)*qtmp(iu) |
|---|
| 1042 | END DO |
|---|
| 1043 | ELSE |
|---|
| 1044 | call xmist(IMR,IML,Qtmp,DC) |
|---|
| 1045 | C |
|---|
| 1046 | do i=-IML,0 |
|---|
| 1047 | DC(i) = DC(IMR+i) |
|---|
| 1048 | DC(IMP-i) = DC(1-i) |
|---|
| 1049 | enddo |
|---|
| 1050 | C |
|---|
| 1051 | DO i=1,IMR |
|---|
| 1052 | itmp = INT(uc(i,j)) |
|---|
| 1053 | rut = uc(i,j) - itmp |
|---|
| 1054 | ISAVE(i) = i - itmp |
|---|
| 1055 | iu = i - uc(i,j) |
|---|
| 1056 | fx1(i) = rut*(qtmp(iu) + DC(iu)*(sign(1.,rut) - rut)) |
|---|
| 1057 | END DO |
|---|
| 1058 | ENDIF |
|---|
| 1059 | C |
|---|
| 1060 | do i=1,IMR |
|---|
| 1061 | IF(uc(i,j)>1.) then |
|---|
| 1062 | CDIR$ NOVECTOR |
|---|
| 1063 | do ist = ISAVE(i),i-1 |
|---|
| 1064 | fx1(i) = fx1(i) + qtmp(ist) |
|---|
| 1065 | enddo |
|---|
| 1066 | elseIF(uc(i,j)<-1.) then |
|---|
| 1067 | do ist = i,ISAVE(i)-1 |
|---|
| 1068 | fx1(i) = fx1(i) - qtmp(ist) |
|---|
| 1069 | enddo |
|---|
| 1070 | CDIR$ VECTOR |
|---|
| 1071 | endif |
|---|
| 1072 | END DO |
|---|
| 1073 | do i=1,IMR |
|---|
| 1074 | fx1(i) = PU(i,j)*fx1(i) |
|---|
| 1075 | enddo |
|---|
| 1076 | C |
|---|
| 1077 | C *************************************** |
|---|
| 1078 | C |
|---|
| 1079 | 1309 fx1(IMP) = fx1(1) |
|---|
| 1080 | DO i=1,IMR |
|---|
| 1081 | DQ(i,j) = DQ(i,j) + fx1(i)-fx1(i+1) |
|---|
| 1082 | END DO |
|---|
| 1083 | C |
|---|
| 1084 | C *************************************** |
|---|
| 1085 | C |
|---|
| 1086 | END DO |
|---|
| 1087 | return |
|---|
| 1088 | end |
|---|
| 1089 | C |
|---|
| 1090 | subroutine fxppm(IMR,IML,UT,P,DC,flux,IORD) |
|---|
| 1091 | implicit none |
|---|
| 1092 | integer IMR,IML,IORD |
|---|
| 1093 | real UT,P,DC,flux |
|---|
| 1094 | real,parameter :: R3 = 1./3., R23 = 2./3. |
|---|
| 1095 | DIMENSION UT(*),flux(*),P(-IML:IMR+IML+1),DC(-IML:IMR+IML+1) |
|---|
| 1096 | REAL :: AR(0:IMR),AL(0:IMR),A6(0:IMR) |
|---|
| 1097 | integer LMT,IMP,JLVL,i |
|---|
| 1098 | c logical first |
|---|
| 1099 | c data first /.true./ |
|---|
| 1100 | c SAVE LMT |
|---|
| 1101 | c if(first) then |
|---|
| 1102 | C |
|---|
| 1103 | C correction calcul de LMT a chaque passage pour pouvoir choisir |
|---|
| 1104 | c plusieurs schemas PPM pour differents traceurs |
|---|
| 1105 | c IF (IORD.LE.0) then |
|---|
| 1106 | c if(IMR.GE.144) then |
|---|
| 1107 | c LMT = 0 |
|---|
| 1108 | c elseif(IMR.GE.72) then |
|---|
| 1109 | c LMT = 1 |
|---|
| 1110 | c else |
|---|
| 1111 | c LMT = 2 |
|---|
| 1112 | c endif |
|---|
| 1113 | c else |
|---|
| 1114 | c LMT = IORD - 3 |
|---|
| 1115 | c endif |
|---|
| 1116 | C |
|---|
| 1117 | LMT = IORD - 3 |
|---|
| 1118 | c write(6,*) 'PPM option in E-W direction = ', LMT |
|---|
| 1119 | c first = .false. |
|---|
| 1120 | C endif |
|---|
| 1121 | C |
|---|
| 1122 | DO i=1,IMR |
|---|
| 1123 | AL(i) = 0.5*(p(i-1)+p(i)) + (DC(i-1) - DC(i))*R3 |
|---|
| 1124 | END DO |
|---|
| 1125 | C |
|---|
| 1126 | do i=1,IMR-1 |
|---|
| 1127 | AR(i) = AL(i+1) |
|---|
| 1128 | END DO |
|---|
| 1129 | AR(IMR) = AL(1) |
|---|
| 1130 | C |
|---|
| 1131 | do i=1,IMR |
|---|
| 1132 | A6(i) = 3.*(p(i)+p(i) - (AL(i)+AR(i))) |
|---|
| 1133 | END DO |
|---|
| 1134 | C |
|---|
| 1135 | if(LMT<=2) call lmtppm(DC(1),A6(1),AR(1),AL(1),P(1),IMR,LMT) |
|---|
| 1136 | C |
|---|
| 1137 | AL(0) = AL(IMR) |
|---|
| 1138 | AR(0) = AR(IMR) |
|---|
| 1139 | A6(0) = A6(IMR) |
|---|
| 1140 | C |
|---|
| 1141 | DO i=1,IMR |
|---|
| 1142 | IF(UT(i)>0.) then |
|---|
| 1143 | flux(i) = AR(i-1) + 0.5*UT(i)*(AL(i-1) - AR(i-1) + |
|---|
| 1144 | & A6(i-1)*(1.-R23*UT(i)) ) |
|---|
| 1145 | else |
|---|
| 1146 | flux(i) = AL(i) - 0.5*UT(i)*(AR(i) - AL(i) + |
|---|
| 1147 | & A6(i)*(1.+R23*UT(i))) |
|---|
| 1148 | endif |
|---|
| 1149 | enddo |
|---|
| 1150 | return |
|---|
| 1151 | end |
|---|
| 1152 | C |
|---|
| 1153 | subroutine xmist(IMR,IML,P,DC) |
|---|
| 1154 | implicit none |
|---|
| 1155 | integer IMR,IML |
|---|
| 1156 | real,parameter :: R24 = 1./24. |
|---|
| 1157 | real :: P(-IML:IMR+1+IML),DC(-IML:IMR+1+IML) |
|---|
| 1158 | integer :: i |
|---|
| 1159 | real :: tmp,pmax,pmin |
|---|
| 1160 | C |
|---|
| 1161 | do i=1,IMR |
|---|
| 1162 | tmp = R24*(8.*(p(i+1) - p(i-1)) + p(i-2) - p(i+2)) |
|---|
| 1163 | Pmax = max(P(i-1), p(i), p(i+1)) - p(i) |
|---|
| 1164 | Pmin = p(i) - min(P(i-1), p(i), p(i+1)) |
|---|
| 1165 | DC(i) = sign(min(abs(tmp),Pmax,Pmin), tmp) |
|---|
| 1166 | END DO |
|---|
| 1167 | return |
|---|
| 1168 | end |
|---|
| 1169 | C |
|---|
| 1170 | subroutine ytp(IMR,JNP,j1,j2,acosp,RCAP,DQ,P,VC,DC2 |
|---|
| 1171 | & ,ymass,fx,A6,AR,AL,JORD) |
|---|
| 1172 | implicit none |
|---|
| 1173 | integer :: IMR,JNP,j1,j2,JORD |
|---|
| 1174 | real :: acosp,RCAP,DQ,P,VC,DC2,ymass,fx,A6,AR,AL |
|---|
| 1175 | dimension P(IMR,JNP),VC(IMR,JNP),ymass(IMR,JNP) |
|---|
| 1176 | & ,DC2(IMR,JNP),DQ(IMR,JNP),acosp(JNP) |
|---|
| 1177 | C Work array |
|---|
| 1178 | DIMENSION fx(IMR,JNP),AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP) |
|---|
| 1179 | integer :: JMR,len,i,jt,j |
|---|
| 1180 | real :: sum1,sum2 |
|---|
| 1181 | C |
|---|
| 1182 | JMR = JNP - 1 |
|---|
| 1183 | len = IMR*(J2-J1+2) |
|---|
| 1184 | C |
|---|
| 1185 | if(JORD==1) then |
|---|
| 1186 | DO i=1,len |
|---|
| 1187 | JT = REAL(J1) - VC(i,J1) |
|---|
| 1188 | fx(i,j1) = p(i,JT) |
|---|
| 1189 | END DO |
|---|
| 1190 | else |
|---|
| 1191 | |
|---|
| 1192 | call ymist(IMR,JNP,j1,P,DC2,4) |
|---|
| 1193 | C |
|---|
| 1194 | if(JORD<=0 .or. JORD>=3) then |
|---|
| 1195 | |
|---|
| 1196 | call fyppm(VC,P,DC2,fx,IMR,JNP,j1,j2,A6,AR,AL,JORD) |
|---|
| 1197 | |
|---|
| 1198 | else |
|---|
| 1199 | DO i=1,len |
|---|
| 1200 | JT = REAL(J1) - VC(i,J1) |
|---|
| 1201 | fx(i,j1) = p(i,JT) + (sign(1.,VC(i,j1))-VC(i,j1))*DC2(i,JT) |
|---|
| 1202 | END DO |
|---|
| 1203 | endif |
|---|
| 1204 | endif |
|---|
| 1205 | C |
|---|
| 1206 | DO i=1,len |
|---|
| 1207 | fx(i,j1) = fx(i,j1)*ymass(i,j1) |
|---|
| 1208 | END DO |
|---|
| 1209 | C |
|---|
| 1210 | DO j=j1,j2 |
|---|
| 1211 | DO i=1,IMR |
|---|
| 1212 | DQ(i,j) = DQ(i,j) + (fx(i,j) - fx(i,j+1)) * acosp(j) |
|---|
| 1213 | END DO |
|---|
| 1214 | END DO |
|---|
| 1215 | C |
|---|
| 1216 | C Poles |
|---|
| 1217 | sum1 = fx(IMR,j1 ) |
|---|
| 1218 | sum2 = fx(IMR,J2+1) |
|---|
| 1219 | do i=1,IMR-1 |
|---|
| 1220 | sum1 = sum1 + fx(i,j1 ) |
|---|
| 1221 | sum2 = sum2 + fx(i,J2+1) |
|---|
| 1222 | enddo |
|---|
| 1223 | C |
|---|
| 1224 | sum1 = DQ(1, 1) - sum1 * RCAP |
|---|
| 1225 | sum2 = DQ(1,JNP) + sum2 * RCAP |
|---|
| 1226 | do i=1,IMR |
|---|
| 1227 | DQ(i, 1) = sum1 |
|---|
| 1228 | DQ(i,JNP) = sum2 |
|---|
| 1229 | enddo |
|---|
| 1230 | C |
|---|
| 1231 | if(j1/=2) then |
|---|
| 1232 | do i=1,IMR |
|---|
| 1233 | DQ(i, 2) = sum1 |
|---|
| 1234 | DQ(i,JMR) = sum2 |
|---|
| 1235 | enddo |
|---|
| 1236 | endif |
|---|
| 1237 | C |
|---|
| 1238 | return |
|---|
| 1239 | end |
|---|
| 1240 | C |
|---|
| 1241 | subroutine ymist(IMR,JNP,j1,P,DC,ID) |
|---|
| 1242 | implicit none |
|---|
| 1243 | integer :: IMR,JNP,j1,ID |
|---|
| 1244 | real,parameter :: R24 = 1./24. |
|---|
| 1245 | real :: P(IMR,JNP),DC(IMR,JNP) |
|---|
| 1246 | integer :: iimh,jmr,ijm3,imh,i |
|---|
| 1247 | real :: pmax,pmin,tmp |
|---|
| 1248 | C |
|---|
| 1249 | IMH = IMR / 2 |
|---|
| 1250 | JMR = JNP - 1 |
|---|
| 1251 | IJM3 = IMR*(JMR-3) |
|---|
| 1252 | C |
|---|
| 1253 | IF(ID==2) THEN |
|---|
| 1254 | do i=1,IMR*(JMR-1) |
|---|
| 1255 | tmp = 0.25*(p(i,3) - p(i,1)) |
|---|
| 1256 | Pmax = max(p(i,1),p(i,2),p(i,3)) - p(i,2) |
|---|
| 1257 | Pmin = p(i,2) - min(p(i,1),p(i,2),p(i,3)) |
|---|
| 1258 | DC(i,2) = sign(min(abs(tmp),Pmin,Pmax),tmp) |
|---|
| 1259 | END DO |
|---|
| 1260 | ELSE |
|---|
| 1261 | do i=1,IMH |
|---|
| 1262 | C J=2 |
|---|
| 1263 | tmp = (8.*(p(i,3) - p(i,1)) + p(i+IMH,2) - p(i,4))*R24 |
|---|
| 1264 | Pmax = max(p(i,1),p(i,2),p(i,3)) - p(i,2) |
|---|
| 1265 | Pmin = p(i,2) - min(p(i,1),p(i,2),p(i,3)) |
|---|
| 1266 | DC(i,2) = sign(min(abs(tmp),Pmin,Pmax),tmp) |
|---|
| 1267 | C J=JMR |
|---|
| 1268 | tmp=(8.*(p(i,JNP)-p(i,JMR-1))+p(i,JMR-2)-p(i+IMH,JMR))*R24 |
|---|
| 1269 | Pmax = max(p(i,JMR-1),p(i,JMR),p(i,JNP)) - p(i,JMR) |
|---|
| 1270 | Pmin = p(i,JMR) - min(p(i,JMR-1),p(i,JMR),p(i,JNP)) |
|---|
| 1271 | DC(i,JMR) = sign(min(abs(tmp),Pmin,Pmax),tmp) |
|---|
| 1272 | END DO |
|---|
| 1273 | do i=IMH+1,IMR |
|---|
| 1274 | C J=2 |
|---|
| 1275 | tmp = (8.*(p(i,3) - p(i,1)) + p(i-IMH,2) - p(i,4))*R24 |
|---|
| 1276 | Pmax = max(p(i,1),p(i,2),p(i,3)) - p(i,2) |
|---|
| 1277 | Pmin = p(i,2) - min(p(i,1),p(i,2),p(i,3)) |
|---|
| 1278 | DC(i,2) = sign(min(abs(tmp),Pmin,Pmax),tmp) |
|---|
| 1279 | C J=JMR |
|---|
| 1280 | tmp=(8.*(p(i,JNP)-p(i,JMR-1))+p(i,JMR-2)-p(i-IMH,JMR))*R24 |
|---|
| 1281 | Pmax = max(p(i,JMR-1),p(i,JMR),p(i,JNP)) - p(i,JMR) |
|---|
| 1282 | Pmin = p(i,JMR) - min(p(i,JMR-1),p(i,JMR),p(i,JNP)) |
|---|
| 1283 | DC(i,JMR) = sign(min(abs(tmp),Pmin,Pmax),tmp) |
|---|
| 1284 | END DO |
|---|
| 1285 | C |
|---|
| 1286 | do i=1,IJM3 |
|---|
| 1287 | tmp = (8.*(p(i,4) - p(i,2)) + p(i,1) - p(i,5))*R24 |
|---|
| 1288 | Pmax = max(p(i,2),p(i,3),p(i,4)) - p(i,3) |
|---|
| 1289 | Pmin = p(i,3) - min(p(i,2),p(i,3),p(i,4)) |
|---|
| 1290 | DC(i,3) = sign(min(abs(tmp),Pmin,Pmax),tmp) |
|---|
| 1291 | END DO |
|---|
| 1292 | ENDIF |
|---|
| 1293 | C |
|---|
| 1294 | if(j1/=2) then |
|---|
| 1295 | do i=1,IMR |
|---|
| 1296 | DC(i,1) = 0. |
|---|
| 1297 | DC(i,JNP) = 0. |
|---|
| 1298 | enddo |
|---|
| 1299 | else |
|---|
| 1300 | C Determine slopes in polar caps for scalars! |
|---|
| 1301 | C |
|---|
| 1302 | do i=1,IMH |
|---|
| 1303 | C South |
|---|
| 1304 | tmp = 0.25*(p(i,2) - p(i+imh,2)) |
|---|
| 1305 | Pmax = max(p(i,2),p(i,1), p(i+imh,2)) - p(i,1) |
|---|
| 1306 | Pmin = p(i,1) - min(p(i,2),p(i,1), p(i+imh,2)) |
|---|
| 1307 | DC(i,1)=sign(min(abs(tmp),Pmax,Pmin),tmp) |
|---|
| 1308 | C North. |
|---|
| 1309 | tmp = 0.25*(p(i+imh,JMR) - p(i,JMR)) |
|---|
| 1310 | Pmax = max(p(i+imh,JMR),p(i,jnp), p(i,JMR)) - p(i,JNP) |
|---|
| 1311 | Pmin = p(i,JNP) - min(p(i+imh,JMR),p(i,jnp), p(i,JMR)) |
|---|
| 1312 | DC(i,JNP) = sign(min(abs(tmp),Pmax,pmin),tmp) |
|---|
| 1313 | END DO |
|---|
| 1314 | C |
|---|
| 1315 | do i=imh+1,IMR |
|---|
| 1316 | DC(i, 1) = - DC(i-imh, 1) |
|---|
| 1317 | DC(i,JNP) = - DC(i-imh,JNP) |
|---|
| 1318 | END DO |
|---|
| 1319 | endif |
|---|
| 1320 | return |
|---|
| 1321 | end |
|---|
| 1322 | C |
|---|
| 1323 | subroutine fyppm(VC,P,DC,flux,IMR,JNP,j1,j2,A6,AR,AL,JORD) |
|---|
| 1324 | implicit none |
|---|
| 1325 | integer IMR,JNP,j1,j2,JORD |
|---|
| 1326 | real,parameter :: R3 = 1./3., R23 = 2./3. |
|---|
| 1327 | real VC(IMR,*),flux(IMR,*),P(IMR,*),DC(IMR,*) |
|---|
| 1328 | C Local work arrays. |
|---|
| 1329 | real AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP) |
|---|
| 1330 | integer LMT,i |
|---|
| 1331 | integer IMH,JMR,j11,IMJM1,len |
|---|
| 1332 | c logical first |
|---|
| 1333 | C data first /.true./ |
|---|
| 1334 | C SAVE LMT |
|---|
| 1335 | C |
|---|
| 1336 | IMH = IMR / 2 |
|---|
| 1337 | JMR = JNP - 1 |
|---|
| 1338 | j11 = j1-1 |
|---|
| 1339 | IMJM1 = IMR*(J2-J1+2) |
|---|
| 1340 | len = IMR*(J2-J1+3) |
|---|
| 1341 | C if(first) then |
|---|
| 1342 | C IF(JORD.LE.0) then |
|---|
| 1343 | C if(JMR.GE.90) then |
|---|
| 1344 | C LMT = 0 |
|---|
| 1345 | C elseif(JMR.GE.45) then |
|---|
| 1346 | C LMT = 1 |
|---|
| 1347 | C else |
|---|
| 1348 | C LMT = 2 |
|---|
| 1349 | C endif |
|---|
| 1350 | C else |
|---|
| 1351 | C LMT = JORD - 3 |
|---|
| 1352 | C endif |
|---|
| 1353 | C |
|---|
| 1354 | C first = .false. |
|---|
| 1355 | C endif |
|---|
| 1356 | C |
|---|
| 1357 | c modifs pour pouvoir choisir plusieurs schemas PPM |
|---|
| 1358 | LMT = JORD - 3 |
|---|
| 1359 | C |
|---|
| 1360 | DO i=1,IMR*JMR |
|---|
| 1361 | AL(i,2) = 0.5*(p(i,1)+p(i,2)) + (DC(i,1) - DC(i,2))*R3 |
|---|
| 1362 | AR(i,1) = AL(i,2) |
|---|
| 1363 | END DO |
|---|
| 1364 | C |
|---|
| 1365 | CPoles: |
|---|
| 1366 | C |
|---|
| 1367 | DO i=1,IMH |
|---|
| 1368 | AL(i,1) = AL(i+IMH,2) |
|---|
| 1369 | AL(i+IMH,1) = AL(i,2) |
|---|
| 1370 | C |
|---|
| 1371 | AR(i,JNP) = AR(i+IMH,JMR) |
|---|
| 1372 | AR(i+IMH,JNP) = AR(i,JMR) |
|---|
| 1373 | ENDDO |
|---|
| 1374 | |
|---|
| 1375 | ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1376 | c Rajout pour LMDZ.3.3 |
|---|
| 1377 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1378 | AR(IMR,1)=AL(1,1) |
|---|
| 1379 | AR(IMR,JNP)=AL(1,JNP) |
|---|
| 1380 | ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1381 | |
|---|
| 1382 | |
|---|
| 1383 | do i=1,len |
|---|
| 1384 | A6(i,j11) = 3.*(p(i,j11)+p(i,j11) - (AL(i,j11)+AR(i,j11))) |
|---|
| 1385 | END DO |
|---|
| 1386 | C |
|---|
| 1387 | if(LMT<=2) call lmtppm(DC(1,j11),A6(1,j11),AR(1,j11) |
|---|
| 1388 | & ,AL(1,j11),P(1,j11),len,LMT) |
|---|
| 1389 | C |
|---|
| 1390 | |
|---|
| 1391 | DO i=1,IMJM1 |
|---|
| 1392 | IF(VC(i,j1)>0.) then |
|---|
| 1393 | flux(i,j1) = AR(i,j11) + 0.5*VC(i,j1)*(AL(i,j11) - AR(i,j11) + |
|---|
| 1394 | & A6(i,j11)*(1.-R23*VC(i,j1)) ) |
|---|
| 1395 | else |
|---|
| 1396 | flux(i,j1) = AL(i,j1) - 0.5*VC(i,j1)*(AR(i,j1) - AL(i,j1) + |
|---|
| 1397 | & A6(i,j1)*(1.+R23*VC(i,j1))) |
|---|
| 1398 | endif |
|---|
| 1399 | END DO |
|---|
| 1400 | return |
|---|
| 1401 | end |
|---|
| 1402 | C |
|---|
| 1403 | subroutine yadv(IMR,JNP,j1,j2,p,VA,ady,wk,IAD) |
|---|
| 1404 | implicit none |
|---|
| 1405 | integer IMR,JNP,j1,j2,IAD |
|---|
| 1406 | REAL p(IMR,JNP),ady(IMR,JNP),VA(IMR,JNP) |
|---|
| 1407 | REAL WK(IMR,-1:JNP+2) |
|---|
| 1408 | INTEGER JMR,IMH,i,j,jp |
|---|
| 1409 | REAL rv,a1,b1,sum1,sum2 |
|---|
| 1410 | C |
|---|
| 1411 | JMR = JNP-1 |
|---|
| 1412 | IMH = IMR/2 |
|---|
| 1413 | do j=1,JNP |
|---|
| 1414 | do i=1,IMR |
|---|
| 1415 | wk(i,j) = p(i,j) |
|---|
| 1416 | enddo |
|---|
| 1417 | enddo |
|---|
| 1418 | C Poles: |
|---|
| 1419 | do i=1,IMH |
|---|
| 1420 | wk(i, -1) = p(i+IMH,3) |
|---|
| 1421 | wk(i+IMH,-1) = p(i,3) |
|---|
| 1422 | wk(i, 0) = p(i+IMH,2) |
|---|
| 1423 | wk(i+IMH,0) = p(i,2) |
|---|
| 1424 | wk(i,JNP+1) = p(i+IMH,JMR) |
|---|
| 1425 | wk(i+IMH,JNP+1) = p(i,JMR) |
|---|
| 1426 | wk(i,JNP+2) = p(i+IMH,JNP-2) |
|---|
| 1427 | wk(i+IMH,JNP+2) = p(i,JNP-2) |
|---|
| 1428 | enddo |
|---|
| 1429 | c write(*,*) 'toto 1' |
|---|
| 1430 | C -------------------------------- |
|---|
| 1431 | IF(IAD==2) then |
|---|
| 1432 | do j=j1-1,j2+1 |
|---|
| 1433 | do i=1,IMR |
|---|
| 1434 | c write(*,*) 'avt NINT','i=',i,'j=',j |
|---|
| 1435 | JP = NINT(VA(i,j)) |
|---|
| 1436 | rv = JP - VA(i,j) |
|---|
| 1437 | c write(*,*) 'VA=',VA(i,j), 'JP1=',JP,'rv=',rv |
|---|
| 1438 | JP = j - JP |
|---|
| 1439 | c write(*,*) 'JP2=',JP |
|---|
| 1440 | a1 = 0.5*(wk(i,jp+1)+wk(i,jp-1)) - wk(i,jp) |
|---|
| 1441 | b1 = 0.5*(wk(i,jp+1)-wk(i,jp-1)) |
|---|
| 1442 | c write(*,*) 'a1=',a1,'b1=',b1 |
|---|
| 1443 | ady(i,j) = wk(i,jp) + rv*(a1*rv + b1) - wk(i,j) |
|---|
| 1444 | enddo |
|---|
| 1445 | enddo |
|---|
| 1446 | c write(*,*) 'toto 2' |
|---|
| 1447 | C |
|---|
| 1448 | ELSEIF(IAD==1) then |
|---|
| 1449 | do j=j1-1,j2+1 |
|---|
| 1450 | do i=1,imr |
|---|
| 1451 | JP = REAL(j)-VA(i,j) |
|---|
| 1452 | ady(i,j) = VA(i,j)*(wk(i,jp)-wk(i,jp+1)) |
|---|
| 1453 | enddo |
|---|
| 1454 | enddo |
|---|
| 1455 | ENDIF |
|---|
| 1456 | C |
|---|
| 1457 | if(j1/=2) then |
|---|
| 1458 | sum1 = 0. |
|---|
| 1459 | sum2 = 0. |
|---|
| 1460 | do i=1,imr |
|---|
| 1461 | sum1 = sum1 + ady(i,2) |
|---|
| 1462 | sum2 = sum2 + ady(i,JMR) |
|---|
| 1463 | enddo |
|---|
| 1464 | sum1 = sum1 / IMR |
|---|
| 1465 | sum2 = sum2 / IMR |
|---|
| 1466 | C |
|---|
| 1467 | do i=1,imr |
|---|
| 1468 | ady(i, 2) = sum1 |
|---|
| 1469 | ady(i,JMR) = sum2 |
|---|
| 1470 | ady(i, 1) = sum1 |
|---|
| 1471 | ady(i,JNP) = sum2 |
|---|
| 1472 | enddo |
|---|
| 1473 | else |
|---|
| 1474 | C Poles: |
|---|
| 1475 | sum1 = 0. |
|---|
| 1476 | sum2 = 0. |
|---|
| 1477 | do i=1,imr |
|---|
| 1478 | sum1 = sum1 + ady(i,1) |
|---|
| 1479 | sum2 = sum2 + ady(i,JNP) |
|---|
| 1480 | enddo |
|---|
| 1481 | sum1 = sum1 / IMR |
|---|
| 1482 | sum2 = sum2 / IMR |
|---|
| 1483 | C |
|---|
| 1484 | do i=1,imr |
|---|
| 1485 | ady(i, 1) = sum1 |
|---|
| 1486 | ady(i,JNP) = sum2 |
|---|
| 1487 | enddo |
|---|
| 1488 | endif |
|---|
| 1489 | C |
|---|
| 1490 | return |
|---|
| 1491 | end |
|---|
| 1492 | C |
|---|
| 1493 | subroutine xadv(IMR,JNP,j1,j2,p,UA,JS,JN,IML,adx,IAD) |
|---|
| 1494 | implicit none |
|---|
| 1495 | INTEGER IMR,JNP,j1,j2,JS,JN,IML,IAD |
|---|
| 1496 | REAL p(IMR,JNP),adx(IMR,JNP),qtmp(-IMR:IMR+IMR),UA(IMR,JNP) |
|---|
| 1497 | INTEGER JMR,j,i,ip,iu,iiu |
|---|
| 1498 | REAL ru,a1,b1 |
|---|
| 1499 | C |
|---|
| 1500 | JMR = JNP-1 |
|---|
| 1501 | do j=j1,j2 |
|---|
| 1502 | if(J>JS .and. J<JN) GO TO 1309 |
|---|
| 1503 | C |
|---|
| 1504 | do i=1,IMR |
|---|
| 1505 | qtmp(i) = p(i,j) |
|---|
| 1506 | enddo |
|---|
| 1507 | C |
|---|
| 1508 | do i=-IML,0 |
|---|
| 1509 | qtmp(i) = p(IMR+i,j) |
|---|
| 1510 | qtmp(IMR+1-i) = p(1-i,j) |
|---|
| 1511 | enddo |
|---|
| 1512 | C |
|---|
| 1513 | IF(IAD==2) THEN |
|---|
| 1514 | DO i=1,IMR |
|---|
| 1515 | IP = NINT(UA(i,j)) |
|---|
| 1516 | ru = IP - UA(i,j) |
|---|
| 1517 | IP = i - IP |
|---|
| 1518 | a1 = 0.5*(qtmp(ip+1)+qtmp(ip-1)) - qtmp(ip) |
|---|
| 1519 | b1 = 0.5*(qtmp(ip+1)-qtmp(ip-1)) |
|---|
| 1520 | adx(i,j) = qtmp(ip) + ru*(a1*ru + b1) |
|---|
| 1521 | enddo |
|---|
| 1522 | ELSEIF(IAD==1) then |
|---|
| 1523 | DO i=1,IMR |
|---|
| 1524 | iu = UA(i,j) |
|---|
| 1525 | ru = UA(i,j) - iu |
|---|
| 1526 | iiu = i-iu |
|---|
| 1527 | if(UA(i,j)>=0.) then |
|---|
| 1528 | adx(i,j) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu)) |
|---|
| 1529 | else |
|---|
| 1530 | adx(i,j) = qtmp(iiu)+ru*(qtmp(iiu)-qtmp(iiu+1)) |
|---|
| 1531 | endif |
|---|
| 1532 | enddo |
|---|
| 1533 | ENDIF |
|---|
| 1534 | C |
|---|
| 1535 | do i=1,IMR |
|---|
| 1536 | adx(i,j) = adx(i,j) - p(i,j) |
|---|
| 1537 | enddo |
|---|
| 1538 | 1309 continue |
|---|
| 1539 | END DO |
|---|
| 1540 | C |
|---|
| 1541 | C Eulerian upwind |
|---|
| 1542 | C |
|---|
| 1543 | do j=JS+1,JN-1 |
|---|
| 1544 | C |
|---|
| 1545 | do i=1,IMR |
|---|
| 1546 | qtmp(i) = p(i,j) |
|---|
| 1547 | enddo |
|---|
| 1548 | C |
|---|
| 1549 | qtmp(0) = p(IMR,J) |
|---|
| 1550 | qtmp(IMR+1) = p(1,J) |
|---|
| 1551 | C |
|---|
| 1552 | IF(IAD==2) THEN |
|---|
| 1553 | qtmp(-1) = p(IMR-1,J) |
|---|
| 1554 | qtmp(IMR+2) = p(2,J) |
|---|
| 1555 | do i=1,imr |
|---|
| 1556 | IP = NINT(UA(i,j)) |
|---|
| 1557 | ru = IP - UA(i,j) |
|---|
| 1558 | IP = i - IP |
|---|
| 1559 | a1 = 0.5*(qtmp(ip+1)+qtmp(ip-1)) - qtmp(ip) |
|---|
| 1560 | b1 = 0.5*(qtmp(ip+1)-qtmp(ip-1)) |
|---|
| 1561 | adx(i,j) = qtmp(ip)- p(i,j) + ru*(a1*ru + b1) |
|---|
| 1562 | enddo |
|---|
| 1563 | ELSEIF(IAD==1) then |
|---|
| 1564 | C 1st order |
|---|
| 1565 | DO i=1,IMR |
|---|
| 1566 | IP = i - UA(i,j) |
|---|
| 1567 | adx(i,j) = UA(i,j)*(qtmp(ip)-qtmp(ip+1)) |
|---|
| 1568 | enddo |
|---|
| 1569 | ENDIF |
|---|
| 1570 | enddo |
|---|
| 1571 | C |
|---|
| 1572 | if(j1/=2) then |
|---|
| 1573 | do i=1,IMR |
|---|
| 1574 | adx(i, 2) = 0. |
|---|
| 1575 | adx(i,JMR) = 0. |
|---|
| 1576 | enddo |
|---|
| 1577 | endif |
|---|
| 1578 | C set cross term due to x-adv at the poles to zero. |
|---|
| 1579 | do i=1,IMR |
|---|
| 1580 | adx(i, 1) = 0. |
|---|
| 1581 | adx(i,JNP) = 0. |
|---|
| 1582 | enddo |
|---|
| 1583 | return |
|---|
| 1584 | end |
|---|
| 1585 | C |
|---|
| 1586 | subroutine lmtppm(DC,A6,AR,AL,P,IM,LMT) |
|---|
| 1587 | implicit none |
|---|
| 1588 | C |
|---|
| 1589 | C A6 = CURVATURE OF THE TEST PARABOLA |
|---|
| 1590 | C AR = RIGHT EDGE VALUE OF THE TEST PARABOLA |
|---|
| 1591 | C AL = LEFT EDGE VALUE OF THE TEST PARABOLA |
|---|
| 1592 | C DC = 0.5 * MISMATCH |
|---|
| 1593 | C P = CELL-AVERAGED VALUE |
|---|
| 1594 | C IM = VECTOR LENGTH |
|---|
| 1595 | C |
|---|
| 1596 | C OPTIONS: |
|---|
| 1597 | C |
|---|
| 1598 | C LMT = 0: FULL MONOTONICITY |
|---|
| 1599 | C LMT = 1: SEMI-MONOTONIC CONSTRAINT (NO UNDERSHOOTS) |
|---|
| 1600 | C LMT = 2: POSITIVE-DEFINITE CONSTRAINT |
|---|
| 1601 | C |
|---|
| 1602 | real,parameter :: R12 = 1./12. |
|---|
| 1603 | real :: A6(IM),AR(IM),AL(IM),P(IM),DC(IM) |
|---|
| 1604 | integer :: IM,LMT |
|---|
| 1605 | INTEGER i |
|---|
| 1606 | REAL da1,da2,a6da,fmin |
|---|
| 1607 | C |
|---|
| 1608 | if(LMT==0) then |
|---|
| 1609 | C Full constraint |
|---|
| 1610 | do i=1,IM |
|---|
| 1611 | if(DC(i)==0.) then |
|---|
| 1612 | AR(i) = p(i) |
|---|
| 1613 | AL(i) = p(i) |
|---|
| 1614 | A6(i) = 0. |
|---|
| 1615 | else |
|---|
| 1616 | da1 = AR(i) - AL(i) |
|---|
| 1617 | da2 = da1**2 |
|---|
| 1618 | A6DA = A6(i)*da1 |
|---|
| 1619 | if(A6DA < -da2) then |
|---|
| 1620 | A6(i) = 3.*(AL(i)-p(i)) |
|---|
| 1621 | AR(i) = AL(i) - A6(i) |
|---|
| 1622 | elseif(A6DA > da2) then |
|---|
| 1623 | A6(i) = 3.*(AR(i)-p(i)) |
|---|
| 1624 | AL(i) = AR(i) - A6(i) |
|---|
| 1625 | endif |
|---|
| 1626 | endif |
|---|
| 1627 | END DO |
|---|
| 1628 | elseif(LMT==1) then |
|---|
| 1629 | C Semi-monotonic constraint |
|---|
| 1630 | do i=1,IM |
|---|
| 1631 | if(abs(AR(i)-AL(i)) >= -A6(i)) go to 150 |
|---|
| 1632 | if(p(i)<AR(i) .and. p(i)<AL(i)) then |
|---|
| 1633 | AR(i) = p(i) |
|---|
| 1634 | AL(i) = p(i) |
|---|
| 1635 | A6(i) = 0. |
|---|
| 1636 | elseif(AR(i) > AL(i)) then |
|---|
| 1637 | A6(i) = 3.*(AL(i)-p(i)) |
|---|
| 1638 | AR(i) = AL(i) - A6(i) |
|---|
| 1639 | else |
|---|
| 1640 | A6(i) = 3.*(AR(i)-p(i)) |
|---|
| 1641 | AL(i) = AR(i) - A6(i) |
|---|
| 1642 | endif |
|---|
| 1643 | 150 continue |
|---|
| 1644 | END DO |
|---|
| 1645 | elseif(LMT==2) then |
|---|
| 1646 | do i=1,IM |
|---|
| 1647 | if(abs(AR(i)-AL(i)) >= -A6(i)) go to 250 |
|---|
| 1648 | fmin = p(i) + 0.25*(AR(i)-AL(i))**2/A6(i) + A6(i)*R12 |
|---|
| 1649 | if(fmin>=0.) go to 250 |
|---|
| 1650 | if(p(i)<AR(i) .and. p(i)<AL(i)) then |
|---|
| 1651 | AR(i) = p(i) |
|---|
| 1652 | AL(i) = p(i) |
|---|
| 1653 | A6(i) = 0. |
|---|
| 1654 | elseif(AR(i) > AL(i)) then |
|---|
| 1655 | A6(i) = 3.*(AL(i)-p(i)) |
|---|
| 1656 | AR(i) = AL(i) - A6(i) |
|---|
| 1657 | else |
|---|
| 1658 | A6(i) = 3.*(AR(i)-p(i)) |
|---|
| 1659 | AL(i) = AR(i) - A6(i) |
|---|
| 1660 | endif |
|---|
| 1661 | 250 continue |
|---|
| 1662 | END DO |
|---|
| 1663 | endif |
|---|
| 1664 | return |
|---|
| 1665 | end |
|---|
| 1666 | C |
|---|
| 1667 | subroutine A2C(U,V,IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5) |
|---|
| 1668 | implicit none |
|---|
| 1669 | integer IMR,JMR,j1,j2 |
|---|
| 1670 | real :: U(IMR,*),V(IMR,*),CRX(IMR,*),CRY(IMR,*),DTDX5(*),DTDY5 |
|---|
| 1671 | integer i,j |
|---|
| 1672 | C |
|---|
| 1673 | do j=j1,j2 |
|---|
| 1674 | do i=2,IMR |
|---|
| 1675 | CRX(i,J) = dtdx5(j)*(U(i,j)+U(i-1,j)) |
|---|
| 1676 | END DO |
|---|
| 1677 | END DO |
|---|
| 1678 | C |
|---|
| 1679 | do j=j1,j2 |
|---|
| 1680 | CRX(1,J) = dtdx5(j)*(U(1,j)+U(IMR,j)) |
|---|
| 1681 | END DO |
|---|
| 1682 | C |
|---|
| 1683 | do i=1,IMR*JMR |
|---|
| 1684 | CRY(i,2) = DTDY5*(V(i,2)+V(i,1)) |
|---|
| 1685 | END DO |
|---|
| 1686 | return |
|---|
| 1687 | end |
|---|
| 1688 | C |
|---|
| 1689 | subroutine cosa(cosp,cose,JNP,PI,DP) |
|---|
| 1690 | implicit none |
|---|
| 1691 | integer JNP |
|---|
| 1692 | real :: cosp(*),cose(*),PI,DP |
|---|
| 1693 | integer JMR,j,jeq |
|---|
| 1694 | real ph5 |
|---|
| 1695 | JMR = JNP-1 |
|---|
| 1696 | do j=2,JNP |
|---|
| 1697 | ph5 = -0.5*PI + (REAL(J-1)-0.5)*DP |
|---|
| 1698 | cose(j) = cos(ph5) |
|---|
| 1699 | END DO |
|---|
| 1700 | C |
|---|
| 1701 | JEQ = (JNP+1) / 2 |
|---|
| 1702 | if(JMR == 2*(JMR/2) ) then |
|---|
| 1703 | do j=JNP, JEQ+1, -1 |
|---|
| 1704 | cose(j) = cose(JNP+2-j) |
|---|
| 1705 | enddo |
|---|
| 1706 | else |
|---|
| 1707 | C cell edge at equator. |
|---|
| 1708 | cose(JEQ+1) = 1. |
|---|
| 1709 | do j=JNP, JEQ+2, -1 |
|---|
| 1710 | cose(j) = cose(JNP+2-j) |
|---|
| 1711 | enddo |
|---|
| 1712 | endif |
|---|
| 1713 | C |
|---|
| 1714 | do j=2,JMR |
|---|
| 1715 | cosp(j) = 0.5*(cose(j)+cose(j+1)) |
|---|
| 1716 | END DO |
|---|
| 1717 | cosp(1) = 0. |
|---|
| 1718 | cosp(JNP) = 0. |
|---|
| 1719 | return |
|---|
| 1720 | end |
|---|
| 1721 | C |
|---|
| 1722 | subroutine cosc(cosp,cose,JNP,PI,DP) |
|---|
| 1723 | implicit none |
|---|
| 1724 | integer JNP |
|---|
| 1725 | real :: cosp(*),cose(*),PI,DP |
|---|
| 1726 | real phi |
|---|
| 1727 | integer j |
|---|
| 1728 | C |
|---|
| 1729 | phi = -0.5*PI |
|---|
| 1730 | do j=2,JNP-1 |
|---|
| 1731 | phi = phi + DP |
|---|
| 1732 | cosp(j) = cos(phi) |
|---|
| 1733 | END DO |
|---|
| 1734 | cosp( 1) = 0. |
|---|
| 1735 | cosp(JNP) = 0. |
|---|
| 1736 | C |
|---|
| 1737 | do j=2,JNP |
|---|
| 1738 | cose(j) = 0.5*(cosp(j)+cosp(j-1)) |
|---|
| 1739 | END DO |
|---|
| 1740 | C |
|---|
| 1741 | do j=2,JNP-1 |
|---|
| 1742 | cosp(j) = 0.5*(cose(j)+cose(j+1)) |
|---|
| 1743 | END DO |
|---|
| 1744 | return |
|---|
| 1745 | end |
|---|
| 1746 | C |
|---|
| 1747 | SUBROUTINE qckxyz (Q,qtmp,IMR,JNP,NLAY,j1,j2,cosp,acosp, |
|---|
| 1748 | & cross,IC,NSTEP) |
|---|
| 1749 | C |
|---|
| 1750 | real,parameter :: tiny = 1.E-60 |
|---|
| 1751 | INTEGER :: IMR,JNP,NLAY,j1,j2,IC,NSTEP |
|---|
| 1752 | REAL :: Q(IMR,JNP,NLAY),qtmp(IMR,JNP),cosp(*),acosp(*) |
|---|
| 1753 | logical cross |
|---|
| 1754 | INTEGER :: NLAYM1,len,ip,L,icr,ipy,ipx,i |
|---|
| 1755 | real :: qup,qly,dup,sum |
|---|
| 1756 | C |
|---|
| 1757 | NLAYM1 = NLAY-1 |
|---|
| 1758 | len = IMR*(j2-j1+1) |
|---|
| 1759 | ip = 0 |
|---|
| 1760 | C |
|---|
| 1761 | C Top layer |
|---|
| 1762 | L = 1 |
|---|
| 1763 | icr = 1 |
|---|
| 1764 | call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny) |
|---|
| 1765 | if(ipy==0) goto 50 |
|---|
| 1766 | call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny) |
|---|
| 1767 | if(ipx==0) goto 50 |
|---|
| 1768 | C |
|---|
| 1769 | if(cross) then |
|---|
| 1770 | call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny) |
|---|
| 1771 | endif |
|---|
| 1772 | if(icr==0) goto 50 |
|---|
| 1773 | C |
|---|
| 1774 | C Vertical filling... |
|---|
| 1775 | do i=1,len |
|---|
| 1776 | IF( Q(i,j1,1)<0.) THEN |
|---|
| 1777 | ip = ip + 1 |
|---|
| 1778 | Q(i,j1,2) = Q(i,j1,2) + Q(i,j1,1) |
|---|
| 1779 | Q(i,j1,1) = 0. |
|---|
| 1780 | endif |
|---|
| 1781 | enddo |
|---|
| 1782 | C |
|---|
| 1783 | 50 continue |
|---|
| 1784 | DO L = 2,NLAYM1 |
|---|
| 1785 | icr = 1 |
|---|
| 1786 | C |
|---|
| 1787 | call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny) |
|---|
| 1788 | if(ipy==0) goto 225 |
|---|
| 1789 | call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny) |
|---|
| 1790 | if(ipx==0) go to 225 |
|---|
| 1791 | if(cross) then |
|---|
| 1792 | call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny) |
|---|
| 1793 | endif |
|---|
| 1794 | if(icr==0) goto 225 |
|---|
| 1795 | C |
|---|
| 1796 | do i=1,len |
|---|
| 1797 | IF( Q(I,j1,L)<0.) THEN |
|---|
| 1798 | C |
|---|
| 1799 | ip = ip + 1 |
|---|
| 1800 | C From above |
|---|
| 1801 | qup = Q(I,j1,L-1) |
|---|
| 1802 | qly = -Q(I,j1,L) |
|---|
| 1803 | dup = min(qly,qup) |
|---|
| 1804 | Q(I,j1,L-1) = qup - dup |
|---|
| 1805 | Q(I,j1,L ) = dup-qly |
|---|
| 1806 | C Below |
|---|
| 1807 | Q(I,j1,L+1) = Q(I,j1,L+1) + Q(I,j1,L) |
|---|
| 1808 | Q(I,j1,L) = 0. |
|---|
| 1809 | ENDIF |
|---|
| 1810 | ENDDO |
|---|
| 1811 | 225 CONTINUE |
|---|
| 1812 | END DO |
|---|
| 1813 | C |
|---|
| 1814 | C BOTTOM LAYER |
|---|
| 1815 | sum = 0. |
|---|
| 1816 | L = NLAY |
|---|
| 1817 | C |
|---|
| 1818 | call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny) |
|---|
| 1819 | if(ipy==0) goto 911 |
|---|
| 1820 | call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny) |
|---|
| 1821 | if(ipx==0) goto 911 |
|---|
| 1822 | C |
|---|
| 1823 | call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny) |
|---|
| 1824 | if(icr==0) goto 911 |
|---|
| 1825 | C |
|---|
| 1826 | DO I=1,len |
|---|
| 1827 | IF( Q(I,j1,L)<0.) THEN |
|---|
| 1828 | ip = ip + 1 |
|---|
| 1829 | c |
|---|
| 1830 | C From above |
|---|
| 1831 | C |
|---|
| 1832 | qup = Q(I,j1,NLAYM1) |
|---|
| 1833 | qly = -Q(I,j1,L) |
|---|
| 1834 | dup = min(qly,qup) |
|---|
| 1835 | Q(I,j1,NLAYM1) = qup - dup |
|---|
| 1836 | C From "below" the surface. |
|---|
| 1837 | sum = sum + qly-dup |
|---|
| 1838 | Q(I,j1,L) = 0. |
|---|
| 1839 | ENDIF |
|---|
| 1840 | ENDDO |
|---|
| 1841 | C |
|---|
| 1842 | 911 continue |
|---|
| 1843 | C |
|---|
| 1844 | if(ip>IMR) then |
|---|
| 1845 | write(6,*) 'IC=',IC,' STEP=',NSTEP, |
|---|
| 1846 | & ' Vertical filling pts=',ip |
|---|
| 1847 | endif |
|---|
| 1848 | C |
|---|
| 1849 | if(sum>1.e-25) then |
|---|
| 1850 | write(6,*) IC,NSTEP,' Mass source from the ground=',sum |
|---|
| 1851 | endif |
|---|
| 1852 | RETURN |
|---|
| 1853 | END |
|---|
| 1854 | C |
|---|
| 1855 | subroutine filcr(q,IMR,JNP,j1,j2,cosp,acosp,icr,tiny) |
|---|
| 1856 | implicit none |
|---|
| 1857 | integer :: IMR,JNP,j1,j2,icr |
|---|
| 1858 | real :: q(IMR,*),cosp(*),acosp(*),tiny |
|---|
| 1859 | integer :: i,j |
|---|
| 1860 | real :: dq,dn,d0,d1,ds,d2 |
|---|
| 1861 | icr = 0 |
|---|
| 1862 | do j=j1+1,j2-1 |
|---|
| 1863 | DO i=1,IMR-1 |
|---|
| 1864 | IF(q(i,j)<0.) THEN |
|---|
| 1865 | icr = 1 |
|---|
| 1866 | dq = - q(i,j)*cosp(j) |
|---|
| 1867 | C N-E |
|---|
| 1868 | dn = q(i+1,j+1)*cosp(j+1) |
|---|
| 1869 | d0 = max(0.,dn) |
|---|
| 1870 | d1 = min(dq,d0) |
|---|
| 1871 | q(i+1,j+1) = (dn - d1)*acosp(j+1) |
|---|
| 1872 | dq = dq - d1 |
|---|
| 1873 | C S-E |
|---|
| 1874 | ds = q(i+1,j-1)*cosp(j-1) |
|---|
| 1875 | d0 = max(0.,ds) |
|---|
| 1876 | d2 = min(dq,d0) |
|---|
| 1877 | q(i+1,j-1) = (ds - d2)*acosp(j-1) |
|---|
| 1878 | q(i,j) = (d2 - dq)*acosp(j) + tiny |
|---|
| 1879 | endif |
|---|
| 1880 | END DO |
|---|
| 1881 | if(icr==0 .and. q(IMR,j)>=0.) goto 65 |
|---|
| 1882 | DO i=2,IMR |
|---|
| 1883 | IF(q(i,j)<0.) THEN |
|---|
| 1884 | icr = 1 |
|---|
| 1885 | dq = - q(i,j)*cosp(j) |
|---|
| 1886 | C N-W |
|---|
| 1887 | dn = q(i-1,j+1)*cosp(j+1) |
|---|
| 1888 | d0 = max(0.,dn) |
|---|
| 1889 | d1 = min(dq,d0) |
|---|
| 1890 | q(i-1,j+1) = (dn - d1)*acosp(j+1) |
|---|
| 1891 | dq = dq - d1 |
|---|
| 1892 | C S-W |
|---|
| 1893 | ds = q(i-1,j-1)*cosp(j-1) |
|---|
| 1894 | d0 = max(0.,ds) |
|---|
| 1895 | d2 = min(dq,d0) |
|---|
| 1896 | q(i-1,j-1) = (ds - d2)*acosp(j-1) |
|---|
| 1897 | q(i,j) = (d2 - dq)*acosp(j) + tiny |
|---|
| 1898 | endif |
|---|
| 1899 | END DO |
|---|
| 1900 | C ***************************************** |
|---|
| 1901 | C i=1 |
|---|
| 1902 | i=1 |
|---|
| 1903 | IF(q(i,j)<0.) THEN |
|---|
| 1904 | icr = 1 |
|---|
| 1905 | dq = - q(i,j)*cosp(j) |
|---|
| 1906 | C N-W |
|---|
| 1907 | dn = q(IMR,j+1)*cosp(j+1) |
|---|
| 1908 | d0 = max(0.,dn) |
|---|
| 1909 | d1 = min(dq,d0) |
|---|
| 1910 | q(IMR,j+1) = (dn - d1)*acosp(j+1) |
|---|
| 1911 | dq = dq - d1 |
|---|
| 1912 | C S-W |
|---|
| 1913 | ds = q(IMR,j-1)*cosp(j-1) |
|---|
| 1914 | d0 = max(0.,ds) |
|---|
| 1915 | d2 = min(dq,d0) |
|---|
| 1916 | q(IMR,j-1) = (ds - d2)*acosp(j-1) |
|---|
| 1917 | q(i,j) = (d2 - dq)*acosp(j) + tiny |
|---|
| 1918 | endif |
|---|
| 1919 | C ***************************************** |
|---|
| 1920 | C i=IMR |
|---|
| 1921 | i=IMR |
|---|
| 1922 | IF(q(i,j)<0.) THEN |
|---|
| 1923 | icr = 1 |
|---|
| 1924 | dq = - q(i,j)*cosp(j) |
|---|
| 1925 | C N-E |
|---|
| 1926 | dn = q(1,j+1)*cosp(j+1) |
|---|
| 1927 | d0 = max(0.,dn) |
|---|
| 1928 | d1 = min(dq,d0) |
|---|
| 1929 | q(1,j+1) = (dn - d1)*acosp(j+1) |
|---|
| 1930 | dq = dq - d1 |
|---|
| 1931 | C S-E |
|---|
| 1932 | ds = q(1,j-1)*cosp(j-1) |
|---|
| 1933 | d0 = max(0.,ds) |
|---|
| 1934 | d2 = min(dq,d0) |
|---|
| 1935 | q(1,j-1) = (ds - d2)*acosp(j-1) |
|---|
| 1936 | q(i,j) = (d2 - dq)*acosp(j) + tiny |
|---|
| 1937 | endif |
|---|
| 1938 | C ***************************************** |
|---|
| 1939 | 65 continue |
|---|
| 1940 | END DO |
|---|
| 1941 | C |
|---|
| 1942 | do i=1,IMR |
|---|
| 1943 | if(q(i,j1)<0. .or. q(i,j2)<0.) then |
|---|
| 1944 | icr = 1 |
|---|
| 1945 | goto 80 |
|---|
| 1946 | endif |
|---|
| 1947 | enddo |
|---|
| 1948 | C |
|---|
| 1949 | 80 continue |
|---|
| 1950 | C |
|---|
| 1951 | if(q(1,1)<0. .or. q(1,jnp)<0.) then |
|---|
| 1952 | icr = 1 |
|---|
| 1953 | endif |
|---|
| 1954 | C |
|---|
| 1955 | return |
|---|
| 1956 | end |
|---|
| 1957 | C |
|---|
| 1958 | subroutine filns(q,IMR,JNP,j1,j2,cosp,acosp,ipy,tiny) |
|---|
| 1959 | implicit none |
|---|
| 1960 | integer :: IMR,JNP,j1,j2,ipy |
|---|
| 1961 | real :: q(IMR,*),cosp(*),acosp(*),tiny |
|---|
| 1962 | real :: DP,CAP1,dq,dn,d0,d1,ds,d2 |
|---|
| 1963 | INTEGER :: i,j |
|---|
| 1964 | c logical first |
|---|
| 1965 | c data first /.true./ |
|---|
| 1966 | c save cap1 |
|---|
| 1967 | C |
|---|
| 1968 | c if(first) then |
|---|
| 1969 | DP = 4.*ATAN(1.)/REAL(JNP-1) |
|---|
| 1970 | CAP1 = IMR*(1.-COS((j1-1.5)*DP))/DP |
|---|
| 1971 | c first = .false. |
|---|
| 1972 | c endif |
|---|
| 1973 | C |
|---|
| 1974 | ipy = 0 |
|---|
| 1975 | do j=j1+1,j2-1 |
|---|
| 1976 | DO i=1,IMR |
|---|
| 1977 | IF(q(i,j)<0.) THEN |
|---|
| 1978 | ipy = 1 |
|---|
| 1979 | dq = - q(i,j)*cosp(j) |
|---|
| 1980 | C North |
|---|
| 1981 | dn = q(i,j+1)*cosp(j+1) |
|---|
| 1982 | d0 = max(0.,dn) |
|---|
| 1983 | d1 = min(dq,d0) |
|---|
| 1984 | q(i,j+1) = (dn - d1)*acosp(j+1) |
|---|
| 1985 | dq = dq - d1 |
|---|
| 1986 | C South |
|---|
| 1987 | ds = q(i,j-1)*cosp(j-1) |
|---|
| 1988 | d0 = max(0.,ds) |
|---|
| 1989 | d2 = min(dq,d0) |
|---|
| 1990 | q(i,j-1) = (ds - d2)*acosp(j-1) |
|---|
| 1991 | q(i,j) = (d2 - dq)*acosp(j) + tiny |
|---|
| 1992 | endif |
|---|
| 1993 | END DO |
|---|
| 1994 | END DO |
|---|
| 1995 | C |
|---|
| 1996 | do i=1,imr |
|---|
| 1997 | IF(q(i,j1)<0.) THEN |
|---|
| 1998 | ipy = 1 |
|---|
| 1999 | dq = - q(i,j1)*cosp(j1) |
|---|
| 2000 | C North |
|---|
| 2001 | dn = q(i,j1+1)*cosp(j1+1) |
|---|
| 2002 | d0 = max(0.,dn) |
|---|
| 2003 | d1 = min(dq,d0) |
|---|
| 2004 | q(i,j1+1) = (dn - d1)*acosp(j1+1) |
|---|
| 2005 | q(i,j1) = (d1 - dq)*acosp(j1) + tiny |
|---|
| 2006 | endif |
|---|
| 2007 | enddo |
|---|
| 2008 | C |
|---|
| 2009 | j = j2 |
|---|
| 2010 | do i=1,imr |
|---|
| 2011 | IF(q(i,j)<0.) THEN |
|---|
| 2012 | ipy = 1 |
|---|
| 2013 | dq = - q(i,j)*cosp(j) |
|---|
| 2014 | C South |
|---|
| 2015 | ds = q(i,j-1)*cosp(j-1) |
|---|
| 2016 | d0 = max(0.,ds) |
|---|
| 2017 | d2 = min(dq,d0) |
|---|
| 2018 | q(i,j-1) = (ds - d2)*acosp(j-1) |
|---|
| 2019 | q(i,j) = (d2 - dq)*acosp(j) + tiny |
|---|
| 2020 | endif |
|---|
| 2021 | enddo |
|---|
| 2022 | C |
|---|
| 2023 | C Check Poles. |
|---|
| 2024 | if(q(1,1)<0.) then |
|---|
| 2025 | dq = q(1,1)*cap1/REAL(IMR)*acosp(j1) |
|---|
| 2026 | do i=1,imr |
|---|
| 2027 | q(i,1) = 0. |
|---|
| 2028 | q(i,j1) = q(i,j1) + dq |
|---|
| 2029 | if(q(i,j1)<0.) ipy = 1 |
|---|
| 2030 | enddo |
|---|
| 2031 | endif |
|---|
| 2032 | C |
|---|
| 2033 | if(q(1,JNP)<0.) then |
|---|
| 2034 | dq = q(1,JNP)*cap1/REAL(IMR)*acosp(j2) |
|---|
| 2035 | do i=1,imr |
|---|
| 2036 | q(i,JNP) = 0. |
|---|
| 2037 | q(i,j2) = q(i,j2) + dq |
|---|
| 2038 | if(q(i,j2)<0.) ipy = 1 |
|---|
| 2039 | enddo |
|---|
| 2040 | endif |
|---|
| 2041 | C |
|---|
| 2042 | return |
|---|
| 2043 | end |
|---|
| 2044 | C |
|---|
| 2045 | subroutine filew(q,qtmp,IMR,JNP,j1,j2,ipx,tiny) |
|---|
| 2046 | implicit none |
|---|
| 2047 | integer :: IMR,JNP,j1,j2,ipx |
|---|
| 2048 | real :: q(IMR,*),qtmp(JNP,IMR),tiny |
|---|
| 2049 | integer :: i,j |
|---|
| 2050 | real :: d0,d1,d2 |
|---|
| 2051 | C |
|---|
| 2052 | ipx = 0 |
|---|
| 2053 | C Copy & swap direction for vectorization. |
|---|
| 2054 | do i=1,imr |
|---|
| 2055 | do j=j1,j2 |
|---|
| 2056 | qtmp(j,i) = q(i,j) |
|---|
| 2057 | END DO |
|---|
| 2058 | END DO |
|---|
| 2059 | C |
|---|
| 2060 | do i=2,imr-1 |
|---|
| 2061 | do j=j1,j2 |
|---|
| 2062 | if(qtmp(j,i)<0.) then |
|---|
| 2063 | ipx = 1 |
|---|
| 2064 | c west |
|---|
| 2065 | d0 = max(0.,qtmp(j,i-1)) |
|---|
| 2066 | d1 = min(-qtmp(j,i),d0) |
|---|
| 2067 | qtmp(j,i-1) = qtmp(j,i-1) - d1 |
|---|
| 2068 | qtmp(j,i) = qtmp(j,i) + d1 |
|---|
| 2069 | c east |
|---|
| 2070 | d0 = max(0.,qtmp(j,i+1)) |
|---|
| 2071 | d2 = min(-qtmp(j,i),d0) |
|---|
| 2072 | qtmp(j,i+1) = qtmp(j,i+1) - d2 |
|---|
| 2073 | qtmp(j,i) = qtmp(j,i) + d2 + tiny |
|---|
| 2074 | endif |
|---|
| 2075 | END DO |
|---|
| 2076 | END DO |
|---|
| 2077 | c |
|---|
| 2078 | i=1 |
|---|
| 2079 | do j=j1,j2 |
|---|
| 2080 | if(qtmp(j,i)<0.) then |
|---|
| 2081 | ipx = 1 |
|---|
| 2082 | c west |
|---|
| 2083 | d0 = max(0.,qtmp(j,imr)) |
|---|
| 2084 | d1 = min(-qtmp(j,i),d0) |
|---|
| 2085 | qtmp(j,imr) = qtmp(j,imr) - d1 |
|---|
| 2086 | qtmp(j,i) = qtmp(j,i) + d1 |
|---|
| 2087 | c east |
|---|
| 2088 | d0 = max(0.,qtmp(j,i+1)) |
|---|
| 2089 | d2 = min(-qtmp(j,i),d0) |
|---|
| 2090 | qtmp(j,i+1) = qtmp(j,i+1) - d2 |
|---|
| 2091 | c |
|---|
| 2092 | qtmp(j,i) = qtmp(j,i) + d2 + tiny |
|---|
| 2093 | endif |
|---|
| 2094 | END DO |
|---|
| 2095 | i=IMR |
|---|
| 2096 | do j=j1,j2 |
|---|
| 2097 | if(qtmp(j,i)<0.) then |
|---|
| 2098 | ipx = 1 |
|---|
| 2099 | c west |
|---|
| 2100 | d0 = max(0.,qtmp(j,i-1)) |
|---|
| 2101 | d1 = min(-qtmp(j,i),d0) |
|---|
| 2102 | qtmp(j,i-1) = qtmp(j,i-1) - d1 |
|---|
| 2103 | qtmp(j,i) = qtmp(j,i) + d1 |
|---|
| 2104 | c east |
|---|
| 2105 | d0 = max(0.,qtmp(j,1)) |
|---|
| 2106 | d2 = min(-qtmp(j,i),d0) |
|---|
| 2107 | qtmp(j,1) = qtmp(j,1) - d2 |
|---|
| 2108 | c |
|---|
| 2109 | qtmp(j,i) = qtmp(j,i) + d2 + tiny |
|---|
| 2110 | endif |
|---|
| 2111 | END DO |
|---|
| 2112 | C |
|---|
| 2113 | if(ipx/=0) then |
|---|
| 2114 | do j=j1,j2 |
|---|
| 2115 | do i=1,imr |
|---|
| 2116 | q(i,j) = qtmp(j,i) |
|---|
| 2117 | END DO |
|---|
| 2118 | END DO |
|---|
| 2119 | else |
|---|
| 2120 | C |
|---|
| 2121 | C Poles. |
|---|
| 2122 | if(q(1,1)<0 .or. q(1,JNP)<0.) ipx = 1 |
|---|
| 2123 | endif |
|---|
| 2124 | return |
|---|
| 2125 | end |
|---|
| 2126 | C |
|---|
| 2127 | subroutine zflip(q,im,km,nc) |
|---|
| 2128 | implicit none |
|---|
| 2129 | C This routine flip the array q (in the vertical). |
|---|
| 2130 | integer :: im,km,nc |
|---|
| 2131 | real q(im,km,nc) |
|---|
| 2132 | C local dynamic array |
|---|
| 2133 | real qtmp(im,km) |
|---|
| 2134 | integer IC,k,i |
|---|
| 2135 | C |
|---|
| 2136 | do IC = 1, nc |
|---|
| 2137 | C |
|---|
| 2138 | do k=1,km |
|---|
| 2139 | do i=1,im |
|---|
| 2140 | qtmp(i,k) = q(i,km+1-k,IC) |
|---|
| 2141 | END DO |
|---|
| 2142 | END DO |
|---|
| 2143 | C |
|---|
| 2144 | do i=1,im*km |
|---|
| 2145 | q(i,1,IC) = qtmp(i,1) |
|---|
| 2146 | END DO |
|---|
| 2147 | END DO |
|---|
| 2148 | return |
|---|
| 2149 | end |
|---|
