source: LMDZ5/trunk/libf/phylmd/phytrac.F90 @ 1605

Last change on this file since 1605 was 1579, checked in by jghattas, 13 years ago

Added latitudes and longitudes in traclmdz_init to simplify initialization of new tracers.

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