Sources/phylmd/phytrac.f

module phytrac_m

  IMPLICIT none

  private
  public phytrac

contains

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

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

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

    ! 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 param\'etrisations dans phytrac
    ! - stockage des moyennes des champs n\'ecessaires en mode traceur off-line 

    use abort_gcm_m, only: abort_gcm
    use clesphys2, only: conv_emanuel
    use cltrac_m, only: cltrac
    use cltracrn_m, only: cltracrn
    USE conf_gcm_m, ONLY: lmt_pas
    use ctherm, only: iflag_thermals
    use cvltr_m, only: cvltr
    use dimens_m, only: llm, nqmx
    use dimphy, only: klon
    use histwrite_phy_m, only: histwrite_phy
    use indicesol, only: nbsrf
    use iniadvtrac_m, only: tname
    use initrrnpb_m, only: initrrnpb
    use minmaxqfi_m, only: minmaxqfi
    use netcdf, only: NF90_FILL_float
    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
    use time_phylmdz, only: itap

    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(:, 2:) ! (klon, 2:llm) coeff melange couche limite
    real cdragh(:) ! (klon)
    real fm_therm(klon, llm+1), entr_therm(klon, llm) ! thermiques
    REAL, intent(in):: yu1(:), yv1(:) ! (klon) vent au premier niveau

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

    ! Lessivage pour le on-line
    REAL, intent(in):: frac_impa(klon, llm) ! fraction d'aerosols impactes
    REAL, intent(in):: 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

    ! 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

    logical, save:: aerosol(nqmx - 2) ! Nature du traceur
    ! ! aerosol(it) = true => aerosol 
    ! ! aerosol(it) = false => gaz 

    logical, save:: clsol(nqmx - 2) ! couche limite sol flux
    ! calcul\'ee, sinon prescrit
    logical, save:: radio(nqmx - 2) ! d\'ecroisssance radioactive

    ! 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)
    ! tendance du lessivage par impaction

    REAL d_tr_lessi_nucl(klon, llm, nqmx - 2)
    ! 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:: convection = .true.

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

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

    if (firstcal) then
       ! 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))
       if (any(trs(:, 1) == NF90_FILL_float)) call abort_gcm("phytrac", &
            "some missing values in 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

       aerosol = .FALSE. ! Tous les traceurs sont des gaz par defaut
       radio = .FALSE. ! par d\'efaut pas de passage par "radiornpb"

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

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

    if (convection) then
       ! Calcul de l'effet de la convection
       DO it=1, nqmx - 2
          if (conv_emanuel) then
             call cvltr(pdtphys, da, phi, mp, paprs, tr_seri(:, :, it), upwd, &
                  dnwd, d_tr_cv(:, :, it))
          else
             CALL nflxtr(pdtphys, pmfu, pmfd, pde_u, pen_d, paprs, &
                  tr_seri(:, :, it), 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

    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, cdragh, 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

    ! 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, paprs, pplay)
       end if
       call o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, tr_seri(:, :, 3))
    end if

    ! Calcul de l'effet de la precipitation

    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

    ! Ecriture des sorties
    CALL histwrite_phy("zmasse", zmasse)
    DO it=1, nqmx - 2
       CALL histwrite_phy(tname(it+2), tr_seri(:, :, it))
       CALL histwrite_phy("fl"//tname(it+2), flestottr(:, :, it))
       CALL histwrite_phy("d_tr_th_"//tname(it+2), d_tr_th(:, :, it))
       CALL histwrite_phy("d_tr_cv_"//tname(it+2), d_tr_cv(:, :, it))
       CALL histwrite_phy("d_tr_cl_"//tname(it+2), d_tr_cl(:, :, it))
    ENDDO

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

  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(julien, gmtime, firstcal, lafin, pdtphys, t_seri, paprs, &
11	pplay, pmfu, pmfd, pde_u, pen_d, coefh, cdragh, fm_therm, entr_therm, &
12	yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, da, phi, mp, upwd, dnwd, &
13	tr_seri, zmasse, ncid_startphy)
14
15	! From phylmd/phytrac.F, version 1.15, 2006/02/21 08:08:30 (SVN
16	! revision 679) and phylmd/write_histrac.h, version 1.9,
17	! 2006/02/21 08:08:30
18
19	! Authors: Fr\'ed\'eric Hourdin, Abderrahmane Idelkadi, Marie-Alice
20	! Foujols
21
22	! Objet : moniteur g\'en\'eral des tendances des traceurs
23
24	! L'appel de "phytrac" se fait avec "nqmx - 2" donc nous avons
25	! bien les vrais traceurs, sans la vapeur d'eau ni l'eau liquide.
26
27	! Modifications pour les traceurs :
28	! - uniformisation des param\'etrisations dans phytrac
29	! - stockage des moyennes des champs n\'ecessaires en mode traceur off-line
30
31	use abort_gcm_m, only: abort_gcm
32	use clesphys2, only: conv_emanuel
33	use cltrac_m, only: cltrac
34	use cltracrn_m, only: cltracrn
35	USE conf_gcm_m, ONLY: lmt_pas
36	use ctherm, only: iflag_thermals
37	use cvltr_m, only: cvltr
38	use dimens_m, only: llm, nqmx
39	use dimphy, only: klon
40	use histwrite_phy_m, only: histwrite_phy
41	use indicesol, only: nbsrf
42	use iniadvtrac_m, only: tname
43	use initrrnpb_m, only: initrrnpb
44	use minmaxqfi_m, only: minmaxqfi
45	use netcdf, only: NF90_FILL_float
46	use netcdf95, only: nf95_inq_varid, nf95_get_var, nf95_put_var
47	use nflxtr_m, only: nflxtr
48	use nr_util, only: assert
49	use o3_chem_m, only: o3_chem
50	use phyetat0_m, only: rlat
51	use phyredem0_m, only: ncid_restartphy
52	use press_coefoz_m, only: press_coefoz
53	use radiornpb_m, only: radiornpb
54	use regr_pr_comb_coefoz_m, only: regr_pr_comb_coefoz
55	use SUPHEC_M, only: rg
56	use time_phylmdz, only: itap
57
58	integer, intent(in):: julien !jour julien, 1 <= julien <= 360
59	real, intent(in):: gmtime ! heure de la journ\'ee en fraction de jour
60	logical, intent(in):: firstcal ! first call to "calfis"
61	logical, intent(in):: lafin ! fin de la physique
62	real, intent(in):: pdtphys ! pas d'integration pour la physique (s)
63	real, intent(in):: t_seri(klon, llm) ! temperature, in K
64
65	real, intent(in):: paprs(klon, llm+1)
66	! (pression pour chaque inter-couche, en Pa)
67
68	real, intent(in):: pplay(klon, llm)
69	! (pression pour le mileu de chaque couche, en Pa)
70
71	! convection:
72
73	REAL, intent(in):: pmfu(klon, llm) ! flux de masse dans le panache montant
74
75	REAL, intent(in):: pmfd(klon, llm)
76	! flux de masse dans le panache descendant
77
78	REAL pde_u(klon, llm) ! flux detraine dans le panache montant
79	REAL pen_d(klon, llm) ! flux entraine dans le panache descendant
80	REAL coefh(:, 2:) ! (klon, 2:llm) coeff melange couche limite
81	real cdragh(:) ! (klon)
82	real fm_therm(klon, llm+1), entr_therm(klon, llm) ! thermiques
83	REAL, intent(in):: yu1(:), yv1(:) ! (klon) vent au premier niveau
84
85	! Arguments n\'ecessaires pour les sources et puits de traceur :
86	real, intent(in):: ftsol(:, :) ! (klon, nbsrf) surface temperature (K)
87	real, intent(in):: pctsrf(klon, nbsrf) ! Pourcentage de sol f(nature du sol)
88
89	! Lessivage pour le on-line
90	REAL, intent(in):: frac_impa(klon, llm) ! fraction d'aerosols impactes
91	REAL, intent(in):: 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
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	logical, save:: aerosol(nqmx - 2) ! Nature du traceur
140	! ! aerosol(it) = true => aerosol
141	! ! aerosol(it) = false => gaz
142
143	logical, save:: clsol(nqmx - 2) ! couche limite sol flux
144	! calcul\'ee, sinon prescrit
145	logical, save:: radio(nqmx - 2) ! d\'ecroisssance radioactive
146
147	! convection tiedtke
148	INTEGER i, k, it
149	REAL delp(klon, llm)
150
151	! Variables liees a l'ecriture de la bande histoire physique
152
153	! Variables locales pour effectuer les appels en serie
154
155	REAL d_tr(klon, llm), d_trs(klon) ! tendances de traceurs
156	REAL d_tr_cl(klon, llm, nqmx - 2) ! tendance de traceurs couche limite
157
158	REAL d_tr_cv(klon, llm, nqmx - 2)
159	! tendance de traceurs conv pour chq traceur
160
161	REAL d_tr_th(klon, llm, nqmx - 2) ! la tendance des thermiques
162	REAL d_tr_dec(klon, llm, 2) ! la tendance de la decroissance
163	! ! radioactive du rn - > pb
164
165	REAL d_tr_lessi_impa(klon, llm, nqmx - 2)
166	! tendance du lessivage par impaction
167
168	REAL d_tr_lessi_nucl(klon, llm, nqmx - 2)
169	! tendance du lessivage par nucleation
170
171	REAL flestottr(klon, llm, nqmx - 2) ! flux de lessivage dans chaque couche
172
173	real ztra_th(klon, llm)
174	integer isplit, varid
175
176	! Controls:
177	logical:: convection = .true.
178
179	!--------------------------------------
180
181	call assert(shape(zmasse) == (/klon, llm/), "phytrac zmasse")
182	call assert(shape(tr_seri) == (/klon, llm, nqmx - 2/), "phytrac tr_seri")
183
184	if (firstcal) then
185	! Initialisation de certaines variables pour le radon et le plomb
186	! Initialisation du traceur dans le sol (couche limite radonique)
187	trs(:, 2:) = 0.
188
189	call nf95_inq_varid(ncid_startphy, "trs", varid)
190	call nf95_get_var(ncid_startphy, varid, trs(:, 1))
191	if (any(trs(:, 1) == NF90_FILL_float)) call abort_gcm("phytrac", &
192	"some missing values in trs(:, 1)")
193
194	! Initialisation de la fraction d'aerosols lessivee
195
196	d_tr_lessi_impa = 0.
197	d_tr_lessi_nucl = 0.
198
199	! Initialisation de la nature des traceurs
200
201	aerosol = .FALSE. ! Tous les traceurs sont des gaz par defaut
202	radio = .FALSE. ! par d\'efaut pas de passage par "radiornpb"
203
204	if (nqmx >= 5) then
205	call press_coefoz ! read input pressure levels for ozone coefficients
206	end if
207
208	! Initialisation du traceur dans le sol (couche limite radonique)
209	radio(1)= .true.
210	radio(2)= .true.
211	clsol(:2)= .true.
212	clsol(3:)= .false.
213	aerosol(2) = .TRUE. ! le Pb est un aerosol
214	call initrrnpb(pctsrf, masktr, fshtr, hsoltr, tautr, vdeptr, scavtr)
215	endif
216
217	if (convection) then
218	! Calcul de l'effet de la convection
219	DO it=1, nqmx - 2
220	if (conv_emanuel) then
221	call cvltr(pdtphys, da, phi, mp, paprs, tr_seri(:, :, it), upwd, &
222	dnwd, d_tr_cv(:, :, it))
223	else
224	CALL nflxtr(pdtphys, pmfu, pmfd, pde_u, pen_d, paprs, &
225	tr_seri(:, :, it), d_tr_cv(:, :, it))
226	endif
227
228	DO k = 1, llm
229	DO i = 1, klon
230	tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cv(i, k, it)
231	ENDDO
232	ENDDO
233	WRITE(unit=itn, fmt='(i1)') it
234	CALL minmaxqfi(tr_seri(:, :, it), 0., 1e33, &
235	'convection, tracer index = ' // itn)
236	ENDDO
237	endif
238
239	! Calcul de l'effet des thermiques
240
241	do it=1, nqmx - 2
242	do k=1, llm
243	do i=1, klon
244	d_tr_th(i, k, it)=0.
245	tr_seri(i, k, it)=max(tr_seri(i, k, it), 0.)
246	tr_seri(i, k, it)=min(tr_seri(i, k, it), 1e10)
247	enddo
248	enddo
249	enddo
250
251	if (iflag_thermals > 0) then
252	nsplit=10
253	DO it=1, nqmx - 2
254	do isplit=1, nsplit
255	! Thermiques
256	call dqthermcell(klon, llm, pdtphys/nsplit &
257	, fm_therm, entr_therm, zmasse &
258	, tr_seri(1:klon, 1:llm, it), d_tr, ztra_th)
259
260	do k=1, llm
261	do i=1, klon
262	d_tr(i, k)=pdtphys*d_tr(i, k)/nsplit
263	d_tr_th(i, k, it)=d_tr_th(i, k, it)+d_tr(i, k)
264	tr_seri(i, k, it)=max(tr_seri(i, k, it)+d_tr(i, k), 0.)
265	enddo
266	enddo
267	enddo
268	ENDDO
269	endif
270
271	! Calcul de l'effet de la couche limite
272
273	DO k = 1, llm
274	DO i = 1, klon
275	delp(i, k) = paprs(i, k)-paprs(i, k+1)
276	ENDDO
277	ENDDO
278
279	! MAF modif pour tenir compte du cas traceur
280	DO it=1, nqmx - 2
281	if (clsol(it)) then
282	! couche limite avec quantite dans le sol calculee
283	CALL cltracrn(it, pdtphys, yu1, yv1, coefh, cdragh, t_seri, ftsol, &
284	pctsrf, tr_seri(:, :, it), trs(:, it), paprs, pplay, delp, &
285	masktr(1, it), fshtr(1, it), hsoltr(it), tautr(it), &
286	vdeptr(it), rlat, d_tr_cl(1, 1, it), d_trs)
287	DO k = 1, llm
288	DO i = 1, klon
289	tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
290	ENDDO
291	ENDDO
292
293	trs(:, it) = trs(:, it) + d_trs
294	else
295	! couche limite avec flux prescrit
296	!MAF provisoire source / traceur a creer
297	DO i=1, klon
298	source(i) = 0. ! pas de source, pour l'instant
299	ENDDO
300
301	CALL cltrac(pdtphys, coefh, t_seri, tr_seri(:, :, it), source, &
302	paprs, pplay, delp, d_tr_cl(1, 1, it))
303	DO k = 1, llm
304	DO i = 1, klon
305	tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
306	ENDDO
307	ENDDO
308	endif
309	ENDDO
310
311	! Calcul de l'effet du puits radioactif
312
313	! MAF il faudrait faire une modification pour passer dans radiornpb
314	! si radio=true
315	d_tr_dec = radiornpb(tr_seri, pdtphys, tautr)
316	DO it = 1, nqmx - 2
317	if (radio(it)) then
318	tr_seri(:, :, it) = tr_seri(:, :, it) + d_tr_dec(:, :, it)
319	WRITE(unit=itn, fmt='(i1)') it
320	CALL minmaxqfi(tr_seri(:, :, it), 0., 1e33, 'puits rn it='//itn)
321	endif
322	ENDDO
323
324	if (nqmx >= 5) then
325	! Ozone as a tracer:
326	if (mod(itap - 1, lmt_pas) == 0) then
327	! Once per day, update the coefficients for ozone chemistry:
328	call regr_pr_comb_coefoz(julien, paprs, pplay)
329	end if
330	call o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, tr_seri(:, :, 3))
331	end if
332
333	! Calcul de l'effet de la precipitation
334
335	d_tr_lessi_nucl = 0.
336	d_tr_lessi_impa = 0.
337	flestottr = 0.
338
339	! tendance des aerosols nuclees et impactes
340
341	DO it = 1, nqmx - 2
342	IF (aerosol(it)) THEN
343	DO k = 1, llm
344	DO i = 1, klon
345	d_tr_lessi_nucl(i, k, it) = d_tr_lessi_nucl(i, k, it) + &
346	(1 - frac_nucl(i, k))*tr_seri(i, k, it)
347	d_tr_lessi_impa(i, k, it) = d_tr_lessi_impa(i, k, it) + &
348	(1 - frac_impa(i, k))*tr_seri(i, k, it)
349	ENDDO
350	ENDDO
351	ENDIF
352	ENDDO
353
354	! Mises a jour des traceurs + calcul des flux de lessivage
355	! Mise a jour due a l'impaction et a la nucleation
356
357	DO it = 1, nqmx - 2
358	IF (aerosol(it)) THEN
359	DO k = 1, llm
360	DO i = 1, klon
361	tr_seri(i, k, it) = tr_seri(i, k, it) * frac_impa(i, k) &
362	* frac_nucl(i, k)
363	ENDDO
364	ENDDO
365	ENDIF
366	ENDDO
367
368	! Flux lessivage total
369	DO it = 1, nqmx - 2
370	DO k = 1, llm
371	DO i = 1, klon
372	flestottr(i, k, it) = flestottr(i, k, it) &
373	- (d_tr_lessi_nucl(i, k, it) + d_tr_lessi_impa(i, k, it)) &
374	* (paprs(i, k)-paprs(i, k+1)) / (RG * pdtphys)
375	ENDDO
376	ENDDO
377	ENDDO
378
379	! Ecriture des sorties
380	CALL histwrite_phy("zmasse", zmasse)
381	DO it=1, nqmx - 2
382	CALL histwrite_phy(tname(it+2), tr_seri(:, :, it))
383	CALL histwrite_phy("fl"//tname(it+2), flestottr(:, :, it))
384	CALL histwrite_phy("d_tr_th_"//tname(it+2), d_tr_th(:, :, it))
385	CALL histwrite_phy("d_tr_cv_"//tname(it+2), d_tr_cv(:, :, it))
386	CALL histwrite_phy("d_tr_cl_"//tname(it+2), d_tr_cl(:, :, it))
387	ENDDO
388
389	if (lafin) then
390	call nf95_inq_varid(ncid_restartphy, "trs", varid)
391	call nf95_put_var(ncid_restartphy, varid, trs(:, 1))
392	endif
393
394	END SUBROUTINE phytrac
395
396	end module phytrac_m