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

Diff of /trunk/phylmd/pbl_surface.f

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

trunk/libf/phylmd/clmain.f90 revision 15 by guez, Fri Aug 1 15:24:12 2008 UTC trunk/Sources/phylmd/clmain.f revision 250 by guez, Fri Jan 5 18:18:53 2018 UTC
# Line 1  Line 1 
1  SUBROUTINE clmain(dtime, itap, date0, pctsrf, pctsrf_new, t, q, u, v,&  module clmain_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, & !IM "slab" ocean  
      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 ioipsl  
   USE interface_surf  
   USE dimens_m  
   USE indicesol  
   USE dimphy  
   USE dimsoil  
   USE temps  
   USE iniprint  
   USE yomcst  
   USE yoethf  
   USE fcttre  
   USE conf_phys_m  
   USE gath_cpl, ONLY : gath2cpl  
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, hbtm2 (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'  
   LOGICAL check  
   PARAMETER (check=.FALSE.)  
   
   !------------------------------------------------------------  
   
   ! initialisation Anne  
   ytherm = 0.  
   
   IF (check) THEN  
      print *, modname, '  klon=', klon  
   END IF  
   
   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  
   
      IF (check) THEN  
         print *, 'CLMAIN, nsrf, knon =', nsrf, knon  
      END IF  
   
      ! 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, & !IM 261103  
           ksta, ksta_ter, & !IM 261103  
           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  
   
      !IM cf AM : pbl, HBTM  
      DO i = 1, knon  
         y_cd_h(i) = ycoefh(i, 1)  
         y_cd_m(i) = ycoefm(i, 1)  
      END DO  
      !     print*, 'appel hbtm2'  
      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)  
      !     print*, 'fin hbtm2'  
   
      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 clmain(dtime, 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, trmb1, trmb2, trmb3, plcl, fqcalving, &
13           ffonte, run_off_lic_0)
14    
15        ! From phylmd/clmain.F, version 1.6, 2005/11/16 14:47:19
16        ! Author: Z. X. Li (LMD/CNRS), date: 1993/08/18
17        ! Objet : interface de couche limite (diffusion verticale)
18    
19        ! Tout ce qui a trait aux traceurs est dans "phytrac". Le calcul
20        ! de la couche limite pour les traceurs se fait avec "cltrac" et
21        ! ne tient pas compte de la diff\'erentiation des sous-fractions
22        ! de sol.
23    
24        use clcdrag_m, only: clcdrag
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, zmasq
31        USE dimsoil, ONLY: nsoilmx
32        use hbtm_m, only: hbtm
33        USE indicesol, ONLY: epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
34        USE interfoce_lim_m, ONLY: interfoce_lim
35        use stdlevvar_m, only: stdlevvar
36        USE suphec_m, ONLY: rd, rg
37        use time_phylmdz, only: itap
38    
39        REAL, INTENT(IN):: dtime ! interval du temps (secondes)
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 d_t(klon, klev), d_q(klon, klev)
78        ! d_t------output-R- le changement pour "t"
79        ! d_q------output-R- le changement pour "q"
80    
81        REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
82        ! changement pour "u" et "v"
83    
84        REAL, intent(out):: d_ts(:, :) ! (klon, nbsrf) variation of ftsol
85    
86        REAL, intent(out):: flux_t(klon, nbsrf)
87        ! flux de chaleur sensible (Cp T) (W / m2) (orientation positive vers
88        ! le bas) à la surface
89    
90        REAL, intent(out):: flux_q(klon, nbsrf)
91        ! flux de vapeur d'eau (kg / m2 / s) à la surface
92    
93        REAL, intent(out):: flux_u(klon, nbsrf), flux_v(klon, nbsrf)
94        ! tension du vent (flux turbulent de vent) à la surface, en Pa
95    
96        REAL, INTENT(out):: cdragh(klon), cdragm(klon)
97        real q2(klon, klev + 1, nbsrf)
98    
99        REAL, INTENT(out):: dflux_t(klon), dflux_q(klon)
100        ! dflux_t derive du flux sensible
101        ! dflux_q derive du flux latent
102        ! IM "slab" ocean
103    
104        REAL, intent(out):: coefh(:, 2:) ! (klon, 2:klev)
105        ! Pour pouvoir extraire les coefficients d'\'echange, le champ
106        ! "coefh" a \'et\'e cr\'e\'e. Nous avons moyenn\'e les valeurs de
107        ! ce champ sur les quatre sous-surfaces du mod\`ele.
108    
109        REAL, INTENT(inout):: t2m(klon, nbsrf), q2m(klon, nbsrf)
110    
111        REAL, INTENT(inout):: u10m_srf(:, :), v10m_srf(:, :) ! (klon, nbsrf)
112        ! composantes du vent \`a 10m sans spirale d'Ekman
113    
114        ! Ionela Musat. Cf. Anne Mathieu : planetary boundary layer, hbtm.
115        ! Comme les autres diagnostics on cumule dans physiq ce qui permet
116        ! de sortir les grandeurs par sous-surface.
117        REAL pblh(klon, nbsrf) ! height of planetary boundary layer
118        REAL capcl(klon, nbsrf)
119        REAL oliqcl(klon, nbsrf)
120        REAL cteicl(klon, nbsrf)
121        REAL, INTENT(inout):: pblt(klon, nbsrf) ! T au nveau HCL
122        REAL therm(klon, nbsrf)
123        REAL trmb1(klon, nbsrf)
124        ! trmb1-------deep_cape
125        REAL trmb2(klon, nbsrf)
126        ! trmb2--------inhibition
127        REAL trmb3(klon, nbsrf)
128        ! trmb3-------Point Omega
129        REAL plcl(klon, nbsrf)
130        REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)
131        ! ffonte----Flux thermique utilise pour fondre la neige
132        ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
133        !           hauteur de neige, en kg / m2 / s
134        REAL run_off_lic_0(klon)
135    
136        ! Local:
137    
138        LOGICAL:: firstcal = .true.
139    
140        ! la nouvelle repartition des surfaces sortie de l'interface
141        REAL, save:: pctsrf_new_oce(klon)
142        REAL, save:: pctsrf_new_sic(klon)
143    
144        REAL y_fqcalving(klon), y_ffonte(klon)
145        real y_run_off_lic_0(klon)
146        REAL rugmer(klon)
147        REAL ytsoil(klon, nsoilmx)
148        REAL yts(klon), ypct(klon), yz0_new(klon)
149        real yrugos(klon) ! longeur de rugosite (en m)
150        REAL yalb(klon)
151        REAL snow(klon), yqsurf(klon), yagesno(klon)
152        real yqsol(klon) ! column-density of water in soil, in kg m-2
153        REAL yrain_f(klon) ! liquid water mass flux (kg / m2 / s), positive down
154        REAL ysnow_f(klon) ! solid water mass flux (kg / m2 / s), positive down
155        REAL yrugm(klon), yrads(klon), yrugoro(klon)
156        REAL yfluxlat(klon)
157        REAL y_d_ts(klon)
158        REAL y_d_t(klon, klev), y_d_q(klon, klev)
159        REAL y_d_u(klon, klev), y_d_v(klon, klev)
160        REAL y_flux_t(klon), y_flux_q(klon)
161        REAL y_flux_u(klon), y_flux_v(klon)
162        REAL y_dflux_t(klon), y_dflux_q(klon)
163        REAL ycoefh(klon, 2:klev), ycoefm(klon, 2:klev)
164        real ycdragh(klon), ycdragm(klon)
165        REAL yu(klon, klev), yv(klon, klev)
166        REAL yt(klon, klev), yq(klon, klev)
167        REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)
168        REAL yq2(klon, klev + 1)
169        REAL delp(klon, klev)
170        INTEGER i, k, nsrf
171        INTEGER ni(klon), knon, j
172    
173        REAL pctsrf_pot(klon, nbsrf)
174        ! "pourcentage potentiel" pour tenir compte des \'eventuelles
175        ! apparitions ou disparitions de la glace de mer
176    
177        REAL yt2m(klon), yq2m(klon), wind10m(klon)
178        REAL ustar(klon)
179    
180        REAL yt10m(klon), yq10m(klon)
181        REAL ypblh(klon)
182        REAL ylcl(klon)
183        REAL ycapcl(klon)
184        REAL yoliqcl(klon)
185        REAL ycteicl(klon)
186        REAL ypblt(klon)
187        REAL ytherm(klon)
188        REAL ytrmb1(klon)
189        REAL ytrmb2(klon)
190        REAL ytrmb3(klon)
191        REAL u1(klon), v1(klon)
192        REAL tair1(klon), qair1(klon), tairsol(klon)
193        REAL psfce(klon), patm(klon)
194    
195        REAL qairsol(klon), zgeo1(klon)
196        REAL rugo1(klon)
197        REAL zgeop(klon, klev)
198    
199        !------------------------------------------------------------
200    
201        ytherm = 0.
202    
203        DO k = 1, klev ! epaisseur de couche
204           DO i = 1, klon
205              delp(i, k) = paprs(i, k) - paprs(i, k + 1)
206           END DO
207        END DO
208    
209        ! Initialization:
210        rugmer = 0.
211        cdragh = 0.
212        cdragm = 0.
213        dflux_t = 0.
214        dflux_q = 0.
215        ypct = 0.
216        yqsurf = 0.
217        yrain_f = 0.
218        ysnow_f = 0.
219        yrugos = 0.
220        ypaprs = 0.
221        ypplay = 0.
222        ydelp = 0.
223        yu = 0.
224        yv = 0.
225        yt = 0.
226        yq = 0.
227        y_dflux_t = 0.
228        y_dflux_q = 0.
229        yrugoro = 0.
230        d_ts = 0.
231        flux_t = 0.
232        flux_q = 0.
233        flux_u = 0.
234        flux_v = 0.
235        fluxlat = 0.
236        d_t = 0.
237        d_q = 0.
238        d_u = 0.
239        d_v = 0.
240        coefh = 0.
241    
242        ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
243        ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
244        ! (\`a affiner)
245    
246        pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
247        pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
248        pctsrf_pot(:, is_oce) = 1. - zmasq
249        pctsrf_pot(:, is_sic) = 1. - zmasq
250    
251        ! Tester si c'est le moment de lire le fichier:
252        if (mod(itap - 1, lmt_pas) == 0) then
253           CALL interfoce_lim(julien, pctsrf_new_oce, pctsrf_new_sic)
254        endif
255    
256        ! Boucler sur toutes les sous-fractions du sol:
257    
258        loop_surface: DO nsrf = 1, nbsrf
259           ! Chercher les indices :
260           ni = 0
261           knon = 0
262           DO i = 1, klon
263              ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
264              ! "potentielles"
265              IF (pctsrf_pot(i, nsrf) > epsfra) THEN
266                 knon = knon + 1
267                 ni(knon) = i
268              END IF
269           END DO
270    
271           if_knon: IF (knon /= 0) then
272              DO j = 1, knon
273                 i = ni(j)
274                 ypct(j) = pctsrf(i, nsrf)
275                 yts(j) = ftsol(i, nsrf)
276                 snow(j) = fsnow(i, nsrf)
277                 yqsurf(j) = qsurf(i, nsrf)
278                 yalb(j) = falbe(i, nsrf)
279                 yrain_f(j) = rain_fall(i)
280                 ysnow_f(j) = snow_f(i)
281                 yagesno(j) = agesno(i, nsrf)
282                 yrugos(j) = frugs(i, nsrf)
283                 yrugoro(j) = rugoro(i)
284                 yrads(j) = fsolsw(i, nsrf) + fsollw(i, nsrf)
285                 ypaprs(j, klev + 1) = paprs(i, klev + 1)
286                 y_run_off_lic_0(j) = run_off_lic_0(i)
287              END DO
288    
289              ! For continent, copy soil water content
290              IF (nsrf == is_ter) yqsol(:knon) = qsol(ni(:knon))
291    
292              ytsoil(:knon, :) = ftsoil(ni(:knon), :, nsrf)
293    
294              DO k = 1, klev
295                 DO j = 1, knon
296                    i = ni(j)
297                    ypaprs(j, k) = paprs(i, k)
298                    ypplay(j, k) = pplay(i, k)
299                    ydelp(j, k) = delp(i, k)
300                    yu(j, k) = u(i, k)
301                    yv(j, k) = v(i, k)
302                    yt(j, k) = t(i, k)
303                    yq(j, k) = q(i, k)
304                 END DO
305              END DO
306    
307              ! Calculer les géopotentiels de chaque couche:
308    
309              zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
310                   + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
311    
312              DO k = 2, klev
313                 zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
314                      * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
315                      * (ypplay(:knon, k - 1) - ypplay(:knon, k))
316              ENDDO
317    
318              CALL clcdrag(nsrf, yu(:knon, 1), yv(:knon, 1), yt(:knon, 1), &
319                   yq(:knon, 1), zgeop(:knon, 1), yts(:knon), yqsurf(:knon), &
320                   yrugos(:knon), ycdragm(:knon), ycdragh(:knon))
321    
322              IF (iflag_pbl == 1) THEN
323                 ycdragm(:knon) = max(ycdragm(:knon), 0.)
324                 ycdragh(:knon) = max(ycdragh(:knon), 0.)
325              end IF
326    
327              ! on met un seuil pour ycdragm et ycdragh
328              IF (nsrf == is_oce) THEN
329                 ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
330                 ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
331              END IF
332    
333              IF (iflag_pbl >= 6) then
334                 DO k = 1, klev + 1
335                    DO j = 1, knon
336                       i = ni(j)
337                       yq2(j, k) = q2(i, k, nsrf)
338                    END DO
339                 END DO
340              end IF
341    
342              call coef_diff_turb(dtime, nsrf, ni(:knon), ypaprs, ypplay, yu, yv, &
343                   yq, yt, yts, ycdragm, zgeop(:knon, :), ycoefm(:knon, :), &
344                   ycoefh(:knon, :), yq2)
345    
346              CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
347                   ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
348                   ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
349                   y_flux_u(:knon))
350              CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
351                   ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
352                   ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
353                   y_flux_v(:knon))
354    
355              ! calculer la diffusion de "q" et de "h"
356              CALL clqh(dtime, julien, firstcal, nsrf, ni(:knon), &
357                   ytsoil(:knon, :), yqsol(:knon), mu0, yrugos, yrugoro, &
358                   yu(:knon, 1), yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), &
359                   yt, yq, yts(:knon), ypaprs, ypplay, ydelp, yrads(:knon), &
360                   yalb(:knon), snow(:knon), yqsurf, yrain_f, ysnow_f, &
361                   yfluxlat(:knon), pctsrf_new_sic, yagesno(:knon), y_d_t, y_d_q, &
362                   y_d_ts(:knon), yz0_new, y_flux_t(:knon), y_flux_q(:knon), &
363                   y_dflux_t(:knon), y_dflux_q(:knon), y_fqcalving, y_ffonte, &
364                   y_run_off_lic_0)
365    
366              ! calculer la longueur de rugosite sur ocean
367              yrugm = 0.
368              IF (nsrf == is_oce) THEN
369                 DO j = 1, knon
370                    yrugm(j) = 0.018 * ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2) &
371                         / rg + 0.11 * 14E-6 &
372                         / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2))
373                    yrugm(j) = max(1.5E-05, yrugm(j))
374                 END DO
375              END IF
376              DO j = 1, knon
377                 y_dflux_t(j) = y_dflux_t(j) * ypct(j)
378                 y_dflux_q(j) = y_dflux_q(j) * ypct(j)
379              END DO
380    
381              DO k = 1, klev
382                 DO j = 1, knon
383                    i = ni(j)
384                    y_d_t(j, k) = y_d_t(j, k) * ypct(j)
385                    y_d_q(j, k) = y_d_q(j, k) * ypct(j)
386                    y_d_u(j, k) = y_d_u(j, k) * ypct(j)
387                    y_d_v(j, k) = y_d_v(j, k) * ypct(j)
388                 END DO
389              END DO
390    
391              flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
392              flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
393              flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
394              flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
395    
396              evap(:, nsrf) = -flux_q(:, nsrf)
397    
398              falbe(:, nsrf) = 0.
399              fsnow(:, nsrf) = 0.
400              qsurf(:, nsrf) = 0.
401              frugs(:, nsrf) = 0.
402              DO j = 1, knon
403                 i = ni(j)
404                 d_ts(i, nsrf) = y_d_ts(j)
405                 falbe(i, nsrf) = yalb(j)
406                 fsnow(i, nsrf) = snow(j)
407                 qsurf(i, nsrf) = yqsurf(j)
408                 frugs(i, nsrf) = yz0_new(j)
409                 fluxlat(i, nsrf) = yfluxlat(j)
410                 IF (nsrf == is_oce) THEN
411                    rugmer(i) = yrugm(j)
412                    frugs(i, nsrf) = yrugm(j)
413                 END IF
414                 agesno(i, nsrf) = yagesno(j)
415                 fqcalving(i, nsrf) = y_fqcalving(j)
416                 ffonte(i, nsrf) = y_ffonte(j)
417                 cdragh(i) = cdragh(i) + ycdragh(j) * ypct(j)
418                 cdragm(i) = cdragm(i) + ycdragm(j) * ypct(j)
419                 dflux_t(i) = dflux_t(i) + y_dflux_t(j)
420                 dflux_q(i) = dflux_q(i) + y_dflux_q(j)
421              END DO
422              IF (nsrf == is_ter) THEN
423                 qsol(ni(:knon)) = yqsol(:knon)
424              else IF (nsrf == is_lic) THEN
425                 DO j = 1, knon
426                    i = ni(j)
427                    run_off_lic_0(i) = y_run_off_lic_0(j)
428                 END DO
429              END IF
430    
431              ftsoil(:, :, nsrf) = 0.
432              ftsoil(ni(:knon), :, nsrf) = ytsoil(:knon, :)
433    
434              DO j = 1, knon
435                 i = ni(j)
436                 DO k = 1, klev
437                    d_t(i, k) = d_t(i, k) + y_d_t(j, k)
438                    d_q(i, k) = d_q(i, k) + y_d_q(j, k)
439                    d_u(i, k) = d_u(i, k) + y_d_u(j, k)
440                    d_v(i, k) = d_v(i, k) + y_d_v(j, k)
441                 END DO
442              END DO
443    
444              forall (k = 2:klev) coefh(ni(:knon), k) &
445                   = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypct(:knon)
446    
447              ! diagnostic t, q a 2m et u, v a 10m
448    
449              DO j = 1, knon
450                 i = ni(j)
451                 u1(j) = yu(j, 1) + y_d_u(j, 1)
452                 v1(j) = yv(j, 1) + y_d_v(j, 1)
453                 tair1(j) = yt(j, 1) + y_d_t(j, 1)
454                 qair1(j) = yq(j, 1) + y_d_q(j, 1)
455                 zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
456                      1))) * (ypaprs(j, 1)-ypplay(j, 1))
457                 tairsol(j) = yts(j) + y_d_ts(j)
458                 rugo1(j) = yrugos(j)
459                 IF (nsrf == is_oce) THEN
460                    rugo1(j) = frugs(i, nsrf)
461                 END IF
462                 psfce(j) = ypaprs(j, 1)
463                 patm(j) = ypplay(j, 1)
464    
465                 qairsol(j) = yqsurf(j)
466              END DO
467    
468              CALL stdlevvar(klon, knon, nsrf, u1(:knon), v1(:knon), tair1(:knon), &
469                   qair1, zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, &
470                   yq2m, yt10m, yq10m, wind10m(:knon), ustar(:knon))
471    
472              DO j = 1, knon
473                 i = ni(j)
474                 t2m(i, nsrf) = yt2m(j)
475                 q2m(i, nsrf) = yq2m(j)
476    
477                 u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
478                      / sqrt(u1(j)**2 + v1(j)**2)
479                 v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
480                      / sqrt(u1(j)**2 + v1(j)**2)
481              END DO
482    
483              CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
484                   y_flux_q(:knon), yu, yv, yt, yq, ypblh(:knon), ycapcl, &
485                   yoliqcl, ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)
486    
487              DO j = 1, knon
488                 i = ni(j)
489                 pblh(i, nsrf) = ypblh(j)
490                 plcl(i, nsrf) = ylcl(j)
491                 capcl(i, nsrf) = ycapcl(j)
492                 oliqcl(i, nsrf) = yoliqcl(j)
493                 cteicl(i, nsrf) = ycteicl(j)
494                 pblt(i, nsrf) = ypblt(j)
495                 therm(i, nsrf) = ytherm(j)
496                 trmb1(i, nsrf) = ytrmb1(j)
497                 trmb2(i, nsrf) = ytrmb2(j)
498                 trmb3(i, nsrf) = ytrmb3(j)
499              END DO
500    
501              DO j = 1, knon
502                 DO k = 1, klev + 1
503                    i = ni(j)
504                    q2(i, k, nsrf) = yq2(j, k)
505                 END DO
506              END DO
507           else
508              fsnow(:, nsrf) = 0.
509           end IF if_knon
510        END DO loop_surface
511    
512        ! On utilise les nouvelles surfaces
513        frugs(:, is_oce) = rugmer
514        pctsrf(:, is_oce) = pctsrf_new_oce
515        pctsrf(:, is_sic) = pctsrf_new_sic
516    
517    rugos(:, is_oce) = rugmer      firstcal = .false.
   pctsrf = pctsrf_new  
518    
519  END SUBROUTINE clmain    END SUBROUTINE clmain
520    
521    end module clmain_m

Legend:
Removed from v.15  
changed lines
  Added in v.250

  ViewVC Help
Powered by ViewVC 1.1.21