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 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, 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, &
               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	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	70
75			REAL fluxlat(klon, nbsrf)
76
77	guez	101	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	guez	72	REAL, INTENT(IN):: solsw(klon, nbsrf), sollw(klon, nbsrf)
84	guez	154	REAL, intent(in):: fder(klon)
85	guez	191	REAL, intent(inout):: rugos(klon, nbsrf) ! longueur de rugosit\'e (en m)
86	guez	70	real agesno(klon, nbsrf)
87			REAL, INTENT(IN):: rugoro(klon)
88
89	guez	38	REAL d_t(klon, klev), d_q(klon, klev)
90	guez	49	! d_t------output-R- le changement pour "t"
91			! d_q------output-R- le changement pour "q"
92	guez	62
93			REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
94			! changement pour "u" et "v"
95
96	guez	208	REAL, intent(out):: d_ts(klon, nbsrf) ! le changement pour ftsol
97	guez	70
98	guez	206	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	guez	70
102	guez	206	REAL, intent(out):: flux_q(klon, nbsrf)
103			! flux de vapeur d'eau (kg/m2/s) à la surface
104	guez	70
105	guez	206	REAL, intent(out):: flux_u(klon, nbsrf), flux_v(klon, nbsrf)
106			! tension du vent à la surface, en Pa
107
108	guez	70	REAL, INTENT(out):: cdragh(klon), cdragm(klon)
109			real q2(klon, klev+1, nbsrf)
110
111	guez	99	REAL, INTENT(out):: dflux_t(klon), dflux_q(klon)
112	guez	49	! dflux_t derive du flux sensible
113			! dflux_q derive du flux latent
114	guez	191	! IM "slab" ocean
115	guez	70
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	guez	191	! 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	guez	70	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	guez	15
147	guez	202	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	guez	70	REAL y_fqcalving(klon), y_ffonte(klon)
154			real y_run_off_lic_0(klon)
155			REAL rugmer(klon)
156	guez	38	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	guez	49	! on rajoute en output yu1 et yv1 qui sont les vents dans
161			! la premiere couche
162	guez	101	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	guez	99	REAL yfder(klon)
174	guez	38	REAL yrugm(klon), yrads(klon), yrugoro(klon)
175	guez	15
176	guez	38	REAL yfluxlat(klon)
177	guez	15
178	guez	38	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	guez	206	REAL y_flux_t(klon), y_flux_q(klon)
182			REAL y_flux_u(klon), y_flux_v(klon)
183	guez	38	REAL y_dflux_t(klon), y_dflux_q(klon)
184	guez	62	REAL coefh(klon, klev), coefm(klon, klev)
185	guez	38	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	guez	15
189	guez	38	REAL ycoefm0(klon, klev), ycoefh0(klon, klev)
190	guez	15
191	guez	38	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	guez	70	REAL yq2(klon, klev+1)
195	guez	38	REAL q2diag(klon, klev+1)
196	guez	15
197	guez	38	REAL u1lay(klon), v1lay(klon)
198			REAL delp(klon, klev)
199			INTEGER i, k, nsrf
200	guez	15
201	guez	38	INTEGER ni(klon), knon, j
202	guez	40
203	guez	38	REAL pctsrf_pot(klon, nbsrf)
204	guez	145	! "pourcentage potentiel" pour tenir compte des \'eventuelles
205	guez	40	! apparitions ou disparitions de la glace de mer
206	guez	15
207	guez	191	REAL zx_alf1, zx_alf2 ! valeur ambiante par extrapolation
208	guez	15
209	guez	38	REAL yt2m(klon), yq2m(klon), yu10m(klon)
210			REAL yustar(klon)
211	guez	15
212	guez	38	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	guez	15
227	guez	38	REAL qairsol(klon), zgeo1(klon)
228			REAL rugo1(klon)
229	guez	15
230	guez	38	! utiliser un jeu de fonctions simples
231			LOGICAL zxli
232			PARAMETER (zxli=.FALSE.)
233	guez	15
234	guez	38	!------------------------------------------------------------
235	guez	15
236	guez	38	ytherm = 0.
237	guez	15
238	guez	38	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	guez	15
250	guez	40	! 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	guez	38	yrugoro = 0.
279	guez	40	d_ts = 0.
280	guez	38	yfluxlat = 0.
281			flux_t = 0.
282			flux_q = 0.
283			flux_u = 0.
284			flux_v = 0.
285	guez	40	d_t = 0.
286			d_q = 0.
287			d_u = 0.
288			d_v = 0.
289	guez	70	ycoefh = 0.
290	guez	15
291	guez	145	! 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	guez	15
295	guez	202	pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
296			pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
297	guez	38	pctsrf_pot(:, is_oce) = 1. - zmasq
298			pctsrf_pot(:, is_sic) = 1. - zmasq
299	guez	15
300	guez	202	! 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	guez	99	! Boucler sur toutes les sous-fractions du sol:
306
307	guez	49	loop_surface: DO nsrf = 1, nbsrf
308			! Chercher les indices :
309	guez	38	ni = 0
310			knon = 0
311			DO i = 1, klon
312	guez	145	! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
313	guez	38	! "potentielles"
314			IF (pctsrf_pot(i, nsrf) > epsfra) THEN
315			knon = knon + 1
316			ni(knon) = i
317			END IF
318			END DO
319	guez	15
320	guez	62	if_knon: IF (knon /= 0) then
321	guez	38	DO j = 1, knon
322			i = ni(j)
323	guez	62	ypct(j) = pctsrf(i, nsrf)
324	guez	207	yts(j) = ftsol(i, nsrf)
325	guez	62	ysnow(j) = snow(i, nsrf)
326			yqsurf(j) = qsurf(i, nsrf)
327	guez	155	yalb(j) = falbe(i, nsrf)
328	guez	62	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	guez	175	yrads(j) = solsw(i, nsrf) + sollw(i, nsrf)
337	guez	62	ypaprs(j, klev+1) = paprs(i, klev+1)
338			y_run_off_lic_0(j) = run_off_lic_0(i)
339	guez	38	END DO
340	guez	3
341	guez	99	! For continent, copy soil water content
342			IF (nsrf == is_ter) THEN
343			yqsol(:knon) = qsol(ni(:knon))
344	guez	62	ELSE
345			yqsol = 0.
346			END IF
347	guez	3
348	guez	208	ytsoil(:knon, :) = ftsoil(ni(:knon), :, nsrf)
349	guez	3
350	guez	38	DO k = 1, klev
351			DO j = 1, knon
352			i = ni(j)
353	guez	62	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	guez	38	END DO
361			END DO
362	guez	3
363	guez	62	! calculer Cdrag et les coefficients d'echange
364	guez	208	CALL coefkz(nsrf, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, yu, &
365			yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))
366	guez	62	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	guez	3
372	guez	70	! on met un seuil pour coefm et coefh
373	guez	62	IF (nsrf == is_oce) THEN
374			coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)
375			coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)
376	guez	38	END IF
377	guez	3
378	guez	62	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	guez	70	coefm(:knon, 1), ycoefm0, ycoefh0)
382	guez	62	coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
383			coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
384	guez	98	END IF
385	guez	3
386	guez	62	IF (iflag_pbl >= 3) THEN
387	guez	145	! Mellor et Yamada adapt\'e \`a Mars, Richard Fournier et
388			! Fr\'ed\'eric Hourdin
389	guez	62	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	guez	99	IF (prt_level > 9) PRINT *, 'USTAR = ', yustar
417	guez	62
418	guez	145	! iflag_pbl peut \^etre utilis\'e comme longueur de m\'elange
419	guez	62
420			IF (iflag_pbl >= 11) THEN
421	guez	145	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	guez	62	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	guez	38	END IF
432	guez	3
433	guez	62	! calculer la diffusion des vitesses "u" et "v"
434	guez	70	CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &
435	guez	206	ypplay, ydelp, y_d_u, y_flux_u(:knon))
436	guez	70	CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &
437	guez	206	ypplay, ydelp, y_d_v, y_flux_v(:knon))
438	guez	3
439	guez	62	! calculer la diffusion de "q" et de "h"
440	guez	209	CALL clqh(dtime, jour, firstcal, nsrf, ni(:knon), ytsoil(:knon, :), &
441	guez	213	yqsol, mu0, yrugos, yrugoro, yu1, yv1, coefh(:knon, :), yt, &
442	guez	209	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	guez	3
448	guez	62	! 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	guez	38	DO j = 1, knon
458	guez	62	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	guez	38	END DO
463	guez	3
464	guez	62	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	guez	38	END DO
475	guez	3
476	guez	206	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	guez	15
484	guez	206	evap(:, nsrf) = -flux_q(:, nsrf)
485
486	guez	155	falbe(:, nsrf) = 0.
487	guez	62	snow(:, nsrf) = 0.
488			qsurf(:, nsrf) = 0.
489			rugos(:, nsrf) = 0.
490			fluxlat(:, nsrf) = 0.
491	guez	38	DO j = 1, knon
492			i = ni(j)
493	guez	62	d_ts(i, nsrf) = y_d_ts(j)
494	guez	155	falbe(i, nsrf) = yalb(j)
495	guez	62	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	guez	38	END DO
513	guez	62	IF (nsrf == is_ter) THEN
514	guez	99	qsol(ni(:knon)) = yqsol(:knon)
515			else IF (nsrf == is_lic) THEN
516	guez	62	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	guez	118
522	guez	62	ftsoil(:, :, nsrf) = 0.
523	guez	208	ftsoil(ni(:knon), :, nsrf) = ytsoil(:knon, :)
524	guez	62
525	guez	38	DO j = 1, knon
526			i = ni(j)
527	guez	62	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	guez	70	ycoefh(i, k) = ycoefh(i, k) + coefh(j, k)
533	guez	62	END DO
534	guez	38	END DO
535	guez	62
536	guez	99	! diagnostic t, q a 2m et u, v a 10m
537	guez	62
538	guez	38	DO j = 1, knon
539			i = ni(j)
540	guez	62	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	guez	15
554	guez	62	qairsol(j) = yqsurf(j)
555	guez	38	END DO
556	guez	15
557	guez	62	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	guez	3
561	guez	62	DO j = 1, knon
562			i = ni(j)
563			t2m(i, nsrf) = yt2m(j)
564			q2m(i, nsrf) = yq2m(j)
565	guez	3
566	guez	62	! 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	guez	15
571	guez	206	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	guez	15
575	guez	38	DO j = 1, knon
576			i = ni(j)
577	guez	62	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	guez	38	END DO
588	guez	3
589	guez	38	DO j = 1, knon
590	guez	62	DO k = 1, klev + 1
591			i = ni(j)
592			q2(i, k, nsrf) = yq2(j, k)
593			END DO
594	guez	38	END DO
595	guez	62	end IF if_knon
596	guez	49	END DO loop_surface
597	guez	15
598	guez	38	! On utilise les nouvelles surfaces
599			rugos(:, is_oce) = rugmer
600	guez	202	pctsrf(:, is_oce) = pctsrf_new_oce
601			pctsrf(:, is_sic) = pctsrf_new_sic
602	guez	15
603	guez	202	firstcal = .false.
604
605	guez	38	END SUBROUTINE clmain
606	guez	15
607	guez	38	end module clmain_m