trunk/phylmd/clmain.f

module clmain_m

  IMPLICIT NONE

contains

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

    ! From phylmd/clmain.F, version 1.6, 2005/11/16 14:47:19
    ! Author: Z. X. Li (LMD/CNRS), date: 1993/08/18
    ! 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 clcdrag_m, only: clcdrag
    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, zmasq
    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 stdlevvar_m, only: stdlevvar
    USE suphec_m, ONLY: rd, rg
    use time_phylmdz, only: itap

    REAL, INTENT(IN):: dtime ! interval du temps (secondes)

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

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

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

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

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

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

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

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

    REAL d_t(klon, klev), d_q(klon, klev)
    ! d_t------output-R- le changement pour "t"
    ! d_q------output-R- le changement pour "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 fqcalving(klon, nbsrf), ffonte(klon, nbsrf)
    ! ffonte----Flux thermique utilise pour fondre la neige
    ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
    !           hauteur de neige, en kg / m2 / s
    REAL run_off_lic_0(klon)

    ! Local:

    LOGICAL:: firstcal = .true.

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

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

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

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

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

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

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

    ytherm = 0.

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

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

    ! 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 clcdrag(nsrf, yu(:knon, 1), yv(:knon, 1), yt(:knon, 1), &
               yq(:knon, 1), zgeop(:knon, 1), yts(:knon), yqsurf(:knon), &
               yrugos(:knon), ycdragm(:knon), ycdragh(:knon)) 

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

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

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

          call coef_diff_turb(dtime, nsrf, ni(:knon), ypaprs(:knon, :), &
               ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
               yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
               ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))

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

          ! calculer la diffusion de "q" et de "h"
          CALL clqh(dtime, julien, firstcal, nsrf, ni(:knon), &
               ytsoil(:knon, :), yqsol(:knon), mu0, yrugos, yrugoro, &
               yu(:knon, 1), yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), &
               yt, yq, yts(:knon), ypaprs, ypplay, ydelp, yrads(:knon), &
               yalb(:knon), snow(:knon), yqsurf, yrain_f, ysnow_f, &
               yfluxlat(:knon), pctsrf_new_sic, yagesno(:knon), y_d_t, y_d_q, &
               y_d_ts(:knon), yz0_new, y_flux_t(:knon), y_flux_q(:knon), &
               y_dflux_t(:knon), y_dflux_q(:knon), y_fqcalving, y_ffonte, &
               y_run_off_lic_0)

          ! calculer la longueur de rugosite sur ocean
          yrugm = 0.
          IF (nsrf == is_oce) THEN
             DO j = 1, knon
                yrugm(j) = 0.018 * ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2) &
                     / rg + 0.11 * 14E-6 &
                     / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2))
                yrugm(j) = max(1.5E-05, yrugm(j))
             END DO
          END IF
          DO j = 1, knon
             y_dflux_t(j) = y_dflux_t(j) * ypct(j)
             y_dflux_q(j) = y_dflux_q(j) * ypct(j)
          END DO

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

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

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

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

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

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

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

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

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

             qairsol(j) = yqsurf(j)
          END DO

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

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

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

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

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

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

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

    firstcal = .false.

  END SUBROUTINE clmain

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