source: LMDZ5/branches/LMDZ5-DOFOCO/libf/phylmd/cvltr.F90 @ 4400

Last change on this file since 4400 was 1759, checked in by Ehouarn Millour, 11 years ago

IOIPSL routine getin is not threadsafe, so when running in OpenMP, it should be called by only one thread (and the result copied to other threads in the case of threadprivate variables).
EM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 25.1 KB
Line 
1!
2! $Id $
3!
4SUBROUTINE 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
679END SUBROUTINE cvltr
Note: See TracBrowser for help on using the repository browser.