/[lmdze]/trunk/phylmd/Interface_surf/pbl_surface.f
ViewVC logotype

Diff of /trunk/phylmd/Interface_surf/pbl_surface.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/phylmd/clmain.f90 revision 37 by guez, Tue Dec 21 15:45:48 2010 UTC trunk/phylmd/Interface_surf/pbl_surface.f revision 303 by guez, Thu Sep 6 14:25:07 2018 UTC
# Line 1  Line 1 
1  SUBROUTINE clmain(dtime, itap, date0, pctsrf, pctsrf_new, t, q, u, v,&  module pbl_surface_m
      jour, rmu0, co2_ppm, ok_veget, ocean, npas, nexca, ts,&  
      soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil,&  
      qsol, paprs, pplay, snow, qsurf, evap, albe, alblw, fluxlat,&  
      rain_f, snow_f, solsw, sollw, sollwdown, fder, rlon, rlat, cufi,&  
      cvfi, rugos, debut, lafin, 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, zcoefh, zu1, zv1, t2m, q2m, u10m, v10m, pblh,&  
      capcl, oliqcl, cteicl, pblt, therm, trmb1, trmb2, trmb3, plcl,&  
      fqcalving, ffonte, run_off_lic_0, flux_o, flux_g, tslab, seaice)  
   
   ! From phylmd/clmain.F, v 1.6 2005/11/16 14:47:19  
   
   !AA Tout ce qui a trait au traceurs est dans phytrac maintenant  
   !AA pour l'instant le calcul de la couche limite pour les traceurs  
   !AA se fait avec cltrac et ne tient pas compte de la differentiation  
   !AA des sous-fraction de sol.  
   
   !AA Pour pouvoir extraire les coefficient d'echanges et le vent  
   !AA dans la premiere couche, 3 champs supplementaires ont ete crees  
   !AA zcoefh, zu1 et zv1. Pour l'instant nous avons moyenne les valeurs  
   !AA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir  
   !AA si les informations des subsurfaces doivent etre prises en compte  
   !AA il faudra sortir ces memes champs en leur ajoutant une dimension,  
   !AA c'est a dire nbsrf (nbre de subsurface).  
   
   ! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818  
   ! Objet: interface de "couche limite" (diffusion verticale)  
   
   ! Arguments:  
   ! dtime----input-R- interval du temps (secondes)  
   ! itap-----input-I- numero du pas de temps  
   ! date0----input-R- jour initial  
   ! t--------input-R- temperature (K)  
   ! q--------input-R- vapeur d'eau (kg/kg)  
   ! u--------input-R- vitesse u  
   ! v--------input-R- vitesse v  
   ! ts-------input-R- temperature du sol (en Kelvin)  
   ! paprs----input-R- pression a intercouche (Pa)  
   ! pplay----input-R- pression au milieu de couche (Pa)  
   ! radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2  
   ! rlat-----input-R- latitude en degree  
   ! rugos----input-R- longeur de rugosite (en m)  
   ! cufi-----input-R- resolution des mailles en x (m)  
   ! cvfi-----input-R- resolution des mailles en y (m)  
   
   ! d_t------output-R- le changement pour "t"  
   ! d_q------output-R- le changement pour "q"  
   ! d_u------output-R- le changement pour "u"  
   ! d_v------output-R- le changement pour "v"  
   ! d_ts-----output-R- le changement pour "ts"  
   ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)  
   !                    (orientation positive vers le bas)  
   ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)  
   ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal  
   ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal  
   ! dflux_t derive du flux sensible  
   ! dflux_q derive du flux latent  
   !IM "slab" ocean  
   ! flux_g---output-R-  flux glace (pour OCEAN='slab  ')  
   ! flux_o---output-R-  flux ocean (pour OCEAN='slab  ')  
   ! tslab-in/output-R temperature du slab ocean (en Kelvin) ! uniqmnt pour slab  
   ! seaice---output-R-  glace de mer (kg/m2) (pour OCEAN='slab  ')  
   !cc  
   ! 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  
   !AA on rajoute en output yu1 et yv1 qui sont les vents dans  
   !AA la premiere couche  
   !AA ces 4 variables sont maintenant traites dans phytrac  
   ! itr--------input-I- nombre de traceurs  
   ! tr---------input-R- q. de traceurs  
   ! flux_surf--input-R- flux de traceurs a la surface  
   ! d_tr-------output-R tendance de traceurs  
   !IM cf. AM : PBL  
   ! trmb1-------deep_cape  
   ! trmb2--------inhibition  
   ! trmb3-------Point Omega  
   ! Cape(klon)-------Cape du thermique  
   ! EauLiq(klon)-------Eau liqu integr du thermique  
   ! ctei(klon)-------Critere d'instab d'entrainmt des nuages de CL  
   ! lcl------- Niveau de condensation  
   ! pblh------- HCL  
   ! pblT------- T au nveau HCL  
   
   !$$$ PB ajout pour soil  
   
   USE histcom, ONLY : histbeg_totreg, histdef, histend, histsync  
   use histwrite_m, only: histwrite  
   use calendar, ONLY : ymds2ju  
   USE dimens_m, ONLY : iim, jjm  
   USE indicesol, ONLY : epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf  
   USE dimphy, ONLY : klev, klon, zmasq  
   USE dimsoil, ONLY : nsoilmx  
   USE temps, ONLY : annee_ref, itau_phy  
   USE dynetat0_m, ONLY : day_ini  
   USE iniprint, ONLY : prt_level  
   USE yomcst, ONLY : rd, rg, rkappa  
   USE conf_phys_m, ONLY : iflag_pbl  
   USE gath_cpl, ONLY : gath2cpl  
   use hbtm_m, only: hbtm  
