/[lmdze]/trunk/Sources/dyn3d/prather.f
ViewVC logotype

Contents of /trunk/Sources/dyn3d/prather.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 157 - (show annotations)
Mon Jul 20 16:01:49 2015 UTC (8 years, 10 months ago) by guez
File size: 10565 byte(s) 
Just encapsulated SUBROUTINE vlsplt in a module and cleaned it.

In procedure vlx, local variables dxqu and adxqu only need indices
iip2:ip1jm. Otherwise, just cleaned vlx.

Procedures dynredem0 and dynredem1 no longer have argument fichnom,
they just operate on a file named "restart.nc". The programming
guideline here is that gcm should not be more complex than it needs by
itself, other programs (ce0l etc.) just have to adapt to gcm. So ce0l
now creates files "restart.nc" and "restartphy.nc".

In order to facilitate decentralizing the writing of "restartphy.nc",
created a procedure phyredem0 out of phyredem. phyredem0 creates the
NetCDF header of "restartphy.nc" while phyredem writes the NetCDF
variables. As the global attribute itau_phy needs to be filled in
phyredem0, at the beginnig of the run, we must compute its value
instead of just using itap. So we have a dummy argument lmt_pas of
phyredem0. Also, the ncid of "startphy.nc" is upgraded from local
variable of phyetat0 to dummy argument. phyetat0 no longer closes
"startphy.nc".

Following the same decentralizing objective, the ncid of "restart.nc"
is upgraded from local variable of dynredem0 to module variable of
dynredem0_m. "restart.nc" is not closed at the end of dynredem0 nor
opened at the beginning of dynredem1.

In procedure etat0, instead of creating many vectors of size klon
which will be filled with zeroes, just create one array null_array.

In procedure phytrac, instead of writing trs(: 1) to a text file,
write it to "restartphy.nc" (following LMDZ). This is better because
now trs(: 1) is next to its coordinates. We can write to
"restartphy.nc" from phytrac directly, and not add trs(: 1) to the
long list of variables in physiq, thanks to the decentralizing of
"restartphy.nc".

