source: LMDZ5/trunk/libf/phylmd/traclmdz_mod.F90 @ 3373

Last change on this file since 3373 was 2320, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation: make an infotrac_phy module, which should be used from within the physics, and is initialized from infotrac (dynamics) via iniphysiq.
EM

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 23.4 KB
RevLine 
[1191]1!$Id $
2!
3MODULE traclmdz_mod
[1742]4
[1191]5!
6! In this module all tracers specific to LMDZ are treated. This module is used
7! only if running without any other chemestry model as INCA or REPROBUS. 
8!
[1403]9  IMPLICIT NONE
10
[1191]11  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: masktr   ! Masque reservoir de sol traceur
12!$OMP THREADPRIVATE(masktr)                        ! Masque de l'echange avec la surface (1 = reservoir) ou (possible >= 1 )
13  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: fshtr    ! Flux surfacique dans le reservoir de sol
14!$OMP THREADPRIVATE(fshtr)
15  REAL,DIMENSION(:),ALLOCATABLE,SAVE   :: hsoltr   ! Epaisseur equivalente du reservoir de sol
16!$OMP THREADPRIVATE(hsoltr)
17!
18!Radioelements:
19!--------------
20!
21  REAL,DIMENSION(:),ALLOCATABLE,SAVE   :: tautr    ! Constante de decroissance radioactive
22!$OMP THREADPRIVATE(tautr)
23  REAL,DIMENSION(:),ALLOCATABLE,SAVE   :: vdeptr   ! Vitesse de depot sec dans la couche Brownienne
24!$OMP THREADPRIVATE(vdeptr)
25  REAL,DIMENSION(:),ALLOCATABLE,SAVE   :: scavtr   ! Coefficient de lessivage
26!$OMP THREADPRIVATE(scavtr)
27  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: srcbe    ! Production du beryllium7 dans l atmosphere (U/s/kgA)
28!$OMP THREADPRIVATE(srcbe)
29
30  LOGICAL,DIMENSION(:),ALLOCATABLE,SAVE :: radio    ! radio(it)   = true  => decroisssance radioactive
31!$OMP THREADPRIVATE(radio) 
32
33  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: trs     ! Conc. radon ds le sol
34!$OMP THREADPRIVATE(trs)
35
[1409]36  INTEGER,SAVE :: id_aga      ! Identification number for tracer : Age of stratospheric air
37!$OMP THREADPRIVATE(id_aga)
38  INTEGER,SAVE :: lev_1p5km   ! Approximative vertical layer number at 1.5km above surface, used for calculation of the age of air. The result shouldn't be that sensible to the exactness of this value as long as it is in the lower troposphere.
39!$OMP THREADPRIVATE(lev_1p5km)
40
41  INTEGER,SAVE :: id_rn, id_pb ! Identification number for tracer : radon (Rn222), lead (Pb210)
42!$OMP THREADPRIVATE(id_rn, id_pb)
43
[1191]44  INTEGER,SAVE :: id_be       ! Activation et position du traceur Be7 [ id_be=0 -> desactive ]
45!$OMP THREADPRIVATE(id_be)
46
[1403]47  INTEGER,SAVE :: id_pcsat, id_pcocsat, id_pcq ! traceurs pseudo-vapeur CL qsat, qsat_oc, q
48!$OMP THREADPRIVATE(id_pcsat, id_pcocsat, id_pcq)
[1421]49  INTEGER,SAVE :: id_pcs0, id_pcos0, id_pcq0   ! traceurs pseudo-vapeur CL qsat, qsat_oc, q
50!                                              ! qui ne sont pas transportes par la convection
[1403]51!$OMP THREADPRIVATE(id_pcs0, id_pcos0, id_pcq0)
52
53  INTEGER, SAVE:: id_o3
[1421]54!$OMP THREADPRIVATE(id_o3)
55! index of ozone tracer with Cariolle parameterization
56! 0 means no ozone tracer
[1403]57
[1409]58  LOGICAL,SAVE :: rnpb=.FALSE. ! Presence du couple Rn222, Pb210
[1191]59!$OMP THREADPRIVATE(rnpb)
60
61
62CONTAINS
63
64  SUBROUTINE traclmdz_from_restart(trs_in)
65    ! This subroutine initialize the module saved variable trs with values from restart file (startphy.nc).
66    ! This subroutine is called from phyetat0 after the field trs_in has been read.
67   
68    USE dimphy
[2320]69    USE infotrac_phy
[1191]70   
71    ! Input argument
72    REAL,DIMENSION(klon,nbtr), INTENT(IN) :: trs_in
73   
74    ! Local variables
75    INTEGER :: ierr
76   
77    ! Allocate restart variables trs
78    ALLOCATE( trs(klon,nbtr), stat=ierr)
[2311]79    IF (ierr /= 0) CALL abort_physic('traclmdz_from_restart', 'pb in allocation 1',1)
[1191]80   
81    ! Initialize trs with values read from restart file
82    trs(:,:) = trs_in(:,:)
83   
84  END SUBROUTINE traclmdz_from_restart
85
86
[1579]87  SUBROUTINE traclmdz_init(pctsrf, xlat, xlon, ftsol, tr_seri, t_seri, pplay, sh, pdtphys, aerosol, lessivage)
[1191]88    ! This subroutine allocates and initialize module variables and control variables.
[1421]89    ! Initialization of the tracers should be done here only for those not found in the restart file.
[1191]90    USE dimphy
[2320]91    USE infotrac_phy
[1403]92    USE regr_pr_comb_coefoz_m, ONLY: alloc_coefoz
93    USE press_coefoz_m, ONLY: press_coefoz
[1227]94    USE carbon_cycle_mod, ONLY : carbon_cycle_init, carbon_cycle_tr, carbon_cycle_cpl
[1421]95    USE mod_grid_phy_lmdz
96    USE mod_phys_lmdz_para
[1785]97    USE indice_sol_mod
[2311]98    USE print_control_mod, ONLY: lunout
[1191]99
100! Input variables
[1227]101    REAL,DIMENSION(klon,nbsrf),INTENT(IN)     :: pctsrf ! Pourcentage de sol f(nature du sol)
[1579]102    REAL,DIMENSION(klon),INTENT(IN)           :: xlat   ! latitudes en degres pour chaque point
103    REAL,DIMENSION(klon),INTENT(IN)           :: xlon   ! longitudes en degres pour chaque point
[1227]104    REAL,DIMENSION(klon,nbsrf),INTENT(IN)     :: ftsol  ! Temperature du sol (surf)(Kelvin)
[1403]105    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri! Concentration Traceur [U/KgA] 
106    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
107    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
108    REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique
[1454]109    REAL,INTENT(IN)                        :: pdtphys ! Pas d'integration pour la physique (seconde) 
[1227]110
[1191]111! Output variables
112    LOGICAL,DIMENSION(nbtr), INTENT(OUT) :: aerosol
113    LOGICAL,INTENT(OUT)                  :: lessivage
114       
115! Local variables   
[1403]116    INTEGER :: ierr, it, iiq, i, k
[1421]117    REAL, DIMENSION(klon_glo,klev) :: varglo ! variable temporaire sur la grille global   
118    REAL, DIMENSION(klev)          :: mintmp, maxtmp
119    LOGICAL                        :: zero
[1742]120! RomP >>> profil initial Be7
121      integer ilesfil
122      parameter (ilesfil=1)
123      integer  irr,kradio
124      real     beryllium(klon,klev)
125! profil initial Pb210
126      integer ilesfil2
127      parameter (ilesfil2=1)
128      integer  irr2,kradio2
129      real     plomb(klon,klev)
130!! RomP <<<
[1191]131! --------------------------------------------
132! Allocation
133! --------------------------------------------
134    ALLOCATE( scavtr(nbtr), stat=ierr)
[2311]135    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 9',1)
[1191]136    scavtr(:)=1.
137   
138    ALLOCATE( radio(nbtr), stat=ierr)
[2311]139    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 11',1)
[1191]140    radio(:) = .false.    ! Par defaut pas decroissance radioactive
141   
142    ALLOCATE( masktr(klon,nbtr), stat=ierr)
[2311]143    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 2',1)
[1191]144   
145    ALLOCATE( fshtr(klon,nbtr), stat=ierr)
[2311]146    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 3',1)
[1191]147   
148    ALLOCATE( hsoltr(nbtr), stat=ierr)
[2311]149    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 4',1)
[1191]150   
151    ALLOCATE( tautr(nbtr), stat=ierr)
[2311]152    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 5',1)
[1191]153    tautr(:)  = 0.
154   
155    ALLOCATE( vdeptr(nbtr), stat=ierr)
[2311]156    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 6',1)
[1191]157    vdeptr(:) = 0.
158
159
160    lessivage  = .TRUE.
[1742]161!!jyg(20130206) : le choix d activation du lessivage est fait dans phytrac avec iflag_lscav
162!!    call getin('lessivage',lessivage)
163!!    if(lessivage) then
164!!     print*,'lessivage lsc ON'
165!!    else
166!!     print*,'lessivage lsc OFF'
167!!    endif
[1191]168    aerosol(:) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
169   
170!
[1403]171! Recherche des traceurs connus : Be7, O3, CO2,...
[1191]172! --------------------------------------------
[1409]173    id_rn=0; id_pb=0; id_aga=0; id_be=0; id_o3=0
[1421]174    id_pcsat=0; id_pcocsat=0; id_pcq=0; id_pcs0=0; id_pcos0=0; id_pcq0=0
[1191]175    DO it=1,nbtr
[2265]176!!       iiq=niadv(it+2)                                                            ! jyg
177       iiq=niadv(it+nqo)                                                            ! jyg
178       IF ( tname(iiq) == "RN" ) THEN                                               
[1409]179          id_rn=it ! radon
180       ELSE IF ( tname(iiq) == "PB") THEN
181          id_pb=it ! plomb
[1742]182! RomP >>> profil initial de PB210
183     open (ilesfil2,file='prof.pb210',status='old',iostat=irr2)
184     IF (irr2 == 0) THEN
185      read(ilesfil2,*) kradio2
186      print*,'number of levels for pb210 profile ',kradio2
187      do k=kradio2,1,-1
188       read (ilesfil2,*) plomb(:,k)
189      enddo
190      close(ilesfil2)
191      do k=1,klev
192       do i=1,klon
193         tr_seri(i,k,id_pb)=plomb(i,k)
194!!        print*, 'tr_seri',i,k,tr_seri(i,k,id_pb)
195        enddo
196      enddo
197     ELSE
198       print *, 'Prof.pb210 does not exist: use restart values'
199     ENDIF
200! RomP <<<
[1409]201       ELSE IF ( tname(iiq) == "Aga" .OR. tname(iiq) == "AGA" ) THEN
202          ! Age of stratospheric air
203          id_aga=it
204          radio(id_aga) = .FALSE.
205          aerosol(id_aga) = .FALSE.
206          pbl_flg(id_aga) = 0
207         
208          ! Find the first model layer above 1.5km from the surface
209          IF (klev>=30) THEN
210             lev_1p5km=6   ! NB! This value is for klev=39
211          ELSE IF (klev>=10) THEN
212             lev_1p5km=5   ! NB! This value is for klev=19
213          ELSE
214             lev_1p5km=klev/2
215          END IF
216       ELSE IF ( tname(iiq) == "BE" .OR. tname(iiq) == "Be" .OR.  &
[1191]217            tname(iiq) == "BE7" .OR. tname(iiq) == "Be7" ) THEN 
[1212]218          ! Recherche du Beryllium 7
[1191]219          id_be=it
220          ALLOCATE( srcbe(klon,klev) )
221          radio(id_be) = .TRUE.
222          aerosol(id_be) = .TRUE. ! le Be est un aerosol
[1742]223!jyg le 13/03/2013 ; ajout de pplay en argument de init_be
224!!!          CALL init_be(pctsrf,masktr(:,id_be),tautr(id_be),vdeptr(id_be),scavtr(id_be),srcbe)
225          CALL init_be(pctsrf,pplay,masktr(:,id_be),tautr(id_be),vdeptr(id_be),scavtr(id_be),srcbe)
[1421]226          WRITE(lunout,*) 'Initialisation srcBe: OK'
[1742]227! RomP >>> profil initial de Be7
228      open (ilesfil,file='prof.be7',status='old',iostat=irr)
229      IF (irr == 0) THEN
230       read(ilesfil,*) kradio
231       print*,'number of levels for Be7 profile ',kradio
232       do k=kradio,1,-1
233        read (ilesfil,*) beryllium(:,k)
234       enddo
235       close(ilesfil)
236       do k=1,klev
237         do i=1,klon
238          tr_seri(i,k,id_be)=beryllium(i,k)
239!!        print*, 'tr_seri',i,k,tr_seri(i,k,id_be)
240         enddo
241       enddo
242     ELSE
243       print *, 'Prof.Be7 does not exist: use restart values'
244     ENDIF
245! RomP <<<
[1403]246       ELSE IF (tname(iiq)=="O3" .OR. tname(iiq)=="o3") THEN
247          ! Recherche de l'ozone : parametrization de la chimie par Cariolle
248          id_o3=it
249          CALL alloc_coefoz   ! allocate ozone coefficients
250          CALL press_coefoz   ! read input pressure levels
[1421]251       ELSE IF ( tname(iiq) == "pcsat" .OR. tname(iiq) == "Pcsat" ) THEN
252          id_pcsat=it
253       ELSE IF ( tname(iiq) == "pcocsat" .OR. tname(iiq) == "Pcocsat" ) THEN
254          id_pcocsat=it
255       ELSE IF ( tname(iiq) == "pcq" .OR. tname(iiq) == "Pcq" ) THEN
256          id_pcq=it
257       ELSE IF ( tname(iiq) == "pcs0" .OR. tname(iiq) == "Pcs0" ) THEN
258          id_pcs0=it
259          conv_flg(it)=0 ! No transport by convection for this tracer
260       ELSE IF ( tname(iiq) == "pcos0" .OR. tname(iiq) == "Pcos0" ) THEN
261          id_pcos0=it
262          conv_flg(it)=0 ! No transport by convection for this tracer
263       ELSE IF ( tname(iiq) == "pcq0" .OR. tname(iiq) == "Pcq0" ) THEN
264          id_pcq0=it
265          conv_flg(it)=0 ! No transport by convection for this tracer
266       ELSE
267          WRITE(lunout,*) 'This is an unknown tracer in LMDZ : ', trim(tname(iiq))
[1403]268       END IF
269    END DO
270
[1191]271!
272! Valeurs specifiques pour les traceurs Rn222 et Pb210
273! ----------------------------------------------
[1409]274    IF ( id_rn/=0 .AND. id_pb/=0 ) THEN
275       rnpb = .TRUE.
276       radio(id_rn)= .TRUE.
277       radio(id_pb)= .TRUE.
278       pbl_flg(id_rn) = 0 ! au lieu de clsol=true ! CL au sol calcule
279       pbl_flg(id_pb) = 0 ! au lieu de clsol=true
280       aerosol(id_rn) = .FALSE.
281       aerosol(id_pb) = .TRUE. ! le Pb est un aerosol
[1191]282       
283       CALL initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr)
284    END IF
285
[1227]286!
287! Initialisation de module carbon_cycle_mod
288! ----------------------------------------------
289    IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
[1454]290       CALL carbon_cycle_init(tr_seri, pdtphys, aerosol, radio)
[1227]291    END IF
292
[1421]293! Check if all tracers have restart values
294! ----------------------------------------------
295    DO it=1,nbtr
[2265]296!!       iiq=niadv(it+2)                                                            ! jyg
297       iiq=niadv(it+nqo)                                                            ! jyg
[1421]298       ! Test if tracer is zero everywhere.
299       ! Done by master process MPI and master thread OpenMP
300       CALL gather(tr_seri(:,:,it),varglo)
301       IF (is_mpi_root .AND. is_omp_root) THEN
302          mintmp=MINVAL(varglo,dim=1)
303          maxtmp=MAXVAL(varglo,dim=1)
304          IF (MINVAL(mintmp,dim=1)==0. .AND. MAXVAL(maxtmp,dim=1)==0.) THEN
305             ! Tracer is zero everywhere
306             zero=.TRUE.
307          ELSE
308             zero=.FALSE.
309          END IF
[1403]310       END IF
311
[1421]312       ! Distribute variable at all processes
313       CALL bcast(zero)
[1403]314
[1421]315       ! Initalize tracer that was not found in restart file.
316       IF (zero) THEN
317          ! The tracer was not found in restart file or it was equal zero everywhere.
318          WRITE(lunout,*) "The tracer ",trim(tname(iiq))," will be initialized"
319          IF (it==id_pcsat .OR. it==id_pcq .OR. &
320               it==id_pcs0 .OR. it==id_pcq0) THEN
321             tr_seri(:,:,it) = 100.
322          ELSE IF (it==id_pcocsat .OR. it==id_pcos0) THEN
323             DO i = 1, klon
324                IF ( pctsrf (i, is_oce) == 0. ) THEN
325                   tr_seri(i,:,it) = 0.
326                ELSE
327                   tr_seri(i,:,it) = 100.
328                END IF
329             END DO
330          ELSE
331             ! No specific initialization exist for this tracer
332             tr_seri(:,:,it) = 0.
333          END IF
[1403]334       END IF
[1421]335    END DO
[1403]336
[1191]337  END SUBROUTINE traclmdz_init
338
[1403]339  SUBROUTINE traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
340       cdragh,  coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon, couchelimite, sh, &
[1813]341       rh, pphi, ustar, wstar, ale_bl, ale_wake,  zu10m, zv10m, &
[2180]342       tr_seri, source, d_tr_cl,d_tr_dec, zmasse)               !RomP
[1191]343   
344    USE dimphy
[2320]345    USE infotrac_phy
[1403]346    USE regr_pr_comb_coefoz_m, ONLY: regr_pr_comb_coefoz
347    USE o3_chem_m, ONLY: o3_chem
[1227]348    USE carbon_cycle_mod, ONLY : carbon_cycle, carbon_cycle_tr, carbon_cycle_cpl
[1785]349    USE indice_sol_mod
350
[1191]351    INCLUDE "YOMCST.h"
352
353!==========================================================================
354!                   -- DESCRIPTION DES ARGUMENTS --
355!==========================================================================
356
357! Input arguments
358!
359!Configuration grille,temps:
[1212]360    INTEGER,INTENT(IN) :: nstep      ! nombre d'appels de la physiq
[1403]361    INTEGER,INTENT(IN) :: julien     ! Jour julien
362    REAL,INTENT(IN)    :: gmtime
[1191]363    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde) 
364    REAL,DIMENSION(klon),INTENT(IN) :: xlat    ! latitudes pour chaque point
[1403]365    REAL, INTENT(IN):: xlon(:) ! dim(klon) longitude
[1191]366
367!
368!Physique:
369!--------
370    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
371    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
372    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
[1403]373    REAL,intent(in):: zmasse (:, :)   ! dim(klon,klev) density of air, in kg/m2
[1191]374
375
376!Couche limite:
377!--------------
378!
[1742]379    REAL,DIMENSION(klon),INTENT(IN)      :: cdragh     ! coeff drag pour T et Q
380    REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh      ! coeff melange CL (m**2/s)
381    REAL,DIMENSION(klon),INTENT(IN)      :: yu1        ! vents au premier niveau
382    REAL,DIMENSION(klon),INTENT(IN)      :: yv1        ! vents au premier niveau
[1191]383    LOGICAL,INTENT(IN)                   :: couchelimite
[1742]384    REAL,DIMENSION(klon,klev),INTENT(IN) :: sh         ! humidite specifique
[1670]385    REAL,DIMENSION(klon,klev),INTENT(IN) :: rh      ! Humidite relative
386    REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi    ! geopotentie
387    REAL,DIMENSION(klon),INTENT(IN)      :: ustar   ! ustar (m/s)
[1813]388    REAL,DIMENSION(klon),INTENT(IN)      :: wstar,ale_bl,ale_wake   ! wstar (m/s) and Avail. Lifti. Energ.
[1670]389    REAL,DIMENSION(klon),INTENT(IN)      :: zu10m   ! vent zonal 10m (m/s)
390    REAL,DIMENSION(klon),INTENT(IN)      :: zv10m   ! vent zonal 10m (m/s)
[1191]391
392! Arguments necessaires pour les sources et puits de traceur:
393    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol  ! Temperature du sol (surf)(Kelvin)
394    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
395
396! InOutput argument
397    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA] 
398
399! Output argument
400    REAL,DIMENSION(klon,nbtr), INTENT(OUT)        :: source  ! a voir lorsque le flux de surface est prescrit
401    REAL,DIMENSION(klon,klev,nbtr), INTENT(OUT)   :: d_tr_cl ! Td couche limite/traceur
402
403!=======================================================================================
404!                        -- VARIABLES LOCALES TRACEURS --
405!=======================================================================================
406
407    INTEGER :: i, k, it
[1421]408    INTEGER :: lmt_pas ! number of time steps of "physics" per day
[1191]409
410    REAL,DIMENSION(klon)           :: d_trs    ! Td dans le reservoir
[1421]411    REAL,DIMENSION(klon,klev)      :: qsat     ! pression de la vapeur a saturation
[1191]412    REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec ! Td radioactive
[1421]413    REAL                           :: zrho     ! Masse Volumique de l'air KgA/m3
414    REAL                           :: amn, amx
[1191]415!
416!=================================================================
417!  Ajout de la production en  Be7 (Beryllium) srcbe U/s/kgA
418!=================================================================
419!
420    IF ( id_be /= 0 ) THEN
421       DO k = 1, klev
422          DO i = 1, klon
423             tr_seri(i,k,id_be) = tr_seri(i,k,id_be)+srcbe(i,k)*pdtphys
424          END DO
425       END DO
426       WRITE(*,*) 'Ajout srcBe dans tr_seri: OK'
427    END IF
428 
[1421]429
430!=================================================================
431! Update pseudo-vapor tracers
432!=================================================================
433
434    CALL q_sat(klon*klev,t_seri,pplay,qsat)
435
[1403]436    IF ( id_pcsat /= 0 ) THEN
437     DO k = 1, klev
438      DO i = 1, klon
[1421]439         IF ( pplay(i,k).GE.85000.) THEN
440            tr_seri(i,k,id_pcsat) = qsat(i,k)
441         ELSE
442            tr_seri(i,k,id_pcsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcsat))           
443         END IF
[1403]444      END DO
445     END DO
446    END IF
[1191]447
[1403]448    IF ( id_pcocsat /= 0 ) THEN
449     DO k = 1, klev
450      DO i = 1, klon
[1421]451         IF ( pplay(i,k).GE.85000.) THEN
452            IF ( pctsrf (i, is_oce) > 0. ) THEN
453               tr_seri(i,k,id_pcocsat) = qsat(i,k)
454            ELSE
455               tr_seri(i,k,id_pcocsat) = 0.
456          END IF
457       ELSE
458          tr_seri(i,k,id_pcocsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcocsat))
[1403]459       END IF
460      END DO
461     END DO
462    END IF
463
464    IF ( id_pcq /= 0 ) THEN
465     DO k = 1, klev
466      DO i = 1, klon
[1421]467         IF ( pplay(i,k).GE.85000.) THEN
468            tr_seri(i,k,id_pcq) = sh(i,k)
469         ELSE
470            tr_seri(i,k,id_pcq) = MIN (qsat(i,k), tr_seri(i,k,id_pcq))   
471         END IF
[1403]472      END DO
473     END DO
474    END IF
475
[1421]476
[1403]477    IF ( id_pcs0 /= 0 ) THEN
478     DO k = 1, klev
479      DO i = 1, klon
[1421]480         IF ( pplay(i,k).GE.85000.) THEN
481            tr_seri(i,k,id_pcs0) = qsat(i,k)
482         ELSE
483            tr_seri(i,k,id_pcs0) = MIN (qsat(i,k), tr_seri(i,k,id_pcs0))   
484         END IF
[1403]485      END DO
486     END DO
487    END IF
488
[1421]489
[1403]490    IF ( id_pcos0 /= 0 ) THEN
491     DO k = 1, klev
492      DO i = 1, klon
[1421]493         IF ( pplay(i,k).GE.85000.) THEN
494            IF ( pctsrf (i, is_oce) > 0. ) THEN
495               tr_seri(i,k,id_pcos0) = qsat(i,k)
496            ELSE
497               tr_seri(i,k,id_pcos0) = 0.
498            END IF
499         ELSE
500            tr_seri(i,k,id_pcos0) = MIN (qsat(i,k), tr_seri(i,k,id_pcos0))
501         END IF
[1403]502      END DO
503     END DO
504    END IF
505
[1421]506
[1403]507    IF ( id_pcq0 /= 0 ) THEN
508     DO k = 1, klev
509      DO i = 1, klon
[1421]510         IF ( pplay(i,k).GE.85000.) THEN
511            tr_seri(i,k,id_pcq0) = sh(i,k)
512         ELSE
513            tr_seri(i,k,id_pcq0) = MIN (qsat(i,k), tr_seri(i,k,id_pcq0))
514         END IF
[1403]515      END DO
516     END DO
517    END IF
518
[1409]519!=================================================================
520! Update tracer : Age of stratospheric air
521!=================================================================
522    IF (id_aga/=0) THEN
523       
524       ! Bottom layers
525       DO k = 1, lev_1p5km
526          tr_seri(:,k,id_aga) = 0.0
527       END DO
528       
529       ! Layers above 1.5km
530       DO k = lev_1p5km+1,klev-1
531          tr_seri(:,k,id_aga) = tr_seri(:,k,id_aga) + pdtphys
532       END DO
533       
534       ! Top layer
535       tr_seri(:,klev,id_aga) = tr_seri(:,klev-1,id_aga)
536       
537    END IF
538
[1191]539!======================================================================
540!     -- Calcul de l'effet de la couche limite --
541!======================================================================
542
543    IF (couchelimite) THEN             
[1212]544       source(:,:) = 0.0
545
546       IF (id_be /=0) THEN
547          DO i=1, klon
548             zrho = pplay(i,1)/t_seri(i,1)/RD
549             source(i,id_be) = - vdeptr(id_be)*tr_seri(i,1,id_be)*zrho
550          END DO
551       END IF
552
[1191]553    END IF
554   
555    DO it=1, nbtr
[1409]556       IF (couchelimite .AND. pbl_flg(it) == 0 .AND. (it==id_rn .OR. it==id_pb)) THEN
557          ! couche limite avec quantite dans le sol calculee
[1191]558          CALL cltracrn(it, pdtphys, yu1, yv1,     &
559               cdragh, coefh,t_seri,ftsol,pctsrf,  &
560               tr_seri(:,:,it),trs(:,it),          &
[1403]561               paprs, pplay, zmasse * rg, &
[1191]562               masktr(:,it),fshtr(:,it),hsoltr(it),&
563               tautr(it),vdeptr(it),               &
564               xlat,d_tr_cl(:,:,it),d_trs)
565         
566          DO k = 1, klev
567             DO i = 1, klon
568                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
569             END DO
570          END DO
571       
572          ! Traceur dans le reservoir sol
573          DO i = 1, klon
574             trs(i,it) = trs(i,it) + d_trs(i)
575          END DO
576       END IF
577    END DO
[1403]578         
579
[1191]580!======================================================================
581!   Calcul de l'effet du puits radioactif
582!======================================================================
583    CALL radio_decay (radio,rnpb,pdtphys,tautr,tr_seri,d_tr_dec)
[1421]584
[1191]585    DO it=1,nbtr
[1421]586       WRITE(solsym(it),'(i2)') it
587    END DO
588
589    DO it=1,nbtr
[1191]590       IF(radio(it)) then     
591          DO k = 1, klev
592             DO i = 1, klon
593                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
594             END DO
595          END DO
596          CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'puits rn it='//solsym(it))
597       END IF
598    END DO
599
[1227]600!======================================================================
[1403]601!   Parameterization of ozone chemistry
602!======================================================================
603
604    IF (id_o3 /= 0) then
605       lmt_pas = NINT(86400./pdtphys)
606       IF (MOD(nstep - 1, lmt_pas) == 0) THEN
607          ! Once per day, update the coefficients for ozone chemistry:
608          CALL regr_pr_comb_coefoz(julien, xlat, paprs, pplay)
609       END IF
610       CALL o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, xlat, &
611            xlon, tr_seri(:, :, id_o3))
612    END IF
613
614!======================================================================
[1227]615!   Calcul de cycle de carbon
616!======================================================================
617    IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
[1454]618       CALL carbon_cycle(nstep, pdtphys, pctsrf, tr_seri, source)
[1227]619    END IF
620
[1191]621  END SUBROUTINE traclmdz
622
623
624  SUBROUTINE traclmdz_to_restart(trs_out)
625    ! This subroutine is called from phyredem.F where the module
626    ! variable trs is written to restart file (restartphy.nc)
627    USE dimphy
[2320]628    USE infotrac_phy
[1191]629   
630    REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: trs_out
[1203]631    INTEGER :: ierr
[1212]632
633    IF ( ALLOCATED(trs) ) THEN
634       trs_out(:,:) = trs(:,:)
[1203]635    ELSE
[1212]636       ! No previous allocate of trs. This is the case for create_etat0_limit.
637       trs_out(:,:) = 0.0
[1203]638    END IF
[1212]639   
[1191]640  END SUBROUTINE traclmdz_to_restart
641 
642
643END MODULE traclmdz_mod
Note: See TracBrowser for help on using the repository browser.