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

    ! 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.
    ! -- 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, &
               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			USE suphec_m, ONLY: rd, rg, rkappa
42			use ustarhb_m, only: ustarhb
43			use vdif_kcay_m, only: vdif_kcay
44	guez	47	use yamada4_m, only: yamada4
45	guez	15
46	guez	62	REAL, INTENT(IN):: dtime ! interval du temps (secondes)
47			INTEGER, INTENT(IN):: itap ! numero du pas de temps
48			REAL, INTENT(inout):: pctsrf(klon, nbsrf)
49
50			! la nouvelle repartition des surfaces sortie de l'interface
51			REAL, INTENT(out):: pctsrf_new(klon, nbsrf)
52
53			REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
54			REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg/kg)
55			REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
56			INTEGER, INTENT(IN):: jour ! jour de l'annee en cours
57			REAL, intent(in):: rmu0(klon) ! cosinus de l'angle solaire zenithal
58	guez	99	REAL, intent(in):: co2_ppm ! taux CO2 atmosphere
59			REAL, INTENT(IN):: ts(klon, nbsrf) ! input-R- temperature du sol (en Kelvin)
60	guez	71	REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh
61	guez	99	REAL, INTENT(IN):: ksta, ksta_ter
62			LOGICAL, INTENT(IN):: ok_kzmin
63	guez	70	REAL ftsoil(klon, nsoilmx, nbsrf)
64	guez	101
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	70	REAL snow(klon, nbsrf)
71			REAL qsurf(klon, nbsrf)
72			REAL evap(klon, nbsrf)
73			REAL albe(klon, nbsrf)
74			REAL alblw(klon, nbsrf)
75
76			REAL fluxlat(klon, nbsrf)
77
78	guez	101	REAL, intent(in):: rain_fall(klon)
79			! liquid water mass flux (kg/m2/s), positive down
80
81			REAL, intent(in):: snow_f(klon)
82			! solid water mass flux (kg/m2/s), positive down
83
84	guez	72	REAL, INTENT(IN):: solsw(klon, nbsrf), sollw(klon, nbsrf)
85	guez	70	REAL fder(klon)
86	guez	62	REAL, INTENT(IN):: rlat(klon) ! latitude en degrés
87	guez	70
88			REAL rugos(klon, nbsrf)
89			! rugos----input-R- longeur de rugosite (en m)
90
91			LOGICAL, INTENT(IN):: debut
92			real agesno(klon, nbsrf)
93			REAL, INTENT(IN):: rugoro(klon)
94
95	guez	38	REAL d_t(klon, klev), d_q(klon, klev)
96	guez	49	! d_t------output-R- le changement pour "t"
97			! d_q------output-R- le changement pour "q"
98	guez	62
99			REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
100			! changement pour "u" et "v"
101
102	guez	70	REAL d_ts(klon, nbsrf)
103			! d_ts-----output-R- le changement pour "ts"
104
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			! -- LOOP
309	guez	40	yu10mx = 0.
310			yu10my = 0.
311			ywindsp = 0.
312	guez	38	! -- LOOP
313	guez	40	d_ts = 0.
314	guez	38	yfluxlat = 0.
315			flux_t = 0.
316			flux_q = 0.
317			flux_u = 0.
318			flux_v = 0.
319	guez	40	d_t = 0.
320			d_q = 0.
321			d_u = 0.
322			d_v = 0.
323	guez	70	ycoefh = 0.
324	guez	15
325	guez	38	! Initialisation des "pourcentages potentiels". On considère ici qu'on
326			! peut avoir potentiellement de la glace sur tout le domaine océanique
327			! (à affiner)
328	guez	15
329	guez	38	pctsrf_pot = pctsrf
330			pctsrf_pot(:, is_oce) = 1. - zmasq
331			pctsrf_pot(:, is_sic) = 1. - zmasq
332	guez	15
333	guez	99	! Boucler sur toutes les sous-fractions du sol:
334
335	guez	49	loop_surface: DO nsrf = 1, nbsrf
336			! Chercher les indices :
337	guez	38	ni = 0
338			knon = 0
339			DO i = 1, klon
340	guez	40	! Pour déterminer le domaine à traiter, on utilise les surfaces
341	guez	38	! "potentielles"
342			IF (pctsrf_pot(i, nsrf) > epsfra) THEN
343			knon = knon + 1
344			ni(knon) = i
345			END IF
346			END DO
347	guez	15
348	guez	62	if_knon: IF (knon /= 0) then
349	guez	38	DO j = 1, knon
350			i = ni(j)
351	guez	62	ypct(j) = pctsrf(i, nsrf)
352			yts(j) = ts(i, nsrf)
353			ytslab(i) = tslab(i)
354			ysnow(j) = snow(i, nsrf)
355			yqsurf(j) = qsurf(i, nsrf)
356			yalb(j) = albe(i, nsrf)
357			yalblw(j) = alblw(i, nsrf)
358			yrain_f(j) = rain_fall(i)
359			ysnow_f(j) = snow_f(i)
360			yagesno(j) = agesno(i, nsrf)
361			yfder(j) = fder(i)
362			ysolsw(j) = solsw(i, nsrf)
363			ysollw(j) = sollw(i, nsrf)
364			yrugos(j) = rugos(i, nsrf)
365			yrugoro(j) = rugoro(i)
366			yu1(j) = u1lay(i)
367			yv1(j) = v1lay(i)
368			yrads(j) = ysolsw(j) + ysollw(j)
369			ypaprs(j, klev+1) = paprs(i, klev+1)
370			y_run_off_lic_0(j) = run_off_lic_0(i)
371			yu10mx(j) = u10m(i, nsrf)
372			yu10my(j) = v10m(i, nsrf)
373			ywindsp(j) = sqrt(yu10mx(j)yu10mx(j)+yu10my(j)yu10my(j))
374	guez	38	END DO
375	guez	3
376	guez	99	! For continent, copy soil water content
377			IF (nsrf == is_ter) THEN
378			yqsol(:knon) = qsol(ni(:knon))
379	guez	62	ELSE
380			yqsol = 0.
381			END IF
382	guez	3
383	guez	62	DO k = 1, nsoilmx
384			DO j = 1, knon
385			i = ni(j)
386			ytsoil(j, k) = ftsoil(i, k, nsrf)
387	guez	38	END DO
388	guez	47	END DO
389	guez	3
390	guez	38	DO k = 1, klev
391			DO j = 1, knon
392			i = ni(j)
393	guez	62	ypaprs(j, k) = paprs(i, k)
394			ypplay(j, k) = pplay(i, k)
395			ydelp(j, k) = delp(i, k)
396			yu(j, k) = u(i, k)
397			yv(j, k) = v(i, k)
398			yt(j, k) = t(i, k)
399			yq(j, k) = q(i, k)
400	guez	38	END DO
401			END DO
402	guez	3
403	guez	62	! calculer Cdrag et les coefficients d'echange
404			CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, &
405			yu, yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))
406			IF (iflag_pbl == 1) THEN
407			CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
408			coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
409			coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
410			END IF
411	guez	3
412	guez	70	! on met un seuil pour coefm et coefh
413	guez	62	IF (nsrf == is_oce) THEN
414			coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)
415			coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)
416	guez	38	END IF
417	guez	3
418	guez	62	IF (ok_kzmin) THEN
419			! Calcul d'une diffusion minimale pour les conditions tres stables
420			CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
421	guez	70	coefm(:knon, 1), ycoefm0, ycoefh0)
422	guez	62	coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
423			coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
424	guez	98	END IF
425	guez	3
426	guez	62	IF (iflag_pbl >= 3) THEN
427	guez	99	! Mellor et Yamada adapté à Mars, Richard Fournier et
428	guez	62	! Frédéric Hourdin
429			yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
430			+ ypplay(:knon, 1))) &
431			* (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg
432			DO k = 2, klev
433			yzlay(1:knon, k) = yzlay(1:knon, k-1) &
434			+ rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
435			/ ypaprs(1:knon, k) &
436			* (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
437			END DO
438			DO k = 1, klev
439			yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
440			/ ypplay(1:knon, k))*rkappa (1.+0.61*yq(1:knon, k))
441			END DO
442			yzlev(1:knon, 1) = 0.
443			yzlev(:knon, klev+1) = 2. * yzlay(:knon, klev) &
444			- yzlay(:knon, klev - 1)
445			DO k = 2, klev
446			yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
447			END DO
448			DO k = 1, klev + 1
449			DO j = 1, knon
450			i = ni(j)
451			yq2(j, k) = q2(i, k, nsrf)
452			END DO
453			END DO
454
455			CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)
456	guez	99	IF (prt_level > 9) PRINT *, 'USTAR = ', yustar
457	guez	62
458			! iflag_pbl peut être utilisé comme longueur de mélange
459
460			IF (iflag_pbl >= 11) THEN
461			CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &
462			yu, yv, yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, &
463			yustar, iflag_pbl)
464			ELSE
465			CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &
466			coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
467			END IF
468
469			coefm(:knon, 2:) = ykmm(:knon, 2:klev)
470			coefh(:knon, 2:) = ykmn(:knon, 2:klev)
471	guez	38	END IF
472	guez	3
473	guez	62	! calculer la diffusion des vitesses "u" et "v"
474	guez	70	CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &
475			ypplay, ydelp, y_d_u, y_flux_u)
476			CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &
477			ypplay, ydelp, y_d_v, y_flux_v)
478	guez	3
479	guez	62	! calculer la diffusion de "q" et de "h"
480	guez	72	CALL clqh(dtime, itap, jour, debut, rlat, knon, nsrf, ni, pctsrf, &
481	guez	101	ytsoil, yqsol, rmu0, co2_ppm, yrugos, yrugoro, &
482	guez	99	yu1, yv1, coefh(:knon, :), yt, yq, yts, ypaprs, ypplay, ydelp, &
483			yrads, yalb, yalblw, ysnow, yqsurf, yrain_f, ysnow_f, yfder, &
484			ysolsw, yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, y_d_ts, &
485			yz0_new, y_flux_t, y_flux_q, y_dflux_t, y_dflux_q, &
486	guez	101	y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, y_flux_g)
487	guez	3
488	guez	62	! calculer la longueur de rugosite sur ocean
489			yrugm = 0.
490			IF (nsrf == is_oce) THEN
491			DO j = 1, knon
492			yrugm(j) = 0.018coefm(j, 1)(yu1(j)2+yv1(j)2)/rg + &
493			0.1114E-6/sqrt(coefm(j, 1)(yu1(j)2+yv1(j)2))
494			yrugm(j) = max(1.5E-05, yrugm(j))
495			END DO
496			END IF
497	guez	38	DO j = 1, knon
498	guez	62	y_dflux_t(j) = y_dflux_t(j)*ypct(j)
499			y_dflux_q(j) = y_dflux_q(j)*ypct(j)
500			yu1(j) = yu1(j)*ypct(j)
501			yv1(j) = yv1(j)*ypct(j)
502	guez	38	END DO
503	guez	3
504	guez	62	DO k = 1, klev
505			DO j = 1, knon
506			i = ni(j)
507			coefh(j, k) = coefh(j, k)*ypct(j)
508			coefm(j, k) = coefm(j, k)*ypct(j)
509			y_d_t(j, k) = y_d_t(j, k)*ypct(j)
510			y_d_q(j, k) = y_d_q(j, k)*ypct(j)
511			flux_t(i, k, nsrf) = y_flux_t(j, k)
512			flux_q(i, k, nsrf) = y_flux_q(j, k)
513			flux_u(i, k, nsrf) = y_flux_u(j, k)
514			flux_v(i, k, nsrf) = y_flux_v(j, k)
515			y_d_u(j, k) = y_d_u(j, k)*ypct(j)
516			y_d_v(j, k) = y_d_v(j, k)*ypct(j)
517			END DO
518	guez	38	END DO
519	guez	3
520	guez	62	evap(:, nsrf) = -flux_q(:, 1, nsrf)
521	guez	15
522	guez	62	albe(:, nsrf) = 0.
523			alblw(:, nsrf) = 0.
524			snow(:, nsrf) = 0.
525			qsurf(:, nsrf) = 0.
526			rugos(:, nsrf) = 0.
527			fluxlat(:, nsrf) = 0.
528	guez	38	DO j = 1, knon
529			i = ni(j)
530	guez	62	d_ts(i, nsrf) = y_d_ts(j)
531			albe(i, nsrf) = yalb(j)
532			alblw(i, nsrf) = yalblw(j)
533			snow(i, nsrf) = ysnow(j)
534			qsurf(i, nsrf) = yqsurf(j)
535			rugos(i, nsrf) = yz0_new(j)
536			fluxlat(i, nsrf) = yfluxlat(j)
537			IF (nsrf == is_oce) THEN
538			rugmer(i) = yrugm(j)
539			rugos(i, nsrf) = yrugm(j)
540			END IF
541			agesno(i, nsrf) = yagesno(j)
542			fqcalving(i, nsrf) = y_fqcalving(j)
543			ffonte(i, nsrf) = y_ffonte(j)
544			cdragh(i) = cdragh(i) + coefh(j, 1)
545			cdragm(i) = cdragm(i) + coefm(j, 1)
546			dflux_t(i) = dflux_t(i) + y_dflux_t(j)
547			dflux_q(i) = dflux_q(i) + y_dflux_q(j)
548			zu1(i) = zu1(i) + yu1(j)
549			zv1(i) = zv1(i) + yv1(j)
550	guez	38	END DO
551	guez	62	IF (nsrf == is_ter) THEN
552	guez	99	qsol(ni(:knon)) = yqsol(:knon)
553			else IF (nsrf == is_lic) THEN
554	guez	62	DO j = 1, knon
555			i = ni(j)
556			run_off_lic_0(i) = y_run_off_lic_0(j)
557			END DO
558			END IF
559			!$$$ PB ajout pour soil
560			ftsoil(:, :, nsrf) = 0.
561			DO k = 1, nsoilmx
562			DO j = 1, knon
563			i = ni(j)
564			ftsoil(i, k, nsrf) = ytsoil(j, k)
565			END DO
566			END DO
567
568	guez	38	DO j = 1, knon
569			i = ni(j)
570	guez	62	DO k = 1, klev
571			d_t(i, k) = d_t(i, k) + y_d_t(j, k)
572			d_q(i, k) = d_q(i, k) + y_d_q(j, k)
573			d_u(i, k) = d_u(i, k) + y_d_u(j, k)
574			d_v(i, k) = d_v(i, k) + y_d_v(j, k)
575	guez	70	ycoefh(i, k) = ycoefh(i, k) + coefh(j, k)
576	guez	62	END DO
577	guez	38	END DO
578	guez	62
579	guez	99	! diagnostic t, q a 2m et u, v a 10m
580	guez	62
581	guez	38	DO j = 1, knon
582			i = ni(j)
583	guez	62	uzon(j) = yu(j, 1) + y_d_u(j, 1)
584			vmer(j) = yv(j, 1) + y_d_v(j, 1)
585			tair1(j) = yt(j, 1) + y_d_t(j, 1)
586			qair1(j) = yq(j, 1) + y_d_q(j, 1)
587			zgeo1(j) = rdtair1(j)/(0.5(ypaprs(j, 1)+ypplay(j, &
588			1)))*(ypaprs(j, 1)-ypplay(j, 1))
589			tairsol(j) = yts(j) + y_d_ts(j)
590			rugo1(j) = yrugos(j)
591			IF (nsrf == is_oce) THEN
592			rugo1(j) = rugos(i, nsrf)
593			END IF
594			psfce(j) = ypaprs(j, 1)
595			patm(j) = ypplay(j, 1)
596	guez	15
597	guez	62	qairsol(j) = yqsurf(j)
598	guez	38	END DO
599	guez	15
600	guez	62	CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, &
601			zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &
602			yt10m, yq10m, yu10m, yustar)
603	guez	3
604	guez	62	DO j = 1, knon
605			i = ni(j)
606			t2m(i, nsrf) = yt2m(j)
607			q2m(i, nsrf) = yq2m(j)
608	guez	3
609	guez	62	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
610			u10m(i, nsrf) = (yu10m(j)uzon(j))/sqrt(uzon(j)2+vmer(j)*2)
611			v10m(i, nsrf) = (yu10m(j)vmer(j))/sqrt(uzon(j)2+vmer(j)*2)
612	guez	15
613	guez	62	END DO
614	guez	15
615	guez	62	CALL hbtm(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, &
616			y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &
617			ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)
618	guez	15
619	guez	38	DO j = 1, knon
620			i = ni(j)
621	guez	62	pblh(i, nsrf) = ypblh(j)
622			plcl(i, nsrf) = ylcl(j)
623			capcl(i, nsrf) = ycapcl(j)
624			oliqcl(i, nsrf) = yoliqcl(j)
625			cteicl(i, nsrf) = ycteicl(j)
626			pblt(i, nsrf) = ypblt(j)
627			therm(i, nsrf) = ytherm(j)
628			trmb1(i, nsrf) = ytrmb1(j)
629			trmb2(i, nsrf) = ytrmb2(j)
630			trmb3(i, nsrf) = ytrmb3(j)
631	guez	38	END DO
632	guez	3
633	guez	38	DO j = 1, knon
634	guez	62	DO k = 1, klev + 1
635			i = ni(j)
636			q2(i, k, nsrf) = yq2(j, k)
637			END DO
638	guez	38	END DO
639	guez	62	!IM "slab" ocean
640	guez	47	IF (nsrf == is_oce) THEN
641	guez	62	DO j = 1, knon
642			! on projette sur la grille globale
643			i = ni(j)
644			IF (pctsrf_new(i, is_oce)>epsfra) THEN
645			flux_o(i) = y_flux_o(j)
646			ELSE
647			flux_o(i) = 0.
648			END IF
649			END DO
650	guez	38	END IF
651	guez	62
652			IF (nsrf == is_sic) THEN
653			DO j = 1, knon
654			i = ni(j)
655			! On pondère lorsque l'on fait le bilan au sol :
656			IF (pctsrf_new(i, is_sic)>epsfra) THEN
657			flux_g(i) = y_flux_g(j)
658			ELSE
659			flux_g(i) = 0.
660			END IF
661			END DO
662
663			END IF
664			end IF if_knon
665	guez	49	END DO loop_surface
666	guez	15
667	guez	38	! On utilise les nouvelles surfaces
668	guez	15
669	guez	38	rugos(:, is_oce) = rugmer
670			pctsrf = pctsrf_new
671	guez	15
672	guez	38	END SUBROUTINE clmain
673	guez	15
674	guez	38	end module clmain_m