source: LMDZ4/trunk/libf/phylmd/phytrac.F90 @ 1357

Last change on this file since 1357 was 1309, checked in by Laurent Fairhead, 15 years ago

The histrac history file is not used at the present. It is disabled for the
CMIP5 project


Le fichier histoire histrac n'est pas utilisé pour l'instant. Il est désactivé
pour CMIP5

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.0 KB
Line 
1!$Id $
2
3SUBROUTINE phytrac(                            &
4     nstep,     julien,   gmtime,   debutphy,  &
5     lafin,     pdtphys,  u, v,     t_seri,     &
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,    &
10     xlat,      frac_impa,frac_nucl,xlon,      &
11     presnivs,  pphis,    pphi,     albsol,    &
12     sh,        rh,       cldfra,   rneb,      &
13     diafra,    cldliq,   itop_con, ibas_con,  &
14     pmflxr,    pmflxs,   prfl,     psfl,      &
15     da,        phi,      mp,       upwd,      &
16     dnwd,      aerosol_couple,     flxmass_w, &
17     tau_aero,  piz_aero,  cg_aero, ccm,       &
18     rfname,                                   &
19     tr_seri)         
20!
21!======================================================================
22! Auteur(s) FH
23! Objet: Moniteur general des tendances traceurs
24!======================================================================
25
26  USE ioipsl
27  USE dimphy
28  USE infotrac
29  USE mod_grid_phy_lmdz
30  USE mod_phys_lmdz_para
31  USE comgeomphy
32  USE iophy
33  USE traclmdz_mod
34  USE tracinca_mod
35
36
37  IMPLICIT NONE
38
39  INCLUDE "YOMCST.h"
40  INCLUDE "dimensions.h"
41  INCLUDE "indicesol.h"
42  INCLUDE "clesphys.h"
43  INCLUDE "temps.h"
44  INCLUDE "paramet.h"
45  INCLUDE "control.h"
46  INCLUDE "thermcell.h"
47!==========================================================================
48!                   -- ARGUMENT DESCRIPTION --
49!==========================================================================
50
51! Input arguments
52!----------------
53!Configuration grille,temps:
54  INTEGER,INTENT(IN) :: nstep      ! Appel physique
55  INTEGER,INTENT(IN) :: julien     ! Jour julien
56  REAL,INTENT(IN)    :: gmtime
57  REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde)
58  LOGICAL,INTENT(IN) :: debutphy   ! le flag de l'initialisation de la physique
59  LOGICAL,INTENT(IN) :: lafin      ! le flag de la fin de la physique
60 
61  REAL,DIMENSION(klon),INTENT(IN) :: xlat    ! latitudes pour chaque point
62  REAL,DIMENSION(klon),INTENT(IN) :: xlon    ! longitudes pour chaque point
63!
64!Physique:
65!--------
66  REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
67  REAL,DIMENSION(klon,klev),INTENT(IN)   :: u       !
68  REAL,DIMENSION(klon,klev),INTENT(IN)   :: v       !
69  REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique
70  REAL,DIMENSION(klon,klev),INTENT(IN)   :: rh      ! humidite relative
71  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
72  REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
73  REAL,DIMENSION(klon,klev),INTENT(IN)   :: pphi    ! geopotentiel
74  REAL,DIMENSION(klon),INTENT(IN)        :: pphis
75  REAL,DIMENSION(klev),INTENT(IN)        :: presnivs
76  REAL,DIMENSION(klon,klev),INTENT(IN)   :: cldliq  ! eau liquide nuageuse
77  REAL,DIMENSION(klon,klev),INTENT(IN)   :: cldfra  ! fraction nuageuse (tous les nuages)
78  REAL,DIMENSION(klon,klev),INTENT(IN)   :: diafra  ! fraction nuageuse (convection ou stratus artificiels)
79  REAL,DIMENSION(klon,klev),INTENT(IN)   :: rneb    ! fraction nuageuse (grande echelle)
80  INTEGER,DIMENSION(klon),INTENT(IN)     :: itop_con
81  INTEGER,DIMENSION(klon),INTENT(IN)     :: ibas_con
82  REAL,DIMENSION(klon),INTENT(IN)        :: albsol  ! albedo surface
83!
84!Convection:
85!----------
86  REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfu  ! flux de masse dans le panache montant
87  REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfd  ! flux de masse dans le panache descendant
88  REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_u ! flux entraine dans le panache montant
89  REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_u ! flux detraine dans le panache montant
90  REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_d ! flux entraine dans le panache descendant
91  REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_d ! flux detraine dans le panache descendant
92
93!...Tiedke     
94  REAL,DIMENSION(klon,klev+1),INTENT(IN)   :: pmflxr, pmflxs ! Flux precipitant de pluie, neige aux interfaces [convection]
95  REAL,DIMENSION(klon,klev+1),INTENT(IN)   :: prfl, psfl ! Flux precipitant de pluie, neige aux interfaces [large-scale]
96
97  LOGICAL,INTENT(IN)                       :: aerosol_couple
98  REAL,DIMENSION(klon,klev),INTENT(IN)     :: flxmass_w
99  REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: tau_aero
100  REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: piz_aero
101  REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: cg_aero
102  CHARACTER(len=4),DIMENSION(9),INTENT(IN) :: rfname
103  REAL,DIMENSION(klon,klev,2),INTENT(IN)   :: ccm
104!... K.Emanuel
105  REAL,DIMENSION(klon,klev),INTENT(IN)     :: da
106  REAL,DIMENSION(klon,klev,klev),INTENT(IN):: phi
107  REAL,DIMENSION(klon,klev),INTENT(IN)     :: mp
108  REAL,DIMENSION(klon,klev),INTENT(IN)     :: upwd      ! saturated updraft mass flux
109  REAL,DIMENSION(klon,klev),INTENT(IN)     :: dnwd      ! saturated downdraft mass flux
110!
111!Thermiques:
112!----------
113  REAL,DIMENSION(klon,klev+1),INTENT(IN)   :: fm_therm
114  REAL,DIMENSION(klon,klev),INTENT(IN)     :: entr_therm
115!
116!Couche limite:
117!--------------
118!
119  REAL,DIMENSION(klon,klev),INTENT(IN) :: cdragh ! coeff drag pour T et Q
120  REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh  ! coeff melange CL (m**2/s)
121  REAL,DIMENSION(klon),INTENT(IN)      :: yu1    ! vents au premier niveau
122  REAL,DIMENSION(klon),INTENT(IN)      :: yv1    ! vents au premier niveau
123!
124!Lessivage:
125!----------
126!
127! pour le ON-LINE
128!
129  REAL,DIMENSION(klon,klev),INTENT(IN) :: frac_impa ! fraction d'aerosols non impactes
130  REAL,DIMENSION(klon,klev),INTENT(IN) :: frac_nucl ! fraction d'aerosols non nuclees
131
132! Arguments necessaires pour les sources et puits de traceur:
133  REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol  ! Temperature du sol (surf)(Kelvin)
134  REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol (nature du sol)
135
136
137! Output argument
138!----------------
139  REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] 
140
141!=======================================================================================
142!                        -- LOCAL VARIABLES --
143!=======================================================================================
144
145  INTEGER :: i, k, it
146  INTEGER :: nsplit
147
148!Sources et Reservoirs de traceurs (ex:Radon):
149!--------------------------------------------
150!
151  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: source  ! a voir lorsque le flux de surface est prescrit
152!$OMP THREADPRIVATE(source)
153
154!
155!Entrees/Sorties: (cf ini_histrac.h et write_histrac.h) 
156!---------------
157  INTEGER                   :: iiq, ierr
158  INTEGER                   :: nhori, nvert
159  REAL                      :: zsto, zout, zjulian
160  INTEGER,SAVE              :: nid_tra     ! pointe vers le fichier histrac.nc         
161!$OMP THREADPRIVATE(nid_tra)
162  REAL,DIMENSION(klon)      :: zx_tmp_fi2d ! variable temporaire grille physique
163  INTEGER                   :: itau_w      ! pas de temps ecriture = nstep + itau_phy
164  LOGICAL,PARAMETER :: ok_sync=.TRUE.
165
166!
167! Nature du traceur
168!------------------
169  LOGICAL,DIMENSION(:),ALLOCATABLE,SAVE :: aerosol  ! aerosol(it) = true  => aerosol => lessivage
170!$OMP THREADPRIVATE(aerosol)                        ! aerosol(it) = false => gaz
171  REAL,DIMENSION(klon,klev)             :: delp     ! epaisseur de couche (Pa)
172!
173! Tendances de traceurs (Td):
174!------------------------
175!
176  REAL,DIMENSION(klon,klev)      :: d_tr     ! Td dans l'atmosphere
177  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_cl  ! Td couche limite/traceur
178  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_cv  ! Td convection/traceur
179  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_th  ! Td thermique
180  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_lessi_impa ! Td du lessivage par impaction
181  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_lessi_nucl ! Td du lessivage par nucleation
182!
183! Physique
184!----------   
185  REAL,DIMENSION(klon,klev,nbtr) :: flestottr ! flux de lessivage dans chaque couche
186  REAL,DIMENSION(klon,klev)      :: zmasse    ! densité atmosphérique Kg/m2
187  REAL,DIMENSION(klon,klev)      :: ztra_th
188 
189!Controles:
190!---------
191  LOGICAL,SAVE :: couchelimite=.TRUE.
192  LOGICAL,SAVE :: convection=.TRUE.
193  LOGICAL,SAVE :: lessivage
194!$OMP THREADPRIVATE(couchelimite,convection,lessivage)
195
196  CHARACTER(len=8),DIMENSION(nbtr) :: solsym
197
198
199!######################################################################
200!                    -- INITIALIZATION --
201!######################################################################
202  IF (debutphy) THEN
203     WRITE(*,*) 'FIRST TIME IN PHYTRAC : pdtphys(sec) = ',pdtphys,'ecrit_tra (sec) = ',ecrit_tra
204     ALLOCATE( source(klon,nbtr), stat=ierr)
205     IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 1',1)
206     
207     ALLOCATE( aerosol(nbtr), stat=ierr)
208     IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 2',1)
209     
210
211     ! Initialize module for specific tracers
212     SELECT CASE(type_trac)
213     CASE('lmdz')
214        CALL traclmdz_init(pctsrf, ftsol, tr_seri, aerosol, lessivage)
215     CASE('inca')
216        source(:,:)=0.
217        CALL tracinca_init(aerosol,lessivage)
218     END SELECT
219!
220! Initialize diagnostic output
221! ----------------------------
222#ifdef CPP_IOIPSL
223!     INCLUDE "ini_histrac.h"
224#endif
225  END IF
226!############################################ END INITIALIZATION #######
227
228!===============================================================================
229!    -- Do specific treatment according to chemestry model or local LMDZ tracers
230!     
231!===============================================================================
232  SELECT CASE(type_trac)
233  CASE('lmdz')
234     !    -- Traitement des traceurs avec traclmdz
235     
236     CALL traclmdz(&
237          nstep,    pdtphys,      t_seri,           &
238          paprs,    pplay,        cdragh,  coefh,   &
239          yu1,      yv1,          ftsol,   pctsrf,  &
240          xlat,     couchelimite,                   &
241          tr_seri,  source,       solsym,  d_tr_cl)
242     
243  CASE('inca')
244     !    -- CHIMIE INCA  config_inca = aero or chem --
245
246     CALL tracinca(&
247          nstep,    julien,   gmtime,         lafin,     &
248          pdtphys,  t_seri,   paprs,          pplay,     &
249          pmfu,     ftsol,    pctsrf,         pphis,     &
250          pphi,     albsol,   sh,             rh,        &
251          cldfra,   rneb,     diafra,         cldliq,    &
252          itop_con, ibas_con, pmflxr,         pmflxs,    &
253          prfl,     psfl,     aerosol_couple, flxmass_w, &
254          tau_aero, piz_aero, cg_aero,        ccm,       &
255          rfname,                                        &
256          tr_seri,  source,   solsym)     
257  END SELECT
258
259!======================================================================
260!       -- Calcul de l'effet de la convection --
261!======================================================================
262  IF (convection) THEN
263     DO it=1, nbtr
264        IF ( conv_flg(it) == 0 ) CYCLE
265       
266        IF (iflag_con.LT.2) THEN
267           d_tr_cv(:,:,:)=0.
268        ELSE IF (iflag_con.EQ.2) THEN
269!..Tiedke
270           CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
271                pplay, paprs, tr_seri(:,:,it), d_tr_cv(:,:,it))
272        ELSE
273!..K.Emanuel
274           CALL cvltr(pdtphys, da, phi, mp, paprs,pplay, tr_seri(:,:,it),&
275                upwd,dnwd,d_tr_cv(:,:,it))
276        END IF
277
278        DO k = 1, klev
279           DO i = 1, klon       
280              tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
281           END DO
282        END DO
283
284        CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'convection it = '//solsym(it))
285             
286     END DO ! nbtr
287  END IF ! convection
288
289!======================================================================
290!    -- Calcul de l'effet des thermiques --
291!======================================================================
292
293  DO k=1,klev
294     DO i=1,klon
295        zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
296     END DO
297  END DO
298
299  DO it=1,nbtr
300     DO k=1,klev
301        DO i=1,klon
302           d_tr_th(i,k,it)=0.
303           tr_seri(i,k,it)=MAX(tr_seri(i,k,it),0.)
304           tr_seri(i,k,it)=MIN(tr_seri(i,k,it),1.e10)
305        END DO
306     END DO
307  END DO
308 
309  IF (iflag_thermals.GT.0) THEN   
310     nsplit=10
311     DO it=1, nbtr
312        DO isplit=1,nsplit
313
314           CALL dqthermcell(klon,klev,pdtphys/nsplit, &
315                fm_therm,entr_therm,zmasse, &
316                tr_seri(1:klon,1:klev,it),d_tr,ztra_th)
317
318           DO k=1,klev
319              DO i=1,klon
320                 d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit
321                 d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k)
322                 tr_seri(i,k,it)=MAX(tr_seri(i,k,it)+d_tr(i,k),0.)
323              END DO
324           END DO
325        END DO ! nsplit
326     END DO ! it
327  END IF ! Thermiques
328
329!======================================================================
330!     -- Calcul de l'effet de la couche limite --
331!======================================================================
332
333  IF (couchelimite) THEN
334
335     DO k = 1, klev
336        DO i = 1, klon
337           delp(i,k) = paprs(i,k)-paprs(i,k+1)
338        END DO
339     END DO
340
341     DO it=1, nbtr
342       
343        IF( pbl_flg(it) /= 0 ) THEN
344       
345           CALL cltrac(pdtphys, coefh,t_seri,       &
346                tr_seri(:,:,it), source(:,it),      &
347                paprs, pplay, delp,                 &
348                d_tr_cl(:,:,it))
349           
350           DO k = 1, klev
351              DO i = 1, klon
352                 tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
353              END DO
354           END DO
355        END IF
356
357     END DO
358     
359  END IF ! couche limite
360
361
362!======================================================================
363!   Calcul de l'effet de la precipitation
364!======================================================================
365
366  IF (lessivage) THEN
367     
368     d_tr_lessi_nucl(:,:,:) = 0.
369     d_tr_lessi_impa(:,:,:) = 0.
370     flestottr(:,:,:) = 0.
371!=========================
372! LESSIVAGE LARGE SCALE :
373!=========================
374
375! Tendance des aerosols nuclees et impactes
376! -----------------------------------------
377     DO it = 1, nbtr
378        IF (aerosol(it)) THEN
379           DO k = 1, klev
380              DO i = 1, klon
381                 d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +    &
382                      ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
383                 d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +    &
384                      ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
385
386!
387! Flux lessivage total
388! ------------------------------------------------------------
389                 flestottr(i,k,it) = flestottr(i,k,it) -   &
390                      ( d_tr_lessi_nucl(i,k,it)   +        &
391                      d_tr_lessi_impa(i,k,it) ) *          &
392                      ( paprs(i,k)-paprs(i,k+1) ) /        &
393                      (RG * pdtphys)
394!
395! Mise a jour des traceurs due a l'impaction,nucleation
396! ----------------------------------------------------------------------
397                 tr_seri(i,k,it)=tr_seri(i,k,it)*frac_impa(i,k)*frac_nucl(i,k)
398              END DO
399           END DO
400        END IF
401     END DO
402     
403  END IF ! lessivage
404
405!=============================================================
406!   Ecriture des sorties
407!=============================================================
408#ifdef CPP_IOIPSL
409!  INCLUDE "write_histrac.h"
410#endif
411
412END SUBROUTINE phytrac
Note: See TracBrowser for help on using the repository browser.