LMDZ
cvltr_spl.F90
Go to the documentation of this file.
1 !
2 ! $Id $
3 !
4 SUBROUTINE 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
15  USE infotrac_phy, ONLY : nbtr,tname
16  IMPLICIT NONE
17 !=====================================================================
18 ! Objet : convection des traceurs / KE
19 ! Auteurs: M-A Filiberti and J-Y Grandpeix
20 ! modifiee par R Pilon : lessivage des traceurs / KE
21 !=====================================================================
22 
23  include "YOMCST.h"
24  include "YOECUMF.h"
25  include "conema3.h"
26  include "chem.h"
27 
28 ! Entree
29  REAL,INTENT(IN) :: pdtime
30  REAL,DIMENSION(klon,klev),INTENT(IN) :: da
31  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi
32 ! RomP
33  REAL,DIMENSION(klon,klev),INTENT(IN) :: d1a,dam ! matrices pour simplifier
34  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi2 ! l'ecriture des tendances
35 !
36  REAL,DIMENSION(klon,klev),INTENT(IN) :: mpIN
37  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression aux 1/2 couches (bas en haut)
38 ! 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 
342 IF (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')
357 ELSE
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 
367 ENDIF
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 
841 END SUBROUTINE cvltr_spl
infotrac_phy::nbtr
integer, save nbtr
Definition: infotrac_phy.F90:20
r
!$Id mode_top_bound COMMON comconstr r
Definition: comconst.h:7
dimphy::klon
integer, save klon
Definition: dimphy.F90:3
infotrac_phy
Definition: infotrac_phy.F90:4
rho_water
!$Id!INTEGER ih2o2 REAL rho_water
Definition: chem.h:4
cvltr_spl
subroutine cvltr_spl(pdtime, da, phi, phi2, d1a, dam, mpIN, epIN, sigd, sij, wght_cvfd, clw, elij, epmlmMm, eplaMm, pmflxrIN, pmflxsIN, ev, te, wdtrainA, wdtrainM, paprs, it, tr, upd, dnd, inb, icb, kk, henry, zrho, ccntrAA_spla, ccntrENV_spla, coefcoli_spla, id_prec, id_fine, id_coss, id_codu, id_scdu, dtrcv, trsptd, dtrSscav, dtrsat, dtrUscav, qDi, qPr, qPa, qMel, qTrdi, dtrcvMA, Mint, zmfd1a, zmfphi2, zmfdam)
Definition: cvltr_spl.F90:13
dimphy::klev
integer, save klev
Definition: dimphy.F90:7
it
!$Id Turb_fcg_gcssold get_uvd it
Definition: 1D_interp_cases.h:16
paprs
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL &zphi geo500!IM on interpole a chaque pas de temps le paprs
Definition: calcul_STDlev.h:144
infotrac_phy::tname
character(len=20), dimension(:), allocatable, save tname
Definition: infotrac_phy.F90:27
epmax
!$Id sig2feed!common comconema2 iflag_cvl_sigd common comconema1 epmax
Definition: conema3.h:15
dimphy
Definition: dimphy.F90:1
sigd
!$Id!Parameters for nlm real sigd
Definition: cv30param.h:5
rg
real rg
Definition: comcstphy.h:1