source: LMDZ6/trunk/libf/phylmd/phytrac_mod.F90 @ 4350

Last change on this file since 4350 was 4298, checked in by pcadule, 2 years ago

modifications to ensure restartability of the model when CO2 tracer is passed to radiative code, and to add diagnostics variables

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:keywords set to Id
File size: 45.7 KB
Line 
1!$Id: phytrac_mod.F90 4298 2022-10-17 08:15:06Z fhourdin $
2MODULE phytrac_mod
3!=================================================================================
4! Interface between the LMDZ physical package and tracer computation.
5! Chemistry modules (INCA, Reprobus or the more specific traclmdz routine)
6! are called from phytrac.
7!
8!======================================================================
9! Auteur(s) FH
10! Objet: Moniteur general des tendances traceurs
11!
12! iflag_vdf_trac : Options for activating transport by vertical diffusion :
13!     1. notmal
14!     0. emission is injected in the first layer only, without diffusion
15!    -1  no emission & no diffusion
16! Modification 2013/07/22 : transformed into a module to pass tendencies to
17!     physics outputs. Additional keys for controling activation of sub processes.
18! Modification R. Pilon 10 octobre 2012 large scale scavenging incloud_scav + bc_scav
19! Modification R. Pilon 01 janvier 2012 transport+scavenging KE scheme : cvltr
20!=================================================================================
21
22!
23! Tracer tendencies, for outputs
24!-------------------------------
25  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_cl  ! Td couche limite/traceur
26  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_dec                            !RomP
27  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_cv  ! Td convection/traceur
28! RomP >>>
29  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_insc
30  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_bcscav
31  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_evapls
32  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_ls
33  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_trsp
34  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_sscav
35  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_sat
36  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_uscav
37  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: qPr,qDi ! concentration tra dans pluie,air descente insaturee
38  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: qPa,qMel
39  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: qTrdi,dtrcvMA ! conc traceur descente air insaturee et td convective MA
40! RomP <<<
41  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_th  ! Td thermique
42  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_lessi_impa ! Td du lessivage par impaction
43  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_lessi_nucl ! Td du lessivage par nucleation
44  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE      :: qPrls      !jyg: concentration tra dans pluie LS a la surf.
45  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE      :: d_tr_dry ! Td depot sec/traceur (1st layer),ALLOCATABLE,SAVE  jyg
46  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE      :: flux_tr_dry ! depot sec/traceur (surface),ALLOCATABLE,SAVE    jyg
47
48!$OMP THREADPRIVATE(qPa,qMel,qTrdi,dtrcvMA,d_tr_th,d_tr_lessi_impa,d_tr_lessi_nucl)
49!$OMP THREADPRIVATE(d_tr_trsp,d_tr_sscav,d_tr_sat,d_tr_uscav,qPr,qDi)
50!$OMP THREADPRIVATE(d_tr_insc,d_tr_bcscav,d_tr_evapls,d_tr_ls,qPrls)
51!$OMP THREADPRIVATE(d_tr_cl,d_tr_dry,flux_tr_dry,d_tr_dec,d_tr_cv)
52
53
54CONTAINS
55
56  SUBROUTINE phytrac_init()
57    USE dimphy
58    USE infotrac_phy, ONLY: nbtr, types_trac
59    USE tracco2i_mod, ONLY: tracco2i_init
60    IMPLICIT NONE
61
62       ALLOCATE(d_tr_cl(klon,klev,nbtr),d_tr_dry(klon,nbtr))
63       ALLOCATE(flux_tr_dry(klon,nbtr),d_tr_dec(klon,klev,nbtr),d_tr_cv(klon,klev,nbtr))
64       ALLOCATE(d_tr_insc(klon,klev,nbtr),d_tr_bcscav(klon,klev,nbtr))
65       ALLOCATE(d_tr_evapls(klon,klev,nbtr),d_tr_ls(klon,klev,nbtr))
66       ALLOCATE(qPrls(klon,nbtr),d_tr_trsp(klon,klev,nbtr))
67       ALLOCATE(d_tr_sscav(klon,klev,nbtr),d_tr_sat(klon,klev,nbtr))
68       ALLOCATE(d_tr_uscav(klon,klev,nbtr),qPr(klon,klev,nbtr),qDi(klon,klev,nbtr))
69       ALLOCATE(qPa(klon,klev,nbtr),qMel(klon,klev,nbtr))
70       ALLOCATE(qTrdi(klon,klev,nbtr),dtrcvMA(klon,klev,nbtr))
71       ALLOCATE(d_tr_th(klon,klev,nbtr))
72       ALLOCATE(d_tr_lessi_impa(klon,klev,nbtr),d_tr_lessi_nucl(klon,klev,nbtr))
73
74
75
76    !===============================================================================
77    !    -- Do specific treatment according to chemestry model or local LMDZ tracers
78    !     
79    !===============================================================================
80    !   -- CO2 interactif --
81    IF(ANY(types_trac == 'co2i') .OR. ANY(types_trac == 'inco')) CALL tracco2i_init()
82
83       !   -- types_trac == 'co2i' ! PC
84       !   -- CO2 interactif --
85       !   -- source is updated with FF and BB emissions
86       !   -- and net fluxes from ocean and orchidee
87       !   -- sign convention : positive into the atmosphere
88
89  END SUBROUTINE phytrac_init
90
91  SUBROUTINE phytrac(                                 &
92       nstep,     julien,   gmtime,   debutphy,       &
93       lafin,     pdtphys,  u, v,     t_seri,         &
94       paprs,     pplay,    pmfu,     pmfd,           &
95       pen_u,     pde_u,    pen_d,    pde_d,          &
96       cdragh,    coefh,    fm_therm, entr_therm,     &
97       yu1,       yv1,      ftsol,    pctsrf,         &
98       ustar,     u10m,      v10m,                    &
99       wstar,     ale_bl,      ale_wake,              &
100       xlat,      xlon,                               &
101       frac_impa,frac_nucl,beta_fisrt,beta_v1,        &
102       presnivs,  pphis,    pphi,     albsol,         &
103       sh,        ch, rh,   cldfra,   rneb,           &
104       diafra,    cldliq,   itop_con, ibas_con,       &
105       pmflxr,    pmflxs,   prfl,     psfl,           &
106       da,        phi,      mp,       upwd,           &
107       phi2,      d1a,      dam,      sij, wght_cvfd, &   ! RomP +RL
108       wdtrainA,  wdtrainM, sigd,     clw, elij,      &   ! RomP
109       evap,      ep,       epmlmMm,  eplaMm,         &   ! RomP
110       dnwd,      aerosol_couple,     flxmass_w,      &
111       tau_aero,  piz_aero,  cg_aero, ccm,            &
112       rfname,                                        &
113       d_tr_dyn,                                      &   ! RomP
114       tr_seri, init_source)         
115    !
116    !======================================================================
117    ! Auteur(s) FH
118    ! Objet: Moniteur general des tendances traceurs
119    ! Modification R. Pilon 01 janvier 2012 transport+scavenging KE scheme : cvltr
120    ! Modification R. Pilon 10 octobre 2012 large scale scavenging incloud_scav + bc_scav
121    !======================================================================
122
123    USE ioipsl
124    USE phys_cal_mod, only : hour
125    USE dimphy
126    USE infotrac_phy, ONLY: nbtr, nqCO2, types_trac, type_trac, conv_flg, pbl_flg
127    USE strings_mod,  ONLY: int2str
128    USE mod_grid_phy_lmdz
129    USE mod_phys_lmdz_para
130    USE iophy
131    USE traclmdz_mod
132    USE tracinca_mod
133    USE tracreprobus_mod
134    USE indice_sol_mod
135    USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
136    USE print_control_mod, ONLY: lunout
137    USE aero_mod, ONLY : naero_grp
138
139    USE tracco2i_mod
140
141#ifdef CPP_StratAer
142    USE traccoag_mod
143    USE phys_local_var_mod, ONLY: mdw
144    USE phys_local_var_mod, ONLY: budg_dep_dry_ocs,   budg_dep_wet_ocs
145    USE phys_local_var_mod, ONLY: budg_dep_dry_so2,   budg_dep_wet_so2
146    USE phys_local_var_mod, ONLY: budg_dep_dry_h2so4, budg_dep_wet_h2so4
147    USE phys_local_var_mod, ONLY: budg_dep_dry_part,  budg_dep_wet_part
148    USE infotrac_phy, ONLY: nbtr_sulgas, id_OCS_strat, id_SO2_strat, id_H2SO4_strat
149    USE aerophys
150#endif
151
152    IMPLICIT NONE
153
154    INCLUDE "YOMCST.h"
155    INCLUDE "clesphys.h"
156    !==========================================================================
157    !                   -- ARGUMENT DESCRIPTION --
158    !==========================================================================
159
160    ! Input arguments
161    !----------------
162    !Configuration grille,temps:
163    INTEGER,INTENT(IN) :: nstep      ! Appel physique
164    INTEGER,INTENT(IN) :: julien     ! Jour julien
165    REAL,INTENT(IN)    :: gmtime     ! Heure courante
166    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde)
167    LOGICAL,INTENT(IN) :: debutphy   ! le flag de l'initialisation de la physique
168    LOGICAL,INTENT(IN) :: lafin      ! le flag de la fin de la physique
169
170    REAL,DIMENSION(klon),INTENT(IN) :: xlat    ! latitudes pour chaque point
171    REAL,DIMENSION(klon),INTENT(IN) :: xlon    ! longitudes pour chaque point
172    !
173    !Physique:
174    !--------
175    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
176    REAL,DIMENSION(klon,klev),INTENT(IN)   :: u       ! variable not used
177    REAL,DIMENSION(klon,klev),INTENT(IN)   :: v       ! variable not used
178    REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique
179    REAL,DIMENSION(klon,klev),INTENT(IN)   :: rh      ! humidite relative
180    REAL,DIMENSION(klon,klev),INTENT(IN)   :: ch      ! eau liquide (+ glace si le traceur existe)
181    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
182    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
183    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pphi    ! geopotentiel
184    REAL,DIMENSION(klon),INTENT(IN)        :: pphis
185    REAL,DIMENSION(klev),INTENT(IN)        :: presnivs
186    REAL,DIMENSION(klon,klev),INTENT(IN)   :: cldliq  ! eau liquide nuageuse
187    REAL,DIMENSION(klon,klev),INTENT(IN)   :: cldfra  ! fraction nuageuse (tous les nuages)
188    REAL,DIMENSION(klon,klev),INTENT(IN)   :: diafra  ! fraction nuageuse (convection ou stratus artificiels)
189    REAL,DIMENSION(klon,klev),INTENT(IN)   :: rneb    ! fraction nuageuse (grande echelle)
190    !
191    REAL                                   :: ql_incl ! contenu en eau liquide nuageuse dans le nuage ! ql_incl=oliq/rneb
192    REAL,DIMENSION(klon,klev),INTENT(IN)   :: beta_fisrt ! taux de conversion de l'eau cond (de fisrtilp)
193    REAL,DIMENSION(klon,klev),INTENT(out)  :: beta_v1    ! -- (originale version)
194
195    !
196    INTEGER,DIMENSION(klon),INTENT(IN)     :: itop_con
197    INTEGER,DIMENSION(klon),INTENT(IN)     :: ibas_con
198    REAL,DIMENSION(klon),INTENT(IN)        :: albsol  ! albedo surface
199    !
200    !Dynamique
201    !--------
202    REAL,DIMENSION(klon,klev,nbtr),INTENT(IN)    :: d_tr_dyn
203    !
204    !Convection:
205    !----------
206    REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfu  ! flux de masse dans le panache montant
207    REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfd  ! flux de masse dans le panache descendant
208    REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_u ! flux entraine dans le panache montant
209    REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_u ! flux detraine dans le panache montant
210    REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_d ! flux entraine dans le panache descendant
211    REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_d ! flux detraine dans le panache descendant
212
213    !...Tiedke     
214    REAL,DIMENSION(klon,klev+1),INTENT(IN)   :: pmflxr, pmflxs ! Flux precipitant de pluie, neige aux interfaces [convection]
215    REAL,DIMENSION(klon,klev+1),INTENT(IN)   :: prfl, psfl ! Flux precipitant de pluie, neige aux interfaces [large-scale]
216
217    LOGICAL,INTENT(IN)                       :: aerosol_couple
218    REAL,DIMENSION(klon,klev),INTENT(IN)     :: flxmass_w
219    REAL,DIMENSION(klon,klev,naero_grp,2),INTENT(IN) :: tau_aero
220    REAL,DIMENSION(klon,klev,naero_grp,2),INTENT(IN) :: piz_aero
221    REAL,DIMENSION(klon,klev,naero_grp,2),INTENT(IN) :: cg_aero
222    CHARACTER(len=4),DIMENSION(naero_grp),INTENT(IN) :: rfname
223    REAL,DIMENSION(klon,klev,2),INTENT(IN)   :: ccm
224    !... K.Emanuel
225    REAL,DIMENSION(klon,klev),INTENT(IN)     :: da
226    REAL,DIMENSION(klon,klev,klev),INTENT(IN):: phi
227    ! RomP >>>
228    REAL,DIMENSION(klon,klev),INTENT(IN)      :: d1a,dam
229    REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi2
230    !
231    REAL,DIMENSION(klon,klev),INTENT(IN)      :: wdtrainA
232    REAL,DIMENSION(klon,klev),INTENT(IN)      :: wdtrainM
233    REAL,DIMENSION(klon),INTENT(IN)           :: sigd
234    ! ---- RomP flux entraine, detraine et precipitant kerry Emanuel
235    REAL,DIMENSION(klon,klev),INTENT(IN)      :: evap
236    REAL,DIMENSION(klon,klev),INTENT(IN)      :: ep
237    REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: sij
238    REAL,DIMENSION(klon,klev),INTENT(IN)      :: wght_cvfd          !RL
239    REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: elij
240    REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: epmlmMm
241    REAL,DIMENSION(klon,klev),INTENT(IN)      :: eplaMm
242    REAL,DIMENSION(klon,klev),INTENT(IN)      :: clw
243    ! RomP <<<
244
245    !
246    REAL,DIMENSION(klon,klev),INTENT(IN)     :: mp
247    REAL,DIMENSION(klon,klev),INTENT(IN)     :: upwd      ! saturated updraft mass flux
248    REAL,DIMENSION(klon,klev),INTENT(IN)     :: dnwd      ! saturated downdraft mass flux
249    !
250    !Thermiques:
251    !----------
252    REAL,DIMENSION(klon,klev+1),INTENT(IN)   :: fm_therm
253    REAL,DIMENSION(klon,klev),INTENT(IN)     :: entr_therm
254    !
255    !Couche limite:
256    !--------------
257    !
258    REAL,DIMENSION(:),INTENT(IN)   :: cdragh          ! (klon) coeff drag pour T et Q
259    REAL,DIMENSION(:,:),INTENT(IN) :: coefh           ! (klon,klev) coeff melange CL (m**2/s)
260    REAL,DIMENSION(:),INTENT(IN)   :: ustar,u10m,v10m ! (klon) u* & vent a 10m (m/s)
261    REAL,DIMENSION(:),INTENT(IN)   :: wstar,ale_bl,ale_wake ! (klon) w* and Avail. Lifting Ener.
262    REAL,DIMENSION(:),INTENT(IN)   :: yu1             ! (klon) vents au premier niveau
263    REAL,DIMENSION(:),INTENT(IN)   :: yv1             ! (klon) vents au premier niveau
264    !
265    !Lessivage:
266    !----------
267    !
268    REAL, DIMENSION(:), ALLOCATABLE, SAVE :: ccntrAA
269    REAL, DIMENSION(:), ALLOCATABLE, SAVE :: ccntrENV
270    REAL, DIMENSION(:), ALLOCATABLE, SAVE :: coefcoli
271    LOGICAL, DIMENSION(:), ALLOCATABLE, SAVE :: flag_cvltr
272!$OMP THREADPRIVATE(ccntrAA,ccntrENV,coefcoli,flag_cvltr)
273    REAL, DIMENSION(klon,klev) :: ccntrAA_3d
274    REAL, DIMENSION(klon,klev) :: ccntrENV_3d
275    REAL, DIMENSION(klon,klev) :: coefcoli_3d
276    !
277    ! pour le ON-LINE
278    !
279    REAL,DIMENSION(klon,klev),INTENT(IN) :: frac_impa ! fraction d'aerosols non impactes
280    REAL,DIMENSION(klon,klev),INTENT(IN) :: frac_nucl ! fraction d'aerosols non nuclees
281
282    ! Arguments necessaires pour les sources et puits de traceur:
283    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol  ! Temperature du sol (surf)(Kelvin)
284    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol (nature du sol)
285
286#ifdef CPP_StratAer
287    REAL,DIMENSION(klon)           :: v_dep_dry !dry deposition velocity of stratospheric sulfate in m/s
288#endif
289    ! Output argument
290    !----------------
291    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
292    REAL,DIMENSION(klon,klev)                    :: sourceBE
293    REAL,DIMENSION(klon,nbtr), INTENT(IN) :: init_source
294
295    !=======================================================================================
296    !                        -- LOCAL VARIABLES --
297    !=======================================================================================
298
299    INTEGER :: i, k, it
300    INTEGER :: nsplit
301
302    !Sources et Reservoirs de traceurs (ex:Radon):
303    !--------------------------------------------
304    !
305    REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: source  ! a voir lorsque le flux de surface est prescrit
306!$OMP THREADPRIVATE(source)
307
308    !
309    !Entrees/Sorties:
310    !---------------
311    INTEGER                   :: iiq, ierr
312    INTEGER                   :: nhori, nvert
313    REAL                      :: zsto, zout, zjulian
314    INTEGER,SAVE              :: nid_tra     ! pointe vers le fichier histrac.nc         
315!$OMP THREADPRIVATE(nid_tra)
316    REAL,DIMENSION(klon)      :: zx_tmp_fi2d ! variable temporaire grille physique
317    INTEGER                   :: itau_w      ! pas de temps ecriture = nstep + itau_phy
318    LOGICAL,PARAMETER         :: ok_sync=.TRUE.
319    !
320    ! Nature du traceur
321    !------------------
322    LOGICAL,DIMENSION(:),ALLOCATABLE,SAVE :: aerosol  ! aerosol(it) = true  => aerosol => lessivage
323!$OMP THREADPRIVATE(aerosol)                        ! aerosol(it) = false => gaz
324    REAL,DIMENSION(klon,klev)             :: delp     ! epaisseur de couche (Pa)
325    !
326    ! Tendances de traceurs (Td) et flux de traceurs:
327    !------------------------
328    REAL,DIMENSION(klon,klev)      :: d_tr     ! Td dans l'atmosphere
329    REAL,DIMENSION(klon,klev)      :: Mint
330    REAL,DIMENSION(klon,klev,nbtr) :: zmfd1a
331    REAL,DIMENSION(klon,klev,nbtr) :: zmfdam
332    REAL,DIMENSION(klon,klev,nbtr) :: zmfphi2
333    ! Physique
334    !----------
335    REAL,DIMENSION(klon,klev,nbtr) :: flestottr ! flux de lessivage dans chaque couche
336    REAL,DIMENSION(klon,klev)      :: zmasse    ! densite atmospherique Kg/m2
337    REAL,DIMENSION(klon,klev)      :: ztra_th
338    !PhH
339    REAL,DIMENSION(klon,klev)      :: zrho
340    REAL,DIMENSION(klon,klev)      :: zdz
341    REAL                           :: evaplsc,dx,beta ! variable pour lessivage Genthon
342    REAL,DIMENSION(klon)           :: his_dh          ! ---
343    ! in-cloud scav variables
344    REAL           :: ql_incloud_ref     ! ref value of in-cloud condensed water content
345
346    !Controles:
347    !---------
348    INTEGER,SAVE :: iflag_vdf_trac,iflag_con_trac,iflag_the_trac
349    INTEGER,SAVE  :: iflag_con_trac_omp, iflag_vdf_trac_omp,iflag_the_trac_omp
350!$OMP THREADPRIVATE(iflag_vdf_trac,iflag_con_trac,iflag_the_trac)
351
352    LOGICAL,SAVE :: lessivage
353!$OMP THREADPRIVATE(lessivage)
354
355    !RomP >>>
356    INTEGER,SAVE  :: iflag_lscav_omp,iflag_lscav
357    REAL, SAVE ::   ccntrAA_in,ccntrAA_omp
358    REAL, SAVE ::   ccntrENV_in,ccntrENV_omp
359    REAL, SAVE ::   coefcoli_in,coefcoli_omp
360
361    LOGICAL,SAVE  :: convscav_omp,convscav
362!$OMP THREADPRIVATE(iflag_lscav)
363!$OMP THREADPRIVATE(ccntrAA_in,ccntrENV_in,coefcoli_in)
364!$OMP THREADPRIVATE(convscav)
365    !RomP <<<
366    !######################################################################
367    !                    -- INITIALIZATION --
368    !######################################################################
369
370    DO k=1,klev
371       DO i=1,klon
372          sourceBE(i,k)=0.
373          Mint(i,k)=0.
374          zrho(i,k)=0.
375          zdz(i,k)=0.
376       END DO
377    END DO
378
379    DO it=1, nbtr
380       DO k=1,klev
381          DO i=1,klon
382             d_tr_insc(i,k,it)=0.
383             d_tr_bcscav(i,k,it)=0.
384             d_tr_evapls(i,k,it)=0.
385             d_tr_ls(i,k,it)=0.
386             d_tr_cv(i,k,it)=0.
387             d_tr_cl(i,k,it)=0.
388             d_tr_trsp(i,k,it)=0.
389             d_tr_sscav(i,k,it)=0.
390             d_tr_sat(i,k,it)=0.
391             d_tr_uscav(i,k,it)=0.
392             d_tr_lessi_impa(i,k,it)=0.
393             d_tr_lessi_nucl(i,k,it)=0.
394             qDi(i,k,it)=0.
395             qPr(i,k,it)=0.
396             qPa(i,k,it)=0.
397             qMel(i,k,it)=0.
398             qTrdi(i,k,it)=0.
399             dtrcvMA(i,k,it)=0.
400             zmfd1a(i,k,it)=0.
401             zmfdam(i,k,it)=0.
402             zmfphi2(i,k,it)=0.
403          END DO
404       END DO
405    END DO
406
407    DO it=1, nbtr
408       DO i=1,klon
409          d_tr_dry(i,it)=0.
410          flux_tr_dry(i,it)=0.
411       END DO
412    END DO
413
414    DO k = 1, klev
415       DO i = 1, klon
416          delp(i,k) = paprs(i,k)-paprs(i,k+1)
417       END DO
418    END DO
419
420    IF (debutphy) THEN
421       !!jyg
422!$OMP BARRIER
423       ecrit_tra=86400. ! frequence de stokage en dur
424       ! obsolete car remplace par des ecritures dans phys_output_write
425       !RomP >>>
426       !
427       !Config Key  = convscav
428       !Config Desc = Convective scavenging switch: 0=off, 1=on.
429       !Config Def  = .FALSE.
430       !Config Help =
431       !
432!$OMP MASTER
433       convscav_omp=.FALSE.
434       call getin('convscav', convscav_omp)
435       iflag_vdf_trac_omp=1
436       call getin('iflag_vdf_trac', iflag_vdf_trac_omp)
437       iflag_con_trac_omp=1
438       call getin('iflag_con_trac', iflag_con_trac_omp)
439       iflag_the_trac_omp=1
440       call getin('iflag_the_trac', iflag_the_trac_omp)
441!$OMP END MASTER
442!$OMP BARRIER
443       convscav=convscav_omp
444       iflag_vdf_trac=iflag_vdf_trac_omp
445       iflag_con_trac=iflag_con_trac_omp
446       iflag_the_trac=iflag_the_trac_omp
447       write(lunout,*) 'phytrac passage dans routine conv avec lessivage', convscav
448       !
449       !Config Key  = iflag_lscav
450       !Config Desc = Large scale scavenging parametrization: 0=none, 1=old(Genthon92),
451       !              2=1+PHeinrich, 3=Reddy_Boucher2004, 4=3+RPilon.
452       !Config Def  = 1
453       !Config Help =
454       !
455!$OMP MASTER
456       iflag_lscav_omp=1
457       call getin('iflag_lscav', iflag_lscav_omp)
458       ccntrAA_omp=1
459       ccntrENV_omp=1.
460       coefcoli_omp=0.001
461       call getin('ccntrAA', ccntrAA_omp)
462       call getin('ccntrENV', ccntrENV_omp)
463       call getin('coefcoli', coefcoli_omp)
464!$OMP END MASTER
465!$OMP BARRIER
466       iflag_lscav=iflag_lscav_omp
467       ccntrAA_in=ccntrAA_omp
468       ccntrENV_in=ccntrENV_omp
469       coefcoli_in=coefcoli_omp
470       !
471       SELECT CASE(iflag_lscav)
472       CASE(0)
473          WRITE(lunout,*)  'Large scale scavenging: none'
474       CASE(1)
475          WRITE(lunout,*)  'Large scale scavenging: C. Genthon, Tellus(1992), 44B, 371-389'
476       CASE(2)
477          WRITE(lunout,*)  'Large scale scavenging: C. Genthon, modified P. Heinrich'
478       CASE(3)
479          WRITE(lunout,*)  'Large scale scavenging: M. Shekkar Reddy and O. Boucher, JGR(2004), 109, D14202'
480       CASE(4)
481          WRITE(lunout,*)  'Large scale scavenging: Reddy and Boucher, modified R. Pilon'
482       END SELECT
483       !RomP <<<
484       WRITE(*,*) 'FIRST TIME IN PHYTRAC : pdtphys(sec) = ',pdtphys,'ecrit_tra (sec) = ',ecrit_tra
485       ALLOCATE( source(klon,nbtr), stat=ierr)
486       IF (ierr /= 0) CALL abort_physic('phytrac', 'pb in allocation 1',1)
487
488       ALLOCATE( aerosol(nbtr), stat=ierr)
489       IF (ierr /= 0) CALL abort_physic('phytrac', 'pb in allocation 2',1)
490
491
492       ! Initialize module for specific tracers
493       IF(ANY(types_trac == 'inca')) THEN
494          source(:,:)=init_source(:,:)
495          CALL tracinca_init(aerosol,lessivage)
496       ELSE IF(ANY(types_trac == 'repr')) THEN
497          source(:,:)=0.
498       ELSE IF(ANY(types_trac == 'co2i')) THEN
499          source(:,:)=0.
500          lessivage  = .FALSE.
501          aerosol(:) = .FALSE.
502          pbl_flg(:) = 1
503          iflag_the_trac= 1
504          iflag_vdf_trac= 1
505          iflag_con_trac= 1
506       ELSE IF(ANY(types_trac == 'inco')) THEN
507          source(:,1:nqCO2) = 0.                          ! from CO2i ModThL
508          source(:,nqCO2+1:nbtr)=init_source(:,:)         ! from INCA ModThL
509          aerosol(1:nqCO2) = .FALSE.                      ! from CO2i ModThL
510          CALL tracinca_init(aerosol(nqCO2+1:nbtr),lessivage)     ! from INCA ModThL
511          pbl_flg(1:nqCO2) = 1                            ! From CO2i ModThL
512          iflag_the_trac = 1                              ! From CO2i
513          iflag_vdf_trac = 1                              ! From CO2i
514          iflag_con_trac = 1                              ! From CO2i
515#ifdef CPP_StratAer
516       ELSE IF(ANY(types_trac == 'coag')) THEN
517          source(:,:)=0.
518          DO it= 1, nbtr_sulgas
519            aerosol(it)=.FALSE.
520            IF (it==id_H2SO4_strat) aerosol(it)=.TRUE.
521          ENDDO
522          DO it= nbtr_sulgas+1, nbtr
523            aerosol(it)=.TRUE.
524          ENDDO
525#endif
526       ELSE IF(ANY(types_trac == 'lmdz')) THEN
527          CALL traclmdz_init(pctsrf,xlat,xlon,ftsol,tr_seri,t_seri,pplay,sh,pdtphys,aerosol,lessivage)
528       END IF
529
530       !
531       !--initialising coefficients for scavenging in the case of NP
532       !
533       ALLOCATE(flag_cvltr(nbtr))
534       IF (iflag_con.EQ.3) THEN
535          !
536          ALLOCATE(ccntrAA(nbtr))
537          ALLOCATE(ccntrENV(nbtr))
538          ALLOCATE(coefcoli(nbtr))
539          !
540          DO it=1, nbtr
541             IF(ANY(types_trac == 'repr')) THEN
542                 flag_cvltr(it)=.FALSE.
543             ELSE IF(ANY(types_trac == 'inca')) THEN
544!                IF ((it.EQ.id_Rn222) .OR. ((it.GE.id_SO2) .AND. (it.LE.id_NH3)) ) THEN
545!                   !--gas-phase species
546!                   flag_cvltr(it)=.FALSE.
547!
548!                ELSEIF ( (it.GE.id_CIDUSTM) .AND. (it.LE.id_AIN) ) THEN
549!                   !--insoluble aerosol species
550!                   flag_cvltr(it)=.TRUE.
551!                   ccntrAA(it)=0.7
552!                   ccntrENV(it)=0.7
553!                   coefcoli(it)=0.001
554!                ELSEIF ( (it.EQ.id_Pb210) .OR. ((it.GE.id_CSSSM) .AND. (it.LE.id_SSN))) THEN
555!                   !--soluble aerosol species
556!                   flag_cvltr(it)=.TRUE.
557!                   ccntrAA(it)=0.9
558!                   ccntrENV(it)=0.9
559!                   coefcoli(it)=0.001
560!                ELSE
561!                   WRITE(lunout,*) 'pb it=', it
562!                   CALL abort_physic('phytrac','pb it scavenging',1)
563!                ENDIF
564                !--test OB
565                !--for now we do not scavenge in cvltr
566                flag_cvltr(it)=.FALSE.
567             ELSE IF(ANY(types_trac == 'co2i')) THEN
568                !--co2 tracers are not scavenged
569                flag_cvltr(it)=.FALSE.
570             ELSE IF(ANY(types_trac == 'inco')) THEN     ! Add ThL
571                flag_cvltr(it)=.FALSE.
572#ifdef CPP_StratAer
573             ELSE IF(ANY(types_trac == 'coag')) THEN
574                IF (convscav.and.aerosol(it)) THEN
575                   flag_cvltr(it)=.TRUE.
576                   ccntrAA(it) =ccntrAA_in   
577                   ccntrENV(it)=ccntrENV_in
578                   coefcoli(it)=coefcoli_in
579                ELSE
580                   flag_cvltr(it)=.FALSE.
581                ENDIF
582#endif
583             ELSE IF(ANY(types_trac == 'lmdz')) THEN
584                IF (convscav.and.aerosol(it)) THEN
585                   flag_cvltr(it)=.TRUE.
586                   ccntrAA(it) =ccntrAA_in    !--a modifier par JYG a lire depuis fichier
587                   ccntrENV(it)=ccntrENV_in
588                   coefcoli(it)=coefcoli_in
589                ELSE
590                   flag_cvltr(it)=.FALSE.
591                ENDIF
592             END IF
593          ENDDO
594          !
595       ELSE ! iflag_con .ne. 3
596          flag_cvltr(:) = .FALSE.
597       ENDIF
598       !
599       ! print out all tracer flags
600       !
601       WRITE(lunout,*) 'print out all tracer flags'
602       WRITE(lunout,*) 'type_trac      =', type_trac
603       WRITE(lunout,*) 'config_inca    =', config_inca
604       WRITE(lunout,*) 'iflag_con_trac =', iflag_con_trac
605       WRITE(lunout,*) 'iflag_con      =', iflag_con
606       WRITE(lunout,*) 'convscav       =', convscav
607       WRITE(lunout,*) 'iflag_lscav    =', iflag_lscav
608       WRITE(lunout,*) 'aerosol        =', aerosol
609       WRITE(lunout,*) 'iflag_the_trac =', iflag_the_trac
610       WRITE(lunout,*) 'iflag_thermals =', iflag_thermals
611       WRITE(lunout,*) 'iflag_vdf_trac =', iflag_vdf_trac
612       WRITE(lunout,*) 'pbl_flg        =', pbl_flg
613       WRITE(lunout,*) 'lessivage      =', lessivage
614       write(lunout,*)  'flag_cvltr    = ', flag_cvltr
615
616       IF (lessivage .AND. (ANY(types_trac == 'inca') .OR. ANY(types_trac=='inco'))) THEN     ! Mod ThL
617          CALL abort_physic('phytrac', 'lessivage=T config_inca=inca impossible',1)
618!          STOP
619       ENDIF
620       !
621    ENDIF ! of IF (debutphy)
622    !############################################ END INITIALIZATION #######
623
624    DO k=1,klev
625       DO i=1,klon
626          zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
627       END DO
628    END DO
629    !
630    IF (id_be .GT. 0) THEN
631       DO k=1,klev
632          DO i=1,klon
633             sourceBE(i,k)=srcbe(i,k)       !RomP  -> pour sortie histrac
634          END DO
635       END DO
636    ENDIF
637
638    !===============================================================================
639    !    -- Do specific treatment according to chemestry model or local LMDZ tracers
640    !     
641    !===============================================================================
642    IF(ANY(types_trac == 'inca')) THEN
643       !    -- CHIMIE INCA  config_inca = aero or chem --
644       ! Appel fait en fin de phytrac pour avoir les emissions modifiees par
645       ! la couche limite et la convection avant le calcul de la chimie
646
647    ELSE IF(ANY(types_trac == 'repr')) THEN
648       !   -- CHIMIE REPROBUS --
649       CALL tracreprobus(pdtphys, gmtime, debutphy, julien, &
650            presnivs, xlat, xlon, pphis, pphi, &
651            t_seri, pplay, paprs, sh , &
652            tr_seri)
653
654    ELSE IF(ANY(types_trac == 'co2i')) THEN
655       !   -- CO2 interactif --
656       !   -- source is updated with FF and BB emissions
657       !   -- and net fluxes from ocean and orchidee
658       !   -- sign convention : positive into the atmosphere
659
660       CALL tracco2i(pdtphys, debutphy, &
661            xlat, xlon, pphis, pphi, &
662            t_seri, pplay, paprs, tr_seri, source)
663    ELSE IF(ANY(types_trac == 'inco')) THEN      ! Add ThL
664       CALL tracco2i(pdtphys, debutphy, &
665            xlat, xlon, pphis, pphi, &
666            t_seri, pplay, paprs, tr_seri, source)
667
668#ifdef CPP_StratAer
669    ELSE IF(ANY(types_trac == 'coag')) THEN
670       !   --STRATOSPHERIC AER IN THE STRAT --
671       CALL traccoag(pdtphys, gmtime, debutphy, julien, &
672            presnivs, xlat, xlon, pphis, pphi, &
673            t_seri, pplay, paprs, sh, rh , &
674            tr_seri)
675#endif
676    ELSE IF(ANY(types_trac == 'lmdz')) THEN
677       !    -- Traitement des traceurs avec traclmdz
678       CALL traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
679            cdragh,  coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon,iflag_vdf_trac>=0,sh, &
680            rh, pphi, ustar, wstar, ale_bl, ale_wake,  u10m, v10m, &
681            tr_seri, source, d_tr_cl,d_tr_dec, zmasse)               !RomP
682    END IF
683    !======================================================================
684    !       -- Calcul de l'effet de la convection --
685    !======================================================================
686
687    IF (iflag_con_trac==1) THEN
688
689       DO it=1, nbtr
690          IF ( conv_flg(it) == 0 ) CYCLE
691          IF (iflag_con.LT.2) THEN
692             !--pas de transport convectif
693             d_tr_cv(:,:,it)=0.
694
695          ELSE IF (iflag_con.EQ.2) THEN
696             !--ancien transport convectif de Tiedtke
697
698             CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
699                  pplay, paprs, tr_seri(:,:,it), d_tr_cv(:,:,it))
700          ELSE   
701             !--nouveau transport convectif de Emanuel
702
703             IF (flag_cvltr(it)) THEN
704                !--nouveau transport convectif de Emanuel avec lessivage convectif
705                !
706                !
707                ccntrAA_3d(:,:) =ccntrAA(it)
708                ccntrENV_3d(:,:)=ccntrENV(it)
709                coefcoli_3d(:,:)=coefcoli(it)
710
711                !--beware this interface is a bit weird because it is called for each tracer
712                !--with the full array tr_seri even if only item it is processed
713
714                CALL cvltr_scav(pdtphys, da, phi,phi2,d1a,dam, mp,ep,         &
715                     sigd,sij,wght_cvfd,clw,elij,epmlmMm,eplaMm,              &     
716                     pmflxr,pmflxs,evap,t_seri,wdtrainA,wdtrainM,             &   
717                     paprs,it,tr_seri,upwd,dnwd,itop_con,ibas_con,            &
718                     ccntrAA_3d,ccntrENV_3d,coefcoli_3d,                      &
719                     d_tr_cv,d_tr_trsp,d_tr_sscav,d_tr_sat,d_tr_uscav,qDi,qPr,&
720                     qPa,qMel,qTrdi,dtrcvMA,Mint,                             &
721                     zmfd1a,zmfphi2,zmfdam)
722
723
724             ELSE  !---flag_cvltr(it).EQ.FALSE
725                !--nouveau transport convectif de Emanuel mais pas de lessivage convectif
726
727                !--beware this interface is a bit weird because it is called for each tracer
728                !--with the full array tr_seri even if only item it is processed
729                !
730                CALL cvltr_noscav(it,pdtphys, da, phi,mp,wght_cvfd,paprs,pplay, &  !jyg
731                     tr_seri,upwd,dnwd,d_tr_cv)                                !jyg
732
733             ENDIF
734
735          ENDIF !--iflag
736
737          !--on ajoute les tendances
738
739          DO k = 1, klev
740             DO i = 1, klon       
741                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
742             END DO
743          END DO
744
745          CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'convection it = '//TRIM(int2str(it)))
746
747       END DO ! nbtr
748
749#ifdef CPP_StratAer
750       IF (ANY(types_trac=='coag')) THEN
751         ! initialize wet deposition flux of sulfur
752         budg_dep_wet_ocs(:)=0.0
753         budg_dep_wet_so2(:)=0.0
754         budg_dep_wet_h2so4(:)=0.0
755         budg_dep_wet_part(:)=0.0
756         ! compute wet deposition flux of sulfur (sum over gases and particles)
757         ! and convert to kg(S)/m2/s
758         DO i = 1, klon
759         DO k = 1, klev
760         DO it = 1, nbtr
761         !do not include SO2 because most of it comes trom the troposphere
762           IF (it==id_OCS_strat) THEN
763             budg_dep_wet_ocs(i)=budg_dep_wet_ocs(i)+d_tr_cv(i,k,it)*(mSatom/mOCSmol) &
764                            & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
765           ELSEIF (it==id_SO2_strat) THEN
766             budg_dep_wet_so2(i)=budg_dep_wet_so2(i)+d_tr_cv(i,k,it)*(mSatom/mSO2mol) &
767                            & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
768           ELSEIF (it==id_H2SO4_strat) THEN
769             budg_dep_wet_h2so4(i)=budg_dep_wet_h2so4(i)+d_tr_cv(i,k,it)*(mSatom/mH2SO4mol) &
770                            & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
771           ELSEIF (it.GT.nbtr_sulgas) THEN
772             budg_dep_wet_part(i)=budg_dep_wet_part(i)+d_tr_cv(i,k,it)*(mSatom/mH2SO4mol)  &
773                            & *dens_aer_dry*4./3.*RPI*(mdw(it-nbtr_sulgas)/2.)**3 &
774                            & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
775           ENDIF
776         ENDDO
777         ENDDO
778         ENDDO
779       ENDIF
780#endif
781
782    ENDIF ! convection
783
784    !======================================================================
785    !    -- Calcul de l'effet des thermiques --
786    !======================================================================
787
788    DO it=1,nbtr
789       DO k=1,klev
790          DO i=1,klon
791             d_tr_th(i,k,it)=0.
792             tr_seri(i,k,it)=MAX(tr_seri(i,k,it),0.)
793! the next safeguard causes some problem for stratospheric aerosol tracers (particle number)
794! and there is little justification for it so it is commented out (4 December 2017) by OB
795! if reinstated please keep the ifndef CPP_StratAer
796!#ifndef CPP_StratAer
797!             tr_seri(i,k,it)=MIN(tr_seri(i,k,it),1.e10)
798!#endif
799          END DO
800       END DO
801    END DO
802
803    IF (iflag_thermals.GT.0.AND.iflag_the_trac>0) THEN   
804
805       DO it=1, nbtr
806
807          CALL thermcell_dq(klon,klev,1,pdtphys,fm_therm,entr_therm, &
808               zmasse,tr_seri(1:klon,1:klev,it),        &
809               d_tr_th(1:klon,1:klev,it),ztra_th,0 )
810
811          DO k=1,klev
812             DO i=1,klon
813                d_tr_th(i,k,it)=pdtphys*d_tr_th(i,k,it)
814                tr_seri(i,k,it)=MAX(tr_seri(i,k,it)+d_tr_th(i,k,it),0.)
815             END DO
816          END DO
817
818       END DO ! it
819
820    ENDIF ! Thermiques
821
822    !======================================================================
823    !     -- Calcul de l'effet de la couche limite --
824    !======================================================================
825
826    IF (iflag_vdf_trac==1) THEN
827
828       !  Injection during BL mixing
829       !
830#ifdef CPP_StratAer
831       IF (ANY(types_trac=='coag')) THEN
832
833         ! initialize dry deposition flux of sulfur
834         budg_dep_dry_ocs(:)=0.0
835         budg_dep_dry_so2(:)=0.0
836         budg_dep_dry_h2so4(:)=0.0
837         budg_dep_dry_part(:)=0.0
838
839         ! compute dry deposition velocity as function of surface type (numbers
840         ! from IPSL note 23, 2002)
841         v_dep_dry(:) =  pctsrf(:,is_ter) * 2.5e-3 &
842                     & + pctsrf(:,is_oce) * 0.5e-3 &
843                     & + pctsrf(:,is_lic) * 2.5e-3 &
844                     & + pctsrf(:,is_sic) * 2.5e-3
845
846         ! compute surface dry deposition flux
847         zrho(:,1)=pplay(:,1)/t_seri(:,1)/RD
848
849         DO it=1, nbtr
850          source(:,it) = - v_dep_dry(:) * tr_seri(:,1,it) * zrho(:,1)
851         ENDDO
852
853       ENDIF
854#endif
855
856       DO it=1, nbtr
857          !
858          IF( pbl_flg(it) /= 0 ) THEN
859             !
860             CALL cltrac(pdtphys, coefh,t_seri,       &
861                  tr_seri(:,:,it), source(:,it),      &
862                  paprs, pplay, delp,                 &
863                  d_tr_cl(:,:,it),d_tr_dry(:,it),flux_tr_dry(:,it))
864             !
865             tr_seri(:,:,it)=tr_seri(:,:,it)+d_tr_cl(:,:,it)
866             !
867#ifdef CPP_StratAer
868             IF (ANY(types_trac=='coag')) THEN
869               ! compute dry deposition flux of sulfur (sum over gases and particles)
870               IF (it==id_OCS_strat) THEN
871                 budg_dep_dry_ocs(:)=budg_dep_dry_ocs(:)-source(:,it)*(mSatom/mOCSmol)
872               ELSEIF (it==id_SO2_strat) THEN
873                 budg_dep_dry_so2(:)=budg_dep_dry_so2(:)-source(:,it)*(mSatom/mSO2mol)
874               ELSEIF (it==id_H2SO4_strat) THEN
875                 budg_dep_dry_h2so4(:)=budg_dep_dry_h2so4(:)-source(:,it)*(mSatom/mH2SO4mol)
876               ELSEIF (it.GT.nbtr_sulgas) THEN
877                 budg_dep_dry_part(:)=budg_dep_dry_part(:)-source(:,it)*(mSatom/mH2SO4mol)*dens_aer_dry &
878                                & *4./3.*RPI*(mdw(it-nbtr_sulgas)/2.)**3
879               ENDIF
880             ENDIF
881#endif
882             !
883          ENDIF
884          !
885       ENDDO
886       !
887    ELSE IF (iflag_vdf_trac==0) THEN
888       !
889       !   Injection of source in the first model layer
890       !
891       DO it=1,nbtr
892          d_tr_cl(:,1,it)=source(:,it)*RG/delp(:,1)*pdtphys
893          tr_seri(:,1,it)=tr_seri(:,1,it)+d_tr_cl(:,1,it)
894       ENDDO
895       d_tr_cl(:,2:klev,1:nbtr)=0.
896       !
897    ELSE IF (iflag_vdf_trac==-1) THEN
898       !
899       ! Nothing happens
900       d_tr_cl=0.
901       !
902    ELSE
903       !
904       CALL abort_physic('iflag_vdf_trac', 'cas non prevu',1)
905       !
906    ENDIF ! couche limite
907
908    !======================================================================
909    !   Calcul de l'effet de la precipitation grande echelle
910    !   POUR INCA le lessivage est fait directement dans INCA
911    !======================================================================
912
913    IF (lessivage) THEN
914
915       ql_incloud_ref = 10.e-4
916       ql_incloud_ref =  5.e-4
917
918
919       ! calcul du contenu en eau liquide au sein du nuage
920       ql_incl = ql_incloud_ref
921       ! choix du lessivage
922       !
923       IF (iflag_lscav .EQ. 3 .OR. iflag_lscav .EQ. 4) THEN
924          ! ********  Olivier Boucher version (3) possibly with modified ql_incl (4)
925          !
926          DO it = 1, nbtr
927
928             IF (aerosol(it)) THEN
929             !  incloud scavenging and removal by large scale rain ! orig : ql_incl was replaced by 0.5e-3 kg/kg
930             ! the value 0.5e-3 kg/kg is from Giorgi and Chameides (1986), JGR
931             ! Liu (2001) proposed to use 1.5e-3 kg/kg
932
933!jyg<
934!!             CALL lsc_scav(pdtphys,it,iflag_lscav,ql_incl,prfl,psfl,rneb,beta_fisrt,  &
935             CALL lsc_scav(pdtphys,it,iflag_lscav,aerosol,ql_incl,prfl,psfl,rneb,beta_fisrt,  &
936!>jyg
937                  beta_v1,pplay,paprs,t_seri,tr_seri,d_tr_insc,   &
938                  d_tr_bcscav,d_tr_evapls,qPrls)
939
940             !large scale scavenging tendency
941             DO k = 1, klev
942                DO i = 1, klon
943                   d_tr_ls(i,k,it)=d_tr_insc(i,k,it)+d_tr_bcscav(i,k,it)+d_tr_evapls(i,k,it)
944                   tr_seri(i,k,it)=tr_seri(i,k,it)+d_tr_ls(i,k,it)
945                ENDDO
946             ENDDO
947             CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'lsc scav it = '//TRIM(int2str(it)))
948             ENDIF
949
950          END DO  !tr
951
952#ifdef CPP_StratAer
953         IF (ANY(types_trac=='coag')) THEN
954           ! compute wet deposition flux of sulfur (sum over gases and
955           ! particles) and convert to kg(S)/m2/s
956           ! adding contribution of d_tr_ls to d_tr_cv (above)
957           DO i = 1, klon
958           DO k = 1, klev
959           DO it = 1, nbtr
960             IF (it==id_OCS_strat) THEN
961               budg_dep_wet_ocs(i)=budg_dep_wet_ocs(i)+d_tr_ls(i,k,it)*(mSatom/mOCSmol) &
962                              & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
963             ELSEIF (it==id_SO2_strat) THEN
964               budg_dep_wet_so2(i)=budg_dep_wet_so2(i)+d_tr_ls(i,k,it)*(mSatom/mSO2mol) &
965                              & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
966             ELSEIF (it==id_H2SO4_strat) THEN
967               budg_dep_wet_h2so4(i)=budg_dep_wet_h2so4(i)+d_tr_ls(i,k,it)*(mSatom/mH2SO4mol) &
968                              & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
969             ELSEIF (it.GT.nbtr_sulgas) THEN
970               budg_dep_wet_part(i)=budg_dep_wet_part(i)+d_tr_ls(i,k,it)*(mSatom/mH2SO4mol)  &
971                              & *dens_aer_dry*4./3.*RPI*(mdw(it-nbtr_sulgas)/2.)**3 &
972                              & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
973             ENDIF
974           ENDDO
975           ENDDO
976           ENDDO
977         ENDIF
978#endif
979
980       ELSE IF (iflag_lscav .EQ. 2) THEN ! frac_impa, frac_nucl
981          ! *********   modified  old version
982
983          d_tr_lessi_nucl(:,:,:) = 0.
984          d_tr_lessi_impa(:,:,:) = 0.
985          flestottr(:,:,:) = 0.
986          ! Tendance des aerosols nuclees et impactes
987          DO it = 1, nbtr
988             IF (aerosol(it)) THEN
989                his_dh(:)=0.
990                DO k = 1, klev
991                   DO i = 1, klon
992                      !PhH
993                      zrho(i,k)=pplay(i,k)/t_seri(i,k)/RD
994                      zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho(i,k)/RG
995                      !
996                   ENDDO
997                ENDDO
998
999                DO k=klev-1, 1, -1
1000                   DO i=1, klon
1001                      !             d_tr_ls(i,k,it)=tr_seri(i,k,it)*(frac_impa(i,k)*frac_nucl(i,k)-1.)
1002                      dx=d_tr_ls(i,k,it)
1003                      his_dh(i)=his_dh(i)-dx*zrho(i,k)*zdz(i,k)/pdtphys !  kg/m2/s
1004                      evaplsc = prfl(i,k) - prfl(i,k+1) + psfl(i,k) - psfl(i,k+1)
1005                      ! Evaporation Partielle -> Liberation Partielle 0.5*evap
1006                      IF ( evaplsc .LT.0..and.abs(prfl(i,k+1)+psfl(i,k+1)).gt.1.e-10) THEN
1007                         evaplsc = (-evaplsc)/(prfl(i,k+1)+psfl(i,k+1))
1008                         ! evaplsc est donc positif, his_dh(i) est positif
1009                         !-------------- 
1010                         d_tr_evapls(i,k,it)=0.5*evaplsc*(d_tr_lessi_nucl(i,k+1,it) &
1011                              +d_tr_lessi_impa(i,k+1,it))
1012                         !-------------   d_tr_evapls(i,k,it)=-0.5*evaplsc*(d_tr_lsc(i,k+1,it))
1013                         beta=0.5*evaplsc
1014                         if ((prfl(i,k)+psfl(i,k)).lt.1.e-10) THEN
1015                            beta=1.0*evaplsc
1016                         endif
1017                         dx=beta*his_dh(i)/zrho(i,k)/zdz(i,k)*pdtphys
1018                         his_dh(i)=(1.-beta)*his_dh(i)   ! tracer from
1019                         d_tr_evapls(i,k,it)=dx
1020                      ENDIF
1021                      d_tr_ls(i,k,it)=tr_seri(i,k,it)*(frac_impa(i,k)*frac_nucl(i,k)-1.) &
1022                           +d_tr_evapls(i,k,it)
1023
1024                      !--------------
1025                      d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +    &
1026                           ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
1027                      d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +    &
1028                           ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
1029                      !
1030                      ! Flux lessivage total
1031                      flestottr(i,k,it) = flestottr(i,k,it) -   &
1032                           ( d_tr_lessi_nucl(i,k,it)   +        &
1033                           d_tr_lessi_impa(i,k,it) ) *          &
1034                           ( paprs(i,k)-paprs(i,k+1) ) /        &
1035                           (RG * pdtphys)
1036                      !! Mise a jour des traceurs due a l'impaction,nucleation
1037                      !                 tr_seri(i,k,it)=tr_seri(i,k,it)*frac_impa(i,k)*frac_nucl(i,k)
1038                      !!  calcul de la tendance liee au lessivage stratiforme
1039                      !                 d_tr_ls(i,k,it)=tr_seri(i,k,it)*&
1040                      !                                (1.-1./(frac_impa(i,k)*frac_nucl(i,k)))
1041                      !--------------
1042                   ENDDO
1043                ENDDO
1044             ENDIF
1045          ENDDO
1046          ! *********   end modified old version
1047
1048       ELSE IF (iflag_lscav .EQ. 1) THEN ! frac_impa, frac_nucl
1049          ! *********    old version
1050
1051          d_tr_lessi_nucl(:,:,:) = 0.
1052          d_tr_lessi_impa(:,:,:) = 0.
1053          flestottr(:,:,:) = 0.
1054          !=========================
1055          ! LESSIVAGE LARGE SCALE :
1056          !=========================
1057
1058          ! Tendance des aerosols nuclees et impactes
1059          ! -----------------------------------------
1060          DO it = 1, nbtr
1061             IF (aerosol(it)) THEN
1062                DO k = 1, klev
1063                   DO i = 1, klon
1064                      d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +    &
1065                           ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
1066                      d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +    &
1067                           ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
1068
1069                      !
1070                      ! Flux lessivage total
1071                      ! ------------------------------------------------------------
1072                      flestottr(i,k,it) = flestottr(i,k,it) -   &
1073                           ( d_tr_lessi_nucl(i,k,it)   +        &
1074                           d_tr_lessi_impa(i,k,it) ) *          &
1075                           ( paprs(i,k)-paprs(i,k+1) ) /        &
1076                           (RG * pdtphys)
1077                      !
1078                      ! Mise a jour des traceurs due a l'impaction,nucleation
1079                      ! ----------------------------------------------------------------------
1080                      tr_seri(i,k,it)=tr_seri(i,k,it)*frac_impa(i,k)*frac_nucl(i,k)
1081                   ENDDO
1082                ENDDO
1083             ENDIF
1084          ENDDO
1085
1086          ! *********   end old version
1087       ENDIF  !  iflag_lscav . EQ. 1, 2, 3 or 4
1088       !
1089    ENDIF !  lessivage
1090
1091
1092    !    -- CHIMIE INCA  config_inca = aero or chem --
1093    IF (ANY(types_trac == 'inca') .OR. ANY(types_trac == 'inco')) THEN  ! ModThL
1094
1095       CALL tracinca(&
1096            nstep,    julien,   gmtime,         lafin,     &
1097            pdtphys,  t_seri,   paprs,          pplay,     &
1098            pmfu,     upwd,     ftsol,  pctsrf, pphis,     &
1099            pphi,     albsol,   sh,    ch,     rh,        &
1100            cldfra,   rneb,     diafra,         cldliq,    &
1101            itop_con, ibas_con, pmflxr,         pmflxs,    &
1102            prfl,     psfl,     aerosol_couple, flxmass_w, &
1103            tau_aero, piz_aero, cg_aero,        ccm,       &
1104            rfname,                                        &
1105            tr_seri(:,:,1+nqCO2:nbtr),  source(:,1+nqCO2:nbtr))  ! ModThL 
1106    ENDIF
1107
1108  END SUBROUTINE phytrac
1109
1110END MODULE
Note: See TracBrowser for help on using the repository browser.