Sources/phylmd/clmain.f

module clmain_m

  IMPLICIT NONE

contains

  SUBROUTINE clmain(dtime, pctsrf, t, q, u, v, julien, mu0, ftsol, cdmmax, &
       cdhmax, 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 coef_diff_turb_m, only: coef_diff_turb
    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
    use time_phylmdz, only: itap

    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(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 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

          IF (iflag_pbl >= 6) then
             DO k = 1, klev + 1
                DO j = 1, knon
                   i = ni(j)
                   yq2(j, k) = q2(i, k, nsrf)
                END DO
             END DO
          end IF

          call coef_diff_turb(dtime, nsrf, ni(:knon), ypaprs(:knon, :), &
               ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
               yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
               ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))

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