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

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