Sources/phylmd/phytrac.f

module phytrac_m

  IMPLICIT none

  private
  public phytrac

contains

  SUBROUTINE phytrac(rnpb, itap, lmt_pas, julien, gmtime, firstcal, lafin, &
       nq_phys, pdtphys, u, t_seri, paprs, pplay, pmfu, pmfd, pde_u, pen_d, &
       coefh, fm_therm, entr_therm, yu1, yv1, ftsol, pctsrf, frac_impa, &
       frac_nucl, pphis, albsol, rh, cldfra, rneb, diafra, cldliq, pmflxr, &
       pmflxs, prfl, psfl, da, phi, mp, upwd, dnwd, tr_seri, zmasse)

    ! 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 (en nombre "nbtr", sans la vapeur d'eau
    ! ni l'eau liquide) dans "phytrac".

    ! 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 cltracrn_m, only: cltracrn
    use ctherm, only: iflag_thermals
    use dimens_m, only: llm
    use dimphy, only: klon, nbtr
    use indicesol, only: nbsrf
    use ini_histrac_m, only: ini_histrac
    use minmaxqfi_m, only: minmaxqfi
    use nflxtr_m, only: nflxtr
    use nr_util, only: assert
    use o3_chem_m, only: o3_chem
    use phyetat0_m, only: rlat
    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

    logical, intent(in):: rnpb
    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

    integer, intent(in):: nq_phys
    ! (nombre de traceurs auxquels on applique la physique)

    real, intent(in):: pdtphys ! pas d'integration pour la physique (s)
    real, intent(in):: u(klon, llm)
    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 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

    real, intent(in):: pphis(klon)
    real albsol(klon) ! albedo surface
    real rh(klon, llm) ! humidite relative
    real cldfra(klon, llm) ! fraction nuageuse (tous les nuages)
    real rneb(klon, llm) ! fraction nuageuse (grande echelle)

    real diafra(klon, llm)
    ! (fraction nuageuse (convection ou stratus artificiels))

    real cldliq(klon, llm) ! eau liquide nuageuse
    REAL pmflxr(klon, llm+1), pmflxs(klon, llm+1) !--lessivage convection
    REAL prfl(klon, llm+1), psfl(klon, llm+1) !--lessivage large-scale

    ! Kerry Emanuel
    real 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, nbtr)
    ! (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)

    ! Variables local to the procedure:

    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, nbtr) ! Concentration de radon dans le sol

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

    CHARACTER itn
    INTEGER, save:: nid_tra

    ! nature du traceur

    logical aerosol(nbtr) ! Nature du traceur
    ! ! aerosol(it) = true => aerosol 
    ! ! aerosol(it) = false => gaz 
    logical clsol(nbtr) ! couche limite sol calcul\'ee
    logical radio(nbtr) ! 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, nbtr) ! tendance de traceurs couche limite
    REAL d_tr_cv(klon, llm, nbtr) ! tendance de traceurs conv pour chq traceur
    REAL d_tr_th(klon, llm, nbtr) ! 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, nbtr) ! la tendance du lessivage 
    ! ! par impaction
    REAL d_tr_lessi_nucl(klon, llm, nbtr) ! la tendance du lessivage 
    ! ! par nucleation 
    REAL flestottr(klon, llm, nbtr) ! flux de lessivage 
    ! ! dans chaque couche 

    real ztra_th(klon, llm)
    integer isplit

    ! 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, nbtr/), "phytrac tr_seri")

    if (firstcal) then
       print *, 'phytrac: pdtphys = ', pdtphys
       PRINT *, 'Fr\'equence de sortie des traceurs : ecrit_tra = ', ecrit_tra
       if (nbtr < nq_phys) call abort_gcm('phytrac', 'nbtr < nq_phys', 1)
       inirnpb=rnpb

       ! Initialisation des sorties :
       call ini_histrac(nid_tra, pdtphys, nq_phys, lessivage)

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

       open (unit=99, file='starttrac', status='old', err=999, &
            form='formatted')
       read(unit=99, fmt=*) (trs(i, 1), i=1, klon)
999    continue
       close(unit=99)

       ! 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, nq_phys
          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 (nq_phys >= 3) 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(ftsol, pctsrf, masktr, fshtr, hsoltr, tautr, vdeptr, &
            scavtr)
       inirnpb=.false.
    endif

    if (convection) then
       ! Calcul de l'effet de la convection
       DO it=1, nq_phys
          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(1, 1, it), upwd, &
                  dnwd, d_tr_cv(1, 1, 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, nq_phys
       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, nq_phys
          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 rnpb + traceur
       DO it=1, nq_phys
          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(1, 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

             ! Traceur ds sol

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

             CALL cltrac(pdtphys, coefh, t_seri, &
                  tr_seri(1, 1, 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 ! couche limite

    ! Calcul de l'effet du puits radioactif

    ! MAF il faudrait faire une modification pour passer dans radiornpb 
    ! si radio=true mais pour l'instant radiornpb propre au cas rnpb
    if (rnpb) then
       d_tr_dec(:, :, :) = radiornpb(tr_seri, pdtphys, tautr)
       DO it=1, nq_phys
          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
    endif

    if (nq_phys >= 3) 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, nq_phys
          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, nq_phys
          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, nq_phys
          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, nq_phys, itap, nid_tra)

    if (lafin) then
       print *, "C'est la fin de la physique."
       open(unit=99, file='restarttrac', form='formatted')
       do i=1, klon
          write(unit=99, fmt=*) trs(i, 1)
       enddo
       PRINT *, 'Ecriture du fichier restarttrac'
       close(unit=99)
    endif

  contains

    subroutine write_histrac(lessivage, nq_phys, itap, nid_tra)

      ! 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: tnom
      use comgeomphy, only: airephy
      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):: nq_phys
      ! (nombre de traceurs auxquels on applique la physique)

      integer, intent(in):: itap ! number of calls to "physiq"
      integer, intent(in):: nid_tra

      ! Variables local to the procedure:
      integer it
      integer itau_w ! pas de temps ecriture
      logical, parameter:: ok_sync = .true.

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

      itau_w = itau_phy + itap

      CALL histwrite(nid_tra, "phis", itau_w, gr_phy_write_2d(pphis))
      CALL histwrite(nid_tra, "aire", itau_w, gr_phy_write_2d(airephy))
      CALL histwrite(nid_tra, "zmasse", itau_w, gr_phy_write_3d(zmasse))

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

      CALL histwrite(nid_tra, "pplay", itau_w, gr_phy_write_3d(pplay))
      CALL histwrite(nid_tra, "T", itau_w, gr_phy_write_3d(t_seri))

      if (ok_sync) then
         call histsync(nid_tra)
      endif

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