GCC Code Coverage Report

Directory: ./
File: dyn/advtrac.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 71 124 57.3%
Branches: 79 128 61.7%
Line Branch Exec Source
1 ! $Id: advtrac.F90 2622 2016-09-04 06:12:02Z emillour $
2
3 2596 SUBROUTINE advtrac(pbaru,pbarv , p, masse,q,iapptrac,teta, flxw, pk)
4 ! Auteur : F. Hourdin
5 !
6 ! Modif. P. Le Van (20/12/97)
7 ! F. Codron (10/99)
8 ! D. Le Croller (07/2001)
9 ! M.A Filiberti (04/2002)
10 !
11 USE infotrac, ONLY: nqtot, iadv,nqperes,ok_iso_verif
12 USE control_mod, ONLY: iapp_tracvl, day_step
13 USE comconst_mod, ONLY: dtvr
14
15 IMPLICIT NONE
16 !
17 include "dimensions.h"
18 include "paramet.h"
19 include "comdissip.h"
20 include "comgeom2.h"
21 include "description.h"
22 include "iniprint.h"
23
24 !-------------------------------------------------------------------
25 ! Arguments
26 !-------------------------------------------------------------------
27 INTEGER,INTENT(OUT) :: iapptrac
28 REAL,INTENT(IN) :: pbaru(ip1jmp1,llm)
29 REAL,INTENT(IN) :: pbarv(ip1jm,llm)
30 REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot)
31 REAL,INTENT(IN) :: masse(ip1jmp1,llm)
32 REAL,INTENT(IN) :: p( ip1jmp1,llmp1 )
33 REAL,INTENT(IN) :: teta(ip1jmp1,llm)
34 REAL,INTENT(IN) :: pk(ip1jmp1,llm)
35 REAL,INTENT(OUT) :: flxw(ip1jmp1,llm)
36 !-------------------------------------------------------------------
37 ! Ajout PPM
38 !--------------------------------------------------------
39 REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
40 !-------------------------------------------------------------
41 ! Variables locales
42 !-------------------------------------------------------------
43
44 REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
45 REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
46 REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
47 REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
48 INTEGER iadvtr
49 INTEGER ij,l,iq,iiq
50 REAL zdpmin, zdpmax
51 EXTERNAL minmax
52 SAVE iadvtr, massem, pbaruc, pbarvc
53 DATA iadvtr/0/
54 !----------------------------------------------------------
55 ! Rajouts pour PPM
56 !----------------------------------------------------------
57 INTEGER indice,n
58 REAL dtbon ! Pas de temps adaptatif pour que CFL<1
59 REAL CFLmaxz,aaa,bbb ! CFL maximum
60 REAL psppm(iim,jjp1) ! pression au sol
61 REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
62 4802 REAL qppm(iim*jjp1,llm,nqtot)
63 REAL fluxwppm(iim,jjp1,llm)
64 REAL apppm(llmp1), bpppm(llmp1)
65 LOGICAL dum,fill
66 DATA fill/.true./
67 DATA dum/.true./
68
69 integer,save :: countcfl=0
70 real cflx(ip1jmp1,llm)
71 real cfly(ip1jm,llm)
72 real cflz(ip1jmp1,llm)
73 real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm)
74
75
2/2
✓ Branch 0 taken 481 times.
✓ Branch 1 taken 1920 times.
 2401 IF(iadvtr.EQ.0) THEN
76 481 pbaruc(:,:)=0
77 481 pbarvc(:,:)=0
78 ENDIF
79
80 ! accumulation des flux de masse horizontaux
81
2/2
✓ Branch 0 taken 93639 times.
✓ Branch 1 taken 2401 times.
 96040 DO l=1,llm
82
2/2
✓ Branch 0 taken 93639 times.
✓ Branch 1 taken 101972871 times.
 102066510 DO ij = 1,ip1jmp1
83 102066510 pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
84 ENDDO
85
2/2
✓ Branch 0 taken 98882784 times.
✓ Branch 1 taken 93639 times.
 98978824 DO ij = 1,ip1jm
86 98976423 pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
87 ENDDO
88 ENDDO
89
90 ! selection de la masse instantannee des mailles avant le transport.
91
2/2
✓ Branch 0 taken 481 times.
✓ Branch 1 taken 1920 times.
 2401 IF(iadvtr.EQ.0) THEN
