/[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 251 by guez, Mon Jan 8 14:12:02 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(:knon, :), &
343                   ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
344                   yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
345                   ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))
346    
347              CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
348                   ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
349                   ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
350                   y_flux_u(:knon))
351              CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
352                   ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
353                   ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
354                   y_flux_v(:knon))
355    
356              ! calculer la diffusion de "q" et de "h"
357              CALL clqh(dtime, julien, firstcal, nsrf, ni(:knon), &
358                   ytsoil(:knon, :), yqsol(:knon), mu0, yrugos, yrugoro, &
359                   yu(:knon, 1), yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), &
360                   yt, yq, yts(:knon), ypaprs, ypplay, ydelp, yrads(:knon), &
361                   yalb(:knon), snow(:knon), yqsurf, yrain_f, ysnow_f, &
362                   yfluxlat(:knon), pctsrf_new_sic, yagesno(:knon), y_d_t, y_d_q, &
363                   y_d_ts(:knon), yz0_new, y_flux_t(:knon), y_flux_q(:knon), &
364                   y_dflux_t(:knon), y_dflux_q(:knon), y_fqcalving, y_ffonte, &
365                   y_run_off_lic_0)
366    
367              ! calculer la longueur de rugosite sur ocean
368              yrugm = 0.
369              IF (nsrf == is_oce) THEN
370                 DO j = 1, knon
371                    yrugm(j) = 0.018 * ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2) &
372                         / rg + 0.11 * 14E-6 &
373                         / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2))
374                    yrugm(j) = max(1.5E-05, yrugm(j))
375                 END DO
376              END IF
377              DO j = 1, knon
378                 y_dflux_t(j) = y_dflux_t(j) * ypct(j)
379                 y_dflux_q(j) = y_dflux_q(j) * ypct(j)
380              END DO
381    
382              DO k = 1, klev
383                 DO j = 1, knon
384                    i = ni(j)
385                    y_d_t(j, k) = y_d_t(j, k) * ypct(j)
386                    y_d_q(j, k) = y_d_q(j, k) * ypct(j)
387                    y_d_u(j, k) = y_d_u(j, k) * ypct(j)
388                    y_d_v(j, k) = y_d_v(j, k) * ypct(j)
389                 END DO
390              END DO
391    
392              flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
393              flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
394              flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
395              flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
396    
397              evap(:, nsrf) = -flux_q(:, nsrf)
398    
399              falbe(:, nsrf) = 0.
400              fsnow(:, nsrf) = 0.
401              qsurf(:, nsrf) = 0.
402              frugs(:, nsrf) = 0.
403              DO j = 1, knon
404                 i = ni(j)
405                 d_ts(i, nsrf) = y_d_ts(j)
406                 falbe(i, nsrf) = yalb(j)
407                 fsnow(i, nsrf) = snow(j)
408                 qsurf(i, nsrf) = yqsurf(j)
409                 frugs(i, nsrf) = yz0_new(j)
410                 fluxlat(i, nsrf) = yfluxlat(j)
411                 IF (nsrf == is_oce) THEN
412                    rugmer(i) = yrugm(j)
413                    frugs(i, nsrf) = yrugm(j)
414                 END IF
415                 agesno(i, nsrf) = yagesno(j)
416                 fqcalving(i, nsrf) = y_fqcalving(j)
417                 ffonte(i, nsrf) = y_ffonte(j)
418                 cdragh(i) = cdragh(i) + ycdragh(j) * ypct(j)
419                 cdragm(i) = cdragm(i) + ycdragm(j) * ypct(j)
420                 dflux_t(i) = dflux_t(i) + y_dflux_t(j)
421                 dflux_q(i) = dflux_q(i) + y_dflux_q(j)
422              END DO
423              IF (nsrf == is_ter) THEN
424                 qsol(ni(:knon)) = yqsol(:knon)
425              else IF (nsrf == is_lic) THEN
426                 DO j = 1, knon
427                    i = ni(j)
428                    run_off_lic_0(i) = y_run_off_lic_0(j)
429                 END DO
430              END IF
431    
432              ftsoil(:, :, nsrf) = 0.
433              ftsoil(ni(:knon), :, nsrf) = ytsoil(:knon, :)
434    
435              DO j = 1, knon
436                 i = ni(j)
437                 DO k = 1, klev
438                    d_t(i, k) = d_t(i, k) + y_d_t(j, k)
439                    d_q(i, k) = d_q(i, k) + y_d_q(j, k)
440                    d_u(i, k) = d_u(i, k) + y_d_u(j, k)
441                    d_v(i, k) = d_v(i, k) + y_d_v(j, k)
442                 END DO
443              END DO
444    
445              forall (k = 2:klev) coefh(ni(:knon), k) &
446                   = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypct(:knon)
447    
448              ! diagnostic t, q a 2m et u, v a 10m
449    
450              DO j = 1, knon
451                 i = ni(j)
452                 u1(j) = yu(j, 1) + y_d_u(j, 1)
453                 v1(j) = yv(j, 1) + y_d_v(j, 1)
454                 tair1(j) = yt(j, 1) + y_d_t(j, 1)
455                 qair1(j) = yq(j, 1) + y_d_q(j, 1)
456                 zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
457                      1))) * (ypaprs(j, 1)-ypplay(j, 1))
458                 tairsol(j) = yts(j) + y_d_ts(j)
459                 rugo1(j) = yrugos(j)
460                 IF (nsrf == is_oce) THEN
461                    rugo1(j) = frugs(i, nsrf)
462                 END IF
463                 psfce(j) = ypaprs(j, 1)
464                 patm(j) = ypplay(j, 1)
465    
466                 qairsol(j) = yqsurf(j)
467              END DO
468    
469              CALL stdlevvar(klon, knon, nsrf, u1(:knon), v1(:knon), tair1(:knon), &
470                   qair1, zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, &
471                   yq2m, yt10m, yq10m, wind10m(:knon), ustar(:knon))
472    
473              DO j = 1, knon
474                 i = ni(j)
475                 t2m(i, nsrf) = yt2m(j)
476                 q2m(i, nsrf) = yq2m(j)
477    
478                 u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
479                      / sqrt(u1(j)**2 + v1(j)**2)
480                 v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
481                      / sqrt(u1(j)**2 + v1(j)**2)
482              END DO
483    
484              CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
485                   y_flux_q(:knon), yu, yv, yt, yq, ypblh(:knon), ycapcl, &
486                   yoliqcl, ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)
487    
488              DO j = 1, knon
489                 i = ni(j)
490                 pblh(i, nsrf) = ypblh(j)
491                 plcl(i, nsrf) = ylcl(j)
492                 capcl(i, nsrf) = ycapcl(j)
493                 oliqcl(i, nsrf) = yoliqcl(j)
494                 cteicl(i, nsrf) = ycteicl(j)
495                 pblt(i, nsrf) = ypblt(j)
496                 therm(i, nsrf) = ytherm(j)
497                 trmb1(i, nsrf) = ytrmb1(j)
498                 trmb2(i, nsrf) = ytrmb2(j)
499                 trmb3(i, nsrf) = ytrmb3(j)
500              END DO
501    
502              DO j = 1, knon
503                 DO k = 1, klev + 1
504                    i = ni(j)
505                    q2(i, k, nsrf) = yq2(j, k)
506                 END DO
507              END DO
508           else
509              fsnow(:, nsrf) = 0.
510           end IF if_knon
511        END DO loop_surface
512    
513        ! On utilise les nouvelles surfaces
514        frugs(:, is_oce) = rugmer
515        pctsrf(:, is_oce) = pctsrf_new_oce
516        pctsrf(:, is_sic) = pctsrf_new_sic
517    
518    rugos(:, is_oce) = rugmer      firstcal = .false.
   pctsrf = pctsrf_new  
519    
520  END SUBROUTINE clmain    END SUBROUTINE clmain
521    
522    end module clmain_m

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

  ViewVC Help
Powered by ViewVC 1.1.21