source: LMDZ6/branches/LMDZ-tracers/libf/phylmd/cvltr_scav.F90 @ 3871

Last change on this file since 3871 was 3852, checked in by dcugnet, 4 years ago

Extension of the tracers management.

The tracers files can be:

1) "traceur.def": old format, with:

  • the number of tracers on the first line
  • one line for each tracer: <tracer name> <hadv> <vadv> [<parent name>]

2) "tracer.def": new format with one section each model component.
3) "tracer_<name>.def": new format with a single section.

The formats 2 and 3 reading is driven by the "type_trac" key, which can be a

coma-separated list of components.

  • Format 2: read the sections from the "tracer.def" file.
  • format 3: read one section each "tracer_<section name>.def" file.
  • the first line of a section is "&<section name>
  • the other lines start with a tracer name followed by <key>=<val> pairs.
  • the "default" tracer name is reserved ; the other tracers of the section inherit its <key>=<val>, except for the keys that are redefined locally.

This format helps keeping the tracers files compact, thanks to the "default"
special tracer and the three levels of factorization:

  • on the tracers names: a tracer name can be a coma-separated list of tracers => all the tracers of the list have the same <key>=<val> properties
  • on the parents names: the value of the "parent" property can be a coma-separated list of tracers => only possible for geographic tagging tracers
  • on the phases: the property "phases" is [g](l][s] (gas/liquid/solid)

Read information is stored in the vector "tracers(:)", of derived type "tra".

"isotopes_params.def" is a similar file, with one section each isotopes family.
It contains a database of isotopes properties ; if there are second generation
tracers (isotopes), the corresponding sections are read.

Read information is stored in the vector "isotopes(:)", of derived type "iso".

