Sources/phylmd/clmain.f

module clmain_m

  IMPLICIT NONE

contains

  SUBROUTINE clmain(dtime, pctsrf, t, q, u, v, julien, mu0, ftsol, cdmmax, &
       cdhmax, ksta, ksta_ter, ok_kzmin, 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, ycoefh, t2m, q2m, &
       u10m_srf, v10m_srf, pblh, capcl, oliqcl, cteicl, pblt, therm, trmb1, &
       trmb2, trmb3, 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 clqh_m, only: clqh
    use clvent_m, only: clvent
    use coefkz_m, only: coefkz
    use coefkzmin_m, only: coefkzmin
    use coefkz2_m, only: coefkz2
    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, rkappa
    use time_phylmdz, only: itap
    use ustarhb_m, only: ustarhb
    use yamada4_m, only: yamada4

    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(IN):: ksta, ksta_ter
    LOGICAL, INTENT(IN):: ok_kzmin

    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):: ycoefh(klon, klev)
    ! Pour pouvoir extraire les coefficients d'\'echange, le champ
    ! "ycoefh" 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 trmb1(klon, nbsrf)
    ! trmb1-------deep_cape
    REAL trmb2(klon, nbsrf)
    ! trmb2--------inhibition
    REAL trmb3(klon, nbsrf)
    ! trmb3-------Point Omega
    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), yrugos(klon), ypct(klon), yz0_new(klon)
    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 coefh(klon, klev), coefm(klon, klev)
    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 ycoefm0(klon, klev), ycoefh0(klon, klev)
    REAL yzlay(klon, klev), zlev(klon, klev + 1), yteta(klon, klev)
    REAL ykmm(klon, klev + 1), ykmn(klon, klev + 1)
    REAL ykmq(klon, klev + 1)
    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 ytrmb1(klon)
    REAL ytrmb2(klon)
    REAL ytrmb3(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)

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

    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.
    ycoefh = 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 Cdrag et les coefficients d'echange
          CALL coefkz(nsrf, ypaprs, ypplay, ksta, ksta_ter, yts(:knon), &
               yrugos, yu, yv, yt, yq, yqsurf(:knon), coefm(:knon, 2:), &
               coefh(:knon, 2:), coefm(:knon, 1), coefh(:knon, 1))

          IF (iflag_pbl == 1) THEN
             CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0(:knon, 2:), &
                  ycoefh0(:knon, 2:))
             ycoefm0(:knon, 1) = 0.
             ycoefh0(:knon, 1) = 0.
             coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
             coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
          END IF

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

          IF (ok_kzmin) THEN
             ! Calcul d'une diffusion minimale pour les conditions tres stables
             CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
                  coefm(:knon, 1), ycoefm0(:knon, 2:), ycoefh0(:knon, 2:))
             coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
             coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
          END IF

          IF (iflag_pbl >= 6) THEN
             ! Mellor et Yamada adapt\'e \`a Mars, Richard Fournier et
             ! Fr\'ed\'eric Hourdin
             yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
                  + ypplay(:knon, 1))) &
                  * (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg

             DO k = 2, klev
                yzlay(:knon, k) = yzlay(:knon, k-1) &
                     + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
                     / ypaprs(1:knon, k) &
                     * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
             END DO

             DO k = 1, klev
                yteta(1:knon, k) = yt(1:knon, k) * (ypaprs(1:knon, 1) &
                     / ypplay(1:knon, k))**rkappa * (1. + 0.61 * yq(1:knon, k))
             END DO

             zlev(:knon, 1) = 0.
             zlev(:knon, klev + 1) = 2. * yzlay(:knon, klev) &
                  - yzlay(:knon, klev - 1)

             DO k = 2, klev
                zlev(:knon, k) = 0.5 * (yzlay(:knon, k) + yzlay(:knon, k-1))
             END DO

             DO k = 1, klev + 1
                DO j = 1, knon
                   i = ni(j)
                   yq2(j, k) = q2(i, k, nsrf)
                END DO
             END DO

             ustar(:knon) = ustarhb(yu(:knon, 1), yv(:knon, 1), coefm(:knon, 1))
             CALL yamada4(dtime, rg, zlev(:knon, :), yzlay(:knon, :), &
                  yu(:knon, :), yv(:knon, :), yteta(:knon, :), &
                  coefm(:knon, 1), yq2(:knon, :), ykmm(:knon, :), &
                  ykmn(:knon, :), ykmq(:knon, :), ustar(:knon))
             coefm(:knon, 2:) = ykmm(:knon, 2:klev)
             coefh(:knon, 2:) = ykmn(:knon, 2:klev)
          END IF

          CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), coefm(:knon, 2:), &
               coefm(:knon, 1), 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), coefm(:knon, 2:), &
               coefm(:knon, 1), 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), coefh(:knon, 2:), coefh(:knon, 1), &
               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 * coefm(j, 1) * (yu(j, 1)**2 + yv(j, 1)**2) &
                     / rg + 0.11 * 14E-6 &
                     / sqrt(coefm(j, 1) * (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)
                coefh(j, k) = coefh(j, k) * ypct(j)
                coefm(j, k) = coefm(j, k) * ypct(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) + coefh(j, 1)
             cdragm(i) = cdragm(i) + coefm(j, 1)
             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)
                ycoefh(i, k) = ycoefh(i, k) + coefh(j, k)
             END DO
          END DO

          ! 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)

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