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, pen_u, pde_u, &
       pen_d, pde_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édéric Hourdin, Abderrahmane Idelkadi, Marie-Alice
    ! Foujols, Olivia
    ! Objet : moniteur général 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écessaires en mode traceur off-line 

    use dimens_m, only: llm
    use indicesol, only: nbsrf
    use dimphy, only: klon, nbtr
    use clesphys, only: ecrit_tra
    use clesphys2, only: iflag_con
    use abort_gcm_m, only: abort_gcm
    use SUPHEC_M, only: rg
    use ctherm, only: iflag_thermals
    use regr_pr_comb_coefoz_m, only: regr_pr_comb_coefoz
    use phyetat0_m, only: rlat
    use o3_chem_m, only: o3_chem
    use ini_histrac_m, only: ini_histrac
    use radiornpb_m, only: radiornpb
    use minmaxqfi_m, only: minmaxqfi
    use nr_util, only: assert
    use press_coefoz_m, only: press_coefoz

    logical, intent(in):: rnpb

    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):: 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ée en fraction de jour
    real, intent(in):: pdtphys  ! pas d'integration pour la physique (s)
    real, intent(in):: t_seri(klon, llm) ! temperature, in K

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

    real, intent(in):: u(klon, llm)
    real rh(klon, llm)     ! humidite relative
    real cldliq(klon, llm) ! eau liquide nuageuse
    real cldfra(klon, llm) ! fraction nuageuse (tous les nuages)

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

    real rneb(klon, llm)   ! fraction nuageuse (grande echelle)
    real albsol(klon)  ! albedo surface

    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)

    real, intent(in):: pphis(klon)
    logical, intent(in):: firstcal ! first call to "calfis"
    logical, intent(in):: lafin ! fin de la physique

    REAL pmflxr(klon, llm+1), pmflxs(klon, llm+1)   !--lessivage convection
    REAL prfl(klon, llm+1),   psfl(klon, llm+1)     !--lessivage large-scale

    !   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 pen_u(klon, llm) ! flux entraine dans le panache montant

    !   thermiques:

    real fm_therm(klon, llm+1), entr_therm(klon, llm)

    REAL pde_u(klon, llm) ! flux detraine dans le panache montant
    REAL pen_d(klon, llm) ! flux entraine dans le panache descendant
    REAL pde_d(klon, llm) ! flux detraine dans le panache descendant
    ! 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

    !   Couche limite:

    REAL coefh(klon, llm) ! coeff melange CL
    REAL yu1(klon)        ! vents au premier niveau
    REAL yv1(klon)        ! vents au premier niveau

    !   Lessivage:

    ! pour le ON-LINE

    REAL frac_impa(klon, llm)  ! fraction d'aerosols impactes
    REAL frac_nucl(klon, llm)  ! fraction d'aerosols nuclees

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

    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ée
    logical radio(nbtr)    ! décroisssance 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équence 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éfaut 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, pen_u, pde_u, pen_d, pde_d, &
                  paprs, tr_seri(1, 1, it), d_tr_cv(1, 1, 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., 1.e33, &
               '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), 1.e10)
          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(1, 1, 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., 1.e33, 'puits rn it='//itn)
          endif
       ENDDO
    endif ! rnpb decroissance  radioactive

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