GCC Code Coverage Report

Directory: ./
File: dyn/vlspltqs.f
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 179 259 69.1%
Branches: 127 204 62.3%
Line Branch Exec Source
1 c
2 c $Id: vlspltqs.F 2603 2016-07-25 09:31:56Z emillour $
3 c
4 480 SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt,
5 , p,pk,teta,iq )
6 USE infotrac, ONLY: nqtot,nqdesc,iqfils
7 c
8 c Auteurs: P.Le Van, F.Hourdin, F.Forget, F.Codron
9 c
10 c ********************************************************************
11 c Shema d'advection " pseudo amont " .
12 c + test sur humidite specifique: Q advecte< Qsat aval
13 c (F. Codron, 10/99)
14 c ********************************************************************
15 c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
16 c
17 c pente_max facteur de limitation des pentes: 2 en general
18 c 0 pour un schema amont
19 c pbaru,pbarv,w flux de masse en u ,v ,w
20 c pdt pas de temps
21 c
22 c teta temperature potentielle, p pression aux interfaces,
23 c pk exner au milieu des couches necessaire pour calculer Qsat
24 c --------------------------------------------------------------------
25
26 USE comconst_mod, ONLY: cpp
27
28 IMPLICIT NONE
29 c
30 include "dimensions.h"
31 include "paramet.h"
32
33 c
34 c Arguments:
35 c ----------
36 REAL masse(ip1jmp1,llm),pente_max
37 REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
38 REAL q(ip1jmp1,llm,nqtot)
39 REAL w(ip1jmp1,llm),pdt
40 REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm)
41 INTEGER iq ! CRisi
42 c
43 c Local
44 c ---------
45 c
46 INTEGER i,ij,l,j,ii
47 INTEGER ifils,iq2 ! CRisi
48 c
49 REAL qsat(ip1jmp1,llm)
50 960 REAL zm(ip1jmp1,llm,nqtot)
51 REAL mu(ip1jmp1,llm)
52 REAL mv(ip1jm,llm)
53 REAL mw(ip1jmp1,llm+1)
54 480 REAL zq(ip1jmp1,llm,nqtot)
55 REAL temps1,temps2,temps3
56 REAL zzpbar, zzw
57 LOGICAL testcpu
58 SAVE testcpu
59 SAVE temps1,temps2,temps3
60
61 REAL qmin,qmax
62 DATA qmin,qmax/0.,1.e33/
63 DATA testcpu/.false./
64 DATA temps1,temps2,temps3/0.,0.,0./
65
66 c--pour rapport de melange saturant--
67
68 REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play
69 REAL ptarg,pdelarg,foeew,zdelta
70 REAL tempe(ip1jmp1)
71
72 c fonction psat(T)
73
74 FOEEW ( PTARG,PDELARG ) = EXP (
75 * (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT)
76 * / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )
77
78 r2es = 380.11733
79 r3les = 17.269
80 r3ies = 21.875
81 r4les = 35.86
82 r4ies = 7.66
83 retv = 0.6077667
84 rtt = 273.16
85
86 c-- Calcul de Qsat en chaque point
87 c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
88 c pour eviter une exponentielle.
89
2/2
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
 19200 DO l = 1, llm
90
2/2
✓ Branch 0 taken 20386080 times.
✓ Branch 1 taken 18720 times.
 20404800 DO ij = 1, ip1jmp1
91 20404800 tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
92 ENDDO
93
2/2
✓ Branch 0 taken 20386080 times.
✓ Branch 1 taken 18720 times.
 20405280 DO ij = 1, ip1jmp1
94 20386080 zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) )
95 20386080 play = 0.5*(p(ij,l)+p(ij,l+1))
96 20386080 qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij),zdelta) / play )
97 20404800 qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )
98 ENDDO
99 ENDDO
100
101 c PRINT*,'Debut vlsplt version debug sans vlyqs'
102
103 480 zzpbar = 0.5 * pdt
104 zzw = pdt
105
2/2
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
 19200 DO l=1,llm
106
2/2
✓ Branch 0 taken 19150560 times.
✓ Branch 1 taken 18720 times.
 19169280 DO ij = iip2,ip1jm
107 19169280 mu(ij,l)=pbaru(ij,l) * zzpbar
108 ENDDO
109
2/2
✓ Branch 0 taken 19768320 times.
✓ Branch 1 taken 18720 times.
 19787040 DO ij=1,ip1jm
110 19787040 mv(ij,l)=pbarv(ij,l) * zzpbar
111 ENDDO
112
2/2
✓ Branch 0 taken 20386080 times.
✓ Branch 1 taken 18720 times.
 20405280 DO ij=1,ip1jmp1
