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

Last change on this file since 2172 was 2171, checked in by acozic, 10 years ago

There are some commits that we must not do just before holiday .... so be back to rev 2168

  • 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.0 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    CHARACTER(len=8),DIMENSION(nbtr) :: solsym
311    !RomP >>>
312    INTEGER,SAVE  :: iflag_lscav_omp,iflag_lscav
313    LOGICAL,SAVE  :: convscav_omp,convscav
314!$OMP THREADPRIVATE(iflag_lscav)
315!$OMP THREADPRIVATE(convscav)
316    !RomP <<<
317    !######################################################################
318    !                    -- INITIALIZATION --
319    !######################################################################
320    IF (debutphy) THEN
321       ALLOCATE(d_tr_cl(klon,klev,nbtr),d_tr_dry(klon,nbtr))
322       ALLOCATE(flux_tr_dry(klon,nbtr),d_tr_dec(klon,klev,nbtr),d_tr_cv(klon,klev,nbtr))
323       ALLOCATE(d_tr_insc(klon,klev,nbtr),d_tr_bcscav(klon,klev,nbtr))
324       ALLOCATE(d_tr_evapls(klon,klev,nbtr),d_tr_ls(klon,klev,nbtr))
325       ALLOCATE(qPrls(klon,nbtr),d_tr_trsp(klon,klev,nbtr))
326       ALLOCATE(d_tr_sscav(klon,klev,nbtr),d_tr_sat(klon,klev,nbtr))
327       ALLOCATE(d_tr_uscav(klon,klev,nbtr),qPr(klon,klev,nbtr),qDi(klon,klev,nbtr))
328       ALLOCATE(qPa(klon,klev,nbtr),qMel(klon,klev,nbtr))
329       ALLOCATE(qTrdi(klon,klev,nbtr),dtrcvMA(klon,klev,nbtr))
330       ALLOCATE(d_tr_th(klon,klev,nbtr))
331       ALLOCATE(d_tr_lessi_impa(klon,klev,nbtr),d_tr_lessi_nucl(klon,klev,nbtr))
332    ENDIF
333
334    DO k=1,klev
335       DO i=1,klon
336          sourceBE(i,k)=0.
337          Mint(i,k)=0.
338          zrho(i,k)=0.
339          zdz(i,k)=0.
340       END DO
341    END DO
342
343    DO it=1, nbtr
344       DO k=1,klev
345          DO i=1,klon
346             d_tr_insc(i,k,it)=0.
347             d_tr_bcscav(i,k,it)=0.
348             d_tr_evapls(i,k,it)=0.
349             d_tr_ls(i,k,it)=0.
350             d_tr_cv(i,k,it)=0.
351             d_tr_cl(i,k,it)=0.
352             d_tr_trsp(i,k,it)=0.
353             d_tr_sscav(i,k,it)=0.
354             d_tr_sat(i,k,it)=0.
355             d_tr_uscav(i,k,it)=0.
356             d_tr_lessi_impa(i,k,it)=0.
357             d_tr_lessi_nucl(i,k,it)=0.
358             qDi(i,k,it)=0.
359             qPr(i,k,it)=0.
360             qPa(i,k,it)=0.
361             qMel(i,k,it)=0.
362             qTrdi(i,k,it)=0.
363             dtrcvMA(i,k,it)=0.
364             zmfd1a(i,k,it)=0.
365             zmfdam(i,k,it)=0.
366             zmfphi2(i,k,it)=0.
367          END DO
368       END DO
369    END DO
370
371    DO k = 1, klev
372       DO i = 1, klon
373          delp(i,k) = paprs(i,k)-paprs(i,k+1)
374       END DO
375    END DO
376
377    IF (debutphy) THEN
378       !!jyg
379!$OMP BARRIER
380       ecrit_tra=86400. ! frequence de stokage en dur
381       ! obsolete car remplace par des ecritures dans phys_output_write
382       !RomP >>>
383       !
384       !Config Key  = convscav
385       !Config Desc = Convective scavenging switch: 0=off, 1=on.
386       !Config Def  = .false.
387       !Config Help =
388       !
389!$OMP MASTER
390       convscav_omp=.false.
391       call getin('convscav', convscav_omp)
392       iflag_vdf_trac_omp=1
393       call getin('iflag_vdf_trac', iflag_vdf_trac_omp)
394       iflag_con_trac_omp=1
395       call getin('iflag_con_trac', iflag_con_trac_omp)
396       iflag_the_trac_omp=1
397       call getin('iflag_the_trac', iflag_the_trac_omp)
398!$OMP END MASTER
399!$OMP BARRIER
400       convscav=convscav_omp
401       iflag_vdf_trac=iflag_vdf_trac_omp
402       iflag_con_trac=iflag_con_trac_omp
403       iflag_the_trac=iflag_the_trac_omp
404       write(lunout,*) 'phytrac passage dans routine conv avec lessivage', convscav
405       !
406       !Config Key  = iflag_lscav
407       !Config Desc = Large scale scavenging parametrization: 0=none, 1=old(Genthon92),
408       !              2=1+PHeinrich, 3=Reddy_Boucher2004, 4=3+RPilon.
409       !Config Def  = 1
410       !Config Help =
411       !
412!$OMP MASTER
413       iflag_lscav_omp=1
414       call getin('iflag_lscav', iflag_lscav_omp)
415!$OMP END MASTER
416!$OMP BARRIER
417       iflag_lscav=iflag_lscav_omp
418       !
419       SELECT CASE(iflag_lscav)
420       CASE(0)
421          WRITE(lunout,*)  'Large scale scavenging: none'
422       CASE(1)
423          WRITE(lunout,*)  'Large scale scavenging: C. Genthon, Tellus(1992), 44B, 371-389'
424       CASE(2)
425          WRITE(lunout,*)  'Large scale scavenging: C. Genthon, modified P. Heinrich'
426       CASE(3)
427          WRITE(lunout,*)  'Large scale scavenging: M. Shekkar Reddy and O. Boucher, JGR(2004), 109, D14202'
428       CASE(4)
429          WRITE(lunout,*)  'Large scale scavenging: Reddy and Boucher, modified R. Pilon'
430       END SELECT
431       !RomP <<<
432       WRITE(*,*) 'FIRST TIME IN PHYTRAC : pdtphys(sec) = ',pdtphys,'ecrit_tra (sec) = ',ecrit_tra
433       ALLOCATE( source(klon,nbtr), stat=ierr)
434       IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 1',1)
435
436       ALLOCATE( aerosol(nbtr), stat=ierr)
437       IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 2',1)
438
439
440       ! Initialize module for specific tracers
441       SELECT CASE(type_trac)
442       CASE('lmdz')
443          CALL traclmdz_init(pctsrf,xlat,xlon,ftsol,tr_seri,t_seri,pplay,sh,pdtphys,aerosol,lessivage)
444       CASE('inca')
445          source(:,:)=0.
446          CALL tracinca_init(aerosol,lessivage)
447       CASE('repr')
448          source(:,:)=0.
449       END SELECT
450
451       !
452       !--initialising coefficients for scavenging in the case of NP
453       !
454       ALLOCATE(flag_cvltr(nbtr))
455       IF (iflag_con.EQ.3) THEN
456          !
457          ALLOCATE(ccntrAA(nbtr))
458          ALLOCATE(ccntrENV(nbtr))
459          ALLOCATE(coefcoli(nbtr))
460          !
461          DO it=1, nbtr
462             SELECT CASE(type_trac)
463             CASE('lmdz')
464                IF (convscav.and.aerosol(it)) THEN
465                   flag_cvltr(it)=.true.
466                   ccntrAA(it) =1.0         !--a modifier par JYG a lire depuis fichier
467                   ccntrENV(it)=1.0
468                   coefcoli(it)=0.001
469                ELSE
470                   flag_cvltr(it)=.false.
471                ENDIF
472
473             CASE('inca')
474!                IF ((it.EQ.id_Rn222) .OR. ((it.GE.id_SO2) .AND. (it.LE.id_NH3)) ) THEN
475!                   !--gas-phase species
476!                   flag_cvltr(it)=.false.
477!
478!                ELSEIF ( (it.GE.id_CIDUSTM) .AND. (it.LE.id_AIN) ) THEN
479!                   !--insoluble aerosol species
480!                   flag_cvltr(it)=.true.
481!                   ccntrAA(it)=0.7
482!                   ccntrENV(it)=0.7
483!                   coefcoli(it)=0.001
484!                ELSEIF ( (it.EQ.id_Pb210) .OR. ((it.GE.id_CSSSM) .AND. (it.LE.id_SSN))) THEN
485!                   !--soluble aerosol species
486!                   flag_cvltr(it)=.true.
487!                   ccntrAA(it)=0.9
488!                   ccntrENV(it)=0.9
489!                   coefcoli(it)=0.001
490!                ELSE
491!                   WRITE(lunout,*) 'pb it=', it
492!                   CALL abort_gcm('phytrac','pb it scavenging',1)
493!                ENDIF
494                !--test OB
495                !--for now we do not scavenge in cvltr
496                flag_cvltr(it)=.false.
497             END SELECT
498          ENDDO
499          !
500       ELSE ! iflag_con .ne. 3
501          flag_cvltr(:) = .false.
502       ENDIF
503       !
504       ! Initialize diagnostic output
505       ! ----------------------------
506#ifdef CPP_IOIPSL
507       !     INCLUDE "ini_histrac.h"
508#endif
509       !
510       ! print out all tracer flags
511       !
512       WRITE(lunout,*) 'print out all tracer flags'
513       WRITE(lunout,*) 'type_trac      =', type_trac
514       WRITE(lunout,*) 'config_inca    =', config_inca
515       WRITE(lunout,*) 'iflag_con_trac =', iflag_con_trac
516       WRITE(lunout,*) 'iflag_con      =', iflag_con
517       WRITE(lunout,*) 'convscav       =', convscav
518       WRITE(lunout,*) 'iflag_lscav    =', iflag_lscav
519       WRITE(lunout,*) 'aerosol        =', aerosol
520       WRITE(lunout,*) 'iflag_the_trac =', iflag_the_trac
521       WRITE(lunout,*) 'iflag_thermals =', iflag_thermals
522       WRITE(lunout,*) 'iflag_vdf_trac =', iflag_vdf_trac
523       WRITE(lunout,*) 'pbl_flg        =', pbl_flg
524       WRITE(lunout,*) 'lessivage      =', lessivage
525       write(lunout,*)  'flag_cvltr    = ', flag_cvltr
526
527       IF (lessivage.AND.config_inca.EQ.'inca') THEN
528          CALL abort_gcm('phytrac', 'lessivage=T config_inca=inca impossible',1)
529          STOP
530       ENDIF
531       !
532    END IF ! of IF (debutphy)
533    !############################################ END INITIALIZATION #######
534
535    DO k=1,klev
536       DO i=1,klon
537          zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
538       END DO
539    END DO
540    !
541    IF (id_be .GT. 0) THEN
542       DO k=1,klev
543          DO i=1,klon
544             sourceBE(i,k)=srcbe(i,k)       !RomP  -> pour sortie histrac
545          END DO
546       END DO
547    ENDIF
548
549    !===============================================================================
550    !    -- Do specific treatment according to chemestry model or local LMDZ tracers
551    !     
552    !===============================================================================
553    SELECT CASE(type_trac)
554    CASE('lmdz')
555       !    -- Traitement des traceurs avec traclmdz
556       CALL traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
557            cdragh,  coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon,iflag_vdf_trac>=0,sh, &
558            rh, pphi, ustar, wstar, ale_bl, ale_wake,  u10m, v10m, &
559            tr_seri, source, solsym, d_tr_cl,d_tr_dec, zmasse)               !RomP
560
561    CASE('inca')
562       !    -- CHIMIE INCA  config_inca = aero or chem --
563
564       CALL tracinca(&
565            nstep,    julien,   gmtime,         lafin,     &
566            pdtphys,  t_seri,   paprs,          pplay,     &
567            pmfu,     upwd,     ftsol,  pctsrf, pphis,     &
568            pphi,     albsol,   sh,             rh,        &
569            cldfra,   rneb,     diafra,         cldliq,    &
570            itop_con, ibas_con, pmflxr,         pmflxs,    &
571            prfl,     psfl,     aerosol_couple, flxmass_w, &
572            tau_aero, piz_aero, cg_aero,        ccm,       &
573            rfname,                                        &
574            tr_seri,  source,   solsym)     
575
576    CASE('repr')
577       !   -- CHIMIE REPROBUS --
578
579       CALL tracreprobus(pdtphys, gmtime, debutphy, julien, &
580            presnivs, xlat, xlon, pphis, pphi, &
581            t_seri, pplay, paprs, sh , &
582            tr_seri, solsym)
583
584    END SELECT
585    !======================================================================
586    !       -- Calcul de l'effet de la convection --
587    !======================================================================
588
589    IF (iflag_con_trac==1) THEN
590
591       DO it=1, nbtr
592          IF ( conv_flg(it) == 0 ) CYCLE
593          IF (iflag_con.LT.2) THEN
594             !--pas de transport convectif
595
596             d_tr_cv(:,:,it)=0.
597          ELSE IF (iflag_con.EQ.2) THEN
598             !--ancien transport convectif de Tiedtke
599
600             CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
601                  pplay, paprs, tr_seri(:,:,it), d_tr_cv(:,:,it))
602          ELSE   
603             !--nouveau transport convectif de Emanuel
604
605             IF (flag_cvltr(it)) THEN
606                !--nouveau transport convectif de Emanuel avec lessivage convectif
607                !
608                !
609                ccntrAA_3d(:,:) =ccntrAA(it)
610                ccntrENV_3d(:,:)=ccntrENV(it)
611                coefcoli_3d(:,:)=coefcoli(it)
612
613                !--beware this interface is a bit weird because it is called for each tracer
614                !--with the full array tr_seri even if only item it is processed
615
616                CALL cvltr_scav(pdtphys, da, phi,phi2,d1a,dam, mp,ep,         &
617                     sigd,sij,wght_cvfd,clw,elij,epmlmMm,eplaMm,              &     
618                     pmflxr,pmflxs,evap,t_seri,wdtrainA,wdtrainM,             &   
619                     paprs,it,tr_seri,upwd,dnwd,itop_con,ibas_con,            &
620                     ccntrAA_3d,ccntrENV_3d,coefcoli_3d,                      &
621                     d_tr_cv,d_tr_trsp,d_tr_sscav,d_tr_sat,d_tr_uscav,qDi,qPr,&
622                     qPa,qMel,qTrdi,dtrcvMA,Mint,                             &
623                     zmfd1a,zmfphi2,zmfdam)
624
625
626             ELSE  !---flag_cvltr(it).EQ.FALSE
627                !--nouveau transport convectif de Emanuel mais pas de lessivage convectif
628
629                !--beware this interface is a bit weird because it is called for each tracer
630                !--with the full array tr_seri even if only item it is processed
631                !
632                CALL cvltr_noscav(it,pdtphys, da, phi,mp,wght_cvfd,paprs,pplay, &  !jyg
633                     tr_seri,upwd,dnwd,d_tr_cv)                                !jyg
634
635             ENDIF
636
637          ENDIF !--iflag
638
639          !--on ajoute les tendances
640
641          DO k = 1, klev
642             DO i = 1, klon       
643                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
644             END DO
645          END DO
646
647          CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'convection it = '//solsym(it))
648
649       END DO ! nbtr
650
651    END IF ! convection
652
653    !======================================================================
654    !    -- Calcul de l'effet des thermiques --
655    !======================================================================
656
657    DO it=1,nbtr
658       DO k=1,klev
659          DO i=1,klon
660             d_tr_th(i,k,it)=0.
661             tr_seri(i,k,it)=MAX(tr_seri(i,k,it),0.)
662             tr_seri(i,k,it)=MIN(tr_seri(i,k,it),1.e10)
663          END DO
664       END DO
665    END DO
666
667    IF (iflag_thermals.GT.0.AND.iflag_the_trac>0) THEN   
668
669       DO it=1, nbtr
670
671          CALL thermcell_dq(klon,klev,1,pdtphys,fm_therm,entr_therm, &
672               zmasse,tr_seri(1:klon,1:klev,it),        &
673               d_tr_th(1:klon,1:klev,it),ztra_th,0 )
674
675          DO k=1,klev
676             DO i=1,klon
677                d_tr_th(i,k,it)=pdtphys*d_tr_th(i,k,it)
678                tr_seri(i,k,it)=MAX(tr_seri(i,k,it)+d_tr_th(i,k,it),0.)
679             END DO
680          END DO
681
682       END DO ! it
683
684    END IF ! Thermiques
685
686    !======================================================================
687    !     -- Calcul de l'effet de la couche limite --
688    !======================================================================
689
690    IF (iflag_vdf_trac==1) THEN
691
692       !  Injection during BL mixing
693       !
694       DO it=1, nbtr
695          !
696          IF( pbl_flg(it) /= 0 ) THEN
697             !
698             CALL cltrac(pdtphys, coefh,t_seri,       &
699                  tr_seri(:,:,it), source(:,it),      &
700                  paprs, pplay, delp,                 &
701                  d_tr_cl(:,:,it),d_tr_dry(:,it),flux_tr_dry(:,it))
702             !
703             tr_seri(:,:,it)=tr_seri(:,:,it)+d_tr_cl(:,:,it)
704             !
705          END IF
706          !
707       END DO
708       !
709    ELSE IF (iflag_vdf_trac==0) THEN
710       !
711       !   Injection of source in the first model layer
712       !
713       DO it=1,nbtr
714          d_tr_cl(:,1,it)=source(:,it)*RG/delp(:,1)*pdtphys
715          tr_seri(:,1,it)=tr_seri(:,1,it)+d_tr_cl(:,1,it)
716       ENDDO
717       d_tr_cl(:,2:klev,1:nbtr)=0.
718       !
719    ELSE IF (iflag_vdf_trac==-1) THEN
720       !
721       ! Nothing happens
722       !
723       d_tr_cl=0.
724       !
725    ELSE
726       !
727       CALL abort_gcm('iflag_vdf_trac', 'cas non prevu',1)
728       !
729    END IF ! couche limite
730
731    !======================================================================
732    !   Calcul de l'effet de la precipitation grande echelle
733    !   POUR INCA le lessivage est fait directement dans INCA
734    !======================================================================
735
736    IF (lessivage) THEN
737
738       ql_incloud_ref = 10.e-4
739       ql_incloud_ref =  5.e-4
740
741
742       ! calcul du contenu en eau liquide au sein du nuage
743       ql_incl = ql_incloud_ref
744       ! choix du lessivage
745       !
746       IF (iflag_lscav .EQ. 3 .OR. iflag_lscav .EQ. 4) THEN
747          ! ********  Olivier Boucher version (3) possibly with modified ql_incl (4)
748          !
749          DO it = 1, nbtr
750             !  incloud scavenging and removal by large scale rain ! orig : ql_incl was replaced by 0.5e-3 kg/kg
751             ! the value 0.5e-3 kg/kg is from Giorgi and Chameides (1986), JGR
752             ! Liu (2001) proposed to use 1.5e-3 kg/kg
753
754             CALL lsc_scav(pdtphys,it,iflag_lscav,ql_incl,prfl,psfl,rneb,beta_fisrt,  &
755                  beta_v1,pplay,paprs,t_seri,tr_seri,d_tr_insc,   &
756                  d_tr_bcscav,d_tr_evapls,qPrls)
757
758             !large scale scavenging tendency
759             DO k = 1, klev
760                DO i = 1, klon
761                   d_tr_ls(i,k,it)=d_tr_insc(i,k,it)+d_tr_bcscav(i,k,it)+d_tr_evapls(i,k,it)
762                   tr_seri(i,k,it)=tr_seri(i,k,it)+d_tr_ls(i,k,it)
763                ENDDO
764             ENDDO
765             CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'lsc scav it = '//solsym(it))
766          END DO  !tr
767
768       ELSE IF (iflag_lscav .EQ. 2) THEN ! frac_impa, frac_nucl
769          ! *********   modified  old version
770
771          d_tr_lessi_nucl(:,:,:) = 0.
772          d_tr_lessi_impa(:,:,:) = 0.
773          flestottr(:,:,:) = 0.
774          ! Tendance des aerosols nuclees et impactes
775          DO it = 1, nbtr
776             IF (aerosol(it)) THEN
777                his_dh(:)=0.
778                DO k = 1, klev
779                   DO i = 1, klon
780                      !PhH
781                      zrho(i,k)=pplay(i,k)/t_seri(i,k)/RD
782                      zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho(i,k)/RG
783                      !
784                   END DO
785                END DO
786
787                DO k=klev-1, 1, -1
788                   DO i=1, klon
789                      !             d_tr_ls(i,k,it)=tr_seri(i,k,it)*(frac_impa(i,k)*frac_nucl(i,k)-1.)
790                      dx=d_tr_ls(i,k,it)
791                      his_dh(i)=his_dh(i)-dx*zrho(i,k)*zdz(i,k)/pdtphys !  kg/m2/s
792                      evaplsc = prfl(i,k) - prfl(i,k+1) + psfl(i,k) - psfl(i,k+1)
793                      ! Evaporation Partielle -> Liberation Partielle 0.5*evap
794                      IF ( evaplsc .LT.0..and.abs(prfl(i,k+1)+psfl(i,k+1)).gt.1.e-10) THEN
795                         evaplsc = (-evaplsc)/(prfl(i,k+1)+psfl(i,k+1))
796                         ! evaplsc est donc positif, his_dh(i) est positif
797                         !-------------- 
798                         d_tr_evapls(i,k,it)=0.5*evaplsc*(d_tr_lessi_nucl(i,k+1,it) &
799                              +d_tr_lessi_impa(i,k+1,it))
800                         !-------------   d_tr_evapls(i,k,it)=-0.5*evaplsc*(d_tr_lsc(i,k+1,it))
801                         beta=0.5*evaplsc
802                         if ((prfl(i,k)+psfl(i,k)).lt.1.e-10) THEN
803                            beta=1.0*evaplsc
804                         endif
805                         dx=beta*his_dh(i)/zrho(i,k)/zdz(i,k)*pdtphys
806                         his_dh(i)=(1.-beta)*his_dh(i)   ! tracer from
807                         d_tr_evapls(i,k,it)=dx
808                      ENDIF
809                      d_tr_ls(i,k,it)=tr_seri(i,k,it)*(frac_impa(i,k)*frac_nucl(i,k)-1.) &
810                           +d_tr_evapls(i,k,it)
811
812                      !--------------
813                      d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +    &
814                           ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
815                      d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +    &
816                           ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
817                      !
818                      ! Flux lessivage total
819                      flestottr(i,k,it) = flestottr(i,k,it) -   &
820                           ( d_tr_lessi_nucl(i,k,it)   +        &
821                           d_tr_lessi_impa(i,k,it) ) *          &
822                           ( paprs(i,k)-paprs(i,k+1) ) /        &
823                           (RG * pdtphys)
824                      !! Mise a jour des traceurs due a l'impaction,nucleation
825                      !                 tr_seri(i,k,it)=tr_seri(i,k,it)*frac_impa(i,k)*frac_nucl(i,k)
826                      !!  calcul de la tendance liee au lessivage stratiforme
827                      !                 d_tr_ls(i,k,it)=tr_seri(i,k,it)*&
828                      !                                (1.-1./(frac_impa(i,k)*frac_nucl(i,k)))
829                      !--------------
830                   END DO
831                END DO
832             END IF
833          END DO
834          ! *********   end modified old version
835
836       ELSE IF (iflag_lscav .EQ. 1) THEN ! frac_impa, frac_nucl
837          ! *********    old version
838
839          d_tr_lessi_nucl(:,:,:) = 0.
840          d_tr_lessi_impa(:,:,:) = 0.
841          flestottr(:,:,:) = 0.
842          !=========================
843          ! LESSIVAGE LARGE SCALE :
844          !=========================
845
846          ! Tendance des aerosols nuclees et impactes
847          ! -----------------------------------------
848          DO it = 1, nbtr
849             IF (aerosol(it)) THEN
850                DO k = 1, klev
851                   DO i = 1, klon
852                      d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +    &
853                           ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
854                      d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +    &
855                           ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
856
857                      !
858                      ! Flux lessivage total
859                      ! ------------------------------------------------------------
860                      flestottr(i,k,it) = flestottr(i,k,it) -   &
861                           ( d_tr_lessi_nucl(i,k,it)   +        &
862                           d_tr_lessi_impa(i,k,it) ) *          &
863                           ( paprs(i,k)-paprs(i,k+1) ) /        &
864                           (RG * pdtphys)
865                      !
866                      ! Mise a jour des traceurs due a l'impaction,nucleation
867                      ! ----------------------------------------------------------------------
868                      tr_seri(i,k,it)=tr_seri(i,k,it)*frac_impa(i,k)*frac_nucl(i,k)
869                   END DO
870                END DO
871             END IF
872          END DO
873
874          ! *********   end old version
875       ENDIF  !  iflag_lscav . EQ. 1, 2, 3 or 4
876       !
877    END IF !  lessivage
878
879    !=============================================================
880    !   Ecriture des sorties
881    !=============================================================
882#ifdef CPP_IOIPSL
883    ! INCLUDE "write_histrac.h"
884#endif
885
886  END SUBROUTINE phytrac
887
888END MODULE
Note: See TracBrowser for help on using the repository browser.