GCC Code Coverage Report


Directory: ./
File: dyn3d_common/advx.f
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 224 0.0%
Branches: 0 134 0.0%

Line Branch Exec Source
1 !
2 ! $Header$
3 !
4 SUBROUTINE advx(limit,dtx,pbaru,sm,s0,
5 $ sx,sy,sz,lati,latf)
6 IMPLICIT NONE
7
8 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
9 C C
10 C first-order moments (FOM) advection of tracer in X direction C
11 C C
12 C Source : Pascal Simon (Meteo,CNRM) C
13 C Adaptation : A.Armengaud (LGGE) juin 94 C
14 C C
15 C limit,dtx,pbaru,pbarv,sm,s0,sx,sy,sz C
16 C sont des arguments d'entree pour le s-pg... C
17 C C
18 C sm,s0,sx,sy,sz C
19 C sont les arguments de sortie pour le s-pg C
20 C C
21 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
22 C
23 C parametres principaux du modele
24 C
25 include "dimensions.h"
26 include "paramet.h"
27
28 C Arguments :
29 C -----------
30 C dtx : frequence fictive d'appel du transport
31 C pbaru, pbarv : flux de masse en x et y en Pa.m2.s-1
32
33 INTEGER ntra
34 PARAMETER (ntra = 1)
35
36 C ATTENTION partout ou on trouve ntra, insertion de boucle
37 C possible dans l'avenir.
38
39 REAL dtx
40 REAL pbaru ( iip1,jjp1,llm )
41
42 C moments: SM total mass in each grid box
43 C S0 mass of tracer in each grid box
44 C Si 1rst order moment in i direction
45 C
46 REAL SM(iip1,jjp1,llm),S0(iip1,jjp1,llm,ntra)
47 REAL sx(iip1,jjp1,llm,ntra)
48 $ ,sy(iip1,jjp1,llm,ntra)
49 REAL sz(iip1,jjp1,llm,ntra)
50
51 C Local :
52 C -------
53
54 C mass fluxes across the boundaries (UGRI,VGRI,WGRI)
55 C mass fluxes in kg
56 C declaration :
57
58 REAL UGRI(iip1,jjp1,llm)
59
60 C Rem : VGRI et WGRI ne sont pas utilises dans
61 C cette subroutine ( advection en x uniquement )
62 C
63 C Ti are the moments for the current latitude and level
64 C
65 REAL TM(iim)
66 REAL T0(iim,ntra),TX(iim,ntra)
67 REAL TY(iim,ntra),TZ(iim,ntra)
68 REAL TEMPTM ! just a temporary variable
69 C
70 C the moments F are similarly defined and used as temporary
71 C storage for portions of the grid boxes in transit
72 C
73 REAL FM(iim)
74 REAL F0(iim,ntra),FX(iim,ntra)
75 REAL FY(iim,ntra),FZ(iim,ntra)
76 C
77 C work arrays
78 C
79 REAL ALF(iim),ALF1(iim),ALFQ(iim),ALF1Q(iim)
80 C
81 REAL SMNEW(iim),UEXT(iim)
82 C
83 REAL sqi,sqf
84
85 LOGICAL LIMIT
86 INTEGER NUM(jjp1),LONK,NUMK
87 INTEGER lon,lati,latf,niv
88 INTEGER i,i2,i3,j,jv,l,k,itrac
89
90 lon = iim
91 niv = llm
92
93 C *** Test de passage d'arguments ******
94
95
96 C -------------------------------------
97 DO 300 j = 1,jjp1
98 NUM(j) = 1
99 300 CONTINUE
100 sqi = 0.
101 sqf = 0.
102
103 DO l = 1,llm
104 DO j = 1,jjp1
105 DO i = 1,iim
106 cIM 240305 sqi = sqi + S0(i,j,l,9)
107 sqi = sqi + S0(i,j,l,ntra)
108 ENDDO
109 ENDDO
110 ENDDO
111 PRINT*,'-------- DIAG DANS ADVX - ENTREE ---------'
112 PRINT*,'sqi=',sqi
113
114
115 C Interface : adaptation nouveau modele
116 C -------------------------------------
117 C
118 C ---------------------------------------------------------
119 C Conversion des flux de masses en kg/s
120 C pbaru est en N/s d'ou :
121 C ugri est en kg/s
122
123 DO 500 l = 1,llm
124 DO 500 j = 1,jjm+1
125 DO 500 i = 1,iip1
126 C ugri (i,j,llm+1-l) = pbaru (i,j,l) * ( dsig(l) / g )
127 ugri (i,j,llm+1-l) = pbaru (i,j,l)
128 500 CONTINUE
129
130
131 C ---------------------------------------------------------
132 C ---------------------------------------------------------
133 C ---------------------------------------------------------
134
135 C start here
136 C
137 C boucle principale sur les niveaux et les latitudes
138 C
139 DO 1 L=1,NIV
140 DO 1 K=lati,latf
141 C
142 C initialisation
143 C
144 C program assumes periodic boundaries in X
145 C
146 DO 10 I=2,LON
147 SMNEW(I)=SM(I,K,L)+(UGRI(I-1,K,L)-UGRI(I,K,L))*DTX
148 10 CONTINUE
149 SMNEW(1)=SM(1,K,L)+(UGRI(LON,K,L)-UGRI(1,K,L))*DTX
150 C
151 C modifications for extended polar zones
152 C
153 NUMK=NUM(K)
154 LONK=LON/NUMK
155 C
156 IF(NUMK.GT.1) THEN
157 C
158 DO 111 I=1,LON
159 TM(I)=0.
160 111 CONTINUE
161 DO 112 JV=1,NTRA
162 DO 1120 I=1,LON
163 T0(I,JV)=0.
164 TX(I,JV)=0.
165 TY(I,JV)=0.
166 TZ(I,JV)=0.
167 1120 CONTINUE
168 112 CONTINUE
169 C
170 DO 11 I2=1,NUMK
171 C
172 DO 113 I=1,LONK
173 I3=(I-1)*NUMK+I2
174 TM(I)=TM(I)+SM(I3,K,L)
175 ALF(I)=SM(I3,K,L)/TM(I)
176 ALF1(I)=1.-ALF(I)
177 113 CONTINUE
178 C
179 DO JV=1,NTRA
180 DO I=1,LONK
181 I3=(I-1)*NUMK+I2
182 TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)
183 $ *S0(I3,K,L,JV)
184 T0(I,JV)=T0(I,JV)+S0(I3,K,L,JV)
185 TX(I,JV)=ALF(I) *sx(I3,K,L,JV)+
186 $ ALF1(I)*TX(I,JV) +3.*TEMPTM
187 TY(I,JV)=TY(I,JV)+sy(I3,K,L,JV)
188 TZ(I,JV)=TZ(I,JV)+sz(I3,K,L,JV)
189 ENDDO
190 ENDDO
191 C
192 11 CONTINUE
193 C
194 ELSE
195 C
196 DO 115 I=1,LON
197 TM(I)=SM(I,K,L)
198 115 CONTINUE
199 DO 116 JV=1,NTRA
200 DO 1160 I=1,LON
201 T0(I,JV)=S0(I,K,L,JV)
202 TX(I,JV)=sx(I,K,L,JV)
203 TY(I,JV)=sy(I,K,L,JV)
204 TZ(I,JV)=sz(I,K,L,JV)
205 1160 CONTINUE
206 116 CONTINUE
207 C
208 ENDIF
209 C
210 DO 117 I=1,LONK
211 UEXT(I)=UGRI(I*NUMK,K,L)
212 117 CONTINUE
213 C
214 C place limits on appropriate moments before transport
215 C (if flux-limiting is to be applied)
216 C
217 IF(.NOT.LIMIT) GO TO 13
218 C
219 DO 12 JV=1,NTRA
220 DO 120 I=1,LONK
221 TX(I,JV)=SIGN(AMIN1(AMAX1(T0(I,JV),0.),ABS(TX(I,JV))),TX(I,JV))
222 120 CONTINUE
223 12 CONTINUE
224 C
225 13 CONTINUE
226 C
227 C calculate flux and moments between adjacent boxes
228 C 1- create temporary moments/masses for partial boxes in transit
229 C 2- reajusts moments remaining in the box
230 C
231 C flux from IP to I if U(I).lt.0
232 C
233 DO 140 I=1,LONK-1
234 IF(UEXT(I).LT.0.) THEN
235 FM(I)=-UEXT(I)*DTX
236 ALF(I)=FM(I)/TM(I+1)
237 TM(I+1)=TM(I+1)-FM(I)
238 ENDIF
239 140 CONTINUE
240 C
241 I=LONK
242 IF(UEXT(I).LT.0.) THEN
243 FM(I)=-UEXT(I)*DTX
244 ALF(I)=FM(I)/TM(1)
245 TM(1)=TM(1)-FM(I)
246 ENDIF
247 C
248 C flux from I to IP if U(I).gt.0
249 C
250 DO 141 I=1,LONK
251 IF(UEXT(I).GE.0.) THEN
252 FM(I)=UEXT(I)*DTX
253 ALF(I)=FM(I)/TM(I)
254 TM(I)=TM(I)-FM(I)
255 ENDIF
256 141 CONTINUE
257 C
258 DO 142 I=1,LONK
259 ALFQ(I)=ALF(I)*ALF(I)
260 ALF1(I)=1.-ALF(I)
261 ALF1Q(I)=ALF1(I)*ALF1(I)
262 142 CONTINUE
263 C
264 DO 150 JV=1,NTRA
265 DO 1500 I=1,LONK-1
266 C
267 IF(UEXT(I).LT.0.) THEN
268 C
269 F0(I,JV)=ALF (I)* ( T0(I+1,JV)-ALF1(I)*TX(I+1,JV) )
270 FX(I,JV)=ALFQ(I)*TX(I+1,JV)
271 FY(I,JV)=ALF (I)*TY(I+1,JV)
272 FZ(I,JV)=ALF (I)*TZ(I+1,JV)
273 C
274 T0(I+1,JV)=T0(I+1,JV)-F0(I,JV)
275 TX(I+1,JV)=ALF1Q(I)*TX(I+1,JV)
276 TY(I+1,JV)=TY(I+1,JV)-FY(I,JV)
277 TZ(I+1,JV)=TZ(I+1,JV)-FZ(I,JV)
278 C
279 ENDIF
280 C
281 1500 CONTINUE
282 150 CONTINUE
283 C
284 I=LONK
285 IF(UEXT(I).LT.0.) THEN
286 C
287 DO 151 JV=1,NTRA
288 C
289 F0 (I,JV)=ALF (I)* ( T0(1,JV)-ALF1(I)*TX(1,JV) )
290 FX (I,JV)=ALFQ(I)*TX(1,JV)
291 FY (I,JV)=ALF (I)*TY(1,JV)
292 FZ (I,JV)=ALF (I)*TZ(1,JV)
293 C
294 T0(1,JV)=T0(1,JV)-F0(I,JV)
295 TX(1,JV)=ALF1Q(I)*TX(1,JV)
296 TY(1,JV)=TY(1,JV)-FY(I,JV)
297 TZ(1,JV)=TZ(1,JV)-FZ(I,JV)
298 C
299 151 CONTINUE
300 C
301 ENDIF
302 C
303 DO 152 JV=1,NTRA
304 DO 1520 I=1,LONK
305 C
306 IF(UEXT(I).GE.0.) THEN
307 C
308 F0(I,JV)=ALF (I)* ( T0(I,JV)+ALF1(I)*TX(I,JV) )
309 FX(I,JV)=ALFQ(I)*TX(I,JV)
310 FY(I,JV)=ALF (I)*TY(I,JV)
311 FZ(I,JV)=ALF (I)*TZ(I,JV)
312 C
313 T0(I,JV)=T0(I,JV)-F0(I,JV)
314 TX(I,JV)=ALF1Q(I)*TX(I,JV)
315 TY(I,JV)=TY(I,JV)-FY(I,JV)
316 TZ(I,JV)=TZ(I,JV)-FZ(I,JV)
317 C
318 ENDIF
319 C
320 1520 CONTINUE
321 152 CONTINUE
322 C
323 C puts the temporary moments Fi into appropriate neighboring boxes
324 C
325 DO 160 I=1,LONK
326 IF(UEXT(I).LT.0.) THEN
327 TM(I)=TM(I)+FM(I)
328 ALF(I)=FM(I)/TM(I)
329 ENDIF
330 160 CONTINUE
331 C
332 DO 161 I=1,LONK-1
333 IF(UEXT(I).GE.0.) THEN
334 TM(I+1)=TM(I+1)+FM(I)
335 ALF(I)=FM(I)/TM(I+1)
336 ENDIF
337 161 CONTINUE
338 C
339 I=LONK
340 IF(UEXT(I).GE.0.) THEN
341 TM(1)=TM(1)+FM(I)
342 ALF(I)=FM(I)/TM(1)
343 ENDIF
344 C
345 DO 162 I=1,LONK
346 ALF1(I)=1.-ALF(I)
347 162 CONTINUE
348 C
349 DO 170 JV=1,NTRA
350 DO 1700 I=1,LONK
351 C
352 IF(UEXT(I).LT.0.) THEN
353 C
354 TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)*F0(I,JV)
355 T0(I,JV)=T0(I,JV)+F0(I,JV)
356 TX(I,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(I,JV)+3.*TEMPTM
357 TY(I,JV)=TY(I,JV)+FY(I,JV)
358 TZ(I,JV)=TZ(I,JV)+FZ(I,JV)
359 C
360 ENDIF
361 C
362 1700 CONTINUE
363 170 CONTINUE
364 C
365 DO 171 JV=1,NTRA
366 DO 1710 I=1,LONK-1
367 C
368 IF(UEXT(I).GE.0.) THEN
369 C
370 TEMPTM=ALF(I)*T0(I+1,JV)-ALF1(I)*F0(I,JV)
371 T0(I+1,JV)=T0(I+1,JV)+F0(I,JV)
372 TX(I+1,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(I+1,JV)+3.*TEMPTM
373 TY(I+1,JV)=TY(I+1,JV)+FY(I,JV)
374 TZ(I+1,JV)=TZ(I+1,JV)+FZ(I,JV)
375 C
376 ENDIF
377 C
378 1710 CONTINUE
379 171 CONTINUE
380 C
381 I=LONK
382 IF(UEXT(I).GE.0.) THEN
383 DO 172 JV=1,NTRA
384 TEMPTM=ALF(I)*T0(1,JV)-ALF1(I)*F0(I,JV)
385 T0(1,JV)=T0(1,JV)+F0(I,JV)
386 TX(1,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(1,JV)+3.*TEMPTM
387 TY(1,JV)=TY(1,JV)+FY(I,JV)
388 TZ(1,JV)=TZ(1,JV)+FZ(I,JV)
389 172 CONTINUE
390 ENDIF
391 C
392 C retour aux mailles d'origine (passage des Tij aux Sij)
393 C
394 IF(NUMK.GT.1) THEN
395 C
396 DO 180 I2=1,NUMK
397 C
398 DO 180 I=1,LONK
399 C
400 I3=I2+(I-1)*NUMK
401 SM(I3,K,L)=SMNEW(I3)
402 ALF(I)=SMNEW(I3)/TM(I)
403 TM(I)=TM(I)-SMNEW(I3)
404 C
405 ALFQ(I)=ALF(I)*ALF(I)
406 ALF1(I)=1.-ALF(I)
407 ALF1Q(I)=ALF1(I)*ALF1(I)
408 C
409 180 CONTINUE
410 C
411 DO JV=1,NTRA
412 DO I=1,LONK
413 C
414 I3=I2+(I-1)*NUMK
415 S0(I3,K,L,JV)=ALF (I)
416 $ * (T0(I,JV)-ALF1(I)*TX(I,JV))
417 sx(I3,K,L,JV)=ALFQ(I)*TX(I,JV)
418 sy(I3,K,L,JV)=ALF (I)*TY(I,JV)
419 sz(I3,K,L,JV)=ALF (I)*TZ(I,JV)
420 C
421 C reajusts moments remaining in the box
422 C
423 T0(I,JV)=T0(I,JV)-S0(I3,K,L,JV)
424 TX(I,JV)=ALF1Q(I)*TX(I,JV)
425 TY(I,JV)=TY(I,JV)-sy(I3,K,L,JV)
426 TZ(I,JV)=TZ(I,JV)-sz(I3,K,L,JV)
427 ENDDO
428 ENDDO
429 C
430 C
431 ELSE
432 C
433 DO 190 I=1,LON
434 SM(I,K,L)=TM(I)
435 190 CONTINUE
436 DO 191 JV=1,NTRA
437 DO 1910 I=1,LON
438 S0(I,K,L,JV)=T0(I,JV)
439 sx(I,K,L,JV)=TX(I,JV)
440 sy(I,K,L,JV)=TY(I,JV)
441 sz(I,K,L,JV)=TZ(I,JV)
442 1910 CONTINUE
443 191 CONTINUE
444 C
445 ENDIF
446 C
447 1 CONTINUE
448 C
449 C ----------- AA Test en fin de ADVX ------ Controle des S*
450 c OK
451 c DO 9998 l = 1, llm
452 c DO 9998 j = 1, jjp1
453 c DO 9998 i = 1, iip1
454 c IF (S0(i,j,l,ntra).lt.0..and.LIMIT) THEN
455 c PRINT*, '-------------------'
456 c PRINT*, 'En fin de ADVX'
457 c PRINT*,'SM(',i,j,l,')=',SM(i,j,l)
458 c PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)
459 c print*, 'sx(',i,j,l,')=',sx(i,j,l,ntra)
460 c print*, 'sy(',i,j,l,')=',sy(i,j,l,ntra)
461 c print*, 'sz(',i,j,l,')=',sz(i,j,l,ntra)
462 c WRITE (*,*) 'On arrete !! - pbl en fin de ADVX1'
463 cc STOP
464 c ENDIF
465 c 9998 CONTINUE
466 c
467 C ---------- bouclage cyclique
468 DO itrac=1,ntra
469 DO l = 1,llm
470 DO j = lati,latf
471 SM(iip1,j,l) = SM(1,j,l)
472 S0(iip1,j,l,itrac) = S0(1,j,l,itrac)
473 sx(iip1,j,l,itrac) = sx(1,j,l,itrac)
474 sy(iip1,j,l,itrac) = sy(1,j,l,itrac)
475 sz(iip1,j,l,itrac) = sz(1,j,l,itrac)
476 END DO
477 END DO
478 ENDDO
479
480 c ----------- qqtite totale de traceur dans tte l'atmosphere
481 DO l = 1, llm
482 DO j = 1, jjp1
483 DO i = 1, iim
484 cIM 240405 sqf = sqf + S0(i,j,l,9)
485 sqf = sqf + S0(i,j,l,ntra)
486 END DO
487 END DO
488 END DO
489 c
490 PRINT*,'------ DIAG DANS ADVX - SORTIE -----'
491 PRINT*,'sqf=',sqf
492 c-------------
493
494 RETURN
495 END
496 C_________________________________________________________________
497 C_________________________________________________________________
498