113 20404800 mw(ij,l)=w(ij,l) * zzw
114 ENDDO
115 ENDDO
116
117
2/2
✓ Branch 0 taken 522720 times.
✓ Branch 1 taken 480 times.
 523200 DO ij=1,ip1jmp1
118 523200 mw(ij,llm+1)=0.
119 ENDDO
120
121 480 CALL SCOPY(ijp1llm,q(1,1,iq),1,zq(1,1,iq),1)
122 480 CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1)
123
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 480 times.
 480 if (nqdesc(iq).gt.0) then
124 ✗ do ifils=1,nqdesc(iq)
125 ✗ iq2=iqfils(ifils,iq)
126 ✗ CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1)
127 enddo
128 endif !if (nqfils(iq).gt.0) then
129
130 c call minmaxq(zq,qmin,qmax,'avant vlxqs ')
131 480 call vlxqs(zq,pente_max,zm,mu,qsat,iq)
132
133 c call minmaxq(zq,qmin,qmax,'avant vlyqs ')
134
135 480 call vlyqs(zq,pente_max,zm,mv,qsat,iq)
136
137 c call minmaxq(zq,qmin,qmax,'avant vlz ')
138
139 480 call vlz(zq,pente_max,zm,mw,iq)
140
141 c call minmaxq(zq,qmin,qmax,'avant vlyqs ')
142 c call minmaxq(zm,qmin,qmax,'M avant vlyqs ')
143
144 480 call vlyqs(zq,pente_max,zm,mv,qsat,iq)
145
146 c call minmaxq(zq,qmin,qmax,'avant vlxqs ')
147 c call minmaxq(zm,qmin,qmax,'M avant vlxqs ')
148
149 480 call vlxqs(zq,pente_max,zm,mu,qsat,iq)
150
151 c call minmaxq(zq,qmin,qmax,'apres vlxqs ')
152 c call minmaxq(zm,qmin,qmax,'M apres vlxqs ')
153
154
155
2/2
✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 480 times.
 19200 DO l=1,llm
156
2/2
✓ Branch 0 taken 20386080 times.
✓ Branch 1 taken 18720 times.
 20404800 DO ij=1,ip1jmp1
157 20404800 q(ij,l,iq)=zq(ij,l,iq)
158 ENDDO
159 480 DO ij=1,ip1jm+1,iip1
160
2/2
✓ Branch 0 taken 599040 times.
✓ Branch 1 taken 18720 times.
 617760 q(ij+iim,l,iq)=q(ij,l,iq)
161 ENDDO
162 ENDDO
163 ! CRisi: aussi pour les fils
164
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 480 times.
 480 if (nqdesc(iq).gt.0) then
165 ✗ do ifils=1,nqdesc(iq)
166 ✗ iq2=iqfils(ifils,iq)
167 ✗ DO l=1,llm
168 ✗ DO ij=1,ip1jmp1
169 ✗ q(ij,l,iq2)=zq(ij,l,iq2)
170 ENDDO
171 ✗ DO ij=1,ip1jm+1,iip1
172 ✗ q(ij+iim,l,iq2)=q(ij,l,iq2)
173 ENDDO
174 ENDDO
175 enddo !do ifils=1,nqdesc(iq)
176 endif ! if (nqfils(iq).gt.0) then
177 !write(*,*) 'vlspltqs 183: fin de la routine'
178
179 480 RETURN
180 END
181 960 SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat,iq)
182 USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
183
184 c
185 c Auteurs: P.Le Van, F.Hourdin, F.Forget
186 c
187 c ********************************************************************
188 c Shema d'advection " pseudo amont " .
189 c ********************************************************************
190 c
191 c --------------------------------------------------------------------
192 IMPLICIT NONE
193 c
194 include "dimensions.h"
195 include "paramet.h"
196 c
197 c
198 c Arguments:
199 c ----------
200 REAL masse(ip1jmp1,llm,nqtot),pente_max
201 REAL u_m( ip1jmp1,llm )
202 REAL q(ip1jmp1,llm,nqtot)
203 REAL qsat(ip1jmp1,llm)
204 INTEGER iq ! CRisi
205 c
206 c Local
207 c ---------
208 c
209 INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
210 INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
211 c
212 REAL new_m,zu_m,zdum(ip1jmp1,llm)
213 REAL dxq(ip1jmp1,llm),dxqu(ip1jmp1)
214 REAL zz(ip1jmp1)
215 REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
216 REAL u_mq(ip1jmp1,llm)
217
218 ! CRisi
219 1920 REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot)
220 INTEGER ifils,iq2 ! CRisi
221
222 Logical first,testcpu
223 SAVE first,testcpu
224
225 REAL SSUM
226 REAL temps0,temps1,temps2,temps3,temps4,temps5
227 SAVE temps0,temps1,temps2,temps3,temps4,temps5
228
229
230 DATA first,testcpu/.true.,.false./
231
232
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 959 times.
 960 IF(first) THEN
