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