trunk/phylmd/pbl_surface.f

module clmain_m

  IMPLICIT NONE

contains

  SUBROUTINE clmain(dtime, pctsrf, t, q, u, v, julien, mu0, ftsol, cdmmax, &
       cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, qsol, paprs, pplay, fsnow, &
       qsurf, evap, falbe, fluxlat, rain_fall, snow_f, fsolsw, fsollw, frugs, &
       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, coefh, t2m, q2m, &
       u10m_srf, v10m_srf, 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.

    use clcdrag_m, only: clcdrag
    use clqh_m, only: clqh
    use clvent_m, only: clvent
    use coefkz_m, only: coefkz
    use coefkzmin_m, only: coefkzmin
    use coefkz2_m, only: coefkz2
    USE conf_gcm_m, ONLY: 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 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):: julien ! 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):: fsolsw(klon, nbsrf), fsollw(klon, nbsrf)
    REAL, intent(inout):: frugs(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) variation of 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 (flux turbulent de 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):: coefh(:, 2:) ! (klon, 2:klev)
    ! Pour pouvoir extraire les coefficients d'\'echange, le champ
    ! "coefh" a \'et\'e cr\'e\'e. Nous avons moyenn\'e les valeurs de
    ! ce champ sur les quatre sous-surfaces du mod\`ele.

    REAL, INTENT(inout):: t2m(klon, nbsrf), q2m(klon, nbsrf)

    REAL, INTENT(inout):: u10m_srf(:, :), v10m_srf(:, :) ! (klon, nbsrf)
    ! composantes du vent \`a 10m sans spirale d'Ekman

    ! 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, INTENT(inout):: pblt(klon, nbsrf) ! 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), ypct(klon), yz0_new(klon)
    real yrugos(klon) ! longeur de rugosite (en m)
    REAL yalb(klon)
    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 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 ycoefh(klon, 2:klev), ycoefm(klon, 2:klev)
    real ycdragh(klon), ycdragm(klon)
    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, 2:klev), ycoefh0(klon, 2:klev)
    REAL yzlay(klon, klev), zlev(klon, klev + 1), yteta(klon, klev)
    REAL yq2(klon, klev + 1)
    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 yt2m(klon), yq2m(klon), wind10m(klon)
    REAL ustar(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 u1(klon), v1(klon)
    REAL tair1(klon), qair1(klon), tairsol(klon)
    REAL psfce(klon), patm(klon)

    REAL qairsol(klon), zgeo1(klon)
    REAL rugo1(klon)
    REAL zgeop(klon, klev)

    !------------------------------------------------------------

    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

    ! Initialization:
    rugmer = 0.
    cdragh = 0.
    cdragm = 0.
    dflux_t = 0.
    dflux_q = 0.
    ypct = 0.
    yqsurf = 0.
    yrain_f = 0.
    ysnow_f = 0.
    yrugos = 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.
    coefh = 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(julien, 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)
             yrugos(j) = frugs(i, nsrf)
             yrugoro(j) = rugoro(i)
             yrads(j) = fsolsw(i, nsrf) + fsollw(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) yqsol(:knon) = qsol(ni(:knon))

          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 les géopotentiels de chaque couche:

          zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
               + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))

          DO k = 2, klev
             zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
                  * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
                  * (ypplay(:knon, k - 1) - ypplay(:knon, k))
          ENDDO

          CALL clcdrag(nsrf, yu(:knon, 1), yv(:knon, 1), yt(:knon, 1), &
               yq(:knon, 1), zgeop(:knon, 1), yts(:knon), yqsurf(:knon), &
               yrugos(:knon), ycdragm(:knon), ycdragh(:knon)) 

          IF (iflag_pbl == 1) THEN
             ycdragm(:knon) = max(ycdragm(:knon), 0.)
             ycdragh(:knon) = max(ycdragh(:knon), 0.)
          end IF
          
          ! on met un seuil pour ycdragm et ycdragh
          IF (nsrf == is_oce) THEN
             ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
             ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
          END IF

          CALL coefkz(nsrf, ypaprs(:knon, :), ypplay(:knon, :), ksta, &
               ksta_ter, yts(:knon), yu(:knon, :), yv(:knon, :), yt(:knon, :), &
               yq(:knon, :), zgeop(:knon, :), ycoefm(:knon, :), &
               ycoefh(:knon, :))

          IF (iflag_pbl == 1) THEN
             CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0(:knon, :), &
                  ycoefh0(:knon, :))
             ycoefm(:knon, :) = max(ycoefm(:knon, :), ycoefm0(:knon, :))
             ycoefh(:knon, :) = max(ycoefh(:knon, :), ycoefh0(:knon, :))
          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, &
                  ycdragm(:knon), ycoefh0(:knon, :))
             ycoefm0(:knon, :) = ycoefh0(:knon, :)
             ycoefm(:knon, :) = max(ycoefm(:knon, :), ycoefm0(:knon, :))
             ycoefh(:knon, :) = max(ycoefh(:knon, :), ycoefh0(:knon, :))
          END IF

          IF (iflag_pbl >= 6) 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(:knon, k) = yzlay(: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

             zlev(:knon, 1) = 0.
             zlev(:knon, klev + 1) = 2. * yzlay(:knon, klev) &
                  - yzlay(:knon, klev - 1)

             DO k = 2, klev
                zlev(:knon, k) = 0.5 * (yzlay(:knon, k) + yzlay(: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

             ustar(:knon) = ustarhb(yu(:knon, 1), yv(:knon, 1), ycdragm(:knon))
             CALL yamada4(dtime, rg, zlev(:knon, :), yzlay(:knon, :), &
                  yu(:knon, :), yv(:knon, :), yteta(:knon, :), yq2(:knon, :), &
                  ycoefm(:knon, :), ycoefh(:knon, :), ustar(:knon))
          END IF

          CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
               ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
               ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
               y_flux_u(:knon))
          CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
               ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
               ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
               y_flux_v(:knon))

          ! calculer la diffusion de "q" et de "h"
          CALL clqh(dtime, julien, firstcal, nsrf, ni(:knon), &
               ytsoil(:knon, :), yqsol(:knon), mu0, yrugos, yrugoro, &
               yu(:knon, 1), yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), &
               yt, yq, yts(:knon), ypaprs, ypplay, ydelp, yrads(:knon), &
               yalb(:knon), snow(:knon), yqsurf, yrain_f, ysnow_f, &
               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(:knon), y_dflux_q(:knon), 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 * ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2) &
                     / rg + 0.11 * 14E-6 &
                     / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**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)
          END DO

          DO k = 1, klev
             DO j = 1, knon
                i = ni(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.
          frugs(:, 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)
             frugs(i, nsrf) = yz0_new(j)
             fluxlat(i, nsrf) = yfluxlat(j)
             IF (nsrf == is_oce) THEN
                rugmer(i) = yrugm(j)
                frugs(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) + ycdragh(j) * ypct(j)
             cdragm(i) = cdragm(i) + ycdragm(j) * ypct(j)
             dflux_t(i) = dflux_t(i) + y_dflux_t(j)
             dflux_q(i) = dflux_q(i) + y_dflux_q(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)
             END DO
          END DO

          forall (k = 2:klev) coefh(ni(:knon), k) &
               = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypct(:knon)

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

          DO j = 1, knon
             i = ni(j)
             u1(j) = yu(j, 1) + y_d_u(j, 1)
             v1(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) = frugs(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, u1(:knon), v1(:knon), tair1(:knon), &
               qair1, zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, &
               yq2m, yt10m, yq10m, wind10m(:knon), ustar(:knon))

          DO j = 1, knon
             i = ni(j)
             t2m(i, nsrf) = yt2m(j)
             q2m(i, nsrf) = yq2m(j)

             u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
                  / sqrt(u1(j)**2 + v1(j)**2)
             v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
                  / sqrt(u1(j)**2 + v1(j)**2)
          END DO

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