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, play, fsnow, fqsurf, falbe, fluxlat, &
       rain_fall, snow_fall, 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, albsol, sollw, solsw, &
       tsol)

    ! 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: masque
    use stdlevvar_m, only: stdlevvar
    USE suphec_m, ONLY: rd, rg, rsigma
    use time_phylmdz, only: itap

    REAL, INTENT(inout):: pctsrf(klon, nbsrf)
    ! 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)
    ! skin temperature of surface fraction, in 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):: play(klon, klev) ! pression au milieu de couche (Pa)
    REAL, INTENT(inout):: fsnow(:, :) ! (klon, nbsrf) \'epaisseur neigeuse
    REAL, INTENT(inout):: fqsurf(klon, nbsrf)
    REAL, intent(inout):: falbe(klon, nbsrf)

    REAL, intent(out):: fluxlat(:, :) ! (klon, nbsrf)
    ! flux de chaleur latente, en W m-2

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

    REAL, intent(out):: albsol(:) ! (klon)
    ! albedo du sol total, visible, moyen par maille

    REAL, intent(in):: sollw(:) ! (klon)
    ! surface net downward longwave flux, in W m-2

    REAL, intent(in):: solsw(:) ! (klon)
    ! surface net downward shortwave flux, in W m-2

    REAL, intent(in):: tsol(:) ! (klon)

    ! Local:

    REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous-surface
    REAL fsolsw(klon, nbsrf) ! flux solaire absorb\'e pour chaque sous-surface

    ! 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), ypctsrf(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), radsol(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)

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

    albsol = sum(falbe * pctsrf, dim = 2)

    ! R\'epartition sous maille des flux longwave et shortwave
    ! R\'epartition du longwave par sous-surface lin\'earis\'ee

    forall (nsrf = 1:nbsrf)
       fsollw(:, nsrf) = sollw + 4. * RSIGMA * tsol**3 &
            * (tsol - ftsol(:, nsrf))
       fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
    END forall

    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.
    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. - masque
    pctsrf_pot(:, is_sic) = 1. - masque

    ! 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
          ypctsrf(:knon) = pctsrf(ni(:knon), nsrf)
          yts(:knon) = ftsol(ni(:knon), nsrf)
          snow(:knon) = fsnow(ni(:knon), nsrf)
          yqsurf(:knon) = fqsurf(ni(:knon), nsrf)
          yalb(:knon) = falbe(ni(:knon), nsrf)
          yrain_fall(:knon) = rain_fall(ni(:knon))
          ysnow_fall(:knon) = snow_fall(ni(:knon))
          yagesno(:knon) = agesno(ni(:knon), nsrf)
          yrugos(:knon) = frugs(ni(:knon), nsrf)
          yrugoro(:knon) = rugoro(ni(:knon))
          radsol(:knon) = fsolsw(ni(:knon), nsrf) + fsollw(ni(:knon), nsrf)
          ypaprs(:knon, klev + 1) = paprs(ni(:knon), klev + 1)
          y_run_off_lic_0(:knon) = run_off_lic_0(ni(:knon))

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