GCC Code Coverage Report


Directory: ./
File: phys/cvltr_scav.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 236 0.0%
Branches: 0 232 0.0%

Line Branch Exec Source
1 !
2 ! $Id $
3 !
4 SUBROUTINE cvltr_scav(pdtime, da, phi,phi2,d1a,dam, mpIN,epIN, &
5 sigd,sij,wght_cvfd,clw,elij,epmlmMm,eplaMm, &
6 pmflxrIN,pmflxsIN,ev,te,wdtrainA,wdtrainM, &
7 paprs,it,tr,upd,dnd,inb,icb, &
8 ccntrAA_3d,ccntrENV_3d,coefcoli_3d, &
9 dtrcv,trsptd,dtrSscav,dtrsat,dtrUscav,qDi,qPr, &
10 qPa,qMel,qTrdi,dtrcvMA,Mint, &
11 zmfd1a,zmfphi2,zmfdam)
12 !
13 USE IOIPSL
14 USE dimphy
15 USE infotrac_phy, ONLY : nbtr,tname
16 IMPLICIT NONE
17 !=====================================================================
18 ! Objet : convection des traceurs / KE
19 ! Auteurs: M-A Filiberti and J-Y Grandpeix
20 ! modifiee par R Pilon : lessivage des traceurs / KE
21 !=====================================================================
22
23 include "YOMCST.h"
24 include "YOECUMF.h"
25 include "conema3.h"
26 include "chem.h"
27
28 ! Entree
29 REAL,INTENT(IN) :: pdtime
30 REAL,DIMENSION(klon,klev),INTENT(IN) :: da
31 REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi
32 ! RomP
33 REAL,DIMENSION(klon,klev),INTENT(IN) :: d1a,dam ! matrices pour simplifier
34 REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi2 ! l'ecriture des tendances
35 !
36 REAL,DIMENSION(klon,klev),INTENT(IN) :: mpIN
37 REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression aux 1/2 couches (bas en haut)
38 INTEGER,INTENT(IN) :: it ! numero du traceur
39 REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr ! q de traceur (bas en haut)
40 REAL,DIMENSION(klon,klev),INTENT(IN) :: upd ! saturated updraft mass flux
41 REAL,DIMENSION(klon,klev),INTENT(IN) :: dnd ! saturated downdraft mass flux
42 !
43 REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainA ! masses precipitantes de l'asc adiab
44 REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainM ! masses precipitantes des melanges
45 !JE REAL,DIMENSION(klon,klev),INTENT(IN) :: pmflxrIN ! vprecip: eau
46 REAL,DIMENSION(klon,klev+1),INTENT(IN) :: pmflxrIN ! vprecip: eau
47 !JE REAL,DIMENSION(klon,klev),INTENT(IN) :: pmflxsIN ! vprecip: neige
48 REAL,DIMENSION(klon,klev+1),INTENT(IN) :: pmflxsIN ! vprecip: neige
49 REAL,DIMENSION(klon,klev),INTENT(IN) :: ev ! evaporation cv30_routine
50 REAL,DIMENSION(klon,klev),INTENT(IN) :: epIN
51 REAL,DIMENSION(klon,klev),INTENT(IN) :: te
52 REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: sij ! fraction dair de lenv
53 REAL,DIMENSION(klon,klev),INTENT(IN) :: wght_cvfd ! weights of the layers feeding convection
54 REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: elij ! contenu en eau condensée spécifique/conc deau condensée massique
55 REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: epmlmMm ! eau condensee precipitee dans mel masse dair sat
56 REAL,DIMENSION(klon,klev),INTENT(IN) :: eplaMm ! eau condensee precipitee dans aa masse dair sat
57
58 REAL,DIMENSION(klon,klev),INTENT(IN) :: clw ! contenu en eau condensée dans lasc adiab
59 REAL,DIMENSION(klon),INTENT(IN) :: sigd
60 INTEGER,DIMENSION(klon),INTENT(IN) :: icb,inb
61 !
62 REAL,DIMENSION(klon,klev),INTENT(IN) :: ccntrAA_3d
63 REAL,DIMENSION(klon,klev),INTENT(IN) :: ccntrENV_3d
64 REAL,DIMENSION(klon,klev),INTENT(IN) :: coefcoli_3d
65 !
66 ! Sortie
67 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrcv ! tendance totale (bas en haut)
68 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrcvMA ! M-A Filiberti
69 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: trsptd ! tendance du transport
70 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrSscav ! tendance du lessivage courant sat
71 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrsat ! tendance trsp+sat scav
72 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrUscav ! tendance du lessivage courant unsat
73 !
74 ! Variables locales
75 INTEGER :: i,j,k
76 REAL,DIMENSION(klon,klev) :: dxpres ! difference de pression entre niveau (j+1) et (j)
77 REAL :: pdtimeRG ! pas de temps * gravite
78 REAL,DIMENSION(klon,nbtr) :: qfeed ! tracer concentration feeding convection
79 ! variables pour les courants satures
80 REAL,DIMENSION(klon,klev,klev) :: zmd
81 REAL,DIMENSION(klon,klev,klev) :: za
82 REAL,DIMENSION(klon,klev,nbtr) :: zmfd,zmfa
83 REAL,DIMENSION(klon,klev,nbtr) :: zmfp,zmfu
84
85 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfd1a
86 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfdam
87 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfphi2
88
89 ! RomP ! les variables sont nettoyees des valeurs aberrantes
90 REAL,DIMENSION(klon,klev) :: Pa, Pm ! pluie AA et mélanges, var temporaire
91 REAL,DIMENSION(klon,klev) :: pmflxs,pmflxr ! pmflxrIN,pmflxsIN sans valeur aberante
92 REAL,DIMENSION(klon,klev) :: mp ! flux de masse
93 REAL,DIMENSION(klon,klev) :: ep ! fraction d'eau convertie en precipitation
94 REAL,DIMENSION(klon,klev) :: evap ! evaporation : variable temporaire
95 REAL,DIMENSION(klon,klev) :: rho !environmental density
96
97 REAL,DIMENSION(klon,klev) :: kappa ! denominateur du au calcul de la matrice
98 ! pour obtenir qd et qp
99 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qTrdi ! traceurs descente air insature transport MA
100 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qDi ! traceurs descente insaturees
101 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qPr ! traceurs colonne precipitante
102 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qPa ! traceurs dans les precip issues lasc. adiab.
103 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qMel ! traceurs dans les precip issues des melanges
104 REAL,DIMENSION(klon,klev,nbtr) :: qMeltmp ! variable temporaire
105 REAL,DIMENSION(klon,klev,nbtr) :: qpmMint
106 REAL,DIMENSION(klon,klev),INTENT(OUT) :: Mint
107 ! tendances
108 REAL :: tdcvMA ! terme de transport de traceur (schema Marie Angele)
109 REAL :: trsptrac ! terme de transport de traceur par l'air
110 REAL :: scavtrac ! terme de lessivage courant sature
111 REAL :: uscavtrac ! terme de lessivage courant insature
112 ! impaction
113 !!! Correction apres discussion Romain P. / Olivier B.
114 !!! REAL,PARAMETER :: rdrop=2.5e-3 ! rayon des gouttes d'eau
115 REAL,PARAMETER :: rdrop=1.e-3 ! rayon des gouttes d'eau
116 !!!
117 REAL,DIMENSION(klon,klev) :: imp ! coefficient d'impaction
118 !
119 LOGICAL,DIMENSION(klon,klev) :: NO_precip
120 ! var tmp tests
121 REAL :: conserv
122 real :: conservMA
123
124 !jyg<
125 !! ! ======================================================
126 !! ! calcul de l'impaction
127 !! ! ======================================================
128 !!
129 !! ! impaction sur la surface de la colonne de la descente insaturee
130 !! ! On prend la moyenne des precip entre le niveau i+1 et i
131 !! ! I=3/4* (P(1+1)+P(i))/2 / (sigd*r*rho_l)
132 !! ! 1000kg/m3= densite de l'eau
133 !! ! 0.75e-3 = 3/4 /1000
134 !! ! Par la suite, I est tout le temps multiplie par sig_d pour avoir l'impaction sur la surface de la maille
135 !!!! ! on le neglige ici pour simplifier le code
136 !!
137 !! DO j=1,klev-1
138 !! DO i=1,klon
139 !! imp(i,j) = coefcoli_3d(i,j)*0.75e-3/rdrop *&
140 !! 0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
141 !! ENDDO
142 !! ENDDO
143 !>jyg
144 !
145 ! initialisation pour flux de traceurs, td et autre
146 !
147 trsptrac = 0.
148 scavtrac = 0.
149 uscavtrac = 0.
150 qfeed(:,it) = 0. !RL
151 DO j=1,klev
152 DO i=1,klon
153 zmfd(i,j,it)=0.
154 zmfa(i,j,it)=0.
155 zmfu(i,j,it)=0.
156 zmfp(i,j,it)=0.
157 zmfphi2(i,j,it)=0.
158 zmfd1a(i,j,it)=0.
159 zmfdam(i,j,it)=0.
160 qDi(i,j,it)=0.
161 qPr(i,j,it)=0.
162 qPa(i,j,it)=0.
163 qMel(i,j,it)=0.
164 qMeltmp(i,j,it)=0.
165 qTrdi(i,j,it)=0.
166 kappa(i,j)=0.
167 trsptd(i,j,it)=0.
168 dtrsat(i,j,it)=0.
169 dtrSscav(i,j,it)=0.
170 dtrUscav(i,j,it)=0.
171 dtrcv(i,j,it)=0.
172 dtrcvMA(i,j,it)=0.
173 evap(i,j)=0.
174 dxpres(i,j)=0.
175 qpmMint(i,j,it)=0.
176 Mint(i,j)=0.
177 END DO
178 END DO
179
180 ! suppression des valeurs très faibles (~1e-320)
181 ! multiplication de levaporation pour lavoir par unite de temps
182 ! et par unite de surface de la maille
183 ! -> cv30_unsat : evap : masse evaporee/s/(m2 de la descente)
184 DO j=1,klev
185 DO i=1,klon
186 IF(ev(i,j).lt.1.e-16) THEN
187 evap(i,j)=0.
188 ELSE
189 evap(i,j)=ev(i,j)*sigd(i)
190 ENDIF
191 END DO
192 END DO
193
194 DO j=1,klev
195 DO i=1,klon
196 IF(j.LT.klev) THEN
197 IF(epIN(i,j).LT.1.e-32) THEN
198 ep(i,j)=0.
199 ELSE
200 ep(i,j)=epIN(i,j)
201 ENDIF
202 ELSE
203 ep(i,j)=epmax
204 ENDIF
205 IF(mpIN(i,j).LT.1.e-32) THEN
206 mp(i,j)=0.
207 ELSE
208 mp(i,j)=mpIN(i,j)
209 ENDIF
210 IF(pmflxsIN(i,j).LT.1.e-32) THEN
211 pmflxs(i,j)=0.
212 ELSE
213 pmflxs(i,j)=pmflxsIN(i,j)
214 ENDIF
215 IF(pmflxrIN(i,j).LT.1.e-32) THEN
216 pmflxr(i,j)=0.
217 ELSE
218 pmflxr(i,j)=pmflxrIN(i,j)
219 ENDIF
220 IF(wdtrainA(i,j).LT.1.e-32) THEN
221 Pa(i,j)=0.
222 ELSE
223 Pa(i,j)=wdtrainA(i,j)
224 ENDIF
225 IF(wdtrainM(i,j).LT.1.e-32) THEN
226 Pm(i,j)=0.
227 ELSE
228 Pm(i,j)=wdtrainM(i,j)
229 ENDIF
230 END DO
231 END DO
232
233 !==========================================
234 DO j = klev-1,1,-1
235 DO i = 1,klon
236 NO_precip(i,j) = (pmflxr(i,j+1)+pmflxs(i,j+1)).LT.1.e-10&
237 .AND.Pa(i,j).LT.1.e-10.AND.Pm(i,j).LT.1.e-10
238 END DO
239 END DO
240
241 !jyg<
242 ! ======================================================
243 ! calcul de l'impaction
244 ! ======================================================
245
246 ! impaction sur la surface de la colonne de la descente insaturee
247 ! On prend la moyenne des precip entre le niveau i+1 et i
248 ! I=3/4* (P(1+1)+P(i))/2 / (sigd*r*rho_l)
249 ! 1000kg/m3= densite de l'eau
250 ! 0.75e-3 = 3/4 /1000
251 ! Par la suite, I est tout le temps multiplie par sig_d pour avoir l'impaction sur la surface de la maille
252 ! on le neglige ici pour simplifier le code
253
254 DO j=1,klev-1
255 DO i=1,klon
256 imp(i,j) = coefcoli_3d(i,j)*0.75e-3/rdrop *&
257 0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
258 ENDDO
259 ENDDO
260 !>jyg
261 ! =========================================
262 ! calcul des tendances liees au downdraft
263 ! =========================================
264 !cdir collapse
265 DO k=1,klev
266 DO j=1,klev
267 DO i=1,klon
268 zmd(i,j,k)=0.
269 za (i,j,k)=0.
270 END DO
271 END DO
272 END DO
273 ! calcul de la matrice d echange
274 ! matrice de distribution de la masse entrainee en k
275 ! commmentaire RomP : mp > 0
276 DO k=1,klev-1
277 DO i=1,klon
278 zmd(i,k,k)=max(0.,mp(i,k)-mp(i,k+1)) ! ~ mk(k)
279 END DO
280 END DO
281 DO k=2,klev
282 DO j=k-1,1,-1
283 DO i=1,klon
284 IF(mp(i,j+1).GT.1.e-10) THEN
285 zmd(i,j,k)=zmd(i,j+1,k)*min(1.,mp(i,j)/mp(i,j+1)) !det ~ mk(j)=mk(j+1)*mp(i,j)/mp(i,j+1)
286 ENDIF
287 END DO
288 END DO
289 END DO
290 DO k=1,klev
291 DO j=1,klev-1
292 DO i=1,klon
293 za(i,j,k)=max(0.,zmd(i,j+1,k)-zmd(i,j,k))
294 END DO
295 END DO
296 END DO
297 !!!!! quantite de traceur dans la descente d'air insaturee : 4 juin 2012
298 DO k=1,klev
299 DO j=1,klev-1
300 DO i=1,klon
301 IF(mp(i,j+1).GT.1.e-10) THEN
302 qTrdi(i,j+1,it)=qTrdi(i,j+1,it)+(zmd(i,j+1,k)/mp(i,j+1))*tr(i,k,it)
303 ELSE
304 qTrdi(i,j,it)=0.!tr(i,j,it)
305 ENDIF
306 ENDDO
307 ENDDO
308 ENDDO
309 !!!!!
310 !
311 ! rajout du terme lie a l ascendance induite
312 !
313 DO j=2,klev
314 DO i=1,klon
315 za(i,j,j-1)=za(i,j,j-1)+mp(i,j)
316 END DO
317 END DO
318 !
319 ! tendance courants insatures ! sans lessivage ancien schema
320 !
321 DO k=1,klev
322 DO j=1,klev
323 DO i=1,klon
324 zmfd(i,j,it)=zmfd(i,j,it)+za(i,j,k)*(tr(i,k,it)-tr(i,j,it))
325 END DO
326 END DO
327 END DO
328 !
329 ! =========================================
330 ! calcul des tendances liees aux courants satures j <-> z ; k <-> z'
331 ! =========================================
332 !RL
333 ! Feeding concentrations
334 DO j=1,klev
335 DO i=1,klon
336 qfeed(i,it)=qfeed(i,it)+wght_cvfd(i,j)*tr(i,j,it)
337 END DO
338 END DO
339 !RL
340 !
341 DO j=1,klev
342 DO i=1,klon
343 !RL
344 !! zmfa(i,j,it)=da(i,j)*(tr(i,1,it)-tr(i,j,it)) ! da
345 zmfa(i,j,it)=da(i,j)*(qfeed(i,it)-tr(i,j,it)) ! da
346 !RL
347 END DO
348 END DO
349 !
350 DO k=1,klev
351 DO j=1,klev
352 DO i=1,klon
353 zmfp(i,j,it)=zmfp(i,j,it)+phi(i,j,k)*(tr(i,k,it)-tr(i,j,it)) ! phi
354 END DO
355 END DO
356 END DO
357 ! RomP ajout des matrices liees au lessivage
358 DO j=1,klev
359 DO i=1,klon
360 zmfd1a(i,j,it)=d1a(i,j)*tr(i,1,it) ! da1
361 zmfdam(i,j,it)=dam(i,j)*tr(i,1,it) ! dam
362 END DO
363 END DO
364 DO k=1,klev
365 DO j=1,klev
366 DO i=1,klon
367 zmfphi2(i,j,it)=zmfphi2(i,j,it)+phi2(i,j,k)*tr(i,k,it) ! psi
368 END DO
369 END DO
370 END DO
371 DO j=1,klev-1
372 DO i=1,klon
373 zmfu(i,j,it)=max(0.,upd(i,j+1)+dnd(i,j+1))*(tr(i,j+1,it)-tr(i,j,it))
374 END DO
375 END DO
376 DO j=2,klev
377 DO i=1,klon
378 zmfu(i,j,it)=zmfu(i,j,it)+min(0.,upd(i,j)+dnd(i,j))*(tr(i,j,it)-tr(i,j-1,it))
379 END DO
380 END DO
381 ! ===================================================
382 ! calcul des tendances liees aux courants insatures
383 ! ===================================================
384 ! pression
385 DO k=1, klev
386 DO i=1, klon
387 dxpres(i,k)=paprs(i,k)-paprs(i,k+1)
388 ENDDO
389 ENDDO
390 pdtimeRG=pdtime*RG
391
392 ! q_pa et q_pm traceurs issues des courants satures se retrouvant dans les precipitations
393 DO j=1,klev
394 DO i=1,klon
395 IF(j.GE.icb(i).AND.j.LE.inb(i)) THEN
396 IF(clw(i,j).GT.1.e-16) THEN
397 !JE qPa(i,j,it)=ccntrAA_coef*tr(i,1,it)/clw(i,j)
398 qPa(i,j,it)=ccntrAA_3d(i,j)*tr(i,1,it)/clw(i,j)
399 ELSE
400 qPa(i,j,it)=0.
401 ENDIF
402 ENDIF
403 END DO
404 END DO
405
406 ! calcul de q_pm en 2 parties :
407 ! 1) calcul de sa valeur pour un niveau z' donne
408 ! 2) integration sur la verticale sur z'
409 DO j=1,klev
410 DO k=1,j-1
411 DO i=1,klon
412 IF(k.GE.icb(i).AND.k.LE.inb(i).AND.&
413 j.LE.inb(i)) THEN
414 IF(elij(i,k,j).GT.1.e-16) THEN
415 !JE qMeltmp(i,j,it)=((1-ep(i,k))*ccntrAA_coef*tr(i,1,it)&
416 !JE *(1.-sij(i,k,j)) +ccntrENV_coef&
417 !JE *tr(i,k,it)*sij(i,k,j)) / elij(i,k,j)
418 qMeltmp(i,j,it)=((1-ep(i,k))*ccntrAA_3d(i,k)*tr(i,1,it)&
419 *(1.-sij(i,k,j)) +ccntrENV_3d(i,k)&
420 *tr(i,k,it)*sij(i,k,j)) / elij(i,k,j)
421 ELSE
422 qMeltmp(i,j,it)=0.
423 ENDIF
424 qpmMint(i,j,it)=qpmMint(i,j,it) + qMeltmp(i,j,it)*epmlmMm(i,j,k)
425 Mint(i,j)=Mint(i,j) + epmlmMm(i,j,k)
426 ENDIF ! end if dans nuage
427 END DO
428 END DO
429 END DO
430
431 DO j=1,klev
432 DO i=1,klon
433 IF(Mint(i,j).GT.1.e-16) THEN
434 qMel(i,j,it)=qpmMint(i,j,it)/Mint(i,j)
435 ELSE
436 qMel(i,j,it)=0.
437 ENDIF
438 END DO
439 END DO
440
441 ! calcul de q_d et q_p traceurs de la descente precipitante
442 DO j=klev-1,1,-1
443 DO i=1,klon
444 IF(mp(i,j+1).GT.mp(i,j).AND.mp(i,j+1).GT.1.e-10) THEN ! detrainement
445 kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
446 (-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))&
447 + (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j)))
448
449 ELSEIF(mp(i,j).GT.mp(i,j+1).AND.mp(i,j).GT.1.e-10) THEN! entrainement
450 IF(j.eq.1) THEN
451 kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
452 (-mp(i,2)-imp(i,j)/RG*dxpres(i,j))&
453 + (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j)))
454 ELSE
455 kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
456 (-mp(i,j)-imp(i,j)/RG*dxpres(i,j))&
457 + (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j)))
458 ENDIF
459 ELSE
460 kappa(i,j)=1.
461 ENDIF
462 ENDDO
463 ENDDO
464
465 DO j=klev-1,1,-1
466 DO i=1,klon
467 IF (abs(kappa(i,j)).LT.1.e-25) THEN !si denominateur nul (il peut y avoir des mp!=0)
468 kappa(i,j)=1.
469 IF(j.eq.1) THEN
470 qDi(i,j,it)=qDi(i,j+1,it) !orig tr(i,j,it) ! mp(1)=0 donc tout vient de la couche supérieure
471 ELSEIF(mp(i,j+1).GT.mp(i,j).AND.mp(i,j+1).GT.1.e-10) THEN
472 qDi(i,j,it)=qDi(i,j+1,it)
473 ELSEIF(mp(i,j).GT.mp(i,j+1).AND.mp(i,j).GT.1.e-10) THEN! entrainement
474 qDi(i,j,it)=(-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it))/(-mp(i,j))
475 ELSE ! si mp (i)=0 et mp(j+1)=0
476 qDi(i,j,it)=tr(i,j,it) ! orig 0.
477 ENDIF
478
479 IF(NO_precip(i,j)) THEN
480 qPr(i,j,it)=0.
481 ELSE
482 qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
483 Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)&
484 +imp(i,j)/RG*dxpres(i,j)*qDi(i,j,it))/&
485 (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
486 ENDIF
487 ELSE ! denominateur non nul
488 kappa(i,j)=1./kappa(i,j)
489 ! calcul de qd et qp
490 !!jyg (20130119) correction pour le sommet du nuage
491 !! if(j.GE.inb(i)) THEN !au-dessus du nuage, sommet inclu
492 if(j.GT.inb(i)) THEN !au-dessus du nuage
493 qDi(i,j,it)=tr(i,j,it) ! pas de descente => environnement = descente insaturee
494 qPr(i,j,it)=0.
495
496 ! vvv premiere couche du modele ou mp(1)=0 ! det tout le temps vvv
497 ELSEIF(j.eq.1) THEN
498 if(mp(i,2).GT.1.e-10) THEN !mp(2) non nul -> detrainement (car mp(1) = 0) !ent pas possible
499 if(NO_precip(i,j)) THEN !pas de precip en (i)
500 qDi(i,j,it)=qDi(i,j+1,it)
501 qPr(i,j,it)=0.
502 ELSE
503 qDi(i,j,it)=kappa(i,j)*(&
504 (-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
505 Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +&
506 (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
507 (-mp(i,j+1)*qDi(i,j+1,it)))
508
509 qPr(i,j,it)=kappa(i,j)*(&
510 (-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))*&
511 ((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
512 Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))&
513 +(-mp(i,j+1)*qDi(i,j+1,it)) * (imp(i,j)/RG*dxpres(i,j)))
514 ENDIF
515
516 ELSE !mp(2) nul -> plus de descente insaturee -> pluie agit sur environnement
517 qDi(i,j,it)=tr(i,j,it) ! orig 0.
518 if(NO_precip(i,j)) THEN
519 qPr(i,j,it)=0.
520 ELSE
521 qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
522 Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)&
523 +imp(i,j)/RG*dxpres(i,j)*tr(i,j,it))/&
524 (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
525 ENDIF
526
527 ENDIF !mp(2) nul ou non
528
529 ! vvv (j!=1.AND.j.LT.inb(i)) en-dessous du sommet nuage vvv
530 ELSE
531 !------------------------------------------------------------- detrainement
532 if(mp(i,j+1).GT.mp(i,j).AND.mp(i,j+1).GT.1.e-10) THEN !mp(i,j).GT.1.e-10) THEN
533 if(NO_precip(i,j)) THEN
534 qDi(i,j,it)=qDi(i,j+1,it)
535 qPr(i,j,it)=0.
536 ELSE
537 qDi(i,j,it)=kappa(i,j)*(&
538 (-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
539 Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +&
540 (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
541 (-mp(i,j+1)*qDi(i,j+1,it)))
542 !
543 qPr(i,j,it)=kappa(i,j)*(&
544 (-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))*&
545 ((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
546 Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))&
547 +(-mp(i,j+1)*qDi(i,j+1,it)) * (imp(i,j)/RG*dxpres(i,j)))
548 ENDIF !precip
549 !------------------------------------------------------------- entrainement
550 ELSEIF(mp(i,j).GT.mp(i,j+1).AND.mp(i,j).GT.1.e-10) THEN
551 if(NO_precip(i,j)) THEN
552 qDi(i,j,it)=(-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it))/(-mp(i,j))
553 qPr(i,j,it)=0.
554 ELSE
555 qDi(i,j,it)=kappa(i,j)*(&
556 (-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
557 Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +&
558 (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
559 (-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it)))
560 !
561 qPr(i,j,it)=kappa(i,j)*(&
562 (-mp(i,j)-imp(i,j)/RG*dxpres(i,j))*&
563 ((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
564 Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))&
565 +(-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it))*&
566 (imp(i,j)/RG*dxpres(i,j)))
567 ENDIF !precip
568 !------------------------------------------------------------- ENDIF ! ent/det
569 ELSE !mp nul
570 qDi(i,j,it)=tr(i,j,it) ! orig 0.
571 if(NO_precip(i,j)) THEN
572 qPr(i,j,it)=0.
573 ELSE
574 qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
575 Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)&
576 +imp(i,j)/RG*dxpres(i,j)*tr(i,j,it))/&
577 (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
578 ENDIF
579 ENDIF ! mp nul ou non
580 ENDIF ! condition sur j
581 ENDIF ! kappa
582 ENDDO
583 ENDDO
584
585 !! print test descente insaturee
586 ! DO j=klev,1,-1
587 ! DO i=1,klon
588 ! if(it.eq.3) THEN
589 ! write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)') j,&
590 !! 'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),&
591 ! 'zmfdam+zmfpsi',zmfdam(i,j,it)+zmfphi2(i,j,it),'qpmMint',qpmMint(i,j,it),&
592 ! 'Pm',Pm(i,j),'Mint',Mint(i,j),&
593 !! 'zmfa',zmfa(i,j,it),'zmfp',zmfp(i,j,it),&
594 ! 'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),'zmfd1a',zmfd1a(i,j,it)
595 !! 'Pa',Pa(i,j),'eplaMm',eplaMm(i,j)
596 !! 'zmfd1a=da1*qa',zmfd1a(i,j,it),'Pa*qPa',wdtrainA(i,j)*qPa(i,j,it),'da1',d1a(i,j)
597 ! ENDIF
598 ! ENDDO
599 ! ENDDO
600
601
602 ! ===================================================
603 ! calcul final des tendances
604 ! ===================================================
605
606 DO k=klev-1,1,-1
607 DO i=1, klon
608 ! transport
609 tdcvMA=zmfd(i,k,it)+zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it) ! double comptage des downdraft insatures
610 trsptrac=zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it)
611 ! lessivage courants satures
612 !JE scavtrac=-ccntrAA_coef*zmfd1a(i,k,it)&
613 !JE -zmfphi2(i,k,it)*ccntrENV_coef&
614 !JE -zmfdam(i,k,it)*ccntrAA_coef
615 scavtrac=-ccntrAA_3d(i,k)*zmfd1a(i,k,it)&
616 -zmfphi2(i,k,it)*ccntrENV_3d(i,k)&
617 -zmfdam(i,k,it)*ccntrAA_3d(i,k)
618 ! lessivage courants insatures
619 if(k.LE.inb(i).AND.k.GT.1) THEN ! tendances dans le nuage
620 !------------------------------------------------------------- detrainement
621 if(mp(i,k+1).GT.mp(i,k).AND.mp(i,k+1).GT.1.e-10) THEN
622 uscavtrac= (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))&
623 + mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
624 !
625 ! if(it.eq.3) write(*,'(I2,1X,a,5X,e20.12,82X,a,e20.12)')k,' det incloud',&
626 ! (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+&
627 ! mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k),&
628 ! 'mp',mp(i,k)
629 !------------------------------------------------------------- entrainement
630 ELSEIF(mp(i,k).GT.mp(i,k+1).AND.mp(i,k).GT.1.e-10) THEN
631 uscavtrac= mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
632 !
633 ! if(it.eq.3) write(*,'(I2,1X,a,5X,e20.12,82X,a,e20.12)')k,' ent incloud',uscavtrac*pdtimeRG/dxpres(i,k), 'mp',mp(i,k)
634 !=!------------------------------------------------------------- end ent/det
635 ELSE ! mp(i,k+1)=0. et mp(i,k)=0. pluie directement sur l environnement
636
637 if(NO_precip(i,k)) THEN
638 uscavtrac=0.
639 ! if(it.eq.3) write(*,'(I2,1X,a,e20.12,82X,a,e20.12)')k,' no P ent incloud',uscavtrac*pdtimeRG/dxpres(i,k), 'mp',mp(i,k)
640 ELSE
641 uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/RG+evap(i,k)*qPr(i,k,it)*dxpres(i,k)/RG
642 ! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac*pdtimeRG/dxpres(i,k), 'mp',mp(i,k)
643 !!JE adds
644 ! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'imp',imp(i,k)
645 ! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'tr',tr(i,k,it)
646 ! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'evap',evap(i,k)
647 ! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'qPr',qPr(i,k,it)
648 ! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'dxpres',dxpres(i,k)
649 !!Je end
650
651 ENDIF
652 ENDIF ! mp/det/ent
653 !------------------------------------------------------------- premiere couche
654 ELSEIF(k.eq.1) THEN ! mp(1)=0.
655 if(mp(i,2).GT.1.e-10) THEN !detrainement
656 uscavtrac= (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it)) !&
657 ! + mp(i,2)*(0.-tr(i,k,it))
658 !
659 ! if(it.eq.3) write(*,'(I2,1X,a,e20.12,84X,a,e20.12)')k,' 1 det',&
660 ! (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+&
661 ! mp(i,2)*(0.-tr(i,k,it))*pdtimeRG/dxpres(i,k),&
662 ! 'mp',mp(i,k)
663 ELSE ! mp(2) = 0 = mp(1) pas de descente insaturee, rien ne se passe s'il ne pleut pas, sinon pluie->env
664 if(NO_precip(i,1)) THEN
665 uscavtrac=0.
666 ELSE
667 uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/RG+evap(i,k)*qPr(i,k,it)*dxpres(i,k)/RG
668 ENDIF
669 ! if(it.eq.3) write(*,'(I2,1X,a,2X,e20.12,82X,a,e20.12)')k,'1 P env incloud',uscavtrac*pdtimeRG/dxpres(i,k), 'mp',mp(i,k)
670 ENDIF
671
672 ELSE ! k > INB au-dessus du nuage
673 uscavtrac=0.
674 ENDIF
675
676 ! ===== tendances finales ======
677 trsptd(i,k,it)=trsptrac*pdtimeRG/dxpres(i,k) ! td transport sans eau dans courants satures
678 dtrSscav(i,k,it)=scavtrac*pdtimeRG/dxpres(i,k) ! td du lessivage dans courants satures
679 dtrUscav(i,k,it)=uscavtrac*pdtimeRG/dxpres(i,k) ! td courant insat
680 dtrsat(i,k,it)=(trsptrac+scavtrac)*pdtimeRG/dxpres(i,k) ! td courant sat
681 dtrcv(i,k,it)=(trsptrac+scavtrac+uscavtrac)*pdtimeRG/dxpres(i,k)!dtrsat(i,k,it)+dtrUscav(i,k,it) td conv
682 !!!!!!
683 dtrcvMA(i,k,it)=tdcvMA*pdtimeRG/dxpres(i,k) ! MA tendance convection
684 ENDDO
685 ENDDO
686
687 ! test de conservation du traceur
688 !print*,"_____________________________________________________________"
689 !print*," "
690 ! conserv=0.
691 ! conservMA=0.
692 ! DO k= klev-1,1,-1
693 ! DO i=1, klon
694 ! conserv=conserv+dtrcv(i,k,it)* &
695 ! (paprs(i,k)-paprs(i,k+1))/RG
696 ! conservMA=conservMA+dtrcvMA(i,k,it)* &
697 ! (paprs(i,k)-paprs(i,k+1))/RG
698 !
699 ! if(it.eq.3) write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)') k,&
700 ! 'MA td ',dtrcvMA(i,k,it)*dxpres(i,k)/RG,&
701 ! ' td',dtrcv(i,k,it)*dxpres(i,k)/RG,' conservMA ',conservMA,'conserv ',conserv
702 !!
703 ! ENDDO
704 ! ENDDO
705 ! if(it.eq.3) print *,'it',it,'conserv ',conserv,'conservMA ',conservMA
706
707 END SUBROUTINE cvltr_scav
708