Sources/phylmd/clmain.f

module clmain_m

  IMPLICIT NONE

contains

  SUBROUTINE clmain(dtime, pctsrf, t, q, u, v, jour, mu0, ftsol, cdmmax, &
       cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, qsol, paprs, pplay, fsnow, &
       qsurf, evap, falbe, fluxlat, rain_fall, snow_f, solsw, sollw, fder, &
       rugos, agesno, rugoro, d_t, d_q, d_u, d_v, d_ts, flux_t, flux_q, &
       flux_u, flux_v, cdragh, cdragm, q2, dflux_t, dflux_q, ycoefh, zu1, &
       zv1, t2m, q2m, u10m, v10m, pblh, capcl, oliqcl, cteicl, pblt, therm, &
       trmb1, trmb2, trmb3, plcl, fqcalving, ffonte, run_off_lic_0)

    ! From phylmd/clmain.F, version 1.6, 2005/11/16 14:47:19
    ! Author: Z. X. Li (LMD/CNRS), date: 1993/08/18
    ! Objet : interface de couche limite (diffusion verticale)

    ! Tout ce qui a trait aux traceurs est dans "phytrac". Le calcul
    ! de la couche limite pour les traceurs se fait avec "cltrac" et
    ! ne tient pas compte de la diff\'erentiation des sous-fractions
    ! de sol.

    ! Pour pouvoir extraire les coefficients d'\'echanges et le vent
    ! dans la premi\`ere couche, trois champs ont \'et\'e cr\'e\'es : "ycoefh",
    ! "zu1" et "zv1". Nous avons moyenn\'e les valeurs de ces trois
    ! champs sur les quatre sous-surfaces du mod\`ele.

    use clqh_m, only: clqh
    use clvent_m, only: clvent
    use coefkz_m, only: coefkz
    use coefkzmin_m, only: coefkzmin
    USE conf_gcm_m, ONLY: prt_level, lmt_pas
    USE conf_phys_m, ONLY: iflag_pbl
    USE dimphy, ONLY: klev, klon, zmasq
    USE dimsoil, ONLY: nsoilmx
    use hbtm_m, only: hbtm
    USE indicesol, ONLY: epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
    USE interfoce_lim_m, ONLY: interfoce_lim
    use stdlevvar_m, only: stdlevvar
    USE suphec_m, ONLY: rd, rg, rkappa
    use time_phylmdz, only: itap
    use ustarhb_m, only: ustarhb
    use vdif_kcay_m, only: vdif_kcay
    use yamada4_m, only: yamada4

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    ! Ionela Musat cf. Anne Mathieu : planetary boundary layer, hbtm
    ! (Comme les autres diagnostics on cumule dans physiq ce qui
    ! permet de sortir les grandeurs par sous-surface)
    REAL pblh(klon, nbsrf) ! height of planetary boundary layer
    REAL capcl(klon, nbsrf)
    REAL oliqcl(klon, nbsrf)
    REAL cteicl(klon, nbsrf)
    REAL pblt(klon, nbsrf)
    ! pblT------- T au nveau HCL
    REAL therm(klon, nbsrf)
    REAL trmb1(klon, nbsrf)
    ! trmb1-------deep_cape
    REAL trmb2(klon, nbsrf)
    ! trmb2--------inhibition
    REAL trmb3(klon, nbsrf)
    ! trmb3-------Point Omega
    REAL plcl(klon, nbsrf)
    REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)
    ! ffonte----Flux thermique utilise pour fondre la neige
    ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
    !           hauteur de neige, en kg/m2/s
    REAL run_off_lic_0(klon)

    ! Local:

    LOGICAL:: firstcal = .true.

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

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

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

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

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

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

    REAL yfder(klon)
    REAL yrugm(klon), yrads(klon), yrugoro(klon)
    REAL yfluxlat(klon)
    REAL y_d_ts(klon)
    REAL y_d_t(klon, klev), y_d_q(klon, klev)
    REAL y_d_u(klon, klev), y_d_v(klon, klev)
    REAL y_flux_t(klon), y_flux_q(klon)
    REAL y_flux_u(klon), y_flux_v(klon)
    REAL y_dflux_t(klon), y_dflux_q(klon)
    REAL coefh(klon, klev), coefm(klon, klev)
    REAL yu(klon, klev), yv(klon, klev)
    REAL yt(klon, klev), yq(klon, klev)
    REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)

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

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

    REAL u1lay(klon), v1lay(klon)
    REAL delp(klon, klev)
    INTEGER i, k, nsrf

    INTEGER ni(klon), knon, j

    REAL pctsrf_pot(klon, nbsrf)
    ! "pourcentage potentiel" pour tenir compte des \'eventuelles
    ! apparitions ou disparitions de la glace de mer

    REAL zx_alf1, zx_alf2 ! valeur ambiante par extrapolation

    REAL yt2m(klon), yq2m(klon), yu10m(klon)
    REAL yustar(klon)

    REAL yt10m(klon), yq10m(klon)
    REAL ypblh(klon)
    REAL ylcl(klon)
    REAL ycapcl(klon)
    REAL yoliqcl(klon)
    REAL ycteicl(klon)
    REAL ypblt(klon)
    REAL ytherm(klon)
    REAL ytrmb1(klon)
    REAL ytrmb2(klon)
    REAL ytrmb3(klon)
    REAL uzon(klon), vmer(klon)
    REAL tair1(klon), qair1(klon), tairsol(klon)
    REAL psfce(klon), patm(klon)

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

    ! utiliser un jeu de fonctions simples               
    LOGICAL zxli
    PARAMETER (zxli=.FALSE.)

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

    ytherm = 0.

    DO k = 1, klev ! epaisseur de couche
       DO i = 1, klon
          delp(i, k) = paprs(i, k) - paprs(i, k+1)
       END DO
    END DO
    DO i = 1, klon ! vent de la premiere couche
       zx_alf1 = 1.0
       zx_alf2 = 1.0 - zx_alf1
       u1lay(i) = u(i, 1)*zx_alf1 + u(i, 2)*zx_alf2
       v1lay(i) = v(i, 1)*zx_alf1 + v(i, 2)*zx_alf2
    END DO

    ! Initialization:
    rugmer = 0.
    cdragh = 0.
    cdragm = 0.
    dflux_t = 0.
    dflux_q = 0.
    zu1 = 0.
    zv1 = 0.
    ypct = 0.
    yts = 0.
    yqsurf = 0.
    yrain_f = 0.
    ysnow_f = 0.
    yfder = 0.
    yrugos = 0.
    yu1 = 0.
    yv1 = 0.
    yrads = 0.
    ypaprs = 0.
    ypplay = 0.
    ydelp = 0.
    yu = 0.
    yv = 0.
    yt = 0.
    yq = 0.
    y_dflux_t = 0.
    y_dflux_q = 0.
    yrugoro = 0.
    d_ts = 0.
    flux_t = 0.
    flux_q = 0.
    flux_u = 0.
    flux_v = 0.
    fluxlat = 0.
    d_t = 0.
    d_q = 0.
    d_u = 0.
    d_v = 0.
    ycoefh = 0.

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

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

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

    ! Boucler sur toutes les sous-fractions du sol:

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

       if_knon: IF (knon /= 0) then
          DO j = 1, knon
             i = ni(j)
             ypct(j) = pctsrf(i, nsrf)
             yts(j) = ftsol(i, nsrf)
             snow(j) = fsnow(i, nsrf)
             yqsurf(j) = qsurf(i, nsrf)
             yalb(j) = falbe(i, nsrf)
             yrain_f(j) = rain_fall(i)
             ysnow_f(j) = snow_f(i)
             yagesno(j) = agesno(i, nsrf)
             yfder(j) = fder(i)
             yrugos(j) = rugos(i, nsrf)
             yrugoro(j) = rugoro(i)
             yu1(j) = u1lay(i)
             yv1(j) = v1lay(i)
             yrads(j) = solsw(i, nsrf) + sollw(i, nsrf)
             ypaprs(j, klev+1) = paprs(i, klev+1)
             y_run_off_lic_0(j) = run_off_lic_0(i)
          END DO

          ! For continent, copy soil water content
          IF (nsrf == is_ter) THEN
             yqsol(:knon) = qsol(ni(:knon))
          ELSE
             yqsol = 0.
          END IF

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

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

          ! calculer Cdrag et les coefficients d'echange
          CALL coefkz(nsrf, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, yu, &
               yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))
          IF (iflag_pbl == 1) THEN
             CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
             coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
             coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
          END IF

          ! on met un seuil pour coefm et coefh
          IF (nsrf == is_oce) THEN
             coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)
             coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)
          END IF

          IF (ok_kzmin) THEN
             ! Calcul d'une diffusion minimale pour les conditions tres stables
             CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
                  coefm(:knon, 1), ycoefm0, ycoefh0)
             coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
             coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
          END IF

          IF (iflag_pbl >= 3) THEN
             ! Mellor et Yamada adapt\'e \`a Mars, Richard Fournier et
             ! Fr\'ed\'eric Hourdin
             yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
                  + ypplay(:knon, 1))) &
                  * (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg
             DO k = 2, klev
                yzlay(1:knon, k) = yzlay(1:knon, k-1) &
                     + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
                     / ypaprs(1:knon, k) &
                     * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
             END DO
             DO k = 1, klev
                yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
                     / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))
             END DO
             yzlev(1:knon, 1) = 0.
             yzlev(:knon, klev+1) = 2. * yzlay(:knon, klev) &
                  - yzlay(:knon, klev - 1)
             DO k = 2, klev
                yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
             END DO
             DO k = 1, klev + 1
                DO j = 1, knon
                   i = ni(j)
                   yq2(j, k) = q2(i, k, nsrf)
                END DO
             END DO

             CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)
             IF (prt_level > 9) PRINT *, 'USTAR = ', yustar

             ! iflag_pbl peut \^etre utilis\'e comme longueur de m\'elange

             IF (iflag_pbl >= 11) THEN
                CALL vdif_kcay(knon, dtime, rg, ypaprs, yzlev, yzlay, yu, yv, &
                     yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, yustar, &
                     iflag_pbl)
             ELSE
                CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &
                     coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
             END IF

             coefm(:knon, 2:) = ykmm(:knon, 2:klev)
             coefh(:knon, 2:) = ykmn(:knon, 2:klev)
          END IF

          ! calculer la diffusion des vitesses "u" et "v"
          CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &
               ypplay, ydelp, y_d_u, y_flux_u(:knon))
          CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &
               ypplay, ydelp, y_d_v, y_flux_v(:knon))

          ! calculer la diffusion de "q" et de "h"
          CALL clqh(dtime, jour, firstcal, nsrf, ni(:knon), ytsoil(:knon, :), &
               yqsol, mu0, yrugos, yrugoro, yu1, yv1, coefh(:knon, :), yt, &
               yq, yts(:knon), ypaprs, ypplay, ydelp, yrads, yalb(:knon), &
               snow(:knon), yqsurf, yrain_f, ysnow_f, yfder, yfluxlat(:knon), &
               pctsrf_new_sic, yagesno(:knon), y_d_t, y_d_q, y_d_ts(:knon), &
               yz0_new, y_flux_t(:knon), y_flux_q(:knon), y_dflux_t, &
               y_dflux_q, y_fqcalving, y_ffonte, y_run_off_lic_0)

          ! calculer la longueur de rugosite sur ocean
          yrugm = 0.
          IF (nsrf == is_oce) THEN
             DO j = 1, knon
                yrugm(j) = 0.018*coefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &
                     0.11*14E-6/sqrt(coefm(j, 1)*(yu1(j)**2+yv1(j)**2))
                yrugm(j) = max(1.5E-05, yrugm(j))
             END DO
          END IF
          DO j = 1, knon
             y_dflux_t(j) = y_dflux_t(j)*ypct(j)
             y_dflux_q(j) = y_dflux_q(j)*ypct(j)
             yu1(j) = yu1(j)*ypct(j)
             yv1(j) = yv1(j)*ypct(j)
          END DO

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

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

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

          falbe(:, nsrf) = 0.
          fsnow(:, nsrf) = 0.
          qsurf(:, nsrf) = 0.
          rugos(:, nsrf) = 0.
          DO j = 1, knon
             i = ni(j)
             d_ts(i, nsrf) = y_d_ts(j)
             falbe(i, nsrf) = yalb(j)
             fsnow(i, nsrf) = snow(j)
             qsurf(i, nsrf) = yqsurf(j)
             rugos(i, nsrf) = yz0_new(j)
             fluxlat(i, nsrf) = yfluxlat(j)
             IF (nsrf == is_oce) THEN
                rugmer(i) = yrugm(j)
                rugos(i, nsrf) = yrugm(j)
             END IF
             agesno(i, nsrf) = yagesno(j)
             fqcalving(i, nsrf) = y_fqcalving(j)
             ffonte(i, nsrf) = y_ffonte(j)
             cdragh(i) = cdragh(i) + coefh(j, 1)
             cdragm(i) = cdragm(i) + coefm(j, 1)
             dflux_t(i) = dflux_t(i) + y_dflux_t(j)
             dflux_q(i) = dflux_q(i) + y_dflux_q(j)
             zu1(i) = zu1(i) + yu1(j)
             zv1(i) = zv1(i) + yv1(j)
          END DO
          IF (nsrf == is_ter) THEN
             qsol(ni(:knon)) = yqsol(:knon)
          else IF (nsrf == is_lic) THEN
             DO j = 1, knon
                i = ni(j)
                run_off_lic_0(i) = y_run_off_lic_0(j)
             END DO
          END IF

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

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

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

          DO j = 1, knon
             i = ni(j)
             uzon(j) = yu(j, 1) + y_d_u(j, 1)
             vmer(j) = yv(j, 1) + y_d_v(j, 1)
             tair1(j) = yt(j, 1) + y_d_t(j, 1)
             qair1(j) = yq(j, 1) + y_d_q(j, 1)
             zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &
                  1)))*(ypaprs(j, 1)-ypplay(j, 1))
             tairsol(j) = yts(j) + y_d_ts(j)
             rugo1(j) = yrugos(j)
             IF (nsrf == is_oce) THEN
                rugo1(j) = rugos(i, nsrf)
             END IF
             psfce(j) = ypaprs(j, 1)
             patm(j) = ypplay(j, 1)

             qairsol(j) = yqsurf(j)
          END DO

          CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, &
               zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &
               yt10m, yq10m, yu10m, yustar)

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

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

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

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

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

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

    firstcal = .false.

  END SUBROUTINE clmain

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