Sources/phylmd/clmain.f

module clmain_m

  IMPLICIT NONE

contains

  SUBROUTINE clmain(dtime, pctsrf, t, q, u, v, jour, mu0, ftsol, cdmmax, &
       cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, qsol, paprs, pplay, fsnow, &
       qsurf, evap, falbe, fluxlat, rain_fall, snow_f, solsw, sollw, fder, &
       rugos, 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, lmt_pas
    USE conf_phys_m, ONLY: iflag_pbl
    USE dimphy, ONLY: klev, klon, zmasq
    USE dimsoil, ONLY: nsoilmx
    use hbtm_m, only: hbtm
    USE indicesol, ONLY: epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
    USE interfoce_lim_m, ONLY: interfoce_lim
    use stdlevvar_m, only: stdlevvar
    USE suphec_m, ONLY: rd, rg, rkappa
    use time_phylmdz, only: itap
    use ustarhb_m, only: ustarhb
    use vdif_kcay_m, only: vdif_kcay
    use yamada4_m, only: yamada4

    REAL, INTENT(IN):: dtime ! interval du temps (secondes)

    REAL, INTENT(inout):: pctsrf(klon, nbsrf)
    ! tableau des pourcentages de surface de chaque maille

    REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
    REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg/kg)
    REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
    INTEGER, INTENT(IN):: jour ! jour de l'annee en cours
    REAL, intent(in):: mu0(klon) ! cosinus de l'angle solaire zenithal     
    REAL, INTENT(IN):: ftsol(klon, nbsrf) ! temp\'erature du sol (en K)
    REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh
    REAL, INTENT(IN):: ksta, ksta_ter
    LOGICAL, INTENT(IN):: ok_kzmin

    REAL, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)
    ! soil temperature of surface fraction

    REAL, INTENT(inout):: qsol(klon)
    ! column-density of water in soil, in kg m-2

    REAL, INTENT(IN):: paprs(klon, klev+1) ! pression a intercouche (Pa)
    REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
    REAL, INTENT(inout):: fsnow(:, :) ! (klon, nbsrf) \'epaisseur neigeuse
    REAL qsurf(klon, nbsrf)
    REAL evap(klon, nbsrf)
    REAL, intent(inout):: falbe(klon, nbsrf)
    REAL, intent(out):: fluxlat(:, :) ! (klon, nbsrf)

    REAL, intent(in):: rain_fall(klon)
    ! liquid water mass flux (kg/m2/s), positive down

    REAL, intent(in):: snow_f(klon)
    ! solid water mass flux (kg/m2/s), positive down

    REAL, INTENT(IN):: solsw(klon, nbsrf), sollw(klon, nbsrf)
    REAL, intent(in):: fder(klon)
    REAL, intent(inout):: rugos(klon, nbsrf) ! longueur de rugosit\'e (en m)
    real agesno(klon, nbsrf)
    REAL, INTENT(IN):: rugoro(klon)

    REAL d_t(klon, klev), d_q(klon, klev)
    ! d_t------output-R- le changement pour "t"
    ! d_q------output-R- le changement pour "q"

    REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
    ! changement pour "u" et "v"

    REAL, intent(out):: d_ts(klon, nbsrf) ! le changement pour ftsol

    REAL, intent(out):: flux_t(klon, nbsrf)
    ! flux de chaleur sensible (Cp T) (W/m2) (orientation positive vers
    ! le bas) à la surface

    REAL, intent(out):: flux_q(klon, nbsrf) 
    ! flux de vapeur d'eau (kg/m2/s) à la surface

    REAL, intent(out):: flux_u(klon, nbsrf), flux_v(klon, nbsrf)
    ! tension du vent à la surface, en Pa

    REAL, INTENT(out):: cdragh(klon), cdragm(klon)
    real q2(klon, klev+1, nbsrf)

    REAL, INTENT(out):: dflux_t(klon), dflux_q(klon)
    ! dflux_t derive du flux sensible
    ! dflux_q derive du flux latent
    ! IM "slab" ocean

    REAL, intent(out):: ycoefh(klon, klev)
    REAL, intent(out):: zu1(klon)
    REAL zv1(klon)
    REAL t2m(klon, nbsrf), q2m(klon, nbsrf)
    REAL u10m(klon, nbsrf), v10m(klon, nbsrf)

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

    LOGICAL:: firstcal = .true.

    ! la nouvelle repartition des surfaces sortie de l'interface
    REAL, save:: pctsrf_new_oce(klon)
    REAL, save:: pctsrf_new_sic(klon)

    REAL y_fqcalving(klon), y_ffonte(klon)
    real y_run_off_lic_0(klon)
    REAL rugmer(klon)
    REAL ytsoil(klon, nsoilmx)
    REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
    REAL yalb(klon)

    REAL yu1(klon), yv1(klon)
    ! On ajoute en output yu1 et yv1 qui sont les vents dans
    ! la premi\`ere couche.
    
    REAL snow(klon), yqsurf(klon), yagesno(klon)

    real yqsol(klon)
    ! column-density of water in soil, in kg m-2

    REAL yrain_f(klon)
    ! liquid water mass flux (kg/m2/s), positive down

    REAL ysnow_f(klon)
    ! solid water mass flux (kg/m2/s), positive down

    REAL 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), y_flux_q(klon)
    REAL y_flux_u(klon), y_flux_v(klon)
    REAL y_dflux_t(klon), y_dflux_q(klon)
    REAL coefh(klon, klev), coefm(klon, klev)
    REAL yu(klon, klev), yv(klon, klev)
    REAL yt(klon, klev), yq(klon, klev)
    REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)

    REAL ycoefm0(klon, klev), ycoefh0(klon, klev)

    REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)
    REAL ykmm(klon, klev+1), ykmn(klon, klev+1)
    REAL ykmq(klon, klev+1)
    REAL yq2(klon, klev+1)
    REAL q2diag(klon, klev+1)

    REAL 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 extrapolation

    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.
    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.
    y_dflux_t = 0.
    y_dflux_q = 0.
    yrugoro = 0.
    d_ts = 0.
    flux_t = 0.
    flux_q = 0.
    flux_u = 0.
    flux_v = 0.
    fluxlat = 0.
    d_t = 0.
    d_q = 0.
    d_u = 0.
    d_v = 0.
    ycoefh = 0.

    ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
    ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
    ! (\`a affiner)

    pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
    pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
    pctsrf_pot(:, is_oce) = 1. - zmasq
    pctsrf_pot(:, is_sic) = 1. - zmasq

    ! Tester si c'est le moment de lire le fichier:
    if (mod(itap - 1, lmt_pas) == 0) then
       CALL interfoce_lim(jour, pctsrf_new_oce, pctsrf_new_sic)
    endif

    ! Boucler sur toutes les sous-fractions du sol:

    loop_surface: DO nsrf = 1, nbsrf
       ! Chercher les indices :
       ni = 0
       knon = 0
       DO i = 1, klon
          ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
          ! "potentielles"
          IF (pctsrf_pot(i, nsrf) > epsfra) THEN
             knon = knon + 1
             ni(knon) = i
          END IF
       END DO

       if_knon: IF (knon /= 0) then
          DO j = 1, knon
             i = ni(j)
             ypct(j) = pctsrf(i, nsrf)
             yts(j) = ftsol(i, nsrf)
             snow(j) = fsnow(i, nsrf)
             yqsurf(j) = qsurf(i, nsrf)
             yalb(j) = falbe(i, nsrf)
             yrain_f(j) = rain_fall(i)
             ysnow_f(j) = snow_f(i)
             yagesno(j) = agesno(i, nsrf)
             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

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

          DO k = 1, klev
             DO j = 1, knon
                i = ni(j)
                ypaprs(j, k) = paprs(i, k)
                ypplay(j, k) = pplay(i, k)
                ydelp(j, k) = delp(i, k)
                yu(j, k) = u(i, k)
                yv(j, k) = v(i, k)
                yt(j, k) = t(i, k)
                yq(j, k) = q(i, k)
             END DO
          END DO

          ! calculer Cdrag et les coefficients d'echange
          CALL coefkz(nsrf, ypaprs, ypplay, ksta, ksta_ter, yts, 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(:knon))
          CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &
               ypplay, ydelp, y_d_v, y_flux_v(:knon))

          ! calculer la diffusion de "q" et de "h"
          CALL clqh(dtime, jour, firstcal, nsrf, ni(:knon), ytsoil(:knon, :), &
               yqsol, mu0, yrugos, yrugoro, yu1, yv1, coefh(:knon, :), yt, &
               yq, yts(:knon), ypaprs, ypplay, ydelp, yrads, yalb(:knon), &
               snow(:knon), yqsurf, yrain_f, ysnow_f, yfder, yfluxlat(:knon), &
               pctsrf_new_sic, yagesno(:knon), y_d_t, y_d_q, y_d_ts(:knon), &
               yz0_new, y_flux_t(:knon), y_flux_q(:knon), y_dflux_t, &
               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)
                y_d_u(j, k) = y_d_u(j, k)*ypct(j)
                y_d_v(j, k) = y_d_v(j, k)*ypct(j)
             END DO
          END DO

          flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
          flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
          flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
          flux_v(ni(:knon), nsrf) = y_flux_v(:knon)

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

          falbe(:, nsrf) = 0.
          fsnow(:, nsrf) = 0.
          qsurf(:, nsrf) = 0.
          rugos(:, nsrf) = 0.
          DO j = 1, knon
             i = ni(j)
             d_ts(i, nsrf) = y_d_ts(j)
             falbe(i, nsrf) = yalb(j)
             fsnow(i, nsrf) = snow(j)
             qsurf(i, nsrf) = yqsurf(j)
             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.
          ftsoil(ni(:knon), :, nsrf) = ytsoil(:knon, :)

          DO j = 1, knon
             i = ni(j)
             DO k = 1, klev
                d_t(i, k) = d_t(i, k) + y_d_t(j, k)
                d_q(i, k) = d_q(i, k) + y_d_q(j, k)
                d_u(i, k) = d_u(i, k) + y_d_u(j, k)
                d_v(i, k) = d_v(i, k) + y_d_v(j, k)
                ycoefh(i, k) = ycoefh(i, k) + coefh(j, k)
             END DO
          END DO

          ! diagnostic t, q a 2m et u, v a 10m

          DO j = 1, knon
             i = ni(j)
             uzon(j) = yu(j, 1) + y_d_u(j, 1)
             vmer(j) = yv(j, 1) + y_d_v(j, 1)
             tair1(j) = yt(j, 1) + y_d_t(j, 1)
             qair1(j) = yq(j, 1) + y_d_q(j, 1)
             zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &
                  1)))*(ypaprs(j, 1)-ypplay(j, 1))
             tairsol(j) = yts(j) + y_d_ts(j)
             rugo1(j) = yrugos(j)
             IF (nsrf == is_oce) THEN
                rugo1(j) = 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(ypaprs, ypplay, yt2m, yq2m, yustar, y_flux_t(:knon), &
               y_flux_q(:knon), yu, yv, yt, yq, ypblh(:knon), ycapcl, &
               yoliqcl, ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)

          DO j = 1, knon
             i = ni(j)
             pblh(i, nsrf) = ypblh(j)
             plcl(i, nsrf) = ylcl(j)
             capcl(i, nsrf) = ycapcl(j)
             oliqcl(i, nsrf) = yoliqcl(j)
             cteicl(i, nsrf) = ycteicl(j)
             pblt(i, nsrf) = ypblt(j)
             therm(i, nsrf) = ytherm(j)
             trmb1(i, nsrf) = ytrmb1(j)
             trmb2(i, nsrf) = ytrmb2(j)
             trmb3(i, nsrf) = ytrmb3(j)
          END DO

          DO j = 1, knon
             DO k = 1, klev + 1
                i = ni(j)
                q2(i, k, nsrf) = yq2(j, k)
             END DO
          END DO
       else
          fsnow(:, nsrf) = 0.
       end IF if_knon
    END DO loop_surface

    ! On utilise les nouvelles surfaces
    rugos(:, is_oce) = rugmer
    pctsrf(:, is_oce) = pctsrf_new_oce
    pctsrf(:, is_sic) = pctsrf_new_sic

    firstcal = .false.

  END SUBROUTINE clmain

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