In procedure phyetat0, we no longer write to standard output the
minimum and maximum values of read arrays. It is ok to check input and
abort on invalid values but just printing statistics on input seems too
much useless computation and out of place clutter.
1
2 ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/prather.F,v 1.1.1.1 2004/05/19
3 ! 12:53:07 lmdzadmin Exp $
4
5 SUBROUTINE prather(q, w, masse, pbaru, pbarv, nt, dt)
6
7 USE comconst
8 USE dimens_m
9 USE disvert_m
10 USE dynetat0_m, only: rlonv, rlonu
11 USE nr_util, ONLY: pi
12 USE paramet_m
13
14 IMPLICIT NONE
15
16 ! =======================================================================
17 ! Adaptation LMDZ: A.Armengaud (LGGE)
18 ! ----------------
19
20 ! ************************************************
21 ! Transport des traceurs par la methode de prather
22 ! Ref :
23
24 ! ************************************************
25 ! q,w,pext,pbaru et pbarv : arguments d'entree pour le s-pg
26
27 ! =======================================================================
28
29
30
31 ! Arguments:
32 ! ----------
33 INTEGER nt
34 REAL, INTENT (IN) :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
35 REAL masse(iip1, jjp1, llm)
36 REAL q(iip1, jjp1, llm, 0:9)
37 REAL w(ip1jmp1, llm)
38
39 ! Local:
40 ! ------
41 LOGICAL limit
42 REAL sm(iip1, jjp1, llm)
43 REAL s0(iip1, jjp1, llm), sx(iip1, jjp1, llm)
44 REAL sy(iip1, jjp1, llm), sz(iip1, jjp1, llm)
45 REAL sxx(iip1, jjp1, llm)
46 REAL sxy(iip1, jjp1, llm)
47 REAL sxz(iip1, jjp1, llm)
48 REAL syy(iip1, jjp1, llm)
49 REAL syz(iip1, jjp1, llm)
50 REAL szz(iip1, jjp1, llm), zz
51 INTEGER i, j, l, indice
52 REAL sxn(iip1), sxs(iip1)
53
54 REAL sinlon(iip1), sinlondlon(iip1)
55 REAL coslon(iip1), coslondlon(iip1)
56 SAVE sinlon, coslon, sinlondlon, coslondlon
57 REAL dyn1, dyn2, dys1, dys2, qpn, qps, dqzpn, dqzps
58 REAL masn, mass
59
60 REAL ssum
61 INTEGER ismax, ismin
62 EXTERNAL ssum, convflu, ismin, ismax
63 LOGICAL first
64 SAVE first
65 EXTERNAL advxp, advyp, advzp
66
67
68 DATA first/.TRUE./
69
70 ! ==========================================================================
71 ! ==========================================================================
72 ! MODIFICATION POUR PAS DE TEMPS ADAPTATIF, dtvr remplace par dt
73 ! ==========================================================================
74 ! ==========================================================================
75 REAL dt
76 ! ==========================================================================
77 limit = .TRUE.
78
79 IF (first) THEN
80 PRINT *, 'SCHEMA PRATHER'
81 first = .FALSE.
82 DO i = 2, iip1
83 coslon(i) = cos(rlonv(i))
84 sinlon(i) = sin(rlonv(i))
85 coslondlon(i) = coslon(i)*(rlonu(i)-rlonu(i-1))/pi
86 sinlondlon(i) = sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
87 END DO
88 coslon(1) = coslon(iip1)
89 coslondlon(1) = coslondlon(iip1)
90 sinlon(1) = sinlon(iip1)
91 sinlondlon(1) = sinlondlon(iip1)
92
93 DO l = 1, llm
94 DO j = 1, jjp1
95 DO i = 1, iip1
96 q(i, j, l, 1) = 0.
97 q(i, j, l, 2) = 0.
98 q(i, j, l, 3) = 0.
99 q(i, j, l, 4) = 0.
100 q(i, j, l, 5) = 0.
101 q(i, j, l, 6) = 0.
102 q(i, j, l, 7) = 0.
103 q(i, j, l, 8) = 0.
104 q(i, j, l, 9) = 0.
105 END DO
106 END DO
107 END DO
108 END IF
109 ! Fin modif Fred
110
111 ! *** On calcule la masse d'air en kg
112
113 DO l = 1, llm
114 DO j = 1, jjp1
115 DO i = 1, iip1
116 sm(i, j, llm+1-l) = masse(i, j, l)
117 END DO
118 END DO
119 END DO
120
121 ! *** q contient les qqtes de traceur avant l'advection
122
123 ! *** Affectation des tableaux S a partir de Q
124
125 DO l = 1, llm
126 DO j = 1, jjp1
127 DO i = 1, iip1
128 s0(i, j, l) = q(i, j, llm+1-l, 0)*sm(i, j, l)
129 sx(i, j, l) = q(i, j, llm+1-l, 1)*sm(i, j, l)
130 sy(i, j, l) = q(i, j, llm+1-l, 2)*sm(i, j, l)
131 sz(i, j, l) = q(i, j, llm+1-l, 3)*sm(i, j, l)
132 sxx(i, j, l) = q(i, j, llm+1-l, 4)*sm(i, j, l)
133 sxy(i, j, l) = q(i, j, llm+1-l, 5)*sm(i, j, l)
134 sxz(i, j, l) = q(i, j, llm+1-l, 6)*sm(i, j, l)
135 syy(i, j, l) = q(i, j, llm+1-l, 7)*sm(i, j, l)
136 syz(i, j, l) = q(i, j, llm+1-l, 8)*sm(i, j, l)
137 szz(i, j, l) = q(i, j, llm+1-l, 9)*sm(i, j, l)
138 END DO
139 END DO
140 END DO
141 ! *** Appel des subroutines d'advection en X, en Y et en Z
142 ! *** Advection avec "time-splitting"
143
144 ! -----------------------------------------------------------
145 DO indice = 1, nt
146 CALL advxp(limit, 0.5*dt, pbaru, sm, s0, sx, sy, sz, sxx, sxy, sxz, syy, &
147 syz, szz, 1)
148 END DO
149 DO l = 1, llm
150 DO i = 1, iip1
151 sy(i, 1, l) = 0.
152 sy(i, jjp1, l) = 0.
153 END DO
154 END DO
155 ! ---------------------------------------------------------
156 CALL advyp(limit, .5*dt*nt, pbarv, sm, s0, sx, sy, sz, sxx, sxy, sxz, syy, &
157 syz, szz, 1)
158 ! ---------------------------------------------------------
159
160 ! ---------------------------------------------------------
161 DO j = 1, jjp1
162 DO i = 1, iip1
163 sz(i, j, 1) = 0.
164 sz(i, j, llm) = 0.
165 sxz(i, j, 1) = 0.
166 sxz(i, j, llm) = 0.
167 syz(i, j, 1) = 0.
168 syz(i, j, llm) = 0.
169 szz(i, j, 1) = 0.
170 szz(i, j, llm) = 0.
171 END DO
172 END DO
173 CALL advzp(limit, dt*nt, w, sm, s0, sx, sy, sz, sxx, sxy, sxz, syy, syz, &
174 szz, 1)
175 DO l = 1, llm
176 DO i = 1, iip1
177 sy(i, 1, l) = 0.
178 sy(i, jjp1, l) = 0.
179 END DO
180 END DO
181
182 ! ---------------------------------------------------------
183
184 ! ---------------------------------------------------------
185 CALL advyp(limit, .5*dt*nt, pbarv, sm, s0, sx, sy, sz, sxx, sxy, sxz, syy, &
186 syz, szz, 1)
187 ! ---------------------------------------------------------
188 DO l = 1, llm
189 DO j = 1, jjp1
190 s0(iip1, j, l) = s0(1, j, l)
191 sx(iip1, j, l) = sx(1, j, l)
192 sy(iip1, j, l) = sy(1, j, l)
193 sz(iip1, j, l) = sz(1, j, l)
194 sxx(iip1, j, l) = sxx(1, j, l)
195 sxy(iip1, j, l) = sxy(1, j, l)
196 sxz(iip1, j, l) = sxz(1, j, l)
197 syy(iip1, j, l) = syy(1, j, l)
198 syz(iip1, j, l) = syz(1, j, l)
199 szz(iip1, j, l) = szz(1, j, l)
200 END DO
201 END DO
202 DO indice = 1, nt
203 CALL advxp(limit, 0.5*dt, pbaru, sm, s0, sx, sy, sz, sxx, sxy, sxz, syy, &
204 syz, szz, 1)
205 END DO
206 ! ---------------------------------------------------------
207 ! ---------------------------------------------------------
208 ! *** On repasse les S dans la variable qpr
209 ! *** On repasse les S dans la variable q directement 14/10/94
210
211 DO l = 1, llm
212 DO j = 1, jjp1
213 DO i = 1, iip1
214 q(i, j, llm+1-l, 0) = s0(i, j, l)/sm(i, j, l)
215 q(i, j, llm+1-l, 1) = sx(i, j, l)/sm(i, j, l)
216 q(i, j, llm+1-l, 2) = sy(i, j, l)/sm(i, j, l)
217 q(i, j, llm+1-l, 3) = sz(i, j, l)/sm(i, j, l)
218 q(i, j, llm+1-l, 4) = sxx(i, j, l)/sm(i, j, l)
219 q(i, j, llm+1-l, 5) = sxy(i, j, l)/sm(i, j, l)
220 q(i, j, llm+1-l, 6) = sxz(i, j, l)/sm(i, j, l)
221 q(i, j, llm+1-l, 7) = syy(i, j, l)/sm(i, j, l)
222 q(i, j, llm+1-l, 8) = syz(i, j, l)/sm(i, j, l)
223 q(i, j, llm+1-l, 9) = szz(i, j, l)/sm(i, j, l)
224 END DO
225 END DO
226 END DO
227
228 ! ---------------------------------------------------------
229 ! go to 777
230 ! filtrages aux poles
231
232 ! Traitements specifiques au pole
233
234 ! filtrages aux poles
235 DO l = 1, llm
236 ! filtrages aux poles
237 masn = ssum(iim, sm(1,1,l), 1)
238 mass = ssum(iim, sm(1,jjp1,l), 1)
239 qpn = ssum(iim, s0(1,1,l), 1)/masn
240 qps = ssum(iim, s0(1,jjp1,l), 1)/mass
241 dqzpn = ssum(iim, sz(1,1,l), 1)/masn
242 dqzps = ssum(iim, sz(1,jjp1,l), 1)/mass
243 DO i = 1, iip1
244 q(i, 1, llm+1-l, 3) = dqzpn
245 q(i, jjp1, llm+1-l, 3) = dqzps
246 q(i, 1, llm+1-l, 0) = qpn
247 q(i, jjp1, llm+1-l, 0) = qps
248 END DO
249 dyn1 = 0.
250 dys1 = 0.
251 dyn2 = 0.
252 dys2 = 0.
253 DO i = 1, iim
254 zz = s0(i, 2, l)/sm(i, 2, l) - q(i, 1, llm+1-l, 0)
255 dyn1 = dyn1 + sinlondlon(i)*zz
256 dyn2 = dyn2 + coslondlon(i)*zz
257 zz = q(i, jjp1, llm+1-l, 0) - s0(i, jjm, l)/sm(i, jjm, l)
258 dys1 = dys1 + sinlondlon(i)*zz
259 dys2 = dys2 + coslondlon(i)*zz
260 END DO
261 DO i = 1, iim
262 q(i, 1, llm+1-l, 2) = (sinlon(i)*dyn1+coslon(i)*dyn2)/2.
263 q(i, 1, llm+1-l, 0) = q(i, 1, llm+1-l, 0) + q(i, 1, llm+1-l, 2)
264 q(i, jjp1, llm+1-l, 2) = (sinlon(i)*dys1+coslon(i)*dys2)/2.
265 q(i, jjp1, llm+1-l, 0) = q(i, jjp1, llm+1-l, 0) - &
266 q(i, jjp1, llm+1-l, 2)
267 END DO
268 q(iip1, 1, llm+1-l, 0) = q(1, 1, llm+1-l, 0)
269 q(iip1, jjp1, llm+1-l, 0) = q(1, jjp1, llm+1-l, 0)
270 DO i = 1, iim
271 sxn(i) = q(i+1, 1, llm+1-l, 0) - q(i, 1, llm+1-l, 0)
272 sxs(i) = q(i+1, jjp1, llm+1-l, 0) - q(i, jjp1, llm+1-l, 0)
273 END DO
274 sxn(iip1) = sxn(1)
275 sxs(iip1) = sxs(1)
276 DO i = 1, iim
277 q(i+1, 1, llm+1-l, 1) = 0.25*(sxn(i)+sxn(i+1))
278 q(i+1, jjp1, llm+1-l, 1) = 0.25*(sxs(i)+sxs(i+1))
279 END DO
280 q(1, 1, llm+1-l, 1) = q(iip1, 1, llm+1-l, 1)
281 q(1, jjp1, llm+1-l, 1) = q(iip1, jjp1, llm+1-l, 1)
282 END DO
283 DO l = 1, llm
284 DO i = 1, iim
285 q(i, 1, llm+1-l, 4) = 0.
286 q(i, jjp1, llm+1-l, 4) = 0.
287 q(i, 1, llm+1-l, 5) = 0.
288 q(i, jjp1, llm+1-l, 5) = 0.
289 q(i, 1, llm+1-l, 6) = 0.
290 q(i, jjp1, llm+1-l, 6) = 0.
291 q(i, 1, llm+1-l, 7) = 0.
292 q(i, jjp1, llm+1-l, 7) = 0.
293 q(i, 1, llm+1-l, 8) = 0.
294 q(i, jjp1, llm+1-l, 8) = 0.
295 q(i, 1, llm+1-l, 9) = 0.
296 q(i, jjp1, llm+1-l, 9) = 0.
297 END DO
298 END DO
299
300 ! bouclage en longitude
301 DO l = 1, llm
302 DO j = 1, jjp1
303 q(iip1, j, l, 0) = q(1, j, l, 0)
304 q(iip1, j, llm+1-l, 0) = q(1, j, llm+1-l, 0)
305 q(iip1, j, llm+1-l, 1) = q(1, j, llm+1-l, 1)
306 q(iip1, j, llm+1-l, 2) = q(1, j, llm+1-l, 2)
307 q(iip1, j, llm+1-l, 3) = q(1, j, llm+1-l, 3)
308 q(iip1, j, llm+1-l, 4) = q(1, j, llm+1-l, 4)
309 q(iip1, j, llm+1-l, 5) = q(1, j, llm+1-l, 5)
310 q(iip1, j, llm+1-l, 6) = q(1, j, llm+1-l, 6)
311 q(iip1, j, llm+1-l, 7) = q(1, j, llm+1-l, 7)
312 q(iip1, j, llm+1-l, 8) = q(1, j, llm+1-l, 8)
313 q(iip1, j, llm+1-l, 9) = q(1, j, llm+1-l, 9)
314 END DO
315 END DO
316 DO l = 1, llm
317 DO j = 2, jjm
318 DO i = 1, iip1
319 IF (q(i,j,l,0)<0.) THEN
320 PRINT *, '------------ BIP-----------'
321 PRINT *, 'S0(', i, j, l, ')=', q(i, j, l, 0), q(i, j-1, l, 0)
322 PRINT *, 'SX(', i, j, l, ')=', q(i, j, l, 1)
323 PRINT *, 'SY(', i, j, l, ')=', q(i, j, l, 2), q(i, j-1, l, 2)
324 PRINT *, 'SZ(', i, j, l, ')=', q(i, j, l, 3)
325 q(i, j, l, 0) = 0.
326 END IF
327 END DO
328 END DO
329 DO j = 1, jjp1, jjm
330 DO i = 1, iip1
331 IF (q(i,j,l,0)<0.) THEN
332 PRINT *, '------------ BIP 2-----------'
333 PRINT *, 'S0(', i, j, l, ')=', q(i, j, l, 0)
334 PRINT *, 'SX(', i, j, l, ')=', q(i, j, l, 1)
335 PRINT *, 'SY(', i, j, l, ')=', q(i, j, l, 2)
336 PRINT *, 'SZ(', i, j, l, ')=', q(i, j, l, 3)
337
338 q(i, j, l, 0) = 0.
339 ! STOP
340 END IF
341 END DO
342 END DO
343 END DO
344 RETURN
345 END SUBROUTINE prather
  ViewVC Help
Powered by ViewVC 1.1.21