source: LMDZ6/branches/cirrus/libf/phylmd/phytrac_mod.F90 @ 5409

Last change on this file since 5409 was 4601, checked in by dcugnet, 18 months ago

StratAer? commit, N. Lebas

  • 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.2 KB
Line 
1!$Id: phytrac_mod.F90 4601 2023-06-30 22:07:30Z evignon $
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
53CONTAINS
54
55  SUBROUTINE phytrac_init()
56
57    USE dimphy
58    USE infotrac_phy, ONLY: nbtr, type_trac
59    USE tracco2i_mod, ONLY: tracco2i_init
60   
61    IMPLICIT NONE
62
63    ALLOCATE(d_tr_cl(klon,klev,nbtr),d_tr_dry(klon,nbtr))
64    ALLOCATE(flux_tr_dry(klon,nbtr),d_tr_dec(klon,klev,nbtr),d_tr_cv(klon,klev,nbtr))
65    ALLOCATE(d_tr_insc(klon,klev,nbtr),d_tr_bcscav(klon,klev,nbtr))
66    ALLOCATE(d_tr_evapls(klon,klev,nbtr),d_tr_ls(klon,klev,nbtr))
67    ALLOCATE(qPrls(klon,nbtr),d_tr_trsp(klon,klev,nbtr))
68    ALLOCATE(d_tr_sscav(klon,klev,nbtr),d_tr_sat(klon,klev,nbtr))
69    ALLOCATE(d_tr_uscav(klon,klev,nbtr),qPr(klon,klev,nbtr),qDi(klon,klev,nbtr))
70    ALLOCATE(qPa(klon,klev,nbtr),qMel(klon,klev,nbtr))
71    ALLOCATE(qTrdi(klon,klev,nbtr),dtrcvMA(klon,klev,nbtr))
72    ALLOCATE(d_tr_th(klon,klev,nbtr))
73    ALLOCATE(d_tr_lessi_impa(klon,klev,nbtr),d_tr_lessi_nucl(klon,klev,nbtr))
74
75    !===============================================================================
76    !    -- Do specific treatment according to chemestry model or local LMDZ tracers
77    !     
78    !===============================================================================
79    !   -- CO2 interactif --
80    IF (ANY(type_trac == ['co2i','inco'])) CALL tracco2i_init()
81
82       !   -- type_trac == 'co2i' ! PC
83       !   -- CO2 interactif --
84       !   -- source is updated with FF and BB emissions
85       !   -- and net fluxes from ocean and orchidee
86       !   -- sign convention : positive into the atmosphere
87
88  END SUBROUTINE phytrac_init
89
90  SUBROUTINE phytrac(                                 &
91       nstep,     julien,   gmtime,   debutphy,       &
92       lafin,     pdtphys,  u, v,     t_seri,         &
93       paprs,     pplay,    pmfu,     pmfd,           &
94       pen_u,     pde_u,    pen_d,    pde_d,          &
95       cdragh,    coefh,    fm_therm, entr_therm,     &
96       yu1,       yv1,      ftsol,    pctsrf,         &
97       ustar,     u10m,      v10m,                    &
98       wstar,     ale_bl,      ale_wake,              &
99       xlat,      xlon,                               &
100       frac_impa,frac_nucl,beta_fisrt,beta_v1,        &
101       presnivs,  pphis,    pphi,     albsol,         &
102       sh,        ch, rh,   cldfra,   rneb,           &
103       diafra,    cldliq,   itop_con, ibas_con,       &
104       pmflxr,    pmflxs,   prfl,     psfl,           &
105       da,        phi,      mp,       upwd,           &
106       phi2,      d1a,      dam,      sij, wght_cvfd, &   ! RomP +RL
107       wdtrainA,  wdtrainM, sigd,     clw, elij,      &   ! RomP
108       evap,      ep,       epmlmMm,  eplaMm,         &   ! RomP
109       dnwd,      aerosol_couple,     flxmass_w,      &
110       tau_aero,  piz_aero,  cg_aero, ccm,            &
111       rfname,                                        &
112       d_tr_dyn,                                      &   ! RomP
113       tr_seri, init_source)         
114    !
115    !======================================================================
116    ! Auteur(s) FH
117    ! Objet: Moniteur general des tendances traceurs
118    ! Modification R. Pilon 01 janvier 2012 transport+scavenging KE scheme : cvltr
119    ! Modification R. Pilon 10 octobre 2012 large scale scavenging incloud_scav + bc_scav
120    !======================================================================
121
122    USE ioipsl
123    USE phys_cal_mod, only : hour
124    USE dimphy
125    USE infotrac_phy, ONLY: nbtr, nqCO2, type_trac, conv_flg, pbl_flg
126    USE strings_mod,  ONLY: int2str
127    USE mod_grid_phy_lmdz
128    USE mod_phys_lmdz_para
129    USE iophy
130    USE traclmdz_mod
131    USE tracinca_mod
132    USE tracreprobus_mod
133    USE indice_sol_mod
134    USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
135    USE print_control_mod, ONLY: lunout
136    USE aero_mod, ONLY : naero_grp
137    USE lmdz_thermcell_dq, ONLY : thermcell_dq
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 strataer_nuc_mod, ONLY : tracstrataer_init
150    USE aerophys
151#endif
152
153    IMPLICIT NONE
154
155    INCLUDE "YOMCST.h"
156    INCLUDE "clesphys.h"
157    !==========================================================================
158    !                   -- ARGUMENT DESCRIPTION --
159    !==========================================================================
160
161    ! Input arguments
162    !----------------
163    !Configuration grille,temps:
164    INTEGER,INTENT(IN) :: nstep      ! Appel physique
165    INTEGER,INTENT(IN) :: julien     ! Jour julien
166    REAL,INTENT(IN)    :: gmtime     ! Heure courante
167    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde)
168    LOGICAL,INTENT(IN) :: debutphy   ! le flag de l'initialisation de la physique
169    LOGICAL,INTENT(IN) :: lafin      ! le flag de la fin de la physique
170
171    REAL,DIMENSION(klon),INTENT(IN) :: xlat    ! latitudes pour chaque point
172    REAL,DIMENSION(klon),INTENT(IN) :: xlon    ! longitudes pour chaque point
173    !
174    !Physique:
175    !--------
176    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
177    REAL,DIMENSION(klon,klev),INTENT(IN)   :: u       ! variable not used
178    REAL,DIMENSION(klon,klev),INTENT(IN)   :: v       ! variable not used
179    REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique
180    REAL,DIMENSION(klon,klev),INTENT(IN)   :: rh      ! humidite relative
181    REAL,DIMENSION(klon,klev),INTENT(IN)   :: ch      ! eau liquide (+ glace si le traceur existe)
182    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
183    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
184    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pphi    ! geopotentiel
185    REAL,DIMENSION(klon),INTENT(IN)        :: pphis
186    REAL,DIMENSION(klev),INTENT(IN)        :: presnivs
187    REAL,DIMENSION(klon,klev),INTENT(IN)   :: cldliq  ! eau condensee totale
188    REAL,DIMENSION(klon,klev),INTENT(IN)   :: cldfra  ! fraction nuageuse (tous les nuages)
189    REAL,DIMENSION(klon,klev),INTENT(IN)   :: diafra  ! fraction nuageuse (convection ou stratus artificiels)
190    REAL,DIMENSION(klon,klev),INTENT(IN)   :: rneb    ! fraction nuageuse (grande echelle)
191    !
192    REAL                                   :: ql_incl ! contenu en eau liquide nuageuse dans le nuage ! ql_incl=oliq/rneb
193    REAL,DIMENSION(klon,klev),INTENT(IN)   :: beta_fisrt ! taux de conversion de l'eau cond (de fisrtilp)
194    REAL,DIMENSION(klon,klev),INTENT(out)  :: beta_v1    ! -- (originale version)
195
196    !
197    INTEGER,DIMENSION(klon),INTENT(IN)     :: itop_con
198    INTEGER,DIMENSION(klon),INTENT(IN)     :: ibas_con
199    REAL,DIMENSION(klon),INTENT(IN)        :: albsol  ! albedo surface
200    !
201    !Dynamique
202    !--------
203    REAL,DIMENSION(klon,klev,nbtr),INTENT(IN)    :: d_tr_dyn
204    !
205    !Convection:
206    !----------
207    REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfu  ! flux de masse dans le panache montant
208    REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfd  ! flux de masse dans le panache descendant
209    REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_u ! flux entraine dans le panache montant
210    REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_u ! flux detraine dans le panache montant
211    REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_d ! flux entraine dans le panache descendant
212    REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_d ! flux detraine dans le panache descendant
213
214    !...Tiedke     
215    REAL,DIMENSION(klon,klev+1),INTENT(IN)   :: pmflxr, pmflxs ! Flux precipitant de pluie, neige aux interfaces [convection]
216    REAL,DIMENSION(klon,klev+1),INTENT(IN)   :: prfl, psfl ! Flux precipitant de pluie, neige aux interfaces [large-scale]
217
218    LOGICAL,INTENT(IN)                       :: aerosol_couple
219    REAL,DIMENSION(klon,klev),INTENT(IN)     :: flxmass_w
220    REAL,DIMENSION(klon,klev,naero_grp,2),INTENT(IN) :: tau_aero
221    REAL,DIMENSION(klon,klev,naero_grp,2),INTENT(IN) :: piz_aero
222    REAL,DIMENSION(klon,klev,naero_grp,2),INTENT(IN) :: cg_aero
223    CHARACTER(len=4),DIMENSION(naero_grp),INTENT(IN) :: rfname
224    REAL,DIMENSION(klon,klev,2),INTENT(IN)   :: ccm
225    !... K.Emanuel
226    REAL,DIMENSION(klon,klev),INTENT(IN)     :: da
227    REAL,DIMENSION(klon,klev,klev),INTENT(IN):: phi
228    ! RomP >>>
229    REAL,DIMENSION(klon,klev),INTENT(IN)      :: d1a,dam
230    REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi2
231    !
232    REAL,DIMENSION(klon,klev),INTENT(IN)      :: wdtrainA
233    REAL,DIMENSION(klon,klev),INTENT(IN)      :: wdtrainM
234    REAL,DIMENSION(klon),INTENT(IN)           :: sigd
235    ! ---- RomP flux entraine, detraine et precipitant kerry Emanuel
236    REAL,DIMENSION(klon,klev),INTENT(IN)      :: evap
237    REAL,DIMENSION(klon,klev),INTENT(IN)      :: ep
238    REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: sij
239    REAL,DIMENSION(klon,klev),INTENT(IN)      :: wght_cvfd          !RL
240    REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: elij
241    REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: epmlmMm
242    REAL,DIMENSION(klon,klev),INTENT(IN)      :: eplaMm
243    REAL,DIMENSION(klon,klev),INTENT(IN)      :: clw
244    ! RomP <<<
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(INOUT)     :: 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       ENDDO
377    ENDDO
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          ENDDO
404       ENDDO
405    ENDDO
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       ENDDO
412    ENDDO
413
414    DO k = 1, klev
415       DO i = 1, klon
416          delp(i,k) = paprs(i,k)-paprs(i,k+1)
417       ENDDO
418    ENDDO
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 (type_trac == 'inca') THEN
494          source(:,:)=init_source(:,:)
495          CALL tracinca_init(aerosol,lessivage)
496       ELSE IF (type_trac == 'repr') THEN
497          source(:,:)=0.
498       ELSE IF (type_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 (type_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 (type_trac == 'coag') THEN
517          source(:,:)=0.
518          CALL tracstrataer_init(aerosol,lessivage)       ! init aerosols and lessivage param
519#endif
520       ELSE IF (type_trac == 'lmdz') THEN
521          CALL traclmdz_init(pctsrf,xlat,xlon,ftsol,tr_seri,t_seri,pplay,sh,pdtphys,aerosol,lessivage)
522       ENDIF
523
524       !
525       !--initialising coefficients for scavenging in the case of NP
526       !
527       ALLOCATE(flag_cvltr(nbtr))
528       IF (iflag_con.EQ.3) THEN
529          !
530          ALLOCATE(ccntrAA(nbtr))
531          ALLOCATE(ccntrENV(nbtr))
532          ALLOCATE(coefcoli(nbtr))
533          !
534          DO it=1, nbtr
535             IF (type_trac == 'repr') THEN
536                 flag_cvltr(it)=.FALSE.
537             ELSE IF (type_trac == 'inca') THEN
538!                IF ((it.EQ.id_Rn222) .OR. ((it.GE.id_SO2) .AND. (it.LE.id_NH3)) ) THEN
539!                   !--gas-phase species
540!                   flag_cvltr(it)=.FALSE.
541!
542!                ELSEIF ( (it.GE.id_CIDUSTM) .AND. (it.LE.id_AIN) ) THEN
543!                   !--insoluble aerosol species
544!                   flag_cvltr(it)=.TRUE.
545!                   ccntrAA(it)=0.7
546!                   ccntrENV(it)=0.7
547!                   coefcoli(it)=0.001
548!                ELSEIF ( (it.EQ.id_Pb210) .OR. ((it.GE.id_CSSSM) .AND. (it.LE.id_SSN))) THEN
549!                   !--soluble aerosol species
550!                   flag_cvltr(it)=.TRUE.
551!                   ccntrAA(it)=0.9
552!                   ccntrENV(it)=0.9
553!                   coefcoli(it)=0.001
554!                ELSE
555!                   WRITE(lunout,*) 'pb it=', it
556!                   CALL abort_physic('phytrac','pb it scavenging',1)
557!                ENDIF
558                !--test OB
559                !--for now we do not scavenge in cvltr
560                flag_cvltr(it)=.FALSE.
561             ELSE IF (type_trac == 'co2i') THEN
562                !--co2 tracers are not scavenged
563                flag_cvltr(it)=.FALSE.
564             ELSE IF (type_trac == 'inco') THEN     ! Add ThL
565                flag_cvltr(it)=.FALSE.
566#ifdef CPP_StratAer
567             ELSE IF (type_trac == 'coag') THEN
568                IF (convscav.and.aerosol(it)) THEN
569                   flag_cvltr(it)=.TRUE.
570                   ccntrAA(it) =ccntrAA_in   
571                   ccntrENV(it)=ccntrENV_in
572                   coefcoli(it)=coefcoli_in
573                ELSE
574                   flag_cvltr(it)=.FALSE.
575                ENDIF
576#endif
577             ELSE IF (type_trac == 'lmdz') THEN
578                IF (convscav.and.aerosol(it)) THEN
579                   flag_cvltr(it)=.TRUE.
580                   ccntrAA(it) =ccntrAA_in    !--a modifier par JYG a lire depuis fichier
581                   ccntrENV(it)=ccntrENV_in
582                   coefcoli(it)=coefcoli_in
583                ELSE
584                   flag_cvltr(it)=.FALSE.
585                ENDIF
586             ENDIF
587          ENDDO
588          !
589       ELSE ! iflag_con .ne. 3
590          flag_cvltr(:) = .FALSE.
591       ENDIF
592       !
593       ! print out all tracer flags
594       !
595       WRITE(lunout,*) 'print out all tracer flags'
596       WRITE(lunout,*) 'type_trac      =', type_trac
597       WRITE(lunout,*) 'config_inca    =', config_inca
598       WRITE(lunout,*) 'iflag_con_trac =', iflag_con_trac
599       WRITE(lunout,*) 'iflag_con      =', iflag_con
600       WRITE(lunout,*) 'convscav       =', convscav
601       WRITE(lunout,*) 'iflag_lscav    =', iflag_lscav
602       WRITE(lunout,*) 'aerosol        =', aerosol
603       WRITE(lunout,*) 'iflag_the_trac =', iflag_the_trac
604       WRITE(lunout,*) 'iflag_thermals =', iflag_thermals
605       WRITE(lunout,*) 'iflag_vdf_trac =', iflag_vdf_trac
606       WRITE(lunout,*) 'pbl_flg        =', pbl_flg
607       WRITE(lunout,*) 'lessivage      =', lessivage
608       write(lunout,*)  'flag_cvltr    = ', flag_cvltr
609
610       IF (lessivage .AND. ANY(type_trac == ['inca','inco'])) &
611          CALL abort_physic('phytrac', 'lessivage=T config_inca=inca impossible',1)
612       !
613    ENDIF ! of IF (debutphy)
614    !############################################ END INITIALIZATION #######
615
616    DO k=1,klev
617       DO i=1,klon
618          zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
619       ENDDO
620    ENDDO
621    !
622    IF (id_be .GT. 0) THEN
623       DO k=1,klev
624          DO i=1,klon
625             sourceBE(i,k)=srcbe(i,k)       !RomP  -> pour sortie histrac
626          ENDDO
627       ENDDO
628    ENDIF
629
630    !===============================================================================
631    !    -- Do specific treatment according to chemestry model or local LMDZ tracers
632    !     
633    !===============================================================================
634    IF (type_trac == 'inca') THEN
635       !    -- CHIMIE INCA  config_inca = aero or chem --
636       ! Appel fait en fin de phytrac pour avoir les emissions modifiees par
637       ! la couche limite et la convection avant le calcul de la chimie
638
639    ELSE IF (type_trac == 'repr') THEN
640       !   -- CHIMIE REPROBUS --
641       CALL tracreprobus(pdtphys, gmtime, debutphy, julien, &
642            presnivs, xlat, xlon, pphis, pphi, &
643            t_seri, pplay, paprs, sh , &
644            tr_seri)
645
646    ELSE IF (type_trac == 'co2i') THEN
647       !   -- CO2 interactif --
648       !   -- source is updated with FF and BB emissions
649       !   -- and net fluxes from ocean and orchidee
650       !   -- sign convention : positive into the atmosphere
651
652       CALL tracco2i(pdtphys, debutphy, &
653            xlat, xlon, pphis, pphi, &
654            t_seri, pplay, paprs, tr_seri, source)
655    ELSE IF (type_trac == 'inco') THEN      ! Add ThL
656       CALL tracco2i(pdtphys, debutphy, &
657            xlat, xlon, pphis, pphi, &
658            t_seri, pplay, paprs, tr_seri, source)
659
660#ifdef CPP_StratAer
661    ELSE IF (type_trac == 'coag') THEN
662       !   --STRATOSPHERIC AER IN THE STRAT --
663       CALL traccoag(pdtphys, gmtime, debutphy, julien, &
664            presnivs, xlat, xlon, pphis, pphi, &
665            t_seri, pplay, paprs, sh, rh , &
666            tr_seri)
667#endif
668    ELSE IF (type_trac == 'lmdz') THEN
669       !    -- Traitement des traceurs avec traclmdz
670       CALL traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
671            cdragh,  coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon,iflag_vdf_trac>=0,sh, &
672            rh, pphi, ustar, wstar, ale_bl, ale_wake,  u10m, v10m, &
673            tr_seri, source, d_tr_cl,d_tr_dec, zmasse)               !RomP
674    ENDIF
675    !======================================================================
676    !       -- Calcul de l'effet de la convection --
677    !======================================================================
678
679    IF (iflag_con_trac==1) THEN
680
681       DO it=1, nbtr
682          IF ( conv_flg(it) == 0 ) CYCLE
683          IF (iflag_con.LT.2) THEN
684             !--pas de transport convectif
685             d_tr_cv(:,:,it)=0.
686
687          ELSE IF (iflag_con.EQ.2) THEN
688             !--ancien transport convectif de Tiedtke
689
690             CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
691                  pplay, paprs, tr_seri(:,:,it), d_tr_cv(:,:,it))
692          ELSE   
693             !--nouveau transport convectif de Emanuel
694
695             IF (flag_cvltr(it)) THEN
696                !--nouveau transport convectif de Emanuel avec lessivage convectif
697                !
698                !
699                ccntrAA_3d(:,:) =ccntrAA(it)
700                ccntrENV_3d(:,:)=ccntrENV(it)
701                coefcoli_3d(:,:)=coefcoli(it)
702
703                !--beware this interface is a bit weird because it is called for each tracer
704                !--with the full array tr_seri even if only item it is processed
705
706                CALL cvltr_scav(pdtphys, da, phi,phi2,d1a,dam, mp,ep,         &
707                     sigd,sij,wght_cvfd,clw,elij,epmlmMm,eplaMm,              &     
708                     pmflxr,pmflxs,evap,t_seri,wdtrainA,wdtrainM,             &   
709                     paprs,it,tr_seri,upwd,dnwd,itop_con,ibas_con,            &
710                     ccntrAA_3d,ccntrENV_3d,coefcoli_3d,                      &
711                     d_tr_cv,d_tr_trsp,d_tr_sscav,d_tr_sat,d_tr_uscav,qDi,qPr,&
712                     qPa,qMel,qTrdi,dtrcvMA,Mint,                             &
713                     zmfd1a,zmfphi2,zmfdam)
714
715
716             ELSE  !---flag_cvltr(it).EQ.FALSE
717                !--nouveau transport convectif de Emanuel mais pas de lessivage convectif
718
719                !--beware this interface is a bit weird because it is called for each tracer
720                !--with the full array tr_seri even if only item it is processed
721                !
722                CALL cvltr_noscav(it,pdtphys, da, phi,mp,wght_cvfd,paprs,pplay, &  !jyg
723                     tr_seri,upwd,dnwd,d_tr_cv)                                !jyg
724
725             ENDIF
726
727          ENDIF !--iflag
728
729          !--on ajoute les tendances
730
731          DO k = 1, klev
732             DO i = 1, klon       
733                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
734             ENDDO
735          ENDDO
736
737          CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'convection it = '//TRIM(int2str(it)))
738
739       ENDDO ! nbtr
740
741#ifdef CPP_StratAer
742       IF (type_trac=='coag') THEN
743         ! initialize wet deposition flux of sulfur
744         budg_dep_wet_ocs(:)=0.0
745         budg_dep_wet_so2(:)=0.0
746         budg_dep_wet_h2so4(:)=0.0
747         budg_dep_wet_part(:)=0.0
748         ! compute wet deposition flux of sulfur (sum over gases and particles)
749         ! and convert to kg(S)/m2/s
750         DO i = 1, klon
751         DO k = 1, klev
752         DO it = 1, nbtr
753         !do not include SO2 because most of it comes trom the troposphere
754           IF (it==id_OCS_strat) THEN
755             budg_dep_wet_ocs(i)=budg_dep_wet_ocs(i)+d_tr_cv(i,k,it)*(mSatom/mOCSmol) &
756                            & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
757           ELSEIF (it==id_SO2_strat) THEN
758             budg_dep_wet_so2(i)=budg_dep_wet_so2(i)+d_tr_cv(i,k,it)*(mSatom/mSO2mol) &
759                            & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
760           ELSEIF (it==id_H2SO4_strat) THEN
761             budg_dep_wet_h2so4(i)=budg_dep_wet_h2so4(i)+d_tr_cv(i,k,it)*(mSatom/mH2SO4mol) &
762                            & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
763           ELSEIF (it.GT.nbtr_sulgas) THEN
764             budg_dep_wet_part(i)=budg_dep_wet_part(i)+d_tr_cv(i,k,it)*(mSatom/mH2SO4mol)  &
765                            & *dens_aer_dry*4./3.*RPI*(mdw(it-nbtr_sulgas)/2.)**3 &
766                            & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
767           ENDIF
768         ENDDO
769         ENDDO
770         ENDDO
771       ENDIF
772#endif
773
774    ENDIF ! convection
775
776    !======================================================================
777    !    -- Calcul de l'effet des thermiques --
778    !======================================================================
779
780    DO it=1,nbtr
781       DO k=1,klev
782          DO i=1,klon
783             d_tr_th(i,k,it)=0.
784             tr_seri(i,k,it)=MAX(tr_seri(i,k,it),0.)
785! the next safeguard causes some problem for stratospheric aerosol tracers (particle number)
786! and there is little justification for it so it is commented out (4 December 2017) by OB
787! if reinstated please keep the ifndef CPP_StratAer
788!#ifndef CPP_StratAer
789!             tr_seri(i,k,it)=MIN(tr_seri(i,k,it),1.e10)
790!#endif
791          ENDDO
792       ENDDO
793    ENDDO
794
795    IF (iflag_thermals.GT.0.AND.iflag_the_trac>0) THEN   
796
797       DO it=1, nbtr
798
799          CALL thermcell_dq(klon,klev,1,pdtphys,fm_therm,entr_therm, &
800               zmasse,tr_seri(1:klon,1:klev,it),        &
801               d_tr_th(1:klon,1:klev,it),ztra_th,0 )
802
803          DO k=1,klev
804             DO i=1,klon
805                d_tr_th(i,k,it)=pdtphys*d_tr_th(i,k,it)
806                tr_seri(i,k,it)=MAX(tr_seri(i,k,it)+d_tr_th(i,k,it),0.)
807             ENDDO
808          ENDDO
809
810       ENDDO ! it
811
812    ENDIF ! Thermiques
813
814    !======================================================================
815    !     -- Calcul de l'effet de la couche limite --
816    !======================================================================
817
818    IF (iflag_vdf_trac==1) THEN
819
820       !  Injection during BL mixing
821       !
822#ifdef CPP_StratAer
823       IF (type_trac=='coag') THEN
824
825         ! initialize dry deposition flux of sulfur
826         budg_dep_dry_ocs(:)=0.0
827         budg_dep_dry_so2(:)=0.0
828         budg_dep_dry_h2so4(:)=0.0
829         budg_dep_dry_part(:)=0.0
830
831         ! compute dry deposition velocity as function of surface type (numbers
832         ! from IPSL note 23, 2002)
833         v_dep_dry(:) =  pctsrf(:,is_ter) * 2.5e-3 &
834                     & + pctsrf(:,is_oce) * 0.5e-3 &
835                     & + pctsrf(:,is_lic) * 2.5e-3 &
836                     & + pctsrf(:,is_sic) * 2.5e-3
837
838         ! compute surface dry deposition flux
839         zrho(:,1)=pplay(:,1)/t_seri(:,1)/RD
840
841         DO it=1, nbtr
842          source(:,it) = - v_dep_dry(:) * tr_seri(:,1,it) * zrho(:,1)
843         ENDDO
844
845       ENDIF
846#endif
847
848       DO it=1, nbtr
849          !
850          IF (pbl_flg(it) /= 0) THEN
851             !
852             CALL cltrac(pdtphys, coefh,t_seri,       &
853                  tr_seri(:,:,it), source(:,it),      &
854                  paprs, pplay, delp,                 &
855                  d_tr_cl(:,:,it),d_tr_dry(:,it),flux_tr_dry(:,it))
856             !
857             tr_seri(:,:,it)=tr_seri(:,:,it)+d_tr_cl(:,:,it)
858             !
859#ifdef CPP_StratAer
860             IF (type_trac=='coag') THEN
861               ! compute dry deposition flux of sulfur (sum over gases and particles)
862               IF (it==id_OCS_strat) THEN
863                 budg_dep_dry_ocs(:)=budg_dep_dry_ocs(:)-source(:,it)*(mSatom/mOCSmol)
864               ELSEIF (it==id_SO2_strat) THEN
865                 budg_dep_dry_so2(:)=budg_dep_dry_so2(:)-source(:,it)*(mSatom/mSO2mol)
866               ELSEIF (it==id_H2SO4_strat) THEN
867                 budg_dep_dry_h2so4(:)=budg_dep_dry_h2so4(:)-source(:,it)*(mSatom/mH2SO4mol)
868               ELSEIF (it.GT.nbtr_sulgas) THEN
869                 budg_dep_dry_part(:)=budg_dep_dry_part(:)-source(:,it)*(mSatom/mH2SO4mol)*dens_aer_dry &
870                                & *4./3.*RPI*(mdw(it-nbtr_sulgas)/2.)**3
871               ENDIF
872             ENDIF
873#endif
874             !
875          ENDIF
876          !
877       ENDDO
878       !
879    ELSE IF (iflag_vdf_trac==0) THEN
880       !
881       !   Injection of source in the first model layer
882       !
883       DO it=1,nbtr
884          d_tr_cl(:,1,it)=source(:,it)*RG/delp(:,1)*pdtphys
885          tr_seri(:,1,it)=tr_seri(:,1,it)+d_tr_cl(:,1,it)
886       ENDDO
887       d_tr_cl(:,2:klev,1:nbtr)=0.
888       !
889    ELSE IF (iflag_vdf_trac==-1) THEN
890       !
891       ! Nothing happens
892       d_tr_cl=0.
893       !
894    ELSE
895       !
896       CALL abort_physic('iflag_vdf_trac', 'cas non prevu',1)
897       !
898    ENDIF ! couche limite
899
900    !======================================================================
901    !   Calcul de l'effet de la precipitation grande echelle
902    !   POUR INCA le lessivage est fait directement dans INCA
903    !======================================================================
904
905    IF (lessivage) THEN
906
907       ql_incloud_ref = 10.e-4
908       ql_incloud_ref =  5.e-4
909
910
911       ! calcul du contenu en eau liquide au sein du nuage
912       ql_incl = ql_incloud_ref
913       ! choix du lessivage
914       !
915       IF (iflag_lscav .EQ. 3 .OR. iflag_lscav .EQ. 4) THEN
916          ! ********  Olivier Boucher version (3) possibly with modified ql_incl (4)
917          !
918          DO it = 1, nbtr
919
920             IF (aerosol(it)) THEN
921             !  incloud scavenging and removal by large scale rain ! orig : ql_incl was replaced by 0.5e-3 kg/kg
922             ! the value 0.5e-3 kg/kg is from Giorgi and Chameides (1986), JGR
923             ! Liu (2001) proposed to use 1.5e-3 kg/kg
924
925             CALL lsc_scav(pdtphys,it,iflag_lscav,aerosol,ql_incl,prfl,psfl,rneb,beta_fisrt,  &
926                           beta_v1,pplay,paprs,t_seri,tr_seri,d_tr_insc,d_tr_bcscav,d_tr_evapls,qPrls)
927
928             !large scale scavenging tendency
929             DO k = 1, klev
930                DO i = 1, klon
931                   d_tr_ls(i,k,it)=d_tr_insc(i,k,it)+d_tr_bcscav(i,k,it)+d_tr_evapls(i,k,it)
932                   tr_seri(i,k,it)=tr_seri(i,k,it)+d_tr_ls(i,k,it)
933                ENDDO
934             ENDDO
935             CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'lsc scav it = '//TRIM(int2str(it)))
936             ENDIF
937
938          ENDDO  !tr
939
940#ifdef CPP_StratAer
941         IF (type_trac=='coag') THEN
942           ! compute wet deposition flux of sulfur (sum over gases and
943           ! particles) and convert to kg(S)/m2/s
944           ! adding contribution of d_tr_ls to d_tr_cv (above)
945           DO i = 1, klon
946           DO k = 1, klev
947           DO it = 1, nbtr
948             IF (it==id_OCS_strat) THEN
949               budg_dep_wet_ocs(i)=budg_dep_wet_ocs(i)+d_tr_ls(i,k,it)*(mSatom/mOCSmol) &
950                              & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
951             ELSEIF (it==id_SO2_strat) THEN
952               budg_dep_wet_so2(i)=budg_dep_wet_so2(i)+d_tr_ls(i,k,it)*(mSatom/mSO2mol) &
953                              & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
954             ELSEIF (it==id_H2SO4_strat) THEN
955               budg_dep_wet_h2so4(i)=budg_dep_wet_h2so4(i)+d_tr_ls(i,k,it)*(mSatom/mH2SO4mol) &
956                              & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
957             ELSEIF (it.GT.nbtr_sulgas) THEN
958               budg_dep_wet_part(i)=budg_dep_wet_part(i)+d_tr_ls(i,k,it)*(mSatom/mH2SO4mol)  &
959                              & *dens_aer_dry*4./3.*RPI*(mdw(it-nbtr_sulgas)/2.)**3 &
960                              & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
961             ENDIF
962           ENDDO
963           ENDDO
964           ENDDO
965         ENDIF
966#endif
967
968       ELSE IF (iflag_lscav .EQ. 2) THEN ! frac_impa, frac_nucl
969          ! *********   modified  old version
970
971          d_tr_lessi_nucl(:,:,:) = 0.
972          d_tr_lessi_impa(:,:,:) = 0.
973          flestottr(:,:,:) = 0.
974          ! Tendance des aerosols nuclees et impactes
975          DO it = 1, nbtr
976             IF (aerosol(it)) THEN
977                his_dh(:)=0.
978                DO k = 1, klev
979                   DO i = 1, klon
980                      !PhH
981                      zrho(i,k)=pplay(i,k)/t_seri(i,k)/RD
982                      zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho(i,k)/RG
983                      !
984                   ENDDO
985                ENDDO
986
987                DO k=klev-1, 1, -1
988                   DO i=1, klon
989                      !             d_tr_ls(i,k,it)=tr_seri(i,k,it)*(frac_impa(i,k)*frac_nucl(i,k)-1.)
990                      dx=d_tr_ls(i,k,it)
991                      his_dh(i)=his_dh(i)-dx*zrho(i,k)*zdz(i,k)/pdtphys !  kg/m2/s
992                      evaplsc = prfl(i,k) - prfl(i,k+1) + psfl(i,k) - psfl(i,k+1)
993                      ! Evaporation Partielle -> Liberation Partielle 0.5*evap
994                      IF ( evaplsc .LT.0..and.abs(prfl(i,k+1)+psfl(i,k+1)).gt.1.e-10) THEN
995                         evaplsc = (-evaplsc)/(prfl(i,k+1)+psfl(i,k+1))
996                         ! evaplsc est donc positif, his_dh(i) est positif
997                         !-------------- 
998                         d_tr_evapls(i,k,it)=0.5*evaplsc*(d_tr_lessi_nucl(i,k+1,it) &
999                              +d_tr_lessi_impa(i,k+1,it))
1000                         !-------------   d_tr_evapls(i,k,it)=-0.5*evaplsc*(d_tr_lsc(i,k+1,it))
1001                         beta=0.5*evaplsc
1002                         if ((prfl(i,k)+psfl(i,k)).lt.1.e-10) THEN
1003                            beta=1.0*evaplsc
1004                         endif
1005                         dx=beta*his_dh(i)/zrho(i,k)/zdz(i,k)*pdtphys
1006                         his_dh(i)=(1.-beta)*his_dh(i)   ! tracer from
1007                         d_tr_evapls(i,k,it)=dx
1008                      ENDIF
1009                      d_tr_ls(i,k,it)=tr_seri(i,k,it)*(frac_impa(i,k)*frac_nucl(i,k)-1.) &
1010                           +d_tr_evapls(i,k,it)
1011
1012                      !--------------
1013                      d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +    &
1014                           ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
1015                      d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +    &
1016                           ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
1017                      !
1018                      ! Flux lessivage total
1019                      flestottr(i,k,it) = flestottr(i,k,it) -   &
1020                           ( d_tr_lessi_nucl(i,k,it)   +        &
1021                           d_tr_lessi_impa(i,k,it) ) *          &
1022                           ( paprs(i,k)-paprs(i,k+1) ) /        &
1023                           (RG * pdtphys)
1024                      !! Mise a jour des traceurs due a l'impaction,nucleation
1025                      !                 tr_seri(i,k,it)=tr_seri(i,k,it)*frac_impa(i,k)*frac_nucl(i,k)
1026                      !!  calcul de la tendance liee au lessivage stratiforme
1027                      !                 d_tr_ls(i,k,it)=tr_seri(i,k,it)*&
1028                      !                                (1.-1./(frac_impa(i,k)*frac_nucl(i,k)))
1029                      !--------------
1030                   ENDDO
1031                ENDDO
1032             ENDIF
1033          ENDDO
1034          ! *********   end modified old version
1035
1036       ELSE IF (iflag_lscav .EQ. 1) THEN ! frac_impa, frac_nucl
1037          ! *********    old version
1038
1039          d_tr_lessi_nucl(:,:,:) = 0.
1040          d_tr_lessi_impa(:,:,:) = 0.
1041          flestottr(:,:,:) = 0.
1042          !=========================
1043          ! LESSIVAGE LARGE SCALE :
1044          !=========================
1045
1046          ! Tendance des aerosols nuclees et impactes
1047          ! -----------------------------------------
1048          DO it = 1, nbtr
1049             IF (aerosol(it)) THEN
1050                DO k = 1, klev
1051                   DO i = 1, klon
1052                      d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +    &
1053                           ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
1054                      d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +    &
1055                           ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
1056
1057                      !
1058                      ! Flux lessivage total
1059                      ! ------------------------------------------------------------
1060                      flestottr(i,k,it) = flestottr(i,k,it) -   &
1061                           ( d_tr_lessi_nucl(i,k,it)   +        &
1062                           d_tr_lessi_impa(i,k,it) ) *          &
1063                           ( paprs(i,k)-paprs(i,k+1) ) /        &
1064                           (RG * pdtphys)
1065                      !
1066                      ! Mise a jour des traceurs due a l'impaction,nucleation
1067                      ! ----------------------------------------------------------------------
1068                      tr_seri(i,k,it)=tr_seri(i,k,it)*frac_impa(i,k)*frac_nucl(i,k)
1069                   ENDDO
1070                ENDDO
1071             ENDIF
1072          ENDDO
1073
1074          ! *********   end old version
1075       ENDIF  !  iflag_lscav . EQ. 1, 2, 3 or 4
1076       !
1077    ENDIF !  lessivage
1078
1079
1080    !    -- CHIMIE INCA  config_inca = aero or chem --
1081    IF (ANY(type_trac == ['inca','inco'])) THEN  ! ModThL
1082
1083       CALL tracinca(&
1084            nstep,    julien,   gmtime,         lafin,     &
1085            pdtphys,  t_seri,   paprs,          pplay,     &
1086            pmfu,     upwd,     ftsol,  pctsrf, pphis,     &
1087            pphi,     albsol,   sh,    ch,     rh,        &
1088            cldfra,   rneb,     diafra,         cldliq,    &
1089            itop_con, ibas_con, pmflxr,         pmflxs,    &
1090            prfl,     psfl,     aerosol_couple, flxmass_w, &
1091            tau_aero, piz_aero, cg_aero,        ccm,       &
1092            rfname,                                        &
1093            tr_seri(:,:,1+nqCO2:nbtr),  source(:,1+nqCO2:nbtr))  ! ModThL 
1094    ENDIF
1095
1096  END SUBROUTINE phytrac
1097
1098END MODULE
Note: See TracBrowser for help on using the repository browser.