GCC Code Coverage Report

Directory: ./
File: dyn/vlsplt.f
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 205 304 67.4%
Branches: 145 240 60.4%
Line Branch Exec Source
1 c
2 c $Id: vlsplt.F 2603 2016-07-25 09:31:56Z emillour $
3 c
4
5 1920 SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt,iq)
6 USE infotrac, ONLY: nqtot,nqdesc,iqfils
7 c
8 c Auteurs: P.Le Van, F.Hourdin, F.Forget
9 c
10 c ********************************************************************
11 c Shema d'advection " pseudo amont " .
12 c ********************************************************************
13 c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
14 c
15 c pente_max facteur de limitation des pentes: 2 en general
16 c 0 pour un schema amont
17 c pbaru,pbarv,w flux de masse en u ,v ,w
18 c pdt pas de temps
19 c
20 c --------------------------------------------------------------------
21 IMPLICIT NONE
22 c
23 include "dimensions.h"
24 include "paramet.h"
25
26 c
27 c Arguments:
28 c ----------
29 REAL masse(ip1jmp1,llm),pente_max
30 c REAL masse(iip1,jjp1,llm),pente_max
31 REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
32 REAL q(ip1jmp1,llm,nqtot)
33 c REAL q(iip1,jjp1,llm)
34 REAL w(ip1jmp1,llm),pdt
35 INTEGER iq ! CRisi
36 c
37 c Local
38 c ---------
39 c
40 INTEGER i,ij,l,j,ii
41 INTEGER ijlqmin,iqmin,jqmin,lqmin
42 c
43 3840 REAL zm(ip1jmp1,llm,nqtot),newmasse
44 REAL mu(ip1jmp1,llm)
45 REAL mv(ip1jm,llm)
46 REAL mw(ip1jmp1,llm+1)
47 3840 REAL zq(ip1jmp1,llm,nqtot),zz
48 REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm)
49 REAL second,temps0,temps1,temps2,temps3
50 REAL ztemps1,ztemps2,ztemps3
51 REAL zzpbar, zzw
52 LOGICAL testcpu
53 SAVE testcpu
54 SAVE temps1,temps2,temps3
55 INTEGER iminn,imaxx
56 INTEGER ifils,iq2 ! CRisi
57
58 REAL qmin,qmax
59 DATA qmin,qmax/0.,1.e33/
60 DATA testcpu/.false./
61 DATA temps1,temps2,temps3/0.,0.,0./
62
63
64 1920 zzpbar = 0.5 * pdt
65 zzw = pdt
66
2/2
✓ Branch 0 taken 74880 times.
✓ Branch 1 taken 1920 times.
 76800 DO l=1,llm
67
2/2
✓ Branch 0 taken 74880 times.
✓ Branch 1 taken 76602240 times.
 76677120 DO ij = iip2,ip1jm
68 76677120 mu(ij,l)=pbaru(ij,l) * zzpbar
69 ENDDO
70
2/2
✓ Branch 0 taken 74880 times.
✓ Branch 1 taken 79073280 times.
 79148160 DO ij=1,ip1jm
71 79148160 mv(ij,l)=pbarv(ij,l) * zzpbar
72 ENDDO
73
2/2
✓ Branch 0 taken 81544320 times.
✓ Branch 1 taken 74880 times.
 81621120 DO ij=1,ip1jmp1
74 81619200 mw(ij,l)=w(ij,l) * zzw
75 ENDDO
76 ENDDO
77
78
2/2
✓ Branch 0 taken 2090880 times.
✓ Branch 1 taken 1920 times.
 2092800 DO ij=1,ip1jmp1
79 2092800 mw(ij,llm+1)=0.
80 ENDDO
81
82 1920 CALL SCOPY(ijp1llm,q(1,1,iq),1,zq(1,1,iq),1)
83 1920 CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1)
84
85
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
 1920 if (nqdesc(iq).gt.0) then
86 ✗ do ifils=1,nqdesc(iq)
87 ✗ iq2=iqfils(ifils,iq)
88 ✗ CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1)
89 enddo
90 endif !if (nqfils(iq).gt.0) then
91
92 cprint*,'Entree vlx1'
93 c call minmaxq(zq,qmin,qmax,'avant vlx ')
94 1920 call vlx(zq,pente_max,zm,mu,iq)
95 cprint*,'Sortie vlx1'
96 c call minmaxq(zq,qmin,qmax,'apres vlx1 ')
97
98 c print*,'Entree vly1'
99
100 1920 call vly(zq,pente_max,zm,mv,iq)
101 c call minmaxq(zq,qmin,qmax,'apres vly1 ')
102 cprint*,'Sortie vly1'
103 1920 call vlz(zq,pente_max,zm,mw,iq)
104 c call minmaxq(zq,qmin,qmax,'apres vlz ')
105
106
107 1920 call vly(zq,pente_max,zm,mv,iq)
108 c call minmaxq(zq,qmin,qmax,'apres vly ')
109
110
111 1920 call vlx(zq,pente_max,zm,mu,iq)
112 c call minmaxq(zq,qmin,qmax,'apres vlx2 ')
113
114
115
2/2
✓ Branch 0 taken 74880 times.
✓ Branch 1 taken 1920 times.
 76800 DO l=1,llm
116
2/2
✓ Branch 0 taken 81544320 times.
✓ Branch 1 taken 74880 times.
 81619200 DO ij=1,ip1jmp1
117 81619200 q(ij,l,iq)=zq(ij,l,iq)
118 ENDDO
119 1920 DO ij=1,ip1jm+1,iip1
120
2/2
✓ Branch 0 taken 2396160 times.
✓ Branch 1 taken 74880 times.
 2471040 q(ij+iim,l,iq)=q(ij,l,iq)
121 ENDDO
122 ENDDO
123 ! CRisi: aussi pour les fils
124
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
 1920 if (nqdesc(iq).gt.0) then
