/[lmdze]/trunk/dyn3d/advy.f
ViewVC logotype

Contents of /trunk/dyn3d/advy.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 81 - (show annotations)
Wed Mar 5 14:38:41 2014 UTC (10 years, 2 months ago) by guez
Original Path: trunk/dyn3d/advy.f90
File size: 10458 byte(s) 
 Converted to free source form files which were still in fixed source
form. The conversion was done using the polish mode of the NAG Fortran
Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).
1
2 ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advy.F,v 1.1.1.1 2004/05/19
3 ! 12:53:06 lmdzadmin Exp $
4
5 SUBROUTINE advy(limit, dty, pbarv, sm, s0, sx, sy, sz)
6 USE dimens_m
7 USE paramet_m
8 USE comconst
9 USE disvert_m
10 USE comgeom
11 IMPLICIT NONE
12
13 ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
14 ! C
15 ! first-order moments (SOM) advection of tracer in Y direction C
16 ! C
17 ! Source : Pascal Simon ( Meteo, CNRM ) C
18 ! Adaptation : A.A. (LGGE) C
19 ! Derniere Modif : 15/12/94 LAST
20 ! C
21 ! sont les arguments d'entree pour le s-pg C
22 ! C
23 ! argument de sortie du s-pg C
24 ! C
25 ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
26 ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
27
28 ! Rem : Probleme aux poles il faut reecrire ce cas specifique
29 ! Attention au sens de l'indexation
30
31 ! parametres principaux du modele
32
33
34
35 ! Arguments :
36 ! ----------
37 ! dty : frequence fictive d'appel du transport
38 ! parbu,pbarv : flux de masse en x et y en Pa.m2.s-1
39
40 INTEGER lon, lat, niv
41 INTEGER i, j, jv, k, kp, l
42 INTEGER ntra
43 PARAMETER (ntra=1)
44
45 REAL dty
46 REAL, INTENT (IN) :: pbarv(iip1, jjm, llm)
47
48 ! moments: SM total mass in each grid box
49 ! S0 mass of tracer in each grid box
50 ! Si 1rst order moment in i direction
51
52 REAL sm(iip1, jjp1, llm), s0(iip1, jjp1, llm, ntra)
53 REAL sx(iip1, jjp1, llm, ntra), sy(iip1, jjp1, llm, ntra), &
54 sz(iip1, jjp1, llm, ntra)
55
56
57 ! Local :
58 ! -------
59
60 ! mass fluxes across the boundaries (UGRI,VGRI,WGRI)
61 ! mass fluxes in kg
62 ! declaration :
63
64 REAL vgri(iip1, 0:jjp1, llm)
65
66 ! Rem : UGRI et WGRI ne sont pas utilises dans
67 ! cette subroutine ( advection en y uniquement )
68 ! Rem 2 :le dimensionnement de VGRI depend de celui de pbarv
69
70 ! the moments F are similarly defined and used as temporary
71 ! storage for portions of the grid boxes in transit
72
73 REAL f0(iim, 0:jjp1, ntra), fm(iim, 0:jjp1)
74 REAL fx(iim, jjm, ntra), fy(iim, jjm, ntra)
75 REAL fz(iim, jjm, ntra)
76 REAL s00(ntra)
77 REAL sm0 ! Just temporal variable
78
79 ! work arrays
80
81 REAL alf(iim, 0:jjp1), alf1(iim, 0:jjp1)
82 REAL alfq(iim, 0:jjp1), alf1q(iim, 0:jjp1)
83 REAL temptm ! Just temporal variable
84
85 ! Special pour poles
86
87 REAL sbms, sfms, sfzs, sbmn, sfmn, sfzn
88 REAL sns0(ntra), snsz(ntra), snsm
89 REAL s1v(llm), slatv(llm)
90 REAL qy1(iim, llm, ntra), qylat(iim, llm, ntra)
91 REAL cx1(llm, ntra), cxlat(llm, ntra)
92 REAL cy1(llm, ntra), cylat(llm, ntra)
93 REAL z1(iim), zcos(iim), zsin(iim)
94 REAL smpn, smps, s0pn, s0ps
95 REAL ssum
96 EXTERNAL ssum
97
98 REAL sqi, sqf
99 LOGICAL limit
100
101 lon = iim ! rem : Il est possible qu'un pbl. arrive ici
102 lat = jjp1 ! a cause des dim. differentes entre les
103 niv = llm
104
105
106 ! the moments Fi are used as temporary storage for
107 ! portions of the grid boxes in transit at the current level
108
109 ! work arrays
110
111
112 DO l = 1, llm
113 DO j = 1, jjm
114 DO i = 1, iip1
115 vgri(i, j, llm+1-l) = -1.*pbarv(i, j, l)
116 END DO
117 END DO
118 DO i = 1, iip1
119 vgri(i, 0, l) = 0.
120 vgri(i, jjp1, l) = 0.
121 END DO
122 END DO
123
124 DO l = 1, niv
125
126 ! place limits on appropriate moments before transport
127 ! (if flux-limiting is to be applied)
128
129 IF (.NOT. limit) GO TO 11
130
131 DO jv = 1, ntra
132 DO k = 1, lat
133 DO i = 1, lon
134 sy(i, k, l, jv) = sign(amin1(amax1(s0(i,k,l,jv), &
135 0.),abs(sy(i,k,l,jv))), sy(i,k,l,jv))
136 END DO
137 END DO
138 END DO
139
140 11 CONTINUE
141
142 ! le flux a travers le pole Nord est traite separement
143
144 sm0 = 0.
145 DO jv = 1, ntra
146 s00(jv) = 0.
147 END DO
148
149 DO i = 1, lon
150
151 IF (vgri(i,0,l)<=0.) THEN
152 fm(i, 0) = -vgri(i, 0, l)*dty
153 alf(i, 0) = fm(i, 0)/sm(i, 1, l)
154 sm(i, 1, l) = sm(i, 1, l) - fm(i, 0)
155 sm0 = sm0 + fm(i, 0)
156 END IF
157
158 alfq(i, 0) = alf(i, 0)*alf(i, 0)
159 alf1(i, 0) = 1. - alf(i, 0)
160 alf1q(i, 0) = alf1(i, 0)*alf1(i, 0)
161
162 END DO
163
164 DO jv = 1, ntra
165 DO i = 1, lon
166
167 IF (vgri(i,0,l)<=0.) THEN
168
169 f0(i, 0, jv) = alf(i, 0)*(s0(i,1,l,jv)-alf1(i,0)*sy(i,1,l,jv))
170
171 s00(jv) = s00(jv) + f0(i, 0, jv)
172 s0(i, 1, l, jv) = s0(i, 1, l, jv) - f0(i, 0, jv)
173 sy(i, 1, l, jv) = alf1q(i, 0)*sy(i, 1, l, jv)
174 sx(i, 1, l, jv) = alf1(i, 0)*sx(i, 1, l, jv)
175 sz(i, 1, l, jv) = alf1(i, 0)*sz(i, 1, l, jv)
176
177 END IF
178
179 END DO
180 END DO
181
182 DO i = 1, lon
183 IF (vgri(i,0,l)>0.) THEN
184 fm(i, 0) = vgri(i, 0, l)*dty
185 alf(i, 0) = fm(i, 0)/sm0
186 END IF
187 END DO
188
189 DO jv = 1, ntra
190 DO i = 1, lon
191 IF (vgri(i,0,l)>0.) THEN
192 f0(i, 0, jv) = alf(i, 0)*s00(jv)
193 END IF
194 END DO
195 END DO
196
197 ! puts the temporary moments Fi into appropriate neighboring boxes
198
199 DO i = 1, lon
200
201 IF (vgri(i,0,l)>0.) THEN
202 sm(i, 1, l) = sm(i, 1, l) + fm(i, 0)
203 alf(i, 0) = fm(i, 0)/sm(i, 1, l)
204 END IF
205
206 alf1(i, 0) = 1. - alf(i, 0)
207
208 END DO
209
210 DO jv = 1, ntra
211 DO i = 1, lon
212
213 IF (vgri(i,0,l)>0.) THEN
214
215 temptm = alf(i, 0)*s0(i, 1, l, jv) - alf1(i, 0)*f0(i, 0, jv)
216 s0(i, 1, l, jv) = s0(i, 1, l, jv) + f0(i, 0, jv)
217 sy(i, 1, l, jv) = alf1(i, 0)*sy(i, 1, l, jv) + 3.*temptm
218
219 END IF
220
221 END DO
222 END DO
223
224 ! calculate flux and moments between adjacent boxes
225 ! 1- create temporary moments/masses for partial boxes in transit
226 ! 2- reajusts moments remaining in the box
227
228 ! flux from KP to K if V(K).lt.0 and from K to KP if V(K).gt.0
229
230 DO k = 1, lat - 1
231 kp = k + 1
232 DO i = 1, lon
233
234 IF (vgri(i,k,l)<0.) THEN
235 fm(i, k) = -vgri(i, k, l)*dty
236 alf(i, k) = fm(i, k)/sm(i, kp, l)
237 sm(i, kp, l) = sm(i, kp, l) - fm(i, k)
238 ELSE
239 fm(i, k) = vgri(i, k, l)*dty
240 alf(i, k) = fm(i, k)/sm(i, k, l)
241 sm(i, k, l) = sm(i, k, l) - fm(i, k)
242 END IF
243
244 alfq(i, k) = alf(i, k)*alf(i, k)
245 alf1(i, k) = 1. - alf(i, k)
246 alf1q(i, k) = alf1(i, k)*alf1(i, k)
247
248 END DO
249 END DO
250
251 DO jv = 1, ntra
252 DO k = 1, lat - 1
253 kp = k + 1
254 DO i = 1, lon
255
256 IF (vgri(i,k,l)<0.) THEN
257
258 f0(i, k, jv) = alf(i, k)*(s0(i,kp,l,jv)-alf1(i,k)*sy(i,kp,l,jv))
259 fy(i, k, jv) = alfq(i, k)*sy(i, kp, l, jv)
260 fx(i, k, jv) = alf(i, k)*sx(i, kp, l, jv)
261 fz(i, k, jv) = alf(i, k)*sz(i, kp, l, jv)
262
263 s0(i, kp, l, jv) = s0(i, kp, l, jv) - f0(i, k, jv)
264 sy(i, kp, l, jv) = alf1q(i, k)*sy(i, kp, l, jv)
265 sx(i, kp, l, jv) = sx(i, kp, l, jv) - fx(i, k, jv)
266 sz(i, kp, l, jv) = sz(i, kp, l, jv) - fz(i, k, jv)
267
268 ELSE
269
270 f0(i, k, jv) = alf(i, k)*(s0(i,k,l,jv)+alf1(i,k)*sy(i,k,l,jv))
271 fy(i, k, jv) = alfq(i, k)*sy(i, k, l, jv)
272 fx(i, k, jv) = alf(i, k)*sx(i, k, l, jv)
273 fz(i, k, jv) = alf(i, k)*sz(i, k, l, jv)
274
275 s0(i, k, l, jv) = s0(i, k, l, jv) - f0(i, k, jv)
276 sy(i, k, l, jv) = alf1q(i, k)*sy(i, k, l, jv)
277 sx(i, k, l, jv) = sx(i, k, l, jv) - fx(i, k, jv)
278 sz(i, k, l, jv) = sz(i, k, l, jv) - fz(i, k, jv)
279
280 END IF
281
282 END DO
283 END DO
284 END DO
285
286 ! puts the temporary moments Fi into appropriate neighboring boxes
287
288 DO k = 1, lat - 1
289 kp = k + 1
290 DO i = 1, lon
291
292 IF (vgri(i,k,l)<0.) THEN
293 sm(i, k, l) = sm(i, k, l) + fm(i, k)
294 alf(i, k) = fm(i, k)/sm(i, k, l)
295 ELSE
296 sm(i, kp, l) = sm(i, kp, l) + fm(i, k)
297 alf(i, k) = fm(i, k)/sm(i, kp, l)
298 END IF
299
300 alf1(i, k) = 1. - alf(i, k)
301
302 END DO
303 END DO
304
305 DO jv = 1, ntra
306 DO k = 1, lat - 1
307 kp = k + 1
308 DO i = 1, lon
309
310 IF (vgri(i,k,l)<0.) THEN
311
312 temptm = -alf(i, k)*s0(i, k, l, jv) + alf1(i, k)*f0(i, k, jv)
313 s0(i, k, l, jv) = s0(i, k, l, jv) + f0(i, k, jv)
314 sy(i, k, l, jv) = alf(i, k)*fy(i, k, jv) + &
315 alf1(i, k)*sy(i, k, l, jv) + 3.*temptm
316 sx(i, k, l, jv) = sx(i, k, l, jv) + fx(i, k, jv)
317 sz(i, k, l, jv) = sz(i, k, l, jv) + fz(i, k, jv)
318
319 ELSE
320
321 temptm = alf(i, k)*s0(i, kp, l, jv) - alf1(i, k)*f0(i, k, jv)
322 s0(i, kp, l, jv) = s0(i, kp, l, jv) + f0(i, k, jv)
323 sy(i, kp, l, jv) = alf(i, k)*fy(i, k, jv) + &
324 alf1(i, k)*sy(i, kp, l, jv) + 3.*temptm
325 sx(i, kp, l, jv) = sx(i, kp, l, jv) + fx(i, k, jv)
326 sz(i, kp, l, jv) = sz(i, kp, l, jv) + fz(i, k, jv)
327
328 END IF
329
330 END DO
331 END DO
332 END DO
333
334 ! traitement special pour le pole Sud (idem pole Nord)
335
336 k = lat
337
338 sm0 = 0.
339 DO jv = 1, ntra
340 s00(jv) = 0.
341 END DO
342
343 DO i = 1, lon
344
345 IF (vgri(i,k,l)>=0.) THEN
346 fm(i, k) = vgri(i, k, l)*dty
347 alf(i, k) = fm(i, k)/sm(i, k, l)
348 sm(i, k, l) = sm(i, k, l) - fm(i, k)
349 sm0 = sm0 + fm(i, k)
350 END IF
351
352 alfq(i, k) = alf(i, k)*alf(i, k)
353 alf1(i, k) = 1. - alf(i, k)
354 alf1q(i, k) = alf1(i, k)*alf1(i, k)
355
356 END DO
357
358 DO jv = 1, ntra
359 DO i = 1, lon
360
361 IF (vgri(i,k,l)>=0.) THEN
362 f0(i, k, jv) = alf(i, k)*(s0(i,k,l,jv)+alf1(i,k)*sy(i,k,l,jv))
363 s00(jv) = s00(jv) + f0(i, k, jv)
364
365 s0(i, k, l, jv) = s0(i, k, l, jv) - f0(i, k, jv)
366 sy(i, k, l, jv) = alf1q(i, k)*sy(i, k, l, jv)
367 sx(i, k, l, jv) = alf1(i, k)*sx(i, k, l, jv)
368 sz(i, k, l, jv) = alf1(i, k)*sz(i, k, l, jv)
369 END IF
370
371 END DO
372 END DO
373
374 DO i = 1, lon
375 IF (vgri(i,k,l)<0.) THEN
376 fm(i, k) = -vgri(i, k, l)*dty
377 alf(i, k) = fm(i, k)/sm0
378 END IF
379 END DO
380
381 DO jv = 1, ntra
382 DO i = 1, lon
383 IF (vgri(i,k,l)<0.) THEN
384 f0(i, k, jv) = alf(i, k)*s00(jv)
385 END IF
386 END DO
387 END DO
388
389 ! puts the temporary moments Fi into appropriate neighboring boxes
390
391 DO i = 1, lon
392
393 IF (vgri(i,k,l)<0.) THEN
394 sm(i, k, l) = sm(i, k, l) + fm(i, k)
395 alf(i, k) = fm(i, k)/sm(i, k, l)
396 END IF
397
398 alf1(i, k) = 1. - alf(i, k)
399
400 END DO
401
402 DO jv = 1, ntra
403 DO i = 1, lon
404
405 IF (vgri(i,k,l)<0.) THEN
406
407 temptm = -alf(i, k)*s0(i, k, l, jv) + alf1(i, k)*f0(i, k, jv)
408 s0(i, k, l, jv) = s0(i, k, l, jv) + f0(i, k, jv)
409 sy(i, k, l, jv) = alf1(i, k)*sy(i, k, l, jv) + 3.*temptm
410
411 END IF
412
413 END DO
414 END DO
415
416 END DO
417
418 RETURN
419 END SUBROUTINE advy
420
  ViewVC Help
Powered by ViewVC 1.1.21