source: LMDZ6/trunk/libf/phylmd/cvltr_spl.f90 @ 5367

Last change on this file since 5367 was 5292, checked in by abarral, 7 weeks ago

Move academic.h chem.h chem_spla.h to module

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