2    
3    IMPLICIT NONE    IMPLICIT NONE
4    
5    REAL, INTENT (IN) :: dtime  contains
   REAL date0  
   INTEGER, INTENT (IN) :: itap  
   REAL t(klon, klev), q(klon, klev)  
   REAL u(klon, klev), v(klon, klev)  
   REAL, INTENT (IN) :: paprs(klon, klev+1)  
   REAL, INTENT (IN) :: pplay(klon, klev)  
   REAL, INTENT (IN) :: rlon(klon), rlat(klon)  
   REAL cufi(klon), cvfi(klon)  
   REAL d_t(klon, klev), d_q(klon, klev)  
   REAL d_u(klon, klev), d_v(klon, klev)  
   REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)  
   REAL dflux_t(klon), dflux_q(klon)  
   !IM "slab" ocean  
   REAL flux_o(klon), flux_g(klon)  
   REAL y_flux_o(klon), y_flux_g(klon)  
   REAL tslab(klon), ytslab(klon)  
   REAL seaice(klon), y_seaice(klon)  
   REAL y_fqcalving(klon), y_ffonte(klon)  
   REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)  
   REAL run_off_lic_0(klon), y_run_off_lic_0(klon)  
   
   REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)  
   REAL rugmer(klon), agesno(klon, nbsrf)  
   REAL, INTENT (IN) :: rugoro(klon)  
   REAL cdragh(klon), cdragm(klon)  
   ! jour de l'annee en cours                  
   INTEGER jour  
   REAL rmu0(klon) ! cosinus de l'angle solaire zenithal      
   ! taux CO2 atmosphere                      
   REAL co2_ppm  
   LOGICAL, INTENT (IN) :: debut  
   LOGICAL, INTENT (IN) :: lafin  
   LOGICAL ok_veget  
   CHARACTER (len=*), INTENT (IN) :: ocean  
   INTEGER npas, nexca  
   
   REAL pctsrf(klon, nbsrf)  
   REAL ts(klon, nbsrf)  
   REAL d_ts(klon, nbsrf)  
   REAL snow(klon, nbsrf)  
   REAL qsurf(klon, nbsrf)  
   REAL evap(klon, nbsrf)  
   REAL albe(klon, nbsrf)  
   REAL alblw(klon, nbsrf)  
   
   REAL fluxlat(klon, nbsrf)  
   
   REAL rain_f(klon), snow_f(klon)  
   REAL fder(klon)  
   
   REAL sollw(klon, nbsrf), solsw(klon, nbsrf), sollwdown(klon)  
   REAL rugos(klon, nbsrf)  
   ! la nouvelle repartition des surfaces sortie de l'interface  
   REAL pctsrf_new(klon, nbsrf)  
   
   REAL zcoefh(klon, klev)  
   REAL zu1(klon)  
   REAL zv1(klon)  
   
   !$$$ PB ajout pour soil  
   LOGICAL, INTENT (IN) :: soil_model  
   !IM ajout seuils cdrm, cdrh  
   REAL cdmmax, cdhmax  
   
   REAL ksta, ksta_ter  
   LOGICAL ok_kzmin  
   
   REAL ftsoil(klon, nsoilmx, nbsrf)  
   REAL ytsoil(klon, nsoilmx)  
   REAL qsol(klon)  
   
   EXTERNAL clqh, clvent, coefkz, calbeta, cltrac  
   
   REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)  
   REAL yalb(klon)  
   REAL yalblw(klon)  
   REAL yu1(klon), yv1(klon)  
   REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)  
   REAL yrain_f(klon), ysnow_f(klon)  
   REAL ysollw(klon), ysolsw(klon), ysollwdown(klon)  
   REAL yfder(klon), ytaux(klon), ytauy(klon)  
   REAL yrugm(klon), yrads(klon), yrugoro(klon)  
   
   REAL yfluxlat(klon)  
   
   REAL y_d_ts(klon)  
   REAL y_d_t(klon, klev), y_d_q(klon, klev)  
   REAL y_d_u(klon, klev), y_d_v(klon, klev)  
   REAL y_flux_t(klon, klev), y_flux_q(klon, klev)  
   REAL y_flux_u(klon, klev), y_flux_v(klon, klev)  
   REAL y_dflux_t(klon), y_dflux_q(klon)  
   REAL ycoefh(klon, klev), ycoefm(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)  
   
   LOGICAL ok_nonloc  
   PARAMETER (ok_nonloc=.FALSE.)  
   REAL ycoefm0(klon, klev), ycoefh0(klon, klev)  
   
   !IM 081204 hcl_Anne ? BEG  
   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), q2(klon, klev+1, nbsrf)  
   REAL q2diag(klon, klev+1)  
   !IM 081204 hcl_Anne ? END  
   
   REAL u1lay(klon), v1lay(klon)  
   REAL delp(klon, klev)  
   INTEGER i, k, nsrf  
   
   INTEGER ni(klon), knon, j  
   ! Introduction d'une variable "pourcentage potentiel" pour tenir compte  
   ! des eventuelles apparitions et/ou disparitions de la glace de mer  
   REAL pctsrf_pot(klon, nbsrf)  
   
   REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.  
   
   ! maf pour sorties IOISPL en cas de debugagage  
   
   CHARACTER (80) cldebug  
   SAVE cldebug  
   CHARACTER (8) cl_surf(nbsrf)  
   SAVE cl_surf  
   INTEGER nhoridbg, nidbg  
   SAVE nhoridbg, nidbg  
   INTEGER ndexbg(iim*(jjm+1))  
   REAL zx_lon(iim, jjm+1), zx_lat(iim, jjm+1), zjulian  
   REAL tabindx(klon)  
   REAL debugtab(iim, jjm+1)  
   LOGICAL first_appel  
   SAVE first_appel  
   DATA first_appel/ .TRUE./  
   LOGICAL :: debugindex = .FALSE.  
   INTEGER idayref  
   REAL t2m(klon, nbsrf), q2m(klon, nbsrf)  
   REAL u10m(klon, nbsrf), v10m(klon, nbsrf)  
   
   REAL yt2m(klon), yq2m(klon), yu10m(klon)  
   REAL yustar(klon)  
   ! -- LOOP  
   REAL yu10mx(klon)  
   REAL yu10my(klon)  
   REAL ywindsp(klon)  
   ! -- LOOP  
   
   REAL yt10m(klon), yq10m(klon)  
   !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds  
   ! physiq ce qui permet de sortir les grdeurs par sous surface)  
   REAL pblh(klon, nbsrf)  
   REAL plcl(klon, nbsrf)  
   REAL capcl(klon, nbsrf)  
   REAL oliqcl(klon, nbsrf)  
   REAL cteicl(klon, nbsrf)  
   REAL pblt(klon, nbsrf)  
   REAL therm(klon, nbsrf)  
   REAL trmb1(klon, nbsrf)  
   REAL trmb2(klon, nbsrf)  
   REAL trmb3(klon, nbsrf)  
   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 y_cd_h(klon), y_cd_m(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.)  
   
   REAL zt, zqs, zdelta, zcor  
   REAL t_coup  
   PARAMETER (t_coup=273.15)  
   
   CHARACTER (len=20) :: modname = 'clmain'  
   
   !------------------------------------------------------------  
   
   ! initialisation Anne  
   ytherm = 0.  
   
   IF (debugindex .AND. first_appel) THEN  
      first_appel = .FALSE.  
   
      ! initialisation sorties netcdf  
   
      idayref = day_ini  
      CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)  
      CALL gr_fi_ecrit(1, klon, iim, jjm+1, rlon, zx_lon)  
      DO i = 1, iim  
         zx_lon(i, 1) = rlon(i+1)  
         zx_lon(i, jjm+1) = rlon(i+1)  
      END DO  
      CALL gr_fi_ecrit(1, klon, iim, jjm+1, rlat, zx_lat)  
      cldebug = 'sous_index'  
      CALL histbeg_totreg(cldebug, zx_lon(:, 1), zx_lat(1, :), 1, &  
           iim, 1, jjm+1, itau_phy, zjulian, dtime, nhoridbg, nidbg)  
      ! no vertical axis  
      cl_surf(1) = 'ter'  
      cl_surf(2) = 'lic'  
      cl_surf(3) = 'oce'  
      cl_surf(4) = 'sic'  
      DO nsrf = 1, nbsrf  
         CALL histdef(nidbg, cl_surf(nsrf), cl_surf(nsrf), '-', iim, jjm+1, &  
              nhoridbg, 1, 1, 1, -99, 'inst', dtime, dtime)  
      END DO  
      CALL histend(nidbg)  
      CALL histsync(nidbg)  
   END IF  
   
   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  
   
   ! initialisation:  
   
   DO i = 1, klon  
      rugmer(i) = 0.0  
      cdragh(i) = 0.0  
      cdragm(i) = 0.0  
      dflux_t(i) = 0.0  
      dflux_q(i) = 0.0  
      zu1(i) = 0.0  
      zv1(i) = 0.0  
   END DO  
   ypct = 0.0  
   yts = 0.0  
   ysnow = 0.0  
   yqsurf = 0.0  
   yalb = 0.0  
   yalblw = 0.0  
   yrain_f = 0.0  
   ysnow_f = 0.0  
   yfder = 0.0  
   ytaux = 0.0  
   ytauy = 0.0  
   ysolsw = 0.0  
   ysollw = 0.0  
   ysollwdown = 0.0  
   yrugos = 0.0  
   yu1 = 0.0  
   yv1 = 0.0  
   yrads = 0.0  
   ypaprs = 0.0  
   ypplay = 0.0  
   ydelp = 0.0  
   yu = 0.0  
   yv = 0.0  
   yt = 0.0  
   yq = 0.0  
   pctsrf_new = 0.0  
   y_flux_u = 0.0  
   y_flux_v = 0.0  
   !$$ PB  
   y_dflux_t = 0.0  
   y_dflux_q = 0.0  
   ytsoil = 999999.  
   yrugoro = 0.  
   ! -- LOOP  
   yu10mx = 0.0  
   yu10my = 0.0  
   ywindsp = 0.0  
   ! -- LOOP  
   DO nsrf = 1, nbsrf  
      DO i = 1, klon  
         d_ts(i, nsrf) = 0.0  
      END DO  
   END DO  
   !§§§ PB  
   yfluxlat = 0.  
   flux_t = 0.  
   flux_q = 0.  
   flux_u = 0.  
   flux_v = 0.  
   DO k = 1, klev  
      DO i = 1, klon  
         d_t(i, k) = 0.0  
         d_q(i, k) = 0.0  
         !$$$         flux_t(i, k) = 0.0  
         !$$$         flux_q(i, k) = 0.0  
         d_u(i, k) = 0.0  
         d_v(i, k) = 0.0  
         !$$$         flux_u(i, k) = 0.0  
         !$$$         flux_v(i, k) = 0.0  
         zcoefh(i, k) = 0.0  
      END DO  
   END DO  
   !AA      IF (itr.GE.1) THEN  
   !AA      DO it = 1, itr  
   !AA      DO k = 1, klev  
   !AA      DO i = 1, klon  
   !AA         d_tr(i, k, it) = 0.0  
   !AA      ENDDO  
   !AA      ENDDO  
   !AA      ENDDO  
   !AA      ENDIF  
   
   
   ! Boucler sur toutes les sous-fractions du sol:  
   
   ! Initialisation des "pourcentages potentiels". On considere ici qu'on  
   ! peut avoir potentiellementdela glace sur tout le domaine oceanique  
   ! (a affiner)  
   
   pctsrf_pot = pctsrf  
   pctsrf_pot(:, is_oce) = 1. - zmasq  
   pctsrf_pot(:, is_sic) = 1. - zmasq  
   
   DO nsrf = 1, nbsrf  
      ! chercher les indices:  
      ni = 0  
      knon = 0  
      DO i = 1, klon  
         ! pour determiner 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  
   
      ! variables pour avoir une sortie IOIPSL des INDEX  
      IF (debugindex) THEN  
         tabindx = 0.  
         DO i = 1, knon  
            tabindx(i) = real(i)  
         END DO  
         debugtab = 0.  
         ndexbg = 0  
         CALL gath2cpl(tabindx, debugtab, klon, knon, iim, jjm, ni)  
         CALL histwrite(nidbg, cl_surf(nsrf), itap, debugtab)  
      END IF  
   
      IF (knon==0) CYCLE  
   
      DO j = 1, knon  
         i = ni(j)  
         ypct(j) = pctsrf(i, nsrf)  
         yts(j) = ts(i, nsrf)  
         ytslab(i) = tslab(i)  
         ysnow(j) = snow(i, nsrf)  
         yqsurf(j) = qsurf(i, nsrf)  
         yalb(j) = albe(i, nsrf)  
         yalblw(j) = alblw(i, nsrf)  
         yrain_f(j) = rain_f(i)  
         ysnow_f(j) = snow_f(i)  
         yagesno(j) = agesno(i, nsrf)  
         yfder(j) = fder(i)  
         ytaux(j) = flux_u(i, 1, nsrf)  
         ytauy(j) = flux_v(i, 1, nsrf)  
         ysolsw(j) = solsw(i, nsrf)  
         ysollw(j) = sollw(i, nsrf)  
         ysollwdown(j) = sollwdown(i)  
         yrugos(j) = rugos(i, nsrf)  
         yrugoro(j) = rugoro(i)  
         yu1(j) = u1lay(i)  
         yv1(j) = v1lay(i)  
         yrads(j) = ysolsw(j) + ysollw(j)  
         ypaprs(j, klev+1) = paprs(i, klev+1)  
         y_run_off_lic_0(j) = run_off_lic_0(i)  
         yu10mx(j) = u10m(i, nsrf)  
         yu10my(j) = v10m(i, nsrf)  
         ywindsp(j) = sqrt(yu10mx(j)*yu10mx(j)+yu10my(j)*yu10my(j))  
      END DO  
   
      !     IF bucket model for continent, copy soil water content  
      IF (nsrf==is_ter .AND. .NOT. ok_veget) THEN  
         DO j = 1, knon  
            i = ni(j)  
            yqsol(j) = qsol(i)  
         END DO  
      ELSE  
         yqsol = 0.  
      END IF  
      !$$$ PB ajour pour soil  
      DO k = 1, nsoilmx  
         DO j = 1, knon  
            i = ni(j)  
            ytsoil(j, k) = ftsoil(i, k, nsrf)  
         END DO  
      END DO  
      DO k = 1, klev  
         DO j = 1, knon  
            i = ni(j)  
            ypaprs(j, k) = paprs(i, k)  
            ypplay(j, k) = pplay(i, k)  
            ydelp(j, k) = delp(i, k)  
            yu(j, k) = u(i, k)  
            yv(j, k) = v(i, k)  
            yt(j, k) = t(i, k)  
            yq(j, k) = q(i, k)  
         END DO  
      END DO  
   
      ! calculer Cdrag et les coefficients d'echange  
      CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts,&  
           yrugos, yu, yv, yt, yq, yqsurf, ycoefm, ycoefh)  
      !IM 081204 BEG  
      !CR test  
      IF (iflag_pbl==1) THEN  
         !IM 081204 END  
         CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)  
         DO k = 1, klev  
            DO i = 1, knon  
               ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))  
               ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))  
            END DO  
         END DO  
      END IF  
   
      !IM cf JLD : on seuille ycoefm et ycoefh  
      IF (nsrf==is_oce) THEN  
         DO j = 1, knon  
            !           ycoefm(j, 1)=min(ycoefm(j, 1), 1.1E-3)  
            ycoefm(j, 1) = min(ycoefm(j, 1), cdmmax)  
            !           ycoefh(j, 1)=min(ycoefh(j, 1), 1.1E-3)  
            ycoefh(j, 1) = min(ycoefh(j, 1), cdhmax)  
         END DO  
      END IF  
   
   
      !IM: 261103  
      IF (ok_kzmin) THEN  
         !IM cf FH: 201103 BEG  
         !   Calcul d'une diffusion minimale pour les conditions tres stables.  
         CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, ycoefm, ycoefm0, &  
              ycoefh0)  
         !      call dump2d(iim, jjm-1, ycoefm(2:klon-1, 2), 'KZ         ')  
         !      call dump2d(iim, jjm-1, ycoefm0(2:klon-1, 2), 'KZMIN      ')  
   
         IF (1==1) THEN  
            DO k = 1, klev  
               DO i = 1, knon  
                  ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))  
                  ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))  
               END DO  
            END DO  
         END IF  
         !IM cf FH: 201103 END  
         !IM: 261103  
      END IF !ok_kzmin  
   
      IF (iflag_pbl>=3) THEN  
   
         !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
         ! MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin  
         !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
   
         yzlay(1:knon, 1) = rd*yt(1:knon, 1)/(0.5*(ypaprs(1:knon, &  
              1)+ypplay(1:knon, 1)))*(ypaprs(1:knon, 1)-ypplay(1: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(1:knon, klev+1) = 2.*yzlay(1:knon, klev) - yzlay(1: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  
   
   
         !   Bug introduit volontairement pour converger avec les resultats  
         !  du papier sur les thermiques.  
         IF (1==1) THEN  
            y_cd_m(1:knon) = ycoefm(1:knon, 1)  
            y_cd_h(1:knon) = ycoefh(1:knon, 1)  
         ELSE  
            y_cd_h(1:knon) = ycoefm(1:knon, 1)  
            y_cd_m(1:knon) = ycoefh(1:knon, 1)  
         END IF  
         CALL ustarhb(knon, yu, yv, y_cd_m, yustar)  
   
         IF (prt_level>9) THEN  
            PRINT *, 'USTAR = ', yustar  
         END IF  
   
         !   iflag_pbl peut etre utilise comme longuer de melange  
   
         IF (iflag_pbl>=11) THEN  
            CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, yu, yv, yteta, &  
                 y_cd_m, yq2, q2diag, ykmm, ykmn, yustar, iflag_pbl)  
         ELSE  
            CALL yamada4(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, yu, yv, yteta, &  
                 y_cd_m, yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)  
         END IF  
   
         ycoefm(1:knon, 1) = y_cd_m(1:knon)  
         ycoefh(1:knon, 1) = y_cd_h(1:knon)  
         ycoefm(1:knon, 2:klev) = ykmm(1:knon, 2:klev)  
         ycoefh(1:knon, 2:klev) = ykmn(1:knon, 2:klev)  
   
   
      END IF  
   
      !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
      ! calculer la diffusion des vitesses "u" et "v"  
      !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
   
      CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yu, ypaprs, ypplay, &  
           ydelp, y_d_u, y_flux_u)  
      CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yv, ypaprs, ypplay, &  
           ydelp, y_d_v, y_flux_v)  
   
      ! pour le couplage  
      ytaux = y_flux_u(:, 1)  
      ytauy = y_flux_v(:, 1)  
   
      ! FH modif sur le cdrag temperature  
      !$$$PB : déplace dans clcdrag  
      !$$$      do i=1, knon  
      !$$$         ycoefh(i, 1)=ycoefm(i, 1)*0.8  
      !$$$      enddo  
   
      ! calculer la diffusion de "q" et de "h"  
      CALL clqh(dtime, itap, date0, jour, debut, lafin, rlon, rlat,&  
           cufi, cvfi, knon, nsrf, ni, pctsrf, soil_model, ytsoil,&  
           yqsol, ok_veget, ocean, npas, nexca, rmu0, co2_ppm, yrugos,&  
           yrugoro, yu1, yv1, ycoefh, yt, yq, yts, ypaprs, ypplay,&  
           ydelp, yrads, yalb, yalblw, ysnow, yqsurf, yrain_f, ysnow_f, &  
           yfder, ytaux, ytauy, ywindsp, ysollw, ysollwdown, ysolsw,&  
           yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, y_d_ts,&  
           yz0_new, y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,&  
           y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, y_flux_g,&  
           ytslab, y_seaice)  
   
      ! calculer la longueur de rugosite sur ocean  
      yrugm = 0.  
      IF (nsrf==is_oce) THEN  
         DO j = 1, knon  
            yrugm(j) = 0.018*ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &  
                 0.11*14E-6/sqrt(ycoefm(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)  
            ycoefh(j, k) = ycoefh(j, k)*ypct(j)  
            ycoefm(j, k) = ycoefm(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)  
            !§§§ PB  
            flux_t(i, k, nsrf) = y_flux_t(j, k)  
            flux_q(i, k, nsrf) = y_flux_q(j, k)  
            flux_u(i, k, nsrf) = y_flux_u(j, k)  
            flux_v(i, k, nsrf) = y_flux_v(j, k)  
            !$$$ PB        y_flux_t(j, k) = y_flux_t(j, k) * ypct(j)  
            !$$$ PB        y_flux_q(j, k) = y_flux_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)  
            !$$$ PB        y_flux_u(j, k) = y_flux_u(j, k) * ypct(j)  
            !$$$ PB        y_flux_v(j, k) = y_flux_v(j, k) * ypct(j)  
         END DO  
      END DO  
   
   
      evap(:, nsrf) = -flux_q(:, 1, nsrf)  
   
      albe(:, nsrf) = 0.  
      alblw(:, 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)  
         albe(i, nsrf) = yalb(j)  
         alblw(i, nsrf) = yalblw(j)  
         snow(i, nsrf) = ysnow(j)  
         qsurf(i, nsrf) = yqsurf(j)  
         rugos(i, nsrf) = yz0_new(j)  
         fluxlat(i, nsrf) = yfluxlat(j)  
         !$$$ pb         rugmer(i) = yrugm(j)  
         IF (nsrf==is_oce) THEN  
            rugmer(i) = yrugm(j)  
            rugos(i, nsrf) = yrugm(j)  
         END IF  
         !IM cf JLD ??  
         agesno(i, nsrf) = yagesno(j)  
         fqcalving(i, nsrf) = y_fqcalving(j)  
         ffonte(i, nsrf) = y_ffonte(j)  
         cdragh(i) = cdragh(i) + ycoefh(j, 1)  
         cdragm(i) = cdragm(i) + ycoefm(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  
         DO j = 1, knon  
            i = ni(j)  
            qsol(i) = yqsol(j)  
         END DO  
      END IF  
      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  
      !$$$ PB ajout pour soil  
      ftsoil(:, :, nsrf) = 0.  
      DO k = 1, nsoilmx  
         DO j = 1, knon  
            i = ni(j)  
            ftsoil(i, k, nsrf) = ytsoil(j, k)  
         END DO  
      END DO  
   
      DO j = 1, knon  
         i = ni(j)  
         DO k = 1, klev  
            d_t(i, k) = d_t(i, k) + y_d_t(j, k)  
            d_q(i, k) = d_q(i, k) + y_d_q(j, k)  
            !$$$ PB        flux_t(i, k) = flux_t(i, k) + y_flux_t(j, k)  
            !$$$         flux_q(i, k) = flux_q(i, k) + y_flux_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)  
            !$$$  PB       flux_u(i, k) = flux_u(i, k) + y_flux_u(j, k)  
            !$$$         flux_v(i, k) = flux_v(i, k) + y_flux_v(j, k)  
            zcoefh(i, k) = zcoefh(i, k) + ycoefh(j, k)  
         END DO  
      END DO  
   
   
      !cc 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)  
      !IM 081204 END  
   
      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  
   
      DO i = 1, knon  
         y_cd_h(i) = ycoefh(i, 1)  
         y_cd_m(i) = ycoefm(i, 1)  
      END DO  
      CALL hbtm(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, &  
           y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &  
           ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)  
   
      DO j = 1, knon  
         i = ni(j)  
         pblh(i, nsrf) = ypblh(j)  
         plcl(i, nsrf) = ylcl(j)  
         capcl(i, nsrf) = ycapcl(j)  
         oliqcl(i, nsrf) = yoliqcl(j)  
         cteicl(i, nsrf) = ycteicl(j)  
         pblt(i, nsrf) = ypblt(j)  
         therm(i, nsrf) = ytherm(j)  
         trmb1(i, nsrf) = ytrmb1(j)  
         trmb2(i, nsrf) = ytrmb2(j)  
         trmb3(i, nsrf) = ytrmb3(j)  
      END DO  
   
   
      DO j = 1, knon  
         DO k = 1, klev + 1  
            i = ni(j)  
            q2(i, k, nsrf) = yq2(j, k)  
         END DO  
      END DO  
      !IM "slab" ocean  
      IF (nsrf==is_oce) THEN  
         DO j = 1, knon  
            ! on projette sur la grille globale  
            i = ni(j)  
            IF (pctsrf_new(i, is_oce)>epsfra) THEN  
               flux_o(i) = y_flux_o(j)  
            ELSE  
               flux_o(i) = 0.  
            END IF  
         END DO  
      END IF  
   
      IF (nsrf==is_sic) THEN  
         DO j = 1, knon  
            i = ni(j)  
            !IM 230604 on pondere lorsque l'on fait le bilan au sol :  flux_g(i) = y_flux_g(j)*ypct(j)  
            IF (pctsrf_new(i, is_sic)>epsfra) THEN  
               flux_g(i) = y_flux_g(j)  
            ELSE  
               flux_g(i) = 0.  
            END IF  
         END DO  
   
      END IF  
      !nsrf.EQ.is_sic                                              
      IF (ocean=='slab  ') THEN  
         IF (nsrf==is_oce) THEN  
            tslab(1:klon) = ytslab(1:klon)  
            seaice(1:klon) = y_seaice(1:klon)  
            !nsrf                                                        
         END IF  
         !OCEAN                                                        
      END IF  
   END DO  
