GCC Code Coverage Report


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

Line Branch Exec Source
1 !
2 ! $Id $
3 !
4 SUBROUTINE cvltr(pdtime, da, phi,phi2,d1a,dam, mpIN,epIN, &
5 !! sigd,sij,clw,elij,epmlmMm,eplaMm, & !RL
6 sigd,sij,wght_cvfd,clw,elij,epmlmMm,eplaMm, & !RL
7 pmflxrIN,pmflxsIN,ev,te,wdtrainA,wdtrainM, &
8 paprs,it,tr,upd,dnd,inb,icb, &
9 dtrcv,trsptd,dtrSscav,dtrsat,dtrUscav,qDi,qPr, &
10 qPa,qMel,qTrdi,dtrcvMA,Mint, &
11 zmfd1a,zmfphi2,zmfdam)
12 USE IOIPSL
13 USE dimphy
14 USE infotrac_phy, ONLY : nbtr,tname
15 IMPLICIT NONE
16 !=====================================================================
17 ! Objet : convection des traceurs / KE
18 ! Auteurs: M-A Filiberti and J-Y Grandpeix
19 ! modifiee par R Pilon : lessivage des traceurs / KE
20 !=====================================================================
21
22 include "YOMCST.h"
23 include "YOECUMF.h"
24 include "conema3.h"
25
26 ! Entree
27 REAL,INTENT(IN) :: pdtime
28 REAL,DIMENSION(klon,klev),INTENT(IN) :: da
29 REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi
30 ! RomP
31 REAL,DIMENSION(klon,klev),INTENT(IN) :: d1a,dam ! matrices pour simplifier
32 REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi2 ! l'ecriture des tendances
33 !
34 REAL,DIMENSION(klon,klev),INTENT(IN) :: mpIN
35 REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression aux 1/2 couches (bas en haut)
36 ! REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression aux 1/2 couches (bas en haut)
37 REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr ! q de traceur (bas en haut)
38 INTEGER,INTENT(IN) :: it
39 REAL,DIMENSION(klon,klev),INTENT(IN) :: upd ! saturated updraft mass flux
40 REAL,DIMENSION(klon,klev),INTENT(IN) :: dnd ! saturated downdraft mass flux
41 !
42 REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainA ! masses precipitantes de l'asc adiab
43 REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainM ! masses precipitantes des melanges
44 REAL,DIMENSION(klon,klev),INTENT(IN) :: pmflxrIN ! vprecip: eau
45 REAL,DIMENSION(klon,klev),INTENT(IN) :: pmflxsIN ! vprecip: neige
46 REAL,DIMENSION(klon,klev),INTENT(IN) :: ev ! evaporation cv30_routine
47 REAL,DIMENSION(klon,klev),INTENT(IN) :: epIN
48 REAL,DIMENSION(klon,klev),INTENT(IN) :: te
49 REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: sij ! fraction dair de lenv
50 REAL,DIMENSION(klon,klev),INTENT(IN) :: wght_cvfd ! weights of the layers feeding convection
51 REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: elij ! contenu en eau condens�e sp�cifique/conc deau condens�e massique
52 REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: epmlmMm ! eau condensee precipitee dans mel masse dair sat
53 REAL,DIMENSION(klon,klev),INTENT(IN) :: eplaMm ! eau condensee precipitee dans aa masse dair sat
54
55 REAL,DIMENSION(klon,klev),INTENT(IN) :: clw ! contenu en eau condens�e dans lasc adiab
56 REAL,DIMENSION(klon),INTENT(IN) :: sigd
57 INTEGER,DIMENSION(klon),INTENT(IN) :: icb,inb
58 ! Sortie
59 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrcv ! tendance totale (bas en haut)
60 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrcvMA ! M-A Filiberti
61 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: trsptd ! tendance du transport
62 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrSscav ! tendance du lessivage courant sat
63 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrsat ! tendance trsp+sat scav
64 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrUscav ! tendance du lessivage courant unsat
65 !
66 ! Variables locales
67 INTEGER :: i,j,k
68 REAL,DIMENSION(klon,klev) :: dxpres ! difference de pression entre niveau (j+1) et (j)
69 REAL :: pdtimeRG ! pas de temps * gravite
70 ! variables pour les courants satures
71 REAL,DIMENSION(klon,klev,klev) :: zmd
72 REAL,DIMENSION(klon,klev,klev) :: za
73 REAL,DIMENSION(klon,klev,nbtr) :: zmfd,zmfa
74 REAL,DIMENSION(klon,klev,nbtr) :: zmfp,zmfu
75 REAL,DIMENSION(klon,nbtr) :: qfeed ! tracer concentration feeding convection
76
77 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfd1a
78 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfdam
79 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfphi2
80
81 ! RomP ! les variables sont nettoyees des valeurs aberrantes
82 REAL,DIMENSION(klon,klev) :: Pa, Pm ! pluie AA et m�langes, var temporaire
83 REAL,DIMENSION(klon,klev) :: pmflxs,pmflxr ! pmflxrIN,pmflxsIN sans valeur aberante
84 REAL,DIMENSION(klon,klev) :: mp ! flux de masse
85 REAL,DIMENSION(klon,klev) :: ep ! fraction d'eau convertie en precipitation
86 REAL,DIMENSION(klon,klev) :: evap ! evaporation : variable temporaire
87 REAL,DIMENSION(klon,klev) :: rho !environmental density
88
89 REAL,DIMENSION(klon,klev) :: kappa ! denominateur du au calcul de la matrice
90 ! pour obtenir qd et qp
91 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qTrdi ! traceurs descente air insature transport MA
92 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qDi ! traceurs descente insaturees
93 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qPr ! traceurs colonne precipitante
94 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qPa ! traceurs dans les precip issues lasc. adiab.
95 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qMel ! traceurs dans les precip issues des melanges
96 REAL,DIMENSION(klon,klev,nbtr) :: qMeltmp ! variable temporaire
97 REAL,DIMENSION(klon,klev,nbtr) :: qpmMint
98 REAL,DIMENSION(klon,klev),INTENT(OUT) :: Mint
99 ! tendances
100 REAL :: tdcvMA ! terme de transport de traceur (schema Marie Angele)
101 REAL :: trsptrac ! terme de transport de traceur par l'air
102 REAL :: scavtrac ! terme de lessivage courant sature
103 REAL :: uscavtrac ! terme de lessivage courant insature
104 ! impaction
105 !!! Correction apres discussion Romain P. / Olivier B.
106 !!! REAL,PARAMETER :: rdrop=2.5e-3 ! rayon des gouttes d'eau
107 REAL,PARAMETER :: rdrop=1.e-3 ! rayon des gouttes d'eau
108 !!!
109 REAL,DIMENSION(klon,klev) :: imp ! coefficient d'impaction
110 ! parametres lessivage
111 REAL :: ccntrAA_coef ! \alpha_a : fract aerosols de l'AA convertis en CCN
112 REAL :: ccntrENV_coef ! \beta_m : fract aerosols de l'env convertis en CCN
113 REAL :: coefcoli ! coefficient de collision des gouttes sur les aerosols
114 !
115 LOGICAL,DIMENSION(klon,klev) :: NO_precip
116 ! LOGICAL :: scavON
117 ! var tmp tests
118 REAL :: conserv
119 real :: conservMA
120
121 ! coefficient lessivage
122 ccntrAA_coef = 0.
123 ccntrENV_coef = 0.
124 coefcoli = 0.
125
126 !$OMP MASTER
127 call getin('ccntrAA_coef',ccntrAA_coef)
128 call getin('ccntrENV_coef',ccntrENV_coef)
129 call getin('coefcoli',coefcoli)
130 !$OMP END MASTER
131 !$OMP BARRIER
132 print*,'cvltr coef lessivage convectif', ccntrAA_coef,ccntrENV_coef,coefcoli
133
134 ! scavON=.TRUE.
135 ! if(scavON) then
136 ! ccntrAA_coef = 1.
137 ! ccntrENV_coef = 1.
138 ! coefcoli = 1.
139 ! else
140 ! ccntrAA_coef = 0.
141 ! ccntrENV_coef = 0.
142 ! coefcoli = 0.
143 ! endif
144
145 ! ======================================================
146 ! calcul de l'impaction
147 ! ======================================================
148 !initialisation
149 do j=1,klev
150 do i=1,klon
151 imp(i,j)=0.
152 enddo
153 enddo
154 ! impaction sur la surface de la colonne de la descente insaturee
155 ! On prend la moyenne des precip entre le niveau i+1 et i
156 ! I=3/4* (P(1+1)+P(i))/2 / (sigd*r*rho_l)
157 ! 1000kg/m3= densit� de l'eau
158 ! 0.75e-3 = 3/4 /1000
159 ! Par la suite, I est tout le temps multipli� par sig_d pour avoir l'impaction sur la surface de la maille
160 ! on le n�glige ici pour simplifier le code
161 do j=1,klev-1
162 do i=1,klon
163 imp(i,j) = coefcoli*0.75e-3/rdrop *&
164 0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
165 ! rho(i,j)=pplay(i,j)/(rd*te(i,j))
166 enddo
167 enddo
168 !
169 ! initialisation pour flux de traceurs, td et autre
170 trsptrac = 0.
171 scavtrac = 0.
172 uscavtrac = 0.
173 qfeed(:,it) = 0. !RL
174 DO j=1,klev
175 DO i=1,klon
176 zmfd(i,j,it)=0.
177 zmfa(i,j,it)=0.
178 zmfu(i,j,it)=0.
179 zmfp(i,j,it)=0.
180 zmfphi2(i,j,it)=0.
181 zmfd1a(i,j,it)=0.
182 zmfdam(i,j,it)=0.
183 qDi(i,j,it)=0.
184 qPr(i,j,it)=0.
185 qPa(i,j,it)=0.
186 qMel(i,j,it)=0.
187 qMeltmp(i,j,it)=0.
188 qTrdi(i,j,it)=0.
189 kappa(i,j)=0.
190 trsptd(i,j,it)=0.
191 dtrsat(i,j,it)=0.
192 dtrSscav(i,j,it)=0.
193 dtrUscav(i,j,it)=0.
194 dtrcv(i,j,it)=0.
195 dtrcvMA(i,j,it)=0.
196 evap(i,j)=0.
197 dxpres(i,j)=0.
198 qpmMint(i,j,it)=0.
199 Mint(i,j)=0.
200 END DO
201 END DO
202
203 ! suppression des valeurs tr�s faibles (~1e-320)
204 ! multiplication de levaporation pour lavoir par unite de temps
205 ! et par unite de surface de la maille
206 ! -> cv30_unsat : evap : masse evaporee/s/(m2 de la descente)
207 DO j=1,klev
208 DO i=1,klon
209 if(ev(i,j).lt.1.e-16) then
210 evap(i,j)=0.
211 else
212 evap(i,j)=ev(i,j)*sigd(i)
213 endif
214 END DO
215 END DO
216
217 DO j=1,klev
218 DO i=1,klon
219 if(j.lt.klev) then
220 if(epIN(i,j).lt.1.e-32) then
221 ep(i,j)=0.
222 else
223 ep(i,j)=epIN(i,j)
224 endif
225 else
226 ep(i,j)=epmax
227 endif
228 if(mpIN(i,j).lt.1.e-32) then
229 mp(i,j)=0.
230 else
231 mp(i,j)=mpIN(i,j)
232 endif
233 if(pmflxsIN(i,j).lt.1.e-32) then
234 pmflxs(i,j)=0.
235 else
236 pmflxs(i,j)=pmflxsIN(i,j)
237 endif
238 if(pmflxrIN(i,j).lt.1.e-32) then
239 pmflxr(i,j)=0.
240 else
241 pmflxr(i,j)=pmflxrIN(i,j)
242 endif
243 if(wdtrainA(i,j).lt.1.e-32) then
244 Pa(i,j)=0.
245 else
246 Pa(i,j)=wdtrainA(i,j)
247 endif
248 if(wdtrainM(i,j).lt.1.e-32) then
249 Pm(i,j)=0.
250 else
251 Pm(i,j)=wdtrainM(i,j)
252 endif
253 END DO
254 END DO
255
256 !==========================================
257 DO j = klev-1,1,-1
258 DO i = 1,klon
259 NO_precip(i,j) = (pmflxr(i,j+1)+pmflxs(i,j+1)).lt.1.e-10&
260 .and.Pa(i,j).lt.1.e-10.and.Pm(i,j).lt.1.e-10
261 END DO
262 END DO
263
264 ! =========================================
265 ! calcul des tendances liees au downdraft
266 ! =========================================
267 !cdir collapse
268 DO k=1,klev
269 DO j=1,klev
270 DO i=1,klon
271 zmd(i,j,k)=0.
272 za (i,j,k)=0.
273 END DO
274 END DO
275 END DO
276 ! calcul de la matrice d echange
277 ! matrice de distribution de la masse entrainee en k
278 ! commmentaire RomP : mp > 0
279 DO k=1,klev-1
280 DO i=1,klon
281 zmd(i,k,k)=max(0.,mp(i,k)-mp(i,k+1)) ! ~ mk(k)
282 END DO
283 END DO
284 DO k=2,klev
285 DO j=k-1,1,-1
286 DO i=1,klon
287 if(mp(i,j+1).gt.1.e-10) then
288 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)
289 ENDif
290 END DO
291 END DO
292 END DO
293 DO k=1,klev
294 DO j=1,klev-1
295 DO i=1,klon
296 za(i,j,k)=max(0.,zmd(i,j+1,k)-zmd(i,j,k))
297 END DO
298 END DO
299 END DO
300 !!!!! quantite de traceur dans la descente d'air insaturee : 4 juin 2012
301 DO k=1,klev
302 DO j=1,klev-1
303 DO i=1,klon
304 if(mp(i,j+1).gt.1.e-10) then
305 qTrdi(i,j+1,it)=qTrdi(i,j+1,it)+(zmd(i,j+1,k)/mp(i,j+1))*tr(i,k,it)
306 else
307 qTrdi(i,j,it)=0.!tr(i,j,it)
308 endif
309 ENDDO
310 ENDDO
311 ENDDO
312 !!!!!
313 !
314 ! rajout du terme lie a l ascendance induite
315 !
316 DO j=2,klev
317 DO i=1,klon
318 za(i,j,j-1)=za(i,j,j-1)+mp(i,j)
319 END DO
320 END DO
321 !
322 ! tendance courants insatures ! sans lessivage ancien schema
323 !
324 DO k=1,klev
325 DO j=1,klev
326 DO i=1,klon
327 zmfd(i,j,it)=zmfd(i,j,it)+za(i,j,k)*(tr(i,k,it)-tr(i,j,it))
328 END DO
329 END DO
330 END DO
331 !
332 ! =========================================
333 ! calcul des tendances liees aux courants satures j <-> z ; k <-> z'
334 ! =========================================
335 !
336 !RL
337 ! Feeding concentrations
338 DO j=1,klev
339 DO i=1,klon
340 qfeed(i,it)=qfeed(i,it)+wght_cvfd(i,j)*tr(i,j,it)
341 END DO
342 END DO
343 !RL
344 !
345 DO j=1,klev
346 DO i=1,klon
347 !RL
348 !! zmfa(i,j,it)=da(i,j)*(tr(i,1,it)-tr(i,j,it)) ! da
349 zmfa(i,j,it)=da(i,j)*(qfeed(i,it)-tr(i,j,it)) ! da
350 !RL
351 END DO
352 END DO
353 !
354 DO k=1,klev
355 DO j=1,klev
356 DO i=1,klon
357 zmfp(i,j,it)=zmfp(i,j,it)+phi(i,j,k)*(tr(i,k,it)-tr(i,j,it)) ! phi
358 END DO
359 END DO
360 END DO
361 ! RomP ajout des matrices liees au lessivage
362 DO j=1,klev
363 DO i=1,klon
364 zmfd1a(i,j,it)=d1a(i,j)*tr(i,1,it) ! da1
365 zmfdam(i,j,it)=dam(i,j)*tr(i,1,it) ! dam
366 END DO
367 END DO
368 DO k=1,klev
369 DO j=1,klev
370 DO i=1,klon
371 zmfphi2(i,j,it)=zmfphi2(i,j,it)+phi2(i,j,k)*tr(i,k,it) ! psi
372 END DO
373 END DO
374 END DO
375 DO j=1,klev-1
376 DO i=1,klon
377 zmfu(i,j,it)=max(0.,upd(i,j+1)+dnd(i,j+1))*(tr(i,j+1,it)-tr(i,j,it))
378 END DO
379 END DO
380 DO j=2,klev
381 DO i=1,klon
382 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))
383 END DO
384 END DO
385 ! ===================================================
386 ! calcul des tendances liees aux courants insatures
387 ! ===================================================
388 ! pression
389 DO k=1, klev
390 DO i=1, klon
391 dxpres(i,k)=paprs(i,k)-paprs(i,k+1)
392 ENDDO
393 ENDDO
394 pdtimeRG=pdtime*RG
395
396 ! q_pa et q_pm traceurs issues des courants satures se retrouvant dans les precipitations
397 DO j=1,klev
398 DO i=1,klon
399 if(j.ge.icb(i).and.j.le.inb(i)) then
400 if(clw(i,j).gt.1.e-16) then
401 qPa(i,j,it)=ccntrAA_coef*tr(i,1,it)/clw(i,j)
402 else
403 qPa(i,j,it)=0.
404 endif
405 endif
406 END DO
407 END DO
408
409 ! calcul de q_pm en 2 parties :
410 ! 1) calcul de sa valeur pour un niveau z' donne
411 ! 2) integration sur la verticale sur z'
412 DO j=1,klev
413 DO k=1,j-1
414 DO i=1,klon
415 if(k.ge.icb(i).and.k.le.inb(i).and.&
416 j.le.inb(i)) then
417 if(elij(i,k,j).gt.1.e-16) then
418 qMeltmp(i,j,it)=((1-ep(i,k))*ccntrAA_coef*tr(i,1,it)&
419 *(1.-sij(i,k,j)) +ccntrENV_coef&
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 scavtrac=-ccntrAA_coef*zmfd1a(i,k,it)&
613 -zmfphi2(i,k,it)*ccntrENV_coef&
614 -zmfdam(i,k,it)*ccntrAA_coef
615 ! lessivage courants insatures
616 if(k.le.inb(i).and.k.gt.1) then ! tendances dans le nuage
617 !------------------------------------------------------------- detrainement
618 if(mp(i,k+1).gt.mp(i,k).and.mp(i,k+1).gt.1.e-10) then
619 uscavtrac= (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))&
620 + mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
621 !
622 ! if(it.eq.3) write(*,'(I2,1X,a,5X,e20.12,82X,a,e20.12)')k,' det incloud',&
623 ! (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+&
624 ! mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k),&
625 ! 'mp',mp(i,k)
626 !------------------------------------------------------------- entrainement
627 elseif(mp(i,k).gt.mp(i,k+1).and.mp(i,k).gt.1.e-10) then
628 uscavtrac= mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
629 !
630 ! 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)
631 !=!------------------------------------------------------------- end ent/det
632 else ! mp(i,k+1)=0. et mp(i,k)=0. pluie directement sur l environnement
633
634 if(NO_precip(i,k)) then
635 uscavtrac=0.
636 ! 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)
637 else
638 uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/RG+evap(i,k)*qPr(i,k,it)*dxpres(i,k)/RG
639 ! 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)
640 endif
641 endif ! mp/det/ent
642 !------------------------------------------------------------- premiere couche
643 elseif(k.eq.1) then ! mp(1)=0.
644 if(mp(i,2).gt.1.e-10) then !detrainement
645 uscavtrac= (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it)) !&
646 ! + mp(i,2)*(0.-tr(i,k,it))
647 !
648 ! if(it.eq.3) write(*,'(I2,1X,a,e20.12,84X,a,e20.12)')k,' 1 det',&
649 ! (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+&
650 ! mp(i,2)*(0.-tr(i,k,it))*pdtimeRG/dxpres(i,k),&
651 ! 'mp',mp(i,k)
652 else ! mp(2) = 0 = mp(1) pas de descente insaturee, rien ne se passe s'il ne pleut pas, sinon pluie->env
653 if(NO_precip(i,1)) then
654 uscavtrac=0.
655 else
656 uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/RG+evap(i,k)*qPr(i,k,it)*dxpres(i,k)/RG
657 endif
658 ! 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)
659 endif
660
661 else ! k > INB au-dessus du nuage
662 uscavtrac=0.
663 endif
664
665 ! ===== tendances finales ======
666 trsptd(i,k,it)=trsptrac*pdtimeRG/dxpres(i,k) ! td transport sans eau dans courants satures
667 dtrSscav(i,k,it)=scavtrac*pdtimeRG/dxpres(i,k) ! td du lessivage dans courants satures
668 dtrUscav(i,k,it)=uscavtrac*pdtimeRG/dxpres(i,k) ! td courant insat
669 dtrsat(i,k,it)=(trsptrac+scavtrac)*pdtimeRG/dxpres(i,k) ! td courant sat
670 dtrcv(i,k,it)=(trsptrac+scavtrac+uscavtrac)*pdtimeRG/dxpres(i,k)!dtrsat(i,k,it)+dtrUscav(i,k,it) td conv
671 !!!!!!
672 dtrcvMA(i,k,it)=tdcvMA*pdtimeRG/dxpres(i,k) ! MA tendance convection
673 ENDDO
674 ENDDO
675
676 ! test de conservation du traceur
677 !print*,"_____________________________________________________________"
678 !print*," "
679 ! conserv=0.
680 ! conservMA=0.
681 ! DO k= klev-1,1,-1
682 ! DO i=1, klon
683 ! conserv=conserv+dtrcv(i,k,it)* &
684 ! (paprs(i,k)-paprs(i,k+1))/RG
685 ! conservMA=conservMA+dtrcvMA(i,k,it)* &
686 ! (paprs(i,k)-paprs(i,k+1))/RG
687 !
688 ! if(it.eq.3) write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)') k,&
689 ! 'MA td ',dtrcvMA(i,k,it)*dxpres(i,k)/RG,&
690 ! ' td',dtrcv(i,k,it)*dxpres(i,k)/RG,' conservMA ',conservMA,'conserv ',conserv
691 !!
692 ! ENDDO
693 ! ENDDO
694 ! if(it.eq.3) print *,'it',it,'conserv ',conserv,'conservMA ',conservMA
695
696 END SUBROUTINE cvltr
697