| GCC Code Coverage Report | |||||||||||||||||||||
|
|||||||||||||||||||||
| Line | Branch | Exec | Source |
1 |
! $Id: advtrac.F90 4143 2022-05-09 10:35:40Z dcugnet $ |
||
2 |
|||
3 |
#define DEBUG_IO |
||
4 |
#undef DEBUG_IO |
||
5 |
1558 |
SUBROUTINE advtrac(pbaru, pbarv, p, masse,q,iapptrac,teta, flxw, pk) |
|
6 |
! Auteur : F. Hourdin |
||
7 |
! |
||
8 |
! Modif. P. Le Van (20/12/97) |
||
9 |
! F. Codron (10/99) |
||
10 |
! D. Le Croller (07/2001) |
||
11 |
! M.A Filiberti (04/2002) |
||
12 |
! |
||
13 |
USE infotrac, ONLY: nqtot, tracers, isoCheck |
||
14 |
USE control_mod, ONLY: iapp_tracvl, day_step |
||
15 |
USE comconst_mod, ONLY: dtvr |
||
16 |
|||
17 |
IMPLICIT NONE |
||
18 |
! |
||
19 |
include "dimensions.h" |
||
20 |
include "paramet.h" |
||
21 |
include "comdissip.h" |
||
22 |
include "comgeom2.h" |
||
23 |
include "description.h" |
||
24 |
include "iniprint.h" |
||
25 |
|||
26 |
!--------------------------------------------------------------------------- |
||
27 |
! Arguments |
||
28 |
!--------------------------------------------------------------------------- |
||
29 |
INTEGER, INTENT(OUT) :: iapptrac |
||
30 |
REAL, INTENT(IN) :: pbaru(ip1jmp1,llm) |
||
31 |
REAL, INTENT(IN) :: pbarv(ip1jm, llm) |
||
32 |
REAL, INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) |
||
33 |
REAL, INTENT(IN) :: masse(ip1jmp1,llm) |
||
34 |
REAL, INTENT(IN) :: p(ip1jmp1,llmp1 ) |
||
35 |
REAL, INTENT(IN) :: teta(ip1jmp1,llm) |
||
36 |
REAL, INTENT(IN) :: pk(ip1jmp1,llm) |
||
37 |
REAL, INTENT(OUT) :: flxw(ip1jmp1,llm) |
||
38 |
!--------------------------------------------------------------------------- |
||
39 |
! Ajout PPM |
||
40 |
!--------------------------------------------------------------------------- |
||
41 |
REAL :: massebx(ip1jmp1,llm), masseby(ip1jm,llm) |
||
42 |
!--------------------------------------------------------------------------- |
||
43 |
! Variables locales |
||
44 |
!--------------------------------------------------------------------------- |
||
45 |
INTEGER :: ij, l, iq, iadv |
||
46 |
! REAL(KIND=KIND(1.d0)) :: t_initial, t_final, tps_cpu |
||
47 |
REAL :: zdp(ip1jmp1), zdpmin, zdpmax |
||
48 |
INTEGER, SAVE :: iadvtr=0 |
||
49 |
REAL, DIMENSION(ip1jmp1,llm) :: pbaruc, pbarug, massem, wg |
||
50 |
REAL, DIMENSION(ip1jm, llm) :: pbarvc, pbarvg |
||
51 |
EXTERNAL minmax |
||
52 |
SAVE massem, pbaruc, pbarvc |
||
53 |
!--------------------------------------------------------------------------- |
||
54 |
! Rajouts pour PPM |
||
55 |
!--------------------------------------------------------------------------- |
||
56 |
INTEGER indice, n |
||
57 |
REAL :: dtbon ! Pas de temps adaptatif pour que CFL<1 |
||
58 |
REAL :: CFLmaxz, aaa, bbb ! CFL maximum |
||
59 |
REAL, DIMENSION(iim,jjp1,llm) :: unatppm, vnatppm, fluxwppm |
||
60 |
1153 |
REAL :: qppm(iim*jjp1,llm,nqtot) |
|
61 |
REAL :: psppm(iim,jjp1) ! pression au sol |
||
62 |
REAL, DIMENSION(llmp1) :: apppm, bpppm |
||
63 |
LOGICAL, SAVE :: dum=.TRUE., fill=.TRUE. |
||
64 |
|||
65 |
INTEGER, SAVE :: countcfl=0 |
||
66 |
REAL, DIMENSION(ip1jmp1,llm) :: cflx, cflz |
||
67 |
REAL, DIMENSION(ip1jm ,llm) :: cfly |
||
68 |
REAL, DIMENSION(llm), SAVE :: cflxmax, cflymax, cflzmax |
||
69 |
|||
70 |
✓✓ | 1441 |
IF(iadvtr == 0) THEN |
71 |
289 |
pbaruc(:,:)=0 |
|
72 |
289 |
pbarvc(:,:)=0 |
|
73 |
END IF |
||
74 |
|||
75 |
!--- Accumulation des flux de masse horizontaux |
||
76 |
✓✓ | 57640 |
DO l=1,llm |
77 |
✓✓ | 61256910 |
DO ij = 1,ip1jmp1 |
78 |
61256910 |
pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l) |
|
79 |
END DO |
||
80 |
✓✓ | 59403784 |
DO ij = 1,ip1jm |
81 |
59402343 |
pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l) |
|
82 |
END DO |
||
83 |
END DO |
||
84 |
|||
85 |
!--- Selection de la masse instantannee des mailles avant le transport. |
||
86 |
✓✓ | 1441 |
IF(iadvtr == 0) THEN |
87 |
289 |
CALL SCOPY(ip1jmp1*llm,masse,1,massem,1) |
|
88 |
! CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 ) |
||
89 |
END IF |
||
90 |
|||
91 |
1441 |
iadvtr = iadvtr+1 |
|
92 |
1441 |
iapptrac = iadvtr |
|
93 |
|||
94 |
!--- Test pour savoir si on advecte a ce pas de temps |
||
95 |
✓✓ | 1441 |
IF(iadvtr /= iapp_tracvl) RETURN |
96 |
|||
97 |
! .. Modif P.Le Van ( 20/12/97 ) .... |
||
98 |
! |
||
99 |
! traitement des flux de masse avant advection. |
||
100 |
! 1. calcul de w |
||
101 |
! 2. groupement des mailles pres du pole. |
||
102 |
|||
103 |
288 |
CALL groupe(massem, pbaruc, pbarvc, pbarug, pbarvg, wg) |
|
104 |
|||
105 |
!--- Flux de masse diaganostiques traceurs |
||
106 |
✓✓✓✓ |
12243168 |
flxw = wg / REAL(iapp_tracvl) |
107 |
|||
108 |
!--- Test sur l'eventuelle creation de valeurs negatives de la masse |
||
109 |
✓✓ | 11232 |
DO l=1,llm-1 |
110 |
✓✓ | 11195712 |
DO ij = iip2+1,ip1jm |
111 |
zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l) & |
||
112 |
- pbarvg(ij-iip1,l) + pbarvg(ij,l) & |
||
113 |
11195712 |
+ wg(ij,l+1) - wg(ij,l) |
|
114 |
END DO |
||
115 |
! ym ---> pourquoi jjm-1 et non jjm ? a cause du pole ? |
||
116 |
10944 |
CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 ) |
|
117 |
✓✓ | 11206656 |
DO ij = iip2,ip1jm |
118 |
11206656 |
zdp(ij)= zdp(ij)*dtvr/ massem(ij,l) |
|
119 |
END DO |
||
120 |
|||
121 |
10944 |
CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax ) |
|
122 |
|||
123 |
✗✓ | 10944 |
IF(MAX(ABS(zdpmin),ABS(zdpmax)) > 0.5) & |
124 |
288 |
WRITE(*,*)'WARNING DP/P l=',l,' MIN:',zdpmin,' MAX:', zdpmax |
|
125 |
|||
126 |
END DO |
||
127 |
|||
128 |
!------------------------------------------------------------------------- |
||
129 |
! Calcul des criteres CFL en X, Y et Z |
||
130 |
!------------------------------------------------------------------------- |
||
131 |
✓✓ | 288 |
IF(countcfl == 0. ) then |
132 |
3 |
cflxmax(:)=0. |
|
133 |
3 |
cflymax(:)=0. |
|
134 |
3 |
cflzmax(:)=0. |
|
135 |
END IF |
||
136 |
|||
137 |
288 |
countcfl=countcfl+iapp_tracvl |
|
138 |
288 |
cflx(:,:)=0. |
|
139 |
288 |
cfly(:,:)=0. |
|
140 |
288 |
cflz(:,:)=0. |
|
141 |
✓✓ | 11520 |
DO l=1,llm |
142 |
✓✓ | 11490624 |
DO ij=iip2,ip1jm-1 |
143 |
✓✓ | 11490336 |
IF(pbarug(ij,l)>=0.) then |
144 |
8186556 |
cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l) |
|
145 |
ELSE |
||
146 |
3292548 |
cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l) |
|
147 |
END IF |
||
148 |
END DO |
||
149 |
END DO |
||
150 |
|||
151 |
✓✓ | 11520 |
DO l=1,llm |
152 |
288 |
DO ij=iip2,ip1jm-1,iip1 |
|
153 |
✓✓ | 348192 |
cflx(ij+iip1,l)=cflx(ij,l) |
154 |
END DO |
||
155 |
END DO |
||
156 |
|||
157 |
✓✓ | 11520 |
DO l=1,llm |
158 |
✓✓ | 11872512 |
DO ij=1,ip1jm |
159 |
✓✓ | 11872224 |
IF(pbarvg(ij,l)>=0.) then |
160 |
5972213 |
cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l) |
|
161 |
ELSE |
||
162 |
5888779 |
cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l) |
|
163 |
END IF |
||
164 |
END DO |
||
165 |
END DO |
||
166 |
|||
167 |
✓✓ | 11232 |
DO l=2,llm |
168 |
✓✓ | 11568096 |
DO ij=1,ip1jm |
169 |
✓✓ | 11567808 |
IF(wg(ij,l) >= 0.) THEN |
170 |
5844813 |
cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l) |
|
171 |
ELSE |
||
172 |
5712051 |
cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1) |
|
173 |
END IF |
||
174 |
END DO |
||
175 |
END DO |
||
176 |
|||
177 |
✓✓ | 11520 |
DO l=1,llm |
178 |
✗✓✓✗ ✓✓✓✓ |
12254112 |
cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l))) |
179 |
✗✓✓✗ ✓✓✓✓ |
11883456 |
cflymax(l)=max(cflymax(l),maxval(cfly(:,l))) |
180 |
✗✓✓✗ ✓✓✓✓ |
12254400 |
cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l))) |
181 |
END DO |
||
182 |
|||
183 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
||
184 |
! Par defaut, on sort le diagnostic des CFL tous les jours. |
||
185 |
! Si on veut le sortir a chaque pas d'advection en cas de plantage |
||
186 |
! IF(countcfl==iapp_tracvl) then |
||
187 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
||
188 |
✓✓ | 288 |
IF(countcfl==day_step) then |
189 |
✓✓ | 120 |
DO l=1,llm |
190 |
120 |
WRITE(lunout,*) 'L, CFL[xyz]max:', l, cflxmax(l), cflymax(l), cflzmax(l) |
|
191 |
END DO |
||
192 |
3 |
countcfl=0 |
|
193 |
END IF |
||
194 |
|||
195 |
!--------------------------------------------------------------------------- |
||
196 |
! Advection proprement dite (Modification Le Croller (07/2001) |
||
197 |
!--------------------------------------------------------------------------- |
||
198 |
|||
199 |
!--------------------------------------------------------------------------- |
||
200 |
! Calcul des moyennes basees sur la masse |
||
201 |
!--------------------------------------------------------------------------- |
||
202 |
288 |
CALL massbar(massem,massebx,masseby) |
|
203 |
|||
204 |
#ifdef DEBUG_IO |
||
205 |
CALL WriteField_u('massem',massem) |
||
206 |
CALL WriteField_u('wg',wg) |
||
207 |
CALL WriteField_u('pbarug',pbarug) |
||
208 |
CALL WriteField_v('pbarvg',pbarvg) |
||
209 |
CALL WriteField_u('p_tmp',p) |
||
210 |
CALL WriteField_u('pk_tmp',pk) |
||
211 |
CALL WriteField_u('teta_tmp',teta) |
||
212 |
DO iq=1,nqtot |
||
213 |
CALL WriteField_u('q_adv'//trim(int2str(iq)),q(:,:,iq)) |
||
214 |
END DO |
||
215 |
#endif |
||
216 |
|||
217 |
✗✓ | 288 |
IF(isoCheck) WRITE(*,*) 'advtrac 227' |
218 |
288 |
CALL check_isotopes_seq(q,ip1jmp1,'advtrac 162') |
|
219 |
|||
220 |
!------------------------------------------------------------------------- |
||
221 |
! Appel des sous programmes d'advection |
||
222 |
!------------------------------------------------------------------------- |
||
223 |
✓✓ | 1728 |
DO iq = 1, nqtot |
224 |
! CALL clock(t_initial) |
||
225 |
✓✗ | 1440 |
IF(tracers(iq)%parent /= 'air') CYCLE |
226 |
1440 |
iadv = tracers(iq)%iadv |
|
227 |
!----------------------------------------------------------------------- |
||
228 |
288 |
SELECT CASE(iadv) |
|
229 |
!----------------------------------------------------------------------- |
||
230 |
1152 |
CASE(0); CYCLE |
|
231 |
!-------------------------------------------------------------------- |
||
232 |
CASE(10) !--- Schema de Van Leer I MUSCL |
||
233 |
!-------------------------------------------------------------------- |
||
234 |
! WRITE(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:) |
||
235 |
1152 |
CALL vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq) |
|
236 |
|||
237 |
!-------------------------------------------------------------------- |
||
238 |
CASE(14) !--- Schema "pseuDO amont" + test sur humidite specifique |
||
239 |
!--- pour la vapeur d'eau. F. Codron |
||
240 |
!-------------------------------------------------------------------- |
||
241 |
! WRITE(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:) |
||
242 |
288 |
CALL vlspltqs(q,2.,massem,wg,pbarug,pbarvg,dtvr,p,pk,teta,iq) |
|
243 |
|||
244 |
!-------------------------------------------------------------------- |
||
245 |
CASE(12) !--- Schema de Frederic Hourdin |
||
246 |
!-------------------------------------------------------------------- |
||
247 |
CALL adaptdt(iadv,dtbon,n,pbarug,massem) ! pas de temps adaptatif |
||
248 |
IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n |
||
249 |
DO indice=1,n |
||
250 |
CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1) |
||
251 |
END DO |
||
252 |
|||
253 |
!-------------------------------------------------------------------- |
||
254 |
CASE(13) !--- Pas de temps adaptatif |
||
255 |
!-------------------------------------------------------------------- |
||
256 |
CALL adaptdt(iadv,dtbon,n,pbarug,massem) |
||
257 |
IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n |
||
258 |
DO indice=1,n |
||
259 |
CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2) |
||
260 |
END DO |
||
261 |
|||
262 |
!-------------------------------------------------------------------- |
||
263 |
CASE(20) !--- Schema de pente SLOPES |
||
264 |
!-------------------------------------------------------------------- |
||
265 |
CALL pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0) |
||
266 |
|||
267 |
!-------------------------------------------------------------------- |
||
268 |
CASE(30) !--- Schema de Prather |
||
269 |
!-------------------------------------------------------------------- |
||
270 |
! Pas de temps adaptatif |
||
271 |
CALL adaptdt(iadv,dtbon,n,pbarug,massem) |
||
272 |
IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n |
||
273 |
CALL prather(q(1,1,iq),wg,massem,pbarug,pbarvg,n,dtbon) |
||
274 |
|||
275 |
!-------------------------------------------------------------------- |
||
276 |
CASE(11,16,17,18) !--- Schemas PPM Lin et Rood |
||
277 |
!-------------------------------------------------------------------- |
||
278 |
! Test sur le flux horizontal |
||
279 |
CALL adaptdt(iadv,dtbon,n,pbarug,massem) ! pas de temps adaptatif |
||
280 |
IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n |
||
281 |
! Test sur le flux vertical |
||
282 |
CFLmaxz=0. |
||
283 |
DO l=2,llm |
||
284 |
DO ij=iip2,ip1jm |
||
285 |
aaa=wg(ij,l)*dtvr/massem(ij,l) |
||
286 |
CFLmaxz=max(CFLmaxz,aaa) |
||
287 |
bbb=-wg(ij,l)*dtvr/massem(ij,l-1) |
||
288 |
CFLmaxz=max(CFLmaxz,bbb) |
||
289 |
END DO |
||
290 |
END DO |
||
291 |
IF(CFLmaxz.GE.1) WRITE(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz |
||
292 |
!---------------------------------------------------------------- |
||
293 |
! Ss-prg interface LMDZ.4->PPM3d (ss-prg de Lin) |
||
294 |
!---------------------------------------------------------------- |
||
295 |
CALL interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, & |
||
296 |
apppm,bpppm,massebx,masseby,pbarug,pbarvg, & |
||
297 |
unatppm,vnatppm,psppm) |
||
298 |
|||
299 |
!---------------------------------------------------------------- |
||
300 |
DO indice=1,n !--- VL (version PPM) horiz. et PPM vert. |
||
301 |
!---------------------------------------------------------------- |
||
302 |
SELECT CASE(iadv) |
||
303 |
!---------------------------------------------------------- |
||
304 |
CASE(11) |
||
305 |
!---------------------------------------------------------- |
||
306 |
CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, & |
||
307 |
2,2,2,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.) |
||
308 |
!---------------------------------------------------------- |
||
309 |
CASE(16) !--- Monotonic PPM |
||
310 |
!---------------------------------------------------------- |
||
311 |
CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, & |
||
312 |
3,3,3,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.) |
||
313 |
!---------------------------------------------------------- |
||
314 |
CASE(17) !--- Semi monotonic PPM |
||
315 |
!---------------------------------------------------------- |
||
316 |
CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, & |
||
317 |
4,4,4,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, fill,dum,220.) |
||
318 |
!---------------------------------------------------------- |
||
319 |
CASE(18) !--- Positive Definite PPM |
||
320 |
!---------------------------------------------------------- |
||
321 |
CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, & |
||
322 |
5,5,5,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.) |
||
323 |
END SELECT |
||
324 |
!---------------------------------------------------------------- |
||
325 |
END DO |
||
326 |
!---------------------------------------------------------------- |
||
327 |
! Ss-prg interface PPM3d-LMDZ.4 |
||
328 |
!---------------------------------------------------------------- |
||
329 |
✓✓✗✗ ✗✗✗✗ |
1440 |
CALL interpost(q(1,1,iq),qppm(1,1,iq)) |
330 |
!---------------------------------------------------------------------- |
||
331 |
END SELECT |
||
332 |
!---------------------------------------------------------------------- |
||
333 |
|||
334 |
!---------------------------------------------------------------------- |
||
335 |
! On impose une seule valeur du traceur au pole Sud j=jjm+1=jjp1 et Nord j=1 |
||
336 |
!---------------------------------------------------------------------- |
||
337 |
! CALL traceurpole(q(1,1,iq),massem) |
||
338 |
|||
339 |
!--- Calcul du temps cpu pour un schema donne |
||
340 |
! CALL clock(t_final) |
||
341 |
!ym tps_cpu=t_final-t_initial |
||
342 |
!ym cpuadv(iq)=cpuadv(iq)+tps_cpu |
||
343 |
|||
344 |
END DO |
||
345 |
|||
346 |
✗✓ | 288 |
IF(isoCheck) WRITE(*,*) 'advtrac 402' |
347 |
288 |
CALL check_isotopes_seq(q,ip1jmp1,'advtrac 397') |
|
348 |
|||
349 |
!------------------------------------------------------------------------- |
||
350 |
! on reinitialise a zero les flux de masse cumules |
||
351 |
!------------------------------------------------------------------------- |
||
352 |
288 |
iadvtr=0 |
|
353 |
|||
354 |
END SUBROUTINE advtrac |
| Generated by: GCOVR (Version 4.2) |