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

    ! 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
       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))
       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

       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 (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

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

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