source: LMDZ6/branches/Ocean_skin/libf/phylmd/phytrac_mod.F90 @ 3458

Last change on this file since 3458 was 3418, checked in by acozic, 6 years ago

Add modification to take into account value read in INCA restart file

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