Sources/phylmd/clmain.f

module clmain_m

  IMPLICIT NONE

contains

  SUBROUTINE clmain(dtime, pctsrf, 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, 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):: 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, INTENT(inout):: 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, 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 "ts"

    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 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), 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.
    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.
    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(:, 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) = 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(: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, rlat, nsrf, ni(:knon), 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_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

          DO j = 1, knon
             i = ni(j)
             flux_t(i, nsrf) = y_flux_t(j)
             flux_q(i, nsrf) = y_flux_q(j)
             flux_u(i, nsrf) = y_flux_u(j)
             flux_v(i, nsrf) = y_flux_v(j)
          END DO

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