125 ✗ do ifils=1,nqdesc(iq)
126 ✗ iq2=iqfils(ifils,iq)
127 ✗ DO l=1,llm
128 ✗ DO ij=1,ip1jmp1
129 ✗ q(ij,l,iq2)=zq(ij,l,iq2)
130 ENDDO
131 ✗ DO ij=1,ip1jm+1,iip1
132 ✗ q(ij+iim,l,iq2)=q(ij,l,iq2)
133 ENDDO
134 ENDDO
135 enddo !do ifils=1,nqdesc(iq)
136 endif ! if (nqdesc(iq).gt.0) then
137
138 1920 RETURN
139 END
140 3840 RECURSIVE SUBROUTINE vlx(q,pente_max,masse,u_m,iq)
141 USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
142
143 c Auteurs: P.Le Van, F.Hourdin, F.Forget
144 c
145 c ********************************************************************
146 c Shema d'advection " pseudo amont " .
147 c ********************************************************************
148 c nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
149 c
150 c
151 c --------------------------------------------------------------------
152 IMPLICIT NONE
153 c
154 include "dimensions.h"
155 include "paramet.h"
156 include "iniprint.h"
157 c
158 c
159 c Arguments:
160 c ----------
161 REAL masse(ip1jmp1,llm,nqtot),pente_max
162 REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)
163 REAL q(ip1jmp1,llm,nqtot)
164 REAL w(ip1jmp1,llm)
165 INTEGER iq ! CRisi
166 c
167 c Local
168 c ---------
169 c
170 INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
171 INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
172 c
173 REAL new_m,zu_m,zdum(ip1jmp1,llm)
174 REAL sigu(ip1jmp1),dxq(ip1jmp1,llm),dxqu(ip1jmp1)
175 REAL zz(ip1jmp1)
176 REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
177 REAL u_mq(ip1jmp1,llm)
178
179 ! CRisi
180 7680 REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot)
181 INTEGER ifils,iq2 ! CRisi
182
183 Logical extremum,first,testcpu
184 SAVE first,testcpu
185
186 REAL SSUM
187 REAL temps0,temps1,temps2,temps3,temps4,temps5,second
188 SAVE temps0,temps1,temps2,temps3,temps4,temps5
189
190 REAL z1,z2,z3
191
192 DATA first,testcpu/.true.,.false./
193
194
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3839 times.
 3840 IF(first) THEN
195 1 temps1=0.
196 1 temps2=0.
197 1 temps3=0.
198 1 temps4=0.
199 1 temps5=0.
200 1 first=.false.
201 ENDIF
202
203 c calcul de la pente a droite et a gauche de la maille
204
205
206
1/2
✓ Branch 0 taken 3840 times.
✗ Branch 1 not taken.
 3840 IF (pente_max.gt.-1.e-5) THEN
207 c IF (pente_max.gt.10) THEN
208
209 c calcul des pentes avec limitation, Van Leer scheme I:
210 c -----------------------------------------------------
211
212 c calcul de la pente aux points u
213
2/2
✓ Branch 0 taken 3840 times.
✓ Branch 1 taken 149760 times.
 153600 DO l = 1, llm
214
2/2
✓ Branch 0 taken 153054720 times.
✓ Branch 1 taken 149760 times.
 153204480 DO ij=iip2,ip1jm-1
215 153204480 dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
216 c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
217 c sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
218 ENDDO
219 4642560 DO ij=iip1+iip1,ip1jm,iip1
220
2/2
✓ Branch 0 taken 4492800 times.
✓ Branch 1 taken 149760 times.
 4642560 dxqu(ij)=dxqu(ij-iim)
221 c sigu(ij)=sigu(ij-iim)
222 ENDDO
223
224
2/2
✓ Branch 0 taken 153204480 times.
✓ Branch 1 taken 149760 times.
 153354240 DO ij=iip2,ip1jm
225 153354240 adxqu(ij)=abs(dxqu(ij))
226 ENDDO
227
228 c calcul de la pente maximum dans la maille en valeur absolue
229
230
2/2
✓ Branch 0 taken 153054720 times.
✓ Branch 1 taken 149760 times.
 153204480 DO ij=iip2+1,ip1jm
231 dxqmax(ij,l)=pente_max*
232 153204480 , min(adxqu(ij-1),adxqu(ij))
233 c limitation subtile
234 c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
235
236
237 ENDDO
238
239 4642560 DO ij=iip1+iip1,ip1jm,iip1
240
2/2
✓ Branch 0 taken 4492800 times.
✓ Branch 1 taken 149760 times.
 4642560 dxqmax(ij-iim,l)=dxqmax(ij,l)
241 ENDDO
242
243
2/2
✓ Branch 0 taken 153054720 times.
✓ Branch 1 taken 149760 times.
 153208320 DO ij=iip2+1,ip1jm
244
2/2
✓ Branch 0 taken 65216412 times.
✓ Branch 1 taken 87838308 times.
 153054720 IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
245 65216412 dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
246 ELSE
247 c extremum local
248 87838308 dxq(ij,l)=0.
249 ENDIF
250 153054720 dxq(ij,l)=0.5*dxq(ij,l)
251 dxq(ij,l)=
252 153204480 , sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
253 ENDDO
254
255 ENDDO ! l=1,llm
256 cprint*,'Ok calcul des pentes'
257
258 ELSE ! (pente_max.lt.-1.e-5)
259
260 c Pentes produits:
261 c ----------------
262
263 ✗ DO l = 1, llm
264 ✗ DO ij=iip2,ip1jm-1
265 ✗ dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
266 ENDDO
267 ✗ DO ij=iip1+iip1,ip1jm,iip1
268 ✗ dxqu(ij)=dxqu(ij-iim)
269 ENDDO
270
271 ✗ DO ij=iip2+1,ip1jm
272 ✗ zz(ij)=dxqu(ij-1)*dxqu(ij)
273 ✗ zz(ij)=zz(ij)+zz(ij)
274 ✗ IF(zz(ij).gt.0) THEN
275 ✗ dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
276 ELSE
277 c extremum local
278 ✗ dxq(ij,l)=0.
279 ENDIF
280 ENDDO
281
282 ENDDO
283
284 ENDIF ! (pente_max.lt.-1.e-5)
285
286 c bouclage de la pente en iip1:
287 c -----------------------------
288
289
2/2
✓ Branch 0 taken 149760 times.
✓ Branch 1 taken 3840 times.
 153600 DO l=1,llm
