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, &
       flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2, dflux_t, dflux_q, &
       coefh, t2m, q2m, u10m_srf, v10m_srf, pblh, capcl, oliqcl, cteicl, pblt, &
       therm, plcl, fqcalving, ffonte, run_off_lic_0, albsol, sollw, solsw, &
       tsol)

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

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

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

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

    REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
    REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg / kg)
    REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
    INTEGER, INTENT(IN):: julien ! jour de l'annee en cours
    REAL, intent(in):: mu0(klon) ! cosinus de l'angle solaire zenithal     

    REAL, INTENT(INout):: ftsol(:, :) ! (klon, nbsrf)
    ! skin temperature of surface fraction, in K

    REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh

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

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

    REAL, INTENT(IN):: paprs(klon, klev + 1) ! pression a intercouche (Pa)
    REAL, INTENT(IN):: play(klon, klev) ! pression au milieu de couche (Pa)
    REAL, INTENT(inout):: fsnow(:, :) ! (klon, nbsrf) \'epaisseur neigeuse
    REAL, INTENT(inout):: fqsurf(klon, nbsrf)
    REAL, intent(inout):: falbe(klon, nbsrf)

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

    REAL, intent(in):: rain_fall(klon)
    ! liquid water mass flux (kg / m2 / s), positive down

    REAL, intent(in):: snow_fall(klon)
    ! solid water mass flux (kg / m2 / s), positive down

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

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

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

    REAL, intent(out):: 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 d_ts(klon, nbsrf) ! variation of ftsol
    REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous-surface
    REAL fsolsw(klon, nbsrf) ! flux solaire absorb\'e pour chaque sous-surface

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

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

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

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

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

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

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

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

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

    ytherm = 0.

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

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

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

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

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

    ! Boucler sur toutes les sous-fractions du sol:

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

       ni = 0
       knon = 0

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

       if_knon: IF (knon /= 0) then
          ypctsrf(:knon) = pctsrf(ni(:knon), nsrf)
          yts(:knon) = ftsol(ni(:knon), nsrf)
          snow(:knon) = fsnow(ni(:knon), nsrf)
          yqsurf(:knon) = fqsurf(ni(:knon), nsrf)
          yalb(:knon) = falbe(ni(:knon), nsrf)
          yrain_fall(:knon) = rain_fall(ni(:knon))
          ysnow_fall(:knon) = snow_fall(ni(:knon))
          yagesno(:knon) = agesno(ni(:knon), nsrf)
          yrugos(:knon) = frugs(ni(:knon), nsrf)
          yrugoro(:knon) = rugoro(ni(:knon))
          radsol(:knon) = fsolsw(ni(:knon), nsrf) + fsollw(ni(:knon), nsrf)
          ypaprs(:knon, klev + 1) = paprs(ni(:knon), klev + 1)
          y_run_off_lic_0(:knon) = run_off_lic_0(ni(:knon))

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

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

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

          ! Calculer les géopotentiels de chaque couche:

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

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

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

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

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

          IF (iflag_pbl >= 6) yq2(:knon, :) = q2(ni(:knon), :, nsrf)
          call coef_diff_turb(nsrf, ni(:knon), ypaprs(:knon, :), &
               ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
               yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
               ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))

          CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
               ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
               ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
               y_flux_u(:knon))
          CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
               ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
               ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
               y_flux_v(:knon))

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

          ! calculer la longueur de rugosite sur ocean

          yrugm = 0.

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

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

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

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

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

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

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

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

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

          CALL stdlevvar(nsrf, u1(:knon), v1(:knon), tair1(:knon), qair1, &
               zgeo1, tairsol, yqsurf(:knon), rugo1, psfce, patm, yt2m, yq2m, &
               yt10m, yq10m, wind10m(:knon), ustar(:knon))

          DO j = 1, knon
             i = ni(j)
             t2m(i, nsrf) = yt2m(j)
             q2m(i, nsrf) = yq2m(j)

             u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
                  / sqrt(u1(j)**2 + v1(j)**2)
             v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
                  / sqrt(u1(j)**2 + v1(j)**2)
          END DO

          CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
               y_flux_q(:knon), yu(:knon, :), yv(:knon, :), yt(:knon, :), &
               yq(:knon, :), ypblh(:knon), ycapcl, yoliqcl, ycteicl, ypblt, &
               ytherm, ylcl)

          DO j = 1, knon
             i = ni(j)
             pblh(i, nsrf) = ypblh(j)
             plcl(i, nsrf) = ylcl(j)
             capcl(i, nsrf) = ycapcl(j)
             oliqcl(i, nsrf) = yoliqcl(j)
             cteicl(i, nsrf) = ycteicl(j)
             pblt(i, nsrf) = ypblt(j)
             therm(i, nsrf) = ytherm(j)
          END DO

          IF (iflag_pbl >= 6) q2(ni(:knon), :, nsrf) = yq2(:knon, :)
       else
          fsnow(:, nsrf) = 0.
       end IF if_knon
    END DO loop_surface

    ! On utilise les nouvelles surfaces
    frugs(:, is_oce) = rugmer
    pctsrf(:, is_oce) = pctsrf_new_oce
    pctsrf(:, is_sic) = pctsrf_new_sic

    CALL histwrite_phy("run_off_lic", run_off_lic)
    ftsol = ftsol + d_ts ! update surface temperature
    CALL histwrite_phy("dtsvdfo", d_ts(:, is_oce))
    CALL histwrite_phy("dtsvdft", d_ts(:, is_ter))
    CALL histwrite_phy("dtsvdfg", d_ts(:, is_lic))
    CALL histwrite_phy("dtsvdfi", d_ts(:, is_sic))

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