233 1 temps1=0.
234 1 temps2=0.
235 1 temps3=0.
236 1 temps4=0.
237 1 temps5=0.
238 1 first=.false.
239 ENDIF
240
241 c calcul de la pente a droite et a gauche de la maille
242
243
244
1/2
✓ Branch 0 taken 960 times.
✗ Branch 1 not taken.
 960 IF (pente_max.gt.-1.e-5) THEN
245 c IF (pente_max.gt.10) THEN
246
247 c calcul des pentes avec limitation, Van Leer scheme I:
248 c -----------------------------------------------------
249
250 c calcul de la pente aux points u
251
2/2
✓ Branch 0 taken 960 times.
✓ Branch 1 taken 37440 times.
 38400 DO l = 1, llm
252
2/2
✓ Branch 0 taken 38263680 times.
✓ Branch 1 taken 37440 times.
 38301120 DO ij=iip2,ip1jm-1
253 38301120 dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
254 c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
255 c sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
256 ENDDO
257 1160640 DO ij=iip1+iip1,ip1jm,iip1
258
2/2
✓ Branch 0 taken 1123200 times.
✓ Branch 1 taken 37440 times.
 1160640 dxqu(ij)=dxqu(ij-iim)
259 c sigu(ij)=sigu(ij-iim)
260 ENDDO
261
262
2/2
✓ Branch 0 taken 38301120 times.
✓ Branch 1 taken 37440 times.
 38338560 DO ij=iip2,ip1jm
263 38338560 adxqu(ij)=abs(dxqu(ij))
264 ENDDO
265
266 c calcul de la pente maximum dans la maille en valeur absolue
267
268
2/2
✓ Branch 0 taken 38263680 times.
✓ Branch 1 taken 37440 times.
 38301120 DO ij=iip2+1,ip1jm
269 dxqmax(ij,l)=pente_max*
270 38301120 , min(adxqu(ij-1),adxqu(ij))
271 c limitation subtile
272 c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
273
274
275 ENDDO
276
277 1160640 DO ij=iip1+iip1,ip1jm,iip1
278
2/2
✓ Branch 0 taken 1123200 times.
✓ Branch 1 taken 37440 times.
 1160640 dxqmax(ij-iim,l)=dxqmax(ij,l)
279 ENDDO
280
281
2/2
✓ Branch 0 taken 38263680 times.
✓ Branch 1 taken 37440 times.
 38302080 DO ij=iip2+1,ip1jm
282
2/2
✓ Branch 0 taken 27096747 times.
✓ Branch 1 taken 11166933 times.
 38263680 IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
283 27096747 dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
284 ELSE
285 c extremum local
286 11166933 dxq(ij,l)=0.
287 ENDIF
288 38263680 dxq(ij,l)=0.5*dxq(ij,l)
289 dxq(ij,l)=
290 38301120 , sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
291 ENDDO
292
293 ENDDO ! l=1,llm
294
295 ELSE ! (pente_max.lt.-1.e-5)
296
297 c Pentes produits:
298 c ----------------
299
300 ✗ DO l = 1, llm
301 ✗ DO ij=iip2,ip1jm-1
302 ✗ dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
303 ENDDO
304 ✗ DO ij=iip1+iip1,ip1jm,iip1
305 ✗ dxqu(ij)=dxqu(ij-iim)
306 ENDDO
307
308 ✗ DO ij=iip2+1,ip1jm
309 ✗ zz(ij)=dxqu(ij-1)*dxqu(ij)
310 ✗ zz(ij)=zz(ij)+zz(ij)
311 ✗ IF(zz(ij).gt.0) THEN
312 ✗ dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
313 ELSE
314 c extremum local
315 ✗ dxq(ij,l)=0.
316 ENDIF
317 ENDDO
318
319 ENDDO
320
321 ENDIF ! (pente_max.lt.-1.e-5)
322
323 c bouclage de la pente en iip1:
324 c -----------------------------
325
326
2/2
✓ Branch 0 taken 37440 times.
✓ Branch 1 taken 960 times.
 38400 DO l=1,llm
