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, zu1, zv1, 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.

    ! Pour pouvoir extraire les coefficients d'\'echanges et le vent
    ! dans la premi\`ere couche, trois champs ont \'et\'e cr\'e\'es : "ycoefh",
    ! "zu1" et "zv1". Nous avons moyenn\'e les valeurs de ces trois
    ! champs sur les quatre sous-surfaces du mod\`ele.

    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: prt_level, 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 vdif_kcay_m, only: vdif_kcay
    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 à 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)
    REAL, intent(out):: zu1(klon), zv1(klon)
    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 u1lay(klon), v1lay(klon) ! vent dans la premi\`ere couche, pour
                              ! une sous-surface donnée
    
    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), yzlev(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 q2diag(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), yu10m(klon)
    REAL yustar(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 uzon(klon), vmer(klon)
    REAL tair1(klon), qair1(klon), tairsol(klon)
    REAL psfce(klon), patm(klon)

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

    ! utiliser un jeu de fonctions simples               
    LOGICAL zxli
    PARAMETER (zxli=.FALSE.)

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

    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.
    zu1 = 0.
    zv1 = 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)
             u1lay(j) = u(i, 1)
             v1lay(j) = v(i, 1)
             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 >= 3) 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(1:knon, k) = yzlay(1: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
             yzlev(1:knon, 1) = 0.
             yzlev(:knon, klev + 1) = 2. * yzlay(:knon, klev) &
                  - yzlay(:knon, klev - 1)
             DO k = 2, klev
                yzlev(1:knon, k) = 0.5 * (yzlay(1:knon, k) + yzlay(1: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

             CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)
             IF (prt_level > 9) PRINT *, 'USTAR = ', yustar

             ! iflag_pbl peut \^etre utilis\'e comme longueur de m\'elange

             IF (iflag_pbl >= 11) THEN
                CALL vdif_kcay(knon, dtime, rg, ypaprs, yzlev, yzlay, yu, yv, &
                     yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, yustar, &
                     iflag_pbl)
             ELSE
                CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &
                     coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
             END IF

             coefm(:knon, 2:) = ykmm(:knon, 2:klev)
             coefh(:knon, 2:) = ykmn(:knon, 2:klev)
          END IF

          ! calculer la diffusion des vitesses "u" et "v"
          CALL clvent(knon, dtime, u1lay(:knon), v1lay(:knon), &
               coefm(:knon, :), yt, yu, ypaprs, ypplay, ydelp, y_d_u, &
               y_flux_u(:knon))
          CALL clvent(knon, dtime, u1lay(:knon), v1lay(:knon), &
               coefm(:knon, :), yt, yv, ypaprs, ypplay, ydelp, y_d_v, &
               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, &
               u1lay(:knon), v1lay(:knon), 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) * (u1lay(j)**2 + v1lay(j)**2) &
                     / rg + 0.11 * 14E-6 &
                     / sqrt(coefm(j, 1) * (u1lay(j)**2 + v1lay(j)**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)
             zu1(i) = zu1(i) + u1lay(j) * ypct(j)
             zv1(i) = zv1(i) + v1lay(j) * ypct(j)
          END DO
          IF (nsrf == is_ter) THEN
             qsol(ni(:knon)) = yqsol(:knon)
          else IF (nsrf == is_lic) THEN
             DO j = 1, knon
                i = ni(j)
                run_off_lic_0(i) = y_run_off_lic_0(j)
             END DO
          END IF

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

          DO j = 1, knon
             i = ni(j)
             DO k = 1, klev
                d_t(i, k) = d_t(i, k) + y_d_t(j, k)
                d_q(i, k) = d_q(i, k) + y_d_q(j, k)
                d_u(i, k) = d_u(i, k) + y_d_u(j, k)
                d_v(i, k) = d_v(i, k) + y_d_v(j, k)
                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)
             uzon(j) = yu(j, 1) + y_d_u(j, 1)
             vmer(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, zxli, uzon(:knon), vmer(:knon), &
               tair1, qair1, zgeo1, tairsol, qairsol, rugo1, psfce, patm, &
               yt2m, yq2m, yt10m, yq10m, yu10m, yustar)

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

             u10m_srf(i, nsrf) = (yu10m(j) * uzon(j)) &
                  / sqrt(uzon(j)**2 + vmer(j)**2)
             v10m_srf(i, nsrf) = (yu10m(j) * vmer(j)) &
                  / sqrt(uzon(j)**2 + vmer(j)**2)
          END DO

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