phylmd/Interface_surf/pbl_surface.f

module pbl_surface_m

  IMPLICIT NONE

contains

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

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

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

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

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

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

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

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

    REAL, INTENT(IN):: paprs(klon, klev + 1) ! pression a intercouche (Pa)
    REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
    REAL, INTENT(inout):: fsnow(:, :) ! (klon, nbsrf) \'epaisseur neigeuse
    REAL qsurf(klon, nbsrf)
    REAL evap(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_f(klon)
    ! solid water mass flux (kg / m2 / s), positive down

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    ! Local:

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

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

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

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

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

    REAL qairsol(klon), zgeo1(klon)
    REAL rugo1(klon)
    REAL zgeop(klon, klev)

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

    ytherm = 0.

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

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

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

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

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

    ! Boucler sur toutes les sous-fractions du sol:

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

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

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

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

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

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

          ! Calculer les géopotentiels de chaque couche:

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

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

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

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

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

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

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

          ! calculer la longueur de rugosite sur ocean

          yrugm = 0.

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

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

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

          evap(:, nsrf) = -flux_q(:, nsrf)

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

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

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

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

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

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

             qairsol(j) = yqsurf(j)
          END DO

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