327 1160640 DO ij=iip1+iip1,ip1jm,iip1
328
2/2
✓ Branch 0 taken 1123200 times.
✓ Branch 1 taken 37440 times.
 1160640 dxq(ij-iim,l)=dxq(ij,l)
329 ENDDO
330
331
2/2
✓ Branch 0 taken 40772160 times.
✓ Branch 1 taken 37440 times.
 40810560 DO ij=1,ip1jmp1
332 40809600 iadvplus(ij,l)=0
333 ENDDO
334
335 ENDDO
336
337
338 c calcul des flux a gauche et a droite
339
340 c on cumule le flux correspondant a toutes les mailles dont la masse
341 c au travers de la paroi pENDant le pas de temps.
342 c le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
343
2/2
✓ Branch 0 taken 37440 times.
✓ Branch 1 taken 960 times.
 38400 DO l=1,llm
344
2/2
✓ Branch 0 taken 38263680 times.
✓ Branch 1 taken 37440 times.
 38302080 DO ij=iip2,ip1jm-1
345
2/2
✓ Branch 0 taken 27282878 times.
✓ Branch 1 taken 10980802 times.
 38301120 IF (u_m(ij,l).gt.0.) THEN
346 27282878 zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
347 u_mq(ij,l)=u_m(ij,l)*
348 27282878 $ min(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
349 ELSE
350 10980802 zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
351 u_mq(ij,l)=u_m(ij,l)*
352 10980802 $ min(q(ij+1,l,iq)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
353 ENDIF
354 ENDDO
355 ENDDO
356
357
358 c detection des points ou on advecte plus que la masse de la
359 c maille
360
2/2
✓ Branch 0 taken 37440 times.
✓ Branch 1 taken 960 times.
 38400 DO l=1,llm
361
2/2
✓ Branch 0 taken 38263680 times.
✓ Branch 1 taken 37440 times.
 38302080 DO ij=iip2,ip1jm-1
362
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 38263680 times.
 38301120 IF(zdum(ij,l).lt.0) THEN
363 ✗ iadvplus(ij,l)=1
364 ✗ u_mq(ij,l)=0.
365 ENDIF
366 ENDDO
367 ENDDO
368
2/2
✓ Branch 0 taken 37440 times.
✓ Branch 1 taken 960 times.
 38400 DO l=1,llm
369 960 DO ij=iip1+iip1,ip1jm,iip1
370
2/2
✓ Branch 0 taken 1123200 times.
✓ Branch 1 taken 37440 times.
 1160640 iadvplus(ij,l)=iadvplus(ij-iim,l)
371 ENDDO
372 ENDDO
373
374
375
376 c traitement special pour le cas ou on advecte en longitude plus que le
377 c contenu de la maille.
378 c cette partie est mal vectorisee.
379
380 c pas d'influence de la pression saturante (pour l'instant)
381
382 c calcul du nombre de maille sur lequel on advecte plus que la maille.
383
384 n0=0
385
2/2
✓ Branch 0 taken 37440 times.
✓ Branch 1 taken 960 times.
 38400 DO l=1,llm
386 37440 nl(l)=0
387
2/2
✓ Branch 0 taken 38301120 times.
✓ Branch 1 taken 37440 times.
 38338560 DO ij=iip2,ip1jm
388 38338560 nl(l)=nl(l)+iadvplus(ij,l)
389 ENDDO
390 38400 n0=n0+nl(l)
391 ENDDO
392
393
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 960 times.
 960 IF(n0.gt.0) THEN
394 ccc PRINT*,'Nombre de points pour lesquels on advect plus que le'
395 ccc & ,'contenu de la maille : ',n0
396
397 ✗ DO l=1,llm
398 ✗ IF(nl(l).gt.0) THEN
399 iju=0
400 c indicage des mailles concernees par le traitement special
401 ✗ DO ij=iip2,ip1jm
402 ✗ IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
403 ✗ iju=iju+1
404 ✗ indu(iju)=ij
405 ENDIF
406 ENDDO
407 niju=iju
408 c PRINT*,'niju,nl',niju,nl(l)
409
410 c traitement des mailles
411 ✗ DO iju=1,niju
412 ✗ ij=indu(iju)
413 ✗ j=(ij-1)/iip1+1
414 ✗ zu_m=u_m(ij,l)
415 ✗ u_mq(ij,l)=0.
416 ✗ IF(zu_m.gt.0.) THEN
417 ijq=ij
418 ✗ i=ijq-(j-1)*iip1
419 c accumulation pour les mailles completements advectees
420 ✗ do while(zu_m.gt.masse(ijq,l,iq))
421 u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq)
422 ✗ & *masse(ijq,l,iq)
423 ✗ zu_m=zu_m-masse(ijq,l,iq)
424 ✗ i=mod(i-2+iim,iim)+1
425 ✗ ijq=(j-1)*iip1+i
426 ENDDO
427 c ajout de la maille non completement advectee
428 u_mq(ij,l)=u_mq(ij,l)+zu_m*
429 & (q(ijq,l,iq)+0.5*(1.-zu_m/masse(ijq,l,iq))
430 ✗ & *dxq(ijq,l))
431 ELSE
432 ✗ ijq=ij+1
433 ✗ i=ijq-(j-1)*iip1
434 c accumulation pour les mailles completements advectees
435 ✗ do while(-zu_m.gt.masse(ijq,l,iq))
436 u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
437 ✗ & *masse(ijq,l,iq)
438 ✗ zu_m=zu_m+masse(ijq,l,iq)
439 ✗ i=mod(i,iim)+1
440 ✗ ijq=(j-1)*iip1+i
441 ENDDO
442 c ajout de la maille non completement advectee
443 u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
444 ✗ & 0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l))
445 ENDIF
446 ENDDO
447 ENDIF
448 ENDDO
449 ENDIF ! n0.gt.0
450
451
452
453 c bouclage en latitude
454
455
2/2
✓ Branch 0 taken 37440 times.
✓ Branch 1 taken 960 times.
 38400 DO l=1,llm