290 4642560 DO ij=iip1+iip1,ip1jm,iip1
291
2/2
✓ Branch 0 taken 4492800 times.
✓ Branch 1 taken 149760 times.
 4642560 dxq(ij-iim,l)=dxq(ij,l)
292 ENDDO
293
2/2
✓ Branch 0 taken 163088640 times.
✓ Branch 1 taken 149760 times.
 163242240 DO ij=1,ip1jmp1
294 163238400 iadvplus(ij,l)=0
295 ENDDO
296
297 ENDDO
298
299 c print*,'Bouclage en iip1'
300
301 c calcul des flux a gauche et a droite
302
303 c on cumule le flux correspondant a toutes les mailles dont la masse
304 c au travers de la paroi pENDant le pas de temps.
305 cprint*,'Cumule ....'
306
307
2/2
✓ Branch 0 taken 149760 times.
✓ Branch 1 taken 3840 times.
 153600 DO l=1,llm
308
2/2
✓ Branch 0 taken 153054720 times.
✓ Branch 1 taken 149760 times.
 153208320 DO ij=iip2,ip1jm-1
309 c print*,'masse(',ij,')=',masse(ij,l,iq)
310
2/2
✓ Branch 0 taken 109131512 times.
✓ Branch 1 taken 43923208 times.
 153204480 IF (u_m(ij,l).gt.0.) THEN
311 109131512 zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
312 109131512 u_mq(ij,l)=u_m(ij,l)*(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l))
313 ELSE
314 43923208 zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
315 u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l,iq)
316 43923208 & -0.5*zdum(ij,l)*dxq(ij+1,l))
317 ENDIF
318 ENDDO
319 ENDDO
320 c stop
321
322 c go to 9999
323 c detection des points ou on advecte plus que la masse de la
324 c maille
325
2/2
✓ Branch 0 taken 149760 times.
✓ Branch 1 taken 3840 times.
 153600 DO l=1,llm
326
2/2
✓ Branch 0 taken 153054720 times.
✓ Branch 1 taken 149760 times.
 153208320 DO ij=iip2,ip1jm-1
327
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 153054720 times.
 153204480 IF(zdum(ij,l).lt.0) THEN
328 ✗ iadvplus(ij,l)=1
329 ✗ u_mq(ij,l)=0.
330 ENDIF
331 ENDDO
332 ENDDO
333 cprint*,'Ok test 1'
334
2/2
✓ Branch 0 taken 149760 times.
✓ Branch 1 taken 3840 times.
 153600 DO l=1,llm
335 3840 DO ij=iip1+iip1,ip1jm,iip1
336
2/2
✓ Branch 0 taken 4492800 times.
✓ Branch 1 taken 149760 times.
 4642560 iadvplus(ij,l)=iadvplus(ij-iim,l)
337 ENDDO
338 ENDDO
339 c print*,'Ok test 2'
340
341
342 c traitement special pour le cas ou on advecte en longitude plus que le
343 c contenu de la maille.
344 c cette partie est mal vectorisee.
345
346 c calcul du nombre de maille sur lequel on advecte plus que la maille.
347
348 3840 n0=0
349
2/2
✓ Branch 0 taken 149760 times.
✓ Branch 1 taken 3840 times.
 153600 DO l=1,llm
350 149760 nl(l)=0
351
2/2
✓ Branch 0 taken 153204480 times.
✓ Branch 1 taken 149760 times.
 153354240 DO ij=iip2,ip1jm
352 153354240 nl(l)=nl(l)+iadvplus(ij,l)
353 ENDDO
354 153600 n0=n0+nl(l)
355 ENDDO
356
357
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3840 times.
 3840 IF(n0.gt.0) THEN
358 ✗ if (prt_level > 2) PRINT *,
359 ✗ $ 'Nombre de points pour lesquels on advect plus que le'
360 ✗ & ,'contenu de la maille : ',n0
361
362 ✗ DO l=1,llm
363 ✗ IF(nl(l).gt.0) THEN
364 iju=0
365 c indicage des mailles concernees par le traitement special
366 ✗ DO ij=iip2,ip1jm
367 ✗ IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
368 ✗ iju=iju+1
369 ✗ indu(iju)=ij
370 ENDIF
371 ENDDO
372 niju=iju
373 c PRINT*,'niju,nl',niju,nl(l)
374
375 c traitement des mailles
376 ✗ DO iju=1,niju
377 ✗ ij=indu(iju)
378 ✗ j=(ij-1)/iip1+1
379 ✗ zu_m=u_m(ij,l)
380 ✗ u_mq(ij,l)=0.
381 ✗ IF(zu_m.gt.0.) THEN
382 ijq=ij
383 ✗ i=ijq-(j-1)*iip1
384 c accumulation pour les mailles completements advectees
385 ✗ do while(zu_m.gt.masse(ijq,l,iq))
386 u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq)
387 ✗ & *masse(ijq,l,iq)
388 ✗ zu_m=zu_m-masse(ijq,l,iq)
389 ✗ i=mod(i-2+iim,iim)+1
390 ✗ ijq=(j-1)*iip1+i
391 ENDDO
392 c ajout de la maille non completement advectee
393 u_mq(ij,l)=u_mq(ij,l)+zu_m*
394 & (q(ijq,l,iq)+0.5*(1.-zu_m/masse(ijq,l,iq))
395 ✗ & *dxq(ijq,l))
396 ELSE
397 ✗ ijq=ij+1
398 ✗ i=ijq-(j-1)*iip1
399 c accumulation pour les mailles completements advectees
400 ✗ do while(-zu_m.gt.masse(ijq,l,iq))
401 u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
402 ✗ & *masse(ijq,l,iq)
403 ✗ zu_m=zu_m+masse(ijq,l,iq)
404 ✗ i=mod(i,iim)+1
405 ✗ ijq=(j-1)*iip1+i
406 ENDDO
407 c ajout de la maille non completement advectee
408 u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
409 ✗ & 0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l))
410 ENDIF
411 ENDDO
412 ENDIF
413 ENDDO
414 ENDIF ! n0.gt.0
415 9999 continue
416
417
418 c bouclage en latitude
419 cprint*,'cvant bouclage en latitude'
420
2/2
✓ Branch 0 taken 149760 times.
✓ Branch 1 taken 3840 times.
 153600 DO l=1,llm
