phylmd/Interface_surf/pbl_surface.f

module pbl_surface_m

  IMPLICIT NONE

contains

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    ! Local:

    LOGICAL:: firstcal = .true.

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

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

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

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

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

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

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

    ytherm = 0.

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

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

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

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

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

    ! Boucler sur toutes les sous-fractions du sol:

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

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

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

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

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

          ! Calculer les géopotentiels de chaque couche:

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

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

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

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

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

          IF (iflag_pbl >= 6) then
             DO k = 1, klev + 1
                DO j = 1, knon
                   i = ni(j)
                   yq2(j, k) = q2(i, k, nsrf)
                END DO
             END DO
          end IF

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

          ! calculer la longueur de rugosite sur ocean

          yrugm = 0.

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

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

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

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

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

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

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

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

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

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

             qairsol(j) = yqsurf(j)
          END DO

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

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

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

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

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

          DO j = 1, knon
             DO k = 1, klev + 1
                i = ni(j)
                q2(i, k, nsrf) = yq2(j, k)
             END DO
          END DO
       else
          fsnow(:, nsrf) = 0.
       end IF if_knon
    END DO loop_surface

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

    firstcal = .false.

  END SUBROUTINE pbl_surface

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