92
93 481 CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
94 !cc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
95 !
96 ENDIF
97
98 2401 iadvtr = iadvtr+1
99 2401 iapptrac = iadvtr
100
101
102 ! Test pour savoir si on advecte a ce pas de temps
103
2/2
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 1921 times.
 2401 IF ( iadvtr.EQ.iapp_tracvl ) THEN
104
105 !c .. Modif P.Le Van ( 20/12/97 ) ....
106 !c
107
108 ! traitement des flux de masse avant advection.
109 ! 1. calcul de w
110 ! 2. groupement des mailles pres du pole.
111
112 480 CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
113
114 ! ... Flux de masse diaganostiques traceurs
115
4/4
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 20386080 times.
✓ Branch 3 taken 18720 times.
 20405280 flxw = wg / REAL(iapp_tracvl)
116
117 ! test sur l'eventuelle creation de valeurs negatives de la masse
118
2/2
✓ Branch 0 taken 18240 times.
✓ Branch 1 taken 480 times.
 18720 DO l=1,llm-1
119
2/2
✓ Branch 0 taken 18641280 times.
✓ Branch 1 taken 18240 times.
 18659520 DO ij = iip2+1,ip1jm
120 zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l) &
121 - pbarvg(ij-iip1,l) + pbarvg(ij,l) &
122 18659520 + wg(ij,l+1) - wg(ij,l)
123 ENDDO
124 18240 CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
125
2/2
✓ Branch 0 taken 18659520 times.
✓ Branch 1 taken 18240 times.
 18677760 DO ij = iip2,ip1jm
126 18677760 zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
127 ENDDO
128
129
130 18240 CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
131
132
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18240 times.
 18720 IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
133 ✗ PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin, &
134 ✗ ' MAX:', zdpmax
135 ENDIF
136
137 ENDDO
138
139
140 !-------------------------------------------------------------------
141 ! Calcul des criteres CFL en X, Y et Z
142 !-------------------------------------------------------------------
143
144
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 475 times.
 480 if (countcfl == 0. ) then
145 5 cflxmax(:)=0.
146 5 cflymax(:)=0.
147 5 cflzmax(:)=0.
148 endif
149
150 480 countcfl=countcfl+iapp_tracvl
151 480 cflx(:,:)=0.
152 480 cfly(:,:)=0.
153 480 cflz(:,:)=0.
154
2/2
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
 19200 do l=1,llm
155
2/2
✓ Branch 0 taken 19131840 times.
✓ Branch 1 taken 18720 times.
 19151040 do ij=iip2,ip1jm-1
156
2/2
✓ Branch 0 taken 13641439 times.
✓ Branch 1 taken 5490401 times.
 19150560 if (pbarug(ij,l)>=0.) then
157 13641439 cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l)
158 else
159 5490401 cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l)
160 endif
161 enddo
162 enddo
163
2/2
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
 19200 do l=1,llm
164 480 do ij=iip2,ip1jm-1,iip1
165
2/2
✓ Branch 0 taken 561600 times.
✓ Branch 1 taken 18720 times.
 580320 cflx(ij+iip1,l)=cflx(ij,l)
166 enddo
167 enddo
168
169
2/2
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
 19200 do l=1,llm
170
2/2
✓ Branch 0 taken 19768320 times.
✓ Branch 1 taken 18720 times.
 19787520 do ij=1,ip1jm
171
2/2
✓ Branch 0 taken 10029057 times.
✓ Branch 1 taken 9739263 times.
 19787040 if (pbarvg(ij,l)>=0.) then
172 10029057 cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l)
173 else
174 9739263 cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l)
175 endif
176 enddo
177 enddo
178
179
2/2
✓ Branch 0 taken 18240 times.
✓ Branch 1 taken 480 times.
 18720 do l=2,llm
180
2/2
✓ Branch 0 taken 19261440 times.
✓ Branch 1 taken 18240 times.
 19280160 do ij=1,ip1jm
181
2/2
✓ Branch 0 taken 9681935 times.
✓ Branch 1 taken 9579505 times.
 19279680 if (wg(ij,l)>=0.) then
182 9681935 cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l)
183 else
184 9579505 cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1)
185 endif
186 enddo
187 enddo
188
189
2/2
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
 19200 do l=1,llm
190
6/8
✗ Branch 0 not taken.
✓ Branch 1 taken 18720 times.
✓ Branch 2 taken 18720 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 20386080 times.
✓ Branch 5 taken 18720 times.
✓ Branch 6 taken 237236 times.
✓ Branch 7 taken 20148844 times.
 20423520 cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l)))
