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)

    !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds
    ! physiq ce qui permet de sortir les grdeurs 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, 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			!IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds
128			! physiq ce qui permet de sortir les grdeurs par sous surface)
129			REAL pblh(klon, nbsrf)
130			! pblh------- HCL
131			REAL capcl(klon, nbsrf)
132			REAL oliqcl(klon, nbsrf)
133			REAL cteicl(klon, nbsrf)
134			REAL pblt(klon, nbsrf)
135			! pblT------- T au nveau HCL
136			REAL therm(klon, nbsrf)
137			REAL trmb1(klon, nbsrf)
138			! trmb1-------deep_cape
139			REAL trmb2(klon, nbsrf)
140			! trmb2--------inhibition
141			REAL trmb3(klon, nbsrf)
142			! trmb3-------Point Omega
143			REAL plcl(klon, nbsrf)
144			REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)
145			! ffonte----Flux thermique utilise pour fondre la neige
146			! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
147			! hauteur de neige, en kg/m2/s
148			REAL run_off_lic_0(klon)
149
150			! Local:
151	guez	15
152	guez	70	REAL y_fqcalving(klon), y_ffonte(klon)
153			real y_run_off_lic_0(klon)
154	guez	15
155	guez	70	REAL rugmer(klon)
156	guez	15
157	guez	38	REAL ytsoil(klon, nsoilmx)
158	guez	15
159	guez	38	REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
160			REAL yalb(klon)
161			REAL yu1(klon), yv1(klon)
162	guez	49	! on rajoute en output yu1 et yv1 qui sont les vents dans
163			! la premiere couche
164	guez	101	REAL ysnow(klon), yqsurf(klon), yagesno(klon)
165
166			real yqsol(klon)
167			! column-density of water in soil, in kg m-2
168
169			REAL yrain_f(klon)
170			! liquid water mass flux (kg/m2/s), positive down
171
172			REAL ysnow_f(klon)
173			! solid water mass flux (kg/m2/s), positive down
174
175	guez	99	REAL yfder(klon)
176	guez	38	REAL yrugm(klon), yrads(klon), yrugoro(klon)
177	guez	15
178	guez	38	REAL yfluxlat(klon)
179	guez	15
180	guez	38	REAL y_d_ts(klon)
181			REAL y_d_t(klon, klev), y_d_q(klon, klev)
182			REAL y_d_u(klon, klev), y_d_v(klon, klev)
183			REAL y_flux_t(klon, klev), y_flux_q(klon, klev)
184			REAL y_flux_u(klon, klev), y_flux_v(klon, klev)
185			REAL y_dflux_t(klon), y_dflux_q(klon)
186	guez	62	REAL coefh(klon, klev), coefm(klon, klev)
187	guez	38	REAL yu(klon, klev), yv(klon, klev)
188			REAL yt(klon, klev), yq(klon, klev)
189			REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)
190	guez	15
191	guez	38	REAL ycoefm0(klon, klev), ycoefh0(klon, klev)
192	guez	15
193	guez	38	REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)
194			REAL ykmm(klon, klev+1), ykmn(klon, klev+1)
195			REAL ykmq(klon, klev+1)
196	guez	70	REAL yq2(klon, klev+1)
197	guez	38	REAL q2diag(klon, klev+1)
198	guez	15
199	guez	38	REAL u1lay(klon), v1lay(klon)
200			REAL delp(klon, klev)
201			INTEGER i, k, nsrf
202	guez	15
203	guez	38	INTEGER ni(klon), knon, j
204	guez	40
205	guez	38	REAL pctsrf_pot(klon, nbsrf)
206	guez	145	! "pourcentage potentiel" pour tenir compte des \'eventuelles
207	guez	40	! apparitions ou disparitions de la glace de mer
208	guez	15
209	guez	38	REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
210	guez	15
211	guez	38	REAL yt2m(klon), yq2m(klon), yu10m(klon)
212			REAL yustar(klon)
213	guez	15
214	guez	38	REAL yt10m(klon), yq10m(klon)
215			REAL ypblh(klon)
216			REAL ylcl(klon)
217			REAL ycapcl(klon)
218			REAL yoliqcl(klon)
219			REAL ycteicl(klon)
220			REAL ypblt(klon)
221			REAL ytherm(klon)
222			REAL ytrmb1(klon)
223			REAL ytrmb2(klon)
224			REAL ytrmb3(klon)
225			REAL uzon(klon), vmer(klon)
226			REAL tair1(klon), qair1(klon), tairsol(klon)
227			REAL psfce(klon), patm(klon)
228	guez	15
229	guez	38	REAL qairsol(klon), zgeo1(klon)
230			REAL rugo1(klon)
231	guez	15
232	guez	38	! utiliser un jeu de fonctions simples
233			LOGICAL zxli
234			PARAMETER (zxli=.FALSE.)
235	guez	15
236	guez	38	!------------------------------------------------------------
237	guez	15
238	guez	38	ytherm = 0.
239	guez	15
240	guez	38	DO k = 1, klev ! epaisseur de couche
241			DO i = 1, klon
242			delp(i, k) = paprs(i, k) - paprs(i, k+1)
243			END DO
244			END DO
245			DO i = 1, klon ! vent de la premiere couche
246			zx_alf1 = 1.0
247			zx_alf2 = 1.0 - zx_alf1
248			u1lay(i) = u(i, 1)zx_alf1 + u(i, 2)zx_alf2
249			v1lay(i) = v(i, 1)zx_alf1 + v(i, 2)zx_alf2
250			END DO
251	guez	15
252	guez	40	! Initialization:
253			rugmer = 0.
254			cdragh = 0.
255			cdragm = 0.
256			dflux_t = 0.
257			dflux_q = 0.
258			zu1 = 0.
259			zv1 = 0.
260			ypct = 0.
261			yts = 0.
262			ysnow = 0.
263			yqsurf = 0.
264			yrain_f = 0.
265			ysnow_f = 0.
266			yfder = 0.
267			yrugos = 0.
268			yu1 = 0.
269			yv1 = 0.
270			yrads = 0.
271			ypaprs = 0.
272			ypplay = 0.
273			ydelp = 0.
274			yu = 0.
275			yv = 0.
276			yt = 0.
277			yq = 0.
278			pctsrf_new = 0.
279			y_flux_u = 0.
280			y_flux_v = 0.
281			y_dflux_t = 0.
282			y_dflux_q = 0.
283	guez	38	ytsoil = 999999.
284			yrugoro = 0.
285	guez	40	d_ts = 0.
286	guez	38	yfluxlat = 0.
287			flux_t = 0.
288			flux_q = 0.
289			flux_u = 0.
290			flux_v = 0.
291	guez	40	d_t = 0.
292			d_q = 0.
293			d_u = 0.
294			d_v = 0.
295	guez	70	ycoefh = 0.
296	guez	15
297	guez	145	! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
298			! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
299			! (\`a affiner)
300	guez	15
301	guez	38	pctsrf_pot = pctsrf
302			pctsrf_pot(:, is_oce) = 1. - zmasq
303			pctsrf_pot(:, is_sic) = 1. - zmasq
304	guez	15
305	guez	99	! Boucler sur toutes les sous-fractions du sol:
306
307	guez	49	loop_surface: DO nsrf = 1, nbsrf
308			! Chercher les indices :
309	guez	38	ni = 0
310			knon = 0
311			DO i = 1, klon
312	guez	145	! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
313	guez	38	! "potentielles"
314			IF (pctsrf_pot(i, nsrf) > epsfra) THEN
315			knon = knon + 1
316			ni(knon) = i
317			END IF
318			END DO
319	guez	15
320	guez	62	if_knon: IF (knon /= 0) then
321	guez	38	DO j = 1, knon
322			i = ni(j)
323	guez	62	ypct(j) = pctsrf(i, nsrf)
324			yts(j) = ts(i, nsrf)
325			ysnow(j) = snow(i, nsrf)
326			yqsurf(j) = qsurf(i, nsrf)
327	guez	155	yalb(j) = falbe(i, nsrf)
328	guez	62	yrain_f(j) = rain_fall(i)
329			ysnow_f(j) = snow_f(i)
330			yagesno(j) = agesno(i, nsrf)
331			yfder(j) = fder(i)
332			yrugos(j) = rugos(i, nsrf)
333			yrugoro(j) = rugoro(i)
334			yu1(j) = u1lay(i)
335			yv1(j) = v1lay(i)
336	guez	175	yrads(j) = solsw(i, nsrf) + sollw(i, nsrf)
337	guez	62	ypaprs(j, klev+1) = paprs(i, klev+1)
338			y_run_off_lic_0(j) = run_off_lic_0(i)
339	guez	38	END DO
340	guez	3
341	guez	99	! For continent, copy soil water content
342			IF (nsrf == is_ter) THEN
343			yqsol(:knon) = qsol(ni(:knon))
344	guez	62	ELSE
345			yqsol = 0.
346			END IF
347	guez	3
348	guez	62	DO k = 1, nsoilmx
349			DO j = 1, knon
350			i = ni(j)
351			ytsoil(j, k) = ftsoil(i, k, nsrf)
352	guez	38	END DO
353	guez	47	END DO
354	guez	3
355	guez	38	DO k = 1, klev
356			DO j = 1, knon
357			i = ni(j)
358	guez	62	ypaprs(j, k) = paprs(i, k)
359			ypplay(j, k) = pplay(i, k)
360			ydelp(j, k) = delp(i, k)
361			yu(j, k) = u(i, k)
362			yv(j, k) = v(i, k)
363			yt(j, k) = t(i, k)
364			yq(j, k) = q(i, k)
365	guez	38	END DO
366			END DO
367	guez	3
368	guez	62	! calculer Cdrag et les coefficients d'echange
369			CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, &
370			yu, yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))
371			IF (iflag_pbl == 1) THEN
372			CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
373			coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
374			coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
375			END IF
376	guez	3
377	guez	70	! on met un seuil pour coefm et coefh
378	guez	62	IF (nsrf == is_oce) THEN
379			coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)
380			coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)
381	guez	38	END IF
382	guez	3
383	guez	62	IF (ok_kzmin) THEN
384			! Calcul d'une diffusion minimale pour les conditions tres stables
385			CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
386	guez	70	coefm(:knon, 1), ycoefm0, ycoefh0)
387	guez	62	coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
388			coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
389	guez	98	END IF
390	guez	3
391	guez	62	IF (iflag_pbl >= 3) THEN
392	guez	145	! Mellor et Yamada adapt\'e \`a Mars, Richard Fournier et
393			! Fr\'ed\'eric Hourdin
394	guez	62	yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
395			+ ypplay(:knon, 1))) &
396			* (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg
397			DO k = 2, klev
398			yzlay(1:knon, k) = yzlay(1:knon, k-1) &
399			+ rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
400			/ ypaprs(1:knon, k) &
401			* (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
402			END DO
403			DO k = 1, klev
404			yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
405			/ ypplay(1:knon, k))*rkappa (1.+0.61*yq(1:knon, k))
406			END DO
407			yzlev(1:knon, 1) = 0.
408			yzlev(:knon, klev+1) = 2. * yzlay(:knon, klev) &
409			- yzlay(:knon, klev - 1)
410			DO k = 2, klev
411			yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
412			END DO
413			DO k = 1, klev + 1
414			DO j = 1, knon
415			i = ni(j)
416			yq2(j, k) = q2(i, k, nsrf)
417			END DO
418			END DO
419
420			CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)
421	guez	99	IF (prt_level > 9) PRINT *, 'USTAR = ', yustar
422	guez	62
423	guez	145	! iflag_pbl peut \^etre utilis\'e comme longueur de m\'elange
424	guez	62
425			IF (iflag_pbl >= 11) THEN
426	guez	145	CALL vdif_kcay(knon, dtime, rg, ypaprs, yzlev, yzlay, yu, yv, &
427			yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, yustar, &
428			iflag_pbl)
429	guez	62	ELSE
430			CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &
431			coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
432			END IF
433
434			coefm(:knon, 2:) = ykmm(:knon, 2:klev)
435			coefh(:knon, 2:) = ykmn(:knon, 2:klev)
436	guez	38	END IF
437	guez	3
438	guez	62	! calculer la diffusion des vitesses "u" et "v"
439	guez	70	CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &
440			ypplay, ydelp, y_d_u, y_flux_u)
441			CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &
442			ypplay, ydelp, y_d_v, y_flux_v)
443	guez	3
444	guez	62	! calculer la diffusion de "q" et de "h"
445	guez	118	CALL clqh(dtime, itap, jour, debut, rlat, knon, nsrf, ni(:knon), &
446	guez	178	pctsrf, ytsoil, yqsol, rmu0, yrugos, yrugoro, yu1, &
447	guez	118	yv1, coefh(:knon, :), yt, yq, yts, ypaprs, ypplay, ydelp, &
448	guez	155	yrads, yalb(:knon), ysnow, yqsurf, yrain_f, ysnow_f, yfder, &
449	guez	175	yfluxlat, pctsrf_new, yagesno(:knon), y_d_t, y_d_q, &
450	guez	118	y_d_ts(:knon), yz0_new, y_flux_t, y_flux_q, y_dflux_t, &
451	guez	175	y_dflux_q, y_fqcalving, y_ffonte, y_run_off_lic_0)
452	guez	3
453	guez	62	! calculer la longueur de rugosite sur ocean
454			yrugm = 0.
455			IF (nsrf == is_oce) THEN
456			DO j = 1, knon
457			yrugm(j) = 0.018coefm(j, 1)(yu1(j)2+yv1(j)2)/rg + &
458			0.1114E-6/sqrt(coefm(j, 1)(yu1(j)2+yv1(j)2))
459			yrugm(j) = max(1.5E-05, yrugm(j))
460			END DO
461			END IF
462	guez	38	DO j = 1, knon
463	guez	62	y_dflux_t(j) = y_dflux_t(j)*ypct(j)
464			y_dflux_q(j) = y_dflux_q(j)*ypct(j)
465			yu1(j) = yu1(j)*ypct(j)
466			yv1(j) = yv1(j)*ypct(j)
467	guez	38	END DO
468	guez	3
469	guez	62	DO k = 1, klev
470			DO j = 1, knon
471			i = ni(j)
472			coefh(j, k) = coefh(j, k)*ypct(j)
473			coefm(j, k) = coefm(j, k)*ypct(j)
474			y_d_t(j, k) = y_d_t(j, k)*ypct(j)
475			y_d_q(j, k) = y_d_q(j, k)*ypct(j)
476			flux_t(i, k, nsrf) = y_flux_t(j, k)
477			flux_q(i, k, nsrf) = y_flux_q(j, k)
478			flux_u(i, k, nsrf) = y_flux_u(j, k)
479			flux_v(i, k, nsrf) = y_flux_v(j, k)
480			y_d_u(j, k) = y_d_u(j, k)*ypct(j)
481			y_d_v(j, k) = y_d_v(j, k)*ypct(j)
482			END DO
483	guez	38	END DO
484	guez	3
485	guez	62	evap(:, nsrf) = -flux_q(:, 1, nsrf)
486	guez	15
487	guez	155	falbe(:, nsrf) = 0.
488	guez	62	snow(:, nsrf) = 0.
489			qsurf(:, nsrf) = 0.
490			rugos(:, nsrf) = 0.
491			fluxlat(:, nsrf) = 0.
492	guez	38	DO j = 1, knon
493			i = ni(j)
494	guez	62	d_ts(i, nsrf) = y_d_ts(j)
495	guez	155	falbe(i, nsrf) = yalb(j)
496	guez	62	snow(i, nsrf) = ysnow(j)
497			qsurf(i, nsrf) = yqsurf(j)
498			rugos(i, nsrf) = yz0_new(j)
499			fluxlat(i, nsrf) = yfluxlat(j)
500			IF (nsrf == is_oce) THEN
501			rugmer(i) = yrugm(j)
502			rugos(i, nsrf) = yrugm(j)
503			END IF
504			agesno(i, nsrf) = yagesno(j)
505			fqcalving(i, nsrf) = y_fqcalving(j)
506			ffonte(i, nsrf) = y_ffonte(j)
507			cdragh(i) = cdragh(i) + coefh(j, 1)
508			cdragm(i) = cdragm(i) + coefm(j, 1)
509			dflux_t(i) = dflux_t(i) + y_dflux_t(j)
510			dflux_q(i) = dflux_q(i) + y_dflux_q(j)
511			zu1(i) = zu1(i) + yu1(j)
512			zv1(i) = zv1(i) + yv1(j)
513	guez	38	END DO
514	guez	62	IF (nsrf == is_ter) THEN
515	guez	99	qsol(ni(:knon)) = yqsol(:knon)
516			else IF (nsrf == is_lic) THEN
517	guez	62	DO j = 1, knon
518			i = ni(j)
519			run_off_lic_0(i) = y_run_off_lic_0(j)
520			END DO
521			END IF
522	guez	118
523	guez	62	ftsoil(:, :, nsrf) = 0.
524			DO k = 1, nsoilmx
525			DO j = 1, knon
526			i = ni(j)
527			ftsoil(i, k, nsrf) = ytsoil(j, k)
528			END DO
529			END DO
530
531	guez	38	DO j = 1, knon
532			i = ni(j)
533	guez	62	DO k = 1, klev
534			d_t(i, k) = d_t(i, k) + y_d_t(j, k)
535			d_q(i, k) = d_q(i, k) + y_d_q(j, k)
536			d_u(i, k) = d_u(i, k) + y_d_u(j, k)
537			d_v(i, k) = d_v(i, k) + y_d_v(j, k)
538	guez	70	ycoefh(i, k) = ycoefh(i, k) + coefh(j, k)
539	guez	62	END DO
540	guez	38	END DO
541	guez	62
542	guez	99	! diagnostic t, q a 2m et u, v a 10m
543	guez	62
544	guez	38	DO j = 1, knon
545			i = ni(j)
546	guez	62	uzon(j) = yu(j, 1) + y_d_u(j, 1)
547			vmer(j) = yv(j, 1) + y_d_v(j, 1)
548			tair1(j) = yt(j, 1) + y_d_t(j, 1)
549			qair1(j) = yq(j, 1) + y_d_q(j, 1)
550			zgeo1(j) = rdtair1(j)/(0.5(ypaprs(j, 1)+ypplay(j, &
551			1)))*(ypaprs(j, 1)-ypplay(j, 1))
552			tairsol(j) = yts(j) + y_d_ts(j)
553			rugo1(j) = yrugos(j)
554			IF (nsrf == is_oce) THEN
555			rugo1(j) = rugos(i, nsrf)
556			END IF
557			psfce(j) = ypaprs(j, 1)
558			patm(j) = ypplay(j, 1)
559	guez	15
560	guez	62	qairsol(j) = yqsurf(j)
561	guez	38	END DO
562	guez	15
563	guez	62	CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, &
564			zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &
565			yt10m, yq10m, yu10m, yustar)
566	guez	3
567	guez	62	DO j = 1, knon
568			i = ni(j)
569			t2m(i, nsrf) = yt2m(j)
570			q2m(i, nsrf) = yq2m(j)
571	guez	3
572	guez	62	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
573			u10m(i, nsrf) = (yu10m(j)uzon(j))/sqrt(uzon(j)2+vmer(j)*2)
574			v10m(i, nsrf) = (yu10m(j)vmer(j))/sqrt(uzon(j)2+vmer(j)*2)
575	guez	15
576	guez	62	END DO
577	guez	15
578	guez	149	CALL hbtm(knon, ypaprs, ypplay, yt2m, yq2m, yustar, &
579	guez	62	y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &
580			ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)
581	guez	15
582	guez	38	DO j = 1, knon
583			i = ni(j)
584	guez	62	pblh(i, nsrf) = ypblh(j)
585			plcl(i, nsrf) = ylcl(j)
586			capcl(i, nsrf) = ycapcl(j)
587			oliqcl(i, nsrf) = yoliqcl(j)
588			cteicl(i, nsrf) = ycteicl(j)
589			pblt(i, nsrf) = ypblt(j)
590			therm(i, nsrf) = ytherm(j)
591			trmb1(i, nsrf) = ytrmb1(j)
592			trmb2(i, nsrf) = ytrmb2(j)
593			trmb3(i, nsrf) = ytrmb3(j)
594	guez	38	END DO
595	guez	3
596	guez	38	DO j = 1, knon
597	guez	62	DO k = 1, klev + 1
598			i = ni(j)
599			q2(i, k, nsrf) = yq2(j, k)
600			END DO
601	guez	38	END DO
602	guez	62	end IF if_knon
603	guez	49	END DO loop_surface
604	guez	15
605	guez	38	! On utilise les nouvelles surfaces
606	guez	15
607	guez	38	rugos(:, is_oce) = rugmer
608			pctsrf = pctsrf_new
609	guez	15
610	guez	38	END SUBROUTINE clmain
611	guez	15
612	guez	38	end module clmain_m