source: LMDZ5/trunk/libf/phylmd/cvltr_scav.F90 @ 2201

Last change on this file since 2201 was 2147, checked in by idelkadi, 10 years ago

Les modifications introduites ont pour but :
1/ d'autoriser le couplage entre INCA-aerosol et les parametrisations de
la nouvelle physique (NP) de LMDZ, en particulier les thermiques et le
transport convectif,
2/ generaliser les routines de calcul de proprietes optiques des
aerosols pour RRTM au cas ou les aerosols sont interactifs
3/ d'inclure les effets LW des aerosols stratospheriques pour RRTM

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