source: LMDZ5/trunk/libf/phylmd/phytrac_mod.F90 @ 2225

Last change on this file since 2225 was 2210, checked in by fhourdin, 10 years ago

Control of convective scavenging by .def files
Conrôle du lessivage convectif dans les fichiers .def
ccntrAA= 1. # fraction d'aerosol en phase condensee dans l'ascendance adiab
ccntrENV=1. # fraction d'aerosol en phase condensee dans les melange (Emanuel)
coefcoli=1.e-3 # efficacité de collision

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