456 960 DO ij=iip1+iip1,ip1jm,iip1
457
2/2
✓ Branch 0 taken 1123200 times.
✓ Branch 1 taken 37440 times.
 1160640 u_mq(ij,l)=u_mq(ij-iim,l)
458 ENDDO
459 ENDDO
460
461 ! CRisi: appel récursif de l'advection sur les fils.
462 ! Il faut faire ça avant d'avoir mis à jour q et masse
463 !write(*,*) 'vlspltqs 326: iq,nqfils(iq)=',iq,nqfils(iq)
464
465
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 960 times.
 960 if (nqfils(iq).gt.0) then
466 ✗ do ifils=1,nqdesc(iq)
467 ✗ iq2=iqfils(ifils,iq)
468 ✗ DO l=1,llm
469 ✗ DO ij=iip2,ip1jm
470 ! On a besoin de q et masse seulement entre iip2 et ip1jm
471 ✗ masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
472 ✗ Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
473 enddo
474 enddo
475 enddo !do ifils=1,nqdesc(iq)
476 ✗ do ifils=1,nqfils(iq)
477 ✗ iq2=iqfils(ifils,iq)
478 ✗ call vlx(Ratio,pente_max,masseq,u_mq,iq2)
479 enddo !do ifils=1,nqfils(iq)
480 endif !if (nqfils(iq).gt.0) then
481 ! end CRisi
482
483 c calcul des tendances
484
485
2/2
✓ Branch 0 taken 37440 times.
✓ Branch 1 taken 960 times.
 38400 DO l=1,llm
486
2/2
✓ Branch 0 taken 38263680 times.
✓ Branch 1 taken 37440 times.
 38301120 DO ij=iip2+1,ip1jm
487 38263680 new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l)
488 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
489 & u_mq(ij-1,l)-u_mq(ij,l))
490 38263680 & /new_m
491 38301120 masse(ij,l,iq)=new_m
492 ENDDO
493 c Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
494 960 DO ij=iip1+iip1,ip1jm,iip1
495 1160640 q(ij-iim,l,iq)=q(ij,l,iq)
496
2/2
✓ Branch 0 taken 1123200 times.
✓ Branch 1 taken 37440 times.
 1160640 masse(ij-iim,l,iq)=masse(ij,l,iq)
497 ENDDO
498 ENDDO
499
500 ! retablir les fils en rapport de melange par rapport a l'air:
501 ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio
502 ! puis on boucle en longitude
503
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 960 times.
 960 if (nqdesc(iq).gt.0) then
