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 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, :), &
               coefh(:knon, :))

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