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, &
       rlat, 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(IN):: rlat(klon) ! latitude en degr\'es

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