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, 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)

    ! 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 stdlevvar_m, only: stdlevvar
    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) ! 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 ftsoil(klon, nsoilmx, nbsrf)

    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 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)
    ! 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 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, intent(out):: d_ts(klon, nbsrf) ! 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

    ! Local:

    REAL y_flux_o(klon), y_flux_g(klon)
    real ytslab(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)

    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 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.
    yu10mx = 0.
    yu10my = 0.
    ywindsp = 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è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, &
               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)

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