191
6/8
✗ Branch 0 not taken.
✓ Branch 1 taken 18720 times.
✓ Branch 2 taken 18720 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 19768320 times.
✓ Branch 5 taken 18720 times.
✓ Branch 6 taken 82538 times.
✓ Branch 7 taken 19685782 times.
 19805760 cflymax(l)=max(cflymax(l),maxval(cfly(:,l)))
192
6/8
✗ Branch 0 not taken.
✓ Branch 1 taken 18720 times.
✓ Branch 2 taken 18720 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 20386080 times.
✓ Branch 5 taken 18720 times.
✓ Branch 6 taken 294638 times.
✓ Branch 7 taken 20091442 times.
 20424000 cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l)))
193 enddo
194
195 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
196 ! Par defaut, on sort le diagnostic des CFL tous les jours.
197 ! Si on veut le sortir a chaque pas d'advection en cas de plantage
198 ! if (countcfl==iapp_tracvl) then
199 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
200
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 475 times.
 480 if (countcfl==day_step) then
201
2/2
✓ Branch 0 taken 195 times.
✓ Branch 1 taken 5 times.
 200 do l=1,llm
202 195 write(lunout,*) 'L, CFL[xyz]max:', l, cflxmax(l), cflymax(l), &
203 395 cflzmax(l)
204 enddo
205 5 countcfl=0
206 endif
207
208 !-------------------------------------------------------------------
209 ! Advection proprement dite (Modification Le Croller (07/2001)
210 !-------------------------------------------------------------------
211
212 !----------------------------------------------------
213 ! Calcul des moyennes bas�es sur la masse
214 !----------------------------------------------------
215 480 call massbar(massem,massebx,masseby)
216
217 !-----------------------------------------------------------
218 ! Appel des sous programmes d'advection
219 !-----------------------------------------------------------
220
221
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 480 times.
 480 if (ok_iso_verif) then
222 ✗ write(*,*) 'advtrac 227'
223 ✗ call check_isotopes_seq(q,ip1jmp1,'advtrac 162')
224 endif !if (ok_iso_verif) then
225
226
2/2
✓ Branch 0 taken 2400 times.
✓ Branch 1 taken 480 times.
 2880 do iq=1,nqperes
227 ! call clock(t_initial)
228
1/2
✓ Branch 0 taken 2400 times.
✗ Branch 1 not taken.
 2400 if(iadv(iq) == 0) cycle
229 ! ----------------------------------------------------------------
230 ! Schema de Van Leer I MUSCL
231 ! ----------------------------------------------------------------
232
2/2
✓ Branch 0 taken 1920 times.
✓ Branch 1 taken 480 times.
 2880 if(iadv(iq).eq.10) THEN
233 ! CRisi: on fait passer tout q pour avoir acces aux fils
234
235 !write(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)
236 1920 call vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq)
237
238 ! ----------------------------------------------------------------
239 ! Schema "pseudo amont" + test sur humidite specifique
240 ! pour la vapeur d'eau. F. Codron
241 ! ----------------------------------------------------------------
242
1/2
✓ Branch 0 taken 480 times.
✗ Branch 1 not taken.
 480 else if(iadv(iq).eq.14) then
