Sources/phylmd/clmain.f

module clmain_m

  IMPLICIT NONE

contains

  SUBROUTINE clmain(dtime, itap, pctsrf, pctsrf_new, t, q, u, v, jour, rmu0, &
       ts, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, qsol, &
       paprs, pplay, snow, qsurf, evap, falbe, fluxlat, rain_fall, snow_f, &
       solsw, sollw, fder, rlat, rugos, debut, 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, v10m, 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
    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 stdlevvar_m, only: stdlevvar
    USE suphec_m, ONLY: rd, rg, rkappa
    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)
    INTEGER, INTENT(IN):: itap ! numero du pas de temps
    REAL, INTENT(inout):: pctsrf(klon, nbsrf)

    ! la nouvelle repartition des surfaces sortie de l'interface
    REAL, INTENT(out):: pctsrf_new(klon, nbsrf)

    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):: jour ! jour de l'annee en cours
    REAL, intent(in):: rmu0(klon) ! cosinus de l'angle solaire zenithal     
    REAL, INTENT(IN):: ts(klon, nbsrf) ! temperature du sol (en Kelvin)
    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 snow(klon, nbsrf)
    REAL qsurf(klon, nbsrf)
    REAL evap(klon, nbsrf)
    REAL, intent(inout):: falbe(klon, nbsrf)

    REAL 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):: solsw(klon, nbsrf), sollw(klon, nbsrf)
    REAL, intent(in):: fder(klon)
    REAL, INTENT(IN):: rlat(klon) ! latitude en degr\'es

    REAL rugos(klon, nbsrf)
    ! rugos----input-R- longeur de rugosite (en m)

    LOGICAL, INTENT(IN):: debut
    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) ! le changement pour "ts"

    REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)
    ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
    !                    (orientation positive vers le bas)
    ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)

    REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)
    ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
    ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal

    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)
    REAL zv1(klon)
    REAL t2m(klon, nbsrf), q2m(klon, nbsrf)
    REAL u10m(klon, nbsrf), v10m(klon, nbsrf)

    ! Ionela Musat cf. Anne Mathieu : pbl, hbtm (Comme les autres
    ! diagnostics on cumule dans physiq ce qui permet de sortir les
    ! grandeurs par sous-surface)
    REAL pblh(klon, nbsrf)
    ! pblh------- HCL
    REAL capcl(klon, nbsrf)
    REAL oliqcl(klon, nbsrf)
    REAL cteicl(klon, nbsrf)
    REAL pblt(klon, nbsrf)
    ! pblT------- 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:

    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 yu1(klon), yv1(klon)
    ! on rajoute en output yu1 et yv1 qui sont les vents dans
    ! la premiere couche
    REAL ysnow(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 yfder(klon)
    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, klev), y_flux_q(klon, klev)
    REAL y_flux_u(klon, klev), y_flux_v(klon, klev)
    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 u1lay(klon), v1lay(klon)
    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 zx_alf1, zx_alf2 !valeur ambiante par extrapola.

    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
    DO i = 1, klon ! vent de la premiere couche
       zx_alf1 = 1.0
       zx_alf2 = 1.0 - zx_alf1
       u1lay(i) = u(i, 1)*zx_alf1 + u(i, 2)*zx_alf2
       v1lay(i) = v(i, 1)*zx_alf1 + v(i, 2)*zx_alf2
    END DO

    ! Initialization:
    rugmer = 0.
    cdragh = 0.
    cdragm = 0.
    dflux_t = 0.
    dflux_q = 0.
    zu1 = 0.
    zv1 = 0.
    ypct = 0.
    yts = 0.
    ysnow = 0.
    yqsurf = 0.
    yrain_f = 0.
    ysnow_f = 0.
    yfder = 0.
    yrugos = 0.
    yu1 = 0.
    yv1 = 0.
    yrads = 0.
    ypaprs = 0.
    ypplay = 0.
    ydelp = 0.
    yu = 0.
    yv = 0.
    yt = 0.
    yq = 0.
    pctsrf_new = 0.
    y_flux_u = 0.
    y_flux_v = 0.
    y_dflux_t = 0.
    y_dflux_q = 0.
    ytsoil = 999999.
    yrugoro = 0.
    d_ts = 0.
    yfluxlat = 0.
    flux_t = 0.
    flux_q = 0.
    flux_u = 0.
    flux_v = 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 = pctsrf
    pctsrf_pot(:, is_oce) = 1. - zmasq
    pctsrf_pot(:, is_sic) = 1. - zmasq

    ! 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) = ts(i, nsrf)
             ysnow(j) = snow(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)
             yfder(j) = fder(i)
             yrugos(j) = rugos(i, nsrf)
             yrugoro(j) = rugoro(i)
             yu1(j) = u1lay(i)
             yv1(j) = v1lay(i)
             yrads(j) = solsw(i, nsrf) + sollw(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) THEN
             yqsol(:knon) = qsol(ni(:knon))
          ELSE
             yqsol = 0.
          END IF

          DO k = 1, nsoilmx
             DO j = 1, knon
                i = ni(j)
                ytsoil(j, k) = ftsoil(i, k, nsrf)
             END DO
          END DO

          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, knon, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, &
               yu, yv, yt, yq, yqsurf, 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, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &
               ypplay, ydelp, y_d_u, y_flux_u)
          CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &
               ypplay, ydelp, y_d_v, y_flux_v)

          ! calculer la diffusion de "q" et de "h"
          CALL clqh(dtime, itap, jour, debut, rlat, knon, nsrf, ni(:knon), &
               pctsrf, ytsoil, yqsol, rmu0, yrugos, yrugoro, yu1, &
               yv1, coefh(:knon, :), yt, yq, yts, ypaprs, ypplay, ydelp, &
               yrads, yalb(:knon), ysnow, yqsurf, yrain_f, ysnow_f, yfder, &
               yfluxlat, pctsrf_new, yagesno(:knon), y_d_t, y_d_q, &
               y_d_ts(:knon), yz0_new, y_flux_t, y_flux_q, y_dflux_t, &
               y_dflux_q, 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)*(yu1(j)**2+yv1(j)**2)/rg + &
                     0.11*14E-6/sqrt(coefm(j, 1)*(yu1(j)**2+yv1(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)
             yu1(j) = yu1(j)*ypct(j)
             yv1(j) = yv1(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)
                flux_t(i, k, nsrf) = y_flux_t(j, k)
                flux_q(i, k, nsrf) = y_flux_q(j, k)
                flux_u(i, k, nsrf) = y_flux_u(j, k)
                flux_v(i, k, nsrf) = y_flux_v(j, k)
                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

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

          falbe(:, nsrf) = 0.
          snow(:, nsrf) = 0.
          qsurf(:, nsrf) = 0.
          rugos(:, nsrf) = 0.
          fluxlat(:, nsrf) = 0.
          DO j = 1, knon
             i = ni(j)
             d_ts(i, nsrf) = y_d_ts(j)
             falbe(i, nsrf) = yalb(j)
             snow(i, nsrf) = ysnow(j)
             qsurf(i, nsrf) = yqsurf(j)
             rugos(i, nsrf) = yz0_new(j)
             fluxlat(i, nsrf) = yfluxlat(j)
             IF (nsrf == is_oce) THEN
                rugmer(i) = yrugm(j)
                rugos(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) + yu1(j)
             zv1(i) = zv1(i) + yv1(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.
          DO k = 1, nsoilmx
             DO j = 1, knon
                i = ni(j)
                ftsoil(i, k, nsrf) = ytsoil(j, k)
             END DO
          END DO

          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) = rugos(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, vmer, 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, v10m : composantes du vent a 10m sans spirale de Ekman
             u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)
             v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)

          END DO

          CALL hbtm(knon, ypaprs, ypplay, yt2m, yq2m, yustar, y_flux_t, &
               y_flux_q, 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
       end IF if_knon
    END DO loop_surface

    ! On utilise les nouvelles surfaces

    rugos(:, is_oce) = rugmer
    pctsrf = pctsrf_new

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