421 3840 DO ij=iip1+iip1,ip1jm,iip1
422
2/2
✓ Branch 0 taken 4492800 times.
✓ Branch 1 taken 149760 times.
 4642560 u_mq(ij,l)=u_mq(ij-iim,l)
423 ENDDO
424 ENDDO
425
426 ! CRisi: appel récursif de l'advection sur les fils.
427 ! Il faut faire ça avant d'avoir mis à jour q et masse
428 !write(*,*) 'vlsplt 326: iq,nqfils(iq)=',iq,nqfils(iq)
429
430
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3840 times.
 3840 if (nqdesc(iq).gt.0) then
431 ✗ do ifils=1,nqdesc(iq)
432 ✗ iq2=iqfils(ifils,iq)
433 ✗ DO l=1,llm
434 ✗ DO ij=iip2,ip1jm
435 ! On a besoin de q et masse seulement entre iip2 et ip1jm
436 ✗ masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
437 ✗ Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
438 enddo
439 enddo
440 enddo !do ifils=1,nqdesc(iq)
441 ✗ do ifils=1,nqfils(iq)
442 ✗ iq2=iqfils(ifils,iq)
443 ✗ call vlx(Ratio,pente_max,masseq,u_mq,iq2)
444 enddo !do ifils=1,nqfils(iq)
445 endif !if (nqfils(iq).gt.0) then
446 ! end CRisi
447
448
449 c calcul des tENDances
450
451
2/2
✓ Branch 0 taken 149760 times.
✓ Branch 1 taken 3840 times.
 153600 DO l=1,llm
452
2/2
✓ Branch 0 taken 153054720 times.
✓ Branch 1 taken 149760 times.
 153204480 DO ij=iip2+1,ip1jm
453 153054720 new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l)
454 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
455 & u_mq(ij-1,l)-u_mq(ij,l))
456 153054720 & /new_m
457 153204480 masse(ij,l,iq)=new_m
458 ENDDO
459 c ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
460 3840 DO ij=iip1+iip1,ip1jm,iip1
461 4642560 q(ij-iim,l,iq)=q(ij,l,iq)
462
2/2
✓ Branch 0 taken 4492800 times.
✓ Branch 1 taken 149760 times.
 4642560 masse(ij-iim,l,iq)=masse(ij,l,iq)
463 ENDDO
464 ENDDO
465
466 ! retablir les fils en rapport de melange par rapport a l'air:
467 ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio
468 ! puis on boucle en longitude
469
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3840 times.
 3840 if (nqdesc(iq).gt.0) then
470 ✗ do ifils=1,nqdesc(iq)
471 ✗ iq2=iqfils(ifils,iq)
472 ✗ DO l=1,llm
473 ✗ DO ij=iip2+1,ip1jm
474 ✗ q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)
475 enddo
476 ✗ DO ij=iip1+iip1,ip1jm,iip1
477 ✗ q(ij-iim,l,iq2)=q(ij,l,iq2)
478 enddo ! DO ij=ijb+iip1-1,ije,iip1
479 enddo !DO l=1,llm
480 enddo !do ifils=1,nqdesc(iq)
481 endif !if (nqfils(iq).gt.0) then
482
483 c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
484 c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
485
486
487 3840 RETURN
488 END
489 3840 RECURSIVE SUBROUTINE vly(q,pente_max,masse,masse_adv_v,iq)
490 USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
491 c
492 c Auteurs: P.Le Van, F.Hourdin, F.Forget
493 c
494 c ********************************************************************
495 c Shema d'advection " pseudo amont " .
496 c ********************************************************************
497 c q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
498 c dq sont des arguments de sortie pour le s-pg ....
499 c
500 c
501 c --------------------------------------------------------------------
502 USE comconst_mod, ONLY: pi
503 IMPLICIT NONE
504 c
505 include "dimensions.h"
506 include "paramet.h"
507 include "comgeom.h"
508 c
509 c
510 c Arguments:
511 c ----------
512 REAL masse(ip1jmp1,llm,nqtot),pente_max
513 REAL masse_adv_v( ip1jm,llm)
514 REAL q(ip1jmp1,llm,nqtot), dq( ip1jmp1,llm)
515 INTEGER iq ! CRisi
516 c
517 c Local
518 c ---------
519 c
520 INTEGER i,ij,l
521 c
522 REAL airej2,airejjm,airescb(iim),airesch(iim)
523 REAL dyq(ip1jmp1,llm),dyqv(ip1jm),zdvm(ip1jmp1,llm)
524 REAL adyqv(ip1jm),dyqmax(ip1jmp1)
525 REAL qbyv(ip1jm,llm)
526
527 REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
528 c REAL newq,oldmasse
529 Logical extremum,first,testcpu
530 REAL temps0,temps1,temps2,temps3,temps4,temps5,second
531 SAVE temps0,temps1,temps2,temps3,temps4,temps5
532 SAVE first,testcpu
533
534 REAL convpn,convps,convmpn,convmps
535 real massepn,masseps,qpn,qps
536 REAL sinlon(iip1),sinlondlon(iip1)
537 REAL coslon(iip1),coslondlon(iip1)
538 SAVE sinlon,coslon,sinlondlon,coslondlon
539 SAVE airej2,airejjm
540
541 7680 REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
542 INTEGER ifils,iq2 ! CRisi
543
544 c
545 c
546 REAL SSUM
547
548 DATA first,testcpu/.true.,.false./
549 DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
550
551 !write(*,*) 'vly 578: entree, iq=',iq
552
553
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3839 times.
 3840 IF(first) THEN
