trunk/phylmd/pbl_surface.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, falbe, 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)

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

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

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

    use clqh_m, only: clqh
    use clvent_m, only: clvent
    use coefkz_m, only: coefkz
    use coefkzmin_m, only: coefkzmin
    USE conf_gcm_m, ONLY: prt_level
    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, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)
    ! soil temperature of surface fraction

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

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

    REAL fluxlat(klon, nbsrf)

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

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

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

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

    ! Local:

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

    REAL rugmer(klon)

    REAL ytsoil(klon, nsoilmx)

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

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

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

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

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

    REAL yfluxlat(klon)

    REAL y_d_ts(klon)
    REAL y_d_t(klon, klev), y_d_q(klon, klev)
    REAL y_d_u(klon, klev), y_d_v(klon, klev)
    REAL y_flux_t(klon, 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 \'eventuelles
    ! 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.
    yrain_f = 0.
    ysnow_f = 0.
    yfder = 0.
    yrugos = 0.
    yu1 = 0.
    yv1 = 0.
    yrads = 0.
    ypaprs = 0.
    ypplay = 0.
    ydelp = 0.
    yu = 0.
    yv = 0.
    yt = 0.
    yq = 0.
    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\`ere ici qu'on
    ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
    ! (\`a 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\'eterminer le domaine \`a traiter, on utilise les surfaces
          ! "potentielles"
          IF (pctsrf_pot(i, nsrf) > epsfra) THEN
             knon = knon + 1
             ni(knon) = i
          END IF
       END DO

       if_knon: IF (knon /= 0) then
          DO j = 1, knon
             i = ni(j)
             ypct(j) = pctsrf(i, nsrf)
             yts(j) = ts(i, nsrf)
             ysnow(j) = snow(i, nsrf)
             yqsurf(j) = qsurf(i, nsrf)
             yalb(j) = falbe(i, nsrf)
             yrain_f(j) = rain_fall(i)
             ysnow_f(j) = snow_f(i)
             yagesno(j) = agesno(i, nsrf)
             yfder(j) = fder(i)
             yrugos(j) = rugos(i, nsrf)
             yrugoro(j) = rugoro(i)
             yu1(j) = u1lay(i)
             yv1(j) = v1lay(i)
             yrads(j) = solsw(i, nsrf) + sollw(i, nsrf)
             ypaprs(j, klev+1) = paprs(i, klev+1)
             y_run_off_lic_0(j) = run_off_lic_0(i)
             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\'e \`a Mars, Richard Fournier et
             ! Fr\'ed\'eric Hourdin
             yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
                  + ypplay(:knon, 1))) &
                  * (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg
             DO k = 2, klev
                yzlay(1:knon, k) = yzlay(1:knon, k-1) &
                     + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
                     / ypaprs(1:knon, k) &
                     * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
             END DO
             DO k = 1, klev
                yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
                     / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))
             END DO
             yzlev(1:knon, 1) = 0.
             yzlev(:knon, klev+1) = 2. * yzlay(:knon, klev) &
                  - yzlay(:knon, klev - 1)
             DO k = 2, klev
                yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
             END DO
             DO k = 1, klev + 1
                DO j = 1, knon
                   i = ni(j)
                   yq2(j, k) = q2(i, k, nsrf)
                END DO
             END DO

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

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

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

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

          ! calculer la diffusion des vitesses "u" et "v"
          CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &
               ypplay, ydelp, y_d_u, y_flux_u)
          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(:knon), &
               pctsrf, ytsoil, yqsol, rmu0, co2_ppm, yrugos, yrugoro, yu1, &
               yv1, coefh(:knon, :), yt, yq, yts, ypaprs, ypplay, ydelp, &
               yrads, yalb(:knon), ysnow, yqsurf, yrain_f, ysnow_f, yfder, &
               yfluxlat, pctsrf_new, yagesno(:knon), y_d_t, y_d_q, &
               y_d_ts(:knon), yz0_new, y_flux_t, y_flux_q, y_dflux_t, &
               y_dflux_q, y_fqcalving, y_ffonte, y_run_off_lic_0)

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

          DO k = 1, klev
             DO j = 1, knon
                i = ni(j)
                coefh(j, k) = coefh(j, k)*ypct(j)
                coefm(j, k) = coefm(j, k)*ypct(j)
                y_d_t(j, k) = y_d_t(j, k)*ypct(j)
                y_d_q(j, k) = y_d_q(j, k)*ypct(j)
                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)

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

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

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

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

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

             qairsol(j) = yqsurf(j)
          END DO

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

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

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

          END DO

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