504 ✗ do ifils=1,nqdesc(iq)
505 ✗ iq2=iqfils(ifils,iq)
506 ✗ DO l=1,llm
507 ✗ DO ij=iip2+1,ip1jm
508 ✗ q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)
509 enddo
510 ✗ DO ij=iip1+iip1,ip1jm,iip1
511 ✗ q(ij-iim,l,iq2)=q(ij,l,iq2)
512 enddo ! DO ij=ijb+iip1-1,ije,iip1
513 enddo !DO l=1,llm
514 enddo !do ifils=1,nqdesc(iq)
515 endif !if (nqfils(iq).gt.0) then
516
517 c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
518 c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
519
520
521 960 RETURN
522 END
523 150720 SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat,iq)
524 USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
525 c
526 c Auteurs: P.Le Van, F.Hourdin, F.Forget
527 c
528 c ********************************************************************
529 c Shema d'advection " pseudo amont " .
530 c ********************************************************************
531 c q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
532 c qsat est un argument de sortie pour le s-pg ....
533 c
534 c
535 c --------------------------------------------------------------------
536
537 USE comconst_mod, ONLY: pi
538
539 IMPLICIT NONE
540 c
541 include "dimensions.h"
542 include "paramet.h"
543 include "comgeom.h"
544 c
545 c
546 c Arguments:
547 c ----------
548 REAL masse(ip1jmp1,llm,nqtot),pente_max
549 REAL masse_adv_v( ip1jm,llm)
550 REAL q(ip1jmp1,llm,nqtot)
551 REAL qsat(ip1jmp1,llm)
552 INTEGER iq ! CRisi
553 c
554 c Local
555 c ---------
556 c
557 INTEGER i,ij,l
558 c
559 REAL airej2,airejjm,airescb(iim),airesch(iim)
560 REAL dyq(ip1jmp1,llm),dyqv(ip1jm)
561 REAL adyqv(ip1jm),dyqmax(ip1jmp1)
562 REAL qbyv(ip1jm,llm)
563
564 REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
565 c REAL newq,oldmasse
566 Logical first,testcpu
567 REAL temps0,temps1,temps2,temps3,temps4,temps5
568 SAVE temps0,temps1,temps2,temps3,temps4,temps5
569 SAVE first,testcpu
570
571 REAL convpn,convps,convmpn,convmps
572 REAL sinlon(iip1),sinlondlon(iip1)
573 REAL coslon(iip1),coslondlon(iip1)
574 SAVE sinlon,coslon,sinlondlon,coslondlon
575 SAVE airej2,airejjm
576
577 1920 REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
578 INTEGER ifils,iq2 ! CRisi
579 c
580 c
581 REAL SSUM
582
583 DATA first,testcpu/.true.,.false./
584 DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
585
586
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 959 times.
 960 IF(first) THEN
587 1 PRINT*,'Shema Amont nouveau appele dans Vanleer '
588 1 first=.false.
589
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 1 times.
 33 do i=2,iip1
590 32 coslon(i)=cos(rlonv(i))
591 32 sinlon(i)=sin(rlonv(i))
592 32 coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
593 33 sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
594 ENDDO
595 1 coslon(1)=coslon(iip1)
596 1 coslondlon(1)=coslondlon(iip1)
597 1 sinlon(1)=sinlon(iip1)
598 1 sinlondlon(1)=sinlondlon(iip1)
599 1 airej2 = SSUM( iim, aire(iip2), 1 )
600 1 airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
601 ENDIF
602
603 c
604
605
606
2/2
✓ Branch 0 taken 37440 times.
✓ Branch 1 taken 960 times.
 38400 DO l = 1, llm
607 c
608 c --------------------------------
609 c CALCUL EN LATITUDE
610 c --------------------------------
611
612 c On commence par calculer la valeur du traceur moyenne sur le premier cercle
613 c de latitude autour du pole (qpns pour le pole nord et qpsn pour
614 c le pole nord) qui sera utilisee pour evaluer les pentes au pole.
615
616
2/2
✓ Branch 0 taken 1198080 times.
✓ Branch 1 taken 37440 times.
 1235520 DO i = 1, iim
617 1198080 airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
618 1235520 airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
619 ENDDO
620 37440 qpns = SSUM( iim, airescb ,1 ) / airej2
621 37440 qpsn = SSUM( iim, airesch ,1 ) / airejjm
622
623 c calcul des pentes aux points v
624
625
2/2
✓ Branch 0 taken 39536640 times.
✓ Branch 1 taken 37440 times.
 39574080 DO ij=1,ip1jm
626 39536640 dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
627 39574080 adyqv(ij)=abs(dyqv(ij))
628 ENDDO
629
630 c calcul des pentes aux points scalaires
631
632
2/2
✓ Branch 0 taken 38301120 times.
✓ Branch 1 taken 37440 times.
 38338560 DO ij=iip2,ip1jm
633 38301120 dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
634 38301120 dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
635 38338560 dyqmax(ij)=pente_max*dyqmax(ij)
636 ENDDO
637
638 c calcul des pentes aux poles
639
640
2/2
✓ Branch 0 taken 1235520 times.
✓ Branch 1 taken 37440 times.
 1272960 DO ij=1,iip1
