My Project
 All Classes Files Functions Variables Macros
phytrac.F90
Go to the documentation of this file.
1 !$Id $
2 
3 SUBROUTINE phytrac( &
4  nstep, julien, gmtime, debutphy, &
5  lafin, pdtphys, u, v, t_seri, &
6  paprs, pplay, pmfu, pmfd, &
7  pen_u, pde_u, pen_d, pde_d, &
8  cdragh, coefh, fm_therm, entr_therm,&
9  yu1, yv1, ftsol, pctsrf, &
10  ustar, u10m, v10m, &
11  xlat, xlon, &
12  frac_impa,frac_nucl,beta_fisrt,beta_v1, &
13  presnivs, pphis, pphi, albsol, &
14  sh, rh, cldfra, rneb, &
15  diafra, cldliq, itop_con, ibas_con, &
16  pmflxr, pmflxs, prfl, psfl, &
17  da, phi, mp, upwd, &
18  phi2, d1a, dam, sij, & ! RomP
19  wdtraina, wdtrainm, sigd, clw,elij, & ! RomP
20  evap, ep, epmlmmm, eplamm, & ! RomP
21  dnwd, aerosol_couple, flxmass_w, &
22  tau_aero, piz_aero, cg_aero, ccm, &
23  rfname, &
24  d_tr_dyn, & ! RomP
25  tr_seri)
26 !
27 !======================================================================
28 ! Auteur(s) FH
29 ! Objet: Moniteur general des tendances traceurs
30 ! Modification R. Pilon 01 janvier 2012 transport+scavenging KE scheme : cvltr
31 ! Modification R. Pilon 10 octobre 2012 large scale scavenging incloud_scav + bc_scav
32 !======================================================================
33 
34  USE ioipsl
35  USE phys_cal_mod, only : hour
36  USE phys_output_mod, only : convers_timesteps
37  USE dimphy
38  USE infotrac
39  USE mod_grid_phy_lmdz
40  USE mod_phys_lmdz_para
41  USE comgeomphy
42  USE iophy
43  USE traclmdz_mod
44  USE tracinca_mod
45  USE tracreprobus_mod
46  USE control_mod
47 
48  IMPLICIT NONE
49 
50  include "YOMCST.h"
51  include "dimensions.h"
52  include "indicesol.h"
53  include "clesphys.h"
54  include "temps.h"
55  include "paramet.h"
56  include "thermcell.h"
57  include "iniprint.h"
58 !==========================================================================
59 ! -- ARGUMENT DESCRIPTION --
60 !==========================================================================
61 
62 ! Input arguments
63 !----------------
64 !Configuration grille,temps:
65  INTEGER,INTENT(IN) :: nstep ! Appel physique
66  INTEGER,INTENT(IN) :: julien ! Jour julien
67  REAL,INTENT(IN) :: gmtime ! Heure courante
68  REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde)
69  LOGICAL,INTENT(IN) :: debutphy ! le flag de l'initialisation de la physique
70  LOGICAL,INTENT(IN) :: lafin ! le flag de la fin de la physique
71 
72  REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point
73  REAL,DIMENSION(klon),INTENT(IN) :: xlon ! longitudes pour chaque point
74 !
75 !Physique:
76 !--------
77  REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
78  REAL,DIMENSION(klon,klev),INTENT(IN) :: u ! variable not used
79  REAL,DIMENSION(klon,klev),INTENT(IN) :: v ! variable not used
80  REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique
81  REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! humidite relative
82  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa)
83  REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
84  REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi ! geopotentiel
85  REAL,DIMENSION(klon),INTENT(IN) :: pphis
86  REAL,DIMENSION(klev),INTENT(IN) :: presnivs
87  REAL,DIMENSION(klon,klev),INTENT(IN) :: cldliq ! eau liquide nuageuse
88  REAL,DIMENSION(klon,klev),INTENT(IN) :: cldfra ! fraction nuageuse (tous les nuages)
89  REAL,DIMENSION(klon,klev),INTENT(IN) :: diafra ! fraction nuageuse (convection ou stratus artificiels)
90  REAL,DIMENSION(klon,klev),INTENT(IN) :: rneb ! fraction nuageuse (grande echelle)
91 !
92  REAL :: ql_incl ! contenu en eau liquide nuageuse dans le nuage ! ql_incl=oliq/rneb
93  REAL,DIMENSION(klon,klev),INTENT(IN) :: beta_fisrt ! taux de conversion de l'eau cond (de fisrtilp)
94  REAL,DIMENSION(klon,klev),INTENT(out) :: beta_v1 ! -- (originale version)
95 
96 !
97  INTEGER,DIMENSION(klon),INTENT(IN) :: itop_con
98  INTEGER,DIMENSION(klon),INTENT(IN) :: ibas_con
99  REAL,DIMENSION(klon),INTENT(IN) :: albsol ! albedo surface
100 !
101 !Dynamique
102 !--------
103  REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: d_tr_dyn
104 !
105 !Convection:
106 !----------
107  REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfu ! flux de masse dans le panache montant
108  REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfd ! flux de masse dans le panache descendant
109  REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_u ! flux entraine dans le panache montant
110  REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_u ! flux detraine dans le panache montant
111  REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_d ! flux entraine dans le panache descendant
112  REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_d ! flux detraine dans le panache descendant
113 
114 !...Tiedke
115  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: pmflxr, pmflxs ! Flux precipitant de pluie, neige aux interfaces [convection]
116  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: prfl, psfl ! Flux precipitant de pluie, neige aux interfaces [large-scale]
117 
118  LOGICAL,INTENT(IN) :: aerosol_couple
119  REAL,DIMENSION(klon,klev),INTENT(IN) :: flxmass_w
120  REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: tau_aero
121  REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: piz_aero
122  REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: cg_aero
123  CHARACTER(len=4),DIMENSION(9),INTENT(IN) :: rfname
124  REAL,DIMENSION(klon,klev,2),INTENT(IN) :: ccm
125 !... K.Emanuel
126  REAL,DIMENSION(klon,klev),INTENT(IN) :: da
127  REAL,DIMENSION(klon,klev,klev),INTENT(IN):: phi
128 ! RomP >>>
129  REAL,DIMENSION(klon,klev),INTENT(IN) :: d1a,dam
130  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi2
131 !
132  REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtraina
133  REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainm
134  REAL,DIMENSION(klon),INTENT(IN) :: sigd
135 ! ---- RomP flux entraine, detraine et precipitant kerry Emanuel
136  REAL,DIMENSION(klon,klev),INTENT(IN) :: evap
137  REAL,DIMENSION(klon,klev),INTENT(IN) :: ep
138  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: sij
139  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: elij
140  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: epmlmmm
141  REAL,DIMENSION(klon,klev),INTENT(IN) :: eplamm
142  REAL,DIMENSION(klon,klev),INTENT(IN) :: clw
143 ! RomP <<<
144 
145 !
146  REAL,DIMENSION(klon,klev),INTENT(IN) :: mp
147  REAL,DIMENSION(klon,klev),INTENT(IN) :: upwd ! saturated updraft mass flux
148  REAL,DIMENSION(klon,klev),INTENT(IN) :: dnwd ! saturated downdraft mass flux
149 !
150 !Thermiques:
151 !----------
152  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: fm_therm
153  REAL,DIMENSION(klon,klev),INTENT(IN) :: entr_therm
154 !
155 !Couche limite:
156 !--------------
157 !
158  REAL,DIMENSION(klon),INTENT(IN) :: cdragh ! coeff drag pour T et Q
159  REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh ! coeff melange CL (m**2/s)
160  REAL,DIMENSION(klon),INTENT(IN) :: ustar,u10m,v10m ! u* & vent a 10m (m/s)
161  REAL,DIMENSION(klon),INTENT(IN) :: yu1 ! vents au premier niveau
162  REAL,DIMENSION(klon),INTENT(IN) :: yv1 ! vents au premier niveau
163 !
164 !Lessivage:
165 !----------
166 !
167 ! pour le ON-LINE
168 !
169  REAL,DIMENSION(klon,klev),INTENT(IN) :: frac_impa ! fraction d'aerosols non impactes
170  REAL,DIMENSION(klon,klev),INTENT(IN) :: frac_nucl ! fraction d'aerosols non nuclees
171 
172 ! Arguments necessaires pour les sources et puits de traceur:
173  REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol ! Temperature du sol (surf)(Kelvin)
174  REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol (nature du sol)
175 
176 
177 ! Output argument
178 !----------------
179  REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
180  REAL,DIMENSION(klon,klev) :: sourcebe
181 !=======================================================================================
182 ! -- LOCAL VARIABLES --
183 !=======================================================================================
184 
185  INTEGER :: i, k, it
186  INTEGER :: nsplit
187 
188 !Sources et Reservoirs de traceurs (ex:Radon):
189 !--------------------------------------------
190 !
191  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: source ! a voir lorsque le flux de surface est prescrit
192 !$OMP THREADPRIVATE(source)
193 
194 !
195 !Entrees/Sorties: (cf ini_histrac.h et write_histrac.h)
196 !---------------
197  INTEGER :: iiq, ierr
198  INTEGER :: nhori, nvert
199  REAL :: zsto, zout, zjulian
200  INTEGER,SAVE :: nid_tra ! pointe vers le fichier histrac.nc
201 !$OMP THREADPRIVATE(nid_tra)
202  REAL,DIMENSION(klon) :: zx_tmp_fi2d ! variable temporaire grille physique
203  INTEGER :: itau_w ! pas de temps ecriture = nstep + itau_phy
204  LOGICAL,PARAMETER :: ok_sync=.true.
205  CHARACTER(len=20),save :: chtratimestep,chtratimestep_omp
206 !$OMP THREADPRIVATE(chtratimestep)
207 !
208 ! Nature du traceur
209 !------------------
210  LOGICAL,DIMENSION(:),ALLOCATABLE,SAVE :: aerosol ! aerosol(it) = true => aerosol => lessivage
211 !$OMP THREADPRIVATE(aerosol) ! aerosol(it) = false => gaz
212  REAL,DIMENSION(klon,klev) :: delp ! epaisseur de couche (Pa)
213 !
214 ! Tendances de traceurs (Td) et flux de traceurs:
215 !------------------------
216  REAL,DIMENSION(klon,klev) :: d_tr ! Td dans l'atmosphere
217  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_cl ! Td couche limite/traceur
218  REAL,DIMENSION(klon,nbtr) :: d_tr_dry ! Td depot sec/traceur (1st layer) jyg
219  REAL,DIMENSION(klon,nbtr) :: flux_tr_dry ! depot sec/traceur (surface) jyg
220  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec !RomP
221  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_cv ! Td convection/traceur
222 ! RomP >>>
223  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_insc
224  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_bcscav
225  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_evapls
226  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_ls
227  REAL,DIMENSION(klon,nbtr) :: qprls !jyg: concentration tra dans pluie LS a la surf.
228  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_trsp
229  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_sscav
230  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_sat
231  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_uscav
232  REAL,DIMENSION(klon,klev,nbtr) :: qpr,qdi ! concentration tra dans pluie,air descente insaturee
233  REAL,DIMENSION(klon,klev,nbtr) :: qpa,qmel
234  REAL,DIMENSION(klon,klev,nbtr) :: qtrdi,dtrcvma ! conc traceur descente air insaturee et td convective MA
235  REAL,DIMENSION(klon,klev) :: mint
236  REAL,DIMENSION(klon,klev,nbtr) :: zmfd1a
237  REAL,DIMENSION(klon,klev,nbtr) :: zmfdam
238  REAL,DIMENSION(klon,klev,nbtr) :: zmfphi2
239 ! RomP <<<
240  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_th ! Td thermique
241  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_lessi_impa ! Td du lessivage par impaction
242  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_lessi_nucl ! Td du lessivage par nucleation
243 !
244 ! Physique
245 !----------
246  REAL,DIMENSION(klon,klev,nbtr) :: flestottr ! flux de lessivage dans chaque couche
247  REAL,DIMENSION(klon,klev) :: zmasse ! densité atmosphérique Kg/m2
248  REAL,DIMENSION(klon,klev) :: ztra_th
249 !PhH
250  REAL,DIMENSION(klon,klev) :: zrho
251  REAL,DIMENSION(klon,klev) :: zdz
252  REAL :: evaplsc,dx,beta ! variable pour lessivage Genthon
253  REAL,DIMENSION(klon) :: his_dh ! ---
254 ! in-cloud scav variables
255  REAL :: ql_incloud_ref ! ref value of in-cloud condensed water content
256 
257 !Controles:
258 !---------
259  LOGICAL,SAVE :: couchelimite=.true.
260  LOGICAL,SAVE :: convection=.true.
261  LOGICAL,SAVE :: lessivage
262 !$OMP THREADPRIVATE(couchelimite,convection,lessivage)
263 
264  CHARACTER(len=8),DIMENSION(nbtr) :: solsym
265 !RomP >>>
266  INTEGER,SAVE :: iflag_lscav_omp,iflag_lscav
267  LOGICAL,SAVE :: convscav_omp,convscav
268 !$OMP THREADPRIVATE(iflag_lscav,convscav)
269 !RomP <<<
270 !######################################################################
271 ! -- INITIALIZATION --
272 !######################################################################
273  DO k=1,klev
274  DO i=1,klon
275  sourcebe(i,k)=0.
276  mint(i,k)=0.
277  zrho(i,k)=0.
278  zdz(i,k)=0.
279  END DO
280  END DO
281 
282  DO it=1, nbtr
283  DO k=1,klev
284  DO i=1,klon
285  d_tr_insc(i,k,it)=0.
286  d_tr_bcscav(i,k,it)=0.
287  d_tr_evapls(i,k,it)=0.
288  d_tr_ls(i,k,it)=0.
289  d_tr_cv(i,k,it)=0.
290  d_tr_cl(i,k,it)=0.
291  d_tr_trsp(i,k,it)=0.
292  d_tr_sscav(i,k,it)=0.
293  d_tr_sat(i,k,it)=0.
294  d_tr_uscav(i,k,it)=0.
295  d_tr_lessi_impa(i,k,it)=0.
296  d_tr_lessi_nucl(i,k,it)=0.
297  qdi(i,k,it)=0.
298  qpr(i,k,it)=0.
299  qpa(i,k,it)=0.
300  qmel(i,k,it)=0.
301  qtrdi(i,k,it)=0.
302  dtrcvma(i,k,it)=0.
303  zmfd1a(i,k,it)=0.
304  zmfdam(i,k,it)=0.
305  zmfphi2(i,k,it)=0.
306  END DO
307  END DO
308  END DO
309  IF (debutphy) THEN
310 !!jyg
311 !$OMP MASTER
312  chtratimestep_omp='DefFreq'
313  CALL getin('tra_time_step',chtratimestep_omp)
314 !$OMP END MASTER
315 !$OMP BARRIER
316  chtratimestep=chtratimestep_omp
317  IF (chtratimestep .NE. 'DefFreq') THEN
318  call convers_timesteps(chtratimestep,pdtphys,ecrit_tra)
319  ENDIF
320 !RomP >>>
321 !
322 !Config Key = convscav
323 !Config Desc = Convective scavenging switch: 0=off, 1=on.
324 !Config Def = .false.
325 !Config Help =
326 !
327 !$OMP MASTER
328  convscav_omp=.false.
329  call getin('convscav', convscav_omp)
330 !$OMP END MASTER
331 !$OMP BARRIER
332  convscav=convscav_omp
333  print*,'phytrac passage dans routine conv avec lessivage', convscav
334 !
335 !Config Key = iflag_lscav
336 !Config Desc = Large scale scavenging parametrization: 0=none, 1=old(Genthon92),
337 ! 2=1+PHeinrich, 3=Reddy_Boucher2004, 4=3+RPilon.
338 !Config Def = 1
339 !Config Help =
340 !
341 !$OMP MASTER
342  iflag_lscav_omp=1
343  call getin('iflag_lscav', iflag_lscav_omp)
344 !$OMP END MASTER
345 !$OMP BARRIER
346  iflag_lscav=iflag_lscav_omp
347 !
348  SELECT CASE(iflag_lscav)
349  CASE(0)
350  print*, 'Large scale scavenging: none'
351  CASE(1)
352  print*, 'Large scale scavenging: C. Genthon, Tellus(1992), 44B, 371-389'
353  CASE(2)
354  print*, 'Large scale scavenging: C. Genthon, modified P. Heinrich'
355  CASE(3)
356  print*, 'Large scale scavenging: M. Shekkar Reddy and O. Boucher, JGR(2004), 109, D14202'
357  CASE(4)
358  print*, 'Large scale scavenging: Reddy and Boucher, modified R. Pilon'
359  END SELECT
360 !RomP <<<
361  WRITE(*,*) 'FIRST TIME IN PHYTRAC : pdtphys(sec) = ',pdtphys,'ecrit_tra (sec) = ',ecrit_tra
362  ALLOCATE( source(klon,nbtr), stat=ierr)
363  IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 1',1)
364 
365  ALLOCATE( aerosol(nbtr), stat=ierr)
366  IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 2',1)
367 
368 
369  ! Initialize module for specific tracers
370  SELECT CASE(type_trac)
371  CASE('lmdz')
372  CALL traclmdz_init(pctsrf, xlat, xlon, ftsol, tr_seri, t_seri, pplay, sh, pdtphys, aerosol, lessivage)
373  CASE('inca')
374  source(:,:)=0.
375  CALL tracinca_init(aerosol,lessivage)
376  CASE('repr')
377  source(:,:)=0.
378  END SELECT
379 !
380 ! Initialize diagnostic output
381 ! ----------------------------
382 #ifdef CPP_IOIPSL
383  include "ini_histrac.h"
384 #endif
385  END IF ! of IF (debutphy)
386 !############################################ END INITIALIZATION #######
387 
388  DO k=1,klev
389  DO i=1,klon
390  zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
391  END DO
392  END DO
393 !
394  IF (id_be .GT. 0) THEN
395  DO k=1,klev
396  DO i=1,klon
397  sourcebe(i,k)=srcbe(i,k) !RomP -> pour sortie histrac
398  END DO
399  END DO
400  ENDIF
401 
402 !===============================================================================
403 ! -- Do specific treatment according to chemestry model or local LMDZ tracers
404 !
405 !===============================================================================
406  SELECT CASE(type_trac)
407  CASE('lmdz')
408  ! -- Traitement des traceurs avec traclmdz
409  CALL traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
410  cdragh, coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon,couchelimite,sh, &
411  rh, pphi, ustar, u10m, v10m, &
412 !! tr_seri, source, solsym, d_tr_cl, zmasse) !RomP
413  tr_seri, source, solsym, d_tr_cl,d_tr_dec, zmasse) !RomP
414  CASE('inca')
415  ! -- CHIMIE INCA config_inca = aero or chem --
416 
417  CALL tracinca(&
418  nstep, julien, gmtime, lafin, &
419  pdtphys, t_seri, paprs, pplay, &
420  pmfu, ftsol, pctsrf, pphis, &
421  pphi, albsol, sh, rh, &
422  cldfra, rneb, diafra, cldliq, &
423  itop_con, ibas_con, pmflxr, pmflxs, &
424  prfl, psfl, aerosol_couple, flxmass_w, &
425  tau_aero, piz_aero, cg_aero, ccm, &
426  rfname, &
427  tr_seri, source, solsym)
428 
429  CASE('repr')
430  ! -- CHIMIE REPROBUS --
431 
432  CALL tracreprobus(pdtphys, gmtime, debutphy, julien, &
433  presnivs, xlat, xlon, pphis, pphi, &
434  t_seri, pplay, paprs, sh , &
435  tr_seri, solsym)
436 
437  END SELECT
438 !======================================================================
439 ! -- Calcul de l'effet de la convection --
440 !======================================================================
441 
442  IF (convection) THEN
443  DO it=1, nbtr
444  IF ( conv_flg(it) == 0 ) cycle
445  IF (iflag_con.LT.2) THEN
446  d_tr_cv(:,:,it)=0.
447  ELSE IF (iflag_con.EQ.2) THEN
448 !..Tiedke
449  CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
450  pplay, paprs, tr_seri(:,:,it), d_tr_cv(:,:,it))
451 ! RomP >>>
452  ELSE
453 !..K.Emanuel !RomP modif arg
454  if (convscav.and.aerosol(it)) then ! lessivage convectif pour aerosol
455 !
456  CALL cvltr(pdtphys, da, phi,phi2,d1a,dam, mp,ep, &
457  sigd,sij,clw,elij,epmlmmm,eplamm, &
458  pmflxr,pmflxs,evap,t_seri,wdtraina,wdtrainm, &
459  paprs,it,tr_seri,upwd,dnwd,itop_con,ibas_con, &
460  d_tr_cv,d_tr_trsp,d_tr_sscav,d_tr_sat,d_tr_uscav,qdi,qpr,&
461  qpa,qmel,qtrdi,dtrcvma,mint, &
462  zmfd1a,zmfphi2,zmfdam)
463  else !pas de lessivage convectif ou n'est pas un aerosol
464  CALL cvltrorig(it,pdtphys, da, phi,mp,paprs,pplay,tr_seri,&
465  upwd,dnwd,d_tr_cv)
466  endif
467  END IF
468 ! RomP <<<
469 
470  DO k = 1, klev
471  DO i = 1, klon
472  tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
473  END DO
474  END DO
475 
476  CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'convection it = '//solsym(it))
477 
478  END DO ! nbtr
479  END IF ! convection
480 
481 !======================================================================
482 ! -- Calcul de l'effet des thermiques --
483 !======================================================================
484 
485  DO it=1,nbtr
486  DO k=1,klev
487  DO i=1,klon
488  d_tr_th(i,k,it)=0.
489  tr_seri(i,k,it)=max(tr_seri(i,k,it),0.)
490  tr_seri(i,k,it)=min(tr_seri(i,k,it),1.e10)
491  END DO
492  END DO
493  END DO
494 
495  IF (iflag_thermals.GT.0) THEN
496  nsplit=10
497  DO it=1, nbtr
498  DO isplit=1,nsplit
499 
500  CALL dqthermcell(klon,klev,pdtphys/nsplit, &
501  fm_therm,entr_therm,zmasse, &
502  tr_seri(1:klon,1:klev,it),d_tr,ztra_th)
503 
504  DO k=1,klev
505  DO i=1,klon
506  d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit
507  d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k)
508  tr_seri(i,k,it)=max(tr_seri(i,k,it)+d_tr(i,k),0.)
509  END DO
510  END DO
511  END DO ! nsplit
512  END DO ! it
513  END IF ! Thermiques
514 
515 !======================================================================
516 ! -- Calcul de l'effet de la couche limite --
517 !======================================================================
518 
519  IF (couchelimite) THEN
520 
521  DO k = 1, klev
522  DO i = 1, klon
523  delp(i,k) = paprs(i,k)-paprs(i,k+1)
524  END DO
525  END DO
526 
527  DO it=1, nbtr
528 
529  IF( pbl_flg(it) /= 0 ) THEN
530 
531  CALL cltrac(pdtphys, coefh,t_seri, &
532  tr_seri(:,:,it), source(:,it), &
533  paprs, pplay, delp, &
534  d_tr_cl(:,:,it),d_tr_dry(:,it),flux_tr_dry(:,it))
535 
536  DO k = 1, klev
537  DO i = 1, klon
538  tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
539  END DO
540  END DO
541  END IF
542 
543  END DO
544 
545  END IF ! couche limite
546 
547 
548 !======================================================================
549 ! Calcul de l'effet de la precipitation grande echelle
550 !======================================================================
551  IF (lessivage) THEN
552 
553  ql_incloud_ref = 10.e-4
554  ql_incloud_ref = 5.e-4
555 
556 
557 ! calcul du contenu en eau liquide au sein du nuage
558  ql_incl = ql_incloud_ref
559 ! choix du lessivage
560 !
561  IF (iflag_lscav .EQ. 3 .OR. iflag_lscav .EQ. 4) THEN
562 ! ******** Olivier Boucher version (3) possibly with modified ql_incl (4)
563 !
564  DO it = 1, nbtr
565 ! incloud scavenging and removal by large scale rain ! orig : ql_incl was replaced by 0.5e-3 kg/kg
566 ! the value 0.5e-3 kg/kg is from Giorgi and Chameides (1986), JGR
567 ! Liu (2001) proposed to use 1.5e-3 kg/kg
568 
569  CALL lsc_scav(pdtphys,it,iflag_lscav,ql_incl,prfl,psfl,rneb,beta_fisrt, &
570  beta_v1,pplay,paprs,t_seri,tr_seri,d_tr_insc, &
571  d_tr_bcscav,d_tr_evapls,qprls)
572 
573 !large scale scavenging tendency
574  DO k = 1, klev
575  DO i = 1, klon
576  d_tr_ls(i,k,it)=d_tr_insc(i,k,it)+d_tr_bcscav(i,k,it)+d_tr_evapls(i,k,it)
577  tr_seri(i,k,it)=tr_seri(i,k,it)+d_tr_ls(i,k,it)
578  ENDDO
579  ENDDO
580  CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'lsc scav it = '//solsym(it))
581  END DO !tr
582 
583  ELSE IF (iflag_lscav .EQ. 2) THEN ! frac_impa, frac_nucl
584 ! ********* modified old version
585 
586  d_tr_lessi_nucl(:,:,:) = 0.
587  d_tr_lessi_impa(:,:,:) = 0.
588  flestottr(:,:,:) = 0.
589 ! Tendance des aerosols nuclees et impactes
590  DO it = 1, nbtr
591  IF (aerosol(it)) THEN
592  his_dh(:)=0.
593  DO k = 1, klev
594  DO i = 1, klon
595 !PhH
596  zrho(i,k)=pplay(i,k)/t_seri(i,k)/rd
597  zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho(i,k)/rg
598 !
599  END DO
600  END DO
601 
602  DO k=klev-1, 1, -1
603  DO i=1, klon
604 ! d_tr_ls(i,k,it)=tr_seri(i,k,it)*(frac_impa(i,k)*frac_nucl(i,k)-1.)
605  dx=d_tr_ls(i,k,it)
606  his_dh(i)=his_dh(i)-dx*zrho(i,k)*zdz(i,k)/pdtphys ! kg/m2/s
607  evaplsc = prfl(i,k) - prfl(i,k+1) + psfl(i,k) - psfl(i,k+1)
608 ! Evaporation Partielle -> Liberation Partielle 0.5*evap
609  IF ( evaplsc .LT.0..and.abs(prfl(i,k+1)+psfl(i,k+1)).gt.1.e-10) THEN
610  evaplsc = (-evaplsc)/(prfl(i,k+1)+psfl(i,k+1))
611 ! evaplsc est donc positif, his_dh(i) est positif
612 !--------------
613  d_tr_evapls(i,k,it)=0.5*evaplsc*(d_tr_lessi_nucl(i,k+1,it) &
614  +d_tr_lessi_impa(i,k+1,it))
615 !------------- d_tr_evapls(i,k,it)=-0.5*evaplsc*(d_tr_lsc(i,k+1,it))
616  beta=0.5*evaplsc
617  if ((prfl(i,k)+psfl(i,k)).lt.1.e-10) THEN
618  beta=1.0*evaplsc
619  endif
620  dx=beta*his_dh(i)/zrho(i,k)/zdz(i,k)*pdtphys
621  his_dh(i)=(1.-beta)*his_dh(i) ! tracer from
622  d_tr_evapls(i,k,it)=dx
623  ENDIF
624  d_tr_ls(i,k,it)=tr_seri(i,k,it)*(frac_impa(i,k)*frac_nucl(i,k)-1.) &
625  +d_tr_evapls(i,k,it)
626 
627 !--------------
628  d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) + &
629  ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
630  d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) + &
631  ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
632 !
633 ! Flux lessivage total
634  flestottr(i,k,it) = flestottr(i,k,it) - &
635  ( d_tr_lessi_nucl(i,k,it) + &
636  d_tr_lessi_impa(i,k,it) ) * &
637  ( paprs(i,k)-paprs(i,k+1) ) / &
638  (rg * pdtphys)
639 !! Mise a jour des traceurs due a l'impaction,nucleation
640 ! tr_seri(i,k,it)=tr_seri(i,k,it)*frac_impa(i,k)*frac_nucl(i,k)
641 !! calcul de la tendance liee au lessivage stratiforme
642 ! d_tr_ls(i,k,it)=tr_seri(i,k,it)*&
643 ! (1.-1./(frac_impa(i,k)*frac_nucl(i,k)))
644 !--------------
645  END DO
646  END DO
647  END IF
648  END DO
649 ! ********* end modified old version
650 
651  ELSE IF (iflag_lscav .EQ. 1) THEN ! frac_impa, frac_nucl
652 ! ********* old version
653 
654  d_tr_lessi_nucl(:,:,:) = 0.
655  d_tr_lessi_impa(:,:,:) = 0.
656  flestottr(:,:,:) = 0.
657 !=========================
658 ! LESSIVAGE LARGE SCALE :
659 !=========================
660 
661 ! Tendance des aerosols nuclees et impactes
662 ! -----------------------------------------
663  DO it = 1, nbtr
664  IF (aerosol(it)) THEN
665  DO k = 1, klev
666  DO i = 1, klon
667  d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) + &
668  ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
669  d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) + &
670  ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
671 
672 !
673 ! Flux lessivage total
674 ! ------------------------------------------------------------
675  flestottr(i,k,it) = flestottr(i,k,it) - &
676  ( d_tr_lessi_nucl(i,k,it) + &
677  d_tr_lessi_impa(i,k,it) ) * &
678  ( paprs(i,k)-paprs(i,k+1) ) / &
679  (rg * pdtphys)
680 !
681 ! Mise a jour des traceurs due a l'impaction,nucleation
682 ! ----------------------------------------------------------------------
683  tr_seri(i,k,it)=tr_seri(i,k,it)*frac_impa(i,k)*frac_nucl(i,k)
684  END DO
685  END DO
686  END IF
687  END DO
688 
689 ! ********* end old version
690  ENDIF ! iflag_lscav . EQ. 1, 2, 3 or 4
691 !
692  END IF ! lessivage
693 
694 !=============================================================
695 ! Ecriture des sorties
696 !=============================================================
697 #ifdef CPP_IOIPSL
698  include "write_histrac.h"
699 #endif
700 
701 END SUBROUTINE phytrac