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

Last change on this file since 1821 was 1820, checked in by Ehouarn Millour, 11 years ago

Correction d'un bug dans phytrac affectant la compilation en openmp.
Attention à l'heure actuelle il n'y a plus convergence entre la version séquentielle et la version omp.
UG
.....................
Bug correction in phytrac, no compiling in openmp as well as in sequential. Be warned that convergence have been lost between sequential and openmp versions.
UG

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