Sources/phylmd/clmain.f

module clmain_m

  IMPLICIT NONE

contains

  SUBROUTINE clmain(dtime, itap, pctsrf, pctsrf_new, t, q, u, v, jour, rmu0, &
       co2_ppm, ts, soil_model, cdmmax, cdhmax, ksta, ksta_ter, &
       ok_kzmin, ftsoil, qsol, paprs, pplay, snow, qsurf, evap, albe, alblw, &
       fluxlat, rain_fall, snow_f, solsw, sollw, fder, rlat, rugos, debut, &
       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, flux_o, flux_g, &
       tslab, seaice)

    ! 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érentiation des sous-fractions de
    ! sol.

    ! Pour pouvoir extraire les coefficients d'échanges et le vent
    ! dans la première couche, trois champs ont été créés : "ycoefh",
    ! "zu1" et "zv1". Nous avons moyenné les valeurs de ces trois
    ! champs sur les quatre sous-surfaces du modèle.

    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
    USE conf_phys_m, ONLY: iflag_pbl
    USE dimens_m, ONLY: iim, jjm
    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 suphec_m, ONLY: rd, rg, rkappa
    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)
    INTEGER, INTENT(IN):: itap ! numero du pas de temps
    REAL, INTENT(inout):: pctsrf(klon, nbsrf)

    ! la nouvelle repartition des surfaces sortie de l'interface
    REAL, INTENT(out):: pctsrf_new(klon, nbsrf)

    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):: co2_ppm ! taux CO2 atmosphere
    REAL, INTENT(IN):: ts(klon, nbsrf) ! input-R- temperature du sol (en Kelvin)
    LOGICAL, INTENT(IN):: soil_model
    REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh
    REAL, INTENT(IN):: ksta, ksta_ter
    LOGICAL, INTENT(IN):: ok_kzmin
    REAL ftsoil(klon, nsoilmx, nbsrf)
    REAL, INTENT(inout):: qsol(klon)
    REAL, INTENT(IN):: paprs(klon, klev+1) ! pression a intercouche (Pa)
    REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
    REAL snow(klon, nbsrf)
    REAL qsurf(klon, nbsrf)
    REAL evap(klon, nbsrf)
    REAL albe(klon, nbsrf)
    REAL alblw(klon, nbsrf)

    REAL fluxlat(klon, nbsrf)

    REAL, intent(in):: rain_fall(klon), snow_f(klon)
    REAL, INTENT(IN):: solsw(klon, nbsrf), sollw(klon, nbsrf)
    REAL fder(klon)
    REAL, INTENT(IN):: rlat(klon) ! latitude en degrés

    REAL rugos(klon, nbsrf)
    ! rugos----input-R- longeur de rugosite (en m)

    LOGICAL, INTENT(IN):: debut
    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 d_ts(klon, nbsrf)
    ! d_ts-----output-R- le changement pour "ts"

    REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)
    ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
    !                    (orientation positive vers le bas)
    ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)

    REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)
    ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
    ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal

    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)

    !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds
    ! physiq ce qui permet de sortir les grdeurs par sous surface)
    REAL pblh(klon, nbsrf)
    ! pblh------- HCL
    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)

    REAL flux_o(klon), flux_g(klon)
    !IM "slab" ocean
    ! flux_g---output-R-  flux glace (pour OCEAN='slab  ')
    ! flux_o---output-R-  flux ocean (pour OCEAN='slab  ')

    REAL tslab(klon)
    ! tslab-in/output-R temperature du slab ocean (en Kelvin) 
    ! uniqmnt pour slab

    REAL seaice(klon)
    ! seaice---output-R-  glace de mer (kg/m2) (pour OCEAN='slab  ')

    ! Local:

    REAL y_flux_o(klon), y_flux_g(klon)
    real ytslab(klon)
    real y_seaice(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 yalblw(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), yqsol(klon)
    REAL yrain_f(klon), ysnow_f(klon)
    REAL ysollw(klon), ysolsw(klon)
    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, klev), y_flux_q(klon, klev)
    REAL y_flux_u(klon, klev), y_flux_v(klon, klev)
    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 éventuelles
    ! apparitions ou disparitions de la glace de mer

    REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.

    REAL yt2m(klon), yq2m(klon), yu10m(klon)
    REAL yustar(klon)
    ! -- LOOP
    REAL yu10mx(klon)
    REAL yu10my(klon)
    REAL ywindsp(klon)
    ! -- LOOP

    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.
    yalb = 0.
    yalblw = 0.
    yrain_f = 0.
    ysnow_f = 0.
    yfder = 0.
    ysolsw = 0.
    ysollw = 0.
    yrugos = 0.
    yu1 = 0.
    yv1 = 0.
    yrads = 0.
    ypaprs = 0.
    ypplay = 0.
    ydelp = 0.
    yu = 0.
    yv = 0.
    yt = 0.
    yq = 0.
    pctsrf_new = 0.
    y_flux_u = 0.
    y_flux_v = 0.
    y_dflux_t = 0.
    y_dflux_q = 0.
    ytsoil = 999999.
    yrugoro = 0.
    ! -- LOOP
    yu10mx = 0.
    yu10my = 0.
    ywindsp = 0.
    ! -- LOOP
    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ère ici qu'on
    ! peut avoir potentiellement de la glace sur tout le domaine océanique
    ! (à affiner)

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

    ! 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éterminer le domaine à 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)
             ytslab(i) = tslab(i)
             ysnow(j) = snow(i, nsrf)
             yqsurf(j) = qsurf(i, nsrf)
             yalb(j) = albe(i, nsrf)
             yalblw(j) = alblw(i, nsrf)
             yrain_f(j) = rain_fall(i)
             ysnow_f(j) = snow_f(i)
             yagesno(j) = agesno(i, nsrf)
             yfder(j) = fder(i)
             ysolsw(j) = solsw(i, nsrf)
             ysollw(j) = sollw(i, nsrf)
             yrugos(j) = rugos(i, nsrf)
             yrugoro(j) = rugoro(i)
             yu1(j) = u1lay(i)
             yv1(j) = v1lay(i)
             yrads(j) = ysolsw(j) + ysollw(j)
             ypaprs(j, klev+1) = paprs(i, klev+1)
             y_run_off_lic_0(j) = run_off_lic_0(i)
             yu10mx(j) = u10m(i, nsrf)
             yu10my(j) = v10m(i, nsrf)
             ywindsp(j) = sqrt(yu10mx(j)*yu10mx(j)+yu10my(j)*yu10my(j))
          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é à Mars, Richard Fournier et
             ! Frédéric 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 être utilisé comme longueur de mélange

             IF (iflag_pbl >= 11) THEN
                CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, 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)
          CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &
               ypplay, ydelp, y_d_v, y_flux_v)

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

          ! 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)
                flux_t(i, k, nsrf) = y_flux_t(j, k)
                flux_q(i, k, nsrf) = y_flux_q(j, k)
                flux_u(i, k, nsrf) = y_flux_u(j, k)
                flux_v(i, k, nsrf) = y_flux_v(j, k)
                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

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

          albe(:, nsrf) = 0.
          alblw(:, 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)
             albe(i, nsrf) = yalb(j)
             alblw(i, nsrf) = yalblw(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
          !$$$ PB ajout pour soil
          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(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, &
               y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, 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
          !IM "slab" ocean
          IF (nsrf == is_oce) THEN
             DO j = 1, knon
                ! on projette sur la grille globale
                i = ni(j)
                IF (pctsrf_new(i, is_oce)>epsfra) THEN
                   flux_o(i) = y_flux_o(j)
                ELSE
                   flux_o(i) = 0.
                END IF
             END DO
          END IF

          IF (nsrf == is_sic) THEN
             DO j = 1, knon
                i = ni(j)
                ! On pondère lorsque l'on fait le bilan au sol :
                IF (pctsrf_new(i, is_sic)>epsfra) THEN
                   flux_g(i) = y_flux_g(j)
                ELSE
                   flux_g(i) = 0.
                END IF
             END DO

          END IF
       end IF if_knon
    END DO loop_surface

    ! On utilise les nouvelles surfaces

    rugos(:, is_oce) = rugmer
    pctsrf = pctsrf_new

  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	99	SUBROUTINE clmain(dtime, itap, pctsrf, pctsrf_new, t, q, u, v, jour, rmu0, &
8			co2_ppm, ts, soil_model, cdmmax, cdhmax, ksta, ksta_ter, &
9			ok_kzmin, ftsoil, qsol, paprs, pplay, snow, qsurf, evap, albe, alblw, &
10			fluxlat, rain_fall, snow_f, solsw, sollw, fder, rlat, rugos, debut, &
11			agesno, rugoro, d_t, d_q, d_u, d_v, d_ts, flux_t, flux_q, flux_u, &
12			flux_v, cdragh, cdragm, q2, dflux_t, dflux_q, ycoefh, zu1, zv1, t2m, &
13			q2m, u10m, v10m, pblh, capcl, oliqcl, cteicl, pblt, therm, trmb1, &
14			trmb2, trmb3, plcl, fqcalving, ffonte, run_off_lic_0, flux_o, flux_g, &
15			tslab, seaice)
16	guez	3
17	guez	99	! From phylmd/clmain.F, version 1.6, 2005/11/16 14:47:19
18	guez	62	! Author: Z. X. Li (LMD/CNRS), date: 1993/08/18
19			! Objet : interface de couche limite (diffusion verticale)
20	guez	3
21	guez	62	! Tout ce qui a trait aux traceurs est dans "phytrac". Le calcul
22			! de la couche limite pour les traceurs se fait avec "cltrac" et
23			! ne tient pas compte de la différentiation des sous-fractions de
24			! sol.
25	guez	3
26	guez	38	! Pour pouvoir extraire les coefficients d'échanges et le vent
27	guez	70	! dans la première couche, trois champs ont été créés : "ycoefh",
28	guez	62	! "zu1" et "zv1". Nous avons moyenné les valeurs de ces trois
29			! champs sur les quatre sous-surfaces du modèle.
30	guez	3
31	guez	49	use clqh_m, only: clqh
32	guez	62	use clvent_m, only: clvent
33	guez	47	use coefkz_m, only: coefkz
34			use coefkzmin_m, only: coefkzmin
35	guez	62	USE conf_gcm_m, ONLY: prt_level
36			USE conf_phys_m, ONLY: iflag_pbl
37			USE dimens_m, ONLY: iim, jjm
38			USE dimphy, ONLY: klev, klon, zmasq
39			USE dimsoil, ONLY: nsoilmx
40	guez	47	use hbtm_m, only: hbtm
41	guez	62	USE indicesol, ONLY: epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
42			USE suphec_m, ONLY: rd, rg, rkappa
43			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			INTEGER, INTENT(IN):: itap ! numero du pas de temps
49			REAL, INTENT(inout):: pctsrf(klon, nbsrf)
50
51			! la nouvelle repartition des surfaces sortie de l'interface
52			REAL, INTENT(out):: pctsrf_new(klon, nbsrf)
53
54			REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
55			REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg/kg)
56			REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
57			INTEGER, INTENT(IN):: jour ! jour de l'annee en cours
58			REAL, intent(in):: rmu0(klon) ! cosinus de l'angle solaire zenithal
59	guez	99	REAL, intent(in):: co2_ppm ! taux CO2 atmosphere
60			REAL, INTENT(IN):: ts(klon, nbsrf) ! input-R- temperature du sol (en Kelvin)
61	guez	70	LOGICAL, INTENT(IN):: soil_model
62	guez	71	REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh
63	guez	99	REAL, INTENT(IN):: ksta, ksta_ter
64			LOGICAL, INTENT(IN):: ok_kzmin
65	guez	70	REAL ftsoil(klon, nsoilmx, nbsrf)
66	guez	99	REAL, INTENT(inout):: qsol(klon)
67	guez	62	REAL, INTENT(IN):: paprs(klon, klev+1) ! pression a intercouche (Pa)
68			REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
69	guez	70	REAL snow(klon, nbsrf)
70			REAL qsurf(klon, nbsrf)
71			REAL evap(klon, nbsrf)
72			REAL albe(klon, nbsrf)
73			REAL alblw(klon, nbsrf)
74
75			REAL fluxlat(klon, nbsrf)
76
77			REAL, intent(in):: rain_fall(klon), snow_f(klon)
78	guez	72	REAL, INTENT(IN):: solsw(klon, nbsrf), sollw(klon, nbsrf)
79	guez	70	REAL fder(klon)
80	guez	62	REAL, INTENT(IN):: rlat(klon) ! latitude en degrés
81	guez	70
82			REAL rugos(klon, nbsrf)
83			! rugos----input-R- longeur de rugosite (en m)
84
85			LOGICAL, INTENT(IN):: debut
86			real agesno(klon, nbsrf)
87			REAL, INTENT(IN):: rugoro(klon)
88
89	guez	38	REAL d_t(klon, klev), d_q(klon, klev)
90	guez	49	! d_t------output-R- le changement pour "t"
91			! d_q------output-R- le changement pour "q"
92	guez	62
93			REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
94			! changement pour "u" et "v"
95
96	guez	70	REAL d_ts(klon, nbsrf)
97			! d_ts-----output-R- le changement pour "ts"
98
99	guez	38	REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)
100	guez	49	! flux_t---output-R- flux de chaleur sensible (CpT) J/m2/s (W/m2)
101			! (orientation positive vers le bas)
102			! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
103	guez	70
104			REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)
105			! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
106			! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
107
108			REAL, INTENT(out):: cdragh(klon), cdragm(klon)
109			real q2(klon, klev+1, nbsrf)
110
111	guez	99	REAL, INTENT(out):: dflux_t(klon), dflux_q(klon)
112	guez	49	! dflux_t derive du flux sensible
113			! dflux_q derive du flux latent
114	guez	38	!IM "slab" ocean
115	guez	70
116			REAL, intent(out):: ycoefh(klon, klev)
117			REAL, intent(out):: zu1(klon)
118			REAL zv1(klon)
119			REAL t2m(klon, nbsrf), q2m(klon, nbsrf)
120			REAL u10m(klon, nbsrf), v10m(klon, nbsrf)
121
122			!IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds
123			! physiq ce qui permet de sortir les grdeurs par sous surface)
124			REAL pblh(klon, nbsrf)
125			! pblh------- HCL
126			REAL capcl(klon, nbsrf)
127			REAL oliqcl(klon, nbsrf)
128			REAL cteicl(klon, nbsrf)
129			REAL pblt(klon, nbsrf)
130			! pblT------- T au nveau HCL
131			REAL therm(klon, nbsrf)
132			REAL trmb1(klon, nbsrf)
133			! trmb1-------deep_cape
134			REAL trmb2(klon, nbsrf)
135			! trmb2--------inhibition
136			REAL trmb3(klon, nbsrf)
137			! trmb3-------Point Omega
138			REAL plcl(klon, nbsrf)
139			REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)
140			! ffonte----Flux thermique utilise pour fondre la neige
141			! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
142			! hauteur de neige, en kg/m2/s
143			REAL run_off_lic_0(klon)
144
145	guez	38	REAL flux_o(klon), flux_g(klon)
146	guez	49	!IM "slab" ocean
147			! flux_g---output-R- flux glace (pour OCEAN='slab ')
148			! flux_o---output-R- flux ocean (pour OCEAN='slab ')
149	guez	70
150			REAL tslab(klon)
151	guez	49	! tslab-in/output-R temperature du slab ocean (en Kelvin)
152			! uniqmnt pour slab
153	guez	70
154			REAL seaice(klon)
155	guez	49	! seaice---output-R- glace de mer (kg/m2) (pour OCEAN='slab ')
156	guez	15
157	guez	70	! Local:
158	guez	15
159	guez	70	REAL y_flux_o(klon), y_flux_g(klon)
160			real ytslab(klon)
161			real y_seaice(klon)
162			REAL y_fqcalving(klon), y_ffonte(klon)
163			real y_run_off_lic_0(klon)
164	guez	15
165	guez	70	REAL rugmer(klon)
166	guez	15
167	guez	38	REAL ytsoil(klon, nsoilmx)
168	guez	15
169	guez	38	REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
170			REAL yalb(klon)
171			REAL yalblw(klon)
172			REAL yu1(klon), yv1(klon)
173	guez	49	! on rajoute en output yu1 et yv1 qui sont les vents dans
174			! la premiere couche
175	guez	38	REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)
176			REAL yrain_f(klon), ysnow_f(klon)
177	guez	72	REAL ysollw(klon), ysolsw(klon)
178	guez	99	REAL yfder(klon)
179	guez	38	REAL yrugm(klon), yrads(klon), yrugoro(klon)
180	guez	15
181	guez	38	REAL yfluxlat(klon)
182	guez	15
183	guez	38	REAL y_d_ts(klon)
184			REAL y_d_t(klon, klev), y_d_q(klon, klev)
185			REAL y_d_u(klon, klev), y_d_v(klon, klev)
186			REAL y_flux_t(klon, klev), y_flux_q(klon, klev)
187			REAL y_flux_u(klon, klev), y_flux_v(klon, klev)
188			REAL y_dflux_t(klon), y_dflux_q(klon)
189	guez	62	REAL coefh(klon, klev), coefm(klon, klev)
190	guez	38	REAL yu(klon, klev), yv(klon, klev)
191			REAL yt(klon, klev), yq(klon, klev)
192			REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)
193	guez	15
194	guez	38	REAL ycoefm0(klon, klev), ycoefh0(klon, klev)
195	guez	15
196	guez	38	REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)
197			REAL ykmm(klon, klev+1), ykmn(klon, klev+1)
198			REAL ykmq(klon, klev+1)
199	guez	70	REAL yq2(klon, klev+1)
200	guez	38	REAL q2diag(klon, klev+1)
201	guez	15
202	guez	38	REAL u1lay(klon), v1lay(klon)
203			REAL delp(klon, klev)
204			INTEGER i, k, nsrf
205	guez	15
206	guez	38	INTEGER ni(klon), knon, j
207	guez	40
208	guez	38	REAL pctsrf_pot(klon, nbsrf)
209	guez	40	! "pourcentage potentiel" pour tenir compte des éventuelles
210			! apparitions ou disparitions de la glace de mer
211	guez	15
212	guez	38	REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
213	guez	15
214	guez	38	REAL yt2m(klon), yq2m(klon), yu10m(klon)
215			REAL yustar(klon)
216			! -- LOOP
217			REAL yu10mx(klon)
218			REAL yu10my(klon)
219			REAL ywindsp(klon)
220			! -- LOOP
221	guez	15
222	guez	38	REAL yt10m(klon), yq10m(klon)
223			REAL ypblh(klon)
224			REAL ylcl(klon)
225			REAL ycapcl(klon)
226			REAL yoliqcl(klon)
227			REAL ycteicl(klon)
228			REAL ypblt(klon)
229			REAL ytherm(klon)
230			REAL ytrmb1(klon)
231			REAL ytrmb2(klon)
232			REAL ytrmb3(klon)
233			REAL uzon(klon), vmer(klon)
234			REAL tair1(klon), qair1(klon), tairsol(klon)
235			REAL psfce(klon), patm(klon)
236	guez	15
237	guez	38	REAL qairsol(klon), zgeo1(klon)
238			REAL rugo1(klon)
239	guez	15
240	guez	38	! utiliser un jeu de fonctions simples
241			LOGICAL zxli
242			PARAMETER (zxli=.FALSE.)
243	guez	15
244	guez	38	!------------------------------------------------------------
245	guez	15
246	guez	38	ytherm = 0.
247	guez	15
248	guez	38	DO k = 1, klev ! epaisseur de couche
249			DO i = 1, klon
250			delp(i, k) = paprs(i, k) - paprs(i, k+1)
251			END DO
252			END DO
253			DO i = 1, klon ! vent de la premiere couche
254			zx_alf1 = 1.0
255			zx_alf2 = 1.0 - zx_alf1
256			u1lay(i) = u(i, 1)zx_alf1 + u(i, 2)zx_alf2
257			v1lay(i) = v(i, 1)zx_alf1 + v(i, 2)zx_alf2
258			END DO
259	guez	15
260	guez	40	! Initialization:
261			rugmer = 0.
262			cdragh = 0.
263			cdragm = 0.
264			dflux_t = 0.
265			dflux_q = 0.
266			zu1 = 0.
267			zv1 = 0.
268			ypct = 0.
269			yts = 0.
270			ysnow = 0.
271			yqsurf = 0.
272			yalb = 0.
273			yalblw = 0.
274			yrain_f = 0.
275			ysnow_f = 0.
276			yfder = 0.
277			ysolsw = 0.
278			ysollw = 0.
279			yrugos = 0.
280			yu1 = 0.
281			yv1 = 0.
282			yrads = 0.
283			ypaprs = 0.
284			ypplay = 0.
285			ydelp = 0.
286			yu = 0.
287			yv = 0.
288			yt = 0.
289			yq = 0.
290			pctsrf_new = 0.
291			y_flux_u = 0.
292			y_flux_v = 0.
293			y_dflux_t = 0.
294			y_dflux_q = 0.
295	guez	38	ytsoil = 999999.
296			yrugoro = 0.
297			! -- LOOP
298	guez	40	yu10mx = 0.
299			yu10my = 0.
300			ywindsp = 0.
301	guez	38	! -- LOOP
302	guez	40	d_ts = 0.
303	guez	38	yfluxlat = 0.
304			flux_t = 0.
305			flux_q = 0.
306			flux_u = 0.
307			flux_v = 0.
308	guez	40	d_t = 0.
309			d_q = 0.
310			d_u = 0.
311			d_v = 0.
312	guez	70	ycoefh = 0.
313	guez	15
314	guez	38	! Initialisation des "pourcentages potentiels". On considère ici qu'on
315			! peut avoir potentiellement de la glace sur tout le domaine océanique
316			! (à affiner)
317	guez	15
318	guez	38	pctsrf_pot = pctsrf
319			pctsrf_pot(:, is_oce) = 1. - zmasq
320			pctsrf_pot(:, is_sic) = 1. - zmasq
321	guez	15
322	guez	99	! Boucler sur toutes les sous-fractions du sol:
323
324	guez	49	loop_surface: DO nsrf = 1, nbsrf
325			! Chercher les indices :
326	guez	38	ni = 0
327			knon = 0
328			DO i = 1, klon
329	guez	40	! Pour déterminer le domaine à traiter, on utilise les surfaces
330	guez	38	! "potentielles"
331			IF (pctsrf_pot(i, nsrf) > epsfra) THEN
332			knon = knon + 1
333			ni(knon) = i
334			END IF
335			END DO
336	guez	15
337	guez	62	if_knon: IF (knon /= 0) then
338	guez	38	DO j = 1, knon
339			i = ni(j)
340	guez	62	ypct(j) = pctsrf(i, nsrf)
341			yts(j) = ts(i, nsrf)
342			ytslab(i) = tslab(i)
343			ysnow(j) = snow(i, nsrf)
344			yqsurf(j) = qsurf(i, nsrf)
345			yalb(j) = albe(i, nsrf)
346			yalblw(j) = alblw(i, nsrf)
347			yrain_f(j) = rain_fall(i)
348			ysnow_f(j) = snow_f(i)
349			yagesno(j) = agesno(i, nsrf)
350			yfder(j) = fder(i)
351			ysolsw(j) = solsw(i, nsrf)
352			ysollw(j) = sollw(i, nsrf)
353			yrugos(j) = rugos(i, nsrf)
354			yrugoro(j) = rugoro(i)
355			yu1(j) = u1lay(i)
356			yv1(j) = v1lay(i)
357			yrads(j) = ysolsw(j) + ysollw(j)
358			ypaprs(j, klev+1) = paprs(i, klev+1)
359			y_run_off_lic_0(j) = run_off_lic_0(i)
360			yu10mx(j) = u10m(i, nsrf)
361			yu10my(j) = v10m(i, nsrf)
362			ywindsp(j) = sqrt(yu10mx(j)yu10mx(j)+yu10my(j)yu10my(j))
363	guez	38	END DO
364	guez	3
365	guez	99	! For continent, copy soil water content
366			IF (nsrf == is_ter) THEN
367			yqsol(:knon) = qsol(ni(:knon))
368	guez	62	ELSE
369			yqsol = 0.
370			END IF
371	guez	3
372	guez	62	DO k = 1, nsoilmx
373			DO j = 1, knon
374			i = ni(j)
375			ytsoil(j, k) = ftsoil(i, k, nsrf)
376	guez	38	END DO
377	guez	47	END DO
378	guez	3
379	guez	38	DO k = 1, klev
380			DO j = 1, knon
381			i = ni(j)
382	guez	62	ypaprs(j, k) = paprs(i, k)
383			ypplay(j, k) = pplay(i, k)
384			ydelp(j, k) = delp(i, k)
385			yu(j, k) = u(i, k)
386			yv(j, k) = v(i, k)
387			yt(j, k) = t(i, k)
388			yq(j, k) = q(i, k)
389	guez	38	END DO
390			END DO
391	guez	3
392	guez	62	! calculer Cdrag et les coefficients d'echange
393			CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, &
394			yu, yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))
395			IF (iflag_pbl == 1) THEN
396			CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
397			coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
398			coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
399			END IF
400	guez	3
401	guez	70	! on met un seuil pour coefm et coefh
402	guez	62	IF (nsrf == is_oce) THEN
403			coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)
404			coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)
405	guez	38	END IF
406	guez	3
407	guez	62	IF (ok_kzmin) THEN
408			! Calcul d'une diffusion minimale pour les conditions tres stables
409			CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
410	guez	70	coefm(:knon, 1), ycoefm0, ycoefh0)
411	guez	62	coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
412			coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
413	guez	98	END IF
414	guez	3
415	guez	62	IF (iflag_pbl >= 3) THEN
416	guez	99	! Mellor et Yamada adapté à Mars, Richard Fournier et
417	guez	62	! Frédéric Hourdin
418			yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
419			+ ypplay(:knon, 1))) &
420			* (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg
421			DO k = 2, klev
422			yzlay(1:knon, k) = yzlay(1:knon, k-1) &
423			+ rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
424			/ ypaprs(1:knon, k) &
425			* (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
426			END DO
427			DO k = 1, klev
428			yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
429			/ ypplay(1:knon, k))*rkappa (1.+0.61*yq(1:knon, k))
430			END DO
431			yzlev(1:knon, 1) = 0.
432			yzlev(:knon, klev+1) = 2. * yzlay(:knon, klev) &
433			- yzlay(:knon, klev - 1)
434			DO k = 2, klev
435			yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
436			END DO
437			DO k = 1, klev + 1
438			DO j = 1, knon
439			i = ni(j)
440			yq2(j, k) = q2(i, k, nsrf)
441			END DO
442			END DO
443
444			CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)
445	guez	99	IF (prt_level > 9) PRINT *, 'USTAR = ', yustar
446	guez	62
447			! iflag_pbl peut être utilisé comme longueur de mélange
448
449			IF (iflag_pbl >= 11) THEN
450			CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &
451			yu, yv, yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, &
452			yustar, iflag_pbl)
453			ELSE
454			CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &
455			coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
456			END IF
457
458			coefm(:knon, 2:) = ykmm(:knon, 2:klev)
459			coefh(:knon, 2:) = ykmn(:knon, 2:klev)
460	guez	38	END IF
461	guez	3
462	guez	62	! calculer la diffusion des vitesses "u" et "v"
463	guez	70	CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &
464			ypplay, ydelp, y_d_u, y_flux_u)
465			CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &
466			ypplay, ydelp, y_d_v, y_flux_v)
467	guez	3
468	guez	62	! calculer la diffusion de "q" et de "h"
469	guez	72	CALL clqh(dtime, itap, jour, debut, rlat, knon, nsrf, ni, pctsrf, &
470	guez	99	soil_model, ytsoil, yqsol, rmu0, co2_ppm, yrugos, yrugoro, &
471			yu1, yv1, coefh(:knon, :), yt, yq, yts, ypaprs, ypplay, ydelp, &
472			yrads, yalb, yalblw, ysnow, yqsurf, yrain_f, ysnow_f, yfder, &
473			ysolsw, yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, y_d_ts, &
474			yz0_new, y_flux_t, y_flux_q, y_dflux_t, y_dflux_q, &
475			y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, y_flux_g, &
476			ytslab, y_seaice)
477	guez	3
478	guez	62	! calculer la longueur de rugosite sur ocean
479			yrugm = 0.
480			IF (nsrf == is_oce) THEN
481			DO j = 1, knon
482			yrugm(j) = 0.018coefm(j, 1)(yu1(j)2+yv1(j)2)/rg + &
483			0.1114E-6/sqrt(coefm(j, 1)(yu1(j)2+yv1(j)2))
484			yrugm(j) = max(1.5E-05, yrugm(j))
485			END DO
486			END IF
487	guez	38	DO j = 1, knon
488	guez	62	y_dflux_t(j) = y_dflux_t(j)*ypct(j)
489			y_dflux_q(j) = y_dflux_q(j)*ypct(j)
490			yu1(j) = yu1(j)*ypct(j)
491			yv1(j) = yv1(j)*ypct(j)
492	guez	38	END DO
493	guez	3
494	guez	62	DO k = 1, klev
495			DO j = 1, knon
496			i = ni(j)
497			coefh(j, k) = coefh(j, k)*ypct(j)
498			coefm(j, k) = coefm(j, k)*ypct(j)
499			y_d_t(j, k) = y_d_t(j, k)*ypct(j)
500			y_d_q(j, k) = y_d_q(j, k)*ypct(j)
501			flux_t(i, k, nsrf) = y_flux_t(j, k)
502			flux_q(i, k, nsrf) = y_flux_q(j, k)
503			flux_u(i, k, nsrf) = y_flux_u(j, k)
504			flux_v(i, k, nsrf) = y_flux_v(j, k)
505			y_d_u(j, k) = y_d_u(j, k)*ypct(j)
506			y_d_v(j, k) = y_d_v(j, k)*ypct(j)
507			END DO
508	guez	38	END DO
509	guez	3
510	guez	62	evap(:, nsrf) = -flux_q(:, 1, nsrf)
511	guez	15
512	guez	62	albe(:, nsrf) = 0.
513			alblw(:, nsrf) = 0.
514			snow(:, nsrf) = 0.
515			qsurf(:, nsrf) = 0.
516			rugos(:, nsrf) = 0.
517			fluxlat(:, nsrf) = 0.
518	guez	38	DO j = 1, knon
519			i = ni(j)
520	guez	62	d_ts(i, nsrf) = y_d_ts(j)
521			albe(i, nsrf) = yalb(j)
522			alblw(i, nsrf) = yalblw(j)
523			snow(i, nsrf) = ysnow(j)
524			qsurf(i, nsrf) = yqsurf(j)
525			rugos(i, nsrf) = yz0_new(j)
526			fluxlat(i, nsrf) = yfluxlat(j)
527			IF (nsrf == is_oce) THEN
528			rugmer(i) = yrugm(j)
529			rugos(i, nsrf) = yrugm(j)
530			END IF
531			agesno(i, nsrf) = yagesno(j)
532			fqcalving(i, nsrf) = y_fqcalving(j)
533			ffonte(i, nsrf) = y_ffonte(j)
534			cdragh(i) = cdragh(i) + coefh(j, 1)
535			cdragm(i) = cdragm(i) + coefm(j, 1)
536			dflux_t(i) = dflux_t(i) + y_dflux_t(j)
537			dflux_q(i) = dflux_q(i) + y_dflux_q(j)
538			zu1(i) = zu1(i) + yu1(j)
539			zv1(i) = zv1(i) + yv1(j)
540	guez	38	END DO
541	guez	62	IF (nsrf == is_ter) THEN
542	guez	99	qsol(ni(:knon)) = yqsol(:knon)
543			else IF (nsrf == is_lic) THEN
544	guez	62	DO j = 1, knon
545			i = ni(j)
546			run_off_lic_0(i) = y_run_off_lic_0(j)
547			END DO
548			END IF
549			!$$$ PB ajout pour soil
550			ftsoil(:, :, nsrf) = 0.
551			DO k = 1, nsoilmx
552			DO j = 1, knon
553			i = ni(j)
554			ftsoil(i, k, nsrf) = ytsoil(j, k)
555			END DO
556			END DO
557
558	guez	38	DO j = 1, knon
559			i = ni(j)
560	guez	62	DO k = 1, klev
561			d_t(i, k) = d_t(i, k) + y_d_t(j, k)
562			d_q(i, k) = d_q(i, k) + y_d_q(j, k)
563			d_u(i, k) = d_u(i, k) + y_d_u(j, k)
564			d_v(i, k) = d_v(i, k) + y_d_v(j, k)
565	guez	70	ycoefh(i, k) = ycoefh(i, k) + coefh(j, k)
566	guez	62	END DO
567	guez	38	END DO
568	guez	62
569	guez	99	! diagnostic t, q a 2m et u, v a 10m
570	guez	62
571	guez	38	DO j = 1, knon
572			i = ni(j)
573	guez	62	uzon(j) = yu(j, 1) + y_d_u(j, 1)
574			vmer(j) = yv(j, 1) + y_d_v(j, 1)
575			tair1(j) = yt(j, 1) + y_d_t(j, 1)
576			qair1(j) = yq(j, 1) + y_d_q(j, 1)
577			zgeo1(j) = rdtair1(j)/(0.5(ypaprs(j, 1)+ypplay(j, &
578			1)))*(ypaprs(j, 1)-ypplay(j, 1))
579			tairsol(j) = yts(j) + y_d_ts(j)
580			rugo1(j) = yrugos(j)
581			IF (nsrf == is_oce) THEN
582			rugo1(j) = rugos(i, nsrf)
583			END IF
584			psfce(j) = ypaprs(j, 1)
585			patm(j) = ypplay(j, 1)
586	guez	15
587	guez	62	qairsol(j) = yqsurf(j)
588	guez	38	END DO
589	guez	15
590	guez	62	CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, &
591			zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &
592			yt10m, yq10m, yu10m, yustar)
593	guez	3
594	guez	62	DO j = 1, knon
595			i = ni(j)
596			t2m(i, nsrf) = yt2m(j)
597			q2m(i, nsrf) = yq2m(j)
598	guez	3
599	guez	62	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
600			u10m(i, nsrf) = (yu10m(j)uzon(j))/sqrt(uzon(j)2+vmer(j)*2)
601			v10m(i, nsrf) = (yu10m(j)vmer(j))/sqrt(uzon(j)2+vmer(j)*2)
602	guez	15
603	guez	62	END DO
604	guez	15
605	guez	62	CALL hbtm(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, &
606			y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &
607			ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)
608	guez	15
609	guez	38	DO j = 1, knon
610			i = ni(j)
611	guez	62	pblh(i, nsrf) = ypblh(j)
612			plcl(i, nsrf) = ylcl(j)
613			capcl(i, nsrf) = ycapcl(j)
614			oliqcl(i, nsrf) = yoliqcl(j)
615			cteicl(i, nsrf) = ycteicl(j)
616			pblt(i, nsrf) = ypblt(j)
617			therm(i, nsrf) = ytherm(j)
618			trmb1(i, nsrf) = ytrmb1(j)
619			trmb2(i, nsrf) = ytrmb2(j)
620			trmb3(i, nsrf) = ytrmb3(j)
621	guez	38	END DO
622	guez	3
623	guez	38	DO j = 1, knon
624	guez	62	DO k = 1, klev + 1
625			i = ni(j)
626			q2(i, k, nsrf) = yq2(j, k)
627			END DO
628	guez	38	END DO
629	guez	62	!IM "slab" ocean
630	guez	47	IF (nsrf == is_oce) THEN
631	guez	62	DO j = 1, knon
632			! on projette sur la grille globale
633			i = ni(j)
634			IF (pctsrf_new(i, is_oce)>epsfra) THEN
635			flux_o(i) = y_flux_o(j)
636			ELSE
637			flux_o(i) = 0.
638			END IF
639			END DO
640	guez	38	END IF
641	guez	62
642			IF (nsrf == is_sic) THEN
643			DO j = 1, knon
644			i = ni(j)
645			! On pondère lorsque l'on fait le bilan au sol :
646			IF (pctsrf_new(i, is_sic)>epsfra) THEN
647			flux_g(i) = y_flux_g(j)
648			ELSE
649			flux_g(i) = 0.
650			END IF
651			END DO
652
653			END IF
654			end IF if_knon
655	guez	49	END DO loop_surface
656	guez	15
657	guez	38	! On utilise les nouvelles surfaces
658	guez	15
659	guez	38	rugos(:, is_oce) = rugmer
660			pctsrf = pctsrf_new
661	guez	15
662	guez	38	END SUBROUTINE clmain
663	guez	15
664	guez	38	end module clmain_m