554 1 PRINT*,'Shema Amont nouveau appele dans Vanleer '
555 1 first=.false.
556
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 1 times.
 33 do i=2,iip1
557 32 coslon(i)=cos(rlonv(i))
558 32 sinlon(i)=sin(rlonv(i))
559 32 coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
560 33 sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
561 ENDDO
562 1 coslon(1)=coslon(iip1)
563 1 coslondlon(1)=coslondlon(iip1)
564 1 sinlon(1)=sinlon(iip1)
565 1 sinlondlon(1)=sinlondlon(iip1)
566 1 airej2 = SSUM( iim, aire(iip2), 1 )
567 1 airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
568 ENDIF
569
570 c
571 cPRINT*,'CALCUL EN LATITUDE'
572
573
2/2
✓ Branch 0 taken 149760 times.
✓ Branch 1 taken 3840 times.
 153600 DO l = 1, llm
574 c
575 c --------------------------------
576 c CALCUL EN LATITUDE
577 c --------------------------------
578
579 c On commence par calculer la valeur du traceur moyenne sur le premier cercle
580 c de latitude autour du pole (qpns pour le pole nord et qpsn pour
581 c le pole nord) qui sera utilisee pour evaluer les pentes au pole.
582
583
2/2
✓ Branch 0 taken 4792320 times.
✓ Branch 1 taken 149760 times.
 4942080 DO i = 1, iim
584 4792320 airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
585 4942080 airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
586 ENDDO
587 149760 qpns = SSUM( iim, airescb ,1 ) / airej2
588 149760 qpsn = SSUM( iim, airesch ,1 ) / airejjm
589
590 c calcul des pentes aux points v
591
592
2/2
✓ Branch 0 taken 158146560 times.
✓ Branch 1 taken 149760 times.
 158296320 DO ij=1,ip1jm
593 158146560 dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
594 158296320 adyqv(ij)=abs(dyqv(ij))
595 ENDDO
596
597 c calcul des pentes aux points scalaires
598
599
2/2
✓ Branch 0 taken 153204480 times.
✓ Branch 1 taken 149760 times.
 153354240 DO ij=iip2,ip1jm
600 153204480 dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
601 153204480 dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
602 153354240 dyqmax(ij)=pente_max*dyqmax(ij)
603 ENDDO
604
605 c calcul des pentes aux poles
606
607
2/2
✓ Branch 0 taken 4942080 times.
✓ Branch 1 taken 149760 times.
 5091840 DO ij=1,iip1
608 4942080 dyq(ij,l)=qpns-q(ij+iip1,l,iq)
609 5091840 dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
610 ENDDO
611
612 c filtrage de la derivee
613 dyn1=0.
614 dys1=0.
615 dyn2=0.
616 dys2=0.
617
2/2
✓ Branch 0 taken 149760 times.
✓ Branch 1 taken 4792320 times.
 4942080 DO ij=1,iim
618 4792320 dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
619 4792320 dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
620 4792320 dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
621 4942080 dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
622 ENDDO
623
2/2
✓ Branch 0 taken 4942080 times.
✓ Branch 1 taken 149760 times.
 5091840 DO ij=1,iip1
624 4942080 dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
625 5091840 dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
626 ENDDO
627
628 c calcul des pentes limites aux poles
629
630 goto 8888
631 fn=1.
632 fs=1.
633 DO ij=1,iim
634 IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
635 fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
636 ENDIF
637 IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
638 fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
639 ENDIF
640 ENDDO
641 DO ij=1,iip1
642 dyq(ij,l)=fn*dyq(ij,l)
643 dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
644 ENDDO
645 8888 continue
646
2/2
✓ Branch 0 taken 4942080 times.
✓ Branch 1 taken 149760 times.
 5091840 DO ij=1,iip1
