phylmd/Interface_surf/pbl_surface.f

module pbl_surface_m

  IMPLICIT NONE

contains

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    ! Local:

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

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

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

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

    REAL yt10m(klon), yq10m(klon)
    REAL ypblh(klon)
    REAL ylcl(klon)
    REAL ycapcl(klon)
    REAL yoliqcl(klon)
    REAL ycteicl(klon)
    REAL ypblt(klon)
    REAL ytherm(klon)
    REAL u1(klon), v1(klon)
    REAL tair1(klon), qair1(klon), tairsol(klon)
    REAL psfce(klon), patm(klon)
    REAL 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.
    yrain_f = 0.
    ysnow_f = 0.
    yrugos = 0.
    ypaprs = 0.
    ypplay = 0.
    ydelp = 0.
    yrugoro = 0.
    d_ts = 0.
    flux_t = 0.
    flux_q = 0.
    flux_u = 0.
    flux_v = 0.
    fluxlat = 0.
    d_t = 0.
    d_q = 0.
    d_u = 0.
    d_v = 0.
    coefh = 0.
    fqcalving = 0.
    run_off_lic = 0.

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

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

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

    ! Boucler sur toutes les sous-fractions du sol:

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

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

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

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

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

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

          ! Calculer les géopotentiels de chaque couche:

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

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

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

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

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

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

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

          ! calculer la longueur de rugosite sur ocean

          yrugm = 0.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    CALL histwrite_phy("run_off_lic", run_off_lic)

  END SUBROUTINE pbl_surface

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