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

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