source: LMDZ6/branches/contrails/libf/phylmd/cvltr.f90 @ 5426

Last change on this file since 5426 was 5289, checked in by abarral, 2 months ago

Turn YOECUMF.h into a module
Fix USE in fxy_new_mod_h.f90

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