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

Last change on this file since 1893 was 1877, checked in by Laurent Fairhead, 11 years ago

Ménage sur le code pour éliminer les calculs spécifiques ISCCP hors COSP
Un certain nombre de variables non utilisées dans physiq.F90 ont aussi
été supprimée. Environ 840 lignes supprimées du code physiq.F90


Code cleanup to eliminate specific references to ISCCP outside the COSP
library. Unused variables have been cleaned up from the physiq.F90 routine
as well.

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