source: LMDZ6/branches/Amaury_dev/libf/phylmd/cvltr_spl.F90 @ 5103

Last change on this file since 5103 was 5103, checked in by abarral, 4 months ago

Handle CPP_INLANDSIS in lmdz_cppkeys_wrapper.F90
Remove obsolete key wrgrads_thermcell, _ADV_HALO, _ADV_HALLO, isminmax
Remove redundant uses of CPPKEY_INCA (thanks acozic)
Remove obsolete misc/write_field.F90
Remove unused ioipsl_* wrappers
Remove calls to WriteField_u with wrong signature
Convert .F -> .[fF]90
(lint) uppercase fortran operators
[note: 1d and iso still broken - working on it]
File size: 30.1 KB
RevLine 
[5099]1
[2147]2! $Id $
[5099]3
[2147]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
[5101]15  USE infotrac_phy, ONLY: nbtr
[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
[5099]35
[2147]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
[5099]43
[2147]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
[5099]69
[2147]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
[5099]118
[2147]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:
[5099]146
147
[2147]148!! coefficient lessivage
149!   ccntrAA_coef     = 0.
150!   ccntrENV_coef    = 0.
151!   coefcoli         = 0.
[5099]152
[2147]153!!$OMP MASTER
[5101]154!  CALL getin('ccntrAA_coef',ccntrAA_coef)
155!  CALL getin('ccntrENV_coef',ccntrENV_coef)
156!  CALL getin('coefcoli',coefcoli)
[2147]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
[5099]168
169
[2147]170!!! for SPLA
171!!!JEend
[5103]172!  PRINT*,'cvltr coef lessivage convectif', ccntrAA_coef,ccntrENV_coef,coefcoli
[5099]173
[2147]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! initialisation pour flux de traceurs, td et autre
214   trsptrac = 0.
215   scavtrac = 0.
216   uscavtrac = 0.
217   qfeed(:,it) = 0.              !RL
218  DO j=1,klev
219   DO i=1,klon
220    zmfd(i,j,it)=0.
221    zmfa(i,j,it)=0.
222    zmfu(i,j,it)=0.
223    zmfp(i,j,it)=0.
224    zmfphi2(i,j,it)=0.
225    zmfd1a(i,j,it)=0.
226    zmfdam(i,j,it)=0.
227    qDi(i,j,it)=0.
228    qPr(i,j,it)=0.
229    qPa(i,j,it)=0.
230    qMel(i,j,it)=0.
231    qMeltmp(i,j,it)=0.
232    qTrdi(i,j,it)=0.
233    kappa(i,j)=0.
234    trsptd(i,j,it)=0.
235    dtrsat(i,j,it)=0.
236    dtrSscav(i,j,it)=0.
237    dtrUscav(i,j,it)=0.
238    dtrcv(i,j,it)=0.
239    dtrcvMA(i,j,it)=0.
240    evap(i,j)=0.
241    dxpres(i,j)=0.
242    qpmMint(i,j,it)=0.
243    Mint(i,j)=0.
244   END DO
245  END DO
246
247! suppression des valeurs très faibles (~1e-320)
248! multiplication de levaporation pour lavoir par unite de temps
249! et par unite de surface de la maille
250! -> cv30_unsat : evap : masse evaporee/s/(m2 de la descente)
251  DO j=1,klev
252   DO i=1,klon
[5082]253    if(ev(i,j)<1.e-16) then
[2147]254     evap(i,j)=0.
255    else
256     evap(i,j)=ev(i,j)*sigd(i)
257    endif
258   END DO
259  END DO
260
261  DO j=1,klev
262   DO i=1,klon
[5082]263   if(j<klev) then
264    if(epIN(i,j)<1.e-32) then
[2147]265     ep(i,j)=0.
266    else
267     ep(i,j)=epIN(i,j)
268    endif
269   else
270    ep(i,j)=epmax
271   endif
[5082]272    if(mpIN(i,j)<1.e-32) then
[2147]273     mp(i,j)=0.
274    else
275     mp(i,j)=mpIN(i,j)
276    endif
[5082]277    if(pmflxsIN(i,j)<1.e-32) then
[2147]278     pmflxs(i,j)=0.
279    else
280     pmflxs(i,j)=pmflxsIN(i,j)
281    endif
[5082]282    if(pmflxrIN(i,j)<1.e-32) then
[2147]283     pmflxr(i,j)=0.
284    else
285     pmflxr(i,j)=pmflxrIN(i,j)
286    endif
[5082]287    if(wdtrainA(i,j)<1.e-32) then
[2147]288     Pa(i,j)=0.
289    else
290     Pa(i,j)=wdtrainA(i,j)
291    endif
[5082]292    if(wdtrainM(i,j)<1.e-32) then
[2147]293     Pm(i,j)=0.
294    else
295     Pm(i,j)=wdtrainM(i,j)
296    endif
297   END DO
298  END DO
299
300!==========================================
301  DO j = klev-1,1,-1
302   DO i = 1,klon
[5082]303     NO_precip(i,j) = (pmflxr(i,j+1)+pmflxs(i,j+1))<1.e-10&
304                    .and.Pa(i,j)<1.e-10.and.Pm(i,j)<1.e-10
[2147]305   END DO
306  END DO
307!==============================================================================
308! JE SPLA: Calc of ccntrAA_coef,ccntrENV_coef, coefcoli for SPLA aerosols and
309! precursors. From SPLA code inscav_spl.F
310!print *,'Using SPLA values for cvltr_spl ccntr and coefcoli params'
[5099]311
312
[2147]313!IF (it.EQ.2) THEN !--fine aerosol
314!  DO j=1,klev
315!   DO i=1,klon
316!    ccntrAA_coef3d(i,j)=0.7
317!    ccntrENV_coef3d(i,j)=0.7
318!    coefcoli3d(i,j)=0.001
319!   ENDDO
320!  ENDDO
321!ELSEIF (it.EQ.3) THEN !-- Coarse Sea salt aerosol
322!  DO j=1,klev
323!   DO i=1,klon
324!    ccntrAA_coef3d(i,j)=1.0
325!    ccntrENV_coef3d(i,j)=1.0
326!    coefcoli3d(i,j)=0.001
327!   ENDDO
328!  ENDDO
[5099]329
[2147]330!ELSEIF (it.EQ.4) THEN !--Coarse Dust aerosol
331!  DO j=1,klev
332 !  DO i=1,klon
333!    ccntrAA_coef3d(i,j)=0.7
334!    ccntrENV_coef3d(i,j)=0.7
335!    coefcoli3d(i,j)=0.001
[5099]336
[2147]337!   ENDDO
338!  ENDDO
339! Gas precursor: Henry's law
340
[5082]341IF (it == id_prec) THEN
[2147]342      DO k=1, klev
343      DO i=1, klon
344        henry_t=henry(it)*exp(-kk(it)*(1./298.-1./te(i,k)))    !--mol/l/atm
345        K1=1.2e-2*exp(-2010*(1/298.-1/te(i,k)))
346        K2=6.6e-8*exp(-1510*(1/298.-1/te(i,k)))
347        henry_t=henry_t*(1. + K1/10.**(-ph) + K1*K2/(10.**(-ph))**2)
348        f_a=henry_t/101.325*R*te(i,k)*qliq*zrho(i,k)/rho_water
349 !       scav(i,k)=f_a/(1.+f_a)
350        ccntrAA_coef3d(i,k)= f_a/(1.+f_a)
351        ccntrENV_coef3d(i,k)= f_a/(1.+f_a)
352        coefcoli3d(i,k)=0.0
353      ENDDO
354      ENDDO
355 !    CALL minmaxqfi2(clw,1.e33,-1.e33,'minmax clw')
356ELSE
357  DO j=1,klev
358   DO i=1,klon
359    ccntrAA_coef3d(i,j)=ccntrAA_spla(it)
360    ccntrENV_coef3d(i,j)=ccntrENV_spla(it)
361    coefcoli3d(i,j)=coefcoli_spla(it)
362   ENDDO
363  ENDDO
364
365
366ENDIF
367
368! JE end SPLA modifs in ccn parameters
369!==============================================================================
370
371
372
373
374!JE init 20140103
375! impaction sur la surface de la colonne de la descente insaturee
376! On prend la moyenne des precip entre le niveau i+1 et i
377! I=3/4* (P(1+1)+P(i))/2 / (sigd*r*rho_l)
378!  1000kg/m3= densité de l'eau
379! 0.75e-3 = 3/4 /1000
380! Par la suite, I est tout le temps multiplié par sig_d pour avoir l'impaction sur la surface de la maille
381! on le néglige ici pour simplifier le code
382  do j=1,klev-1
383   do i=1,klon
384!JE    imp(i,j) = coefcoli*0.75e-3/rdrop *&
385!JE             0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
386    imp(i,j) = coefcoli3d(i,j)*0.75e-3/rdrop *&
387             0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
388!    rho(i,j)=pplay(i,j)/(rd*te(i,j))
389   enddo
390  enddo
391!JE end 20140103
392
393
394! =========================================
395! calcul des tendances liees au downdraft
396! =========================================
397!cdir collapse
398  DO k=1,klev
399   DO j=1,klev
400    DO i=1,klon
401     zmd(i,j,k)=0.
402     za (i,j,k)=0.
403    END DO
404   END DO
405  END DO
406! calcul de la matrice d echange
407! matrice de distribution de la masse entrainee en k
408! commmentaire RomP : mp > 0
409  DO k=1,klev-1
410     DO i=1,klon
411        zmd(i,k,k)=max(0.,mp(i,k)-mp(i,k+1))   ! ~ mk(k)
412     END DO
413  END DO
414  DO k=2,klev
415     DO j=k-1,1,-1
416        DO i=1,klon
[5082]417           if(mp(i,j+1)>1.e-10) then
[2147]418              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)
419           ENDif
420        END DO
421     END DO
422  END DO
423  DO k=1,klev
424     DO j=1,klev-1
425        DO i=1,klon
426           za(i,j,k)=max(0.,zmd(i,j+1,k)-zmd(i,j,k))
427        END DO
428     END DO
429  END DO
430!!!!! quantite  de traceur dans la descente d'air insaturee  :   4 juin 2012
431  DO k=1,klev
432     DO j=1,klev-1
433        DO i=1,klon
[5082]434        if(mp(i,j+1)>1.e-10) then
[2147]435          qTrdi(i,j+1,it)=qTrdi(i,j+1,it)+(zmd(i,j+1,k)/mp(i,j+1))*tr(i,k,it)
436        else
437          qTrdi(i,j,it)=0.!tr(i,j,it)
438        endif
439        ENDDO
440     ENDDO
441  ENDDO
442!!!!!
[5099]443
[2147]444! rajout du terme lie a l ascendance induite
[5099]445
[2147]446  DO j=2,klev
447     DO i=1,klon
448        za(i,j,j-1)=za(i,j,j-1)+mp(i,j)
449     END DO
450  END DO
[5099]451
[2147]452! tendance courants insatures  ! sans lessivage ancien schema
[5099]453
[2147]454  DO k=1,klev
455     DO j=1,klev
456        DO i=1,klon
457           zmfd(i,j,it)=zmfd(i,j,it)+za(i,j,k)*(tr(i,k,it)-tr(i,j,it))
458        END DO
459     END DO
460  END DO
[5099]461
[2147]462! =========================================
463! calcul des tendances liees aux courants satures   j <-> z ; k <-> z'
464! =========================================
465!RL
466!  Feeding concentrations
467  DO j=1,klev
468     DO i=1,klon
469        qfeed(i,it)=qfeed(i,it)+wght_cvfd(i,j)*tr(i,j,it)
470     END DO
471  END DO
472!RL
[5099]473
[2147]474  DO j=1,klev
475     DO i=1,klon
476!RL
477!!        zmfa(i,j,it)=da(i,j)*(tr(i,1,it)-tr(i,j,it))                     ! da
478        zmfa(i,j,it)=da(i,j)*(qfeed(i,it)-tr(i,j,it))                     ! da
479!RL
480     END DO
481  END DO
[5099]482
[2147]483  DO k=1,klev
484     DO j=1,klev
485        DO i=1,klon
486           zmfp(i,j,it)=zmfp(i,j,it)+phi(i,j,k)*(tr(i,k,it)-tr(i,j,it))  ! phi
487        END DO
488     END DO
489  END DO
490! RomP ajout des matrices liees au lessivage
491  DO j=1,klev
492     DO i=1,klon
493        zmfd1a(i,j,it)=d1a(i,j)*tr(i,1,it)   ! da1
494        zmfdam(i,j,it)=dam(i,j)*tr(i,1,it)   ! dam
495     END DO
496  END DO
497  DO k=1,klev
498     DO j=1,klev
499        DO i=1,klon
500          zmfphi2(i,j,it)=zmfphi2(i,j,it)+phi2(i,j,k)*tr(i,k,it)  ! psi
501        END DO
502     END DO
503  END DO
504  DO j=1,klev-1
505     DO i=1,klon
506        zmfu(i,j,it)=max(0.,upd(i,j+1)+dnd(i,j+1))*(tr(i,j+1,it)-tr(i,j,it))
507     END DO
508  END DO
509  DO j=2,klev
510     DO i=1,klon
511        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))
512     END DO
513  END DO
514! ===================================================
515! calcul des tendances liees aux courants insatures
516! ===================================================
517!  pression 
518  DO k=1, klev
519     DO i=1, klon
520        dxpres(i,k)=paprs(i,k)-paprs(i,k+1)
521     ENDDO
522  ENDDO
523  pdtimeRG=pdtime*RG
524
525! q_pa et q_pm    traceurs issues des courants satures se retrouvant dans les precipitations
526  DO j=1,klev
527     DO i=1,klon
[5082]528        if(j>=icb(i).and.j<=inb(i)) then
529          if(clw(i,j)>1.e-16) then
[2147]530!           qPa(i,j,it)=ccntrAA_coef*tr(i,1,it)/clw(i,j)
531           qPa(i,j,it)=ccntrAA_coef3d(i,j)*tr(i,1,it)/clw(i,j)
532          else
533           qPa(i,j,it)=0.
534          endif
535        endif
536     END DO
537  END DO
538
539! calcul de q_pm en 2 parties :
540! 1) calcul de sa valeur pour un niveau z' donne
541! 2) integration sur la verticale sur z'
542     DO j=1,klev
543      DO k=1,j-1
544        DO i=1,klon
[5082]545        if(k>=icb(i).and.k<=inb(i).and.&
546           j<=inb(i)) then
547            if(elij(i,k,j)>1.e-16) then
[2147]548!JE             qMeltmp(i,j,it)=((1-ep(i,k))*ccntrAA_coef*tr(i,1,it)&
549!JE                         *(1.-sij(i,k,j))  +ccntrENV_coef&
550!JE                         *tr(i,k,it)*sij(i,k,j)) / elij(i,k,j)
551             qMeltmp(i,j,it)=((1-ep(i,k))*ccntrAA_coef3d(i,k)*tr(i,1,it)&
552                         *(1.-sij(i,k,j))  +ccntrENV_coef3d(i,k)&
553                         *tr(i,k,it)*sij(i,k,j)) / elij(i,k,j)
554            else
555             qMeltmp(i,j,it)=0.
556            endif
557          qpmMint(i,j,it)=qpmMint(i,j,it) + qMeltmp(i,j,it)*epmlmMm(i,j,k)
558          Mint(i,j)=Mint(i,j) + epmlmMm(i,j,k)
559        endif ! end if dans nuage
560        END DO
561     END DO
562  END DO
563
564     DO j=1,klev
565        DO i=1,klon
[5082]566          if(Mint(i,j)>1.e-16) then
[2147]567           qMel(i,j,it)=qpmMint(i,j,it)/Mint(i,j)
568          else
569           qMel(i,j,it)=0.
570          endif
571     END DO
572  END DO
573
574! calcul de q_d et q_p    traceurs de la descente precipitante
575   DO j=klev-1,1,-1
576    DO i=1,klon
[5082]577     if(mp(i,j+1)>mp(i,j).and.mp(i,j+1)>1.e-10) then  ! detrainement
[2147]578       kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
579                (-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))&
580                + (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j)))
581             
[5082]582     elseif(mp(i,j)>mp(i,j+1).and.mp(i,j)>1.e-10) then! entrainement
583       if(j==1) then
[2147]584        kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
585                 (-mp(i,2)-imp(i,j)/RG*dxpres(i,j))&
586                 + (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j)))
587       else
588        kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
589                 (-mp(i,j)-imp(i,j)/RG*dxpres(i,j))&
590                 + (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j)))
591       endif
592      else
593        kappa(i,j)=1.
594      endif
595    ENDDO
596   ENDDO
597
598  DO j=klev-1,1,-1
599   DO i=1,klon
[5082]600    if (abs(kappa(i,j))<1.e-25) then    !si denominateur nul (il peut y avoir des mp!=0)
[2147]601     kappa(i,j)=1.
[5082]602     if(j==1) then
[2147]603       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
[5082]604     elseif(mp(i,j+1)>mp(i,j).and.mp(i,j+1)>1.e-10) then
[2147]605       qDi(i,j,it)=qDi(i,j+1,it)
[5082]606     elseif(mp(i,j)>mp(i,j+1).and.mp(i,j)>1.e-10) then! entrainement
[2147]607       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))
608     else  ! si  mp (i)=0 et mp(j+1)=0
609       qDi(i,j,it)=tr(i,j,it) ! orig 0.
610     endif
611
612      if(NO_precip(i,j)) then
613       qPr(i,j,it)=0.
614      else
615       qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
616                   Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)&
617                   +imp(i,j)/RG*dxpres(i,j)*qDi(i,j,it))/&
618                  (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
619      endif
620    else   !     denominateur non nul
621     kappa(i,j)=1./kappa(i,j)     
622! calcul de qd et qp
623!!jyg  (20130119) correction pour le sommet du nuage
624!!     if(j.ge.inb(i)) then       !au-dessus du nuage, sommet inclu
[5082]625     if(j>inb(i)) then       !au-dessus du nuage
[2147]626       qDi(i,j,it)=tr(i,j,it)   ! pas de descente => environnement = descente insaturee
627       qPr(i,j,it)=0.
628
629!  vvv premiere couche du modele ou mp(1)=0  ! det tout le temps  vvv
[5082]630     elseif(j==1) then
631      if(mp(i,2)>1.e-10) then !mp(2) non nul -> detrainement (car mp(1) = 0) !ent pas possible
[2147]632       if(NO_precip(i,j)) then !pas de precip en (i)
633        qDi(i,j,it)=qDi(i,j+1,it)
634        qPr(i,j,it)=0.
635       else
636        qDi(i,j,it)=kappa(i,j)*(&
637            (-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
638            Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +&
639            (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
640            (-mp(i,j+1)*qDi(i,j+1,it)))
641
642        qPr(i,j,it)=kappa(i,j)*(&
643            (-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))*&
644            ((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
645            Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))&
646            +(-mp(i,j+1)*qDi(i,j+1,it)) * (imp(i,j)/RG*dxpres(i,j)))
647       endif
648
649      else !mp(2) nul -> plus de descente insaturee -> pluie agit sur environnement
650        qDi(i,j,it)=tr(i,j,it) ! orig 0.
651        if(NO_precip(i,j)) then
652         qPr(i,j,it)=0.
653        else
654         qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
655                   Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)&
656                   +imp(i,j)/RG*dxpres(i,j)*tr(i,j,it))/&
657                  (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
658        endif
659
660      endif  !mp(2) nul ou non
661
662! vvv  (j!=1.and.j.lt.inb(i))  en-dessous du sommet nuage   vvv
663     else
664!------------------------------------------------------------- detrainement
[5082]665      if(mp(i,j+1)>mp(i,j).and.mp(i,j+1)>1.e-10) then !mp(i,j).gt.1.e-10) then
[2147]666         if(NO_precip(i,j)) then
667          qDi(i,j,it)=qDi(i,j+1,it)
668          qPr(i,j,it)=0.
669         else
670          qDi(i,j,it)=kappa(i,j)*(&
671            (-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
672            Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +&
673            (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
674            (-mp(i,j+1)*qDi(i,j+1,it)))
[5099]675
[2147]676          qPr(i,j,it)=kappa(i,j)*(&
677            (-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))*&
678            ((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
679            Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))&
680            +(-mp(i,j+1)*qDi(i,j+1,it)) * (imp(i,j)/RG*dxpres(i,j)))
681         endif !precip
682!------------------------------------------------------------- entrainement
[5082]683      elseif(mp(i,j)>mp(i,j+1).and.mp(i,j)>1.e-10) then
[2147]684         if(NO_precip(i,j)) then
685          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))
686          qPr(i,j,it)=0.
687         else
688          qDi(i,j,it)=kappa(i,j)*(&
689            (-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
690            Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +&
691            (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
692            (-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it)))
[5099]693
[2147]694          qPr(i,j,it)=kappa(i,j)*(&
695            (-mp(i,j)-imp(i,j)/RG*dxpres(i,j))*&
696            ((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
697            Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))&
698            +(-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it))*&
699            (imp(i,j)/RG*dxpres(i,j)))
700         endif !precip
701!------------------------------------------------------------- endif ! ent/det
702      else  !mp nul
703        qDi(i,j,it)=tr(i,j,it) ! orig 0.
704        if(NO_precip(i,j)) then
705         qPr(i,j,it)=0.
706        else
707         qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
708                   Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)&
709                   +imp(i,j)/RG*dxpres(i,j)*tr(i,j,it))/&
710                  (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
711        endif
712      endif ! mp nul ou non
713     endif ! condition sur j
714    endif ! kappa
715   ENDDO
716  ENDDO
717
718!! print test descente insaturee
719!  DO j=klev,1,-1
720!   DO i=1,klon
721!     if(it.eq.3) then
722!   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,&
723!!         'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),&
724!          'zmfdam+zmfpsi',zmfdam(i,j,it)+zmfphi2(i,j,it),'qpmMint',qpmMint(i,j,it),&
725!          'Pm',Pm(i,j),'Mint',Mint(i,j),&
726!!          'zmfa',zmfa(i,j,it),'zmfp',zmfp(i,j,it),&
727!        'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),'zmfd1a',zmfd1a(i,j,it)
728!!          'Pa',Pa(i,j),'eplaMm',eplaMm(i,j)
729!!         'zmfd1a=da1*qa',zmfd1a(i,j,it),'Pa*qPa',wdtrainA(i,j)*qPa(i,j,it),'da1',d1a(i,j)
730!     endif
731!   ENDDO
732!  ENDDO
733
734
735! ===================================================
736! calcul final des tendances
737! ===================================================
738
739  DO k=klev-1,1,-1
740    DO i=1, klon
741! transport
742     tdcvMA=zmfd(i,k,it)+zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it)   ! double comptage des downdraft insatures
743     trsptrac=zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it)
744! lessivage courants satures
745!JE     scavtrac=-ccntrAA_coef*zmfd1a(i,k,it)&
746!JE               -zmfphi2(i,k,it)*ccntrENV_coef&
747!JE               -zmfdam(i,k,it)*ccntrAA_coef
748     scavtrac=-ccntrAA_coef3d(i,k)*zmfd1a(i,k,it)&
749               -zmfphi2(i,k,it)*ccntrENV_coef3d(i,k)&
750               -zmfdam(i,k,it)*ccntrAA_coef3d(i,k)
751! lessivage courants insatures
[5082]752   if(k<=inb(i).and.k>1) then   ! tendances dans le nuage
[2147]753!------------------------------------------------------------- detrainement
[5082]754      if(mp(i,k+1)>mp(i,k).and.mp(i,k+1)>1.e-10) then
[2147]755        uscavtrac= (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))&
756                   + mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
[5099]757
[2147]758!        if(it.eq.3) write(*,'(I2,1X,a,5X,e20.12,82X,a,e20.12)')k,' det incloud',&
759!                    (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+&
760!                    mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k),&
761!                    'mp',mp(i,k)
762!------------------------------------------------------------- entrainement
[5082]763      elseif(mp(i,k)>mp(i,k+1).and.mp(i,k)>1.e-10) then
[2147]764         uscavtrac= mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
[5099]765
[2147]766!         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)
767!=!------------------------------------------------------------- end ent/det
768      else !        mp(i,k+1)=0. et mp(i,k)=0.    pluie directement sur l environnement
769
770        if(NO_precip(i,k)) then
771          uscavtrac=0.
772!         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)
773        else
774          uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/RG+evap(i,k)*qPr(i,k,it)*dxpres(i,k)/RG
775!         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)
776!!JE adds
777!         if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'imp',imp(i,k)
778!         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)
779!         if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'evap',evap(i,k)
780!         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)
781!         if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'dxpres',dxpres(i,k)
782!!Je end
783
784        endif
785      endif  ! mp/det/ent
786!------------------------------------------------------------- premiere couche
[5082]787   elseif(k==1) then  !                                      mp(1)=0.
788      if(mp(i,2)>1.e-10) then  !detrainement
[2147]789         uscavtrac= (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it)) !&
790!                   + mp(i,2)*(0.-tr(i,k,it))
[5099]791
[2147]792!       if(it.eq.3) write(*,'(I2,1X,a,e20.12,84X,a,e20.12)')k,' 1 det',&
793!                  (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+&
794!                  mp(i,2)*(0.-tr(i,k,it))*pdtimeRG/dxpres(i,k),&
795!                   'mp',mp(i,k)
796      else   ! mp(2) = 0 = mp(1) pas de descente insaturee, rien ne se passe s'il ne pleut pas, sinon pluie->env
797        if(NO_precip(i,1)) then
798          uscavtrac=0.
799        else
800          uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/RG+evap(i,k)*qPr(i,k,it)*dxpres(i,k)/RG
801        endif
802!         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)
803      endif
804
805   else   ! k > INB  au-dessus du nuage
806    uscavtrac=0.
807   endif
808
809! =====    tendances finales  ======
810     trsptd(i,k,it)=trsptrac*pdtimeRG/dxpres(i,k)              ! td transport sans eau dans courants satures
811     dtrSscav(i,k,it)=scavtrac*pdtimeRG/dxpres(i,k)            ! td du lessivage dans courants satures
812     dtrUscav(i,k,it)=uscavtrac*pdtimeRG/dxpres(i,k)           ! td courant insat
813     dtrsat(i,k,it)=(trsptrac+scavtrac)*pdtimeRG/dxpres(i,k)   ! td courant sat
814     dtrcv(i,k,it)=(trsptrac+scavtrac+uscavtrac)*pdtimeRG/dxpres(i,k)!dtrsat(i,k,it)+dtrUscav(i,k,it)         td conv
815!!!!!!
816     dtrcvMA(i,k,it)=tdcvMA*pdtimeRG/dxpres(i,k) ! MA tendance convection
817     ENDDO
818  ENDDO
819
820! test de conservation du traceur
[5103]821!PRINT*,"_____________________________________________________________"
822!PRINT*,"                                                             "
[2147]823!      conserv=0.
824!      conservMA=0.
825!      DO k= klev-1,1,-1
826!        DO i=1, klon
827!         conserv=conserv+dtrcv(i,k,it)*   &
828!        (paprs(i,k)-paprs(i,k+1))/RG
829!         conservMA=conservMA+dtrcvMA(i,k,it)*   &
830!        (paprs(i,k)-paprs(i,k+1))/RG
[5099]831
[2147]832!      if(it.eq.3)  write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)') k,&
833!         'MA td ',dtrcvMA(i,k,it)*dxpres(i,k)/RG,&
834!         ' td',dtrcv(i,k,it)*dxpres(i,k)/RG,'         conservMA ',conservMA,'conserv ',conserv
835!!
836!        ENDDO
837!      ENDDO
838!       if(it.eq.3) print *,'it',it,'conserv ',conserv,'conservMA ',conservMA
839
840END SUBROUTINE cvltr_spl
Note: See TracBrowser for help on using the repository browser.