Sources/phylmd/phytrac.f

module phytrac_m

  IMPLICIT none

  private
  public phytrac

contains

  SUBROUTINE phytrac(itap, lmt_pas, julien, gmtime, firstcal, lafin, pdtphys, &
       t_seri, paprs, pplay, pmfu, pmfd, pde_u, pen_d, coefh, fm_therm, &
       entr_therm, yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, da, phi, &
       mp, upwd, dnwd, tr_seri, zmasse, ncid_startphy, nid_ins)

    ! From phylmd/phytrac.F, version 1.15 2006/02/21 08:08:30 (SVN revision 679)

    ! Authors: Fr\'ed\'eric Hourdin, Abderrahmane Idelkadi, Marie-Alice
    ! Foujols, Olivia

    ! Objet : moniteur g\'en\'eral des tendances des traceurs

    ! L'appel de "phytrac" se fait avec "nqmx - 2" donc nous avons
    ! bien les vrais traceurs, sans la vapeur d'eau ni l'eau liquide.

    ! Modifications pour les traceurs :
    ! - uniformisation des parametrisations dans phytrac
    ! - stockage des moyennes des champs n\'ecessaires en mode traceur off-line 

    use abort_gcm_m, only: abort_gcm
    use clesphys, only: ecrit_tra
    use clesphys2, only: iflag_con
    use cltrac_m, only: cltrac
    use cltracrn_m, only: cltracrn
    use ctherm, only: iflag_thermals
    use cvltr_m, only: cvltr
    use dimens_m, only: llm, nqmx
    use dimphy, only: klon
    use indicesol, only: nbsrf
    use initrrnpb_m, only: initrrnpb
    use minmaxqfi_m, only: minmaxqfi
    use netcdf95, only: nf95_inq_varid, nf95_get_var, nf95_put_var
    use nflxtr_m, only: nflxtr
    use nr_util, only: assert
    use o3_chem_m, only: o3_chem
    use phyetat0_m, only: rlat
    use phyredem0_m, only: ncid_restartphy
    use press_coefoz_m, only: press_coefoz
    use radiornpb_m, only: radiornpb
    use regr_pr_comb_coefoz_m, only: regr_pr_comb_coefoz
    use SUPHEC_M, only: rg

    integer, intent(in):: itap ! number of calls to "physiq"
    integer, intent(in):: lmt_pas ! number of time steps of "physics" per day
    integer, intent(in):: julien !jour julien, 1 <= julien <= 360
    real, intent(in):: gmtime ! heure de la journ\'ee en fraction de jour
    logical, intent(in):: firstcal ! first call to "calfis"
    logical, intent(in):: lafin ! fin de la physique
    real, intent(in):: pdtphys ! pas d'integration pour la physique (s)
    real, intent(in):: t_seri(klon, llm) ! temperature, in K

    real, intent(in):: paprs(klon, llm+1)
    ! (pression pour chaque inter-couche, en Pa)

    real, intent(in):: pplay(klon, llm)
    ! (pression pour le mileu de chaque couche, en Pa)

    ! convection:

    REAL, intent(in):: pmfu(klon, llm) ! flux de masse dans le panache montant

    REAL, intent(in):: pmfd(klon, llm)
    ! flux de masse dans le panache descendant

    REAL pde_u(klon, llm) ! flux detraine dans le panache montant
    REAL pen_d(klon, llm) ! flux entraine dans le panache descendant
    REAL coefh(klon, llm) ! coeff melange couche limite

    ! thermiques:
    real fm_therm(klon, llm+1), entr_therm(klon, llm)

    ! Couche limite:
    REAL yu1(klon) ! vents au premier niveau
    REAL yv1(klon) ! vents au premier niveau

    ! Arguments n\'ecessaires pour les sources et puits de traceur :
    real, intent(in):: ftsol(klon, nbsrf) ! Temperature du sol (surf)(Kelvin)
    real pctsrf(klon, nbsrf) ! Pourcentage de sol f(nature du sol)

    ! Lessivage pour le on-line
    REAL frac_impa(klon, llm) ! fraction d'aerosols impactes
    REAL frac_nucl(klon, llm) ! fraction d'aerosols nuclees

    ! Kerry Emanuel
    real, intent(in):: da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
    REAL, intent(in):: upwd(klon, llm) ! saturated updraft mass flux
    REAL, intent(in):: dnwd(klon, llm) ! saturated downdraft mass flux

    real, intent(inout):: tr_seri(:, :, :) ! (klon, llm, nqmx - 2)
    ! (mass fractions of tracers, excluding water, at mid-layers)

    real, intent(in):: zmasse(:, :) ! (klon, llm)
    ! (column-density of mass of air in a cell, in kg m-2)

    integer, intent(in):: ncid_startphy, nid_ins

    ! Local:

    integer nsplit

    ! TRACEURS

    ! Sources et puits des traceurs:

    ! Pour l'instant seuls les cas du rn et du pb ont ete envisages.

    REAL source(klon) ! a voir lorsque le flux est prescrit 
    ! 
    ! Pour la source de radon et son reservoir de sol

    REAL, save:: trs(klon, nqmx - 2) ! Concentration de traceur dans le sol

    REAL masktr(klon, nqmx - 2) ! Masque reservoir de sol traceur
    ! Masque de l'echange avec la surface
    ! (1 = reservoir) ou (possible => 1)
    SAVE masktr
    REAL fshtr(klon, nqmx - 2) ! Flux surfacique dans le reservoir de sol
    SAVE fshtr
    REAL hsoltr(nqmx - 2) ! Epaisseur equivalente du reservoir de sol
    SAVE hsoltr
    REAL tautr(nqmx - 2) ! Constante de decroissance radioactive
    SAVE tautr
    REAL vdeptr(nqmx - 2) ! Vitesse de depot sec dans la couche Brownienne
    SAVE vdeptr
    REAL scavtr(nqmx - 2) ! Coefficient de lessivage
    SAVE scavtr

    CHARACTER itn

    ! nature du traceur

    logical aerosol(nqmx - 2) ! Nature du traceur
    ! ! aerosol(it) = true => aerosol 
    ! ! aerosol(it) = false => gaz 
    logical clsol(nqmx - 2) ! couche limite sol calcul\'ee
    logical radio(nqmx - 2) ! d\'ecroisssance radioactive
    save aerosol, clsol, radio

    ! convection tiedtke
    INTEGER i, k, it
    REAL delp(klon, llm)

    ! Variables liees a l'ecriture de la bande histoire physique

    ! Variables locales pour effectuer les appels en serie

    REAL d_tr(klon, llm), d_trs(klon) ! tendances de traceurs 
    REAL d_tr_cl(klon, llm, nqmx - 2) ! tendance de traceurs couche limite

    REAL d_tr_cv(klon, llm, nqmx - 2) 
    ! tendance de traceurs conv pour chq traceur

    REAL d_tr_th(klon, llm, nqmx - 2) ! la tendance des thermiques
    REAL d_tr_dec(klon, llm, 2) ! la tendance de la decroissance 
    ! ! radioactive du rn - > pb 
    REAL d_tr_lessi_impa(klon, llm, nqmx - 2) ! la tendance du lessivage 
    ! ! par impaction
    REAL d_tr_lessi_nucl(klon, llm, nqmx - 2) ! la tendance du lessivage 
    ! ! par nucleation 
    REAL flestottr(klon, llm, nqmx - 2) ! flux de lessivage 
    ! ! dans chaque couche 

    real ztra_th(klon, llm)
    integer isplit, varid

    ! Controls:
    logical:: couchelimite = .true.
    logical:: convection = .true.
    logical:: lessivage = .true.
    logical, save:: inirnpb

    !--------------------------------------

    call assert(shape(zmasse) == (/klon, llm/), "phytrac zmasse")
    call assert(shape(tr_seri) == (/klon, llm, nqmx - 2/), "phytrac tr_seri")

    if (firstcal) then
       print *, 'phytrac: pdtphys = ', pdtphys
       PRINT *, 'Frequency of tracer output: ecrit_tra = ', ecrit_tra
       inirnpb = .true.

       ! Initialisation de certaines variables pour le radon et le plomb 
       ! Initialisation du traceur dans le sol (couche limite radonique)
       trs(:, 2:) = 0.

       call nf95_inq_varid(ncid_startphy, "trs", varid)
       call nf95_get_var(ncid_startphy, varid, trs(:, 1))

       ! Initialisation de la fraction d'aerosols lessivee

       d_tr_lessi_impa = 0.
       d_tr_lessi_nucl = 0. 

       ! Initialisation de la nature des traceurs

       DO it = 1, nqmx - 2
          aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
          radio(it) = .FALSE. ! par d\'efaut pas de passage par "radiornpb"
          clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
       ENDDO

       if (nqmx >= 5) then
          call press_coefoz ! read input pressure levels for ozone coefficients
       end if
    ENDIF

    if (inirnpb) THEN
       ! Initialisation du traceur dans le sol (couche limite radonique)
       radio(1)= .true.
       radio(2)= .true.
       clsol(1)= .true.
       clsol(2)= .true.
       aerosol(2) = .TRUE. ! le Pb est un aerosol 
       call initrrnpb(pctsrf, masktr, fshtr, hsoltr, tautr, vdeptr, scavtr)
       inirnpb=.false.
    endif

    if (convection) then
       ! Calcul de l'effet de la convection
       DO it=1, nqmx - 2
          if (iflag_con == 2) then
             ! Tiedke
             CALL nflxtr(pdtphys, pmfu, pmfd, pde_u, pen_d, paprs, &
                  tr_seri(:, :, it), d_tr_cv(:, :, it))
          else if (iflag_con == 3) then
             ! Emanuel
             call cvltr(pdtphys, da, phi, mp, paprs, tr_seri(:, :, it), upwd, &
                  dnwd, d_tr_cv(:, :, it))
          endif

          DO k = 1, llm
             DO i = 1, klon
                tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cv(i, k, it)
             ENDDO
          ENDDO
          WRITE(unit=itn, fmt='(i1)') it
          CALL minmaxqfi(tr_seri(:, :, it), 0., 1e33, &
               'convection, tracer index = ' // itn)
       ENDDO
    endif

    ! Calcul de l'effet des thermiques

    do it=1, nqmx - 2
       do k=1, llm
          do i=1, klon
             d_tr_th(i, k, it)=0.
             tr_seri(i, k, it)=max(tr_seri(i, k, it), 0.)
             tr_seri(i, k, it)=min(tr_seri(i, k, it), 1e10)
          enddo
       enddo
    enddo

    if (iflag_thermals > 0) then
       nsplit=10
       DO it=1, nqmx - 2
          do isplit=1, nsplit
             ! Thermiques 
             call dqthermcell(klon, llm, pdtphys/nsplit &
                  , fm_therm, entr_therm, zmasse &
                  , tr_seri(1:klon, 1:llm, it), d_tr, ztra_th)

             do k=1, llm
                do i=1, klon
                   d_tr(i, k)=pdtphys*d_tr(i, k)/nsplit
                   d_tr_th(i, k, it)=d_tr_th(i, k, it)+d_tr(i, k)
                   tr_seri(i, k, it)=max(tr_seri(i, k, it)+d_tr(i, k), 0.)
                enddo
             enddo
          enddo
       ENDDO
    endif

    ! Calcul de l'effet de la couche limite

    if (couchelimite) then
       DO k = 1, llm
          DO i = 1, klon
             delp(i, k) = paprs(i, k)-paprs(i, k+1)
          ENDDO
       ENDDO

       ! MAF modif pour tenir compte du cas traceur
       DO it=1, nqmx - 2
          if (clsol(it)) then 
             ! couche limite avec quantite dans le sol calculee
             CALL cltracrn(it, pdtphys, yu1, yv1, coefh, t_seri, ftsol, &
                  pctsrf, tr_seri(:, :, it), trs(:, it), paprs, pplay, delp, &
                  masktr(1, it), fshtr(1, it), hsoltr(it), tautr(it), &
                  vdeptr(it), rlat, d_tr_cl(1, 1, it), d_trs)
             DO k = 1, llm
                DO i = 1, klon
                   tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
                ENDDO
             ENDDO

             trs(:, it) = trs(:, it) + d_trs
          else
             ! couche limite avec flux prescrit
             !MAF provisoire source / traceur a creer
             DO i=1, klon
                source(i) = 0. ! pas de source, pour l'instant
             ENDDO

             CALL cltrac(pdtphys, coefh, t_seri, tr_seri(:, :, it), source, &
                  paprs, pplay, delp, d_tr_cl(1, 1, it))
             DO k = 1, llm
                DO i = 1, klon
                   tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
                ENDDO
             ENDDO
          endif
       ENDDO
    endif

    ! Calcul de l'effet du puits radioactif

    ! MAF il faudrait faire une modification pour passer dans radiornpb 
    ! si radio=true
    d_tr_dec = radiornpb(tr_seri, pdtphys, tautr)
    DO it = 1, nqmx - 2
       if (radio(it)) then
          tr_seri(:, :, it) = tr_seri(:, :, it) + d_tr_dec(:, :, it)
          WRITE(unit=itn, fmt='(i1)') it
          CALL minmaxqfi(tr_seri(:, :, it), 0., 1e33, 'puits rn it='//itn)
       endif
    ENDDO

    if (nqmx >= 5) then
       ! Ozone as a tracer:
       if (mod(itap - 1, lmt_pas) == 0) then
          ! Once per day, update the coefficients for ozone chemistry:
          call regr_pr_comb_coefoz(julien)
       end if
       call o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, tr_seri(:, :, 3))
    end if

    ! Calcul de l'effet de la precipitation

    IF (lessivage) THEN
       d_tr_lessi_nucl = 0. 
       d_tr_lessi_impa = 0. 
       flestottr = 0. 

       ! tendance des aerosols nuclees et impactes 

       DO it = 1, nqmx - 2
          IF (aerosol(it)) THEN
             DO k = 1, llm
                DO i = 1, klon
                   d_tr_lessi_nucl(i, k, it) = d_tr_lessi_nucl(i, k, it) + &
                        (1 - frac_nucl(i, k))*tr_seri(i, k, it)
                   d_tr_lessi_impa(i, k, it) = d_tr_lessi_impa(i, k, it) + &
                        (1 - frac_impa(i, k))*tr_seri(i, k, it)
                ENDDO
             ENDDO
          ENDIF
       ENDDO

       ! Mises a jour des traceurs + calcul des flux de lessivage 
       ! Mise a jour due a l'impaction et a la nucleation

       DO it = 1, nqmx - 2
          IF (aerosol(it)) THEN
             DO k = 1, llm
                DO i = 1, klon
                   tr_seri(i, k, it) = tr_seri(i, k, it) * frac_impa(i, k) &
                        * frac_nucl(i, k)
                ENDDO
             ENDDO
          ENDIF
       ENDDO

       ! Flux lessivage total 

       DO it = 1, nqmx - 2
          DO k = 1, llm
             DO i = 1, klon
                flestottr(i, k, it) = flestottr(i, k, it) &
                     - (d_tr_lessi_nucl(i, k, it) + d_tr_lessi_impa(i, k, it)) &
                     * (paprs(i, k)-paprs(i, k+1)) / (RG * pdtphys)
             ENDDO
          ENDDO
       ENDDO
    ENDIF

    ! Ecriture des sorties
    call write_histrac(lessivage, itap, nid_ins)

    if (lafin) then
       call nf95_inq_varid(ncid_restartphy, "trs", varid)
       call nf95_put_var(ncid_restartphy, varid, trs(:, 1))
    endif

  contains

    subroutine write_histrac(lessivage, itap, nid_ins)

      ! From phylmd/write_histrac.h, version 1.9 2006/02/21 08:08:30

      use dimens_m, only: iim, jjm, llm
      use histsync_m, only: histsync
      use histwrite_m, only: histwrite
      use temps, only: itau_phy
      use iniadvtrac_m, only: tname
      use dimphy, only: klon
      use grid_change, only: gr_phy_write_2d
      use gr_phy_write_3d_m, only: gr_phy_write_3d

      logical, intent(in):: lessivage
      integer, intent(in):: itap ! number of calls to "physiq"
      integer, intent(in):: nid_ins

      ! Variables local to the procedure:
      integer it
      integer itau_w ! pas de temps ecriture

      !-----------------------------------------------------

      itau_w = itau_phy + itap

      CALL histwrite(nid_ins, "zmasse", itau_w, gr_phy_write_3d(zmasse))

      DO it=1, nqmx - 2
         CALL histwrite(nid_ins, tname(it+2), itau_w, &
              gr_phy_write_3d(tr_seri(:, :, it)))
         if (lessivage) THEN
            CALL histwrite(nid_ins, "fl"//tname(it+2), itau_w, &
                 gr_phy_write_3d(flestottr(:, :, it)))
         endif
         CALL histwrite(nid_ins, "d_tr_th_"//tname(it+2), itau_w, &
              gr_phy_write_3d(d_tr_th(:, :, it)))
         CALL histwrite(nid_ins, "d_tr_cv_"//tname(it+2), itau_w, &
              gr_phy_write_3d(d_tr_cv(:, :, it)))
         CALL histwrite(nid_ins, "d_tr_cl_"//tname(it+2), itau_w, &
              gr_phy_write_3d(d_tr_cl(:, :, it)))
      ENDDO

    end subroutine write_histrac

  END SUBROUTINE phytrac

end module phytrac_m
1	module phytrac_m
2
3	IMPLICIT none
4
5	private
6	public phytrac
7
8	contains
9
10	SUBROUTINE phytrac(itap, lmt_pas, julien, gmtime, firstcal, lafin, pdtphys, &
11	t_seri, paprs, pplay, pmfu, pmfd, pde_u, pen_d, coefh, fm_therm, &
12	entr_therm, yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, da, phi, &
13	mp, upwd, dnwd, tr_seri, zmasse, ncid_startphy, nid_ins)
14
15	! From phylmd/phytrac.F, version 1.15 2006/02/21 08:08:30 (SVN revision 679)
16
17	! Authors: Fr\'ed\'eric Hourdin, Abderrahmane Idelkadi, Marie-Alice
18	! Foujols, Olivia
19
20	! Objet : moniteur g\'en\'eral des tendances des traceurs
21
22	! L'appel de "phytrac" se fait avec "nqmx - 2" donc nous avons
23	! bien les vrais traceurs, sans la vapeur d'eau ni l'eau liquide.
24
25	! Modifications pour les traceurs :
26	! - uniformisation des parametrisations dans phytrac
27	! - stockage des moyennes des champs n\'ecessaires en mode traceur off-line
28
29	use abort_gcm_m, only: abort_gcm
30	use clesphys, only: ecrit_tra
31	use clesphys2, only: iflag_con
32	use cltrac_m, only: cltrac
33	use cltracrn_m, only: cltracrn
34	use ctherm, only: iflag_thermals
35	use cvltr_m, only: cvltr
36	use dimens_m, only: llm, nqmx
37	use dimphy, only: klon
38	use indicesol, only: nbsrf
39	use initrrnpb_m, only: initrrnpb
40	use minmaxqfi_m, only: minmaxqfi
41	use netcdf95, only: nf95_inq_varid, nf95_get_var, nf95_put_var
42	use nflxtr_m, only: nflxtr
43	use nr_util, only: assert
44	use o3_chem_m, only: o3_chem
45	use phyetat0_m, only: rlat
46	use phyredem0_m, only: ncid_restartphy
47	use press_coefoz_m, only: press_coefoz
48	use radiornpb_m, only: radiornpb
49	use regr_pr_comb_coefoz_m, only: regr_pr_comb_coefoz
50	use SUPHEC_M, only: rg
51
52	integer, intent(in):: itap ! number of calls to "physiq"
53	integer, intent(in):: lmt_pas ! number of time steps of "physics" per day
54	integer, intent(in):: julien !jour julien, 1 <= julien <= 360
55	real, intent(in):: gmtime ! heure de la journ\'ee en fraction de jour
56	logical, intent(in):: firstcal ! first call to "calfis"
57	logical, intent(in):: lafin ! fin de la physique
58	real, intent(in):: pdtphys ! pas d'integration pour la physique (s)
59	real, intent(in):: t_seri(klon, llm) ! temperature, in K
60
61	real, intent(in):: paprs(klon, llm+1)
62	! (pression pour chaque inter-couche, en Pa)
63
64	real, intent(in):: pplay(klon, llm)
65	! (pression pour le mileu de chaque couche, en Pa)
66
67	! convection:
68
69	REAL, intent(in):: pmfu(klon, llm) ! flux de masse dans le panache montant
70
71	REAL, intent(in):: pmfd(klon, llm)
72	! flux de masse dans le panache descendant
73
74	REAL pde_u(klon, llm) ! flux detraine dans le panache montant
75	REAL pen_d(klon, llm) ! flux entraine dans le panache descendant
76	REAL coefh(klon, llm) ! coeff melange couche limite
77
78	! thermiques:
79	real fm_therm(klon, llm+1), entr_therm(klon, llm)
80
81	! Couche limite:
82	REAL yu1(klon) ! vents au premier niveau
83	REAL yv1(klon) ! vents au premier niveau
84
85	! Arguments n\'ecessaires pour les sources et puits de traceur :
86	real, intent(in):: ftsol(klon, nbsrf) ! Temperature du sol (surf)(Kelvin)
87	real pctsrf(klon, nbsrf) ! Pourcentage de sol f(nature du sol)
88
89	! Lessivage pour le on-line
90	REAL frac_impa(klon, llm) ! fraction d'aerosols impactes
91	REAL frac_nucl(klon, llm) ! fraction d'aerosols nuclees
92
93	! Kerry Emanuel
94	real, intent(in):: da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
95	REAL, intent(in):: upwd(klon, llm) ! saturated updraft mass flux
96	REAL, intent(in):: dnwd(klon, llm) ! saturated downdraft mass flux
97
98	real, intent(inout):: tr_seri(:, :, :) ! (klon, llm, nqmx - 2)
99	! (mass fractions of tracers, excluding water, at mid-layers)
100
101	real, intent(in):: zmasse(:, :) ! (klon, llm)
102	! (column-density of mass of air in a cell, in kg m-2)
103
104	integer, intent(in):: ncid_startphy, nid_ins
105
106	! Local:
107
108	integer nsplit
109
110	! TRACEURS
111
112	! Sources et puits des traceurs:
113
114	! Pour l'instant seuls les cas du rn et du pb ont ete envisages.
115
116	REAL source(klon) ! a voir lorsque le flux est prescrit
117	!
118	! Pour la source de radon et son reservoir de sol
119
120	REAL, save:: trs(klon, nqmx - 2) ! Concentration de traceur dans le sol
121
122	REAL masktr(klon, nqmx - 2) ! Masque reservoir de sol traceur
123	! Masque de l'echange avec la surface
124	! (1 = reservoir) ou (possible => 1)
125	SAVE masktr
126	REAL fshtr(klon, nqmx - 2) ! Flux surfacique dans le reservoir de sol
127	SAVE fshtr
128	REAL hsoltr(nqmx - 2) ! Epaisseur equivalente du reservoir de sol
129	SAVE hsoltr
130	REAL tautr(nqmx - 2) ! Constante de decroissance radioactive
131	SAVE tautr
132	REAL vdeptr(nqmx - 2) ! Vitesse de depot sec dans la couche Brownienne
133	SAVE vdeptr
134	REAL scavtr(nqmx - 2) ! Coefficient de lessivage
135	SAVE scavtr
136
137	CHARACTER itn
138
139	! nature du traceur
140
141	logical aerosol(nqmx - 2) ! Nature du traceur
142	! ! aerosol(it) = true => aerosol
143	! ! aerosol(it) = false => gaz
144	logical clsol(nqmx - 2) ! couche limite sol calcul\'ee
145	logical radio(nqmx - 2) ! d\'ecroisssance radioactive
146	save aerosol, clsol, radio
147
148	! convection tiedtke
149	INTEGER i, k, it
150	REAL delp(klon, llm)
151
152	! Variables liees a l'ecriture de la bande histoire physique
153
154	! Variables locales pour effectuer les appels en serie
155
156	REAL d_tr(klon, llm), d_trs(klon) ! tendances de traceurs
157	REAL d_tr_cl(klon, llm, nqmx - 2) ! tendance de traceurs couche limite
158
159	REAL d_tr_cv(klon, llm, nqmx - 2)
160	! tendance de traceurs conv pour chq traceur
161
162	REAL d_tr_th(klon, llm, nqmx - 2) ! la tendance des thermiques
163	REAL d_tr_dec(klon, llm, 2) ! la tendance de la decroissance
164	! ! radioactive du rn - > pb
165	REAL d_tr_lessi_impa(klon, llm, nqmx - 2) ! la tendance du lessivage
166	! ! par impaction
167	REAL d_tr_lessi_nucl(klon, llm, nqmx - 2) ! la tendance du lessivage
168	! ! par nucleation
169	REAL flestottr(klon, llm, nqmx - 2) ! flux de lessivage
170	! ! dans chaque couche
171
172	real ztra_th(klon, llm)
173	integer isplit, varid
174
175	! Controls:
176	logical:: couchelimite = .true.
177	logical:: convection = .true.
178	logical:: lessivage = .true.
179	logical, save:: inirnpb
180
181	!--------------------------------------
182
183	call assert(shape(zmasse) == (/klon, llm/), "phytrac zmasse")
184	call assert(shape(tr_seri) == (/klon, llm, nqmx - 2/), "phytrac tr_seri")
185
186	if (firstcal) then
187	print *, 'phytrac: pdtphys = ', pdtphys
188	PRINT *, 'Frequency of tracer output: ecrit_tra = ', ecrit_tra
189	inirnpb = .true.
190
191	! Initialisation de certaines variables pour le radon et le plomb
192	! Initialisation du traceur dans le sol (couche limite radonique)
193	trs(:, 2:) = 0.
194
195	call nf95_inq_varid(ncid_startphy, "trs", varid)
196	call nf95_get_var(ncid_startphy, varid, trs(:, 1))
197
198	! Initialisation de la fraction d'aerosols lessivee
199
200	d_tr_lessi_impa = 0.
201	d_tr_lessi_nucl = 0.
202
203	! Initialisation de la nature des traceurs
204
205	DO it = 1, nqmx - 2
206	aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
207	radio(it) = .FALSE. ! par d\'efaut pas de passage par "radiornpb"
208	clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
209	ENDDO
210
211	if (nqmx >= 5) then
212	call press_coefoz ! read input pressure levels for ozone coefficients
213	end if
214	ENDIF
215
216	if (inirnpb) THEN
217	! Initialisation du traceur dans le sol (couche limite radonique)
218	radio(1)= .true.
219	radio(2)= .true.
220	clsol(1)= .true.
221	clsol(2)= .true.
222	aerosol(2) = .TRUE. ! le Pb est un aerosol
223	call initrrnpb(pctsrf, masktr, fshtr, hsoltr, tautr, vdeptr, scavtr)
224	inirnpb=.false.
225	endif
226
227	if (convection) then
228	! Calcul de l'effet de la convection
229	DO it=1, nqmx - 2
230	if (iflag_con == 2) then
231	! Tiedke
232	CALL nflxtr(pdtphys, pmfu, pmfd, pde_u, pen_d, paprs, &
233	tr_seri(:, :, it), d_tr_cv(:, :, it))
234	else if (iflag_con == 3) then
235	! Emanuel
236	call cvltr(pdtphys, da, phi, mp, paprs, tr_seri(:, :, it), upwd, &
237	dnwd, d_tr_cv(:, :, it))
238	endif
239
240	DO k = 1, llm
241	DO i = 1, klon
242	tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cv(i, k, it)
243	ENDDO
244	ENDDO
245	WRITE(unit=itn, fmt='(i1)') it
246	CALL minmaxqfi(tr_seri(:, :, it), 0., 1e33, &
247	'convection, tracer index = ' // itn)
248	ENDDO
249	endif
250
251	! Calcul de l'effet des thermiques
252
253	do it=1, nqmx - 2
254	do k=1, llm
255	do i=1, klon
256	d_tr_th(i, k, it)=0.
257	tr_seri(i, k, it)=max(tr_seri(i, k, it), 0.)
258	tr_seri(i, k, it)=min(tr_seri(i, k, it), 1e10)
259	enddo
260	enddo
261	enddo
262
263	if (iflag_thermals > 0) then
264	nsplit=10
265	DO it=1, nqmx - 2
266	do isplit=1, nsplit
267	! Thermiques
268	call dqthermcell(klon, llm, pdtphys/nsplit &
269	, fm_therm, entr_therm, zmasse &
270	, tr_seri(1:klon, 1:llm, it), d_tr, ztra_th)
271
272	do k=1, llm
273	do i=1, klon
274	d_tr(i, k)=pdtphys*d_tr(i, k)/nsplit
275	d_tr_th(i, k, it)=d_tr_th(i, k, it)+d_tr(i, k)
276	tr_seri(i, k, it)=max(tr_seri(i, k, it)+d_tr(i, k), 0.)
277	enddo
278	enddo
279	enddo
280	ENDDO
281	endif
282
283	! Calcul de l'effet de la couche limite
284
285	if (couchelimite) then
286	DO k = 1, llm
287	DO i = 1, klon
288	delp(i, k) = paprs(i, k)-paprs(i, k+1)
289	ENDDO
290	ENDDO
291
292	! MAF modif pour tenir compte du cas traceur
293	DO it=1, nqmx - 2
294	if (clsol(it)) then
295	! couche limite avec quantite dans le sol calculee
296	CALL cltracrn(it, pdtphys, yu1, yv1, coefh, t_seri, ftsol, &
297	pctsrf, tr_seri(:, :, it), trs(:, it), paprs, pplay, delp, &
298	masktr(1, it), fshtr(1, it), hsoltr(it), tautr(it), &
299	vdeptr(it), rlat, d_tr_cl(1, 1, it), d_trs)
300	DO k = 1, llm
301	DO i = 1, klon
302	tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
303	ENDDO
304	ENDDO
305
306	trs(:, it) = trs(:, it) + d_trs
307	else
308	! couche limite avec flux prescrit
309	!MAF provisoire source / traceur a creer
310	DO i=1, klon
311	source(i) = 0. ! pas de source, pour l'instant
312	ENDDO
313
314	CALL cltrac(pdtphys, coefh, t_seri, tr_seri(:, :, it), source, &
315	paprs, pplay, delp, d_tr_cl(1, 1, it))
316	DO k = 1, llm
317	DO i = 1, klon
318	tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
319	ENDDO
320	ENDDO
321	endif
322	ENDDO
323	endif
324
325	! Calcul de l'effet du puits radioactif
326
327	! MAF il faudrait faire une modification pour passer dans radiornpb
328	! si radio=true
329	d_tr_dec = radiornpb(tr_seri, pdtphys, tautr)
330	DO it = 1, nqmx - 2
331	if (radio(it)) then
332	tr_seri(:, :, it) = tr_seri(:, :, it) + d_tr_dec(:, :, it)
333	WRITE(unit=itn, fmt='(i1)') it
334	CALL minmaxqfi(tr_seri(:, :, it), 0., 1e33, 'puits rn it='//itn)
335	endif
336	ENDDO
337
338	if (nqmx >= 5) then
339	! Ozone as a tracer:
340	if (mod(itap - 1, lmt_pas) == 0) then
341	! Once per day, update the coefficients for ozone chemistry:
342	call regr_pr_comb_coefoz(julien)
343	end if
344	call o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, tr_seri(:, :, 3))
345	end if
346
347	! Calcul de l'effet de la precipitation
348
349	IF (lessivage) THEN
350	d_tr_lessi_nucl = 0.
351	d_tr_lessi_impa = 0.
352	flestottr = 0.
353
354	! tendance des aerosols nuclees et impactes
355
356	DO it = 1, nqmx - 2
357	IF (aerosol(it)) THEN
358	DO k = 1, llm
359	DO i = 1, klon
360	d_tr_lessi_nucl(i, k, it) = d_tr_lessi_nucl(i, k, it) + &
361	(1 - frac_nucl(i, k))*tr_seri(i, k, it)
362	d_tr_lessi_impa(i, k, it) = d_tr_lessi_impa(i, k, it) + &
363	(1 - frac_impa(i, k))*tr_seri(i, k, it)
364	ENDDO
365	ENDDO
366	ENDIF
367	ENDDO
368
369	! Mises a jour des traceurs + calcul des flux de lessivage
370	! Mise a jour due a l'impaction et a la nucleation
371
372	DO it = 1, nqmx - 2
373	IF (aerosol(it)) THEN
374	DO k = 1, llm
375	DO i = 1, klon
376	tr_seri(i, k, it) = tr_seri(i, k, it) * frac_impa(i, k) &
377	* frac_nucl(i, k)
378	ENDDO
379	ENDDO
380	ENDIF
381	ENDDO
382
383	! Flux lessivage total
384
385	DO it = 1, nqmx - 2
386	DO k = 1, llm
387	DO i = 1, klon
388	flestottr(i, k, it) = flestottr(i, k, it) &
389	- (d_tr_lessi_nucl(i, k, it) + d_tr_lessi_impa(i, k, it)) &
390	* (paprs(i, k)-paprs(i, k+1)) / (RG * pdtphys)
391	ENDDO
392	ENDDO
393	ENDDO
394	ENDIF
395
396	! Ecriture des sorties
397	call write_histrac(lessivage, itap, nid_ins)
398
399	if (lafin) then
400	call nf95_inq_varid(ncid_restartphy, "trs", varid)
401	call nf95_put_var(ncid_restartphy, varid, trs(:, 1))
402	endif
403
404	contains
405
406	subroutine write_histrac(lessivage, itap, nid_ins)
407
408	! From phylmd/write_histrac.h, version 1.9 2006/02/21 08:08:30
409
410	use dimens_m, only: iim, jjm, llm
411	use histsync_m, only: histsync
412	use histwrite_m, only: histwrite
413	use temps, only: itau_phy
414	use iniadvtrac_m, only: tname
415	use dimphy, only: klon
416	use grid_change, only: gr_phy_write_2d
417	use gr_phy_write_3d_m, only: gr_phy_write_3d
418
419	logical, intent(in):: lessivage
420	integer, intent(in):: itap ! number of calls to "physiq"
421	integer, intent(in):: nid_ins
422
423	! Variables local to the procedure:
424	integer it
425	integer itau_w ! pas de temps ecriture
426
427	!-----------------------------------------------------
428
429	itau_w = itau_phy + itap
430
431	CALL histwrite(nid_ins, "zmasse", itau_w, gr_phy_write_3d(zmasse))
432
433	DO it=1, nqmx - 2
434	CALL histwrite(nid_ins, tname(it+2), itau_w, &
435	gr_phy_write_3d(tr_seri(:, :, it)))
436	if (lessivage) THEN
437	CALL histwrite(nid_ins, "fl"//tname(it+2), itau_w, &
438	gr_phy_write_3d(flestottr(:, :, it)))
439	endif
440	CALL histwrite(nid_ins, "d_tr_th_"//tname(it+2), itau_w, &
441	gr_phy_write_3d(d_tr_th(:, :, it)))
442	CALL histwrite(nid_ins, "d_tr_cv_"//tname(it+2), itau_w, &
443	gr_phy_write_3d(d_tr_cv(:, :, it)))
444	CALL histwrite(nid_ins, "d_tr_cl_"//tname(it+2), itau_w, &
445	gr_phy_write_3d(d_tr_cl(:, :, it)))
446	ENDDO
447
448	end subroutine write_histrac
449
450	END SUBROUTINE phytrac
451
452	end module phytrac_m