source: LMDZ5/trunk/libf/phylmd/cvltr_spl.F90 @ 3405

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

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

File size: 30.3 KB
RevLine 
[2147]1!
2! $Id $
3!
4SUBROUTINE cvltr_spl(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           kk,henry,zrho, ccntrAA_spla,ccntrENV_spla,coefcoli_spla, &
9           id_prec,id_fine,id_coss, id_codu, id_scdu, &
10           dtrcv,trsptd,dtrSscav,dtrsat,dtrUscav,qDi,qPr, &
11           qPa,qMel,qTrdi,dtrcvMA,Mint,                   &
12           zmfd1a,zmfphi2,zmfdam)
13  USE IOIPSL
14  USE dimphy
[2320]15  USE infotrac_phy, ONLY : nbtr,tname
[2147]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!  REAL,DIMENSION(klon,klev),INTENT(IN)    :: pplay ! pression aux 1/2 couches (bas en haut)
39  REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr     ! q de traceur (bas en haut)
40  INTEGER,INTENT(IN)                        :: it
41  REAL,DIMENSION(klon,klev),INTENT(IN)      :: upd   ! saturated updraft mass flux
42  REAL,DIMENSION(klon,klev),INTENT(IN)      :: dnd   ! saturated downdraft mass flux
43!
44  REAL,DIMENSION(klon,klev),INTENT(IN)      :: wdtrainA   ! masses precipitantes de l'asc adiab
45  REAL,DIMENSION(klon,klev),INTENT(IN)      :: wdtrainM   ! masses precipitantes des melanges
46!JE  REAL,DIMENSION(klon,klev),INTENT(IN)      :: pmflxrIN   ! vprecip: eau
47  REAL,DIMENSION(klon,klev+1),INTENT(IN)      :: pmflxrIN   ! vprecip: eau
48!JE  REAL,DIMENSION(klon,klev),INTENT(IN)      :: pmflxsIN   ! vprecip: neige
49  REAL,DIMENSION(klon,klev+1),INTENT(IN)      :: pmflxsIN   ! vprecip: neige
50  REAL,DIMENSION(klon,klev),INTENT(IN)      :: ev         ! evaporation cv30_routine
51  REAL,DIMENSION(klon,klev),INTENT(IN)      :: epIN
52  REAL,DIMENSION(klon,klev),INTENT(IN)      :: te
53  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: sij        ! fraction dair de lenv
54  REAL,DIMENSION(klon,klev),INTENT(IN)      :: wght_cvfd  ! weights of the layers feeding convection
55  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: elij       ! contenu en eau condensée spécifique/conc deau condensée massique
56  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: epmlmMm    ! eau condensee precipitee dans mel masse dair sat
57  REAL,DIMENSION(klon,klev),INTENT(IN)      :: eplaMm    ! eau condensee precipitee dans aa masse dair sat
58
59  REAL,DIMENSION(klon,klev),INTENT(IN)      :: clw        ! contenu en eau condensée dans lasc adiab
60  REAL,DIMENSION(klon),INTENT(IN)           :: sigd
61  INTEGER,DIMENSION(klon),INTENT(IN)        :: icb,inb
62! Sortie
63  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrcv     ! tendance totale (bas en haut)
64  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrcvMA ! M-A Filiberti
65  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: trsptd    ! tendance du transport
66  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrSscav  ! tendance du lessivage courant sat
67  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrsat    ! tendance trsp+sat scav
68  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrUscav  ! tendance du lessivage courant unsat
69!
70! Variables locales
71  INTEGER                           :: i,j,k
72  REAL,DIMENSION(klon,klev)         :: dxpres   ! difference de pression entre niveau (j+1) et (j)
73  REAL                              :: pdtimeRG ! pas de temps * gravite
74  REAL,DIMENSION(klon,nbtr)         :: qfeed     ! tracer concentration feeding convection
75! variables pour les courants satures
76  REAL,DIMENSION(klon,klev,klev)    :: zmd
77  REAL,DIMENSION(klon,klev,klev)    :: za
78  REAL,DIMENSION(klon,klev,nbtr)    :: zmfd,zmfa
79  REAL,DIMENSION(klon,klev,nbtr)    :: zmfp,zmfu
80
81  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: zmfd1a
82  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: zmfdam
83  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: zmfphi2
84
85! RomP ! les variables sont nettoyees des valeurs aberrantes
86  REAL,DIMENSION(klon,klev)         :: Pa, Pm  ! pluie AA et mélanges, var temporaire
87  REAL,DIMENSION(klon,klev)         :: pmflxs,pmflxr ! pmflxrIN,pmflxsIN sans valeur aberante
88  REAL,DIMENSION(klon,klev)         :: mp            ! flux de masse
89  REAL,DIMENSION(klon,klev)         :: ep            ! fraction d'eau convertie en precipitation
90  REAL,DIMENSION(klon,klev)         :: evap          ! evaporation : variable temporaire
91  REAL,DIMENSION(klon,klev)         :: rho    !environmental density
92
93  REAL,DIMENSION(klon,klev)         :: kappa ! denominateur du au calcul de la matrice
94                                             ! pour obtenir qd et qp
95  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: qTrdi ! traceurs descente air insature transport MA
96  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: qDi  ! traceurs descente insaturees
97  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: qPr  ! traceurs colonne precipitante
98  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: qPa  ! traceurs dans les precip issues lasc. adiab.
99  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: qMel ! traceurs dans les precip issues des melanges
100  REAL,DIMENSION(klon,klev,nbtr)                :: qMeltmp ! variable temporaire
101  REAL,DIMENSION(klon,klev,nbtr)                :: qpmMint
102  REAL,DIMENSION(klon,klev),INTENT(OUT)         :: Mint
103! tendances
104  REAL                              :: tdcvMA           ! terme de transport de traceur (schema Marie Angele)
105  REAL                              :: trsptrac         ! terme de transport de traceur par l'air
106  REAL                              :: scavtrac         ! terme de lessivage courant sature
107  REAL                              :: uscavtrac        ! terme de lessivage courant insature
108! impaction
109!!!       Correction apres discussion Romain P. / Olivier B.
110!!!  REAL,PARAMETER                    :: rdrop=2.5e-3     ! rayon des gouttes d'eau
111  REAL,PARAMETER                    :: rdrop=1.e-3     ! rayon des gouttes d'eau
112!!!
113  REAL,DIMENSION(klon,klev)         :: imp              ! coefficient d'impaction
114! parametres lessivage
115  REAL                    :: ccntrAA_coef        ! \alpha_a : fract aerosols de l'AA convertis en CCN
116  REAL                    :: ccntrENV_coef       ! \beta_m  : fract aerosols de l'env convertis en CCN
117  REAL                    :: coefcoli            ! coefficient de collision des gouttes sur les aerosols
118!
119  LOGICAL,DIMENSION(klon,klev) :: NO_precip
120!  LOGICAL                      :: scavON
121! var tmp tests
122  REAL                           :: conserv
123  real                           :: conservMA
124! JE SPLA adaptation
125
126   INTEGER :: id_prec,id_fine,id_coss, id_codu, id_scdu
127   REAL,DIMENSION(nbtr) :: ccntrAA_spla,ccntrENV_spla,coefcoli_spla
128   REAL,DIMENSION(klon,klev) :: ccntrAA_coef3d
129   REAL,DIMENSION(klon,klev) :: ccntrENV_coef3d
130   REAL,DIMENSION(klon,klev) :: coefcoli3d
131
132   REAL,DIMENSION(nbtr) :: henry         !--cste de Henry  mol/l/atm
133   REAL,DIMENSION(nbtr) :: kk            !--coefficient de var avec T (K)
134   REAL :: henry_t !--constante de Henry a T t  (mol/l/atm)
135   REAL :: f_a     !--rapport de la phase aqueuse a la phase gazeuse
136
137   REAL, PARAMETER :: ph=5.
138   REAL ::  K1, K2
139   REAL,DIMENSION(klon,klev) :: zrho
140   REAL, PARAMETER :: qliq=1.e-3 ! CONVECTIVE ONLY
141
142! Je end 20140105
143
144!JE20140724<<
145!JE SPLA adapt:
146!
147!
148!! coefficient lessivage
149!   ccntrAA_coef     = 0.
150!   ccntrENV_coef    = 0.
151!   coefcoli         = 0.
152!
153!!$OMP MASTER
154!  call getin('ccntrAA_coef',ccntrAA_coef)
155!  call getin('ccntrENV_coef',ccntrENV_coef)
156!  call getin('coefcoli',coefcoli)
157!!$OMP END MASTER
158!!$OMP BARRIER
159!!! JE
160!  do j=1,klev
161!   do i=1,klon
162!     imp(i,j)=0.
163!     ccntrAA_coef3d(i,j)=ccntrAA_coef
164!     ccntrENV_coef3d(i,j)=ccntrENV_coef
165!     coefcoli3d(i,j)=coefcoli
166!   enddo
167!  enddo
168!
169!
170!!! for SPLA
171!!!JEend
172!  print*,'cvltr coef lessivage convectif', ccntrAA_coef,ccntrENV_coef,coefcoli
173!
174!JE20140724>>
175
176!  scavON=.TRUE.
177!  if(scavON) then
178!   ccntrAA_coef     = 1.
179!   ccntrENV_coef    = 1.
180!   coefcoli         = 1.
181!  else
182!   ccntrAA_coef     = 0.
183!   ccntrENV_coef    = 0.
184!   coefcoli         = 0.
185!  endif
186
187! ======================================================
188! calcul de l'impaction
189! ======================================================
190!initialisation
191  do j=1,klev
192   do i=1,klon
193     imp(i,j)=0.
194   enddo
195  enddo
196!JE init 20140103
197!! impaction sur la surface de la colonne de la descente insaturee
198!! On prend la moyenne des precip entre le niveau i+1 et i
199!! I=3/4* (P(1+1)+P(i))/2 / (sigd*r*rho_l)
200!!  1000kg/m3= densité de l'eau
201!! 0.75e-3 = 3/4 /1000
202!! Par la suite, I est tout le temps multiplié par sig_d pour avoir l'impaction sur la surface de la maille
203!! on le néglige ici pour simplifier le code
204!  do j=1,klev-1
205!   do i=1,klon
206!    imp(i,j) = coefcoli*0.75e-3/rdrop *&
207!             0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
208!!    rho(i,j)=pplay(i,j)/(rd*te(i,j))
209!   enddo
210!  enddo
211!JE end 20140103
212
213!
214! initialisation pour flux de traceurs, td et autre
215   trsptrac = 0.
216   scavtrac = 0.
217   uscavtrac = 0.
218   qfeed(:,it) = 0.              !RL
219  DO j=1,klev
220   DO i=1,klon
221    zmfd(i,j,it)=0.
222    zmfa(i,j,it)=0.
223    zmfu(i,j,it)=0.
224    zmfp(i,j,it)=0.
225    zmfphi2(i,j,it)=0.
226    zmfd1a(i,j,it)=0.
227    zmfdam(i,j,it)=0.
228    qDi(i,j,it)=0.
229    qPr(i,j,it)=0.
230    qPa(i,j,it)=0.
231    qMel(i,j,it)=0.
232    qMeltmp(i,j,it)=0.
233    qTrdi(i,j,it)=0.
234    kappa(i,j)=0.
235    trsptd(i,j,it)=0.
236    dtrsat(i,j,it)=0.
237    dtrSscav(i,j,it)=0.
238    dtrUscav(i,j,it)=0.
239    dtrcv(i,j,it)=0.
240    dtrcvMA(i,j,it)=0.
241    evap(i,j)=0.
242    dxpres(i,j)=0.
243    qpmMint(i,j,it)=0.
244    Mint(i,j)=0.
245   END DO
246  END DO
247
248! suppression des valeurs très faibles (~1e-320)
249! multiplication de levaporation pour lavoir par unite de temps
250! et par unite de surface de la maille
251! -> cv30_unsat : evap : masse evaporee/s/(m2 de la descente)
252  DO j=1,klev
253   DO i=1,klon
254    if(ev(i,j).lt.1.e-16) then
255     evap(i,j)=0.
256    else
257     evap(i,j)=ev(i,j)*sigd(i)
258    endif
259   END DO
260  END DO
261
262  DO j=1,klev
263   DO i=1,klon
264   if(j.lt.klev) then
265    if(epIN(i,j).lt.1.e-32) then
266     ep(i,j)=0.
267    else
268     ep(i,j)=epIN(i,j)
269    endif
270   else
271    ep(i,j)=epmax
272   endif
273    if(mpIN(i,j).lt.1.e-32) then
274     mp(i,j)=0.
275    else
276     mp(i,j)=mpIN(i,j)
277    endif
278    if(pmflxsIN(i,j).lt.1.e-32) then
279     pmflxs(i,j)=0.
280    else
281     pmflxs(i,j)=pmflxsIN(i,j)
282    endif
283    if(pmflxrIN(i,j).lt.1.e-32) then
284     pmflxr(i,j)=0.
285    else
286     pmflxr(i,j)=pmflxrIN(i,j)
287    endif
288    if(wdtrainA(i,j).lt.1.e-32) then
289     Pa(i,j)=0.
290    else
291     Pa(i,j)=wdtrainA(i,j)
292    endif
293    if(wdtrainM(i,j).lt.1.e-32) then
294     Pm(i,j)=0.
295    else
296     Pm(i,j)=wdtrainM(i,j)
297    endif
298   END DO
299  END DO
300
301!==========================================
302  DO j = klev-1,1,-1
303   DO i = 1,klon
304     NO_precip(i,j) = (pmflxr(i,j+1)+pmflxs(i,j+1)).lt.1.e-10&
305                    .and.Pa(i,j).lt.1.e-10.and.Pm(i,j).lt.1.e-10
306   END DO
307  END DO
308!==============================================================================
309! JE SPLA: Calc of ccntrAA_coef,ccntrENV_coef, coefcoli for SPLA aerosols and
310! precursors. From SPLA code inscav_spl.F
311!print *,'Using SPLA values for cvltr_spl ccntr and coefcoli params'
312!
313!
314!IF (it.EQ.2) THEN !--fine aerosol
315!  DO j=1,klev
316!   DO i=1,klon
317!    ccntrAA_coef3d(i,j)=0.7
318!    ccntrENV_coef3d(i,j)=0.7
319!    coefcoli3d(i,j)=0.001
320!   ENDDO
321!  ENDDO
322!ELSEIF (it.EQ.3) THEN !-- Coarse Sea salt aerosol
323!  DO j=1,klev
324!   DO i=1,klon
325!    ccntrAA_coef3d(i,j)=1.0
326!    ccntrENV_coef3d(i,j)=1.0
327!    coefcoli3d(i,j)=0.001
328!   ENDDO
329!  ENDDO
330!
331!ELSEIF (it.EQ.4) THEN !--Coarse Dust aerosol
332!  DO j=1,klev
333 !  DO i=1,klon
334!    ccntrAA_coef3d(i,j)=0.7
335!    ccntrENV_coef3d(i,j)=0.7
336!    coefcoli3d(i,j)=0.001
337!
338!   ENDDO
339!  ENDDO
340! Gas precursor: Henry's law
341
342IF (it .EQ. id_prec) THEN
343      DO k=1, klev
344      DO i=1, klon
345        henry_t=henry(it)*exp(-kk(it)*(1./298.-1./te(i,k)))    !--mol/l/atm
346        K1=1.2e-2*exp(-2010*(1/298.-1/te(i,k)))
347        K2=6.6e-8*exp(-1510*(1/298.-1/te(i,k)))
348        henry_t=henry_t*(1. + K1/10.**(-ph) + K1*K2/(10.**(-ph))**2)
349        f_a=henry_t/101.325*R*te(i,k)*qliq*zrho(i,k)/rho_water
350 !       scav(i,k)=f_a/(1.+f_a)
351        ccntrAA_coef3d(i,k)= f_a/(1.+f_a)
352        ccntrENV_coef3d(i,k)= f_a/(1.+f_a)
353        coefcoli3d(i,k)=0.0
354      ENDDO
355      ENDDO
356 !    CALL minmaxqfi2(clw,1.e33,-1.e33,'minmax clw')
357ELSE
358  DO j=1,klev
359   DO i=1,klon
360    ccntrAA_coef3d(i,j)=ccntrAA_spla(it)
361    ccntrENV_coef3d(i,j)=ccntrENV_spla(it)
362    coefcoli3d(i,j)=coefcoli_spla(it)
363   ENDDO
364  ENDDO
365
366
367ENDIF
368
369! JE end SPLA modifs in ccn parameters
370!==============================================================================
371
372
373
374
375!JE init 20140103
376! impaction sur la surface de la colonne de la descente insaturee
377! On prend la moyenne des precip entre le niveau i+1 et i
378! I=3/4* (P(1+1)+P(i))/2 / (sigd*r*rho_l)
379!  1000kg/m3= densité de l'eau
380! 0.75e-3 = 3/4 /1000
381! Par la suite, I est tout le temps multiplié par sig_d pour avoir l'impaction sur la surface de la maille
382! on le néglige ici pour simplifier le code
383  do j=1,klev-1
384   do i=1,klon
385!JE    imp(i,j) = coefcoli*0.75e-3/rdrop *&
386!JE             0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
387    imp(i,j) = coefcoli3d(i,j)*0.75e-3/rdrop *&
388             0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
389!    rho(i,j)=pplay(i,j)/(rd*te(i,j))
390   enddo
391  enddo
392!JE end 20140103
393
394
395! =========================================
396! calcul des tendances liees au downdraft
397! =========================================
398!cdir collapse
399  DO k=1,klev
400   DO j=1,klev
401    DO i=1,klon
402     zmd(i,j,k)=0.
403     za (i,j,k)=0.
404    END DO
405   END DO
406  END DO
407! calcul de la matrice d echange
408! matrice de distribution de la masse entrainee en k
409! commmentaire RomP : mp > 0
410  DO k=1,klev-1
411     DO i=1,klon
412        zmd(i,k,k)=max(0.,mp(i,k)-mp(i,k+1))   ! ~ mk(k)
413     END DO
414  END DO
415  DO k=2,klev
416     DO j=k-1,1,-1
417        DO i=1,klon
418           if(mp(i,j+1).gt.1.e-10) then
419              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)
420           ENDif
421        END DO
422     END DO
423  END DO
424  DO k=1,klev
425     DO j=1,klev-1
426        DO i=1,klon
427           za(i,j,k)=max(0.,zmd(i,j+1,k)-zmd(i,j,k))
428        END DO
429     END DO
430  END DO
431!!!!! quantite  de traceur dans la descente d'air insaturee  :   4 juin 2012
432  DO k=1,klev
433     DO j=1,klev-1
434        DO i=1,klon
435        if(mp(i,j+1).gt.1.e-10) then
436          qTrdi(i,j+1,it)=qTrdi(i,j+1,it)+(zmd(i,j+1,k)/mp(i,j+1))*tr(i,k,it)
437        else
438          qTrdi(i,j,it)=0.!tr(i,j,it)
439        endif
440        ENDDO
441     ENDDO
442  ENDDO
443!!!!!
444!
445! rajout du terme lie a l ascendance induite
446!
447  DO j=2,klev
448     DO i=1,klon
449        za(i,j,j-1)=za(i,j,j-1)+mp(i,j)
450     END DO
451  END DO
452!
453! tendance courants insatures  ! sans lessivage ancien schema
454!
455  DO k=1,klev
456     DO j=1,klev
457        DO i=1,klon
458           zmfd(i,j,it)=zmfd(i,j,it)+za(i,j,k)*(tr(i,k,it)-tr(i,j,it))
459        END DO
460     END DO
461  END DO
462!
463! =========================================
464! calcul des tendances liees aux courants satures   j <-> z ; k <-> z'
465! =========================================
466!RL
467!  Feeding concentrations
468  DO j=1,klev
469     DO i=1,klon
470        qfeed(i,it)=qfeed(i,it)+wght_cvfd(i,j)*tr(i,j,it)
471     END DO
472  END DO
473!RL
474!
475  DO j=1,klev
476     DO i=1,klon
477!RL
478!!        zmfa(i,j,it)=da(i,j)*(tr(i,1,it)-tr(i,j,it))                     ! da
479        zmfa(i,j,it)=da(i,j)*(qfeed(i,it)-tr(i,j,it))                     ! da
480!RL
481     END DO
482  END DO
483!
484  DO k=1,klev
485     DO j=1,klev
486        DO i=1,klon
487           zmfp(i,j,it)=zmfp(i,j,it)+phi(i,j,k)*(tr(i,k,it)-tr(i,j,it))  ! phi
488        END DO
489     END DO
490  END DO
491! RomP ajout des matrices liees au lessivage
492  DO j=1,klev
493     DO i=1,klon
494        zmfd1a(i,j,it)=d1a(i,j)*tr(i,1,it)   ! da1
495        zmfdam(i,j,it)=dam(i,j)*tr(i,1,it)   ! dam
496     END DO
497  END DO
498  DO k=1,klev
499     DO j=1,klev
500        DO i=1,klon
501          zmfphi2(i,j,it)=zmfphi2(i,j,it)+phi2(i,j,k)*tr(i,k,it)  ! psi
502        END DO
503     END DO
504  END DO
505  DO j=1,klev-1
506     DO i=1,klon
507        zmfu(i,j,it)=max(0.,upd(i,j+1)+dnd(i,j+1))*(tr(i,j+1,it)-tr(i,j,it))
508     END DO
509  END DO
510  DO j=2,klev
511     DO i=1,klon
512        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))
513     END DO
514  END DO
515! ===================================================
516! calcul des tendances liees aux courants insatures
517! ===================================================
518!  pression 
519  DO k=1, klev
520     DO i=1, klon
521        dxpres(i,k)=paprs(i,k)-paprs(i,k+1)
522     ENDDO
523  ENDDO
524  pdtimeRG=pdtime*RG
525
526! q_pa et q_pm    traceurs issues des courants satures se retrouvant dans les precipitations
527  DO j=1,klev
528     DO i=1,klon
529        if(j.ge.icb(i).and.j.le.inb(i)) then
530          if(clw(i,j).gt.1.e-16) then
531!           qPa(i,j,it)=ccntrAA_coef*tr(i,1,it)/clw(i,j)
532           qPa(i,j,it)=ccntrAA_coef3d(i,j)*tr(i,1,it)/clw(i,j)
533          else
534           qPa(i,j,it)=0.
535          endif
536        endif
537     END DO
538  END DO
539
540! calcul de q_pm en 2 parties :
541! 1) calcul de sa valeur pour un niveau z' donne
542! 2) integration sur la verticale sur z'
543     DO j=1,klev
544      DO k=1,j-1
545        DO i=1,klon
546        if(k.ge.icb(i).and.k.le.inb(i).and.&
547           j.le.inb(i)) then
548            if(elij(i,k,j).gt.1.e-16) then
549!JE             qMeltmp(i,j,it)=((1-ep(i,k))*ccntrAA_coef*tr(i,1,it)&
550!JE                         *(1.-sij(i,k,j))  +ccntrENV_coef&
551!JE                         *tr(i,k,it)*sij(i,k,j)) / elij(i,k,j)
552             qMeltmp(i,j,it)=((1-ep(i,k))*ccntrAA_coef3d(i,k)*tr(i,1,it)&
553                         *(1.-sij(i,k,j))  +ccntrENV_coef3d(i,k)&
554                         *tr(i,k,it)*sij(i,k,j)) / elij(i,k,j)
555            else
556             qMeltmp(i,j,it)=0.
557            endif
558          qpmMint(i,j,it)=qpmMint(i,j,it) + qMeltmp(i,j,it)*epmlmMm(i,j,k)
559          Mint(i,j)=Mint(i,j) + epmlmMm(i,j,k)
560        endif ! end if dans nuage
561        END DO
562     END DO
563  END DO
564
565     DO j=1,klev
566        DO i=1,klon
567          if(Mint(i,j).gt.1.e-16) then
568           qMel(i,j,it)=qpmMint(i,j,it)/Mint(i,j)
569          else
570           qMel(i,j,it)=0.
571          endif
572     END DO
573  END DO
574
575! calcul de q_d et q_p    traceurs de la descente precipitante
576   DO j=klev-1,1,-1
577    DO i=1,klon
578     if(mp(i,j+1).gt.mp(i,j).and.mp(i,j+1).gt.1.e-10) then  ! detrainement
579       kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
580                (-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))&
581                + (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j)))
582             
583     elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then! entrainement
584       if(j.eq.1) then
585        kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
586                 (-mp(i,2)-imp(i,j)/RG*dxpres(i,j))&
587                 + (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j)))
588       else
589        kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
590                 (-mp(i,j)-imp(i,j)/RG*dxpres(i,j))&
591                 + (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j)))
592       endif
593      else
594        kappa(i,j)=1.
595      endif
596    ENDDO
597   ENDDO
598
599  DO j=klev-1,1,-1
600   DO i=1,klon
601    if (abs(kappa(i,j)).lt.1.e-25) then    !si denominateur nul (il peut y avoir des mp!=0)
602     kappa(i,j)=1.
603     if(j.eq.1) then
604       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
605     elseif(mp(i,j+1).gt.mp(i,j).and.mp(i,j+1).gt.1.e-10) then
606       qDi(i,j,it)=qDi(i,j+1,it)
607     elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then! entrainement
608       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))
609     else  ! si  mp (i)=0 et mp(j+1)=0
610       qDi(i,j,it)=tr(i,j,it) ! orig 0.
611     endif
612
613      if(NO_precip(i,j)) then
614       qPr(i,j,it)=0.
615      else
616       qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
617                   Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)&
618                   +imp(i,j)/RG*dxpres(i,j)*qDi(i,j,it))/&
619                  (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
620      endif
621    else   !     denominateur non nul
622     kappa(i,j)=1./kappa(i,j)     
623! calcul de qd et qp
624!!jyg  (20130119) correction pour le sommet du nuage
625!!     if(j.ge.inb(i)) then       !au-dessus du nuage, sommet inclu
626     if(j.gt.inb(i)) then       !au-dessus du nuage
627       qDi(i,j,it)=tr(i,j,it)   ! pas de descente => environnement = descente insaturee
628       qPr(i,j,it)=0.
629
630!  vvv premiere couche du modele ou mp(1)=0  ! det tout le temps  vvv
631     elseif(j.eq.1) then
632      if(mp(i,2).gt.1.e-10) then !mp(2) non nul -> detrainement (car mp(1) = 0) !ent pas possible
633       if(NO_precip(i,j)) then !pas de precip en (i)
634        qDi(i,j,it)=qDi(i,j+1,it)
635        qPr(i,j,it)=0.
636       else
637        qDi(i,j,it)=kappa(i,j)*(&
638            (-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
639            Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +&
640            (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
641            (-mp(i,j+1)*qDi(i,j+1,it)))
642
643        qPr(i,j,it)=kappa(i,j)*(&
644            (-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))*&
645            ((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
646            Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))&
647            +(-mp(i,j+1)*qDi(i,j+1,it)) * (imp(i,j)/RG*dxpres(i,j)))
648       endif
649
650      else !mp(2) nul -> plus de descente insaturee -> pluie agit sur environnement
651        qDi(i,j,it)=tr(i,j,it) ! orig 0.
652        if(NO_precip(i,j)) then
653         qPr(i,j,it)=0.
654        else
655         qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
656                   Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)&
657                   +imp(i,j)/RG*dxpres(i,j)*tr(i,j,it))/&
658                  (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
659        endif
660
661      endif  !mp(2) nul ou non
662
663! vvv  (j!=1.and.j.lt.inb(i))  en-dessous du sommet nuage   vvv
664     else
665!------------------------------------------------------------- detrainement
666      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
667         if(NO_precip(i,j)) then
668          qDi(i,j,it)=qDi(i,j+1,it)
669          qPr(i,j,it)=0.
670         else
671          qDi(i,j,it)=kappa(i,j)*(&
672            (-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
673            Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +&
674            (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
675            (-mp(i,j+1)*qDi(i,j+1,it)))
676!
677          qPr(i,j,it)=kappa(i,j)*(&
678            (-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))*&
679            ((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
680            Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))&
681            +(-mp(i,j+1)*qDi(i,j+1,it)) * (imp(i,j)/RG*dxpres(i,j)))
682         endif !precip
683!------------------------------------------------------------- entrainement
684      elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then
685         if(NO_precip(i,j)) then
686          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))
687          qPr(i,j,it)=0.
688         else
689          qDi(i,j,it)=kappa(i,j)*(&
690            (-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
691            Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +&
692            (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
693            (-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it)))
694!
695          qPr(i,j,it)=kappa(i,j)*(&
696            (-mp(i,j)-imp(i,j)/RG*dxpres(i,j))*&
697            ((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
698            Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))&
699            +(-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it))*&
700            (imp(i,j)/RG*dxpres(i,j)))
701         endif !precip
702!------------------------------------------------------------- endif ! ent/det
703      else  !mp nul
704        qDi(i,j,it)=tr(i,j,it) ! orig 0.
705        if(NO_precip(i,j)) then
706         qPr(i,j,it)=0.
707        else
708         qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
709                   Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)&
710                   +imp(i,j)/RG*dxpres(i,j)*tr(i,j,it))/&
711                  (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
712        endif
713      endif ! mp nul ou non
714     endif ! condition sur j
715    endif ! kappa
716   ENDDO
717  ENDDO
718
719!! print test descente insaturee
720!  DO j=klev,1,-1
721!   DO i=1,klon
722!     if(it.eq.3) then
723!   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,&
724!!         'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),&
725!          'zmfdam+zmfpsi',zmfdam(i,j,it)+zmfphi2(i,j,it),'qpmMint',qpmMint(i,j,it),&
726!          'Pm',Pm(i,j),'Mint',Mint(i,j),&
727!!          'zmfa',zmfa(i,j,it),'zmfp',zmfp(i,j,it),&
728!        'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),'zmfd1a',zmfd1a(i,j,it)
729!!          'Pa',Pa(i,j),'eplaMm',eplaMm(i,j)
730!!         'zmfd1a=da1*qa',zmfd1a(i,j,it),'Pa*qPa',wdtrainA(i,j)*qPa(i,j,it),'da1',d1a(i,j)
731!     endif
732!   ENDDO
733!  ENDDO
734
735
736! ===================================================
737! calcul final des tendances
738! ===================================================
739
740  DO k=klev-1,1,-1
741    DO i=1, klon
742! transport
743     tdcvMA=zmfd(i,k,it)+zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it)   ! double comptage des downdraft insatures
744     trsptrac=zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it)
745! lessivage courants satures
746!JE     scavtrac=-ccntrAA_coef*zmfd1a(i,k,it)&
747!JE               -zmfphi2(i,k,it)*ccntrENV_coef&
748!JE               -zmfdam(i,k,it)*ccntrAA_coef
749     scavtrac=-ccntrAA_coef3d(i,k)*zmfd1a(i,k,it)&
750               -zmfphi2(i,k,it)*ccntrENV_coef3d(i,k)&
751               -zmfdam(i,k,it)*ccntrAA_coef3d(i,k)
752! lessivage courants insatures
753   if(k.le.inb(i).and.k.gt.1) then   ! tendances dans le nuage
754!------------------------------------------------------------- detrainement
755      if(mp(i,k+1).gt.mp(i,k).and.mp(i,k+1).gt.1.e-10) then
756        uscavtrac= (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))&
757                   + mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
758!
759!        if(it.eq.3) write(*,'(I2,1X,a,5X,e20.12,82X,a,e20.12)')k,' det incloud',&
760!                    (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+&
761!                    mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k),&
762!                    'mp',mp(i,k)
763!------------------------------------------------------------- entrainement
764      elseif(mp(i,k).gt.mp(i,k+1).and.mp(i,k).gt.1.e-10) then
765         uscavtrac= mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
766!
767!         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)
768!=!------------------------------------------------------------- end ent/det
769      else !        mp(i,k+1)=0. et mp(i,k)=0.    pluie directement sur l environnement
770
771        if(NO_precip(i,k)) then
772          uscavtrac=0.
773!         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)
774        else
775          uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/RG+evap(i,k)*qPr(i,k,it)*dxpres(i,k)/RG
776!         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)
777!!JE adds
778!         if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'imp',imp(i,k)
779!         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)
780!         if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'evap',evap(i,k)
781!         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)
782!         if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'dxpres',dxpres(i,k)
783!!Je end
784
785        endif
786      endif  ! mp/det/ent
787!------------------------------------------------------------- premiere couche
788   elseif(k.eq.1) then  !                                      mp(1)=0.
789      if(mp(i,2).gt.1.e-10) then  !detrainement
790         uscavtrac= (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it)) !&
791!                   + mp(i,2)*(0.-tr(i,k,it))
792!
793!       if(it.eq.3) write(*,'(I2,1X,a,e20.12,84X,a,e20.12)')k,' 1 det',&
794!                  (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+&
795!                  mp(i,2)*(0.-tr(i,k,it))*pdtimeRG/dxpres(i,k),&
796!                   'mp',mp(i,k)
797      else   ! mp(2) = 0 = mp(1) pas de descente insaturee, rien ne se passe s'il ne pleut pas, sinon pluie->env
798        if(NO_precip(i,1)) then
799          uscavtrac=0.
800        else
801          uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/RG+evap(i,k)*qPr(i,k,it)*dxpres(i,k)/RG
802        endif
803!         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)
804      endif
805
806   else   ! k > INB  au-dessus du nuage
807    uscavtrac=0.
808   endif
809
810! =====    tendances finales  ======
811     trsptd(i,k,it)=trsptrac*pdtimeRG/dxpres(i,k)              ! td transport sans eau dans courants satures
812     dtrSscav(i,k,it)=scavtrac*pdtimeRG/dxpres(i,k)            ! td du lessivage dans courants satures
813     dtrUscav(i,k,it)=uscavtrac*pdtimeRG/dxpres(i,k)           ! td courant insat
814     dtrsat(i,k,it)=(trsptrac+scavtrac)*pdtimeRG/dxpres(i,k)   ! td courant sat
815     dtrcv(i,k,it)=(trsptrac+scavtrac+uscavtrac)*pdtimeRG/dxpres(i,k)!dtrsat(i,k,it)+dtrUscav(i,k,it)         td conv
816!!!!!!
817     dtrcvMA(i,k,it)=tdcvMA*pdtimeRG/dxpres(i,k) ! MA tendance convection
818     ENDDO
819  ENDDO
820
821! test de conservation du traceur
822!print*,"_____________________________________________________________"
823!print*,"                                                             "
824!      conserv=0.
825!      conservMA=0.
826!      DO k= klev-1,1,-1
827!        DO i=1, klon
828!         conserv=conserv+dtrcv(i,k,it)*   &
829!        (paprs(i,k)-paprs(i,k+1))/RG
830!         conservMA=conservMA+dtrcvMA(i,k,it)*   &
831!        (paprs(i,k)-paprs(i,k+1))/RG
832!
833!      if(it.eq.3)  write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)') k,&
834!         'MA td ',dtrcvMA(i,k,it)*dxpres(i,k)/RG,&
835!         ' td',dtrcv(i,k,it)*dxpres(i,k)/RG,'         conservMA ',conservMA,'conserv ',conserv
836!!
837!        ENDDO
838!      ENDDO
839!       if(it.eq.3) print *,'it',it,'conserv ',conserv,'conservMA ',conservMA
840
841END SUBROUTINE cvltr_spl
Note: See TracBrowser for help on using the repository browser.