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

Last change on this file since 1667 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
RevLine 
[1191]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)         
[524]20!
[1191]21!======================================================================
22! Auteur(s) FH
23! Objet: Moniteur general des tendances traceurs
24!======================================================================
[524]25
[1191]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
[1565]35  USE tracreprobus_mod
[1403]36  USE control_mod
[524]37
38
[1191]39  IMPLICIT NONE
[959]40
[1191]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"
[1575]48  INCLUDE "iniprint.h"
[1191]49!==========================================================================
50!                   -- ARGUMENT DESCRIPTION --
51!==========================================================================
[541]52
[1191]53! Input arguments
54!----------------
55!Configuration grille,temps:
56  INTEGER,INTENT(IN) :: nstep      ! Appel physique
57  INTEGER,INTENT(IN) :: julien     ! Jour julien
[1565]58  REAL,INTENT(IN)    :: gmtime     ! Heure courante
[1191]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
[1454]69  REAL,DIMENSION(klon,klev),INTENT(IN)   :: u       ! variable not used
70  REAL,DIMENSION(klon,klev),INTENT(IN)   :: v       ! variable not used
[1191]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
[619]94
[1191]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]
[524]98
[1191]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!
[1454]121  REAL,DIMENSION(klon),INTENT(IN)      :: cdragh ! coeff drag pour T et Q
[1191]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
[524]133
[1191]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)
[524]137
138
[1191]139! Output argument
140!----------------
141  REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] 
[524]142
[1191]143!=======================================================================================
144!                        -- LOCAL VARIABLES --
145!=======================================================================================
[766]146
[1191]147  INTEGER :: i, k, it
148  INTEGER :: nsplit
[524]149
[1191]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)
[524]155
[1191]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.
[524]167
[1191]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)
[766]197
[1191]198  CHARACTER(len=8),DIMENSION(nbtr) :: solsym
[524]199
200
[1191]201!######################################################################
202!                    -- INITIALIZATION --
203!######################################################################
204  IF (debutphy) THEN
[1575]205     IF (prt_level >9) WRITE(lunout,*) 'FIRST TIME IN PHYTRAC : pdtphys(sec) = ',pdtphys,'ecrit_tra (sec) = ',ecrit_tra
[1191]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     
[524]212
[1191]213     ! Initialize module for specific tracers
214     SELECT CASE(type_trac)
215     CASE('lmdz')
[1579]216        CALL traclmdz_init(pctsrf, xlat, xlon, ftsol, tr_seri, t_seri, pplay, sh, pdtphys, aerosol, lessivage)
[1191]217     CASE('inca')
218        source(:,:)=0.
219        CALL tracinca_init(aerosol,lessivage)
[1565]220     CASE('repr')
221        source(:,:)=0.
[1191]222     END SELECT
223!
224! Initialize diagnostic output
225! ----------------------------
[1192]226#ifdef CPP_IOIPSL
[1575]227     INCLUDE "ini_histrac.h"
[1192]228#endif
[1191]229  END IF
230!############################################ END INITIALIZATION #######
[524]231
[1403]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
[1191]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
[1403]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)
[1191]248  CASE('inca')
249     !    -- CHIMIE INCA  config_inca = aero or chem --
[776]250
[1191]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)     
[1565]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     
[1191]271  END SELECT
[524]272
[1191]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
[959]291
[1191]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
[959]297
[1191]298        CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'convection it = '//solsym(it))
299             
300     END DO ! nbtr
301  END IF ! convection
[959]302
[524]303!======================================================================
[1191]304!    -- Calcul de l'effet des thermiques --
[524]305!======================================================================
306
[1191]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)
[524]313        END DO
[1191]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
[524]321
[1191]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)
[524]325
[1191]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
[682]336
[1191]337!======================================================================
338!     -- Calcul de l'effet de la couche limite --
339!======================================================================
[682]340
[1191]341  IF (couchelimite) THEN
[524]342
[1191]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
[524]348
[1191]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
[959]364
[1191]365     END DO
366     
367  END IF ! couche limite
[959]368
[619]369
[1191]370!======================================================================
371!   Calcul de l'effet de la precipitation
372!======================================================================
[959]373
[1191]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!=========================
[524]382
[1191]383! Tendance des aerosols nuclees et impactes
384! -----------------------------------------
385     DO it = 1, nbtr
386        IF (aerosol(it)) THEN
387           DO k = 1, klev
[524]388              DO i = 1, klon
[1191]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)
[524]393
[1191]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
[524]412
[1191]413!=============================================================
414!   Ecriture des sorties
415!=============================================================
[1192]416#ifdef CPP_IOIPSL
[1575]417  INCLUDE "write_histrac.h"
[1192]418#endif
[524]419
[1191]420END SUBROUTINE phytrac
Note: See TracBrowser for help on using the repository browser.