243 !
244 !write(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
245 CALL vlspltqs( q, 2., massem, wg , &
246 480 pbarug,pbarvg,dtvr,p,pk,teta,iq)
247
248 ! ----------------------------------------------------------------
249 ! Schema de Frederic Hourdin
250 ! ----------------------------------------------------------------
251 ✗ else if(iadv(iq).eq.12) then
252 ! Pas de temps adaptatif
253 ✗ call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
254 ✗ if (n.GT.1) then
255 ✗ write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
256 ✗ dtvr,'n=',n
257 endif
258 ✗ do indice=1,n
259 ✗ call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
260 end do
261 ✗ else if(iadv(iq).eq.13) then
262 ! Pas de temps adaptatif
263 ✗ call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
264 ✗ if (n.GT.1) then
265 ✗ write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
266 ✗ dtvr,'n=',n
267 endif
268 ✗ do indice=1,n
269 ✗ call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
270 end do
271 ! ----------------------------------------------------------------
272 ! Schema de pente SLOPES
273 ! ----------------------------------------------------------------
274 ✗ else if (iadv(iq).eq.20) then
275 ✗ call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
276
277 ! ----------------------------------------------------------------
278 ! Schema de Prather
279 ! ----------------------------------------------------------------
280 ✗ else if (iadv(iq).eq.30) then
281 ! Pas de temps adaptatif
282 ✗ call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
283 ✗ if (n.GT.1) then
284 ✗ write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
285 ✗ dtvr,'n=',n
286 endif
287 call prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
288 ✗ n,dtbon)
289
290 ! ----------------------------------------------------------------
291 ! Schemas PPM Lin et Rood
292 ! ----------------------------------------------------------------
293 ✗ else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
294 iadv(iq).LE.18)) then
295
296 ! Test sur le flux horizontal
297 ! Pas de temps adaptatif
298 ✗ call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
299 ✗ if (n.GT.1) then
300 ✗ write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
301 ✗ dtvr,'n=',n
302 endif
303 ! Test sur le flux vertical
304 ✗ CFLmaxz=0.
305 ✗ do l=2,llm
306 ✗ do ij=iip2,ip1jm
307 ✗ aaa=wg(ij,l)*dtvr/massem(ij,l)
308 ✗ CFLmaxz=max(CFLmaxz,aaa)
309 ✗ bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
310 ✗ CFLmaxz=max(CFLmaxz,bbb)
311 enddo
312 enddo
313 ✗ if (CFLmaxz.GE.1) then
314 ✗ write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
315 endif
316
317 !-----------------------------------------------------------
318 ! Ss-prg interface LMDZ.4->PPM3d
319 !-----------------------------------------------------------
320
321 call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
322 apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
323 ✗ unatppm,vnatppm,psppm)
324
325 ✗ do indice=1,n
326 !----------------------------------------------------------------
327 ! VL (version PPM) horiz. et PPM vert.
328 !----------------------------------------------------------------
329 ✗ if (iadv(iq).eq.11) then
330 ! Ss-prg PPM3d de Lin
331 call ppm3d(1,qppm(1,1,iq), &
332 psppm,psppm, &
333 unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
334 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
335 ✗ fill,dum,220.)
336
337 !-------------------------------------------------------------
338 ! Monotonic PPM
339 !-------------------------------------------------------------
340 ✗ else if (iadv(iq).eq.16) then
341 ! Ss-prg PPM3d de Lin
342 call ppm3d(1,qppm(1,1,iq), &
343 psppm,psppm, &
344 unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
345 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
346 ✗ fill,dum,220.)
347 !-------------------------------------------------------------
348
349 !-------------------------------------------------------------
350 ! Semi Monotonic PPM
351 !-------------------------------------------------------------
352 ✗ else if (iadv(iq).eq.17) then
353 ! Ss-prg PPM3d de Lin
354 call ppm3d(1,qppm(1,1,iq), &
355 psppm,psppm, &
356 unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
357 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
358 ✗ fill,dum,220.)
359 !-------------------------------------------------------------
360
361 !-------------------------------------------------------------
362 ! Positive Definite PPM
363 !-------------------------------------------------------------
364 ✗ else if (iadv(iq).eq.18) then
365 ! Ss-prg PPM3d de Lin
366 call ppm3d(1,qppm(1,1,iq), &
367 psppm,psppm, &
368 unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
369 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
370 ✗ fill,dum,220.)
371 !-------------------------------------------------------------
372 endif
373 enddo
374 !-----------------------------------------------------------------
375 ! Ss-prg interface PPM3d-LMDZ.4
376 !-----------------------------------------------------------------
377 ✗ call interpost(q(1,1,iq),qppm(1,1,iq))
378 endif
379 !----------------------------------------------------------------------
380
381 !-----------------------------------------------------------------
382 ! On impose une seule valeur du traceur au p�le Sud j=jjm+1=jjp1
383 ! et Nord j=1
384 !-----------------------------------------------------------------
385
386 ! call traceurpole(q(1,1,iq),massem)
387
388 ! calcul du temps cpu pour un schema donne
389
390 ! call clock(t_final)
391 !ym tps_cpu=t_final-t_initial
392 !ym cpuadv(iq)=cpuadv(iq)+tps_cpu
393
394 end DO
395
396
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 480 times.
 480 if (ok_iso_verif) then
397 ✗ write(*,*) 'advtrac 402'
398 ✗ call check_isotopes_seq(q,ip1jmp1,'advtrac 397')
399 endif !if (ok_iso_verif) then
400
401 !------------------------------------------------------------------
402 ! on reinitialise a zero les flux de masse cumules
403 !---------------------------------------------------
404 480 iadvtr=0
405
406 ENDIF ! if iadvtr.EQ.iapp_tracvl
407
408 2401 END SUBROUTINE advtrac
409