trunk/phylmd/pbl_surface.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, 2:klev), coefm(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 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, :), &
               coefh(:knon, :), ycdragm(:knon), ycdragh(:knon))

          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, 2:))
             coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, 2:))
             ycdragm(:knon) = max(ycdragm(:knon), ycoefm0(:knon, 1))
             ycdragh(:knon) = max(ycdragh(:knon), ycoefh0(:knon, 1))
          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 (ok_kzmin) THEN
             ! Calcul d'une diffusion minimale pour les conditions tres stables
             CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
                  ycdragm(:knon), ycoefm0(:knon, 2:), ycoefh0(:knon, 2:))
             coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, 2:))
             coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, 2:))
             ycdragm(:knon) = max(ycdragm(:knon), ycoefm0(:knon, 1))
             ycdragh(:knon) = max(ycdragh(:knon), ycoefh0(:knon, 1))
          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), ycdragm(:knon))
             CALL yamada4(dtime, rg, zlev(:knon, :), yzlay(:knon, :), &
                  yu(:knon, :), yv(:knon, :), yteta(:knon, :), &
                  ycdragm(:knon), 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, :), &
               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), coefm(: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), coefh(: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 = 2, klev
             DO j = 1, knon
                i = ni(j)
                coefh(j, k) = coefh(j, k) * ypct(j)
                coefm(j, k) = coefm(j, k) * ypct(j)
             END DO
          END DO
          DO j = 1, knon
             i = ni(j)
             ycdragh(j) = ycdragh(j) * ypct(j)
             ycdragm(j) = ycdragm(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)
             cdragm(i) = cdragm(i) + ycdragm(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
          
          DO j = 1, knon
             i = ni(j)
             DO k = 2, klev
                ycoefh(i, k) = ycoefh(i, k) + coefh(j, k)
             END DO
          END DO

          DO j = 1, knon
             i = ni(j)
             ycoefh(i, 1) = ycoefh(i, 1) + ycdragh(j)
          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	guez	237	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	237	REAL coefh(klon, 2:klev), coefm(klon, 2:klev)
168			real ycdragh(klon), ycdragm(klon)
169	guez	38	REAL yu(klon, klev), yv(klon, klev)
170			REAL yt(klon, klev), yq(klon, klev)
171	guez	225	REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)
172	guez	38	REAL ycoefm0(klon, klev), ycoefh0(klon, klev)
173	guez	227	REAL yzlay(klon, klev), zlev(klon, klev + 1), yteta(klon, klev)
174	guez	225	REAL ykmm(klon, klev + 1), ykmn(klon, klev + 1)
175			REAL ykmq(klon, klev + 1)
176			REAL yq2(klon, klev + 1)
177	guez	38	REAL delp(klon, klev)
178			INTEGER i, k, nsrf
179			INTEGER ni(klon), knon, j
180	guez	40
181	guez	38	REAL pctsrf_pot(klon, nbsrf)
182	guez	145	! "pourcentage potentiel" pour tenir compte des \'eventuelles
183	guez	40	! apparitions ou disparitions de la glace de mer
184	guez	15
185	guez	227	REAL yt2m(klon), yq2m(klon), wind10m(klon)
186			REAL ustar(klon)
187	guez	15
188	guez	38	REAL yt10m(klon), yq10m(klon)
189			REAL ypblh(klon)
190			REAL ylcl(klon)
191			REAL ycapcl(klon)
192			REAL yoliqcl(klon)
193			REAL ycteicl(klon)
194			REAL ypblt(klon)
195			REAL ytherm(klon)
196			REAL ytrmb1(klon)
197			REAL ytrmb2(klon)
198			REAL ytrmb3(klon)
199	guez	227	REAL u1(klon), v1(klon)
200	guez	38	REAL tair1(klon), qair1(klon), tairsol(klon)
201			REAL psfce(klon), patm(klon)
202	guez	15
203	guez	38	REAL qairsol(klon), zgeo1(klon)
204			REAL rugo1(klon)
205	guez	15
206	guez	38	!------------------------------------------------------------
207	guez	15
208	guez	38	ytherm = 0.
209	guez	15
210	guez	38	DO k = 1, klev ! epaisseur de couche
211			DO i = 1, klon
212	guez	225	delp(i, k) = paprs(i, k) - paprs(i, k + 1)
213	guez	38	END DO
214			END DO
215	guez	15
216	guez	40	! Initialization:
217			rugmer = 0.
218			cdragh = 0.
219			cdragm = 0.
220			dflux_t = 0.
221			dflux_q = 0.
222			ypct = 0.
223			yqsurf = 0.
224			yrain_f = 0.
225			ysnow_f = 0.
226			yrugos = 0.
227			ypaprs = 0.
228			ypplay = 0.
229			ydelp = 0.
230			yu = 0.
231			yv = 0.
232			yt = 0.
233			yq = 0.
234			y_dflux_t = 0.
235			y_dflux_q = 0.
236	guez	38	yrugoro = 0.
237	guez	40	d_ts = 0.
238	guez	38	flux_t = 0.
239			flux_q = 0.
240			flux_u = 0.
241			flux_v = 0.
242	guez	214	fluxlat = 0.
243	guez	40	d_t = 0.
244			d_q = 0.
245			d_u = 0.
246			d_v = 0.
247	guez	70	ycoefh = 0.
248	guez	15
249	guez	145	! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
250			! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
251			! (\`a affiner)
252	guez	15
253	guez	202	pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
254			pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
255	guez	38	pctsrf_pot(:, is_oce) = 1. - zmasq
256			pctsrf_pot(:, is_sic) = 1. - zmasq
257	guez	15
258	guez	202	! Tester si c'est le moment de lire le fichier:
259			if (mod(itap - 1, lmt_pas) == 0) then
260	guez	221	CALL interfoce_lim(julien, pctsrf_new_oce, pctsrf_new_sic)
261	guez	202	endif
262
263	guez	99	! Boucler sur toutes les sous-fractions du sol:
264
265	guez	49	loop_surface: DO nsrf = 1, nbsrf
266			! Chercher les indices :
267	guez	38	ni = 0
268			knon = 0
269			DO i = 1, klon
270	guez	145	! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
271	guez	38	! "potentielles"
272			IF (pctsrf_pot(i, nsrf) > epsfra) THEN
273			knon = knon + 1
274			ni(knon) = i
275			END IF
276			END DO
277	guez	15
278	guez	62	if_knon: IF (knon /= 0) then
279	guez	38	DO j = 1, knon
280			i = ni(j)
281	guez	62	ypct(j) = pctsrf(i, nsrf)
282	guez	207	yts(j) = ftsol(i, nsrf)
283	guez	215	snow(j) = fsnow(i, nsrf)
284	guez	62	yqsurf(j) = qsurf(i, nsrf)
285	guez	155	yalb(j) = falbe(i, nsrf)
286	guez	62	yrain_f(j) = rain_fall(i)
287			ysnow_f(j) = snow_f(i)
288			yagesno(j) = agesno(i, nsrf)
289	guez	222	yrugos(j) = frugs(i, nsrf)
290	guez	62	yrugoro(j) = rugoro(i)
291	guez	222	yrads(j) = fsolsw(i, nsrf) + fsollw(i, nsrf)
292	guez	225	ypaprs(j, klev + 1) = paprs(i, klev + 1)
293	guez	62	y_run_off_lic_0(j) = run_off_lic_0(i)
294	guez	38	END DO
295	guez	3
296	guez	99	! For continent, copy soil water content
297	guez	225	IF (nsrf == is_ter) yqsol(:knon) = qsol(ni(:knon))
298	guez	3
299	guez	208	ytsoil(:knon, :) = ftsoil(ni(:knon), :, nsrf)
300	guez	3
301	guez	38	DO k = 1, klev
302			DO j = 1, knon
303			i = ni(j)
304	guez	62	ypaprs(j, k) = paprs(i, k)
305			ypplay(j, k) = pplay(i, k)
306			ydelp(j, k) = delp(i, k)
307			yu(j, k) = u(i, k)
308			yv(j, k) = v(i, k)
309			yt(j, k) = t(i, k)
310			yq(j, k) = q(i, k)
311	guez	38	END DO
312			END DO
313	guez	3
314	guez	62	! calculer Cdrag et les coefficients d'echange
315	guez	221	CALL coefkz(nsrf, ypaprs, ypplay, ksta, ksta_ter, yts(:knon), &
316	guez	237	yrugos, yu, yv, yt, yq, yqsurf(:knon), coefm(:knon, :), &
317			coefh(:knon, :), ycdragm(:knon), ycdragh(:knon))
318	guez	228
319	guez	62	IF (iflag_pbl == 1) THEN
320	guez	235	CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0(:knon, 2:), &
321			ycoefh0(:knon, 2:))
322			ycoefm0(:knon, 1) = 0.
323			ycoefh0(:knon, 1) = 0.
324	guez	237	coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, 2:))
325			coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, 2:))
326			ycdragm(:knon) = max(ycdragm(:knon), ycoefm0(:knon, 1))
327			ycdragh(:knon) = max(ycdragh(:knon), ycoefh0(:knon, 1))
328	guez	62	END IF
329	guez	3
330	guez	237	! on met un seuil pour ycdragm et ycdragh
331	guez	62	IF (nsrf == is_oce) THEN
332	guez	237	ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
333			ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
334	guez	38	END IF
335	guez	3
336	guez	62	IF (ok_kzmin) THEN
337			! Calcul d'une diffusion minimale pour les conditions tres stables
338			CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
339	guez	237	ycdragm(:knon), ycoefm0(:knon, 2:), ycoefh0(:knon, 2:))
340			coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, 2:))
341			coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, 2:))
342			ycdragm(:knon) = max(ycdragm(:knon), ycoefm0(:knon, 1))
343			ycdragh(:knon) = max(ycdragh(:knon), ycoefh0(:knon, 1))
344	guez	98	END IF
345	guez	3
346	guez	228	IF (iflag_pbl >= 6) THEN
347	guez	145	! Mellor et Yamada adapt\'e \`a Mars, Richard Fournier et
348			! Fr\'ed\'eric Hourdin
349	guez	62	yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
350			+ ypplay(:knon, 1))) &
351			* (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg
352	guez	228
353	guez	62	DO k = 2, klev
354	guez	227	yzlay(:knon, k) = yzlay(:knon, k-1) &
355	guez	62	+ rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
356			/ ypaprs(1:knon, k) &
357			* (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
358			END DO
359	guez	227
360	guez	62	DO k = 1, klev
361	guez	225	yteta(1:knon, k) = yt(1:knon, k) * (ypaprs(1:knon, 1) &
362			/ ypplay(1:knon, k))*rkappa (1. + 0.61 * yq(1:knon, k))
363	guez	62	END DO
364	guez	227
365			zlev(:knon, 1) = 0.
366			zlev(:knon, klev + 1) = 2. * yzlay(:knon, klev) &
367	guez	62	- yzlay(:knon, klev - 1)
368	guez	227
369	guez	62	DO k = 2, klev
370	guez	227	zlev(:knon, k) = 0.5 * (yzlay(:knon, k) + yzlay(:knon, k-1))
371	guez	62	END DO
372	guez	227
373	guez	62	DO k = 1, klev + 1
374			DO j = 1, knon
375			i = ni(j)
376			yq2(j, k) = q2(i, k, nsrf)
377			END DO
378			END DO
379
380	guez	237	ustar(:knon) = ustarhb(yu(:knon, 1), yv(:knon, 1), ycdragm(:knon))
381	guez	228	CALL yamada4(dtime, rg, zlev(:knon, :), yzlay(:knon, :), &
382			yu(:knon, :), yv(:knon, :), yteta(:knon, :), &
383	guez	237	ycdragm(:knon), yq2(:knon, :), ykmm(:knon, :), &
384	guez	229	ykmn(:knon, :), ykmq(:knon, :), ustar(:knon))
385	guez	62	coefm(:knon, 2:) = ykmm(:knon, 2:klev)
386			coefh(:knon, 2:) = ykmn(:knon, 2:klev)
387	guez	38	END IF
388	guez	3
389	guez	237	CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), coefm(:knon, :), &
390			ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
391	guez	229	ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
392	guez	225	y_flux_u(:knon))
393	guez	237	CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), coefm(:knon, :), &
394			ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
395	guez	229	ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
396	guez	225	y_flux_v(:knon))
397	guez	3
398	guez	62	! calculer la diffusion de "q" et de "h"
399	guez	221	CALL clqh(dtime, julien, firstcal, nsrf, ni(:knon), &
400	guez	225	ytsoil(:knon, :), yqsol(:knon), mu0, yrugos, yrugoro, &
401	guez	237	yu(:knon, 1), yv(:knon, 1), coefh(:knon, :), ycdragh(:knon), &
402	guez	236	yt, yq, yts(:knon), ypaprs, ypplay, ydelp, yrads(:knon), &
403			yalb(:knon), snow(:knon), yqsurf, yrain_f, ysnow_f, &
404			yfluxlat(:knon), pctsrf_new_sic, yagesno(:knon), y_d_t, y_d_q, &
405			y_d_ts(:knon), yz0_new, y_flux_t(:knon), y_flux_q(:knon), &
406			y_dflux_t(:knon), y_dflux_q(:knon), y_fqcalving, y_ffonte, &
407			y_run_off_lic_0)
408	guez	3
409	guez	62	! calculer la longueur de rugosite sur ocean
410			yrugm = 0.
411			IF (nsrf == is_oce) THEN
412			DO j = 1, knon
413	guez	237	yrugm(j) = 0.018 * ycdragm(j) * (yu(j, 1)2 + yv(j, 1)2) &
414	guez	225	/ rg + 0.11 * 14E-6 &
415	guez	237	/ sqrt(ycdragm(j) * (yu(j, 1)2 + yv(j, 1)2))
416	guez	62	yrugm(j) = max(1.5E-05, yrugm(j))
417			END DO
418			END IF
419	guez	38	DO j = 1, knon
420	guez	225	y_dflux_t(j) = y_dflux_t(j) * ypct(j)
421			y_dflux_q(j) = y_dflux_q(j) * ypct(j)
422	guez	38	END DO
423	guez	3
424	guez	237	DO k = 2, klev
425	guez	62	DO j = 1, knon
426			i = ni(j)
427	guez	225	coefh(j, k) = coefh(j, k) * ypct(j)
428			coefm(j, k) = coefm(j, k) * ypct(j)
429	guez	237	END DO
430			END DO
431			DO j = 1, knon
432			i = ni(j)
433			ycdragh(j) = ycdragh(j) * ypct(j)
434			ycdragm(j) = ycdragm(j) * ypct(j)
435			END DO
436			DO k = 1, klev
437			DO j = 1, knon
438			i = ni(j)
439	guez	225	y_d_t(j, k) = y_d_t(j, k) * ypct(j)
440			y_d_q(j, k) = y_d_q(j, k) * ypct(j)
441			y_d_u(j, k) = y_d_u(j, k) * ypct(j)
442			y_d_v(j, k) = y_d_v(j, k) * ypct(j)
443	guez	62	END DO
444	guez	38	END DO
445	guez	3
446	guez	214	flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
447			flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
448			flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
449			flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
450	guez	15
451	guez	206	evap(:, nsrf) = -flux_q(:, nsrf)
452
453	guez	155	falbe(:, nsrf) = 0.
454	guez	215	fsnow(:, nsrf) = 0.
455	guez	62	qsurf(:, nsrf) = 0.
456	guez	222	frugs(:, nsrf) = 0.
457	guez	38	DO j = 1, knon
458			i = ni(j)
459	guez	62	d_ts(i, nsrf) = y_d_ts(j)
460	guez	155	falbe(i, nsrf) = yalb(j)
461	guez	215	fsnow(i, nsrf) = snow(j)
462	guez	62	qsurf(i, nsrf) = yqsurf(j)
463	guez	222	frugs(i, nsrf) = yz0_new(j)
464	guez	62	fluxlat(i, nsrf) = yfluxlat(j)
465			IF (nsrf == is_oce) THEN
466			rugmer(i) = yrugm(j)
467	guez	222	frugs(i, nsrf) = yrugm(j)
468	guez	62	END IF
469			agesno(i, nsrf) = yagesno(j)
470			fqcalving(i, nsrf) = y_fqcalving(j)
471			ffonte(i, nsrf) = y_ffonte(j)
472	guez	237	cdragh(i) = cdragh(i) + ycdragh(j)
473			cdragm(i) = cdragm(i) + ycdragm(j)
474	guez	62	dflux_t(i) = dflux_t(i) + y_dflux_t(j)
475			dflux_q(i) = dflux_q(i) + y_dflux_q(j)
476	guez	38	END DO
477	guez	62	IF (nsrf == is_ter) THEN
478	guez	99	qsol(ni(:knon)) = yqsol(:knon)
479			else IF (nsrf == is_lic) THEN
480	guez	62	DO j = 1, knon
481			i = ni(j)
482			run_off_lic_0(i) = y_run_off_lic_0(j)
483			END DO
484			END IF
485	guez	118
486	guez	62	ftsoil(:, :, nsrf) = 0.
487	guez	208	ftsoil(ni(:knon), :, nsrf) = ytsoil(:knon, :)
488	guez	62
489	guez	38	DO j = 1, knon
490			i = ni(j)
491	guez	62	DO k = 1, klev
492			d_t(i, k) = d_t(i, k) + y_d_t(j, k)
493			d_q(i, k) = d_q(i, k) + y_d_q(j, k)
494			d_u(i, k) = d_u(i, k) + y_d_u(j, k)
495			d_v(i, k) = d_v(i, k) + y_d_v(j, k)
496	guez	237	END DO
497			END DO
498
499			DO j = 1, knon
500			i = ni(j)
501			DO k = 2, klev
502	guez	70	ycoefh(i, k) = ycoefh(i, k) + coefh(j, k)
503	guez	62	END DO
504	guez	38	END DO
505	guez	62
506	guez	237	DO j = 1, knon
507			i = ni(j)
508			ycoefh(i, 1) = ycoefh(i, 1) + ycdragh(j)
509			END DO
510
511	guez	99	! diagnostic t, q a 2m et u, v a 10m
512	guez	62
513	guez	38	DO j = 1, knon
514			i = ni(j)
515	guez	227	u1(j) = yu(j, 1) + y_d_u(j, 1)
516			v1(j) = yv(j, 1) + y_d_v(j, 1)
517	guez	62	tair1(j) = yt(j, 1) + y_d_t(j, 1)
518			qair1(j) = yq(j, 1) + y_d_q(j, 1)
519	guez	225	zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
520			1))) * (ypaprs(j, 1)-ypplay(j, 1))
521	guez	62	tairsol(j) = yts(j) + y_d_ts(j)
522			rugo1(j) = yrugos(j)
523			IF (nsrf == is_oce) THEN
524	guez	222	rugo1(j) = frugs(i, nsrf)
525	guez	62	END IF
526			psfce(j) = ypaprs(j, 1)
527			patm(j) = ypplay(j, 1)
528	guez	15
529	guez	62	qairsol(j) = yqsurf(j)
530	guez	38	END DO
531	guez	15
532	guez	227	CALL stdlevvar(klon, knon, nsrf, u1(:knon), v1(:knon), tair1(:knon), &
533			qair1, zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, &
534			yq2m, yt10m, yq10m, wind10m(:knon), ustar)
535	guez	3
536	guez	62	DO j = 1, knon
537			i = ni(j)
538			t2m(i, nsrf) = yt2m(j)
539			q2m(i, nsrf) = yq2m(j)
540	guez	3
541	guez	227	u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
542			/ sqrt(u1(j)2 + v1(j)2)
543			v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
544			/ sqrt(u1(j)2 + v1(j)2)
545	guez	62	END DO
546	guez	15
547	guez	227	CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
548	guez	206	y_flux_q(:knon), yu, yv, yt, yq, ypblh(:knon), ycapcl, &
549			yoliqcl, ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)
550	guez	15
551	guez	38	DO j = 1, knon
552			i = ni(j)
553	guez	62	pblh(i, nsrf) = ypblh(j)
554			plcl(i, nsrf) = ylcl(j)
555			capcl(i, nsrf) = ycapcl(j)
556			oliqcl(i, nsrf) = yoliqcl(j)
557			cteicl(i, nsrf) = ycteicl(j)
558			pblt(i, nsrf) = ypblt(j)
559			therm(i, nsrf) = ytherm(j)
560			trmb1(i, nsrf) = ytrmb1(j)
561			trmb2(i, nsrf) = ytrmb2(j)
562			trmb3(i, nsrf) = ytrmb3(j)
563	guez	38	END DO
564	guez	3
565	guez	38	DO j = 1, knon
566	guez	62	DO k = 1, klev + 1
567			i = ni(j)
568			q2(i, k, nsrf) = yq2(j, k)
569			END DO
570	guez	38	END DO
571	guez	215	else
572			fsnow(:, nsrf) = 0.
573	guez	62	end IF if_knon
574	guez	49	END DO loop_surface
575	guez	15
576	guez	38	! On utilise les nouvelles surfaces
577	guez	222	frugs(:, is_oce) = rugmer
578	guez	202	pctsrf(:, is_oce) = pctsrf_new_oce
579			pctsrf(:, is_sic) = pctsrf_new_sic
580	guez	15
581	guez	202	firstcal = .false.
582
583	guez	38	END SUBROUTINE clmain
584	guez	15
585	guez	38	end module clmain_m