source: LMDZ6/branches/Amaury_dev/libf/phylmd/cltracrn.F90 @ 5218

!$Id $

SUBROUTINE cltracrn(itr, dtime, u1lay, v1lay, &
        cdrag, coef, t, ftsol, pctsrf, &
        tr, trs, paprs, pplay, delp, &
        masktr, fshtr, hsoltr, tautr, vdeptr, &
        lat, d_tr, d_trs)

  USE dimphy
  USE traclmdz_mod, ONLY: id_rn, id_pb
  USE indice_sol_mod
  USE lmdz_yomcst

  IMPLICIT NONE
  !======================================================================
  ! Auteur(s): Alex/LMD) date:  fev 99
  !            inspire de clqh + clvent
  ! Objet: diffusion verticale de traceurs avec quantite de traceur ds
  !        le sol ( reservoir de sol de radon )

  ! note : pour l'instant le traceur dans le sol et le flux sont
  !        calcules mais ils ne servent que de diagnostiques
  !        seule la tendance sur le traceur est sortie (d_tr)
  !---------------------------------------------------------------------
  ! Arguments:
  ! itr......input-R-  le type de traceur : id_rn(radon), id_pb(plomb)
  ! dtime....input-R-  intervalle du temps (en secondes) ~ pdtphys
  ! u1lay....input-R-  vent u de la premiere couche (m/s)
  ! v1lay....input-R-  vent v de la premiere couche (m/s)
  ! cdrag....input-R-  cdrag
  ! coef.....input-R-  le coefficient d'echange (m**2/s) l>1, valable uniquement pour k entre 2 et klev
  ! t........input-R-  temperature (K)
  ! paprs....input-R-  pression a inter-couche (Pa)
  ! pplay....input-R-  pression au milieu de couche (Pa)
  ! delp.....input-R-  epaisseur de couche (Pa)
  ! ftsol....input-R-  temperature du sol (en Kelvin)
  ! tr.......input-R-  traceurs
  ! trs......input-R-  traceurs dans le sol
  ! masktr...input-R-  Masque reservoir de sol traceur (1 = reservoir)
  ! fshtr....input-R-  Flux surfacique de production dans le sol
  ! tautr....input-R-  Constante de decroissance du traceur
  ! vdeptr...input-R-  Vitesse de depot sec dans la couche brownienne
  ! hsoltr...input-R-  Epaisseur equivalente du reservoir de sol
  ! lat......input-R-  latitude en degree
  ! d_tr.....output-R- le changement de "tr"
  ! d_trs....output-R- le changement de "trs"
  !======================================================================

  !Entrees
  INTEGER, INTENT(IN) :: itr
  REAL, INTENT(IN) :: dtime
  REAL, DIMENSION(klon), INTENT(IN) :: u1lay, v1lay
  REAL, DIMENSION(klon), INTENT(IN) :: cdrag
  REAL, DIMENSION(klon, klev), INTENT(IN) :: coef, t
  REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: ftsol, pctsrf
  REAL, DIMENSION(klon, klev), INTENT(IN) :: tr
  REAL, DIMENSION(klon), INTENT(IN) :: trs
  REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs
  REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay, delp
  REAL, DIMENSION(klon), INTENT(IN) :: masktr
  REAL, DIMENSION(klon), INTENT(IN) :: fshtr
  REAL, INTENT(IN) :: hsoltr
  REAL, INTENT(IN) :: tautr
  REAL, INTENT(IN) :: vdeptr
  REAL, DIMENSION(klon), INTENT(IN) :: lat

  !Sorties
  REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_tr
  REAL, DIMENSION(klon), INTENT(OUT) :: d_trs  ! (diagnostic) traceur ds le sol

  !Locales
  REAL, DIMENSION(klon, klev) :: flux_tr  ! (diagnostic) flux de traceur
  INTEGER :: i, k, n, l
  REAL, DIMENSION(klon) :: rotrhi
  REAL, DIMENSION(klon, klev) :: zx_coef
  REAL, DIMENSION(klon) :: zx_buf
  REAL, DIMENSION(klon, klev) :: zx_ctr
  REAL, DIMENSION(klon, klev) :: zx_dtr
  REAL, DIMENSION(klon) :: zx_trs
  REAL :: zx_a, zx_b

  REAL, DIMENSION(klon, klev) :: local_tr
  REAL, DIMENSION(klon) :: local_trs
  REAL, DIMENSION(klon) :: zts      ! champ de temperature du sol
  REAL, DIMENSION(klon) :: zx_alpha1, zx_alpha2
  !======================================================================
  !AA Pour l'instant les 4 types de surface ne sont pas pris en compte
  !AA On fabrique avec zts un champ de temperature de sol
  !AA que le pondere par la fraction de nature de sol.

  DO i = 1, klon
    zts(i) = 0.
  ENDDO

  DO n = 1, nbsrf
    DO i = 1, klon
      zts(i) = zts(i) + ftsol(i, n) * pctsrf(i, n)
    ENDDO
  ENDDO

  DO i = 1, klon
    rotrhi(i) = RD * zts(i) / hsoltr
  ENDDO

  DO k = 1, klev
    DO i = 1, klon
      local_tr(i, k) = tr(i, k)
    ENDDO
  ENDDO

  DO i = 1, klon
    local_trs(i) = trs(i)
  ENDDO
  !======================================================================
  !AA   Attention si dans clmain zx_alf1(i) = 1.0
  !AA   Il doit y avoir coherence (dc la meme chose ici)

  DO i = 1, klon
    !AA         zx_alpha1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
    zx_alpha1(i) = 1.0
    zx_alpha2(i) = 1.0 - zx_alpha1(i)
  ENDDO
  !======================================================================
  DO i = 1, klon
    zx_coef(i, 1) = cdrag(i) * (1.0 + SQRT(u1lay(i)**2 + v1lay(i)**2)) &
            * pplay(i, 1) / (RD * t(i, 1))
    zx_coef(i, 1) = zx_coef(i, 1) * dtime * RG
  ENDDO

  DO k = 2, klev
    DO i = 1, klon
      zx_coef(i, k) = coef(i, k) * RG / (pplay(i, k - 1) - pplay(i, k)) &
              * (paprs(i, k) * 2 / (t(i, k) + t(i, k - 1)) / RD)**2
      zx_coef(i, k) = zx_coef(i, k) * dtime * RG
    ENDDO
  ENDDO
  !======================================================================
  DO i = 1, klon
    zx_buf(i) = delp(i, klev) + zx_coef(i, klev)
    zx_ctr(i, klev) = local_tr(i, klev) * delp(i, klev) / zx_buf(i)
    zx_dtr(i, klev) = zx_coef(i, klev) / zx_buf(i)
  ENDDO

  DO l = klev - 1, 2, -1
    DO i = 1, klon
      zx_buf(i) = delp(i, l) + zx_coef(i, l)      &
              + zx_coef(i, l + 1) * (1. - zx_dtr(i, l + 1))

      zx_ctr(i, l) = (local_tr(i, l) * delp(i, l) &
              + zx_coef(i, l + 1) * zx_ctr(i, l + 1)) / zx_buf(i)
      zx_dtr(i, l) = zx_coef(i, l) / zx_buf(i)
    ENDDO
  ENDDO

  DO i = 1, klon
    zx_buf(i) = delp(i, 1) + zx_coef(i, 2) * (1. - zx_dtr(i, 2))  &
            + masktr(i) * zx_coef(i, 1)                        &
                    * (zx_alpha1(i) + zx_alpha2(i) * zx_dtr(i, 2))

    zx_ctr(i, 1) = (local_tr(i, 1) * delp(i, 1)                &
            + zx_ctr(i, 2)                                     &
                    * (zx_coef(i, 2)                                    &
                            - masktr(i) * zx_coef(i, 1)                        &
                                    * zx_alpha2(i))) / zx_buf(i)
    zx_dtr(i, 1) = masktr(i) * zx_coef(i, 1) / zx_buf(i)
  ENDDO
  !======================================================================
  ! Calculer d'abord local_trs nouvelle quantite dans le reservoir
  ! de sol
  !=====================================================================

  DO i = 1, klon
    !-------------------------
    ! Au dessus des continents
    !--
    ! Le pb peut se deposer partout : vdeptr = 10-3 m/s
    ! Le Rn est traiter commme une couche Brownienne puisque vdeptr = 0.
    !-------------------------------------------------------------------
    IF (NINT(masktr(i)) == 1) THEN
      zx_trs(i) = local_trs(i)
      zx_a = zx_trs(i)                                           &
              + fshtr(i) * dtime * rotrhi(i)                             &
              + rotrhi(i) * masktr(i) * zx_coef(i, 1) / RG                  &
                      * (zx_ctr(i, 1) * (zx_alpha1(i) + zx_alpha2(i) * zx_dtr(i, 2)) &
                              + zx_alpha2(i) * zx_ctr(i, 2))
      ! Pour l'instant, pour aller vite, le depot sec est traite comme une decroissance
      zx_b = 1. + rotrhi(i) * masktr(i) * zx_coef(i, 1) / RG            &
              * (1. - zx_dtr(i, 1)                                     &
                      * (zx_alpha1(i) + zx_alpha2(i) * zx_dtr(i, 2)))             &
              + dtime / tautr                                       &
              + dtime * vdeptr / hsoltr
      zx_trs(i) = zx_a / zx_b
      local_trs(i) = zx_trs(i)
    ENDIF
    !--------------------------------------------------------
    ! Si on est entre 60N et 70N on divise par 2 l'emanation
    !--------------------------------------------------------

    IF ((itr==id_rn.AND.NINT(masktr(i))==1.AND.lat(i)>=60..AND.lat(i)<=70.).OR.      &
            (itr==id_pb.AND.NINT(masktr(i))==1.AND.lat(i)>=60..AND.lat(i)<=70.)) THEN
      zx_trs(i) = local_trs(i)
      zx_a = zx_trs(i)                                           &
              + (fshtr(i) / 2.) * dtime * rotrhi(i)                        &
              + rotrhi(i) * masktr(i) * zx_coef(i, 1) / RG                  &
                      * (zx_ctr(i, 1) * (zx_alpha1(i) + zx_alpha2(i) * zx_dtr(i, 2)) &
                              + zx_alpha2(i) * zx_ctr(i, 2))

      zx_b = 1. + rotrhi(i) * masktr(i) * zx_coef(i, 1) / RG  &
              * (1. - zx_dtr(i, 1)                           &
                      * (zx_alpha1(i) + zx_alpha2(i) * zx_dtr(i, 2)))   &
              + dtime / tautr                             &
              + dtime * vdeptr / hsoltr

      zx_trs(i) = zx_a / zx_b
      local_trs(i) = zx_trs(i)
    ENDIF

    !----------------------------------------------
    ! Au dessus des oceans et aux hautes latitudes
    !--
    ! au dessous de -60S  pas d'emission de radon au dessus
    ! des oceans et des continents
    !---------------------------------------------------------------

    IF ((itr==id_rn.AND.NINT(masktr(i))==0).OR.       &
            (itr==id_rn.AND.NINT(masktr(i))==1.AND.lat(i)<-60.)) THEN
      zx_trs(i) = 0.
      local_trs(i) = 0.
    END IF
    !--
    ! au dessus de 70 N pas d'emission de radon au dessus
    ! des oceans et des continents
    !--------------------------------------------------------------
    IF ((itr==id_rn.AND.NINT(masktr(i))==0).OR.    &
            (itr==id_rn.AND.NINT(masktr(i))==1.AND.lat(i)>70.)) THEN
      zx_trs(i) = 0.
      local_trs(i) = 0.
    END IF
    !---------------------------------------------
    ! Au dessus des oceans la source est nulle
    !--------------------------------------------

    IF (itr==id_rn.AND.NINT(masktr(i))==0) THEN
      zx_trs(i) = 0.
      local_trs(i) = 0.
    END IF

  ENDDO    ! sur le i=1,klon

  !======================================================================
  ! Une fois on a zx_trs, on peut faire l'iteration
  !======================================================================

  DO i = 1, klon
    local_tr(i, 1) = zx_ctr(i, 1) + zx_dtr(i, 1) * zx_trs(i)
  ENDDO
  DO l = 2, klev
    DO i = 1, klon
      local_tr(i, l) = zx_ctr(i, l) + zx_dtr(i, l) * local_tr(i, l - 1)
    ENDDO
  ENDDO
  !======================================================================
  ! Calcul du flux de traceur (flux_tr): UA/(m**2 s)
  !======================================================================
  DO i = 1, klon
    flux_tr(i, 1) = masktr(i) * zx_coef(i, 1) / RG                      &
            * (zx_alpha1(i) * local_tr(i, 1) + zx_alpha2(i) * local_tr(i, 2) &
                    - zx_trs(i)) / dtime
  ENDDO
  DO l = 2, klev
    DO i = 1, klon
      flux_tr(i, l) = zx_coef(i, l) / RG                    &
              * (local_tr(i, l) - local_tr(i, l - 1)) / dtime
    ENDDO
  ENDDO
  !======================================================================
  ! Calcul des tendances du traceur ds le sol et dans l'atmosphere
  !======================================================================
  DO l = 1, klev
    DO i = 1, klon
      d_tr(i, l) = local_tr(i, l) - tr(i, l)
    ENDDO
  ENDDO
  DO i = 1, klon
    d_trs(i) = local_trs(i) - trs(i)
  ENDDO

END SUBROUTINE cltracrn