source: LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cvltr.F90 @ 3758

Last change on this file since 3758 was 2320, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation: make an infotrac_phy module, which should be used from within the physics, and is initialized from infotrac (dynamics) via iniphysiq.
EM

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 25.6 KB
Line 
1!
2! $Id $
3!
4SUBROUTINE 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
696END SUBROUTINE cvltr
Note: See TracBrowser for help on using the repository browser.