phylmd/Interface_surf/pbl_surface.f

module pbl_surface_m

  IMPLICIT NONE

contains

  SUBROUTINE pbl_surface(pctsrf, t, q, u, v, julien, mu0, ftsol, cdmmax, &
       cdhmax, ftsoil, qsol, paprs, pplay, fsnow, qsurf, falbe, fluxlat, &
       rain_fall, snow_fall, fsolsw, fsollw, frugs, agesno, rugoro, d_t, d_q, &
       d_u, d_v, d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2, &
       dflux_t, dflux_q, coefh, t2m, q2m, u10m_srf, v10m_srf, pblh, capcl, &
       oliqcl, cteicl, pblt, therm, plcl, fqcalving, ffonte, run_off_lic_0)

    ! From phylmd/clmain.F, version 1.6, 2005/11/16 14:47:19
    ! Author: Z. X. Li (LMD/CNRS)
    ! Date: Aug. 18th, 1993
    ! Objet : interface de couche limite (diffusion verticale)

    ! Tout ce qui a trait aux traceurs est dans "phytrac". Le calcul
    ! de la couche limite pour les traceurs se fait avec "cltrac" et
    ! ne tient pas compte de la diff\'erentiation des sous-fractions
    ! de sol.

    use cdrag_m, only: cdrag
    use clqh_m, only: clqh
    use clvent_m, only: clvent
    use coef_diff_turb_m, only: coef_diff_turb
    USE conf_gcm_m, ONLY: lmt_pas
    USE conf_phys_m, ONLY: iflag_pbl
    USE dimphy, ONLY: klev, klon
    USE dimsoil, ONLY: nsoilmx
    use hbtm_m, only: hbtm
    USE histwrite_phy_m, ONLY: histwrite_phy
    USE indicesol, ONLY: epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
    USE interfoce_lim_m, ONLY: interfoce_lim
    use phyetat0_m, only: zmasq
    use stdlevvar_m, only: stdlevvar
    USE suphec_m, ONLY: rd, rg
    use time_phylmdz, only: itap

    REAL, INTENT(inout):: pctsrf(klon, nbsrf)
    ! tableau des pourcentages de surface de chaque maille

    REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
    REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg / kg)
    REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
    INTEGER, INTENT(IN):: julien ! jour de l'annee en cours
    REAL, intent(in):: mu0(klon) ! cosinus de l'angle solaire zenithal     
    REAL, INTENT(IN):: ftsol(:, :) ! (klon, nbsrf) temp\'erature du sol (en K)
    REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh

    REAL, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)
    ! soil temperature of surface fraction

    REAL, INTENT(inout):: qsol(:) ! (klon)
    ! column-density of water in soil, in kg m-2

    REAL, INTENT(IN):: paprs(klon, klev + 1) ! pression a intercouche (Pa)
    REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
    REAL, INTENT(inout):: fsnow(:, :) ! (klon, nbsrf) \'epaisseur neigeuse
    REAL, INTENT(inout):: qsurf(klon, nbsrf)
    REAL, intent(inout):: falbe(klon, nbsrf)
    REAL, intent(out):: fluxlat(:, :) ! (klon, nbsrf)

    REAL, intent(in):: rain_fall(klon)
    ! liquid water mass flux (kg / m2 / s), positive down

    REAL, intent(in):: snow_fall(klon)
    ! solid water mass flux (kg / m2 / s), positive down

    REAL, INTENT(IN):: fsolsw(klon, nbsrf), fsollw(klon, nbsrf)
    REAL, intent(inout):: frugs(klon, nbsrf) ! longueur de rugosit\'e (en m)
    real agesno(klon, nbsrf)
    REAL, INTENT(IN):: rugoro(klon)

    REAL, intent(out):: d_t(:, :), d_q(:, :) ! (klon, klev)
    ! changement pour t et q

    REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
    ! changement pour "u" et "v"

    REAL, intent(out):: d_ts(:, :) ! (klon, nbsrf) variation of ftsol

    REAL, intent(out):: flux_t(klon, nbsrf)
    ! flux de chaleur sensible (c_p T) (W / m2) (orientation positive
    ! vers le bas) à la surface

    REAL, intent(out):: flux_q(klon, nbsrf) 
    ! flux de vapeur d'eau (kg / m2 / s) à la surface

    REAL, intent(out):: flux_u(klon, nbsrf), flux_v(klon, nbsrf)
    ! tension du vent (flux turbulent de vent) à la surface, en Pa

    REAL, INTENT(out):: cdragh(klon), cdragm(klon)
    real q2(klon, klev + 1, nbsrf)

    ! Ocean slab:
    REAL, INTENT(out):: dflux_t(klon) ! derive du flux sensible
    REAL, INTENT(out):: dflux_q(klon) ! derive du flux latent

    REAL, intent(out):: coefh(:, 2:) ! (klon, 2:klev)
    ! Pour pouvoir extraire les coefficients d'\'echange, le champ
    ! "coefh" a \'et\'e cr\'e\'e. Nous avons moyenn\'e les valeurs de
    ! ce champ sur les quatre sous-surfaces du mod\`ele.

    REAL, INTENT(inout):: t2m(klon, nbsrf), q2m(klon, nbsrf)

    REAL, INTENT(inout):: u10m_srf(:, :), v10m_srf(:, :) ! (klon, nbsrf)
    ! composantes du vent \`a 10m sans spirale d'Ekman

    ! Ionela Musat. Cf. Anne Mathieu : planetary boundary layer, hbtm.
    ! Comme les autres diagnostics on cumule dans physiq ce qui permet
    ! de sortir les grandeurs par sous-surface.
    REAL pblh(klon, nbsrf) ! height of planetary boundary layer
    REAL capcl(klon, nbsrf)
    REAL oliqcl(klon, nbsrf)
    REAL cteicl(klon, nbsrf)
    REAL, INTENT(inout):: pblt(klon, nbsrf) ! T au nveau HCL
    REAL therm(klon, nbsrf)
    REAL plcl(klon, nbsrf)

    REAL, intent(out):: fqcalving(klon, nbsrf)
    ! flux d'eau "perdue" par la surface et necessaire pour limiter la
    ! hauteur de neige, en kg / m2 / s

    real ffonte(klon, nbsrf) ! flux thermique utilise pour fondre la neige
    REAL, intent(inout):: run_off_lic_0(:) ! (klon)

    ! Local:

    ! la nouvelle repartition des surfaces sortie de l'interface
    REAL, save:: pctsrf_new_oce(klon)
    REAL, save:: pctsrf_new_sic(klon)

    REAL y_fqcalving(klon), y_ffonte(klon)
    real y_run_off_lic_0(klon), y_run_off_lic(klon)
    REAL run_off_lic(klon) ! ruissellement total
    REAL rugmer(klon)
    REAL ytsoil(klon, nsoilmx)
    REAL yts(klon), ypct(klon), yz0_new(klon)
    real yrugos(klon) ! longueur de rugosite (en m)
    REAL yalb(klon)
    REAL snow(klon), yqsurf(klon), yagesno(klon)
    real yqsol(klon) ! column-density of water in soil, in kg m-2
    REAL yrain_fall(klon) ! liquid water mass flux (kg / m2 / s), positive down
    REAL ysnow_fall(klon) ! solid water mass flux (kg / m2 / s), positive down
    REAL yrugm(klon), yrads(klon), yrugoro(klon)
    REAL yfluxlat(klon)
    REAL y_d_ts(klon)
    REAL y_d_t(klon, klev), y_d_q(klon, klev)
    REAL y_d_u(klon, klev), y_d_v(klon, klev)
    REAL y_flux_t(klon), y_flux_q(klon)
    REAL y_flux_u(klon), y_flux_v(klon)
    REAL y_dflux_t(klon), y_dflux_q(klon)
    REAL ycoefh(klon, 2:klev), ycoefm(klon, 2:klev)
    real ycdragh(klon), ycdragm(klon)
    REAL yu(klon, klev), yv(klon, klev)
    REAL yt(klon, klev), yq(klon, klev)
    REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)
    REAL yq2(klon, klev + 1)
    REAL delp(klon, klev)
    INTEGER i, k, nsrf
    INTEGER ni(klon), knon, j

    REAL pctsrf_pot(klon, nbsrf)
    ! "pourcentage potentiel" pour tenir compte des \'eventuelles
    ! apparitions ou disparitions de la glace de mer

    REAL yt2m(klon), yq2m(klon), wind10m(klon)
    REAL ustar(klon)

    REAL yt10m(klon), yq10m(klon)
    REAL ypblh(klon)
    REAL ylcl(klon)
    REAL ycapcl(klon)
    REAL yoliqcl(klon)
    REAL ycteicl(klon)
    REAL ypblt(klon)
    REAL ytherm(klon)
    REAL u1(klon), v1(klon)
    REAL tair1(klon), qair1(klon), tairsol(klon)
    REAL psfce(klon), patm(klon)
    REAL zgeo1(klon)
    REAL rugo1(klon)
    REAL zgeop(klon, klev)

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

    ytherm = 0.

    DO k = 1, klev ! epaisseur de couche
       DO i = 1, klon
          delp(i, k) = paprs(i, k) - paprs(i, k + 1)
       END DO
    END DO

    ! Initialization:
    rugmer = 0.
    cdragh = 0.
    cdragm = 0.
    dflux_t = 0.
    dflux_q = 0.
    ypct = 0.
    yrugos = 0.
    ypaprs = 0.
    ypplay = 0.
    ydelp = 0.
    yrugoro = 0.
    d_ts = 0.
    flux_t = 0.
    flux_q = 0.
    flux_u = 0.
    flux_v = 0.
    fluxlat = 0.
    d_t = 0.
    d_q = 0.
    d_u = 0.
    d_v = 0.
    coefh = 0.
    fqcalving = 0.
    run_off_lic = 0.

    ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
    ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
    ! (\`a affiner).

    pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
    pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
    pctsrf_pot(:, is_oce) = 1. - zmasq
    pctsrf_pot(:, is_sic) = 1. - zmasq

    ! Tester si c'est le moment de lire le fichier:
    if (mod(itap - 1, lmt_pas) == 0) then
       CALL interfoce_lim(julien, pctsrf_new_oce, pctsrf_new_sic)
    endif

    ! Boucler sur toutes les sous-fractions du sol:

    loop_surface: DO nsrf = 1, nbsrf
       ! Define ni and knon:
       
       ni = 0
       knon = 0

       DO i = 1, klon
          ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
          ! "potentielles"
          IF (pctsrf_pot(i, nsrf) > epsfra) THEN
             knon = knon + 1
             ni(knon) = i
          END IF
       END DO

       if_knon: IF (knon /= 0) then
          DO j = 1, knon
             i = ni(j)
             ypct(j) = pctsrf(i, nsrf)
             yts(j) = ftsol(i, nsrf)
             snow(j) = fsnow(i, nsrf)
             yqsurf(j) = qsurf(i, nsrf)
             yalb(j) = falbe(i, nsrf)
             yrain_fall(j) = rain_fall(i)
             ysnow_fall(j) = snow_fall(i)
             yagesno(j) = agesno(i, nsrf)
             yrugos(j) = frugs(i, nsrf)
             yrugoro(j) = rugoro(i)
             yrads(j) = fsolsw(i, nsrf) + fsollw(i, nsrf)
             ypaprs(j, klev + 1) = paprs(i, klev + 1)
             y_run_off_lic_0(j) = run_off_lic_0(i)
          END DO

          ! For continent, copy soil water content
          IF (nsrf == is_ter) yqsol(:knon) = qsol(ni(:knon))

          ytsoil(:knon, :) = ftsoil(ni(:knon), :, nsrf)

          DO k = 1, klev
             DO j = 1, knon
                i = ni(j)
                ypaprs(j, k) = paprs(i, k)
                ypplay(j, k) = pplay(i, k)
                ydelp(j, k) = delp(i, k)
                yu(j, k) = u(i, k)
                yv(j, k) = v(i, k)
                yt(j, k) = t(i, k)
                yq(j, k) = q(i, k)
             END DO
          END DO

          ! Calculer les géopotentiels de chaque couche:

          zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
               + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))

          DO k = 2, klev
             zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
                  * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
                  * (ypplay(:knon, k - 1) - ypplay(:knon, k))
          ENDDO

          CALL cdrag(nsrf, sqrt(yu(:knon, 1)**2 + yv(:knon, 1)**2), &
               yt(:knon, 1), yq(:knon, 1), zgeop(:knon, 1), ypaprs(:knon, 1), &
               yts(:knon), yqsurf(:knon), yrugos(:knon), ycdragm(:knon), &
               ycdragh(:knon)) 

          IF (iflag_pbl == 1) THEN
             ycdragm(:knon) = max(ycdragm(:knon), 0.)
             ycdragh(:knon) = max(ycdragh(:knon), 0.)
          end IF

          ! on met un seuil pour ycdragm et ycdragh
          IF (nsrf == is_oce) THEN
             ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
             ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
          END IF

          IF (iflag_pbl >= 6) yq2(:knon, :) = q2(ni(:knon), :, nsrf)
          call coef_diff_turb(nsrf, ni(:knon), ypaprs(:knon, :), &
               ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
               yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
               ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))
          
          CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
               ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
               ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
               y_flux_u(:knon))
          CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
               ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
               ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
               y_flux_v(:knon))

          CALL clqh(julien, nsrf, ni(:knon), ytsoil(:knon, :), yqsol(:knon), &
               mu0(ni(:knon)), yrugos(:knon), yrugoro(:knon), yu(:knon, 1), &
               yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), yt(:knon, :), &
               yq(:knon, :), yts(:knon), ypaprs(:knon, :), ypplay(:knon, :), &
               ydelp(:knon, :), yrads(:knon), yalb(:knon), snow(:knon), &
               yqsurf(:knon), yrain_fall(:knon), ysnow_fall(:knon), &
               yfluxlat(:knon), pctsrf_new_sic(ni(:knon)), yagesno(:knon), &
               y_d_t(:knon, :), y_d_q(:knon, :), y_d_ts(:knon), &
               yz0_new(:knon), y_flux_t(:knon), y_flux_q(:knon), &
               y_dflux_t(:knon), y_dflux_q(:knon), y_fqcalving(:knon), &
               y_ffonte(:knon), y_run_off_lic_0(:knon), y_run_off_lic(:knon))

          ! calculer la longueur de rugosite sur ocean

          yrugm = 0.

          IF (nsrf == is_oce) THEN
             DO j = 1, knon
                yrugm(j) = 0.018 * ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2) &
                     / rg + 0.11 * 14E-6 &
                     / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2))
                yrugm(j) = max(1.5E-05, yrugm(j))
             END DO
          END IF

          DO k = 1, klev
             DO j = 1, knon
                i = ni(j)
                y_d_t(j, k) = y_d_t(j, k) * ypct(j)
                y_d_q(j, k) = y_d_q(j, k) * ypct(j)
                y_d_u(j, k) = y_d_u(j, k) * ypct(j)
                y_d_v(j, k) = y_d_v(j, k) * ypct(j)
             END DO
          END DO

          flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
          flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
          flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
          flux_v(ni(:knon), nsrf) = y_flux_v(:knon)

          falbe(:, nsrf) = 0.
          fsnow(:, nsrf) = 0.
          qsurf(:, nsrf) = 0.
          frugs(:, nsrf) = 0.
          DO j = 1, knon
             i = ni(j)
             d_ts(i, nsrf) = y_d_ts(j)
             falbe(i, nsrf) = yalb(j)
             fsnow(i, nsrf) = snow(j)
             qsurf(i, nsrf) = yqsurf(j)
             frugs(i, nsrf) = yz0_new(j)
             fluxlat(i, nsrf) = yfluxlat(j)
             IF (nsrf == is_oce) THEN
                rugmer(i) = yrugm(j)
                frugs(i, nsrf) = yrugm(j)
             END IF
             agesno(i, nsrf) = yagesno(j)
             fqcalving(i, nsrf) = y_fqcalving(j)
             ffonte(i, nsrf) = y_ffonte(j)
             cdragh(i) = cdragh(i) + ycdragh(j) * ypct(j)
             cdragm(i) = cdragm(i) + ycdragm(j) * ypct(j)
             dflux_t(i) = dflux_t(i) + y_dflux_t(j) * ypct(j)
             dflux_q(i) = dflux_q(i) + y_dflux_q(j) * ypct(j)
          END DO
          IF (nsrf == is_ter) THEN
             qsol(ni(:knon)) = yqsol(:knon)
          else IF (nsrf == is_lic) THEN
             DO j = 1, knon
                i = ni(j)
                run_off_lic_0(i) = y_run_off_lic_0(j)
                run_off_lic(i) = y_run_off_lic(j)
             END DO
          END IF

          ftsoil(:, :, nsrf) = 0.
          ftsoil(ni(:knon), :, nsrf) = ytsoil(:knon, :)

          DO j = 1, knon
             i = ni(j)
             DO k = 1, klev
                d_t(i, k) = d_t(i, k) + y_d_t(j, k)
                d_q(i, k) = d_q(i, k) + y_d_q(j, k)
                d_u(i, k) = d_u(i, k) + y_d_u(j, k)
                d_v(i, k) = d_v(i, k) + y_d_v(j, k)
             END DO
          END DO

          forall (k = 2:klev) coefh(ni(:knon), k) &
               = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypct(:knon)

          ! diagnostic t, q a 2m et u, v a 10m

          DO j = 1, knon
             i = ni(j)
             u1(j) = yu(j, 1) + y_d_u(j, 1)
             v1(j) = yv(j, 1) + y_d_v(j, 1)
             tair1(j) = yt(j, 1) + y_d_t(j, 1)
             qair1(j) = yq(j, 1) + y_d_q(j, 1)
             zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
                  1))) * (ypaprs(j, 1)-ypplay(j, 1))
             tairsol(j) = yts(j) + y_d_ts(j)
             rugo1(j) = yrugos(j)
             IF (nsrf == is_oce) THEN
                rugo1(j) = frugs(i, nsrf)
             END IF
             psfce(j) = ypaprs(j, 1)
             patm(j) = ypplay(j, 1)
          END DO

          CALL stdlevvar(nsrf, u1(:knon), v1(:knon), tair1(:knon), qair1, &
               zgeo1, tairsol, yqsurf(:knon), rugo1, psfce, patm, yt2m, yq2m, &
               yt10m, yq10m, wind10m(:knon), ustar(:knon))

          DO j = 1, knon
             i = ni(j)
             t2m(i, nsrf) = yt2m(j)
             q2m(i, nsrf) = yq2m(j)

             u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
                  / sqrt(u1(j)**2 + v1(j)**2)
             v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
                  / sqrt(u1(j)**2 + v1(j)**2)
          END DO

          CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
               y_flux_q(:knon), yu(:knon, :), yv(:knon, :), yt(:knon, :), &
               yq(:knon, :), ypblh(:knon), ycapcl, yoliqcl, ycteicl, ypblt, &
               ytherm, ylcl)

          DO j = 1, knon
             i = ni(j)
             pblh(i, nsrf) = ypblh(j)
             plcl(i, nsrf) = ylcl(j)
             capcl(i, nsrf) = ycapcl(j)
             oliqcl(i, nsrf) = yoliqcl(j)
             cteicl(i, nsrf) = ycteicl(j)
             pblt(i, nsrf) = ypblt(j)
             therm(i, nsrf) = ytherm(j)
          END DO

          IF (iflag_pbl >= 6) q2(ni(:knon), :, nsrf) = yq2(:knon, :)
       else
          fsnow(:, nsrf) = 0.
       end IF if_knon
    END DO loop_surface

    ! On utilise les nouvelles surfaces
    frugs(:, is_oce) = rugmer
    pctsrf(:, is_oce) = pctsrf_new_oce
    pctsrf(:, is_sic) = pctsrf_new_sic

    CALL histwrite_phy("run_off_lic", run_off_lic)

  END SUBROUTINE pbl_surface

end module pbl_surface_m
1	module pbl_surface_m
2
3	IMPLICIT NONE
4
5	contains
6
7	SUBROUTINE pbl_surface(pctsrf, t, q, u, v, julien, mu0, ftsol, cdmmax, &
8	cdhmax, ftsoil, qsol, paprs, pplay, fsnow, qsurf, falbe, fluxlat, &
9	rain_fall, snow_fall, fsolsw, fsollw, frugs, agesno, rugoro, d_t, d_q, &
10	d_u, d_v, d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2, &
11	dflux_t, dflux_q, coefh, t2m, q2m, u10m_srf, v10m_srf, pblh, capcl, &
12	oliqcl, cteicl, pblt, therm, plcl, fqcalving, ffonte, run_off_lic_0)
13
14	! From phylmd/clmain.F, version 1.6, 2005/11/16 14:47:19
15	! Author: Z. X. Li (LMD/CNRS)
16	! Date: Aug. 18th, 1993
17	! Objet : interface de couche limite (diffusion verticale)
18
19	! Tout ce qui a trait aux traceurs est dans "phytrac". Le calcul
20	! de la couche limite pour les traceurs se fait avec "cltrac" et
21	! ne tient pas compte de la diff\'erentiation des sous-fractions
22	! de sol.
23
24	use cdrag_m, only: cdrag
25	use clqh_m, only: clqh
26	use clvent_m, only: clvent
27	use coef_diff_turb_m, only: coef_diff_turb
28	USE conf_gcm_m, ONLY: lmt_pas
29	USE conf_phys_m, ONLY: iflag_pbl
30	USE dimphy, ONLY: klev, klon
31	USE dimsoil, ONLY: nsoilmx
32	use hbtm_m, only: hbtm
33	USE histwrite_phy_m, ONLY: histwrite_phy
34	USE indicesol, ONLY: epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
35	USE interfoce_lim_m, ONLY: interfoce_lim
36	use phyetat0_m, only: zmasq
37	use stdlevvar_m, only: stdlevvar
38	USE suphec_m, ONLY: rd, rg
39	use time_phylmdz, only: itap
40
41	REAL, INTENT(inout):: pctsrf(klon, nbsrf)
42	! tableau des pourcentages de surface de chaque maille
43
44	REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
45	REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg / kg)
46	REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
47	INTEGER, INTENT(IN):: julien ! jour de l'annee en cours
48	REAL, intent(in):: mu0(klon) ! cosinus de l'angle solaire zenithal
49	REAL, INTENT(IN):: ftsol(:, :) ! (klon, nbsrf) temp\'erature du sol (en K)
50	REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh
51
52	REAL, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)
53	! soil temperature of surface fraction
54
55	REAL, INTENT(inout):: qsol(:) ! (klon)
56	! column-density of water in soil, in kg m-2
57
58	REAL, INTENT(IN):: paprs(klon, klev + 1) ! pression a intercouche (Pa)
59	REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
60	REAL, INTENT(inout):: fsnow(:, :) ! (klon, nbsrf) \'epaisseur neigeuse
61	REAL, INTENT(inout):: qsurf(klon, nbsrf)
62	REAL, intent(inout):: falbe(klon, nbsrf)
63	REAL, intent(out):: fluxlat(:, :) ! (klon, nbsrf)
64
65	REAL, intent(in):: rain_fall(klon)
66	! liquid water mass flux (kg / m2 / s), positive down
67
68	REAL, intent(in):: snow_fall(klon)
69	! solid water mass flux (kg / m2 / s), positive down
70
71	REAL, INTENT(IN):: fsolsw(klon, nbsrf), fsollw(klon, nbsrf)
72	REAL, intent(inout):: frugs(klon, nbsrf) ! longueur de rugosit\'e (en m)
73	real agesno(klon, nbsrf)
74	REAL, INTENT(IN):: rugoro(klon)
75
76	REAL, intent(out):: d_t(:, :), d_q(:, :) ! (klon, klev)
77	! changement pour t et q
78
79	REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
80	! changement pour "u" et "v"
81
82	REAL, intent(out):: d_ts(:, :) ! (klon, nbsrf) variation of ftsol
83
84	REAL, intent(out):: flux_t(klon, nbsrf)
85	! flux de chaleur sensible (c_p T) (W / m2) (orientation positive
86	! vers le bas) à la surface
87
88	REAL, intent(out):: flux_q(klon, nbsrf)
89	! flux de vapeur d'eau (kg / m2 / s) à la surface
90
91	REAL, intent(out):: flux_u(klon, nbsrf), flux_v(klon, nbsrf)
92	! tension du vent (flux turbulent de vent) à la surface, en Pa
93
94	REAL, INTENT(out):: cdragh(klon), cdragm(klon)
95	real q2(klon, klev + 1, nbsrf)
96
97	! Ocean slab:
98	REAL, INTENT(out):: dflux_t(klon) ! derive du flux sensible
99	REAL, INTENT(out):: dflux_q(klon) ! derive du flux latent
100
101	REAL, intent(out):: coefh(:, 2:) ! (klon, 2:klev)
102	! Pour pouvoir extraire les coefficients d'\'echange, le champ
103	! "coefh" a \'et\'e cr\'e\'e. Nous avons moyenn\'e les valeurs de
104	! ce champ sur les quatre sous-surfaces du mod\`ele.
105
106	REAL, INTENT(inout):: t2m(klon, nbsrf), q2m(klon, nbsrf)
107
108	REAL, INTENT(inout):: u10m_srf(:, :), v10m_srf(:, :) ! (klon, nbsrf)
109	! composantes du vent \`a 10m sans spirale d'Ekman
110
111	! Ionela Musat. Cf. Anne Mathieu : planetary boundary layer, hbtm.
112	! Comme les autres diagnostics on cumule dans physiq ce qui permet
113	! de sortir les grandeurs par sous-surface.
114	REAL pblh(klon, nbsrf) ! height of planetary boundary layer
115	REAL capcl(klon, nbsrf)
116	REAL oliqcl(klon, nbsrf)
117	REAL cteicl(klon, nbsrf)
118	REAL, INTENT(inout):: pblt(klon, nbsrf) ! T au nveau HCL
119	REAL therm(klon, nbsrf)
120	REAL plcl(klon, nbsrf)
121
122	REAL, intent(out):: fqcalving(klon, nbsrf)
123	! flux d'eau "perdue" par la surface et necessaire pour limiter la
124	! hauteur de neige, en kg / m2 / s
125
126	real ffonte(klon, nbsrf) ! flux thermique utilise pour fondre la neige
127	REAL, intent(inout):: run_off_lic_0(:) ! (klon)
128
129	! Local:
130
131	! la nouvelle repartition des surfaces sortie de l'interface
132	REAL, save:: pctsrf_new_oce(klon)
133	REAL, save:: pctsrf_new_sic(klon)
134
135	REAL y_fqcalving(klon), y_ffonte(klon)
136	real y_run_off_lic_0(klon), y_run_off_lic(klon)
137	REAL run_off_lic(klon) ! ruissellement total
138	REAL rugmer(klon)
139	REAL ytsoil(klon, nsoilmx)
140	REAL yts(klon), ypct(klon), yz0_new(klon)
141	real yrugos(klon) ! longueur de rugosite (en m)
142	REAL yalb(klon)
143	REAL snow(klon), yqsurf(klon), yagesno(klon)
144	real yqsol(klon) ! column-density of water in soil, in kg m-2
145	REAL yrain_fall(klon) ! liquid water mass flux (kg / m2 / s), positive down
146	REAL ysnow_fall(klon) ! solid water mass flux (kg / m2 / s), positive down
147	REAL yrugm(klon), yrads(klon), yrugoro(klon)
148	REAL yfluxlat(klon)
149	REAL y_d_ts(klon)
150	REAL y_d_t(klon, klev), y_d_q(klon, klev)
151	REAL y_d_u(klon, klev), y_d_v(klon, klev)
152	REAL y_flux_t(klon), y_flux_q(klon)
153	REAL y_flux_u(klon), y_flux_v(klon)
154	REAL y_dflux_t(klon), y_dflux_q(klon)
155	REAL ycoefh(klon, 2:klev), ycoefm(klon, 2:klev)
156	real ycdragh(klon), ycdragm(klon)
157	REAL yu(klon, klev), yv(klon, klev)
158	REAL yt(klon, klev), yq(klon, klev)
159	REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)
160	REAL yq2(klon, klev + 1)
161	REAL delp(klon, klev)
162	INTEGER i, k, nsrf
163	INTEGER ni(klon), knon, j
164
165	REAL pctsrf_pot(klon, nbsrf)
166	! "pourcentage potentiel" pour tenir compte des \'eventuelles
167	! apparitions ou disparitions de la glace de mer
168
169	REAL yt2m(klon), yq2m(klon), wind10m(klon)
170	REAL ustar(klon)
171
172	REAL yt10m(klon), yq10m(klon)
173	REAL ypblh(klon)
174	REAL ylcl(klon)
175	REAL ycapcl(klon)
176	REAL yoliqcl(klon)
177	REAL ycteicl(klon)
178	REAL ypblt(klon)
179	REAL ytherm(klon)
180	REAL u1(klon), v1(klon)
181	REAL tair1(klon), qair1(klon), tairsol(klon)
182	REAL psfce(klon), patm(klon)
183	REAL zgeo1(klon)
184	REAL rugo1(klon)
185	REAL zgeop(klon, klev)
186
187	!------------------------------------------------------------
188
189	ytherm = 0.
190
191	DO k = 1, klev ! epaisseur de couche
192	DO i = 1, klon
193	delp(i, k) = paprs(i, k) - paprs(i, k + 1)
194	END DO
195	END DO
196
197	! Initialization:
198	rugmer = 0.
199	cdragh = 0.
200	cdragm = 0.
201	dflux_t = 0.
202	dflux_q = 0.
203	ypct = 0.
204	yrugos = 0.
205	ypaprs = 0.
206	ypplay = 0.
207	ydelp = 0.
208	yrugoro = 0.
209	d_ts = 0.
210	flux_t = 0.
211	flux_q = 0.
212	flux_u = 0.
213	flux_v = 0.
214	fluxlat = 0.
215	d_t = 0.
216	d_q = 0.
217	d_u = 0.
218	d_v = 0.
219	coefh = 0.
220	fqcalving = 0.
221	run_off_lic = 0.
222
223	! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
224	! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
225	! (\`a affiner).
226
227	pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
228	pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
229	pctsrf_pot(:, is_oce) = 1. - zmasq
230	pctsrf_pot(:, is_sic) = 1. - zmasq
231
232	! Tester si c'est le moment de lire le fichier:
233	if (mod(itap - 1, lmt_pas) == 0) then
234	CALL interfoce_lim(julien, pctsrf_new_oce, pctsrf_new_sic)
235	endif
236
237	! Boucler sur toutes les sous-fractions du sol:
238
239	loop_surface: DO nsrf = 1, nbsrf
240	! Define ni and knon:
241
242	ni = 0
243	knon = 0
244
245	DO i = 1, klon
246	! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
247	! "potentielles"
248	IF (pctsrf_pot(i, nsrf) > epsfra) THEN
249	knon = knon + 1
250	ni(knon) = i
251	END IF
252	END DO
253
254	if_knon: IF (knon /= 0) then
255	DO j = 1, knon
256	i = ni(j)
257	ypct(j) = pctsrf(i, nsrf)
258	yts(j) = ftsol(i, nsrf)
259	snow(j) = fsnow(i, nsrf)
260	yqsurf(j) = qsurf(i, nsrf)
261	yalb(j) = falbe(i, nsrf)
262	yrain_fall(j) = rain_fall(i)
263	ysnow_fall(j) = snow_fall(i)
264	yagesno(j) = agesno(i, nsrf)
265	yrugos(j) = frugs(i, nsrf)
266	yrugoro(j) = rugoro(i)
267	yrads(j) = fsolsw(i, nsrf) + fsollw(i, nsrf)
268	ypaprs(j, klev + 1) = paprs(i, klev + 1)
269	y_run_off_lic_0(j) = run_off_lic_0(i)
270	END DO
271
272	! For continent, copy soil water content
273	IF (nsrf == is_ter) yqsol(:knon) = qsol(ni(:knon))
274
275	ytsoil(:knon, :) = ftsoil(ni(:knon), :, nsrf)
276
277	DO k = 1, klev
278	DO j = 1, knon
279	i = ni(j)
280	ypaprs(j, k) = paprs(i, k)
281	ypplay(j, k) = pplay(i, k)
282	ydelp(j, k) = delp(i, k)
283	yu(j, k) = u(i, k)
284	yv(j, k) = v(i, k)
285	yt(j, k) = t(i, k)
286	yq(j, k) = q(i, k)
287	END DO
288	END DO
289
290	! Calculer les géopotentiels de chaque couche:
291
292	zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
293	+ ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
294
295	DO k = 2, klev
296	zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
297	* (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
298	* (ypplay(:knon, k - 1) - ypplay(:knon, k))
299	ENDDO
300
301	CALL cdrag(nsrf, sqrt(yu(:knon, 1)2 + yv(:knon, 1)2), &
302	yt(:knon, 1), yq(:knon, 1), zgeop(:knon, 1), ypaprs(:knon, 1), &
303	yts(:knon), yqsurf(:knon), yrugos(:knon), ycdragm(:knon), &
304	ycdragh(:knon))
305
306	IF (iflag_pbl == 1) THEN
307	ycdragm(:knon) = max(ycdragm(:knon), 0.)
308	ycdragh(:knon) = max(ycdragh(:knon), 0.)
309	end IF
310
311	! on met un seuil pour ycdragm et ycdragh
312	IF (nsrf == is_oce) THEN
313	ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
314	ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
315	END IF
316
317	IF (iflag_pbl >= 6) yq2(:knon, :) = q2(ni(:knon), :, nsrf)
318	call coef_diff_turb(nsrf, ni(:knon), ypaprs(:knon, :), &
319	ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
320	yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
321	ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))
322
323	CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
324	ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
325	ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
326	y_flux_u(:knon))
327	CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
328	ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
329	ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
330	y_flux_v(:knon))
331
332	CALL clqh(julien, nsrf, ni(:knon), ytsoil(:knon, :), yqsol(:knon), &
333	mu0(ni(:knon)), yrugos(:knon), yrugoro(:knon), yu(:knon, 1), &
334	yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), yt(:knon, :), &
335	yq(:knon, :), yts(:knon), ypaprs(:knon, :), ypplay(:knon, :), &
336	ydelp(:knon, :), yrads(:knon), yalb(:knon), snow(:knon), &
337	yqsurf(:knon), yrain_fall(:knon), ysnow_fall(:knon), &
338	yfluxlat(:knon), pctsrf_new_sic(ni(:knon)), yagesno(:knon), &
339	y_d_t(:knon, :), y_d_q(:knon, :), y_d_ts(:knon), &
340	yz0_new(:knon), y_flux_t(:knon), y_flux_q(:knon), &
341	y_dflux_t(:knon), y_dflux_q(:knon), y_fqcalving(:knon), &
342	y_ffonte(:knon), y_run_off_lic_0(:knon), y_run_off_lic(:knon))
343
344	! calculer la longueur de rugosite sur ocean
345
346	yrugm = 0.
347
348	IF (nsrf == is_oce) THEN
349	DO j = 1, knon
350	yrugm(j) = 0.018 * ycdragm(j) * (yu(j, 1)2 + yv(j, 1)2) &
351	/ rg + 0.11 * 14E-6 &
352	/ sqrt(ycdragm(j) * (yu(j, 1)2 + yv(j, 1)2))
353	yrugm(j) = max(1.5E-05, yrugm(j))
354	END DO
355	END IF
356
357	DO k = 1, klev
358	DO j = 1, knon
359	i = ni(j)
360	y_d_t(j, k) = y_d_t(j, k) * ypct(j)
361	y_d_q(j, k) = y_d_q(j, k) * ypct(j)
362	y_d_u(j, k) = y_d_u(j, k) * ypct(j)
363	y_d_v(j, k) = y_d_v(j, k) * ypct(j)
364	END DO
365	END DO
366
367	flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
368	flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
369	flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
370	flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
371
372	falbe(:, nsrf) = 0.
373	fsnow(:, nsrf) = 0.
374	qsurf(:, nsrf) = 0.
375	frugs(:, nsrf) = 0.
376	DO j = 1, knon
377	i = ni(j)
378	d_ts(i, nsrf) = y_d_ts(j)
379	falbe(i, nsrf) = yalb(j)
380	fsnow(i, nsrf) = snow(j)
381	qsurf(i, nsrf) = yqsurf(j)
382	frugs(i, nsrf) = yz0_new(j)
383	fluxlat(i, nsrf) = yfluxlat(j)
384	IF (nsrf == is_oce) THEN
385	rugmer(i) = yrugm(j)
386	frugs(i, nsrf) = yrugm(j)
387	END IF
388	agesno(i, nsrf) = yagesno(j)
389	fqcalving(i, nsrf) = y_fqcalving(j)
390	ffonte(i, nsrf) = y_ffonte(j)
391	cdragh(i) = cdragh(i) + ycdragh(j) * ypct(j)
392	cdragm(i) = cdragm(i) + ycdragm(j) * ypct(j)
393	dflux_t(i) = dflux_t(i) + y_dflux_t(j) * ypct(j)
394	dflux_q(i) = dflux_q(i) + y_dflux_q(j) * ypct(j)
395	END DO
396	IF (nsrf == is_ter) THEN
397	qsol(ni(:knon)) = yqsol(:knon)
398	else IF (nsrf == is_lic) THEN
399	DO j = 1, knon
400	i = ni(j)
401	run_off_lic_0(i) = y_run_off_lic_0(j)
402	run_off_lic(i) = y_run_off_lic(j)
403	END DO
404	END IF
405
406	ftsoil(:, :, nsrf) = 0.
407	ftsoil(ni(:knon), :, nsrf) = ytsoil(:knon, :)
408
409	DO j = 1, knon
410	i = ni(j)
411	DO k = 1, klev
412	d_t(i, k) = d_t(i, k) + y_d_t(j, k)
413	d_q(i, k) = d_q(i, k) + y_d_q(j, k)
414	d_u(i, k) = d_u(i, k) + y_d_u(j, k)
415	d_v(i, k) = d_v(i, k) + y_d_v(j, k)
416	END DO
417	END DO
418
419	forall (k = 2:klev) coefh(ni(:knon), k) &
420	= coefh(ni(:knon), k) + ycoefh(:knon, k) * ypct(:knon)
421
422	! diagnostic t, q a 2m et u, v a 10m
423
424	DO j = 1, knon
425	i = ni(j)
426	u1(j) = yu(j, 1) + y_d_u(j, 1)
427	v1(j) = yv(j, 1) + y_d_v(j, 1)
428	tair1(j) = yt(j, 1) + y_d_t(j, 1)
429	qair1(j) = yq(j, 1) + y_d_q(j, 1)
430	zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
431	1))) * (ypaprs(j, 1)-ypplay(j, 1))
432	tairsol(j) = yts(j) + y_d_ts(j)
433	rugo1(j) = yrugos(j)
434	IF (nsrf == is_oce) THEN
435	rugo1(j) = frugs(i, nsrf)
436	END IF
437	psfce(j) = ypaprs(j, 1)
438	patm(j) = ypplay(j, 1)
439	END DO
440
441	CALL stdlevvar(nsrf, u1(:knon), v1(:knon), tair1(:knon), qair1, &
442	zgeo1, tairsol, yqsurf(:knon), rugo1, psfce, patm, yt2m, yq2m, &
443	yt10m, yq10m, wind10m(:knon), ustar(:knon))
444
445	DO j = 1, knon
446	i = ni(j)
447	t2m(i, nsrf) = yt2m(j)
448	q2m(i, nsrf) = yq2m(j)
449
450	u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
451	/ sqrt(u1(j)2 + v1(j)2)
452	v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
453	/ sqrt(u1(j)2 + v1(j)2)
454	END DO
455
456	CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
457	y_flux_q(:knon), yu(:knon, :), yv(:knon, :), yt(:knon, :), &
458	yq(:knon, :), ypblh(:knon), ycapcl, yoliqcl, ycteicl, ypblt, &
459	ytherm, ylcl)
460
461	DO j = 1, knon
462	i = ni(j)
463	pblh(i, nsrf) = ypblh(j)
464	plcl(i, nsrf) = ylcl(j)
465	capcl(i, nsrf) = ycapcl(j)
466	oliqcl(i, nsrf) = yoliqcl(j)
467	cteicl(i, nsrf) = ycteicl(j)
468	pblt(i, nsrf) = ypblt(j)
469	therm(i, nsrf) = ytherm(j)
470	END DO
471
472	IF (iflag_pbl >= 6) q2(ni(:knon), :, nsrf) = yq2(:knon, :)
473	else
474	fsnow(:, nsrf) = 0.
475	end IF if_knon
476	END DO loop_surface
477
478	! On utilise les nouvelles surfaces
479	frugs(:, is_oce) = rugmer
480	pctsrf(:, is_oce) = pctsrf_new_oce
481	pctsrf(:, is_sic) = pctsrf_new_sic
482
483	CALL histwrite_phy("run_off_lic", run_off_lic)
484
485	END SUBROUTINE pbl_surface
486
487	end module pbl_surface_m