My Project
 All Classes Files Functions Variables Macros
cvltr.F90
Go to the documentation of this file.
1 !
2 ! $Id $
3 !
4 SUBROUTINE cvltr(pdtime, da, phi,phi2,d1a,dam, mpIN,epIN, &
5  sigd,sij,clw,elij,epmlmmm,eplamm, &
6  pmflxrin,pmflxsin,ev,te,wdtraina,wdtrainm, &
7  paprs,it,tr,upd,dnd,inb,icb, &
8  dtrcv,trsptd,dtrsscav,dtrsat,dtruscav,qdi,qpr, &
9  qpa,qmel,qtrdi,dtrcvma,mint, &
10  zmfd1a,zmfphi2,zmfdam)
11  USE ioipsl
12  USE dimphy
13  USE infotrac, ONLY : nbtr,tname
14  IMPLICIT NONE
15 !=====================================================================
16 ! Objet : convection des traceurs / KE
17 ! Auteurs: M-A Filiberti and J-Y Grandpeix
18 ! modifiee par R Pilon : lessivage des traceurs / KE
19 !=====================================================================
20 
21  include "YOMCST.h"
22  include "YOECUMF.h"
23  include "conema3.h"
24 
25 ! Entree
26  REAL,INTENT(IN) :: pdtime
27  REAL,DIMENSION(klon,klev),INTENT(IN) :: da
28  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi
29 ! RomP
30  REAL,DIMENSION(klon,klev),INTENT(IN) :: d1a,dam ! matrices pour simplifier
31  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi2 ! l'ecriture des tendances
32 !
33  REAL,DIMENSION(klon,klev),INTENT(IN) :: mpin
34  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression aux 1/2 couches (bas en haut)
35 ! REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression aux 1/2 couches (bas en haut)
36  REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr ! q de traceur (bas en haut)
37  INTEGER,INTENT(IN) :: it
38  REAL,DIMENSION(klon,klev),INTENT(IN) :: upd ! saturated updraft mass flux
39  REAL,DIMENSION(klon,klev),INTENT(IN) :: dnd ! saturated downdraft mass flux
40 !
41  REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtraina ! masses precipitantes de l'asc adiab
42  REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainm ! masses precipitantes des melanges
43  REAL,DIMENSION(klon,klev),INTENT(IN) :: pmflxrin ! vprecip: eau
44  REAL,DIMENSION(klon,klev),INTENT(IN) :: pmflxsin ! vprecip: neige
45  REAL,DIMENSION(klon,klev),INTENT(IN) :: ev ! evaporation cv30_routine
46  REAL,DIMENSION(klon,klev),INTENT(IN) :: epin
47  REAL,DIMENSION(klon,klev),INTENT(IN) :: te
48  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: sij ! fraction dair de lenv
49  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: elij ! contenu en eau condensée spécifique/conc deau condensée massique
50  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: epmlmmm ! eau condensee precipitee dans mel masse dair sat
51  REAL,DIMENSION(klon,klev),INTENT(IN) :: eplamm ! eau condensee precipitee dans aa masse dair sat
52 
53  REAL,DIMENSION(klon,klev),INTENT(IN) :: clw ! contenu en eau condensée dans lasc adiab
54  REAL,DIMENSION(klon),INTENT(IN) :: sigd
55  INTEGER,DIMENSION(klon),INTENT(IN) :: icb,inb
56 ! Sortie
57  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrcv ! tendance totale (bas en haut)
58  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrcvma ! M-A Filiberti
59  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: trsptd ! tendance du transport
60  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrsscav ! tendance du lessivage courant sat
61  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrsat ! tendance trsp+sat scav
62  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtruscav ! tendance du lessivage courant unsat
63 !
64 ! Variables locales
65  INTEGER :: i,j,k
66  REAL,DIMENSION(klon,klev) :: dxpres ! difference de pression entre niveau (j+1) et (j)
67  REAL :: pdtimerg ! pas de temps * gravite
68 ! variables pour les courants satures
69  REAL,DIMENSION(klon,klev,klev) :: zmd
70  REAL,DIMENSION(klon,klev,klev) :: za
71  REAL,DIMENSION(klon,klev,nbtr) :: zmfd,zmfa
72  REAL,DIMENSION(klon,klev,nbtr) :: zmfp,zmfu
73 
74  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfd1a
75  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfdam
76  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfphi2
77 
78 ! RomP ! les variables sont nettoyees des valeurs aberrantes
79  REAL,DIMENSION(klon,klev) :: pa, pm ! pluie AA et mélanges, var temporaire
80  REAL,DIMENSION(klon,klev) :: pmflxs,pmflxr ! pmflxrIN,pmflxsIN sans valeur aberante
81  REAL,DIMENSION(klon,klev) :: mp ! flux de masse
82  REAL,DIMENSION(klon,klev) :: ep ! fraction d'eau convertie en precipitation
83  REAL,DIMENSION(klon,klev) :: evap ! evaporation : variable temporaire
84  REAL,DIMENSION(klon,klev) :: rho !environmental density
85 
86  REAL,DIMENSION(klon,klev) :: kappa ! denominateur du au calcul de la matrice
87  ! pour obtenir qd et qp
88  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qtrdi ! traceurs descente air insature transport MA
89  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qdi ! traceurs descente insaturees
90  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qpr ! traceurs colonne precipitante
91  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qpa ! traceurs dans les precip issues lasc. adiab.
92  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qmel ! traceurs dans les precip issues des melanges
93  REAL,DIMENSION(klon,klev,nbtr) :: qmeltmp ! variable temporaire
94  REAL,DIMENSION(klon,klev,nbtr) :: qpmmint
95  REAL,DIMENSION(klon,klev),INTENT(OUT) :: mint
96 ! tendances
97  REAL :: tdcvma ! terme de transport de traceur (schema Marie Angele)
98  REAL :: trsptrac ! terme de transport de traceur par l'air
99  REAL :: scavtrac ! terme de lessivage courant sature
100  REAL :: uscavtrac ! terme de lessivage courant insature
101 ! impaction
102 !!! Correction apres discussion Romain P. / Olivier B.
103 !!! REAL,PARAMETER :: rdrop=2.5e-3 ! rayon des gouttes d'eau
104  REAL,PARAMETER :: rdrop=1.e-3 ! rayon des gouttes d'eau
105 !!!
106  REAL,DIMENSION(klon,klev) :: imp ! coefficient d'impaction
107 ! parametres lessivage
108  REAL :: ccntraa_coef ! \alpha_a : fract aerosols de l'AA convertis en CCN
109  REAL :: ccntrenv_coef ! \beta_m : fract aerosols de l'env convertis en CCN
110  REAL :: coefcoli ! coefficient de collision des gouttes sur les aerosols
111 !
112  LOGICAL,DIMENSION(klon,klev) :: no_precip
113 ! LOGICAL :: scavON
114 ! var tmp tests
115  REAL :: conserv
116  real :: conservma
117 
118 ! coefficient lessivage
119  ccntraa_coef = 0.
120  ccntrenv_coef = 0.
121  coefcoli = 0.
122 
123 !$OMP MASTER
124  call getin('ccntrAA_coef',ccntraa_coef)
125  call getin('ccntrENV_coef',ccntrenv_coef)
126  call getin('coefcoli',coefcoli)
127 !$OMP END MASTER
128 !$OMP BARRIER
129  print*,'cvltr coef lessivage convectif', ccntraa_coef,ccntrenv_coef,coefcoli
130 
131 ! scavON=.TRUE.
132 ! if(scavON) then
133 ! ccntrAA_coef = 1.
134 ! ccntrENV_coef = 1.
135 ! coefcoli = 1.
136 ! else
137 ! ccntrAA_coef = 0.
138 ! ccntrENV_coef = 0.
139 ! coefcoli = 0.
140 ! endif
141 
142 ! ======================================================
143 ! calcul de l'impaction
144 ! ======================================================
145 !initialisation
146  do j=1,klev
147  do i=1,klon
148  imp(i,j)=0.
149  enddo
150  enddo
151 ! impaction sur la surface de la colonne de la descente insaturee
152 ! On prend la moyenne des precip entre le niveau i+1 et i
153 ! I=3/4* (P(1+1)+P(i))/2 / (sigd*r*rho_l)
154 ! 1000kg/m3= densité de l'eau
155 ! 0.75e-3 = 3/4 /1000
156 ! Par la suite, I est tout le temps multiplié par sig_d pour avoir l'impaction sur la surface de la maille
157 ! on le néglige ici pour simplifier le code
158  do j=1,klev-1
159  do i=1,klon
160  imp(i,j) = coefcoli*0.75e-3/rdrop *&
161  0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
162 ! rho(i,j)=pplay(i,j)/(rd*te(i,j))
163  enddo
164  enddo
165 !
166 ! initialisation pour flux de traceurs, td et autre
167  trsptrac = 0.
168  scavtrac = 0.
169  uscavtrac = 0.
170 
171  DO j=1,klev
172  DO i=1,klon
173  zmfd(i,j,it)=0.
174  zmfa(i,j,it)=0.
175  zmfu(i,j,it)=0.
176  zmfp(i,j,it)=0.
177  zmfphi2(i,j,it)=0.
178  zmfd1a(i,j,it)=0.
179  zmfdam(i,j,it)=0.
180  qdi(i,j,it)=0.
181  qpr(i,j,it)=0.
182  qpa(i,j,it)=0.
183  qmel(i,j,it)=0.
184  qmeltmp(i,j,it)=0.
185  qtrdi(i,j,it)=0.
186  kappa(i,j)=0.
187  trsptd(i,j,it)=0.
188  dtrsat(i,j,it)=0.
189  dtrsscav(i,j,it)=0.
190  dtruscav(i,j,it)=0.
191  dtrcv(i,j,it)=0.
192  dtrcvma(i,j,it)=0.
193  evap(i,j)=0.
194  dxpres(i,j)=0.
195  qpmmint(i,j,it)=0.
196  mint(i,j)=0.
197  END DO
198  END DO
199 
200 ! suppression des valeurs très faibles (~1e-320)
201 ! multiplication de levaporation pour lavoir par unite de temps
202 ! et par unite de surface de la maille
203 ! -> cv30_unsat : evap : masse evaporee/s/(m2 de la descente)
204  DO j=1,klev
205  DO i=1,klon
206  if(ev(i,j).lt.1.e-16) then
207  evap(i,j)=0.
208  else
209  evap(i,j)=ev(i,j)*sigd(i)
210  endif
211  END DO
212  END DO
213 
214  DO j=1,klev
215  DO i=1,klon
216  if(j.lt.klev) then
217  if(epin(i,j).lt.1.e-32) then
218  ep(i,j)=0.
219  else
220  ep(i,j)=epin(i,j)
221  endif
222  else
223  ep(i,j)=epmax
224  endif
225  if(mpin(i,j).lt.1.e-32) then
226  mp(i,j)=0.
227  else
228  mp(i,j)=mpin(i,j)
229  endif
230  if(pmflxsin(i,j).lt.1.e-32) then
231  pmflxs(i,j)=0.
232  else
233  pmflxs(i,j)=pmflxsin(i,j)
234  endif
235  if(pmflxrin(i,j).lt.1.e-32) then
236  pmflxr(i,j)=0.
237  else
238  pmflxr(i,j)=pmflxrin(i,j)
239  endif
240  if(wdtraina(i,j).lt.1.e-32) then
241  pa(i,j)=0.
242  else
243  pa(i,j)=wdtraina(i,j)
244  endif
245  if(wdtrainm(i,j).lt.1.e-32) then
246  pm(i,j)=0.
247  else
248  pm(i,j)=wdtrainm(i,j)
249  endif
250  END DO
251  END DO
252 
253 !==========================================
254  DO j = klev-1,1,-1
255  DO i = 1,klon
256  no_precip(i,j) = (pmflxr(i,j+1)+pmflxs(i,j+1)).lt.1.e-10&
257  .and.pa(i,j).lt.1.e-10.and.pm(i,j).lt.1.e-10
258  END DO
259  END DO
260 
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  DO j=1,klev
333  DO i=1,klon
334  zmfa(i,j,it)=da(i,j)*(tr(i,1,it)-tr(i,j,it)) ! da
335  END DO
336  END DO
337  DO k=1,klev
338  DO j=1,klev
339  DO i=1,klon
340  zmfp(i,j,it)=zmfp(i,j,it)+phi(i,j,k)*(tr(i,k,it)-tr(i,j,it)) ! phi
341  END DO
342  END DO
343  END DO
344 ! RomP ajout des matrices liees au lessivage
345  DO j=1,klev
346  DO i=1,klon
347  zmfd1a(i,j,it)=d1a(i,j)*tr(i,1,it) ! da1
348  zmfdam(i,j,it)=dam(i,j)*tr(i,1,it) ! dam
349  END DO
350  END DO
351  DO k=1,klev
352  DO j=1,klev
353  DO i=1,klon
354  zmfphi2(i,j,it)=zmfphi2(i,j,it)+phi2(i,j,k)*tr(i,k,it) ! psi
355  END DO
356  END DO
357  END DO
358  DO j=1,klev-1
359  DO i=1,klon
360  zmfu(i,j,it)=max(0.,upd(i,j+1)+dnd(i,j+1))*(tr(i,j+1,it)-tr(i,j,it))
361  END DO
362  END DO
363  DO j=2,klev
364  DO i=1,klon
365  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))
366  END DO
367  END DO
368 ! ===================================================
369 ! calcul des tendances liees aux courants insatures
370 ! ===================================================
371 ! pression
372  DO k=1, klev
373  DO i=1, klon
374  dxpres(i,k)=paprs(i,k)-paprs(i,k+1)
375  ENDDO
376  ENDDO
377  pdtimerg=pdtime*rg
378 
379 ! q_pa et q_pm traceurs issues des courants satures se retrouvant dans les precipitations
380  DO j=1,klev
381  DO i=1,klon
382  if(j.ge.icb(i).and.j.le.inb(i)) then
383  if(clw(i,j).gt.1.e-16) then
384  qpa(i,j,it)=ccntraa_coef*tr(i,1,it)/clw(i,j)
385  else
386  qpa(i,j,it)=0.
387  endif
388  endif
389  END DO
390  END DO
391 
392 ! calcul de q_pm en 2 parties :
393 ! 1) calcul de sa valeur pour un niveau z' donne
394 ! 2) integration sur la verticale sur z'
395  DO j=1,klev
396  DO k=1,j-1
397  DO i=1,klon
398  if(k.ge.icb(i).and.k.le.inb(i).and.&
399  j.le.inb(i)) then
400  if(elij(i,k,j).gt.1.e-16) then
401  qmeltmp(i,j,it)=((1-ep(i,k))*ccntraa_coef*tr(i,1,it)&
402  *(1.-sij(i,k,j)) +ccntrenv_coef&
403  *tr(i,k,it)*sij(i,k,j)) / elij(i,k,j)
404  else
405  qmeltmp(i,j,it)=0.
406  endif
407  qpmmint(i,j,it)=qpmmint(i,j,it) + qmeltmp(i,j,it)*epmlmmm(i,j,k)
408  mint(i,j)=mint(i,j) + epmlmmm(i,j,k)
409  endif ! end if dans nuage
410  END DO
411  END DO
412  END DO
413 
414  DO j=1,klev
415  DO i=1,klon
416  if(mint(i,j).gt.1.e-16) then
417  qmel(i,j,it)=qpmmint(i,j,it)/mint(i,j)
418  else
419  qmel(i,j,it)=0.
420  endif
421  END DO
422  END DO
423 
424 ! calcul de q_d et q_p traceurs de la descente precipitante
425  DO j=klev-1,1,-1
426  DO i=1,klon
427  if(mp(i,j+1).gt.mp(i,j).and.mp(i,j+1).gt.1.e-10) then ! detrainement
428  kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+pa(i,j)+pm(i,j))*&
429  (-mp(i,j+1)-imp(i,j)/rg*dxpres(i,j))&
430  + (imp(i,j)/rg*dxpres(i,j))*(evap(i,j)/rg*dxpres(i,j)))
431 
432  elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then! entrainement
433  if(j.eq.1) then
434  kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+pa(i,j)+pm(i,j))*&
435  (-mp(i,2)-imp(i,j)/rg*dxpres(i,j))&
436  + (imp(i,j)/rg*dxpres(i,j))*(evap(i,j)/rg*dxpres(i,j)))
437  else
438  kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+pa(i,j)+pm(i,j))*&
439  (-mp(i,j)-imp(i,j)/rg*dxpres(i,j))&
440  + (imp(i,j)/rg*dxpres(i,j))*(evap(i,j)/rg*dxpres(i,j)))
441  endif
442  else
443  kappa(i,j)=1.
444  endif
445  ENDDO
446  ENDDO
447 
448  DO j=klev-1,1,-1
449  DO i=1,klon
450  if (abs(kappa(i,j)).lt.1.e-25) then !si denominateur nul (il peut y avoir des mp!=0)
451  kappa(i,j)=1.
452  if(j.eq.1) then
453  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
454  elseif(mp(i,j+1).gt.mp(i,j).and.mp(i,j+1).gt.1.e-10) then
455  qdi(i,j,it)=qdi(i,j+1,it)
456  elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then! entrainement
457  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))
458  else ! si mp (i)=0 et mp(j+1)=0
459  qdi(i,j,it)=tr(i,j,it) ! orig 0.
460  endif
461 
462  if(no_precip(i,j)) then
463  qpr(i,j,it)=0.
464  else
465  qpr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qpr(i,j+1,it)+&
466  pa(i,j)*qpa(i,j,it)+pm(i,j)*qmel(i,j,it)&
467  +imp(i,j)/rg*dxpres(i,j)*qdi(i,j,it))/&
468  (pmflxr(i,j+1)+pmflxs(i,j+1)+pa(i,j)+pm(i,j))
469  endif
470  else ! denominateur non nul
471  kappa(i,j)=1./kappa(i,j)
472 ! calcul de qd et qp
473 !!jyg (20130119) correction pour le sommet du nuage
474 !! if(j.ge.inb(i)) then !au-dessus du nuage, sommet inclu
475  if(j.gt.inb(i)) then !au-dessus du nuage
476  qdi(i,j,it)=tr(i,j,it) ! pas de descente => environnement = descente insaturee
477  qpr(i,j,it)=0.
478 
479 ! vvv premiere couche du modele ou mp(1)=0 ! det tout le temps vvv
480  elseif(j.eq.1) then
481  if(mp(i,2).gt.1.e-10) then !mp(2) non nul -> detrainement (car mp(1) = 0) !ent pas possible
482  if(no_precip(i,j)) then !pas de precip en (i)
483  qdi(i,j,it)=qdi(i,j+1,it)
484  qpr(i,j,it)=0.
485  else
486  qdi(i,j,it)=kappa(i,j)*(&
487  (-evap(i,j)/rg*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qpr(i,j+1,it)+&
488  pa(i,j)*qpa(i,j,it)+pm(i,j)*qmel(i,j,it)) +&
489  (pmflxr(i,j+1)+pmflxs(i,j+1)+pa(i,j)+pm(i,j))*&
490  (-mp(i,j+1)*qdi(i,j+1,it)))
491 
492  qpr(i,j,it)=kappa(i,j)*(&
493  (-mp(i,j+1)-imp(i,j)/rg*dxpres(i,j))*&
494  ((pmflxr(i,j+1)+pmflxs(i,j+1))*qpr(i,j+1,it)+&
495  pa(i,j)*qpa(i,j,it)+pm(i,j)*qmel(i,j,it))&
496  +(-mp(i,j+1)*qdi(i,j+1,it)) * (imp(i,j)/rg*dxpres(i,j)))
497  endif
498 
499  else !mp(2) nul -> plus de descente insaturee -> pluie agit sur environnement
500  qdi(i,j,it)=tr(i,j,it) ! orig 0.
501  if(no_precip(i,j)) then
502  qpr(i,j,it)=0.
503  else
504  qpr(i,j,it)=((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  +imp(i,j)/rg*dxpres(i,j)*tr(i,j,it))/&
507  (pmflxr(i,j+1)+pmflxs(i,j+1)+pa(i,j)+pm(i,j))
508  endif
509 
510  endif !mp(2) nul ou non
511 
512 ! vvv (j!=1.and.j.lt.inb(i)) en-dessous du sommet nuage vvv
513  else
514 !------------------------------------------------------------- detrainement
515  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
516  if(no_precip(i,j)) then
517  qdi(i,j,it)=qdi(i,j+1,it)
518  qpr(i,j,it)=0.
519  else
520  qdi(i,j,it)=kappa(i,j)*(&
521  (-evap(i,j)/rg*dxpres(i,j))*((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  (pmflxr(i,j+1)+pmflxs(i,j+1)+pa(i,j)+pm(i,j))*&
524  (-mp(i,j+1)*qdi(i,j+1,it)))
525 !
526  qpr(i,j,it)=kappa(i,j)*(&
527  (-mp(i,j+1)-imp(i,j)/rg*dxpres(i,j))*&
528  ((pmflxr(i,j+1)+pmflxs(i,j+1))*qpr(i,j+1,it)+&
529  pa(i,j)*qpa(i,j,it)+pm(i,j)*qmel(i,j,it))&
530  +(-mp(i,j+1)*qdi(i,j+1,it)) * (imp(i,j)/rg*dxpres(i,j)))
531  endif !precip
532 !------------------------------------------------------------- entrainement
533  elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then
534  if(no_precip(i,j)) then
535  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))
536  qpr(i,j,it)=0.
537  else
538  qdi(i,j,it)=kappa(i,j)*(&
539  (-evap(i,j)/rg*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qpr(i,j+1,it)+&
540  pa(i,j)*qpa(i,j,it)+pm(i,j)*qmel(i,j,it)) +&
541  (pmflxr(i,j+1)+pmflxs(i,j+1)+pa(i,j)+pm(i,j))*&
542  (-mp(i,j+1)*(qdi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it)))
543 !
544  qpr(i,j,it)=kappa(i,j)*(&
545  (-mp(i,j)-imp(i,j)/rg*dxpres(i,j))*&
546  ((pmflxr(i,j+1)+pmflxs(i,j+1))*qpr(i,j+1,it)+&
547  pa(i,j)*qpa(i,j,it)+pm(i,j)*qmel(i,j,it))&
548  +(-mp(i,j+1)*(qdi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it))*&
549  (imp(i,j)/rg*dxpres(i,j)))
550  endif !precip
551 !------------------------------------------------------------- endif ! ent/det
552  else !mp nul
553  qdi(i,j,it)=tr(i,j,it) ! orig 0.
554  if(no_precip(i,j)) then
555  qpr(i,j,it)=0.
556  else
557  qpr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qpr(i,j+1,it)+&
558  pa(i,j)*qpa(i,j,it)+pm(i,j)*qmel(i,j,it)&
559  +imp(i,j)/rg*dxpres(i,j)*tr(i,j,it))/&
560  (pmflxr(i,j+1)+pmflxs(i,j+1)+pa(i,j)+pm(i,j))
561  endif
562  endif ! mp nul ou non
563  endif ! condition sur j
564  endif ! kappa
565  ENDDO
566  ENDDO
567 
568 !! print test descente insaturee
569 ! DO j=klev,1,-1
570 ! DO i=1,klon
571 ! if(it.eq.3) then
572 ! 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,&
573 !! 'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),&
574 ! 'zmfdam+zmfpsi',zmfdam(i,j,it)+zmfphi2(i,j,it),'qpmMint',qpmMint(i,j,it),&
575 ! 'Pm',Pm(i,j),'Mint',Mint(i,j),&
576 !! 'zmfa',zmfa(i,j,it),'zmfp',zmfp(i,j,it),&
577 ! 'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),'zmfd1a',zmfd1a(i,j,it)
578 !! 'Pa',Pa(i,j),'eplaMm',eplaMm(i,j)
579 !! 'zmfd1a=da1*qa',zmfd1a(i,j,it),'Pa*qPa',wdtrainA(i,j)*qPa(i,j,it),'da1',d1a(i,j)
580 ! endif
581 ! ENDDO
582 ! ENDDO
583 
584 
585 ! ===================================================
586 ! calcul final des tendances
587 ! ===================================================
588 
589  DO k=klev-1,1,-1
590  DO i=1, klon
591 ! transport
592  tdcvma=zmfd(i,k,it)+zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it) ! double comptage des downdraft insatures
593  trsptrac=zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it)
594 ! lessivage courants satures
595  scavtrac=-ccntraa_coef*zmfd1a(i,k,it)&
596  -zmfphi2(i,k,it)*ccntrenv_coef&
597  -zmfdam(i,k,it)*ccntraa_coef
598 ! lessivage courants insatures
599  if(k.le.inb(i).and.k.gt.1) then ! tendances dans le nuage
600 !------------------------------------------------------------- detrainement
601  if(mp(i,k+1).gt.mp(i,k).and.mp(i,k+1).gt.1.e-10) then
602  uscavtrac= (-mp(i,k)+mp(i,k+1))*(qdi(i,k,it)-tr(i,k,it))&
603  + mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
604 !
605 ! if(it.eq.3) write(*,'(I2,1X,a,5X,e20.12,82X,a,e20.12)')k,' det incloud',&
606 ! (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+&
607 ! mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k),&
608 ! 'mp',mp(i,k)
609 !------------------------------------------------------------- entrainement
610  elseif(mp(i,k).gt.mp(i,k+1).and.mp(i,k).gt.1.e-10) then
611  uscavtrac= mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
612 !
613 ! 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)
614 !=!------------------------------------------------------------- end ent/det
615  else ! mp(i,k+1)=0. et mp(i,k)=0. pluie directement sur l environnement
616 
617  if(no_precip(i,k)) then
618  uscavtrac=0.
619 ! 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)
620  else
621  uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/rg+evap(i,k)*qpr(i,k,it)*dxpres(i,k)/rg
622 ! 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)
623  endif
624  endif ! mp/det/ent
625 !------------------------------------------------------------- premiere couche
626  elseif(k.eq.1) then ! mp(1)=0.
627  if(mp(i,2).gt.1.e-10) then !detrainement
628  uscavtrac= (-0.+mp(i,2))*(qdi(i,k,it)-tr(i,k,it)) !&
629 ! + mp(i,2)*(0.-tr(i,k,it))
630 !
631 ! if(it.eq.3) write(*,'(I2,1X,a,e20.12,84X,a,e20.12)')k,' 1 det',&
632 ! (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+&
633 ! mp(i,2)*(0.-tr(i,k,it))*pdtimeRG/dxpres(i,k),&
634 ! 'mp',mp(i,k)
635  else ! mp(2) = 0 = mp(1) pas de descente insaturee, rien ne se passe s'il ne pleut pas, sinon pluie->env
636  if(no_precip(i,1)) then
637  uscavtrac=0.
638  else
639  uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/rg+evap(i,k)*qpr(i,k,it)*dxpres(i,k)/rg
640  endif
641 ! 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)
642  endif
643 
644  else ! k > INB au-dessus du nuage
645  uscavtrac=0.
646  endif
647 
648 ! ===== tendances finales ======
649  trsptd(i,k,it)=trsptrac*pdtimerg/dxpres(i,k) ! td transport sans eau dans courants satures
650  dtrsscav(i,k,it)=scavtrac*pdtimerg/dxpres(i,k) ! td du lessivage dans courants satures
651  dtruscav(i,k,it)=uscavtrac*pdtimerg/dxpres(i,k) ! td courant insat
652  dtrsat(i,k,it)=(trsptrac+scavtrac)*pdtimerg/dxpres(i,k) ! td courant sat
653  dtrcv(i,k,it)=(trsptrac+scavtrac+uscavtrac)*pdtimerg/dxpres(i,k)!dtrsat(i,k,it)+dtrUscav(i,k,it) td conv
654 !!!!!!
655  dtrcvma(i,k,it)=tdcvma*pdtimerg/dxpres(i,k) ! MA tendance convection
656  ENDDO
657  ENDDO
658 
659 ! test de conservation du traceur
660 !print*,"_____________________________________________________________"
661 !print*," "
662 ! conserv=0.
663 ! conservMA=0.
664 ! DO k= klev-1,1,-1
665 ! DO i=1, klon
666 ! conserv=conserv+dtrcv(i,k,it)* &
667 ! (paprs(i,k)-paprs(i,k+1))/RG
668 ! conservMA=conservMA+dtrcvMA(i,k,it)* &
669 ! (paprs(i,k)-paprs(i,k+1))/RG
670 !
671 ! if(it.eq.3) write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)') k,&
672 ! 'MA td ',dtrcvMA(i,k,it)*dxpres(i,k)/RG,&
673 ! ' td',dtrcv(i,k,it)*dxpres(i,k)/RG,' conservMA ',conservMA,'conserv ',conserv
674 !!
675 ! ENDDO
676 ! ENDDO
677 ! if(it.eq.3) print *,'it',it,'conserv ',conserv,'conservMA ',conservMA
678 
679 END SUBROUTINE cvltr