641 1235520 dyq(ij,l)=qpns-q(ij+iip1,l,iq)
642 1272960 dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
643 ENDDO
644
645 c filtrage de la derivee
646 dyn1=0.
647 dys1=0.
648 dyn2=0.
649 dys2=0.
650
2/2
✓ Branch 0 taken 37440 times.
✓ Branch 1 taken 1198080 times.
 1235520 DO ij=1,iim
651 1198080 dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
652 1198080 dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
653 1198080 dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
654 1235520 dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
655 ENDDO
656
2/2
✓ Branch 0 taken 1235520 times.
✓ Branch 1 taken 37440 times.
 1272960 DO ij=1,iip1
657 1235520 dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
658 1272960 dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
659 ENDDO
660
661 c calcul des pentes limites aux poles
662
663 fn=1.
664 fs=1.
665
2/2
✓ Branch 0 taken 37440 times.
✓ Branch 1 taken 1198080 times.
 1235520 DO ij=1,iim
666
2/2
✓ Branch 0 taken 202056 times.
✓ Branch 1 taken 996024 times.
 1198080 IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
667 202056 fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
668 ENDIF
669
2/2
✓ Branch 0 taken 245116 times.
✓ Branch 1 taken 952964 times.
 1235520 IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
670 245116 fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
671 ENDIF
672 ENDDO
673
2/2
✓ Branch 0 taken 1235520 times.
✓ Branch 1 taken 37440 times.
 1272960 DO ij=1,iip1
674 1235520 dyq(ij,l)=fn*dyq(ij,l)
675 1272960 dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
676 ENDDO
677
678 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
679 C En memoire de dIFferents tests sur la
680 C limitation des pentes aux poles.
681 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
682 C PRINT*,dyq(1)
683 C PRINT*,dyqv(iip1+1)
684 C appn=abs(dyq(1)/dyqv(iip1+1))
685 C PRINT*,dyq(ip1jm+1)
686 C PRINT*,dyqv(ip1jm-iip1+1)
687 C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
688 C DO ij=2,iim
689 C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
690 C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
691 C ENDDO
692 C appn=min(pente_max/appn,1.)
693 C apps=min(pente_max/apps,1.)
694 C
695 C
696 C cas ou on a un extremum au pole
697 C
698 C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
699 C & appn=0.
700 C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
701 C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
702 C & apps=0.
703 C
704 C limitation des pentes aux poles
705 C DO ij=1,iip1
706 C dyq(ij)=appn*dyq(ij)
707 C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
708 C ENDDO
709 C
710 C test
711 C DO ij=1,iip1
712 C dyq(iip1+ij)=0.
713 C dyq(ip1jm+ij-iip1)=0.
714 C ENDDO
715 C DO ij=1,ip1jmp1
716 C dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
717 C ENDDO
718 C
719 C changement 10 07 96
720 C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
721 C & THEN
722 C DO ij=1,iip1
723 C dyqmax(ij)=0.
724 C ENDDO
725 C ELSE
726 C DO ij=1,iip1
727 C dyqmax(ij)=pente_max*abs(dyqv(ij))
728 C ENDDO
729 C ENDIF
730 C
731 C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
732 C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
733 C &THEN
734 C DO ij=ip1jm+1,ip1jmp1
735 C dyqmax(ij)=0.
736 C ENDDO
737 C ELSE
738 C DO ij=ip1jm+1,ip1jmp1
739 C dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
740 C ENDDO
741 C ENDIF
742 C fin changement 10 07 96
743 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
744
745 c calcul des pentes limitees
746
747
2/2
✓ Branch 0 taken 38301120 times.
✓ Branch 1 taken 37440 times.
 38339520 DO ij=iip2,ip1jm
748
2/2
✓ Branch 0 taken 26560912 times.
✓ Branch 1 taken 11740208 times.
 38338560 IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
749 26560912 dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
750 ELSE
751 11740208 dyq(ij,l)=0.
752 ENDIF
753 ENDDO
754
755 ENDDO
756
757
2/2
✓ Branch 0 taken 37440 times.
✓ Branch 1 taken 960 times.
 38400 DO l=1,llm
758
2/2
✓ Branch 0 taken 39536640 times.
✓ Branch 1 taken 37440 times.
 39575040 DO ij=1,ip1jm
759
2/2
✓ Branch 0 taken 20058114 times.
✓ Branch 1 taken 19478526 times.
 39536640 IF( masse_adv_v(ij,l).GT.0. ) THEN
