Sources/phylmd/phytrac.f

module phytrac_m

  IMPLICIT none

  private
  public phytrac

contains

  SUBROUTINE phytrac(itap, lmt_pas, julien, gmtime, firstcal, lafin, 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, 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 cltracrn_m, only: cltracrn
    use ctherm, only: iflag_thermals
    use dimens_m, only: llm, nqmx
    use dimphy, only: klon
    use indicesol, only: nbsrf
    use ini_histrac_m, only: ini_histrac
    use initrrnpb_m, only: initrrnpb
    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

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

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

    ! 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, nqmx - 2) ! Concentration de radon 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
    INTEGER, save:: nid_tra

    ! 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

    ! 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 des sorties :
       call ini_histrac(nid_tra, pdtphys, nqmx - 2, 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, 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(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, 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(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

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