647 4942080 dyq(ij,l)=0.
648 5091840 dyq(ip1jm+ij,l)=0.
649 ENDDO
650
651 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
652 C En memoire de dIFferents tests sur la
653 C limitation des pentes aux poles.
654 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
655 C PRINT*,dyq(1)
656 C PRINT*,dyqv(iip1+1)
657 C appn=abs(dyq(1)/dyqv(iip1+1))
658 C PRINT*,dyq(ip1jm+1)
659 C PRINT*,dyqv(ip1jm-iip1+1)
660 C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
661 C DO ij=2,iim
662 C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
663 C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
664 C ENDDO
665 C appn=min(pente_max/appn,1.)
666 C apps=min(pente_max/apps,1.)
667 C
668 C
669 C cas ou on a un extremum au pole
670 C
671 C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
672 C & appn=0.
673 C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
674 C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
675 C & apps=0.
676 C
677 C limitation des pentes aux poles
678 C DO ij=1,iip1
679 C dyq(ij)=appn*dyq(ij)
680 C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
681 C ENDDO
682 C
683 C test
684 C DO ij=1,iip1
685 C dyq(iip1+ij)=0.
686 C dyq(ip1jm+ij-iip1)=0.
687 C ENDDO
688 C DO ij=1,ip1jmp1
689 C dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
690 C ENDDO
691 C
692 C changement 10 07 96
693 C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
694 C & THEN
695 C DO ij=1,iip1
696 C dyqmax(ij)=0.
697 C ENDDO
698 C ELSE
699 C DO ij=1,iip1
700 C dyqmax(ij)=pente_max*abs(dyqv(ij))
701 C ENDDO
702 C ENDIF
703 C
704 C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
705 C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
706 C &THEN
707 C DO ij=ip1jm+1,ip1jmp1
708 C dyqmax(ij)=0.
709 C ENDDO
710 C ELSE
711 C DO ij=ip1jm+1,ip1jmp1
712 C dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
713 C ENDDO
714 C ENDIF
715 C fin changement 10 07 96
716 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
717
718 c calcul des pentes limitees
719
720
2/2
✓ Branch 0 taken 153204480 times.
✓ Branch 1 taken 149760 times.
 153358080 DO ij=iip2,ip1jm
721
2/2
✓ Branch 0 taken 66366910 times.
✓ Branch 1 taken 86837570 times.
 153354240 IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
722 66366910 dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
723 ELSE
724 86837570 dyq(ij,l)=0.
725 ENDIF
726 ENDDO
727
728 ENDDO
729
730 !write(*,*) 'vly 756'
731
2/2
✓ Branch 0 taken 149760 times.
✓ Branch 1 taken 3840 times.
 153600 DO l=1,llm
732
2/2
✓ Branch 0 taken 158146560 times.
✓ Branch 1 taken 149760 times.
 158300160 DO ij=1,ip1jm
733
2/2
✓ Branch 0 taken 80232456 times.
✓ Branch 1 taken 77914104 times.
 158146560 IF(masse_adv_v(ij,l).gt.0) THEN
734 qbyv(ij,l)=q(ij+iip1,l,iq)+dyq(ij+iip1,l)*
735 , 0.5*(1.-masse_adv_v(ij,l)
736 80232456 , /masse(ij+iip1,l,iq))
737 ELSE
738 qbyv(ij,l)=q(ij,l,iq)-dyq(ij,l)*
739 , 0.5*(1.+masse_adv_v(ij,l)
740 77914104 , /masse(ij,l,iq))
741 ENDIF
742 158296320 qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
743 ENDDO
744 ENDDO
745
746 ! CRisi: appel récursif de l'advection sur les fils.
747 ! Il faut faire ça avant d'avoir mis à jour q et masse
748 !write(*,*) 'vly 689: iq,nqfils(iq)=',iq,nqfils(iq)
749
750
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3840 times.
 3840 if (nqfils(iq).gt.0) then
751 ✗ do ifils=1,nqdesc(iq)
752 ✗ iq2=iqfils(ifils,iq)
753 ✗ DO l=1,llm
754 ✗ DO ij=1,ip1jmp1
755 ! attention, chaque fils doit avoir son masseq, sinon, le 1er
756 ! fils ecrase le masseq de ses freres.
757 ✗ masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
758 ✗ Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
759 enddo
760 enddo
761 enddo !do ifils=1,nqdesc(iq)
762
763 ✗ do ifils=1,nqfils(iq)
764 ✗ iq2=iqfils(ifils,iq)
765 ✗ call vly(Ratio,pente_max,masseq,qbyv,iq2)
766 enddo !do ifils=1,nqfils(iq)
767 endif !if (nqfils(iq).gt.0) then
768
769
2/2
✓ Branch 0 taken 149760 times.
✓ Branch 1 taken 3840 times.
 153600 DO l=1,llm
770
2/2
✓ Branch 0 taken 153204480 times.
✓ Branch 1 taken 149760 times.
 153354240 DO ij=iip2,ip1jm
771 newmasse=masse(ij,l,iq)
772 153204480 & +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
773 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l)
774 153204480 & -qbyv(ij-iip1,l))/newmasse
775 153354240 masse(ij,l,iq)=newmasse
776 ENDDO
777 c.-. ancienne version
778 c convpn=SSUM(iim,qbyv(1,l),1)/apoln
779 c convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
780
781 149760 convpn=SSUM(iim,qbyv(1,l),1)
782 149760 convmpn=ssum(iim,masse_adv_v(1,l),1)
783 149760 massepn=ssum(iim,masse(1,l,iq),1)
784 qpn=0.
785
2/2
✓ Branch 0 taken 4792320 times.
✓ Branch 1 taken 149760 times.
 4942080 do ij=1,iim
786 4942080 qpn=qpn+masse(ij,l,iq)*q(ij,l,iq)
787 enddo
788 149760 qpn=(qpn+convpn)/(massepn+convmpn)
789
2/2
✓ Branch 0 taken 4942080 times.
✓ Branch 1 taken 149760 times.
 5091840 do ij=1,iip1
790 5091840 q(ij,l,iq)=qpn
791 enddo
792
793 c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
794 c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
795
796 149760 convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
797 149760 convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
798 149760 masseps=ssum(iim, masse(ip1jm+1,l,iq),1)
799 qps=0.
800
2/2
✓ Branch 0 taken 4792320 times.
✓ Branch 1 taken 149760 times.
 4942080 do ij = ip1jm+1,ip1jmp1-1
801 4942080 qps=qps+masse(ij,l,iq)*q(ij,l,iq)
802 enddo
803 149760 qps=(qps+convps)/(masseps+convmps)
804
2/2
✓ Branch 0 taken 4942080 times.
✓ Branch 1 taken 149760 times.
 5095680 do ij=ip1jm+1,ip1jmp1