6    
7    ! On utilise les nouvelles surfaces    SUBROUTINE pbl_surface(pctsrf, t, q, u, v, julien, mu0, ftsol, cdmmax, &
8    ! A rajouter: conservation de l'albedo         cdhmax, ftsoil, qsol, paprs, pplay, fsnow, qsurf, evap, falbe, fluxlat, &
9           rain_fall, snow_f, fsolsw, fsollw, frugs, agesno, rugoro, d_t, d_q, &
10           d_u, d_v, d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2, &
11           dflux_t, dflux_q, coefh, t2m, q2m, u10m_srf, v10m_srf, pblh, capcl, &
12           oliqcl, cteicl, pblt, therm, plcl, fqcalving, ffonte, run_off_lic_0)
13    
14        ! From phylmd/clmain.F, version 1.6, 2005/11/16 14:47:19
15        ! Author: Z. X. Li (LMD/CNRS)
16        ! Date: Aug. 18th, 1993
17        ! Objet : interface de couche limite (diffusion verticale)
18    
19        ! Tout ce qui a trait aux traceurs est dans "phytrac". Le calcul
20        ! de la couche limite pour les traceurs se fait avec "cltrac" et
21        ! ne tient pas compte de la diff\'erentiation des sous-fractions
22        ! de sol.
23    
24        use cdrag_m, only: cdrag
25        use clqh_m, only: clqh
26        use clvent_m, only: clvent
27        use coef_diff_turb_m, only: coef_diff_turb
28        USE conf_gcm_m, ONLY: lmt_pas
29        USE conf_phys_m, ONLY: iflag_pbl
30        USE dimphy, ONLY: klev, klon
31        USE dimsoil, ONLY: nsoilmx
32        use hbtm_m, only: hbtm
33        USE histwrite_phy_m, ONLY: histwrite_phy
34        USE indicesol, ONLY: epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
35        USE interfoce_lim_m, ONLY: interfoce_lim
36        use phyetat0_m, only: zmasq
37        use stdlevvar_m, only: stdlevvar
38        USE suphec_m, ONLY: rd, rg
39        use time_phylmdz, only: itap
40    
41        REAL, INTENT(inout):: pctsrf(klon, nbsrf)
42        ! tableau des pourcentages de surface de chaque maille
43    
44        REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
45        REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg / kg)
46        REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
47        INTEGER, INTENT(IN):: julien ! jour de l'annee en cours
48        REAL, intent(in):: mu0(klon) ! cosinus de l'angle solaire zenithal    
49        REAL, INTENT(IN):: ftsol(:, :) ! (klon, nbsrf) temp\'erature du sol (en K)
50        REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh
51    
52        REAL, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)
53        ! soil temperature of surface fraction
54    
55        REAL, INTENT(inout):: qsol(:) ! (klon)
56        ! column-density of water in soil, in kg m-2
57    
58        REAL, INTENT(IN):: paprs(klon, klev + 1) ! pression a intercouche (Pa)
59        REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
60        REAL, INTENT(inout):: fsnow(:, :) ! (klon, nbsrf) \'epaisseur neigeuse
61        REAL qsurf(klon, nbsrf)
62        REAL evap(klon, nbsrf)
63        REAL, intent(inout):: falbe(klon, nbsrf)
64        REAL, intent(out):: fluxlat(:, :) ! (klon, nbsrf)
65    
66        REAL, intent(in):: rain_fall(klon)
67        ! liquid water mass flux (kg / m2 / s), positive down
68    
69        REAL, intent(in):: snow_f(klon)
70        ! solid water mass flux (kg / m2 / s), positive down
71    
72        REAL, INTENT(IN):: fsolsw(klon, nbsrf), fsollw(klon, nbsrf)
73        REAL, intent(inout):: frugs(klon, nbsrf) ! longueur de rugosit\'e (en m)
74        real agesno(klon, nbsrf)
75        REAL, INTENT(IN):: rugoro(klon)
76    
77        REAL, intent(out):: d_t(:, :), d_q(:, :) ! (klon, klev)
78        ! changement pour t et q
79    
80        REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
81        ! changement pour "u" et "v"
82    
83        REAL, intent(out):: d_ts(:, :) ! (klon, nbsrf) variation of ftsol
84    
85        REAL, intent(out):: flux_t(klon, nbsrf)
86        ! flux de chaleur sensible (c_p T) (W / m2) (orientation positive
87        ! vers le bas) à la surface
88    
89        REAL, intent(out):: flux_q(klon, nbsrf)
90        ! flux de vapeur d'eau (kg / m2 / s) à la surface
91    
92        REAL, intent(out):: flux_u(klon, nbsrf), flux_v(klon, nbsrf)
93        ! tension du vent (flux turbulent de vent) à la surface, en Pa
94    
95        REAL, INTENT(out):: cdragh(klon), cdragm(klon)
96        real q2(klon, klev + 1, nbsrf)
97    
98        ! Ocean slab:
99        REAL, INTENT(out):: dflux_t(klon) ! derive du flux sensible
100        REAL, INTENT(out):: dflux_q(klon) ! derive du flux latent
101    
102        REAL, intent(out):: coefh(:, 2:) ! (klon, 2:klev)
103        ! Pour pouvoir extraire les coefficients d'\'echange, le champ
104        ! "coefh" a \'et\'e cr\'e\'e. Nous avons moyenn\'e les valeurs de
105        ! ce champ sur les quatre sous-surfaces du mod\`ele.
106    
107        REAL, INTENT(inout):: t2m(klon, nbsrf), q2m(klon, nbsrf)
108    
109        REAL, INTENT(inout):: u10m_srf(:, :), v10m_srf(:, :) ! (klon, nbsrf)
110        ! composantes du vent \`a 10m sans spirale d'Ekman
111    
112        ! Ionela Musat. Cf. Anne Mathieu : planetary boundary layer, hbtm.
113        ! Comme les autres diagnostics on cumule dans physiq ce qui permet
114        ! de sortir les grandeurs par sous-surface.
115        REAL pblh(klon, nbsrf) ! height of planetary boundary layer
116        REAL capcl(klon, nbsrf)
117        REAL oliqcl(klon, nbsrf)
118        REAL cteicl(klon, nbsrf)
119        REAL, INTENT(inout):: pblt(klon, nbsrf) ! T au nveau HCL
120        REAL therm(klon, nbsrf)
121        REAL plcl(klon, nbsrf)
122    
123        REAL, intent(out):: fqcalving(klon, nbsrf)
124        ! flux d'eau "perdue" par la surface et necessaire pour limiter la
125        ! hauteur de neige, en kg / m2 / s
126    
127        real ffonte(klon, nbsrf) ! flux thermique utilise pour fondre la neige
128        REAL, intent(inout):: run_off_lic_0(:) ! (klon)
129    
130        ! Local:
131    
132        ! la nouvelle repartition des surfaces sortie de l'interface
133        REAL, save:: pctsrf_new_oce(klon)
134        REAL, save:: pctsrf_new_sic(klon)
135    
136        REAL y_fqcalving(klon), y_ffonte(klon)
137        real y_run_off_lic_0(klon), y_run_off_lic(klon)
138        REAL run_off_lic(klon) ! ruissellement total
139        REAL rugmer(klon)
140        REAL ytsoil(klon, nsoilmx)
141        REAL yts(klon), ypct(klon), yz0_new(klon)
142        real yrugos(klon) ! longueur de rugosite (en m)
143        REAL yalb(klon)
144        REAL snow(klon), yqsurf(klon), yagesno(klon)
145        real yqsol(klon) ! column-density of water in soil, in kg m-2
146        REAL yrain_f(klon) ! liquid water mass flux (kg / m2 / s), positive down
147        REAL ysnow_f(klon) ! solid water mass flux (kg / m2 / s), positive down
148        REAL yrugm(klon), yrads(klon), yrugoro(klon)
149        REAL yfluxlat(klon)
150        REAL y_d_ts(klon)
151        REAL y_d_t(klon, klev), y_d_q(klon, klev)
152        REAL y_d_u(klon, klev), y_d_v(klon, klev)
153        REAL y_flux_t(klon), y_flux_q(klon)
154        REAL y_flux_u(klon), y_flux_v(klon)
155        REAL y_dflux_t(klon), y_dflux_q(klon)
156        REAL ycoefh(klon, 2:klev), ycoefm(klon, 2:klev)
157        real ycdragh(klon), ycdragm(klon)
158        REAL yu(klon, klev), yv(klon, klev)
159        REAL yt(klon, klev), yq(klon, klev)
160        REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)
161        REAL yq2(klon, klev + 1)
162        REAL delp(klon, klev)
163        INTEGER i, k, nsrf
164        INTEGER ni(klon), knon, j
165    
166        REAL pctsrf_pot(klon, nbsrf)
167        ! "pourcentage potentiel" pour tenir compte des \'eventuelles
168        ! apparitions ou disparitions de la glace de mer
169    
170        REAL yt2m(klon), yq2m(klon), wind10m(klon)
171        REAL ustar(klon)
172    
173        REAL yt10m(klon), yq10m(klon)
174        REAL ypblh(klon)
175        REAL ylcl(klon)
176        REAL ycapcl(klon)
177        REAL yoliqcl(klon)
178        REAL ycteicl(klon)
179        REAL ypblt(klon)
180        REAL ytherm(klon)
181        REAL u1(klon), v1(klon)
182        REAL tair1(klon), qair1(klon), tairsol(klon)
183        REAL psfce(klon), patm(klon)
184    
185        REAL qairsol(klon), zgeo1(klon)
186        REAL rugo1(klon)
187        REAL zgeop(klon, klev)
188    
189        !------------------------------------------------------------
190    
191        ytherm = 0.
192    
193        DO k = 1, klev ! epaisseur de couche
194           DO i = 1, klon
195              delp(i, k) = paprs(i, k) - paprs(i, k + 1)
196           END DO
197        END DO
198    
199        ! Initialization:
200        rugmer = 0.
201        cdragh = 0.
202        cdragm = 0.
203        dflux_t = 0.
204        dflux_q = 0.
205        ypct = 0.
206        yqsurf = 0.
207        yrain_f = 0.
208        ysnow_f = 0.
209        yrugos = 0.
210        ypaprs = 0.
211        ypplay = 0.
212        ydelp = 0.
213        yrugoro = 0.
214        d_ts = 0.
215        flux_t = 0.
216        flux_q = 0.
217        flux_u = 0.
218        flux_v = 0.
219        fluxlat = 0.
220        d_t = 0.
221        d_q = 0.
222        d_u = 0.
223        d_v = 0.
224        coefh = 0.
225        fqcalving = 0.
226        run_off_lic = 0.
227    
228        ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
229        ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
230        ! (\`a affiner).
231    
232        pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
233        pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
234        pctsrf_pot(:, is_oce) = 1. - zmasq
235        pctsrf_pot(:, is_sic) = 1. - zmasq
236    
237        ! Tester si c'est le moment de lire le fichier:
238        if (mod(itap - 1, lmt_pas) == 0) then
239           CALL interfoce_lim(julien, pctsrf_new_oce, pctsrf_new_sic)
240        endif
241    
242        ! Boucler sur toutes les sous-fractions du sol:
243    
244        loop_surface: DO nsrf = 1, nbsrf
245           ! Define ni and knon:
246          
247           ni = 0
248           knon = 0
249    
250           DO i = 1, klon
251              ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
252              ! "potentielles"
253              IF (pctsrf_pot(i, nsrf) > epsfra) THEN
254                 knon = knon + 1
255                 ni(knon) = i
256              END IF
257           END DO
258    
259           if_knon: IF (knon /= 0) then
260              DO j = 1, knon
261                 i = ni(j)
262                 ypct(j) = pctsrf(i, nsrf)
263                 yts(j) = ftsol(i, nsrf)
264                 snow(j) = fsnow(i, nsrf)
265                 yqsurf(j) = qsurf(i, nsrf)
266                 yalb(j) = falbe(i, nsrf)
267                 yrain_f(j) = rain_fall(i)
268                 ysnow_f(j) = snow_f(i)
269                 yagesno(j) = agesno(i, nsrf)
270                 yrugos(j) = frugs(i, nsrf)
271                 yrugoro(j) = rugoro(i)
272                 yrads(j) = fsolsw(i, nsrf) + fsollw(i, nsrf)
273                 ypaprs(j, klev + 1) = paprs(i, klev + 1)
274                 y_run_off_lic_0(j) = run_off_lic_0(i)
275              END DO
276    
277              ! For continent, copy soil water content
278              IF (nsrf == is_ter) yqsol(:knon) = qsol(ni(:knon))
279    
280              ytsoil(:knon, :) = ftsoil(ni(:knon), :, nsrf)
281    
282              DO k = 1, klev
283                 DO j = 1, knon
284                    i = ni(j)
285                    ypaprs(j, k) = paprs(i, k)
286                    ypplay(j, k) = pplay(i, k)
287                    ydelp(j, k) = delp(i, k)
288                    yu(j, k) = u(i, k)
289                    yv(j, k) = v(i, k)
290                    yt(j, k) = t(i, k)
291                    yq(j, k) = q(i, k)
292                 END DO
293              END DO
294    
295              ! Calculer les géopotentiels de chaque couche:
296    
297              zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
298                   + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
299    
300              DO k = 2, klev
301                 zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
302                      * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
303                      * (ypplay(:knon, k - 1) - ypplay(:knon, k))
304              ENDDO
305    
306              CALL cdrag(nsrf, sqrt(yu(:knon, 1)**2 + yv(:knon, 1)**2), &
307                   yt(:knon, 1), yq(:knon, 1), zgeop(:knon, 1), ypaprs(:knon, 1), &
308                   yts(:knon), yqsurf(:knon), yrugos(:knon), ycdragm(:knon), &
309                   ycdragh(:knon))
310    
311              IF (iflag_pbl == 1) THEN
312                 ycdragm(:knon) = max(ycdragm(:knon), 0.)
313                 ycdragh(:knon) = max(ycdragh(:knon), 0.)
314              end IF
315    
316              ! on met un seuil pour ycdragm et ycdragh
317              IF (nsrf == is_oce) THEN
318                 ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
319                 ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
320              END IF
321    
322              IF (iflag_pbl >= 6) yq2(:knon, :) = q2(ni(:knon), :, nsrf)
323              call coef_diff_turb(nsrf, ni(:knon), ypaprs(:knon, :), &
324                   ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
325                   yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
326                   ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))
327              
328              CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
329                   ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
330                   ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
331                   y_flux_u(:knon))
332              CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
333                   ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
334                   ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
335                   y_flux_v(:knon))
336    
337              CALL clqh(julien, nsrf, ni(:knon), ytsoil(:knon, :), yqsol(:knon), &
338                   mu0(ni(:knon)), yrugos(:knon), yrugoro(:knon), yu(:knon, 1), &
339                   yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), yt(:knon, :), &
340                   yq(:knon, :), yts(:knon), ypaprs(:knon, :), ypplay(:knon, :), &
341                   ydelp(:knon, :), yrads(:knon), yalb(:knon), snow(:knon), &
342                   yqsurf(:knon), yrain_f(:knon), ysnow_f(:knon), yfluxlat(:knon), &
343                   pctsrf_new_sic(ni(:knon)), yagesno(:knon), y_d_t(:knon, :), &
344                   y_d_q(:knon, :), y_d_ts(:knon), yz0_new(:knon), &
345                   y_flux_t(:knon), y_flux_q(:knon), y_dflux_t(:knon), &
346                   y_dflux_q(:knon), y_fqcalving(:knon), y_ffonte(:knon), &
347                   y_run_off_lic_0(:knon), y_run_off_lic(:knon))
348    
349              ! calculer la longueur de rugosite sur ocean
350    
351              yrugm = 0.
352    
353              IF (nsrf == is_oce) THEN
354                 DO j = 1, knon
355                    yrugm(j) = 0.018 * ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2) &
356                         / rg + 0.11 * 14E-6 &
357                         / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2))
358                    yrugm(j) = max(1.5E-05, yrugm(j))
359                 END DO
360              END IF
361    
362              DO k = 1, klev
363                 DO j = 1, knon
364                    i = ni(j)
365                    y_d_t(j, k) = y_d_t(j, k) * ypct(j)
366                    y_d_q(j, k) = y_d_q(j, k) * ypct(j)
367                    y_d_u(j, k) = y_d_u(j, k) * ypct(j)
368                    y_d_v(j, k) = y_d_v(j, k) * ypct(j)
369                 END DO
370              END DO
371    
372              flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
373              flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
374              flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
375              flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
376    
377              evap(:, nsrf) = -flux_q(:, nsrf)
378    
379              falbe(:, nsrf) = 0.
380              fsnow(:, nsrf) = 0.
381              qsurf(:, nsrf) = 0.
382              frugs(:, nsrf) = 0.
383              DO j = 1, knon
384                 i = ni(j)
385                 d_ts(i, nsrf) = y_d_ts(j)
386                 falbe(i, nsrf) = yalb(j)
387                 fsnow(i, nsrf) = snow(j)
388                 qsurf(i, nsrf) = yqsurf(j)
389                 frugs(i, nsrf) = yz0_new(j)
390                 fluxlat(i, nsrf) = yfluxlat(j)
391                 IF (nsrf == is_oce) THEN
392                    rugmer(i) = yrugm(j)
393                    frugs(i, nsrf) = yrugm(j)
394                 END IF
395                 agesno(i, nsrf) = yagesno(j)
396                 fqcalving(i, nsrf) = y_fqcalving(j)
397                 ffonte(i, nsrf) = y_ffonte(j)
398                 cdragh(i) = cdragh(i) + ycdragh(j) * ypct(j)
399                 cdragm(i) = cdragm(i) + ycdragm(j) * ypct(j)
400                 dflux_t(i) = dflux_t(i) + y_dflux_t(j) * ypct(j)
401                 dflux_q(i) = dflux_q(i) + y_dflux_q(j) * ypct(j)
402              END DO
403              IF (nsrf == is_ter) THEN
404                 qsol(ni(:knon)) = yqsol(:knon)
405              else IF (nsrf == is_lic) THEN
406                 DO j = 1, knon
407                    i = ni(j)
408                    run_off_lic_0(i) = y_run_off_lic_0(j)
409                    run_off_lic(i) = y_run_off_lic(j)
410                 END DO
411              END IF
412    
413              ftsoil(:, :, nsrf) = 0.
414              ftsoil(ni(:knon), :, nsrf) = ytsoil(:knon, :)
415    
416              DO j = 1, knon
417                 i = ni(j)
418                 DO k = 1, klev
419                    d_t(i, k) = d_t(i, k) + y_d_t(j, k)
420                    d_q(i, k) = d_q(i, k) + y_d_q(j, k)
421                    d_u(i, k) = d_u(i, k) + y_d_u(j, k)
422                    d_v(i, k) = d_v(i, k) + y_d_v(j, k)
423                 END DO
424              END DO
425    
426              forall (k = 2:klev) coefh(ni(:knon), k) &
427                   = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypct(:knon)
428    
429              ! diagnostic t, q a 2m et u, v a 10m
430    
431              DO j = 1, knon
432                 i = ni(j)
433                 u1(j) = yu(j, 1) + y_d_u(j, 1)
434                 v1(j) = yv(j, 1) + y_d_v(j, 1)
435                 tair1(j) = yt(j, 1) + y_d_t(j, 1)
436                 qair1(j) = yq(j, 1) + y_d_q(j, 1)
437                 zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
438                      1))) * (ypaprs(j, 1)-ypplay(j, 1))
439                 tairsol(j) = yts(j) + y_d_ts(j)
440                 rugo1(j) = yrugos(j)
441                 IF (nsrf == is_oce) THEN
442                    rugo1(j) = frugs(i, nsrf)
443                 END IF
444                 psfce(j) = ypaprs(j, 1)
445                 patm(j) = ypplay(j, 1)
446    
447                 qairsol(j) = yqsurf(j)
448              END DO
449    
450              CALL stdlevvar(nsrf, u1(:knon), v1(:knon), tair1(:knon), qair1, &
451                   zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, &
452                   yq10m, wind10m(:knon), ustar(:knon))
453    
454              DO j = 1, knon
455                 i = ni(j)
456                 t2m(i, nsrf) = yt2m(j)
457                 q2m(i, nsrf) = yq2m(j)
458    
459                 u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
460                      / sqrt(u1(j)**2 + v1(j)**2)
461                 v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
462                      / sqrt(u1(j)**2 + v1(j)**2)
463              END DO
464    
465              CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
466                   y_flux_q(:knon), yu(:knon, :), yv(:knon, :), yt(:knon, :), &
467                   yq(:knon, :), ypblh(:knon), ycapcl, yoliqcl, ycteicl, ypblt, &
468                   ytherm, ylcl)
469    
470              DO j = 1, knon
471                 i = ni(j)
472                 pblh(i, nsrf) = ypblh(j)
473                 plcl(i, nsrf) = ylcl(j)
474                 capcl(i, nsrf) = ycapcl(j)
475                 oliqcl(i, nsrf) = yoliqcl(j)
476                 cteicl(i, nsrf) = ycteicl(j)
477                 pblt(i, nsrf) = ypblt(j)
478                 therm(i, nsrf) = ytherm(j)
479              END DO
480    
481              IF (iflag_pbl >= 6) q2(ni(:knon), :, nsrf) = yq2(:knon, :)
482           else
483              fsnow(:, nsrf) = 0.
484           end IF if_knon
485        END DO loop_surface
486    
487        ! On utilise les nouvelles surfaces
488        frugs(:, is_oce) = rugmer
489        pctsrf(:, is_oce) = pctsrf_new_oce
490        pctsrf(:, is_sic) = pctsrf_new_sic
491    
492    rugos(:, is_oce) = rugmer      CALL histwrite_phy("run_off_lic", run_off_lic)
   pctsrf = pctsrf_new  
493    
494  END SUBROUTINE clmain    END SUBROUTINE pbl_surface
495    
496    end module pbl_surface_m

Legend:
Removed from v.37  
changed lines
  Added in v.303

  ViewVC Help
Powered by ViewVC 1.1.21