Sources/phylmd/clmain.f

module clmain_m

  IMPLICIT NONE

contains

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

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

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

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

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

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

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

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

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

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

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

    REAL fluxlat(klon, nbsrf)

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

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

    REAL, INTENT(IN):: solsw(klon, nbsrf), sollw(klon, nbsrf)
    REAL, intent(in):: fder(klon)
    REAL, intent(inout):: rugos(klon, nbsrf) ! longueur de rugosit\'e (en m)
    real agesno(klon, nbsrf)
    REAL, INTENT(IN):: rugoro(klon)

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

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

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

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

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

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

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

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

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

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

    ! Local:

    LOGICAL:: firstcal = .true.

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

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

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

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

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

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

    REAL yfluxlat(klon)

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

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

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

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

    INTEGER ni(klon), knon, j

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

    REAL zx_alf1, zx_alf2 ! valeur ambiante par extrapolation

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

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

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

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

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

    ytherm = 0.

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

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

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

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

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

    ! Boucler sur toutes les sous-fractions du sol:

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

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

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

          ytsoil(:knon, :) = ftsoil(ni(:knon), :, nsrf)

          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, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, yu, &
               yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))
          IF (iflag_pbl == 1) THEN
             CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
             coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
             coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
          END IF

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

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

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

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

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

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

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

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

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

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

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

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

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

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

          ftsoil(:, :, nsrf) = 0.
          ftsoil(ni(:knon), :, nsrf) = ytsoil(:knon, :)

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

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

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

             qairsol(j) = yqsurf(j)
          END DO

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

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

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

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

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

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

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

    firstcal = .false.

  END SUBROUTINE clmain

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