805 5091840 q(ij,l,iq)=qps
806 enddo
807
808 c.-. fin ancienne version
809
810 c._. nouvelle version
811 c convpn=SSUM(iim,qbyv(1,l),1)
812 c convmpn=ssum(iim,masse_adv_v(1,l),1)
813 c oldmasse=ssum(iim,masse(1,l),1)
814 c newmasse=oldmasse+convmpn
815 c newq=(q(1,l)*oldmasse+convpn)/newmasse
816 c newmasse=newmasse/apoln
817 c DO ij = 1,iip1
818 c q(ij,l)=newq
819 c masse(ij,l,iq)=newmasse*aire(ij)
820 c ENDDO
821 c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
822 c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
823 c oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
824 c newmasse=oldmasse+convmps
825 c newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
826 c newmasse=newmasse/apols
827 c DO ij = ip1jm+1,ip1jmp1
828 c q(ij,l)=newq
829 c masse(ij,l,iq)=newmasse*aire(ij)
830 c ENDDO
831 c._. fin nouvelle version
832 ENDDO
833
834 ! retablir les fils en rapport de melange par rapport a l'air:
835
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3840 times.
 3840 if (nqfils(iq).gt.0) then
836 ✗ do ifils=1,nqdesc(iq)
837 ✗ iq2=iqfils(ifils,iq)
838 ✗ DO l=1,llm
839 ✗ DO ij=1,ip1jmp1
840 ✗ q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)
841 enddo
842 enddo
843 enddo !do ifils=1,nqdesc(iq)
844 endif !if (nqfils(iq).gt.0) then
845
846 !write(*,*) 'vly 853: sortie'
847
848 3840 RETURN
849 END
850 2400 RECURSIVE SUBROUTINE vlz(q,pente_max,masse,w,iq)
851 USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
852 c
853 c Auteurs: P.Le Van, F.Hourdin, F.Forget
854 c
855 c ********************************************************************
856 c Shema d'advection " pseudo amont " .
857 c ********************************************************************
858 c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
859 c dq sont des arguments de sortie pour le s-pg ....
860 c
861 c
862 c --------------------------------------------------------------------
863 IMPLICIT NONE
864 c
865 include "dimensions.h"
866 include "paramet.h"
867 c
868 c
869 c Arguments:
870 c ----------
871 REAL masse(ip1jmp1,llm,nqtot),pente_max
872 REAL q(ip1jmp1,llm,nqtot)
873 REAL w(ip1jmp1,llm+1)
874 INTEGER iq
875 c
876 c Local
877 c ---------
878 c
879 INTEGER i,ij,l,j,ii
880 c
881 REAL wq(ip1jmp1,llm+1),newmasse
882
883 REAL dzq(ip1jmp1,llm),dzqw(ip1jmp1,llm),adzqw(ip1jmp1,llm),dzqmax
884 REAL sigw
885
886 2400 REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
887 INTEGER ifils,iq2 ! CRisi
888
889 LOGICAL testcpu
890 SAVE testcpu
891
892 REAL temps0,temps1,temps2,temps3,temps4,temps5,second
893 SAVE temps0,temps1,temps2,temps3,temps4,temps5
894 REAL SSUM
895
896 DATA testcpu/.false./
897 DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
898
899 c On oriente tout dans le sens de la pression c'est a dire dans le
900 c sens de W
901
902 !write(*,*) 'vlz 923: entree'
903
904
2/2
✓ Branch 0 taken 91200 times.
✓ Branch 1 taken 2400 times.
 93600 DO l=2,llm
905
2/2
✓ Branch 0 taken 99316800 times.
✓ Branch 1 taken 91200 times.
 99410400 DO ij=1,ip1jmp1
906 99316800 dzqw(ij,l)=q(ij,l-1,iq)-q(ij,l,iq)
907 99408000 adzqw(ij,l)=abs(dzqw(ij,l))
908 ENDDO
909 ENDDO
910
911
2/2
✓ Branch 0 taken 88800 times.
✓ Branch 1 taken 2400 times.
 91200 DO l=2,llm-1
912
2/2
✓ Branch 0 taken 96703200 times.
✓ Branch 1 taken 88800 times.
 96794400 DO ij=1,ip1jmp1
913
2/2
✓ Branch 0 taken 55269341 times.
✓ Branch 1 taken 41433859 times.
 96703200 IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN
914 55269341 dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
915 ELSE
916 41433859 dzq(ij,l)=0.
917 ENDIF
918 96703200 dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
919 96792000 dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
920 ENDDO
921 ENDDO
922
923 !write(*,*) 'vlz 954'
924
2/2
✓ Branch 0 taken 2613600 times.
✓ Branch 1 taken 2400 times.
 2616000 DO ij=1,ip1jmp1
925 2613600 dzq(ij,1)=0.
926 2616000 dzq(ij,llm)=0.
927 ENDDO
928
929 c ---------------------------------------------------------------
930 c .... calcul des termes d'advection verticale .......
931 c ---------------------------------------------------------------
932
933 c calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq
934
935 !write(*,*) 'vlz 969'
936
2/2
✓ Branch 0 taken 91200 times.
✓ Branch 1 taken 2400 times.
 93600 DO l = 1,llm-1
937
2/2
✓ Branch 0 taken 99316800 times.
✓ Branch 1 taken 91200 times.
 99410400 do ij = 1,ip1jmp1
938
2/2
✓ Branch 0 taken 49793035 times.
✓ Branch 1 taken 49523765 times.
 99408000 IF(w(ij,l+1).gt.0.) THEN
