phylmd/Interface_surf/pbl_surface.f

module pbl_surface_m

  IMPLICIT NONE

contains

  SUBROUTINE pbl_surface(dtime, 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: 1993 Aug. 18th
    ! 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 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(IN):: dtime ! interval du temps (secondes)

    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(klon, klev), 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 (Cp 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)

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

    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)
    ! ffonte----Flux thermique utilise pour fondre la neige
    REAL run_off_lic_0(klon)

    ! Local:

    LOGICAL:: firstcal = .true.

    ! 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)
    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.
    yu = 0.
    yv = 0.
    yt = 0.
    yq = 0.
    y_dflux_t = 0.
    y_dflux_q = 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.

    ! 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
       ! Chercher les indices :
       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) then
             DO k = 1, klev + 1
                DO j = 1, knon
                   i = ni(j)
                   yq2(j, k) = q2(i, k, nsrf)
                END DO
             END DO
          end IF

          call coef_diff_turb(dtime, 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(dtime, 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(dtime, 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))

          ! calculer la diffusion de "q" et de "h"
          CALL clqh(dtime, julien, firstcal, nsrf, ni(:knon), &
               ytsoil(:knon, :), yqsol(:knon), mu0, 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, yrads(:knon), &
               yalb(:knon), snow(:knon), yqsurf, yrain_f, ysnow_f, &
               yfluxlat(:knon), pctsrf_new_sic, 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, y_run_off_lic_0)

          ! 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 j = 1, knon
             y_dflux_t(j) = y_dflux_t(j) * ypct(j)
             y_dflux_q(j) = y_dflux_q(j) * ypct(j)
          END DO

          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)
             dflux_q(i) = dflux_q(i) + y_dflux_q(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)
             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, yv, yt, yq, 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

          DO j = 1, knon
             DO k = 1, klev + 1
                i = ni(j)
                q2(i, k, nsrf) = yq2(j, k)
             END DO
          END DO
       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

    firstcal = .false.

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