/[lmdze]/trunk/Sources/phylmd/clmain.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/clmain.f

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

revision 15 by guez, Fri Aug 1 15:24:12 2008 UTC revision 38 by guez, Thu Jan 6 17:52:19 2011 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, itap, date0, pctsrf, pctsrf_new, t, q, u, v,&
8    ! A rajouter: conservation de l'albedo         jour, rmu0, co2_ppm, ok_veget, ocean, npas, nexca, ts,&
9           soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil,&
10           qsol, paprs, pplay, snow, qsurf, evap, albe, alblw, fluxlat,&
11           rain_f, snow_f, solsw, sollw, sollwdown, fder, rlon, rlat, cufi,&
12           cvfi, rugos, debut, lafin, agesno, rugoro, d_t, d_q, d_u, d_v,&
13           d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2,&
14           dflux_t, dflux_q, zcoefh, zu1, zv1, t2m, q2m, u10m, v10m, pblh,&
15           capcl, oliqcl, cteicl, pblt, therm, trmb1, trmb2, trmb3, plcl,&
16           fqcalving, ffonte, run_off_lic_0, flux_o, flux_g, tslab, seaice)
17    
18        ! From phylmd/clmain.F, version 1.6 2005/11/16 14:47:19
19    
20        ! Tout ce qui a trait aux traceurs est dans phytrac maintenant.
21        ! Pour l'instant le calcul de la couche limite pour les traceurs
22        ! se fait avec cltrac et ne tient pas compte de la différentiation
23        ! des sous-fractions de sol.
24    
25        ! Pour pouvoir extraire les coefficients d'échanges et le vent
26        ! dans la première couche, trois champs supplémentaires ont été créés :
27        ! zcoefh, zu1 et zv1. Pour l'instant nous avons moyenné les valeurs
28        ! de ces trois champs sur les 4 sous-surfaces du modèle. Dans l'avenir
29        ! si les informations des sous-surfaces doivent être prises en compte
30        ! il faudra sortir ces mêmes champs en leur ajoutant une dimension,
31        ! c'est a dire nbsrf (nombre de sous-surfaces).
32    
33        ! Auteur Z.X. Li (LMD/CNRS) date: 1993/08/18
34        ! Objet : interface de "couche limite" (diffusion verticale)
35    
36        ! Arguments:
37        ! dtime----input-R- interval du temps (secondes)
38        ! itap-----input-I- numero du pas de temps
39        ! date0----input-R- jour initial
40        ! t--------input-R- temperature (K)
41        ! q--------input-R- vapeur d'eau (kg/kg)
42        ! u--------input-R- vitesse u
43        ! v--------input-R- vitesse v
44        ! ts-------input-R- temperature du sol (en Kelvin)
45        ! paprs----input-R- pression a intercouche (Pa)
46        ! pplay----input-R- pression au milieu de couche (Pa)
47        ! radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2
48        ! rlat-----input-R- latitude en degree
49        ! rugos----input-R- longeur de rugosite (en m)
50        ! cufi-----input-R- resolution des mailles en x (m)
51        ! cvfi-----input-R- resolution des mailles en y (m)
52    
53        ! d_t------output-R- le changement pour "t"
54        ! d_q------output-R- le changement pour "q"
55        ! d_u------output-R- le changement pour "u"
56        ! d_v------output-R- le changement pour "v"
57        ! d_ts-----output-R- le changement pour "ts"
58        ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
59        !                    (orientation positive vers le bas)
60        ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
61        ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
62        ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
63        ! dflux_t derive du flux sensible
64        ! dflux_q derive du flux latent
65        !IM "slab" ocean
66        ! flux_g---output-R-  flux glace (pour OCEAN='slab  ')
67        ! flux_o---output-R-  flux ocean (pour OCEAN='slab  ')
68    
69        ! tslab-in/output-R temperature du slab ocean (en Kelvin)
70        ! uniqmnt pour slab
71    
72        ! seaice---output-R-  glace de mer (kg/m2) (pour OCEAN='slab  ')
73        !cc
74        ! ffonte----Flux thermique utilise pour fondre la neige
75        ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
76        !           hauteur de neige, en kg/m2/s
77        ! on rajoute en output yu1 et yv1 qui sont les vents dans
78        ! la premiere couche
79        ! ces 4 variables sont maintenant traites dans phytrac
80        ! itr--------input-I- nombre de traceurs
81        ! tr---------input-R- q. de traceurs
82        ! flux_surf--input-R- flux de traceurs a la surface
83        ! d_tr-------output-R tendance de traceurs
84        !IM cf. AM : PBL
85        ! trmb1-------deep_cape
86        ! trmb2--------inhibition
87        ! trmb3-------Point Omega
88        ! Cape(klon)-------Cape du thermique
89        ! EauLiq(klon)-------Eau liqu integr du thermique
90        ! ctei(klon)-------Critere d'instab d'entrainmt des nuages de CL
91        ! lcl------- Niveau de condensation
92        ! pblh------- HCL
93        ! pblT------- T au nveau HCL
94    
95        USE histcom, ONLY : histbeg_totreg, histdef, histend, histsync
96        use histwrite_m, only: histwrite
97        use calendar, ONLY : ymds2ju
98        USE dimens_m, ONLY : iim, jjm
99        USE indicesol, ONLY : epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
100        USE dimphy, ONLY : klev, klon, zmasq
101        USE dimsoil, ONLY : nsoilmx
102        USE temps, ONLY : annee_ref, itau_phy
103        USE dynetat0_m, ONLY : day_ini
104        USE iniprint, ONLY : prt_level
105        USE suphec_m, ONLY : rd, rg, rkappa
106        USE conf_phys_m, ONLY : iflag_pbl
107        USE gath_cpl, ONLY : gath2cpl
108        use hbtm_m, only: hbtm
109    
110        REAL, INTENT (IN) :: dtime
111        REAL date0
112        INTEGER, INTENT (IN) :: itap
113        REAL t(klon, klev), q(klon, klev)
114        REAL u(klon, klev), v(klon, klev)
115        REAL, INTENT (IN) :: paprs(klon, klev+1)
116        REAL, INTENT (IN) :: pplay(klon, klev)
117        REAL, INTENT (IN) :: rlon(klon), rlat(klon)
118        REAL cufi(klon), cvfi(klon)
119        REAL d_t(klon, klev), d_q(klon, klev)
120        REAL d_u(klon, klev), d_v(klon, klev)
121        REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)
122        REAL dflux_t(klon), dflux_q(klon)
123        !IM "slab" ocean
124        REAL flux_o(klon), flux_g(klon)
125        REAL y_flux_o(klon), y_flux_g(klon)
126        REAL tslab(klon), ytslab(klon)
127        REAL seaice(klon), y_seaice(klon)
128        REAL y_fqcalving(klon), y_ffonte(klon)
129        REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)
130        REAL run_off_lic_0(klon), y_run_off_lic_0(klon)
131    
132        REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)
133        REAL rugmer(klon), agesno(klon, nbsrf)
134        REAL, INTENT (IN) :: rugoro(klon)
135        REAL cdragh(klon), cdragm(klon)
136        ! jour de l'annee en cours                
137        INTEGER jour
138        REAL rmu0(klon) ! cosinus de l'angle solaire zenithal    
139        ! taux CO2 atmosphere                    
140        REAL co2_ppm
141        LOGICAL, INTENT (IN) :: debut
142        LOGICAL, INTENT (IN) :: lafin
143        LOGICAL ok_veget
144        CHARACTER (len=*), INTENT (IN) :: ocean
145        INTEGER npas, nexca
146    
147        REAL pctsrf(klon, nbsrf)
148        REAL ts(klon, nbsrf)
149        REAL d_ts(klon, nbsrf)
150        REAL snow(klon, nbsrf)
151        REAL qsurf(klon, nbsrf)
152        REAL evap(klon, nbsrf)
153        REAL albe(klon, nbsrf)
154        REAL alblw(klon, nbsrf)
155    
156        REAL fluxlat(klon, nbsrf)
157    
158        REAL rain_f(klon), snow_f(klon)
159        REAL fder(klon)
160    
161        REAL sollw(klon, nbsrf), solsw(klon, nbsrf), sollwdown(klon)
162        REAL rugos(klon, nbsrf)
163        ! la nouvelle repartition des surfaces sortie de l'interface
164        REAL pctsrf_new(klon, nbsrf)
165    
166        REAL zcoefh(klon, klev)
167        REAL zu1(klon)
168        REAL zv1(klon)
169    
170        !$$$ PB ajout pour soil
171        LOGICAL, INTENT (IN) :: soil_model
172        !IM ajout seuils cdrm, cdrh
173        REAL cdmmax, cdhmax
174    
175        REAL ksta, ksta_ter
176        LOGICAL ok_kzmin
177    
178        REAL ftsoil(klon, nsoilmx, nbsrf)
179        REAL ytsoil(klon, nsoilmx)
180        REAL qsol(klon)
181    
182        EXTERNAL clqh, clvent, coefkz, calbeta, cltrac
183    
184        REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
185        REAL yalb(klon)
186        REAL yalblw(klon)
187        REAL yu1(klon), yv1(klon)
188        REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)
189        REAL yrain_f(klon), ysnow_f(klon)
190        REAL ysollw(klon), ysolsw(klon), ysollwdown(klon)
191        REAL yfder(klon), ytaux(klon), ytauy(klon)
192        REAL yrugm(klon), yrads(klon), yrugoro(klon)
193    
194        REAL yfluxlat(klon)
195    
196        REAL y_d_ts(klon)
197        REAL y_d_t(klon, klev), y_d_q(klon, klev)
198        REAL y_d_u(klon, klev), y_d_v(klon, klev)
199        REAL y_flux_t(klon, klev), y_flux_q(klon, klev)
200        REAL y_flux_u(klon, klev), y_flux_v(klon, klev)
201        REAL y_dflux_t(klon), y_dflux_q(klon)
202        REAL ycoefh(klon, klev), ycoefm(klon, klev)
203        REAL yu(klon, klev), yv(klon, klev)
204        REAL yt(klon, klev), yq(klon, klev)
205        REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)
206    
207        LOGICAL ok_nonloc
208        PARAMETER (ok_nonloc=.FALSE.)
209        REAL ycoefm0(klon, klev), ycoefh0(klon, klev)
210    
211        !IM 081204 hcl_Anne ? BEG
212        REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)
213        REAL ykmm(klon, klev+1), ykmn(klon, klev+1)
214        REAL ykmq(klon, klev+1)
215        REAL yq2(klon, klev+1), q2(klon, klev+1, nbsrf)
216        REAL q2diag(klon, klev+1)
217        !IM 081204 hcl_Anne ? END
218    
219        REAL u1lay(klon), v1lay(klon)
220        REAL delp(klon, klev)
221        INTEGER i, k, nsrf
222    
223        INTEGER ni(klon), knon, j
224        ! Introduction d'une variable "pourcentage potentiel" pour tenir compte
225        ! des eventuelles apparitions et/ou disparitions de la glace de mer
226        REAL pctsrf_pot(klon, nbsrf)
227    
228        REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
229    
230        ! maf pour sorties IOISPL en cas de debugagage
231    
232        CHARACTER (80) cldebug
233        SAVE cldebug
234        CHARACTER (8) cl_surf(nbsrf)
235        SAVE cl_surf
236        INTEGER nhoridbg, nidbg
237        SAVE nhoridbg, nidbg
238        INTEGER ndexbg(iim*(jjm+1))
239        REAL zx_lon(iim, jjm+1), zx_lat(iim, jjm+1), zjulian
240        REAL tabindx(klon)
241        REAL debugtab(iim, jjm+1)
242        LOGICAL first_appel
243        SAVE first_appel
244        DATA first_appel/ .TRUE./
245        LOGICAL :: debugindex = .FALSE.
246        INTEGER idayref
247        REAL t2m(klon, nbsrf), q2m(klon, nbsrf)
248        REAL u10m(klon, nbsrf), v10m(klon, nbsrf)
249    
250        REAL yt2m(klon), yq2m(klon), yu10m(klon)
251        REAL yustar(klon)
252        ! -- LOOP
253        REAL yu10mx(klon)
254        REAL yu10my(klon)
255        REAL ywindsp(klon)
256        ! -- LOOP
257    
258        REAL yt10m(klon), yq10m(klon)
259        !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds
260        ! physiq ce qui permet de sortir les grdeurs par sous surface)
261        REAL pblh(klon, nbsrf)
262        REAL plcl(klon, nbsrf)
263        REAL capcl(klon, nbsrf)
264        REAL oliqcl(klon, nbsrf)
265        REAL cteicl(klon, nbsrf)
266        REAL pblt(klon, nbsrf)
267        REAL therm(klon, nbsrf)
268        REAL trmb1(klon, nbsrf)
269        REAL trmb2(klon, nbsrf)
270        REAL trmb3(klon, nbsrf)
271        REAL ypblh(klon)
272        REAL ylcl(klon)
273        REAL ycapcl(klon)
274        REAL yoliqcl(klon)
275        REAL ycteicl(klon)
276        REAL ypblt(klon)
277        REAL ytherm(klon)
278        REAL ytrmb1(klon)
279        REAL ytrmb2(klon)
280        REAL ytrmb3(klon)
281        REAL y_cd_h(klon), y_cd_m(klon)
282        REAL uzon(klon), vmer(klon)
283        REAL tair1(klon), qair1(klon), tairsol(klon)
284        REAL psfce(klon), patm(klon)
285    
286        REAL qairsol(klon), zgeo1(klon)
287        REAL rugo1(klon)
288    
289        ! utiliser un jeu de fonctions simples              
290        LOGICAL zxli
291        PARAMETER (zxli=.FALSE.)
292    
293        REAL zt, zqs, zdelta, zcor
294        REAL t_coup
295        PARAMETER (t_coup=273.15)
296    
297        CHARACTER (len=20) :: modname = 'clmain'
298    
299        !------------------------------------------------------------
300    
301        ! initialisation Anne
302        ytherm = 0.
303    
304        IF (debugindex .AND. first_appel) THEN
305           first_appel = .FALSE.
306    
307           ! initialisation sorties netcdf
308    
309           idayref = day_ini
310           CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
311           CALL gr_fi_ecrit(1, klon, iim, jjm+1, rlon, zx_lon)
312           DO i = 1, iim
313              zx_lon(i, 1) = rlon(i+1)
314              zx_lon(i, jjm+1) = rlon(i+1)
315           END DO
316           CALL gr_fi_ecrit(1, klon, iim, jjm+1, rlat, zx_lat)
317           cldebug = 'sous_index'
318           CALL histbeg_totreg(cldebug, zx_lon(:, 1), zx_lat(1, :), 1, &
319                iim, 1, jjm+1, itau_phy, zjulian, dtime, nhoridbg, nidbg)
320           ! no vertical axis
321           cl_surf(1) = 'ter'
322           cl_surf(2) = 'lic'
323           cl_surf(3) = 'oce'
324           cl_surf(4) = 'sic'
325           DO nsrf = 1, nbsrf
326              CALL histdef(nidbg, cl_surf(nsrf), cl_surf(nsrf), '-', iim, jjm+1, &
327                   nhoridbg, 1, 1, 1, -99, 'inst', dtime, dtime)
328           END DO
329           CALL histend(nidbg)
330           CALL histsync(nidbg)
331        END IF
332    
333        DO k = 1, klev ! epaisseur de couche
334           DO i = 1, klon
335              delp(i, k) = paprs(i, k) - paprs(i, k+1)
336           END DO
337        END DO
338        DO i = 1, klon ! vent de la premiere couche
339           zx_alf1 = 1.0
340           zx_alf2 = 1.0 - zx_alf1
341           u1lay(i) = u(i, 1)*zx_alf1 + u(i, 2)*zx_alf2
342           v1lay(i) = v(i, 1)*zx_alf1 + v(i, 2)*zx_alf2
343        END DO
344    
345        ! initialisation:
346    
347        DO i = 1, klon
348           rugmer(i) = 0.0
349           cdragh(i) = 0.0
350           cdragm(i) = 0.0
351           dflux_t(i) = 0.0
352           dflux_q(i) = 0.0
353           zu1(i) = 0.0
354           zv1(i) = 0.0
355        END DO
356        ypct = 0.0
357        yts = 0.0
358        ysnow = 0.0
359        yqsurf = 0.0
360        yalb = 0.0
361        yalblw = 0.0
362        yrain_f = 0.0
363        ysnow_f = 0.0
364        yfder = 0.0
365        ytaux = 0.0
366        ytauy = 0.0
367        ysolsw = 0.0
368        ysollw = 0.0
369        ysollwdown = 0.0
370        yrugos = 0.0
371        yu1 = 0.0
372        yv1 = 0.0
373        yrads = 0.0
374        ypaprs = 0.0
375        ypplay = 0.0
376        ydelp = 0.0
377        yu = 0.0
378        yv = 0.0
379        yt = 0.0
380        yq = 0.0
381        pctsrf_new = 0.0
382        y_flux_u = 0.0
383        y_flux_v = 0.0
384        !$$ PB
385        y_dflux_t = 0.0
386        y_dflux_q = 0.0
387        ytsoil = 999999.
388        yrugoro = 0.
389        ! -- LOOP
390        yu10mx = 0.0
391        yu10my = 0.0
392        ywindsp = 0.0
393        ! -- LOOP
394        DO nsrf = 1, nbsrf
395           DO i = 1, klon
396              d_ts(i, nsrf) = 0.0
397           END DO
398        END DO
399        !§§§ PB
400        yfluxlat = 0.
401        flux_t = 0.
402        flux_q = 0.
403        flux_u = 0.
404        flux_v = 0.
405        DO k = 1, klev
406           DO i = 1, klon
407              d_t(i, k) = 0.0
408              d_q(i, k) = 0.0
409              d_u(i, k) = 0.0
410              d_v(i, k) = 0.0
411              zcoefh(i, k) = 0.0
412           END DO
413        END DO
414    
415        ! Boucler sur toutes les sous-fractions du sol:
416    
417        ! Initialisation des "pourcentages potentiels". On considère ici qu'on
418        ! peut avoir potentiellement de la glace sur tout le domaine océanique
419        ! (à affiner)
420    
421        pctsrf_pot = pctsrf
422        pctsrf_pot(:, is_oce) = 1. - zmasq
423        pctsrf_pot(:, is_sic) = 1. - zmasq
424    
425        DO nsrf = 1, nbsrf
426           ! chercher les indices:
427           ni = 0
428           knon = 0
429           DO i = 1, klon
430              ! pour determiner le domaine a traiter on utilise les surfaces
431              ! "potentielles"
432              IF (pctsrf_pot(i, nsrf) > epsfra) THEN
433                 knon = knon + 1
434                 ni(knon) = i
435              END IF
436           END DO
437    
438           ! variables pour avoir une sortie IOIPSL des INDEX
439           IF (debugindex) THEN
440              tabindx = 0.
441              DO i = 1, knon
442                 tabindx(i) = real(i)
443              END DO
444              debugtab = 0.
445              ndexbg = 0
446              CALL gath2cpl(tabindx, debugtab, klon, knon, iim, jjm, ni)
447              CALL histwrite(nidbg, cl_surf(nsrf), itap, debugtab)
448           END IF
449    
450           IF (knon==0) CYCLE
451    
452           DO j = 1, knon
453              i = ni(j)
454              ypct(j) = pctsrf(i, nsrf)
455              yts(j) = ts(i, nsrf)
456              ytslab(i) = tslab(i)
457              ysnow(j) = snow(i, nsrf)
458              yqsurf(j) = qsurf(i, nsrf)
459              yalb(j) = albe(i, nsrf)
460              yalblw(j) = alblw(i, nsrf)
461              yrain_f(j) = rain_f(i)
462              ysnow_f(j) = snow_f(i)
463              yagesno(j) = agesno(i, nsrf)
464              yfder(j) = fder(i)
465              ytaux(j) = flux_u(i, 1, nsrf)
466              ytauy(j) = flux_v(i, 1, nsrf)
467              ysolsw(j) = solsw(i, nsrf)
468              ysollw(j) = sollw(i, nsrf)
469              ysollwdown(j) = sollwdown(i)
470              yrugos(j) = rugos(i, nsrf)
471              yrugoro(j) = rugoro(i)
472              yu1(j) = u1lay(i)
473              yv1(j) = v1lay(i)
474              yrads(j) = ysolsw(j) + ysollw(j)
475              ypaprs(j, klev+1) = paprs(i, klev+1)
476              y_run_off_lic_0(j) = run_off_lic_0(i)
477              yu10mx(j) = u10m(i, nsrf)
478              yu10my(j) = v10m(i, nsrf)
479              ywindsp(j) = sqrt(yu10mx(j)*yu10mx(j)+yu10my(j)*yu10my(j))
480           END DO
481    
482           !     IF bucket model for continent, copy soil water content
483           IF (nsrf==is_ter .AND. .NOT. ok_veget) THEN
484              DO j = 1, knon
485                 i = ni(j)
486                 yqsol(j) = qsol(i)
487              END DO
488           ELSE
489              yqsol = 0.
490           END IF
491           !$$$ PB ajour pour soil
492           DO k = 1, nsoilmx
493              DO j = 1, knon
494                 i = ni(j)
495                 ytsoil(j, k) = ftsoil(i, k, nsrf)
496              END DO
497           END DO
498           DO k = 1, klev
499              DO j = 1, knon
500                 i = ni(j)
501                 ypaprs(j, k) = paprs(i, k)
502                 ypplay(j, k) = pplay(i, k)
503                 ydelp(j, k) = delp(i, k)
504                 yu(j, k) = u(i, k)
505                 yv(j, k) = v(i, k)
506                 yt(j, k) = t(i, k)
507                 yq(j, k) = q(i, k)
508              END DO
509           END DO
510    
511           ! calculer Cdrag et les coefficients d'echange
512           CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts,&
513                yrugos, yu, yv, yt, yq, yqsurf, ycoefm, ycoefh)
514           !IM 081204 BEG
515           !CR test
516           IF (iflag_pbl==1) THEN
517              !IM 081204 END
518              CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
519              DO k = 1, klev
520                 DO i = 1, knon
521                    ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))
522                    ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))
523                 END DO
524              END DO
525           END IF
526    
527           !IM cf JLD : on seuille ycoefm et ycoefh
528           IF (nsrf==is_oce) THEN
529              DO j = 1, knon
530                 !           ycoefm(j, 1)=min(ycoefm(j, 1), 1.1E-3)
531                 ycoefm(j, 1) = min(ycoefm(j, 1), cdmmax)
532                 !           ycoefh(j, 1)=min(ycoefh(j, 1), 1.1E-3)
533                 ycoefh(j, 1) = min(ycoefh(j, 1), cdhmax)
534              END DO
535           END IF
536    
537           !IM: 261103
538           IF (ok_kzmin) THEN
539              !IM cf FH: 201103 BEG
540              !   Calcul d'une diffusion minimale pour les conditions tres stables.
541              CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, ycoefm, &
542                   ycoefm0, ycoefh0)
543    
544              IF (1==1) THEN
545                 DO k = 1, klev
546                    DO i = 1, knon
547                       ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))
548                       ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))
549                    END DO
550                 END DO
551              END IF
552              !IM cf FH: 201103 END
553              !IM: 261103
554           END IF !ok_kzmin
555    
556           IF (iflag_pbl>=3) THEN
557              ! MELLOR ET YAMADA adapté à Mars, Richard Fournier et Frédéric Hourdin
558              yzlay(1:knon, 1) = rd*yt(1:knon, 1)/(0.5*(ypaprs(1:knon, &
559                   1)+ypplay(1:knon, 1)))*(ypaprs(1:knon, 1)-ypplay(1:knon, 1))/rg
560              DO k = 2, klev
561                 yzlay(1:knon, k) = yzlay(1:knon, k-1) &
562                      + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
563                      / ypaprs(1:knon, k) &
564                      * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
565              END DO
566              DO k = 1, klev
567                 yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
568                      / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))
569              END DO
570              yzlev(1:knon, 1) = 0.
571              yzlev(1:knon, klev+1) = 2.*yzlay(1:knon, klev) - yzlay(1:knon, klev-1)
572              DO k = 2, klev
573                 yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
574              END DO
575              DO k = 1, klev + 1
576                 DO j = 1, knon
577                    i = ni(j)
578                    yq2(j, k) = q2(i, k, nsrf)
579                 END DO
580              END DO
581    
582              !   Bug introduit volontairement pour converger avec les resultats
583              !  du papier sur les thermiques.
584              IF (1==1) THEN
585                 y_cd_m(1:knon) = ycoefm(1:knon, 1)
586                 y_cd_h(1:knon) = ycoefh(1:knon, 1)
587              ELSE
588                 y_cd_h(1:knon) = ycoefm(1:knon, 1)
589                 y_cd_m(1:knon) = ycoefh(1:knon, 1)
590              END IF
591              CALL ustarhb(knon, yu, yv, y_cd_m, yustar)
592    
593              IF (prt_level>9) THEN
594                 PRINT *, 'USTAR = ', yustar
595              END IF
596    
597              !   iflag_pbl peut etre utilise comme longuer de melange
598    
599              IF (iflag_pbl>=11) THEN
600                 CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &
601                      yu, yv, yteta, y_cd_m, yq2, q2diag, ykmm, ykmn, yustar, &
602                      iflag_pbl)
603              ELSE
604                 CALL yamada4(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, yu, &
605                      yv, yteta, y_cd_m, yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
606              END IF
607    
608              ycoefm(1:knon, 1) = y_cd_m(1:knon)
609              ycoefh(1:knon, 1) = y_cd_h(1:knon)
610              ycoefm(1:knon, 2:klev) = ykmm(1:knon, 2:klev)
611              ycoefh(1:knon, 2:klev) = ykmn(1:knon, 2:klev)
612           END IF
613    
614           ! calculer la diffusion des vitesses "u" et "v"
615           CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yu, ypaprs, ypplay, &
616                ydelp, y_d_u, y_flux_u)
617           CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yv, ypaprs, ypplay, &
618                ydelp, y_d_v, y_flux_v)
619    
620           ! pour le couplage
621           ytaux = y_flux_u(:, 1)
622           ytauy = y_flux_v(:, 1)
623    
624           ! FH modif sur le cdrag temperature
625           !$$$PB : déplace dans clcdrag
626           !$$$      do i=1, knon
627           !$$$         ycoefh(i, 1)=ycoefm(i, 1)*0.8
628           !$$$      enddo
629    
630           ! calculer la diffusion de "q" et de "h"
631           CALL clqh(dtime, itap, date0, jour, debut, lafin, rlon, rlat,&
632                cufi, cvfi, knon, nsrf, ni, pctsrf, soil_model, ytsoil,&
633                yqsol, ok_veget, ocean, npas, nexca, rmu0, co2_ppm, yrugos,&
634                yrugoro, yu1, yv1, ycoefh, yt, yq, yts, ypaprs, ypplay,&
635                ydelp, yrads, yalb, yalblw, ysnow, yqsurf, yrain_f, ysnow_f, &
636                yfder, ytaux, ytauy, ywindsp, ysollw, ysollwdown, ysolsw,&
637                yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, y_d_ts,&
638                yz0_new, y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,&
639                y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, y_flux_g,&
640                ytslab, y_seaice)
641    
642           ! calculer la longueur de rugosite sur ocean
643           yrugm = 0.
644           IF (nsrf==is_oce) THEN
645              DO j = 1, knon
646                 yrugm(j) = 0.018*ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &
647                      0.11*14E-6/sqrt(ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2))
648                 yrugm(j) = max(1.5E-05, yrugm(j))
649              END DO
650           END IF
651           DO j = 1, knon
652              y_dflux_t(j) = y_dflux_t(j)*ypct(j)
653              y_dflux_q(j) = y_dflux_q(j)*ypct(j)
654              yu1(j) = yu1(j)*ypct(j)
655              yv1(j) = yv1(j)*ypct(j)
656           END DO
657    
658           DO k = 1, klev
659              DO j = 1, knon
660                 i = ni(j)
661                 ycoefh(j, k) = ycoefh(j, k)*ypct(j)
662                 ycoefm(j, k) = ycoefm(j, k)*ypct(j)
663                 y_d_t(j, k) = y_d_t(j, k)*ypct(j)
664                 y_d_q(j, k) = y_d_q(j, k)*ypct(j)
665                 !§§§ PB
666                 flux_t(i, k, nsrf) = y_flux_t(j, k)
667                 flux_q(i, k, nsrf) = y_flux_q(j, k)
668                 flux_u(i, k, nsrf) = y_flux_u(j, k)
669                 flux_v(i, k, nsrf) = y_flux_v(j, k)
670                 !$$$ PB        y_flux_t(j, k) = y_flux_t(j, k) * ypct(j)
671                 !$$$ PB        y_flux_q(j, k) = y_flux_q(j, k) * ypct(j)
672                 y_d_u(j, k) = y_d_u(j, k)*ypct(j)
673                 y_d_v(j, k) = y_d_v(j, k)*ypct(j)
674                 !$$$ PB        y_flux_u(j, k) = y_flux_u(j, k) * ypct(j)
675                 !$$$ PB        y_flux_v(j, k) = y_flux_v(j, k) * ypct(j)
676              END DO
677           END DO
678    
679           evap(:, nsrf) = -flux_q(:, 1, nsrf)
680    
681           albe(:, nsrf) = 0.
682           alblw(:, nsrf) = 0.
683           snow(:, nsrf) = 0.
684           qsurf(:, nsrf) = 0.
685           rugos(:, nsrf) = 0.
686           fluxlat(:, nsrf) = 0.
687           DO j = 1, knon
688              i = ni(j)
689              d_ts(i, nsrf) = y_d_ts(j)
690              albe(i, nsrf) = yalb(j)
691              alblw(i, nsrf) = yalblw(j)
692              snow(i, nsrf) = ysnow(j)
693              qsurf(i, nsrf) = yqsurf(j)
694              rugos(i, nsrf) = yz0_new(j)
695              fluxlat(i, nsrf) = yfluxlat(j)
696              !$$$ pb         rugmer(i) = yrugm(j)
697              IF (nsrf==is_oce) THEN
698                 rugmer(i) = yrugm(j)
699                 rugos(i, nsrf) = yrugm(j)
700              END IF
701              !IM cf JLD ??
702              agesno(i, nsrf) = yagesno(j)
703              fqcalving(i, nsrf) = y_fqcalving(j)
704              ffonte(i, nsrf) = y_ffonte(j)
705              cdragh(i) = cdragh(i) + ycoefh(j, 1)
706              cdragm(i) = cdragm(i) + ycoefm(j, 1)
707              dflux_t(i) = dflux_t(i) + y_dflux_t(j)
708              dflux_q(i) = dflux_q(i) + y_dflux_q(j)
709              zu1(i) = zu1(i) + yu1(j)
710              zv1(i) = zv1(i) + yv1(j)
711           END DO
712           IF (nsrf==is_ter) THEN
713              DO j = 1, knon
714                 i = ni(j)
715                 qsol(i) = yqsol(j)
716              END DO
717           END IF
718           IF (nsrf==is_lic) THEN
719              DO j = 1, knon
720                 i = ni(j)
721                 run_off_lic_0(i) = y_run_off_lic_0(j)
722              END DO
723           END IF
724           !$$$ PB ajout pour soil
725           ftsoil(:, :, nsrf) = 0.
726           DO k = 1, nsoilmx
727              DO j = 1, knon
728                 i = ni(j)
729                 ftsoil(i, k, nsrf) = ytsoil(j, k)
730              END DO
731           END DO
732    
733           DO j = 1, knon
734              i = ni(j)
735              DO k = 1, klev
736                 d_t(i, k) = d_t(i, k) + y_d_t(j, k)
737                 d_q(i, k) = d_q(i, k) + y_d_q(j, k)
738                 !$$$ PB        flux_t(i, k) = flux_t(i, k) + y_flux_t(j, k)
739                 !$$$         flux_q(i, k) = flux_q(i, k) + y_flux_q(j, k)
740                 d_u(i, k) = d_u(i, k) + y_d_u(j, k)
741                 d_v(i, k) = d_v(i, k) + y_d_v(j, k)
742                 !$$$  PB       flux_u(i, k) = flux_u(i, k) + y_flux_u(j, k)
743                 !$$$         flux_v(i, k) = flux_v(i, k) + y_flux_v(j, k)
744                 zcoefh(i, k) = zcoefh(i, k) + ycoefh(j, k)
745              END DO
746           END DO
747    
748           !cc diagnostic t, q a 2m et u, v a 10m
749    
750           DO j = 1, knon
751              i = ni(j)
752              uzon(j) = yu(j, 1) + y_d_u(j, 1)
753              vmer(j) = yv(j, 1) + y_d_v(j, 1)
754              tair1(j) = yt(j, 1) + y_d_t(j, 1)
755              qair1(j) = yq(j, 1) + y_d_q(j, 1)
756              zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &
757                   1)))*(ypaprs(j, 1)-ypplay(j, 1))
758              tairsol(j) = yts(j) + y_d_ts(j)
759              rugo1(j) = yrugos(j)
760              IF (nsrf==is_oce) THEN
761                 rugo1(j) = rugos(i, nsrf)
762              END IF
763              psfce(j) = ypaprs(j, 1)
764              patm(j) = ypplay(j, 1)
765    
766              qairsol(j) = yqsurf(j)
767           END DO
768    
769           CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, zgeo1, &
770                tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, yq10m, &
771                yu10m, yustar)
772           !IM 081204 END
773    
774           DO j = 1, knon
775              i = ni(j)
776              t2m(i, nsrf) = yt2m(j)
777              q2m(i, nsrf) = yq2m(j)
778    
779              ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
780              u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)
781              v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)
782    
783           END DO
784    
785           DO i = 1, knon
786              y_cd_h(i) = ycoefh(i, 1)
787              y_cd_m(i) = ycoefm(i, 1)
788           END DO
789           CALL hbtm(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, &
790                y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &
791                ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)
792    
793           DO j = 1, knon
794              i = ni(j)
795              pblh(i, nsrf) = ypblh(j)
796              plcl(i, nsrf) = ylcl(j)
797              capcl(i, nsrf) = ycapcl(j)
798              oliqcl(i, nsrf) = yoliqcl(j)
799              cteicl(i, nsrf) = ycteicl(j)
800              pblt(i, nsrf) = ypblt(j)
801              therm(i, nsrf) = ytherm(j)
802              trmb1(i, nsrf) = ytrmb1(j)
803              trmb2(i, nsrf) = ytrmb2(j)
804              trmb3(i, nsrf) = ytrmb3(j)
805           END DO
806    
807           DO j = 1, knon
808              DO k = 1, klev + 1
809                 i = ni(j)
810                 q2(i, k, nsrf) = yq2(j, k)
811              END DO
812           END DO
813           !IM "slab" ocean
814           IF (nsrf==is_oce) THEN
815              DO j = 1, knon
816                 ! on projette sur la grille globale
817                 i = ni(j)
818                 IF (pctsrf_new(i, is_oce)>epsfra) THEN
819                    flux_o(i) = y_flux_o(j)
820                 ELSE
821                    flux_o(i) = 0.
822                 END IF
823              END DO
824           END IF
825    
826           IF (nsrf==is_sic) THEN
827              DO j = 1, knon
828                 i = ni(j)
829                 ! On pondère lorsque l'on fait le bilan au sol :
830                 ! flux_g(i) = y_flux_g(j)*ypct(j)
831                 IF (pctsrf_new(i, is_sic)>epsfra) THEN
832                    flux_g(i) = y_flux_g(j)
833                 ELSE
834                    flux_g(i) = 0.
835                 END IF
836              END DO
837    
838           END IF
839           !nsrf.EQ.is_sic                                            
840           IF (ocean=='slab  ') THEN
841              IF (nsrf==is_oce) THEN
842                 tslab(1:klon) = ytslab(1:klon)
843                 seaice(1:klon) = y_seaice(1:klon)
844                 !nsrf                                                      
845              END IF
846              !OCEAN                                                      
847           END IF
848        END DO
849    
850        ! On utilise les nouvelles surfaces
851        ! A rajouter: conservation de l'albedo
852    
853    rugos(:, is_oce) = rugmer      rugos(:, is_oce) = rugmer
854    pctsrf = pctsrf_new      pctsrf = pctsrf_new
855    
856  END SUBROUTINE clmain    END SUBROUTINE clmain
857    
858    end module clmain_m

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

  ViewVC Help
Powered by ViewVC 1.1.21