phylmd/Interface_surf/pbl_surface.f90

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)
    ! column-density of mass of snow at the surface, in kg m-2

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