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