Sources/phylmd/clmain.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, 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):: 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 à 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, INTENT(inout):: 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, 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), yrugos(klon), ypct(klon), yz0_new(klon)
    REAL yalb(klon)

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

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

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

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

    REAL 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.
    yqsurf = 0.
    yrain_f = 0.
    ysnow_f = 0.
    yrugos = 0.
    yu1 = 0.
    yv1 = 0.
    ypaprs = 0.
    ypplay = 0.
    ydelp = 0.
    yu = 0.
    yv = 0.
    yt = 0.
    yq = 0.
    y_dflux_t = 0.
    y_dflux_q = 0.
    yrugoro = 0.
    d_ts = 0.
    flux_t = 0.
    flux_q = 0.
    flux_u = 0.
    flux_v = 0.
    fluxlat = 0.
    d_t = 0.
    d_q = 0.
    d_u = 0.
    d_v = 0.
    ycoefh = 0.

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

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

    ! Tester si c'est le moment de lire le fichier:
    if (mod(itap - 1, lmt_pas) == 0) then
       CALL interfoce_lim(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)
             yu1(j) = u1lay(i)
             yv1(j) = v1lay(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) THEN
             yqsol(:knon) = qsol(ni(:knon))
          ELSE
             yqsol = 0.
          END IF

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

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

          ! calculer Cdrag et les coefficients d'echange
          CALL coefkz(nsrf, ypaprs, ypplay, ksta, ksta_ter, yts(:knon), &
               yrugos, yu, yv, yt, yq, yqsurf(:knon), 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, julien, firstcal, nsrf, ni(:knon), &
               ytsoil(:knon, :), yqsol, mu0, yrugos, yrugoro, yu1, yv1, &
               coefh(: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*coefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &
                     0.11*14E-6/sqrt(coefm(j, 1)*(yu1(j)**2+yv1(j)**2))
                yrugm(j) = max(1.5E-05, yrugm(j))
             END DO
          END IF
          DO j = 1, knon
             y_dflux_t(j) = y_dflux_t(j)*ypct(j)
             y_dflux_q(j) = y_dflux_q(j)*ypct(j)
             yu1(j) = yu1(j)*ypct(j)
             yv1(j) = yv1(j)*ypct(j)
          END DO

          DO k = 1, klev
             DO j = 1, knon
                i = ni(j)
                coefh(j, k) = coefh(j, k)*ypct(j)
                coefm(j, k) = coefm(j, k)*ypct(j)
                y_d_t(j, k) = y_d_t(j, k)*ypct(j)
                y_d_q(j, k) = y_d_q(j, k)*ypct(j)
                y_d_u(j, k) = y_d_u(j, k)*ypct(j)
                y_d_v(j, k) = y_d_v(j, k)*ypct(j)
             END DO
          END DO

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

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

          falbe(:, nsrf) = 0.
          fsnow(:, nsrf) = 0.
          qsurf(:, nsrf) = 0.
          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) + coefh(j, 1)
             cdragm(i) = cdragm(i) + coefm(j, 1)
             dflux_t(i) = dflux_t(i) + y_dflux_t(j)
             dflux_q(i) = dflux_q(i) + y_dflux_q(j)
             zu1(i) = zu1(i) + yu1(j)
             zv1(i) = zv1(i) + yv1(j)
          END DO
          IF (nsrf == is_ter) THEN
             qsol(ni(:knon)) = yqsol(:knon)
          else IF (nsrf == is_lic) THEN
             DO j = 1, knon
                i = ni(j)
                run_off_lic_0(i) = y_run_off_lic_0(j)
             END DO
          END IF

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

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

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

          DO j = 1, knon
             i = ni(j)
             uzon(j) = yu(j, 1) + y_d_u(j, 1)
             vmer(j) = yv(j, 1) + y_d_v(j, 1)
             tair1(j) = yt(j, 1) + y_d_t(j, 1)
             qair1(j) = yq(j, 1) + y_d_q(j, 1)
             zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &
                  1)))*(ypaprs(j, 1)-ypplay(j, 1))
             tairsol(j) = yts(j) + y_d_ts(j)
             rugo1(j) = yrugos(j)
             IF (nsrf == is_oce) THEN
                rugo1(j) = 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, zxli, uzon, vmer, tair1, qair1, &
               zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &
               yt10m, yq10m, yu10m, yustar)

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

             ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
             u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)
             v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)
          END DO

          CALL hbtm(ypaprs, ypplay, yt2m, yq2m, yustar, y_flux_t(:knon), &
               y_flux_q(:knon), yu, yv, yt, yq, ypblh(:knon), ycapcl, &
               yoliqcl, ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)

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

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

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