760 qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l,iq ) +
761 , dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)
762 20058114 , /masse(ij+iip1,l,iq)))
763 ELSE
764 qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l,iq) - dyq(ij,l) *
765 19478526 , 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq)) )
766 ENDIF
767 39574080 qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
768 ENDDO
769 ENDDO
770
771
772 ! CRisi: appel récursif de l'advection sur les fils.
773 ! Il faut faire ça avant d'avoir mis à jour q et masse
774 !write(*,*) 'vlyqs 689: iq,nqfils(iq)=',iq,nqfils(iq)
775
776
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 960 times.
 960 if (nqfils(iq).gt.0) then
777 ✗ do ifils=1,nqdesc(iq)
778 ✗ iq2=iqfils(ifils,iq)
779 ✗ DO l=1,llm
780 ✗ DO ij=1,ip1jmp1
781 ✗ masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
782 ✗ Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
783 enddo
784 enddo
785 enddo !do ifils=1,nqdesc(iq)
786
787 ✗ do ifils=1,nqfils(iq)
788 ✗ iq2=iqfils(ifils,iq)
789 !write(*,*) 'vlyqs 783: appel rec de vly, iq2=',iq2
790 ✗ call vly(Ratio,pente_max,masseq,qbyv,iq2)
791 enddo !do ifils=1,nqfils(iq)
792 endif !if (nqfils(iq).gt.0) then
793
794
2/2
✓ Branch 0 taken 37440 times.
✓ Branch 1 taken 960 times.
 38400 DO l=1,llm
795
2/2
✓ Branch 0 taken 38301120 times.
✓ Branch 1 taken 37440 times.
 38338560 DO ij=iip2,ip1jm
796 newmasse=masse(ij,l,iq)
797 38301120 & +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
798 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l)
799 38301120 & -qbyv(ij-iip1,l))/newmasse
800 38338560 masse(ij,l,iq)=newmasse
801 ENDDO
802 c.-. ancienne version
803 37440 convpn=SSUM(iim,qbyv(1,l),1)/apoln
804 37440 convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
805
2/2
✓ Branch 0 taken 1235520 times.
✓ Branch 1 taken 37440 times.
 1272960 DO ij = 1,iip1
806 1235520 newmasse=masse(ij,l,iq)+convmpn*aire(ij)
807 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convpn*aire(ij))/
808 1235520 & newmasse
809 1272960 masse(ij,l,iq)=newmasse
810 ENDDO
811 37440 convps = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
812 37440 convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
813
2/2
✓ Branch 0 taken 1235520 times.
✓ Branch 1 taken 37440 times.
 1273920 DO ij = ip1jm+1,ip1jmp1
814 1235520 newmasse=masse(ij,l,iq)+convmps*aire(ij)
815 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convps*aire(ij))/
816 1235520 & newmasse
817 1272960 masse(ij,l,iq)=newmasse
818 ENDDO
819 c.-. fin ancienne version
820
821 c._. nouvelle version
822 c convpn=SSUM(iim,qbyv(1,l),1)
823 c convmpn=ssum(iim,masse_adv_v(1,l),1)
824 c oldmasse=ssum(iim,masse(1,l),1)
825 c newmasse=oldmasse+convmpn
826 c newq=(q(1,l)*oldmasse+convpn)/newmasse
827 c newmasse=newmasse/apoln
828 c DO ij = 1,iip1
829 c q(ij,l)=newq
830 c masse(ij,l,iq)=newmasse*aire(ij)
831 c ENDDO
832 c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
833 c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
834 c oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
835 c newmasse=oldmasse+convmps
836 c newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
837 c newmasse=newmasse/apols
838 c DO ij = ip1jm+1,ip1jmp1
839 c q(ij,l)=newq
840 c masse(ij,l,iq)=newmasse*aire(ij)
841 c ENDDO
842 c._. fin nouvelle version
843 ENDDO
844
845 !write(*,*) 'vly 866'
846
847 ! retablir les fils en rapport de melange par rapport a l'air:
848
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 960 times.
 960 if (nqdesc(iq).gt.0) then
849 ✗ do ifils=1,nqdesc(iq)
850 ✗ iq2=iqfils(ifils,iq)
851 ✗ DO l=1,llm
852 ✗ DO ij=1,ip1jmp1
853 ✗ q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)
854 enddo
855 enddo
856 enddo !do ifils=1,nqdesc(iq)
857 endif !if (nqfils(iq).gt.0) then
858 !write(*,*) 'vly 879'
859
860 960 RETURN
861 END
862