Sources/phylmd/clmain.f

module clmain_m

  IMPLICIT NONE

contains

  SUBROUTINE clmain(dtime, itap, pctsrf, pctsrf_new, t, q, u, v, jour, rmu0, &
       co2_ppm, 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 dimens_m, ONLY: iim, jjm
    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):: co2_ppm ! taux CO2 atmosphere
    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)
    ! -- LOOP
    REAL yu10mx(klon)
    REAL yu10my(klon)
    REAL ywindsp(klon)
    ! -- LOOP

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