Sources/phylmd/clmain.f

module clmain_m

  IMPLICIT NONE

contains

  SUBROUTINE clmain(dtime, pctsrf, t, q, u, v, jour, mu0, 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):: mu0(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, intent(out):: 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.
    flux_t = 0.
    flux_q = 0.
    flux_u = 0.
    flux_v = 0.
    fluxlat = 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, mu0, yrugos, yrugoro, yu1, yv1, coefh(:knon, :), yt, &
               yq, yts(:knon), ypaprs, ypplay, ydelp, yrads, yalb(:knon), &
               ysnow, yqsurf, yrain_f, ysnow_f, yfder, yfluxlat(:knon), &
               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

          flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
          flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
          flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
          flux_v(ni(:knon), nsrf) = y_flux_v(:knon)

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

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