source: LMDZ5/branches/testing/libf/phylmd/phytrac.F90 @ 1790

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

Version testing basée sur r1745


Testing release based on r1745

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