The "getKey" function helps to get the values of the parameters stored in
"tracers" or "isotopes".
File size: 28.8 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_phy, ONLY : nbtr
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!jyg<
125!!  ! ======================================================
126!!  ! calcul de l'impaction
127!!  ! ======================================================
128!!
129!!  ! impaction sur la surface de la colonne de la descente insaturee
130!!  ! On prend la moyenne des precip entre le niveau i+1 et i
131!!  ! I=3/4* (P(1+1)+P(i))/2 / (sigd*r*rho_l)
132!!  ! 1000kg/m3= densite de l'eau
133!!  ! 0.75e-3 = 3/4 /1000
134!!  ! Par la suite, I est tout le temps multiplie par sig_d pour avoir l'impaction sur la surface de la maille
135!!!!  ! on le neglige ici pour simplifier le code
136!!
137!!  DO j=1,klev-1
138!!     DO i=1,klon
139!!        imp(i,j) = coefcoli_3d(i,j)*0.75e-3/rdrop *&
140!!             0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
141!!     ENDDO
142!!  ENDDO
143!>jyg
144  !
145  ! initialisation pour flux de traceurs, td et autre
146  !
147  trsptrac = 0.
148  scavtrac = 0.
149  uscavtrac = 0.
150  qfeed(:,it) = 0.              !RL
151  DO j=1,klev
152     DO i=1,klon
153        zmfd(i,j,it)=0.
154        zmfa(i,j,it)=0.
155        zmfu(i,j,it)=0.
156        zmfp(i,j,it)=0.
157        zmfphi2(i,j,it)=0.
158        zmfd1a(i,j,it)=0.
159        zmfdam(i,j,it)=0.
160        qDi(i,j,it)=0.
161        qPr(i,j,it)=0.
162        qPa(i,j,it)=0.
163        qMel(i,j,it)=0.
164        qMeltmp(i,j,it)=0.
165        qTrdi(i,j,it)=0.
166        kappa(i,j)=0.
167        trsptd(i,j,it)=0.
168        dtrsat(i,j,it)=0.
169        dtrSscav(i,j,it)=0.
170        dtrUscav(i,j,it)=0.
171        dtrcv(i,j,it)=0.
172        dtrcvMA(i,j,it)=0.
173        evap(i,j)=0.
174        dxpres(i,j)=0.
175        qpmMint(i,j,it)=0.
176        Mint(i,j)=0.
177     END DO
178  END DO
179
180  ! suppression des valeurs très faibles (~1e-320)
181  ! multiplication de levaporation pour lavoir par unite de temps
182  ! et par unite de surface de la maille
183  ! -> cv30_unsat : evap : masse evaporee/s/(m2 de la descente)
184  DO j=1,klev
185     DO i=1,klon
186        IF(ev(i,j).lt.1.e-16) THEN
187           evap(i,j)=0.
188        ELSE
189           evap(i,j)=ev(i,j)*sigd(i)
190        ENDIF
191     END DO
192  END DO
193
194  DO j=1,klev
195     DO i=1,klon
196        IF(j.LT.klev) THEN
197           IF(epIN(i,j).LT.1.e-32) THEN
198              ep(i,j)=0.
199           ELSE
200              ep(i,j)=epIN(i,j)
201           ENDIF
202        ELSE
203           ep(i,j)=epmax
204        ENDIF
205        IF(mpIN(i,j).LT.1.e-32) THEN
206           mp(i,j)=0.
207        ELSE
208           mp(i,j)=mpIN(i,j)
209        ENDIF
210        IF(pmflxsIN(i,j).LT.1.e-32) THEN
211           pmflxs(i,j)=0.
212        ELSE
213           pmflxs(i,j)=pmflxsIN(i,j)
214        ENDIF
215        IF(pmflxrIN(i,j).LT.1.e-32) THEN
216           pmflxr(i,j)=0.
217        ELSE
218           pmflxr(i,j)=pmflxrIN(i,j)
219        ENDIF
220        IF(wdtrainA(i,j).LT.1.e-32) THEN
221           Pa(i,j)=0.
222        ELSE
223           Pa(i,j)=wdtrainA(i,j)
224        ENDIF
225        IF(wdtrainM(i,j).LT.1.e-32) THEN
226           Pm(i,j)=0.
227        ELSE
228           Pm(i,j)=wdtrainM(i,j)
229        ENDIF
230     END DO
231  END DO
232
233  !==========================================
234  DO j = klev-1,1,-1
235     DO i = 1,klon
236        NO_precip(i,j) = (pmflxr(i,j+1)+pmflxs(i,j+1)).LT.1.e-10&
237             .AND.Pa(i,j).LT.1.e-10.AND.Pm(i,j).LT.1.e-10
238     END DO
239  END DO
240
241!jyg<
242  ! ======================================================
243  ! calcul de l'impaction
244  ! ======================================================
245
246  ! impaction sur la surface de la colonne de la descente insaturee
247  ! On prend la moyenne des precip entre le niveau i+1 et i
248  ! I=3/4* (P(1+1)+P(i))/2 / (sigd*r*rho_l)
249  ! 1000kg/m3= densite de l'eau
250  ! 0.75e-3 = 3/4 /1000
251  ! Par la suite, I est tout le temps multiplie par sig_d pour avoir l'impaction sur la surface de la maille
252  ! on le neglige ici pour simplifier le code
253
254  DO j=1,klev-1
255     DO i=1,klon
256        imp(i,j) = coefcoli_3d(i,j)*0.75e-3/rdrop *&
257             0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
258     ENDDO
259  ENDDO
260!>jyg
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  !RL
333  !  Feeding concentrations
334  DO j=1,klev
335     DO i=1,klon
336        qfeed(i,it)=qfeed(i,it)+wght_cvfd(i,j)*tr(i,j,it)
337     END DO
338  END DO
339  !RL
340  !
341  DO j=1,klev
342     DO i=1,klon
343        !RL
344        !!        zmfa(i,j,it)=da(i,j)*(tr(i,1,it)-tr(i,j,it))                     ! da
345        zmfa(i,j,it)=da(i,j)*(qfeed(i,it)-tr(i,j,it))                     ! da
346        !RL
347     END DO
348  END DO
349  !
350  DO k=1,klev
351     DO j=1,klev
352        DO i=1,klon
353           zmfp(i,j,it)=zmfp(i,j,it)+phi(i,j,k)*(tr(i,k,it)-tr(i,j,it))  ! phi
354        END DO
355     END DO
356  END DO
357  ! RomP ajout des matrices liees au lessivage
358  DO j=1,klev
359     DO i=1,klon
360        zmfd1a(i,j,it)=d1a(i,j)*tr(i,1,it)   ! da1
361        zmfdam(i,j,it)=dam(i,j)*tr(i,1,it)   ! dam
362     END DO
363  END DO
364  DO k=1,klev
365     DO j=1,klev
366        DO i=1,klon
367           zmfphi2(i,j,it)=zmfphi2(i,j,it)+phi2(i,j,k)*tr(i,k,it)  ! psi
368        END DO
369     END DO
370  END DO
371  DO j=1,klev-1
372     DO i=1,klon
373        zmfu(i,j,it)=max(0.,upd(i,j+1)+dnd(i,j+1))*(tr(i,j+1,it)-tr(i,j,it))
374     END DO
375  END DO
376  DO j=2,klev
377     DO i=1,klon
378        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))
379     END DO
380  END DO
381  ! ===================================================
382  ! calcul des tendances liees aux courants insatures
383  ! ===================================================
384  !  pression 
385  DO k=1, klev
386     DO i=1, klon
387        dxpres(i,k)=paprs(i,k)-paprs(i,k+1)
388     ENDDO
389  ENDDO
390  pdtimeRG=pdtime*RG
391
392  ! q_pa et q_pm    traceurs issues des courants satures se retrouvant dans les precipitations
393  DO j=1,klev
394     DO i=1,klon
395        IF(j.GE.icb(i).AND.j.LE.inb(i)) THEN
396           IF(clw(i,j).GT.1.e-16) THEN
397              !JE           qPa(i,j,it)=ccntrAA_coef*tr(i,1,it)/clw(i,j)
398              qPa(i,j,it)=ccntrAA_3d(i,j)*tr(i,1,it)/clw(i,j)
399           ELSE
400              qPa(i,j,it)=0.
401           ENDIF
402        ENDIF
403     END DO
404  END DO
405
406  ! calcul de q_pm en 2 parties :
407  ! 1) calcul de sa valeur pour un niveau z' donne
408  ! 2) integration sur la verticale sur z'
409  DO j=1,klev
410     DO k=1,j-1
411        DO i=1,klon
412           IF(k.GE.icb(i).AND.k.LE.inb(i).AND.&
413                j.LE.inb(i)) THEN
414              IF(elij(i,k,j).GT.1.e-16) THEN
415                 !JE             qMeltmp(i,j,it)=((1-ep(i,k))*ccntrAA_coef*tr(i,1,it)&
416                 !JE                         *(1.-sij(i,k,j))  +ccntrENV_coef&
417                 !JE                         *tr(i,k,it)*sij(i,k,j)) / elij(i,k,j)
418                 qMeltmp(i,j,it)=((1-ep(i,k))*ccntrAA_3d(i,k)*tr(i,1,it)&
419                      *(1.-sij(i,k,j))  +ccntrENV_3d(i,k)&
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        !JE     scavtrac=-ccntrAA_coef*zmfd1a(i,k,it)&
613        !JE               -zmfphi2(i,k,it)*ccntrENV_coef&
614        !JE               -zmfdam(i,k,it)*ccntrAA_coef
615        scavtrac=-ccntrAA_3d(i,k)*zmfd1a(i,k,it)&
616             -zmfphi2(i,k,it)*ccntrENV_3d(i,k)&
617             -zmfdam(i,k,it)*ccntrAA_3d(i,k)
618        ! lessivage courants insatures
619        if(k.LE.inb(i).AND.k.GT.1) THEN   ! tendances dans le nuage
620           !------------------------------------------------------------- detrainement
621           if(mp(i,k+1).GT.mp(i,k).AND.mp(i,k+1).GT.1.e-10) THEN
622              uscavtrac= (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))&
623                   + mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
624              !
625              !        if(it.eq.3) write(*,'(I2,1X,a,5X,e20.12,82X,a,e20.12)')k,' det incloud',&
626              !                    (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+&
627              !                    mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k),&
628              !                    'mp',mp(i,k)
629              !------------------------------------------------------------- entrainement
630           ELSEIF(mp(i,k).GT.mp(i,k+1).AND.mp(i,k).GT.1.e-10) THEN
631              uscavtrac= mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
632              !
633              !         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)
634              !=!------------------------------------------------------------- end ent/det
635           ELSE !        mp(i,k+1)=0. et mp(i,k)=0.    pluie directement sur l environnement
636
637              if(NO_precip(i,k)) THEN
638                 uscavtrac=0.
639                 !         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)
640              ELSE
641                 uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/RG+evap(i,k)*qPr(i,k,it)*dxpres(i,k)/RG
642                 !         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)
643                 !!JE adds
644                 !         if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'imp',imp(i,k)
645                 !         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)
646                 !         if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'evap',evap(i,k)
647                 !         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)
648                 !         if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'dxpres',dxpres(i,k)
649                 !!Je end
650
651              ENDIF
652           ENDIF  ! mp/det/ent
653           !------------------------------------------------------------- premiere couche
654        ELSEIF(k.eq.1) THEN  !                                      mp(1)=0.
655           if(mp(i,2).GT.1.e-10) THEN  !detrainement
656              uscavtrac= (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it)) !&
657              !                   + mp(i,2)*(0.-tr(i,k,it))
658              !
659              !       if(it.eq.3) write(*,'(I2,1X,a,e20.12,84X,a,e20.12)')k,' 1 det',&
660              !                  (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+&
661              !                  mp(i,2)*(0.-tr(i,k,it))*pdtimeRG/dxpres(i,k),&
662              !                   'mp',mp(i,k)
663           ELSE   ! mp(2) = 0 = mp(1) pas de descente insaturee, rien ne se passe s'il ne pleut pas, sinon pluie->env
664              if(NO_precip(i,1)) THEN
665                 uscavtrac=0.
666              ELSE
667                 uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/RG+evap(i,k)*qPr(i,k,it)*dxpres(i,k)/RG
668              ENDIF
669              !         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)
670           ENDIF
671
672        ELSE   ! k > INB  au-dessus du nuage
673           uscavtrac=0.
674        ENDIF
675
676        ! =====    tendances finales  ======
677        trsptd(i,k,it)=trsptrac*pdtimeRG/dxpres(i,k)              ! td transport sans eau dans courants satures
678        dtrSscav(i,k,it)=scavtrac*pdtimeRG/dxpres(i,k)            ! td du lessivage dans courants satures
679        dtrUscav(i,k,it)=uscavtrac*pdtimeRG/dxpres(i,k)           ! td courant insat
680        dtrsat(i,k,it)=(trsptrac+scavtrac)*pdtimeRG/dxpres(i,k)   ! td courant sat
681        dtrcv(i,k,it)=(trsptrac+scavtrac+uscavtrac)*pdtimeRG/dxpres(i,k)!dtrsat(i,k,it)+dtrUscav(i,k,it)         td conv
682!!!!!!
683        dtrcvMA(i,k,it)=tdcvMA*pdtimeRG/dxpres(i,k) ! MA tendance convection
684     ENDDO
685  ENDDO
686
687  ! test de conservation du traceur
688  !print*,"_____________________________________________________________"
689  !print*,"                                                             "
690  !      conserv=0.
691  !      conservMA=0.
692  !      DO k= klev-1,1,-1
693  !        DO i=1, klon
694  !         conserv=conserv+dtrcv(i,k,it)*   &
695  !        (paprs(i,k)-paprs(i,k+1))/RG
696  !         conservMA=conservMA+dtrcvMA(i,k,it)*   &
697  !        (paprs(i,k)-paprs(i,k+1))/RG
698  !
699  !      if(it.eq.3)  write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)') k,&
700  !         'MA td ',dtrcvMA(i,k,it)*dxpres(i,k)/RG,&
701  !         ' td',dtrcv(i,k,it)*dxpres(i,k)/RG,'         conservMA ',conservMA,'conserv ',conserv
702  !!
703  !        ENDDO
704  !      ENDDO
705  !       if(it.eq.3) print *,'it',it,'conserv ',conserv,'conservMA ',conservMA
706
707END SUBROUTINE cvltr_scav
Note: See TracBrowser for help on using the repository browser.