source: LMDZ6/trunk/libf/phylmd/phytrac_mod.f90 @ 5473

Last change on this file since 5473 was 5473, checked in by jyg, 24 hours ago

output the total wet deposit of tracers (convective + large scale)

  • 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: 46.0 KB
Line 
1!$Id: phytrac_mod.f90 5473 2025-01-13 10:22:23Z jyg $
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   :: flux_tr_wet    ! tracer wet deposit (surface)                    jyg
38  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE   :: flux_tr_wet_cv ! tracer convective wet deposit (surface)         jyg
39  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: qPr,qDi ! concentration tra dans pluie,air descente insaturee
40  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: qPa,qMel
41  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: qTrdi,dtrcvMA ! conc traceur descente air insaturee et td convective MA
42! RomP <<<
43  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_th  ! Td thermique
44  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_lessi_impa ! Td du lessivage par impaction
45  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_lessi_nucl ! Td du lessivage par nucleation
46  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE   :: qPrls      !jyg: concentration tra dans pluie LS a la surf.
47  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE   :: d_tr_dry ! Td depot sec/traceur (1st layer),ALLOCATABLE,SAVE  jyg
48  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE   :: flux_tr_dry ! depot sec/traceur (surface),ALLOCATABLE,SAVE    jyg
49
50!$OMP THREADPRIVATE(qPa,qMel,qTrdi,dtrcvMA,d_tr_th,d_tr_lessi_impa,d_tr_lessi_nucl)
51!$OMP THREADPRIVATE(d_tr_trsp,d_tr_sscav,d_tr_sat,d_tr_uscav)
52!$OMP THREADPRIVATE(flux_tr_wet,flux_tr_wet_cv,qPr,qDi)
53!$OMP THREADPRIVATE(d_tr_insc,d_tr_bcscav,d_tr_evapls,d_tr_ls,qPrls)
54!$OMP THREADPRIVATE(d_tr_cl,d_tr_dry,flux_tr_dry,d_tr_dec,d_tr_cv)
55
56CONTAINS
57
58  SUBROUTINE phytrac_init()
59
60    USE dimphy
61    USE infotrac_phy, ONLY: nbtr, type_trac
62    USE tracco2i_mod, ONLY: tracco2i_init
63
64    IMPLICIT NONE
65
66    ALLOCATE(d_tr_cl(klon,klev,nbtr),d_tr_dry(klon,nbtr))
67    ALLOCATE(flux_tr_dry(klon,nbtr),d_tr_dec(klon,klev,nbtr),d_tr_cv(klon,klev,nbtr))
68    ALLOCATE(d_tr_insc(klon,klev,nbtr),d_tr_bcscav(klon,klev,nbtr))
69    ALLOCATE(d_tr_evapls(klon,klev,nbtr),d_tr_ls(klon,klev,nbtr))
70    ALLOCATE(qPrls(klon,nbtr),d_tr_trsp(klon,klev,nbtr))
71    ALLOCATE(d_tr_sscav(klon,klev,nbtr),d_tr_sat(klon,klev,nbtr))
72    ALLOCATE(d_tr_uscav(klon,klev,nbtr),qPr(klon,klev,nbtr),qDi(klon,klev,nbtr))
73    ALLOCATE(flux_tr_wet(klon,nbtr),flux_tr_wet_cv(klon,nbtr))
74    ALLOCATE(qPa(klon,klev,nbtr),qMel(klon,klev,nbtr))
75    ALLOCATE(qTrdi(klon,klev,nbtr),dtrcvMA(klon,klev,nbtr))
76    ALLOCATE(d_tr_th(klon,klev,nbtr))
77    ALLOCATE(d_tr_lessi_impa(klon,klev,nbtr),d_tr_lessi_nucl(klon,klev,nbtr))
78
79    !===============================================================================
80    !    -- Do specific treatment according to chemestry model or local LMDZ tracers
81    !
82    !===============================================================================
83    !   -- CO2 interactif --
84    IF (ANY(type_trac == ['co2i','inco'])) CALL tracco2i_init()
85
86       !   -- type_trac == 'co2i' ! PC
87       !   -- CO2 interactif --
88       !   -- source is updated with FF and BB emissions
89       !   -- and net fluxes from ocean and orchidee
90       !   -- sign convention : positive into the atmosphere
91
92  END SUBROUTINE phytrac_init
93
94  SUBROUTINE phytrac(                                 &
95       nstep,     julien,   gmtime,   debutphy,       &
96       lafin,     pdtphys,  u, v,     t_seri,         &
97       paprs,     pplay,    pmfu,     pmfd,           &
98       pen_u,     pde_u,    pen_d,    pde_d,          &
99       cdragh,    coefh,    fm_therm, entr_therm,     &
100       yu1,       yv1,      ftsol,    pctsrf,         &
101       ustar,     u10m,      v10m,                    &
102       wstar,     ale_bl,      ale_wake,              &
103       xlat,      xlon,                               &
104       frac_impa,frac_nucl,beta_fisrt,beta_v1,        &
105       presnivs,  pphis,    pphi,     albsol,         &
106       sh,        ch, rh,   cldfra,   rneb,           &
107       diafra,    cldliq,   itop_con, ibas_con,       &
108       pmflxr,    pmflxs,   prfl,     psfl,           &
109       da,        phi,      mp,       upwd,           &
110       phi2,      d1a,      dam,      sij, wght_cvfd, &   ! RomP +RL
111       wdtrainA,  wdtrainM, sigd,     clw, elij,      &   ! RomP
112       evap,      ep,       epmlmMm,  eplaMm,         &   ! RomP
113       dnwd,      aerosol_couple,     flxmass_w,      &
114       tau_aero,  piz_aero,  cg_aero, ccm,            &
115       rfname,                                        &
116       d_tr_dyn,                                      &   ! RomP
117       tr_seri, init_source)
118    !
119    !======================================================================
120    ! Auteur(s) FH
121    ! Objet: Moniteur general des tendances traceurs
122    ! Modification R. Pilon 01 janvier 2012 transport+scavenging KE scheme : cvltr
123    ! Modification R. Pilon 10 octobre 2012 large scale scavenging incloud_scav + bc_scav
124    !======================================================================
125
126    USE ioipsl
127    USE phys_cal_mod, only : hour
128    USE dimphy
129    USE infotrac_phy, ONLY: nbtr, nqCO2, type_trac, conv_flg, pbl_flg
130    USE strings_mod,  ONLY: int2str
131    USE mod_grid_phy_lmdz
132    USE mod_phys_lmdz_para
133    USE iophy
134    USE traclmdz_mod
135    USE tracinca_mod
136    USE tracreprobus_mod
137    USE indice_sol_mod
138    USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
139    USE print_control_mod, ONLY: lunout
140    USE aero_mod, ONLY : naero_grp
141    USE lmdz_thermcell_dq, ONLY : thermcell_dq
142
143    USE tracco2i_mod
144
145    USE traccoag_mod
146    USE phys_local_var_mod, ONLY: mdw
147    USE phys_local_var_mod, ONLY: budg_dep_dry_ocs,   budg_dep_wet_ocs
148    USE phys_local_var_mod, ONLY: budg_dep_dry_so2,   budg_dep_wet_so2
149    USE phys_local_var_mod, ONLY: budg_dep_dry_h2so4, budg_dep_wet_h2so4
150    USE phys_local_var_mod, ONLY: budg_dep_dry_part,  budg_dep_wet_part
151    USE infotrac_phy, ONLY: nbtr_sulgas, id_OCS_strat, id_SO2_strat, id_H2SO4_strat
152    USE strataer_nuc_mod, ONLY : tracstrataer_init
153    USE aerophys
154
155    USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_STRATAER
156
157    USE yomcst_mod_h
158    USE clesphys_mod_h
159IMPLICIT NONE
160
161
162    !==========================================================================
163    !                   -- ARGUMENT DESCRIPTION --
164    !==========================================================================
165
166    ! Input arguments
167    !----------------
168    !Configuration grille,temps:
169    INTEGER,INTENT(IN) :: nstep      ! Appel physique
170    INTEGER,INTENT(IN) :: julien     ! Jour julien
171    REAL,INTENT(IN)    :: gmtime     ! Heure courante
172    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde)
173    LOGICAL,INTENT(IN) :: debutphy   ! le flag de l'initialisation de la physique
174    LOGICAL,INTENT(IN) :: lafin      ! le flag de la fin de la physique
175
176    REAL,DIMENSION(klon),INTENT(IN) :: xlat    ! latitudes pour chaque point
177    REAL,DIMENSION(klon),INTENT(IN) :: xlon    ! longitudes pour chaque point
178    !
179    !Physique:
180    !--------
181    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
182    REAL,DIMENSION(klon,klev),INTENT(IN)   :: u       ! variable not used
183    REAL,DIMENSION(klon,klev),INTENT(IN)   :: v       ! variable not used
184    REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique
185    REAL,DIMENSION(klon,klev),INTENT(IN)   :: rh      ! humidite relative
186    REAL,DIMENSION(klon,klev),INTENT(IN)   :: ch      ! eau liquide (+ glace si le traceur existe)
187    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
188    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
189    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pphi    ! geopotentiel
190    REAL,DIMENSION(klon),INTENT(IN)        :: pphis
191    REAL,DIMENSION(klev),INTENT(IN)        :: presnivs
192    REAL,DIMENSION(klon,klev),INTENT(IN)   :: cldliq  ! eau condensee totale
193    REAL,DIMENSION(klon,klev),INTENT(IN)   :: cldfra  ! fraction nuageuse (tous les nuages)
194    REAL,DIMENSION(klon,klev),INTENT(IN)   :: diafra  ! fraction nuageuse (convection ou stratus artificiels)
195    REAL,DIMENSION(klon,klev),INTENT(IN)   :: rneb    ! fraction nuageuse (grande echelle)
196    !
197    REAL                                   :: ql_incl ! contenu en eau liquide nuageuse dans le nuage ! ql_incl=oliq/rneb
198    REAL,DIMENSION(klon,klev),INTENT(IN)   :: beta_fisrt ! taux de conversion de l'eau cond (de fisrtilp)
199    REAL,DIMENSION(klon,klev),INTENT(out)  :: beta_v1    ! -- (originale version)
200
201    !
202    INTEGER,DIMENSION(klon),INTENT(IN)     :: itop_con
203    INTEGER,DIMENSION(klon),INTENT(IN)     :: ibas_con
204    REAL,DIMENSION(klon),INTENT(IN)        :: albsol  ! albedo surface
205    !
206    !Dynamique
207    !--------
208    REAL,DIMENSION(klon,klev,nbtr),INTENT(IN)    :: d_tr_dyn
209    !
210    !Convection:
211    !----------
212    REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfu  ! flux de masse dans le panache montant
213    REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfd  ! flux de masse dans le panache descendant
214    REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_u ! flux entraine dans le panache montant
215    REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_u ! flux detraine dans le panache montant
216    REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_d ! flux entraine dans le panache descendant
217    REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_d ! flux detraine dans le panache descendant
218
219    !...Tiedke     
220    REAL,DIMENSION(klon,klev+1),INTENT(IN)   :: pmflxr, pmflxs ! Flux precipitant de pluie, neige aux interfaces [convection]
221    REAL,DIMENSION(klon,klev+1),INTENT(IN)   :: prfl, psfl ! Flux precipitant de pluie, neige aux interfaces [large-scale]
222
223    LOGICAL,INTENT(IN)                       :: aerosol_couple
224    REAL,DIMENSION(klon,klev),INTENT(IN)     :: flxmass_w
225    REAL,DIMENSION(klon,klev,naero_grp,2),INTENT(IN) :: tau_aero
226    REAL,DIMENSION(klon,klev,naero_grp,2),INTENT(IN) :: piz_aero
227    REAL,DIMENSION(klon,klev,naero_grp,2),INTENT(IN) :: cg_aero
228    CHARACTER(len=4),DIMENSION(naero_grp),INTENT(IN) :: rfname
229    REAL,DIMENSION(klon,klev,2),INTENT(IN)   :: ccm
230    !... K.Emanuel
231    REAL,DIMENSION(klon,klev),INTENT(IN)     :: da
232    REAL,DIMENSION(klon,klev,klev),INTENT(IN):: phi
233    ! RomP >>>
234    REAL,DIMENSION(klon,klev),INTENT(IN)      :: d1a,dam
235    REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi2
236    !
237    REAL,DIMENSION(klon,klev),INTENT(IN)      :: wdtrainA
238    REAL,DIMENSION(klon,klev),INTENT(IN)      :: wdtrainM
239    REAL,DIMENSION(klon),INTENT(IN)           :: sigd
240    ! ---- RomP flux entraine, detraine et precipitant kerry Emanuel
241    REAL,DIMENSION(klon,klev),INTENT(IN)      :: evap
242    REAL,DIMENSION(klon,klev),INTENT(IN)      :: ep
243    REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: sij
244    REAL,DIMENSION(klon,klev),INTENT(IN)      :: wght_cvfd          !RL
245    REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: elij
246    REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: epmlmMm
247    REAL,DIMENSION(klon,klev),INTENT(IN)      :: eplaMm
248    REAL,DIMENSION(klon,klev),INTENT(IN)      :: clw
249    ! RomP <<<
250    !
251    REAL,DIMENSION(klon,klev),INTENT(IN)     :: mp
252    REAL,DIMENSION(klon,klev),INTENT(IN)     :: upwd      ! saturated updraft mass flux
253    REAL,DIMENSION(klon,klev),INTENT(IN)     :: dnwd      ! saturated downdraft mass flux
254    !
255    !Thermiques:
256    !----------
257    REAL,DIMENSION(klon,klev+1),INTENT(IN)   :: fm_therm
258    REAL,DIMENSION(klon,klev),INTENT(INOUT)     :: entr_therm
259    !
260    !Couche limite:
261    !--------------
262    !
263    REAL,DIMENSION(:),INTENT(IN)   :: cdragh          ! (klon) coeff drag pour T et Q
264    REAL,DIMENSION(:,:),INTENT(IN) :: coefh           ! (klon,klev) coeff melange CL (m**2/s)
265    REAL,DIMENSION(:),INTENT(IN)   :: ustar,u10m,v10m ! (klon) u* & vent a 10m (m/s)
266    REAL,DIMENSION(:),INTENT(IN)   :: wstar,ale_bl,ale_wake ! (klon) w* and Avail. Lifting Ener.
267    REAL,DIMENSION(:),INTENT(IN)   :: yu1             ! (klon) vents au premier niveau
268    REAL,DIMENSION(:),INTENT(IN)   :: yv1             ! (klon) vents au premier niveau
269    !
270    !Lessivage:
271    !----------
272    !
273    REAL, DIMENSION(:), ALLOCATABLE, SAVE :: ccntrAA
274    REAL, DIMENSION(:), ALLOCATABLE, SAVE :: ccntrENV
275    REAL, DIMENSION(:), ALLOCATABLE, SAVE :: coefcoli
276    LOGICAL, DIMENSION(:), ALLOCATABLE, SAVE :: flag_cvltr
277!$OMP THREADPRIVATE(ccntrAA,ccntrENV,coefcoli,flag_cvltr)
278    REAL, DIMENSION(klon,klev) :: ccntrAA_3d
279    REAL, DIMENSION(klon,klev) :: ccntrENV_3d
280    REAL, DIMENSION(klon,klev) :: coefcoli_3d
281    !
282    ! pour le ON-LINE
283    !
284    REAL,DIMENSION(klon,klev),INTENT(IN) :: frac_impa ! fraction d'aerosols non impactes
285    REAL,DIMENSION(klon,klev),INTENT(IN) :: frac_nucl ! fraction d'aerosols non nuclees
286
287    ! Arguments necessaires pour les sources et puits de traceur:
288    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol  ! Temperature du sol (surf)(Kelvin)
289    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol (nature du sol)
290
291    REAL,DIMENSION(klon)           :: v_dep_dry !dry deposition velocity of stratospheric sulfate in m/s
292    ! Output argument
293    !----------------
294    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
295    REAL,DIMENSION(klon,klev)                    :: sourceBE
296    REAL,DIMENSION(klon,nbtr), INTENT(IN) :: init_source
297
298    !=======================================================================================
299    !                        -- LOCAL VARIABLES --
300    !=======================================================================================
301
302    INTEGER :: i, k, it
303    INTEGER :: nsplit
304
305    !Sources et Reservoirs de traceurs (ex:Radon):
306    !--------------------------------------------
307    !
308    REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: source  ! a voir lorsque le flux de surface est prescrit
309!$OMP THREADPRIVATE(source)
310
311    !
312    !Entrees/Sorties:
313    !---------------
314    INTEGER                   :: iiq, ierr
315    INTEGER                   :: nhori, nvert
316    REAL                      :: zsto, zout, zjulian
317    INTEGER,SAVE              :: nid_tra     ! pointe vers le fichier histrac.nc         
318!$OMP THREADPRIVATE(nid_tra)
319    REAL,DIMENSION(klon)      :: zx_tmp_fi2d ! variable temporaire grille physique
320    INTEGER                   :: itau_w      ! pas de temps ecriture = nstep + itau_phy
321    LOGICAL,PARAMETER         :: ok_sync=.TRUE.
322    !
323    ! Nature du traceur
324    !------------------
325    LOGICAL,DIMENSION(:),ALLOCATABLE,SAVE :: aerosol  ! aerosol(it) = true  => aerosol => lessivage
326!$OMP THREADPRIVATE(aerosol)                        ! aerosol(it) = false => gaz
327    REAL,DIMENSION(klon,klev)             :: delp     ! epaisseur de couche (Pa)
328    !
329    ! Tendances de traceurs (Td) et flux de traceurs:
330    !------------------------
331    REAL,DIMENSION(klon,klev)      :: d_tr     ! Td dans l'atmosphere
332    REAL,DIMENSION(klon,klev)      :: Mint
333    REAL,DIMENSION(klon,klev,nbtr) :: zmfd1a
334    REAL,DIMENSION(klon,klev,nbtr) :: zmfdam
335    REAL,DIMENSION(klon,klev,nbtr) :: zmfphi2
336    ! Physique
337    !----------
338    REAL,DIMENSION(klon,klev,nbtr) :: flestottr ! flux de lessivage dans chaque couche
339    REAL,DIMENSION(klon,klev)      :: zmasse    ! densite atmospherique Kg/m2
340    REAL,DIMENSION(klon,klev)      :: ztra_th
341    !PhH
342    REAL,DIMENSION(klon,klev)      :: zrho
343    REAL,DIMENSION(klon,klev)      :: zdz
344    REAL                           :: evaplsc,dx,beta ! variable pour lessivage Genthon
345    REAL,DIMENSION(klon)           :: his_dh          ! ---
346    ! in-cloud scav variables
347    REAL           :: ql_incloud_ref     ! ref value of in-cloud condensed water content
348
349    !Controles:
350    !---------
351    INTEGER,SAVE :: iflag_vdf_trac,iflag_con_trac,iflag_the_trac
352    INTEGER,SAVE  :: iflag_con_trac_omp, iflag_vdf_trac_omp,iflag_the_trac_omp
353!$OMP THREADPRIVATE(iflag_vdf_trac,iflag_con_trac,iflag_the_trac)
354
355    LOGICAL,SAVE :: lessivage
356!$OMP THREADPRIVATE(lessivage)
357
358    !RomP >>>
359    INTEGER,SAVE  :: iflag_lscav_omp,iflag_lscav
360    REAL, SAVE ::   ccntrAA_in,ccntrAA_omp
361    REAL, SAVE ::   ccntrENV_in,ccntrENV_omp
362    REAL, SAVE ::   coefcoli_in,coefcoli_omp
363
364    LOGICAL,SAVE  :: convscav_omp,convscav
365!$OMP THREADPRIVATE(iflag_lscav)
366!$OMP THREADPRIVATE(ccntrAA_in,ccntrENV_in,coefcoli_in)
367!$OMP THREADPRIVATE(convscav)
368    !RomP <<<
369    !######################################################################
370    !                    -- INITIALIZATION --
371    !######################################################################
372
373    DO k=1,klev
374       DO i=1,klon
375          sourceBE(i,k)=0.
376          Mint(i,k)=0.
377          zrho(i,k)=0.
378          zdz(i,k)=0.
379       ENDDO
380    ENDDO
381
382    DO it=1, nbtr
383       DO k=1,klev
384          DO i=1,klon
385             d_tr_insc(i,k,it)=0.
386             d_tr_bcscav(i,k,it)=0.
387             d_tr_evapls(i,k,it)=0.
388             d_tr_ls(i,k,it)=0.
389             d_tr_cv(i,k,it)=0.
390             d_tr_cl(i,k,it)=0.
391             d_tr_trsp(i,k,it)=0.
392             d_tr_sscav(i,k,it)=0.
393             d_tr_sat(i,k,it)=0.
394             d_tr_uscav(i,k,it)=0.
395             d_tr_lessi_impa(i,k,it)=0.
396             d_tr_lessi_nucl(i,k,it)=0.
397             qDi(i,k,it)=0.
398             qPr(i,k,it)=0.
399             qPa(i,k,it)=0.
400             qMel(i,k,it)=0.
401             qTrdi(i,k,it)=0.
402             dtrcvMA(i,k,it)=0.
403             zmfd1a(i,k,it)=0.
404             zmfdam(i,k,it)=0.
405             zmfphi2(i,k,it)=0.
406          ENDDO
407       ENDDO
408    ENDDO
409
410    DO it=1, nbtr
411       DO i=1,klon
412          d_tr_dry(i,it)=0.
413          flux_tr_dry(i,it)=0.
414          flux_tr_wet(i,it)=0.
415          flux_tr_wet_cv(i,it)=0.
416       ENDDO
417    ENDDO
418
419    DO k = 1, klev
420       DO i = 1, klon
421          delp(i,k) = paprs(i,k)-paprs(i,k+1)
422       ENDDO
423    ENDDO
424
425    IF (debutphy) THEN
426       !!jyg
427!$OMP BARRIER
428       ecrit_tra=86400. ! frequence de stokage en dur
429       ! obsolete car remplace par des ecritures dans phys_output_write
430       !RomP >>>
431       !
432       !Config Key  = convscav
433       !Config Desc = Convective scavenging switch: 0=off, 1=on.
434       !Config Def  = .FALSE.
435       !Config Help =
436       !
437!$OMP MASTER
438       convscav_omp=.FALSE.
439       call getin('convscav', convscav_omp)
440       iflag_vdf_trac_omp=1
441       call getin('iflag_vdf_trac', iflag_vdf_trac_omp)
442       iflag_con_trac_omp=1
443       call getin('iflag_con_trac', iflag_con_trac_omp)
444       iflag_the_trac_omp=1
445       call getin('iflag_the_trac', iflag_the_trac_omp)
446!$OMP END MASTER
447!$OMP BARRIER
448       convscav=convscav_omp
449       iflag_vdf_trac=iflag_vdf_trac_omp
450       iflag_con_trac=iflag_con_trac_omp
451       iflag_the_trac=iflag_the_trac_omp
452       write(lunout,*) 'phytrac passage dans routine conv avec lessivage', convscav
453       !
454       !Config Key  = iflag_lscav
455       !Config Desc = Large scale scavenging parametrization: 0=none, 1=old(Genthon92),
456       !              2=1+PHeinrich, 3=Reddy_Boucher2004, 4=3+RPilon.
457       !Config Def  = 1
458       !Config Help =
459       !
460!$OMP MASTER
461       iflag_lscav_omp=1
462       call getin('iflag_lscav', iflag_lscav_omp)
463       ccntrAA_omp=1
464       ccntrENV_omp=1.
465       coefcoli_omp=0.001
466       call getin('ccntrAA', ccntrAA_omp)
467       call getin('ccntrENV', ccntrENV_omp)
468       call getin('coefcoli', coefcoli_omp)
469!$OMP END MASTER
470!$OMP BARRIER
471       iflag_lscav=iflag_lscav_omp
472       ccntrAA_in=ccntrAA_omp
473       ccntrENV_in=ccntrENV_omp
474       coefcoli_in=coefcoli_omp
475       !
476       SELECT CASE(iflag_lscav)
477       CASE(0)
478          WRITE(lunout,*)  'Large scale scavenging: none'
479       CASE(1)
480          WRITE(lunout,*)  'Large scale scavenging: C. Genthon, Tellus(1992), 44B, 371-389'
481       CASE(2)
482          WRITE(lunout,*)  'Large scale scavenging: C. Genthon, modified P. Heinrich'
483       CASE(3)
484          WRITE(lunout,*)  'Large scale scavenging: M. Shekkar Reddy and O. Boucher, JGR(2004), 109, D14202'
485       CASE(4)
486          WRITE(lunout,*)  'Large scale scavenging: Reddy and Boucher, modified R. Pilon'
487       END SELECT
488       !RomP <<<
489       WRITE(*,*) 'FIRST TIME IN PHYTRAC : pdtphys(sec) = ',pdtphys,'ecrit_tra (sec) = ',ecrit_tra
490       ALLOCATE( source(klon,nbtr), stat=ierr)
491       IF (ierr /= 0) CALL abort_physic('phytrac', 'pb in allocation 1',1)
492
493       ALLOCATE( aerosol(nbtr), stat=ierr)
494       IF (ierr /= 0) CALL abort_physic('phytrac', 'pb in allocation 2',1)
495
496
497       ! Initialize module for specific tracers
498       IF (type_trac == 'inca') THEN
499          source(:,:)=init_source(:,:)
500          CALL tracinca_init(aerosol,lessivage)
501       ELSE IF (type_trac == 'repr') THEN
502          source(:,:)=0.
503       ELSE IF (type_trac == 'co2i') THEN
504          source(:,:)=0.
505          lessivage  = .FALSE.
506          aerosol(:) = .FALSE.
507          pbl_flg(:) = 1
508          iflag_the_trac= 1
509          iflag_vdf_trac= 1
510          iflag_con_trac= 1
511       ELSE IF (type_trac == 'inco') THEN
512          source(:,1:nqCO2) = 0.                          ! from CO2i ModThL
513          source(:,nqCO2+1:nbtr)=init_source(:,:)         ! from INCA ModThL
514          aerosol(1:nqCO2) = .FALSE.                      ! from CO2i ModThL
515          CALL tracinca_init(aerosol(nqCO2+1:nbtr),lessivage)     ! from INCA ModThL
516          pbl_flg(1:nqCO2) = 1                            ! From CO2i ModThL
517          iflag_the_trac = 1                              ! From CO2i
518          iflag_vdf_trac = 1                              ! From CO2i
519          iflag_con_trac = 1                              ! From CO2i
520       ELSE IF (type_trac == 'coag' .AND. CPPKEY_STRATAER) THEN
521          source(:,:)=0.
522          CALL tracstrataer_init(aerosol,lessivage)       ! init aerosols and lessivage param
523       ELSE IF (type_trac == 'lmdz') THEN
524          CALL traclmdz_init(pctsrf,xlat,xlon,ftsol,tr_seri,t_seri,pplay,sh,pdtphys,aerosol,lessivage)
525       ENDIF
526
527       !
528       !--initialising coefficients for scavenging in the case of NP
529       !
530       ALLOCATE(flag_cvltr(nbtr))
531       IF (iflag_con.EQ.3) THEN
532          !
533          ALLOCATE(ccntrAA(nbtr))
534          ALLOCATE(ccntrENV(nbtr))
535          ALLOCATE(coefcoli(nbtr))
536          !
537          DO it=1, nbtr
538             IF (type_trac == 'repr') THEN
539                 flag_cvltr(it)=.FALSE.
540             ELSE IF (type_trac == 'inca') THEN
541!                IF ((it.EQ.id_Rn222) .OR. ((it.GE.id_SO2) .AND. (it.LE.id_NH3)) ) THEN
542!                   !--gas-phase species
543!                   flag_cvltr(it)=.FALSE.
544!
545!                ELSEIF ( (it.GE.id_CIDUSTM) .AND. (it.LE.id_AIN) ) THEN
546!                   !--insoluble aerosol species
547!                   flag_cvltr(it)=.TRUE.
548!                   ccntrAA(it)=0.7
549!                   ccntrENV(it)=0.7
550!                   coefcoli(it)=0.001
551!                ELSEIF ( (it.EQ.id_Pb210) .OR. ((it.GE.id_CSSSM) .AND. (it.LE.id_SSN))) THEN
552!                   !--soluble aerosol species
553!                   flag_cvltr(it)=.TRUE.
554!                   ccntrAA(it)=0.9
555!                   ccntrENV(it)=0.9
556!                   coefcoli(it)=0.001
557!                ELSE
558!                   WRITE(lunout,*) 'pb it=', it
559!                   CALL abort_physic('phytrac','pb it scavenging',1)
560!                ENDIF
561                !--test OB
562                !--for now we do not scavenge in cvltr
563                flag_cvltr(it)=.FALSE.
564             ELSE IF (type_trac == 'co2i') THEN
565                !--co2 tracers are not scavenged
566                flag_cvltr(it)=.FALSE.
567             ELSE IF (type_trac == 'inco') THEN     ! Add ThL
568                flag_cvltr(it)=.FALSE.
569             ELSE IF (type_trac == 'coag' .AND. CPPKEY_STRATAER) THEN
570                IF (convscav.and.aerosol(it)) THEN
571                   flag_cvltr(it)=.TRUE.
572                   ccntrAA(it) =ccntrAA_in
573                   ccntrENV(it)=ccntrENV_in
574                   coefcoli(it)=coefcoli_in
575                ELSE
576                   flag_cvltr(it)=.FALSE.
577                ENDIF
578             ELSE IF (type_trac == 'lmdz') THEN
579                IF (convscav.and.aerosol(it)) THEN
580                   flag_cvltr(it)=.TRUE.
581                   ccntrAA(it) =ccntrAA_in    !--a modifier par JYG a lire depuis fichier
582                   ccntrENV(it)=ccntrENV_in
583                   coefcoli(it)=coefcoli_in
584                ELSE
585                   flag_cvltr(it)=.FALSE.
586                ENDIF
587             ENDIF
588          ENDDO
589          !
590       ELSE ! iflag_con .ne. 3
591          flag_cvltr(:) = .FALSE.
592       ENDIF
593       !
594       ! print out all tracer flags
595       !
596       WRITE(lunout,*) 'print out all tracer flags'
597       WRITE(lunout,*) 'type_trac      =', type_trac
598       WRITE(lunout,*) 'config_inca    =', config_inca
599       WRITE(lunout,*) 'iflag_con_trac =', iflag_con_trac
600       WRITE(lunout,*) 'iflag_con      =', iflag_con
601       WRITE(lunout,*) 'convscav       =', convscav
602       WRITE(lunout,*) 'iflag_lscav    =', iflag_lscav
603       WRITE(lunout,*) 'aerosol        =', aerosol
604       WRITE(lunout,*) 'iflag_the_trac =', iflag_the_trac
605       WRITE(lunout,*) 'iflag_thermals =', iflag_thermals
606       WRITE(lunout,*) 'iflag_vdf_trac =', iflag_vdf_trac
607       WRITE(lunout,*) 'pbl_flg        =', pbl_flg
608       WRITE(lunout,*) 'lessivage      =', lessivage
609       write(lunout,*)  'flag_cvltr    = ', flag_cvltr
610
611       IF (lessivage .AND. ANY(type_trac == ['inca','inco'])) &
612          CALL abort_physic('phytrac', 'lessivage=T config_inca=inca impossible',1)
613       !
614    ENDIF ! of IF (debutphy)
615    !############################################ END INITIALIZATION #######
616
617    DO k=1,klev
618       DO i=1,klon
619          zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
620       ENDDO
621    ENDDO
622    !
623    IF (id_be .GT. 0) THEN
624       DO k=1,klev
625          DO i=1,klon
626             sourceBE(i,k)=srcbe(i,k)       !RomP  -> pour sortie histrac
627          ENDDO
628       ENDDO
629    ENDIF
630
631    !===============================================================================
632    !    -- Do specific treatment according to chemestry model or local LMDZ tracers
633    !     
634    !===============================================================================
635    IF (type_trac == 'inca') THEN
636       !    -- CHIMIE INCA  config_inca = aero or chem --
637       ! Appel fait en fin de phytrac pour avoir les emissions modifiees par
638       ! la couche limite et la convection avant le calcul de la chimie
639
640    ELSE IF (type_trac == 'repr') THEN
641       !   -- CHIMIE REPROBUS --
642       CALL tracreprobus(pdtphys, gmtime, debutphy, julien, &
643            presnivs, xlat, xlon, pphis, pphi, &
644            t_seri, pplay, paprs, sh , &
645            tr_seri)
646
647    ELSE IF (type_trac == 'co2i') THEN
648       !   -- CO2 interactif --
649       !   -- source is updated with FF and BB emissions
650       !   -- and net fluxes from ocean and orchidee
651       !   -- sign convention : positive into the atmosphere
652
653       CALL tracco2i(pdtphys, debutphy, &
654            xlat, xlon, pphis, pphi, &
655            t_seri, pplay, paprs, tr_seri, source)
656    ELSE IF (type_trac == 'inco') THEN      ! Add ThL
657       CALL tracco2i(pdtphys, debutphy, &
658            xlat, xlon, pphis, pphi, &
659            t_seri, pplay, paprs, tr_seri, source)
660
661    ELSE IF (type_trac == 'coag' .AND. CPPKEY_STRATAER) 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    ELSE IF (type_trac == 'lmdz') THEN
668       !    -- Traitement des traceurs avec traclmdz
669       CALL traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
670            cdragh,  coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon,iflag_vdf_trac>=0,sh, &
671            rh, pphi, ustar, wstar, ale_bl, ale_wake,  u10m, v10m, &
672            tr_seri, source, d_tr_cl,d_tr_dec, zmasse)               !RomP
673    ENDIF
674    !======================================================================
675    !       -- Calcul de l'effet de la convection --
676    !======================================================================
677
678    IF (iflag_con_trac==1) THEN
679
680       DO it=1, nbtr
681          IF ( conv_flg(it) == 0 ) CYCLE
682          IF (iflag_con.LT.2) THEN
683             !--pas de transport convectif
684             d_tr_cv(:,:,it)=0.
685
686          ELSE IF (iflag_con.EQ.2) THEN
687             !--ancien transport convectif de Tiedtke
688
689             CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
690                  pplay, paprs, tr_seri(:,:,it), d_tr_cv(:,:,it))
691          ELSE   
692             !--nouveau transport convectif de Emanuel
693
694             IF (flag_cvltr(it)) THEN
695                !--nouveau transport convectif de Emanuel avec lessivage convectif
696                !
697                !
698                ccntrAA_3d(:,:) =ccntrAA(it)
699                ccntrENV_3d(:,:)=ccntrENV(it)
700                coefcoli_3d(:,:)=coefcoli(it)
701
702                !--beware this interface is a bit weird because it is called for each tracer
703                !--with the full array tr_seri even if only item it is processed
704
705                CALL cvltr_scav(pdtphys, da, phi,phi2,d1a,dam, mp,ep,                 &
706                     sigd,sij,wght_cvfd,clw,elij,epmlmMm,eplaMm,                      &   
707                     pmflxr,pmflxs,evap,t_seri,wdtrainA,wdtrainM,                     & 
708                     paprs,it,tr_seri,upwd,dnwd,itop_con,ibas_con,                    &
709                     ccntrAA_3d,ccntrENV_3d,coefcoli_3d,                              &
710                     d_tr_cv,d_tr_trsp,d_tr_sscav,d_tr_sat,d_tr_uscav,flux_tr_wet_cv, &
711                     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
741IF (CPPKEY_STRATAER) THEN
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
772END IF
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       !
822IF (CPPKEY_STRATAER) THEN
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
846END IF
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             !
859IF (CPPKEY_STRATAER) THEN
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
873END IF
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             !total wet deposit = large scale wet deposit + convective wet deposit
929             DO i = 1, klon
930               flux_tr_wet(i, it) = flux_tr_wet_cv(i, it) + &
931                                    qPrls(i, it)*(prfl(i, 1)+psfl(i, 1))*pdtphys
932             ENDDO  ! i = 1, klon
933
934             !large scale scavenging tendency
935             DO k = 1, klev
936                DO i = 1, klon
937                   d_tr_ls(i,k,it)=d_tr_insc(i,k,it)+d_tr_bcscav(i,k,it)+d_tr_evapls(i,k,it)
938                   tr_seri(i,k,it)=tr_seri(i,k,it)+d_tr_ls(i,k,it)
939                ENDDO
940             ENDDO
941             CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'lsc scav it = '//TRIM(int2str(it)))
942             ENDIF
943
944          ENDDO  !tr
945
946IF (CPPKEY_STRATAER) THEN
947         IF (type_trac=='coag') THEN
948           ! compute wet deposition flux of sulfur (sum over gases and
949           ! particles) and convert to kg(S)/m2/s
950           ! adding contribution of d_tr_ls to d_tr_cv (above)
951           DO i = 1, klon
952           DO k = 1, klev
953           DO it = 1, nbtr
954             IF (it==id_OCS_strat) THEN
955               budg_dep_wet_ocs(i)=budg_dep_wet_ocs(i)+d_tr_ls(i,k,it)*(mSatom/mOCSmol) &
956                              & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
957             ELSEIF (it==id_SO2_strat) THEN
958               budg_dep_wet_so2(i)=budg_dep_wet_so2(i)+d_tr_ls(i,k,it)*(mSatom/mSO2mol) &
959                              & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
960             ELSEIF (it==id_H2SO4_strat) THEN
961               budg_dep_wet_h2so4(i)=budg_dep_wet_h2so4(i)+d_tr_ls(i,k,it)*(mSatom/mH2SO4mol) &
962                              & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
963             ELSEIF (it.GT.nbtr_sulgas) THEN
964               budg_dep_wet_part(i)=budg_dep_wet_part(i)+d_tr_ls(i,k,it)*(mSatom/mH2SO4mol)  &
965                              & *dens_aer_dry*4./3.*RPI*(mdw(it-nbtr_sulgas)/2.)**3 &
966                              & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
967             ENDIF
968           ENDDO
969           ENDDO
970           ENDDO
971         ENDIF
972END IF
973
974       ELSE IF (iflag_lscav .EQ. 2) THEN ! frac_impa, frac_nucl
975          ! *********   modified  old version
976
977          d_tr_lessi_nucl(:,:,:) = 0.
978          d_tr_lessi_impa(:,:,:) = 0.
979          flestottr(:,:,:) = 0.
980          ! Tendance des aerosols nuclees et impactes
981          DO it = 1, nbtr
982             IF (aerosol(it)) THEN
983                his_dh(:)=0.
984                DO k = 1, klev
985                   DO i = 1, klon
986                      !PhH
987                      zrho(i,k)=pplay(i,k)/t_seri(i,k)/RD
988                      zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho(i,k)/RG
989                      !
990                   ENDDO
991                ENDDO
992
993                DO k=klev-1, 1, -1
994                   DO i=1, klon
995                      !             d_tr_ls(i,k,it)=tr_seri(i,k,it)*(frac_impa(i,k)*frac_nucl(i,k)-1.)
996                      dx=d_tr_ls(i,k,it)
997                      his_dh(i)=his_dh(i)-dx*zrho(i,k)*zdz(i,k)/pdtphys !  kg/m2/s
998                      evaplsc = prfl(i,k) - prfl(i,k+1) + psfl(i,k) - psfl(i,k+1)
999                      ! Evaporation Partielle -> Liberation Partielle 0.5*evap
1000                      IF ( evaplsc .LT.0..and.abs(prfl(i,k+1)+psfl(i,k+1)).gt.1.e-10) THEN
1001                         evaplsc = (-evaplsc)/(prfl(i,k+1)+psfl(i,k+1))
1002                         ! evaplsc est donc positif, his_dh(i) est positif
1003                         !-------------- 
1004                         d_tr_evapls(i,k,it)=0.5*evaplsc*(d_tr_lessi_nucl(i,k+1,it) &
1005                              +d_tr_lessi_impa(i,k+1,it))
1006                         !-------------   d_tr_evapls(i,k,it)=-0.5*evaplsc*(d_tr_lsc(i,k+1,it))
1007                         beta=0.5*evaplsc
1008                         if ((prfl(i,k)+psfl(i,k)).lt.1.e-10) THEN
1009                            beta=1.0*evaplsc
1010                         endif
1011                         dx=beta*his_dh(i)/zrho(i,k)/zdz(i,k)*pdtphys
1012                         his_dh(i)=(1.-beta)*his_dh(i)   ! tracer from
1013                         d_tr_evapls(i,k,it)=dx
1014                      ENDIF
1015                      d_tr_ls(i,k,it)=tr_seri(i,k,it)*(frac_impa(i,k)*frac_nucl(i,k)-1.) &
1016                           +d_tr_evapls(i,k,it)
1017
1018                      !--------------
1019                      d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +    &
1020                           ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
1021                      d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +    &
1022                           ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
1023                      !
1024                      ! Flux lessivage total
1025                      flestottr(i,k,it) = flestottr(i,k,it) -   &
1026                           ( d_tr_lessi_nucl(i,k,it)   +        &
1027                           d_tr_lessi_impa(i,k,it) ) *          &
1028                           ( paprs(i,k)-paprs(i,k+1) ) /        &
1029                           (RG * pdtphys)
1030                      !! Mise a jour des traceurs due a l'impaction,nucleation
1031                      !                 tr_seri(i,k,it)=tr_seri(i,k,it)*frac_impa(i,k)*frac_nucl(i,k)
1032                      !!  calcul de la tendance liee au lessivage stratiforme
1033                      !                 d_tr_ls(i,k,it)=tr_seri(i,k,it)*&
1034                      !                                (1.-1./(frac_impa(i,k)*frac_nucl(i,k)))
1035                      !--------------
1036                   ENDDO
1037                ENDDO
1038             ENDIF
1039          ENDDO
1040          ! *********   end modified old version
1041
1042       ELSE IF (iflag_lscav .EQ. 1) THEN ! frac_impa, frac_nucl
1043          ! *********    old version
1044
1045          d_tr_lessi_nucl(:,:,:) = 0.
1046          d_tr_lessi_impa(:,:,:) = 0.
1047          flestottr(:,:,:) = 0.
1048          !=========================
1049          ! LESSIVAGE LARGE SCALE :
1050          !=========================
1051
1052          ! Tendance des aerosols nuclees et impactes
1053          ! -----------------------------------------
1054          DO it = 1, nbtr
1055             IF (aerosol(it)) THEN
1056                DO k = 1, klev
1057                   DO i = 1, klon
1058                      d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +    &
1059                           ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
1060                      d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +    &
1061                           ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
1062
1063                      !
1064                      ! Flux lessivage total
1065                      ! ------------------------------------------------------------
1066                      flestottr(i,k,it) = flestottr(i,k,it) -   &
1067                           ( d_tr_lessi_nucl(i,k,it)   +        &
1068                           d_tr_lessi_impa(i,k,it) ) *          &
1069                           ( paprs(i,k)-paprs(i,k+1) ) /        &
1070                           (RG * pdtphys)
1071                      !
1072                      ! Mise a jour des traceurs due a l'impaction,nucleation
1073                      ! ----------------------------------------------------------------------
1074                      tr_seri(i,k,it)=tr_seri(i,k,it)*frac_impa(i,k)*frac_nucl(i,k)
1075                   ENDDO
1076                ENDDO
1077             ENDIF
1078          ENDDO
1079
1080          ! *********   end old version
1081       ENDIF  !  iflag_lscav .EQ. 1, 2, 3 or 4
1082       !
1083    ENDIF !  lessivage
1084
1085
1086    !    -- CHIMIE INCA  config_inca = aero or chem --
1087    IF (ANY(type_trac == ['inca','inco'])) THEN  ! ModThL
1088
1089       CALL tracinca(&
1090            nstep,    julien,   gmtime,         lafin,     &
1091            pdtphys,  t_seri,   paprs,          pplay,     &
1092            pmfu,     upwd,     ftsol,  pctsrf, pphis,     &
1093            pphi,     albsol,   sh,    ch,     rh,        &
1094            cldfra,   rneb,     diafra,         cldliq,    &
1095            itop_con, ibas_con, pmflxr,         pmflxs,    &
1096            prfl,     psfl,     aerosol_couple, flxmass_w, &
1097            tau_aero, piz_aero, cg_aero,        ccm,       &
1098            rfname,                                        &
1099            tr_seri(:,:,1+nqCO2:nbtr),  source(:,1+nqCO2:nbtr))  ! ModThL 
1100    ENDIF
1101
1102  END SUBROUTINE phytrac
1103
1104END MODULE
Note: See TracBrowser for help on using the repository browser.