Sources/phylmd/clmain.f

module clmain_m

  IMPLICIT NONE

contains

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

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

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

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

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

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

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

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

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

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

    REAL, INTENT(IN):: paprs(klon, klev+1) ! pression a intercouche (Pa)
    REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
    REAL, INTENT(inout):: snow(klon, nbsrf)
    REAL qsurf(klon, nbsrf)
    REAL evap(klon, nbsrf)
    REAL, intent(inout):: falbe(klon, nbsrf)

    REAL fluxlat(klon, nbsrf)

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

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

    REAL, INTENT(IN):: solsw(klon, nbsrf), sollw(klon, nbsrf)
    REAL, intent(in):: fder(klon)
    REAL, INTENT(IN):: rlat(klon) ! latitude en degr\'es

    REAL, intent(inout):: rugos(klon, nbsrf) ! longueur de rugosit\'e (en m)

    real agesno(klon, nbsrf)
    REAL, INTENT(IN):: rugoro(klon)

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

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

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

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

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

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

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

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

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

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

    ! Local:

    LOGICAL:: firstcal = .true.

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

    REAL y_fqcalving(klon), y_ffonte(klon)
    real y_run_off_lic_0(klon)

    REAL rugmer(klon)

    REAL ytsoil(klon, nsoilmx)

    REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
    REAL yalb(klon)
    REAL yu1(klon), yv1(klon)
    ! on rajoute en output yu1 et yv1 qui sont les vents dans
    ! la premiere couche
    REAL ysnow(klon), yqsurf(klon), yagesno(klon)

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

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

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

    REAL yfder(klon)
    REAL yrugm(klon), yrads(klon), yrugoro(klon)

    REAL yfluxlat(klon)

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

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

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

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

    INTEGER ni(klon), knon, j

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

    REAL zx_alf1, zx_alf2 ! valeur ambiante par extrapolation

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

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

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

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

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

    ytherm = 0.

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

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

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

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

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

    ! Boucler sur toutes les sous-fractions du sol:

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

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

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

          DO k = 1, nsoilmx
             DO j = 1, knon
                i = ni(j)
                ytsoil(j, k) = ftsoil(i, k, nsrf)
             END DO
          END DO

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

          ftsoil(:, :, nsrf) = 0.
          DO k = 1, nsoilmx
             DO j = 1, knon
                i = ni(j)
                ftsoil(i, k, nsrf) = ytsoil(j, k)
             END DO
          END DO

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

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

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

             qairsol(j) = yqsurf(j)
          END DO

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

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

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

          END DO

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

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

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

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

    firstcal = .false.

  END SUBROUTINE clmain

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