939 49793035 sigw=w(ij,l+1)/masse(ij,l+1,iq)
940 wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1,iq)
941 49793035 & +0.5*(1.-sigw)*dzq(ij,l+1))
942 ELSE
943 49523765 sigw=w(ij,l+1)/masse(ij,l,iq)
944 49523765 wq(ij,l+1)=w(ij,l+1)*(q(ij,l,iq)-0.5*(1.+sigw)*dzq(ij,l))
945 ENDIF
946 ENDDO
947 ENDDO
948
949
2/2
✓ Branch 0 taken 2613600 times.
✓ Branch 1 taken 2400 times.
 2616000 DO ij=1,ip1jmp1
950 2613600 wq(ij,llm+1)=0.
951 2616000 wq(ij,1)=0.
952 ENDDO
953
954 ! CRisi: appel récursif de l'advection sur les fils.
955 ! Il faut faire ça avant d'avoir mis à jour q et masse
956 !write(*,*) 'vlsplt 942: iq,nqfils(iq)=',iq,nqfils(iq)
957
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2400 times.
 2400 if (nqfils(iq).gt.0) then
958 ✗ do ifils=1,nqdesc(iq)
959 ✗ iq2=iqfils(ifils,iq)
960 ✗ DO l=1,llm
961 ✗ DO ij=1,ip1jmp1
962 ✗ masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
963 ✗ Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
964 enddo
965 enddo
966 enddo !do ifils=1,nqdesc(iq)
967
968 ✗ do ifils=1,nqfils(iq)
969 ✗ iq2=iqfils(ifils,iq)
970 ✗ call vlz(Ratio,pente_max,masseq,wq,iq2)
971 enddo !do ifils=1,nqfils(iq)
972 endif !if (nqfils(iq).gt.0) then
973 ! end CRisi
974
975
2/2
✓ Branch 0 taken 93600 times.
✓ Branch 1 taken 2400 times.
 96000 DO l=1,llm
976
2/2
✓ Branch 0 taken 101930400 times.
✓ Branch 1 taken 93600 times.
 102026400 DO ij=1,ip1jmp1
977 101930400 newmasse=masse(ij,l,iq)+w(ij,l+1)-w(ij,l)
978 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+wq(ij,l+1)-wq(ij,l))
979 101930400 & /newmasse
980 102024000 masse(ij,l,iq)=newmasse
981 ENDDO
982 ENDDO
983
984 ! retablir les fils en rapport de melange par rapport a l'air:
985
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2400 times.
 2400 if (nqfils(iq).gt.0) then
986 ✗ do ifils=1,nqdesc(iq)
987 ✗ iq2=iqfils(ifils,iq)
988 ✗ DO l=1,llm
989 ✗ DO ij=1,ip1jmp1
990 ✗ q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)
991 enddo
992 enddo
993 enddo !do ifils=1,nqdesc(iq)
994 endif !if (nqfils(iq).gt.0) then
995 !write(*,*) 'vlsplt 1032'
996
997 2400 RETURN
998 END
999 c SUBROUTINE minmaxq(zq,qmin,qmax,comment)
1000 c
1001 c#include "dimensions.h"
1002 c#include "paramet.h"
1003
1004 c CHARACTER*(*) comment
1005 c real qmin,qmax
1006 c real zq(ip1jmp1,llm)
1007
1008 c INTEGER jadrs(ip1jmp1), jbad, k, i
1009
1010
1011 c DO k = 1, llm
1012 c jbad = 0
1013 c DO i = 1, ip1jmp1
1014 c IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
1015 c jbad = jbad + 1
1016 c jadrs(jbad) = i
1017 c ENDIF
1018 c ENDDO
1019 c IF (jbad.GT.0) THEN
1020 c PRINT*, comment
1021 c DO i = 1, jbad
1022 cc PRINT*, "i,k,zq=", jadrs(i),k,zq(jadrs(i),k)
1023 c ENDDO
1024 c ENDIF
1025 c ENDDO
1026
1027 c return
1028 c end
1029 ✗ subroutine minmaxq(zq,qmin,qmax,comment)
1030
1031 !-----------------------------------------------------------------------
1032 ! INCLUDE 'dimensions.h'
1033 !
1034 ! dimensions.h contient les dimensions du modele
1035 ! ndm est tel que iim=2**ndm
1036 !-----------------------------------------------------------------------
1037
1038 INTEGER iim,jjm,llm,ndm
1039
1040 PARAMETER (iim= 32,jjm=32,llm=39,ndm=1)
1041
1042 !-----------------------------------------------------------------------
1043 !
1044 ! $Header$
1045 !
1046 !
1047 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
1048 ! veillez n'utiliser que des ! pour les commentaires
1049 ! et bien positionner les & des lignes de continuation
1050 ! (les placer en colonne 6 et en colonne 73)
1051 !
1052 !
1053 !-----------------------------------------------------------------------
1054 ! INCLUDE 'paramet.h'
1055
1056 INTEGER iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
1057 INTEGER kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
1058 INTEGER ijmllm,mvar
1059 INTEGER jcfil,jcfllm
1060
1061 PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3 &
1062 & ,jjp1=jjm+1-1/jjm)
1063 PARAMETER( llmp1 = llm+1, llmp2 = llm+2, llmm1 = llm-1 )
1064 PARAMETER( kftd = iim/2 -ndm )
1065 PARAMETER( ip1jm = iip1*jjm, ip1jmp1= iip1*jjp1 )
1066 PARAMETER( ip1jmi1= ip1jm - iip1 )
1067 PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
1068 PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm )
1069 PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
1070
1071 !-----------------------------------------------------------------------
1072
1073 character*20 comment
1074 real qmin,qmax
1075 real zq(ip1jmp1,llm)
1076 real zzq(iip1,jjp1,llm)
1077
1078 integer imin,jmin,lmin,ijlmin
1079 integer imax,jmax,lmax,ijlmax
1080
1081 integer ismin,ismax
1082
1083 ✗ return
1084 9999 format(a20,' q(',i3,',',i2,',',i2,')=',e12.5,e12.5)
1085 end
1086
1087
1088
1089