/[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 30 by guez, Thu Apr 1 09:07:28 2010 UTC revision 40 by guez, Tue Feb 22 13:49:36 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, flux_o, flux_g, tslab, seaice)  
   
   ! From phylmd/clmain.F, v 1.6 2005/11/16 14:47:19  
   
   !AA Tout ce qui a trait au traceurs est dans phytrac maintenant  
   !AA pour l'instant le calcul de la couche limite pour les traceurs  
   !AA se fait avec cltrac et ne tient pas compte de la differentiation  
   !AA des sous-fraction de sol.  
   
   !AA Pour pouvoir extraire les coefficient d'echanges et le vent  
   !AA dans la premiere couche, 3 champs supplementaires ont ete crees  
   !AA zcoefh, zu1 et zv1. Pour l'instant nous avons moyenne les valeurs  
   !AA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir  
   !AA si les informations des subsurfaces doivent etre prises en compte  
   !AA il faudra sortir ces memes champs en leur ajoutant une dimension,  
   !AA c'est a dire nbsrf (nbre de subsurface).  
   
   ! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818  
   ! Objet: interface de "couche limite" (diffusion verticale)  
   
   ! Arguments:  
   ! dtime----input-R- interval du temps (secondes)  
   ! itap-----input-I- numero du pas de temps  
   ! date0----input-R- jour initial  
   ! t--------input-R- temperature (K)  
   ! q--------input-R- vapeur d'eau (kg/kg)  
   ! u--------input-R- vitesse u  
   ! v--------input-R- vitesse v  
   ! ts-------input-R- temperature du sol (en Kelvin)  
   ! paprs----input-R- pression a intercouche (Pa)  
   ! pplay----input-R- pression au milieu de couche (Pa)  
   ! radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2  
   ! rlat-----input-R- latitude en degree  
   ! rugos----input-R- longeur de rugosite (en m)  
   ! cufi-----input-R- resolution des mailles en x (m)  
   ! cvfi-----input-R- resolution des mailles en y (m)  
   
   ! d_t------output-R- le changement pour "t"  
   ! d_q------output-R- le changement pour "q"  
   ! d_u------output-R- le changement pour "u"  
   ! d_v------output-R- le changement pour "v"  
   ! d_ts-----output-R- le changement pour "ts"  
   ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)  
   !                    (orientation positive vers le bas)  
   ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)  
   ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal  
   ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal  
   ! dflux_t derive du flux sensible  
   ! dflux_q derive du flux latent  
   !IM "slab" ocean  
   ! flux_g---output-R-  flux glace (pour OCEAN='slab  ')  
   ! flux_o---output-R-  flux ocean (pour OCEAN='slab  ')  
   ! tslab-in/output-R temperature du slab ocean (en Kelvin) ! uniqmnt pour slab  
   ! seaice---output-R-  glace de mer (kg/m2) (pour OCEAN='slab  ')  
   !cc  
   ! ffonte----Flux thermique utilise pour fondre la neige  
   ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la  
   !           hauteur de neige, en kg/m2/s  
   !AA on rajoute en output yu1 et yv1 qui sont les vents dans  
   !AA la premiere couche  
   !AA ces 4 variables sont maintenant traites dans phytrac  
   ! itr--------input-I- nombre de traceurs  
   ! tr---------input-R- q. de traceurs  
   ! flux_surf--input-R- flux de traceurs a la surface  
   ! d_tr-------output-R tendance de traceurs  
   !IM cf. AM : PBL  
   ! trmb1-------deep_cape  
   ! trmb2--------inhibition  
   ! trmb3-------Point Omega  
   ! Cape(klon)-------Cape du thermique  
   ! EauLiq(klon)-------Eau liqu integr du thermique  
   ! ctei(klon)-------Critere d'instab d'entrainmt des nuages de CL  
   ! lcl------- Niveau de condensation  
   ! pblh------- HCL  
   ! pblT------- T au nveau HCL  
   
   !$$$ PB ajout pour soil  
   
   USE histcom, ONLY : histbeg_totreg, histdef, histend, histsync  
   use histwrite_m, only: histwrite  
   use calendar, ONLY : ymds2ju  
   USE dimens_m, ONLY : iim, jjm  
   USE indicesol, ONLY : epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf  
   USE dimphy, ONLY : klev, klon, zmasq  
   USE dimsoil, ONLY : nsoilmx  
   USE temps, ONLY : annee_ref, itau_phy  
   USE dynetat0_m, ONLY : day_ini  
   USE iniprint, ONLY : prt_level  
   USE yomcst, ONLY : rd, rg, rkappa  
   USE conf_phys_m, ONLY : iflag_pbl  
   USE gath_cpl, ONLY : gath2cpl  
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, ksta, ksta_ter, yts,&  
           yrugos, yu, yv, yt, yq, yqsurf, ycoefm, ycoefh)  
      !IM 081204 BEG  
      !CR test  
      IF (iflag_pbl==1) THEN  
         !IM 081204 END  
         CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)  
         DO k = 1, klev  
            DO i = 1, knon  
               ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))  
               ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))  
            END DO  
         END DO  
      END IF  
   
      !IM cf JLD : on seuille ycoefm et ycoefh  
      IF (nsrf==is_oce) THEN  
         DO j = 1, knon  
            !           ycoefm(j, 1)=min(ycoefm(j, 1), 1.1E-3)  
            ycoefm(j, 1) = min(ycoefm(j, 1), cdmmax)  
            !           ycoefh(j, 1)=min(ycoefh(j, 1), 1.1E-3)  
            ycoefh(j, 1) = min(ycoefh(j, 1), cdhmax)  
         END DO  
      END IF  
   
   
      !IM: 261103  
      IF (ok_kzmin) THEN  
         !IM cf FH: 201103 BEG  
         !   Calcul d'une diffusion minimale pour les conditions tres stables.  
         CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, ycoefm, ycoefm0, &  
              ycoefh0)  
         !      call dump2d(iim, jjm-1, ycoefm(2:klon-1, 2), 'KZ         ')  
         !      call dump2d(iim, jjm-1, ycoefm0(2:klon-1, 2), 'KZMIN      ')  
   
         IF (1==1) THEN  
            DO k = 1, klev  
               DO i = 1, knon  
                  ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))  
                  ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))  
               END DO  
            END DO  
         END IF  
         !IM cf FH: 201103 END  
         !IM: 261103  
      END IF !ok_kzmin  
   
      IF (iflag_pbl>=3) THEN  
   
         !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
         ! MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin  
         !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
   
         yzlay(1:knon, 1) = rd*yt(1:knon, 1)/(0.5*(ypaprs(1:knon, &  
              1)+ypplay(1:knon, 1)))*(ypaprs(1:knon, 1)-ypplay(1:knon, 1))/rg  
         DO k = 2, klev  
            yzlay(1:knon, k) = yzlay(1:knon, k-1) &  
                 + rd*0.5*(yt(1:knon, k-1) +yt(1: knon, k)) &  
                 / ypaprs(1:knon, k) *(ypplay(1:knon, k-1)-ypplay(1:knon, k))/ &  
                 rg  
         END DO  
         DO k = 1, klev  
            yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &  
                 / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))  
         END DO  
         yzlev(1:knon, 1) = 0.  
         yzlev(1:knon, klev+1) = 2.*yzlay(1:knon, klev) - yzlay(1:knon, klev-1)  
         DO k = 2, klev  
            yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))  
         END DO  
         DO k = 1, klev + 1  
            DO j = 1, knon  
               i = ni(j)  
               yq2(j, k) = q2(i, k, nsrf)  
            END DO  
         END DO  
   
   
         !   Bug introduit volontairement pour converger avec les resultats  
         !  du papier sur les thermiques.  
         IF (1==1) THEN  
            y_cd_m(1:knon) = ycoefm(1:knon, 1)  
            y_cd_h(1:knon) = ycoefh(1:knon, 1)  
         ELSE  
            y_cd_h(1:knon) = ycoefm(1:knon, 1)  
            y_cd_m(1:knon) = ycoefh(1:knon, 1)  
         END IF  
         CALL ustarhb(knon, yu, yv, y_cd_m, yustar)  
   
         IF (prt_level>9) THEN  
            PRINT *, 'USTAR = ', yustar  
         END IF  
   
         !   iflag_pbl peut etre utilise comme longuer de melange  
   
         IF (iflag_pbl>=11) THEN  
            CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, yu, yv, yteta, &  
                 y_cd_m, yq2, q2diag, ykmm, ykmn, yustar, iflag_pbl)  
         ELSE  
            CALL yamada4(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, yu, yv, yteta, &  
                 y_cd_m, yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)  
         END IF  
   
         ycoefm(1:knon, 1) = y_cd_m(1:knon)  
         ycoefh(1:knon, 1) = y_cd_h(1:knon)  
         ycoefm(1:knon, 2:klev) = ykmm(1:knon, 2:klev)  
         ycoefh(1:knon, 2:klev) = ykmn(1:knon, 2:klev)  
   
   
      END IF  
   
      !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
      ! calculer la diffusion des vitesses "u" et "v"  
      !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
   
      CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yu, ypaprs, ypplay, &  
           ydelp, y_d_u, y_flux_u)  
      CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yv, ypaprs, ypplay, &  
           ydelp, y_d_v, y_flux_v)  
   
      ! pour le couplage  
      ytaux = y_flux_u(:, 1)  
      ytauy = y_flux_v(:, 1)  
   
      ! FH modif sur le cdrag temperature  
      !$$$PB : déplace dans clcdrag  
      !$$$      do i=1, knon  
      !$$$         ycoefh(i, 1)=ycoefm(i, 1)*0.8  
      !$$$      enddo  
   
      ! calculer la diffusion de "q" et de "h"  
      CALL clqh(dtime, itap, date0, jour, debut, lafin, rlon, rlat,&  
           cufi, cvfi, knon, nsrf, ni, pctsrf, soil_model, ytsoil,&  
           yqsol, ok_veget, ocean, npas, nexca, rmu0, co2_ppm, yrugos,&  
           yrugoro, yu1, yv1, ycoefh, yt, yq, yts, ypaprs, ypplay,&  
           ydelp, yrads, yalb, yalblw, ysnow, yqsurf, yrain_f, ysnow_f, &  
           yfder, ytaux, ytauy, ywindsp, ysollw, ysollwdown, ysolsw,&  
           yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, y_d_ts,&  
           yz0_new, y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,&  
           y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, y_flux_g,&  
           ytslab, y_seaice)  
   
      ! calculer la longueur de rugosite sur ocean  
      yrugm = 0.  
      IF (nsrf==is_oce) THEN  
         DO j = 1, knon  
            yrugm(j) = 0.018*ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &  
                 0.11*14E-6/sqrt(ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2))  
            yrugm(j) = max(1.5E-05, yrugm(j))  
         END DO  
      END IF  
      DO j = 1, knon  
         y_dflux_t(j) = y_dflux_t(j)*ypct(j)  
         y_dflux_q(j) = y_dflux_q(j)*ypct(j)  
         yu1(j) = yu1(j)*ypct(j)  
         yv1(j) = yv1(j)*ypct(j)  
      END DO  
   
      DO k = 1, klev  
         DO j = 1, knon  
            i = ni(j)  
            ycoefh(j, k) = ycoefh(j, k)*ypct(j)  
            ycoefm(j, k) = ycoefm(j, k)*ypct(j)  
            y_d_t(j, k) = y_d_t(j, k)*ypct(j)  
            y_d_q(j, k) = y_d_q(j, k)*ypct(j)  
            !§§§ PB  
            flux_t(i, k, nsrf) = y_flux_t(j, k)  
            flux_q(i, k, nsrf) = y_flux_q(j, k)  
            flux_u(i, k, nsrf) = y_flux_u(j, k)  
            flux_v(i, k, nsrf) = y_flux_v(j, k)  
            !$$$ PB        y_flux_t(j, k) = y_flux_t(j, k) * ypct(j)  
            !$$$ PB        y_flux_q(j, k) = y_flux_q(j, k) * ypct(j)  
            y_d_u(j, k) = y_d_u(j, k)*ypct(j)  
            y_d_v(j, k) = y_d_v(j, k)*ypct(j)  
            !$$$ PB        y_flux_u(j, k) = y_flux_u(j, k) * ypct(j)  
            !$$$ PB        y_flux_v(j, k) = y_flux_v(j, k) * ypct(j)  
         END DO  
      END DO  
   
   
      evap(:, nsrf) = -flux_q(:, 1, nsrf)  
   
      albe(:, nsrf) = 0.  
      alblw(:, nsrf) = 0.  
      snow(:, nsrf) = 0.  
      qsurf(:, nsrf) = 0.  
      rugos(:, nsrf) = 0.  
      fluxlat(:, nsrf) = 0.  
      DO j = 1, knon  
         i = ni(j)  
         d_ts(i, nsrf) = y_d_ts(j)  
         albe(i, nsrf) = yalb(j)  
         alblw(i, nsrf) = yalblw(j)  
         snow(i, nsrf) = ysnow(j)  
         qsurf(i, nsrf) = yqsurf(j)  
         rugos(i, nsrf) = yz0_new(j)  
         fluxlat(i, nsrf) = yfluxlat(j)  
         !$$$ pb         rugmer(i) = yrugm(j)  
         IF (nsrf==is_oce) THEN  
            rugmer(i) = yrugm(j)  
            rugos(i, nsrf) = yrugm(j)  
         END IF  
         !IM cf JLD ??  
         agesno(i, nsrf) = yagesno(j)  
         fqcalving(i, nsrf) = y_fqcalving(j)  
         ffonte(i, nsrf) = y_ffonte(j)  
         cdragh(i) = cdragh(i) + ycoefh(j, 1)  
         cdragm(i) = cdragm(i) + ycoefm(j, 1)  
         dflux_t(i) = dflux_t(i) + y_dflux_t(j)  
         dflux_q(i) = dflux_q(i) + y_dflux_q(j)  
         zu1(i) = zu1(i) + yu1(j)  
         zv1(i) = zv1(i) + yv1(j)  
      END DO  
      IF (nsrf==is_ter) THEN  
         DO j = 1, knon  
            i = ni(j)  
            qsol(i) = yqsol(j)  
         END DO  
      END IF  
      IF (nsrf==is_lic) THEN  
         DO j = 1, knon  
            i = ni(j)  
            run_off_lic_0(i) = y_run_off_lic_0(j)  
         END DO  
      END IF  
      !$$$ PB ajout pour soil  
      ftsoil(:, :, nsrf) = 0.  
      DO k = 1, nsoilmx  
         DO j = 1, knon  
            i = ni(j)  
            ftsoil(i, k, nsrf) = ytsoil(j, k)  
         END DO  
      END DO  
   
      DO j = 1, knon  
         i = ni(j)  
         DO k = 1, klev  
            d_t(i, k) = d_t(i, k) + y_d_t(j, k)  
            d_q(i, k) = d_q(i, k) + y_d_q(j, k)  
            !$$$ PB        flux_t(i, k) = flux_t(i, k) + y_flux_t(j, k)  
            !$$$         flux_q(i, k) = flux_q(i, k) + y_flux_q(j, k)  
            d_u(i, k) = d_u(i, k) + y_d_u(j, k)  
            d_v(i, k) = d_v(i, k) + y_d_v(j, k)  
            !$$$  PB       flux_u(i, k) = flux_u(i, k) + y_flux_u(j, k)  
            !$$$         flux_v(i, k) = flux_v(i, k) + y_flux_v(j, k)  
            zcoefh(i, k) = zcoefh(i, k) + ycoefh(j, k)  
         END DO  
      END DO  
   
   
      !cc diagnostic t, q a 2m et u, v a 10m  
   
      DO j = 1, knon  
         i = ni(j)  
         uzon(j) = yu(j, 1) + y_d_u(j, 1)  
         vmer(j) = yv(j, 1) + y_d_v(j, 1)  
         tair1(j) = yt(j, 1) + y_d_t(j, 1)  
         qair1(j) = yq(j, 1) + y_d_q(j, 1)  
         zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &  
              1)))*(ypaprs(j, 1)-ypplay(j, 1))  
         tairsol(j) = yts(j) + y_d_ts(j)  
         rugo1(j) = yrugos(j)  
         IF (nsrf==is_oce) THEN  
            rugo1(j) = rugos(i, nsrf)  
         END IF  
         psfce(j) = ypaprs(j, 1)  
         patm(j) = ypplay(j, 1)  
   
         qairsol(j) = yqsurf(j)  
      END DO  
   
      CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, zgeo1, &  
           tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, yq10m, &  
           yu10m, yustar)  
      !IM 081204 END  
   
      DO j = 1, knon  
         i = ni(j)  
         t2m(i, nsrf) = yt2m(j)  
         q2m(i, nsrf) = yq2m(j)  
   
         ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman  
         u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)  
         v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)  
   
      END DO  
   
      !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        ! Author: Z.X. Li (LMD/CNRS), date: 1993/08/18
20        ! Objet : interface de "couche limite" (diffusion verticale)
21    
22        ! Tout ce qui a trait aux traceurs est dans "phytrac" maintenant.
23        ! Pour l'instant le calcul de la couche limite pour les traceurs
24        ! se fait avec "cltrac" et ne tient pas compte de la différentiation
25        ! des sous-fractions de sol.
26    
27        ! Pour pouvoir extraire les coefficients d'échanges et le vent
28        ! dans la première couche, trois champs supplémentaires ont été
29        ! créés : "zcoefh", "zu1" et "zv1". Pour l'instant nous avons
30        ! moyenné les valeurs de ces trois champs sur les 4 sous-surfaces
31        ! du modèle. Dans l'avenir, si les informations des sous-surfaces
32        ! doivent être prises en compte, il faudra sortir ces mêmes champs
33        ! en leur ajoutant une dimension, c'est-à-dire "nbsrf" (nombre de
34        ! sous-surfaces).
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    
225        REAL pctsrf_pot(klon, nbsrf)
226        ! "pourcentage potentiel" pour tenir compte des éventuelles
227        ! apparitions ou disparitions de la glace de mer
228    
229        REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
230    
231        ! maf pour sorties IOISPL en cas de debugagage
232    
233        CHARACTER (80) cldebug
234        SAVE cldebug
235        CHARACTER (8) cl_surf(nbsrf)
236        SAVE cl_surf
237        INTEGER nhoridbg, nidbg
238        SAVE nhoridbg, nidbg
239        INTEGER ndexbg(iim*(jjm+1))
240        REAL zx_lon(iim, jjm+1), zx_lat(iim, jjm+1), zjulian
241        REAL tabindx(klon)
242        REAL debugtab(iim, jjm+1)
243        LOGICAL first_appel
244        SAVE first_appel
245        DATA first_appel/ .TRUE./
246        LOGICAL :: debugindex = .FALSE.
247        INTEGER idayref
248        REAL t2m(klon, nbsrf), q2m(klon, nbsrf)
249        REAL u10m(klon, nbsrf), v10m(klon, nbsrf)
250    
251        REAL yt2m(klon), yq2m(klon), yu10m(klon)
252        REAL yustar(klon)
253        ! -- LOOP
254        REAL yu10mx(klon)
255        REAL yu10my(klon)
256        REAL ywindsp(klon)
257        ! -- LOOP
258    
259        REAL yt10m(klon), yq10m(klon)
260        !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds
261        ! physiq ce qui permet de sortir les grdeurs par sous surface)
262        REAL pblh(klon, nbsrf)
263        REAL plcl(klon, nbsrf)
264        REAL capcl(klon, nbsrf)
265        REAL oliqcl(klon, nbsrf)
266        REAL cteicl(klon, nbsrf)
267        REAL pblt(klon, nbsrf)
268        REAL therm(klon, nbsrf)
269        REAL trmb1(klon, nbsrf)
270        REAL trmb2(klon, nbsrf)
271        REAL trmb3(klon, nbsrf)
272        REAL ypblh(klon)
273        REAL ylcl(klon)
274        REAL ycapcl(klon)
275        REAL yoliqcl(klon)
276        REAL ycteicl(klon)
277        REAL ypblt(klon)
278        REAL ytherm(klon)
279        REAL ytrmb1(klon)
280        REAL ytrmb2(klon)
281        REAL ytrmb3(klon)
282        REAL y_cd_h(klon), y_cd_m(klon)
283        REAL uzon(klon), vmer(klon)
284        REAL tair1(klon), qair1(klon), tairsol(klon)
285        REAL psfce(klon), patm(klon)
286    
287        REAL qairsol(klon), zgeo1(klon)
288        REAL rugo1(klon)
289    
290        ! utiliser un jeu de fonctions simples              
291        LOGICAL zxli
292        PARAMETER (zxli=.FALSE.)
293    
294        REAL zt, zqs, zdelta, zcor
295        REAL t_coup
296        PARAMETER (t_coup=273.15)
297    
298        CHARACTER (len=20) :: modname = 'clmain'
299    
300        !------------------------------------------------------------
301    
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., 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        ! Initialization:
346        rugmer = 0.
347        cdragh = 0.
348        cdragm = 0.
349        dflux_t = 0.
350        dflux_q = 0.
351        zu1 = 0.
352        zv1 = 0.
353        ypct = 0.
354        yts = 0.
355        ysnow = 0.
356        yqsurf = 0.
357        yalb = 0.
358        yalblw = 0.
359        yrain_f = 0.
360        ysnow_f = 0.
361        yfder = 0.
362        ytaux = 0.
363        ytauy = 0.
364        ysolsw = 0.
365        ysollw = 0.
366        ysollwdown = 0.
367        yrugos = 0.
368        yu1 = 0.
369        yv1 = 0.
370        yrads = 0.
371        ypaprs = 0.
372        ypplay = 0.
373        ydelp = 0.
374        yu = 0.
375        yv = 0.
376        yt = 0.
377        yq = 0.
378        pctsrf_new = 0.
379        y_flux_u = 0.
380        y_flux_v = 0.
381        !$$ PB
382        y_dflux_t = 0.
383        y_dflux_q = 0.
384        ytsoil = 999999.
385        yrugoro = 0.
386        ! -- LOOP
387        yu10mx = 0.
388        yu10my = 0.
389        ywindsp = 0.
390        ! -- LOOP
391        d_ts = 0.
392        !§§§ PB
393        yfluxlat = 0.
394        flux_t = 0.
395        flux_q = 0.
396        flux_u = 0.
397        flux_v = 0.
398        d_t = 0.
399        d_q = 0.
400        d_u = 0.
401        d_v = 0.
402        zcoefh = 0.
403    
404        ! Boucler sur toutes les sous-fractions du sol:
405    
406        ! Initialisation des "pourcentages potentiels". On considère ici qu'on
407        ! peut avoir potentiellement de la glace sur tout le domaine océanique
408        ! (à affiner)
409    
410        pctsrf_pot = pctsrf
411        pctsrf_pot(:, is_oce) = 1. - zmasq
412        pctsrf_pot(:, is_sic) = 1. - zmasq
413    
414        DO nsrf = 1, nbsrf
415           ! chercher les indices:
416           ni = 0
417           knon = 0
418           DO i = 1, klon
419              ! Pour déterminer le domaine à traiter, on utilise les surfaces
420              ! "potentielles"
421              IF (pctsrf_pot(i, nsrf) > epsfra) THEN
422                 knon = knon + 1
423                 ni(knon) = i
424              END IF
425           END DO
426    
427           ! variables pour avoir une sortie IOIPSL des INDEX
428           IF (debugindex) THEN
429              tabindx = 0.
430              DO i = 1, knon
431                 tabindx(i) = real(i)
432              END DO
433              debugtab = 0.
434              ndexbg = 0
435              CALL gath2cpl(tabindx, debugtab, klon, knon, iim, jjm, ni)
436              CALL histwrite(nidbg, cl_surf(nsrf), itap, debugtab)
437           END IF
438    
439           IF (knon==0) CYCLE
440    
441           DO j = 1, knon
442              i = ni(j)
443              ypct(j) = pctsrf(i, nsrf)
444              yts(j) = ts(i, nsrf)
445              ytslab(i) = tslab(i)
446              ysnow(j) = snow(i, nsrf)
447              yqsurf(j) = qsurf(i, nsrf)
448              yalb(j) = albe(i, nsrf)
449              yalblw(j) = alblw(i, nsrf)
450              yrain_f(j) = rain_f(i)
451              ysnow_f(j) = snow_f(i)
452              yagesno(j) = agesno(i, nsrf)
453              yfder(j) = fder(i)
454              ytaux(j) = flux_u(i, 1, nsrf)
455              ytauy(j) = flux_v(i, 1, nsrf)
456              ysolsw(j) = solsw(i, nsrf)
457              ysollw(j) = sollw(i, nsrf)
458              ysollwdown(j) = sollwdown(i)
459              yrugos(j) = rugos(i, nsrf)
460              yrugoro(j) = rugoro(i)
461              yu1(j) = u1lay(i)
462              yv1(j) = v1lay(i)
463              yrads(j) = ysolsw(j) + ysollw(j)
464              ypaprs(j, klev+1) = paprs(i, klev+1)
465              y_run_off_lic_0(j) = run_off_lic_0(i)
466              yu10mx(j) = u10m(i, nsrf)
467              yu10my(j) = v10m(i, nsrf)
468              ywindsp(j) = sqrt(yu10mx(j)*yu10mx(j)+yu10my(j)*yu10my(j))
469           END DO
470    
471           !     IF bucket model for continent, copy soil water content
472           IF (nsrf==is_ter .AND. .NOT. ok_veget) THEN
473              DO j = 1, knon
474                 i = ni(j)
475                 yqsol(j) = qsol(i)
476              END DO
477           ELSE
478              yqsol = 0.
479           END IF
480           !$$$ PB ajour pour soil
481           DO k = 1, nsoilmx
482              DO j = 1, knon
483                 i = ni(j)
484                 ytsoil(j, k) = ftsoil(i, k, nsrf)
485              END DO
486           END DO
487           DO k = 1, klev
488              DO j = 1, knon
489                 i = ni(j)
490                 ypaprs(j, k) = paprs(i, k)
491                 ypplay(j, k) = pplay(i, k)
492                 ydelp(j, k) = delp(i, k)
493                 yu(j, k) = u(i, k)
494                 yv(j, k) = v(i, k)
495                 yt(j, k) = t(i, k)
496                 yq(j, k) = q(i, k)
497              END DO
498           END DO
499    
500           ! calculer Cdrag et les coefficients d'echange
501           CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts,&
502                yrugos, yu, yv, yt, yq, yqsurf, ycoefm, ycoefh)
503           !IM 081204 BEG
504           !CR test
505           IF (iflag_pbl==1) THEN
506              !IM 081204 END
507              CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
508              DO k = 1, klev
509                 DO i = 1, knon
510                    ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))
511                    ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))
512                 END DO
513              END DO
514           END IF
515    
516           !IM cf JLD : on seuille ycoefm et ycoefh
517           IF (nsrf==is_oce) THEN
518              DO j = 1, knon
519                 !           ycoefm(j, 1)=min(ycoefm(j, 1), 1.1E-3)
520                 ycoefm(j, 1) = min(ycoefm(j, 1), cdmmax)
521                 !           ycoefh(j, 1)=min(ycoefh(j, 1), 1.1E-3)
522                 ycoefh(j, 1) = min(ycoefh(j, 1), cdhmax)
523              END DO
524           END IF
525    
526           !IM: 261103
527           IF (ok_kzmin) THEN
528              !IM cf FH: 201103 BEG
529              !   Calcul d'une diffusion minimale pour les conditions tres stables.
530              CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, ycoefm, &
531                   ycoefm0, ycoefh0)
532    
533              IF (1==1) THEN
534                 DO k = 1, klev
535                    DO i = 1, knon
536                       ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))
537                       ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))
538                    END DO
539                 END DO
540              END IF
541              !IM cf FH: 201103 END
542              !IM: 261103
543           END IF !ok_kzmin
544    
545           IF (iflag_pbl>=3) THEN
546              ! MELLOR ET YAMADA adapté à Mars, Richard Fournier et Frédéric Hourdin
547              yzlay(1:knon, 1) = rd*yt(1:knon, 1)/(0.5*(ypaprs(1:knon, &
548                   1)+ypplay(1:knon, 1)))*(ypaprs(1:knon, 1)-ypplay(1:knon, 1))/rg
549              DO k = 2, klev
550                 yzlay(1:knon, k) = yzlay(1:knon, k-1) &
551                      + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
552                      / ypaprs(1:knon, k) &
553                      * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
554              END DO
555              DO k = 1, klev
556                 yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
557                      / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))
558              END DO
559              yzlev(1:knon, 1) = 0.
560              yzlev(1:knon, klev+1) = 2.*yzlay(1:knon, klev) - yzlay(1:knon, klev-1)
561              DO k = 2, klev
562                 yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
563              END DO
564              DO k = 1, klev + 1
565                 DO j = 1, knon
566                    i = ni(j)
567                    yq2(j, k) = q2(i, k, nsrf)
568                 END DO
569              END DO
570    
571              !   Bug introduit volontairement pour converger avec les resultats
572              !  du papier sur les thermiques.
573              IF (1==1) THEN
574                 y_cd_m(1:knon) = ycoefm(1:knon, 1)
575                 y_cd_h(1:knon) = ycoefh(1:knon, 1)
576              ELSE
577                 y_cd_h(1:knon) = ycoefm(1:knon, 1)
578                 y_cd_m(1:knon) = ycoefh(1:knon, 1)
579              END IF
580              CALL ustarhb(knon, yu, yv, y_cd_m, yustar)
581    
582              IF (prt_level>9) THEN
583                 PRINT *, 'USTAR = ', yustar
584              END IF
585    
586              !   iflag_pbl peut etre utilise comme longuer de melange
587    
588              IF (iflag_pbl>=11) THEN
589                 CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &
590                      yu, yv, yteta, y_cd_m, yq2, q2diag, ykmm, ykmn, yustar, &
591                      iflag_pbl)
592              ELSE
593                 CALL yamada4(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, yu, &
594                      yv, yteta, y_cd_m, yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
595              END IF
596    
597              ycoefm(1:knon, 1) = y_cd_m(1:knon)
598              ycoefh(1:knon, 1) = y_cd_h(1:knon)
599              ycoefm(1:knon, 2:klev) = ykmm(1:knon, 2:klev)
600              ycoefh(1:knon, 2:klev) = ykmn(1:knon, 2:klev)
601           END IF
602    
603           ! calculer la diffusion des vitesses "u" et "v"
604           CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yu, ypaprs, ypplay, &
605                ydelp, y_d_u, y_flux_u)
606           CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yv, ypaprs, ypplay, &
607                ydelp, y_d_v, y_flux_v)
608    
609           ! pour le couplage
610           ytaux = y_flux_u(:, 1)
611           ytauy = y_flux_v(:, 1)
612    
613           ! calculer la diffusion de "q" et de "h"
614           CALL clqh(dtime, itap, date0, jour, debut, lafin, rlon, rlat,&
615                cufi, cvfi, knon, nsrf, ni, pctsrf, soil_model, ytsoil,&
616                yqsol, ok_veget, ocean, npas, nexca, rmu0, co2_ppm, yrugos,&
617                yrugoro, yu1, yv1, ycoefh, yt, yq, yts, ypaprs, ypplay,&
618                ydelp, yrads, yalb, yalblw, ysnow, yqsurf, yrain_f, ysnow_f, &
619                yfder, ytaux, ytauy, ywindsp, ysollw, ysollwdown, ysolsw,&
620                yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, y_d_ts,&
621                yz0_new, y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,&
622                y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, y_flux_g,&
623                ytslab, y_seaice)
624    
625           ! calculer la longueur de rugosite sur ocean
626           yrugm = 0.
627           IF (nsrf==is_oce) THEN
628              DO j = 1, knon
629                 yrugm(j) = 0.018*ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &
630                      0.11*14E-6/sqrt(ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2))
631                 yrugm(j) = max(1.5E-05, yrugm(j))
632              END DO
633           END IF
634           DO j = 1, knon
635              y_dflux_t(j) = y_dflux_t(j)*ypct(j)
636              y_dflux_q(j) = y_dflux_q(j)*ypct(j)
637              yu1(j) = yu1(j)*ypct(j)
638              yv1(j) = yv1(j)*ypct(j)
639           END DO
640    
641           DO k = 1, klev
642              DO j = 1, knon
643                 i = ni(j)
644                 ycoefh(j, k) = ycoefh(j, k)*ypct(j)
645                 ycoefm(j, k) = ycoefm(j, k)*ypct(j)
646                 y_d_t(j, k) = y_d_t(j, k)*ypct(j)
647                 y_d_q(j, k) = y_d_q(j, k)*ypct(j)
648                 !§§§ PB
649                 flux_t(i, k, nsrf) = y_flux_t(j, k)
650                 flux_q(i, k, nsrf) = y_flux_q(j, k)
651                 flux_u(i, k, nsrf) = y_flux_u(j, k)
652                 flux_v(i, k, nsrf) = y_flux_v(j, k)
653                 !$$$ PB        y_flux_t(j, k) = y_flux_t(j, k) * ypct(j)
654                 !$$$ PB        y_flux_q(j, k) = y_flux_q(j, k) * ypct(j)
655                 y_d_u(j, k) = y_d_u(j, k)*ypct(j)
656                 y_d_v(j, k) = y_d_v(j, k)*ypct(j)
657                 !$$$ PB        y_flux_u(j, k) = y_flux_u(j, k) * ypct(j)
658                 !$$$ PB        y_flux_v(j, k) = y_flux_v(j, k) * ypct(j)
659              END DO
660           END DO
661    
662           evap(:, nsrf) = -flux_q(:, 1, nsrf)
663    
664           albe(:, nsrf) = 0.
665           alblw(:, nsrf) = 0.
666           snow(:, nsrf) = 0.
667           qsurf(:, nsrf) = 0.
668           rugos(:, nsrf) = 0.
669           fluxlat(:, nsrf) = 0.
670           DO j = 1, knon
671              i = ni(j)
672              d_ts(i, nsrf) = y_d_ts(j)
673              albe(i, nsrf) = yalb(j)
674              alblw(i, nsrf) = yalblw(j)
675              snow(i, nsrf) = ysnow(j)
676              qsurf(i, nsrf) = yqsurf(j)
677              rugos(i, nsrf) = yz0_new(j)
678              fluxlat(i, nsrf) = yfluxlat(j)
679              !$$$ pb         rugmer(i) = yrugm(j)
680              IF (nsrf==is_oce) THEN
681                 rugmer(i) = yrugm(j)
682                 rugos(i, nsrf) = yrugm(j)
683              END IF
684              !IM cf JLD ??
685              agesno(i, nsrf) = yagesno(j)
686              fqcalving(i, nsrf) = y_fqcalving(j)
687              ffonte(i, nsrf) = y_ffonte(j)
688              cdragh(i) = cdragh(i) + ycoefh(j, 1)
689              cdragm(i) = cdragm(i) + ycoefm(j, 1)
690              dflux_t(i) = dflux_t(i) + y_dflux_t(j)
691              dflux_q(i) = dflux_q(i) + y_dflux_q(j)
692              zu1(i) = zu1(i) + yu1(j)
693              zv1(i) = zv1(i) + yv1(j)
694           END DO
695           IF (nsrf==is_ter) THEN
696              DO j = 1, knon
697                 i = ni(j)
698                 qsol(i) = yqsol(j)
699              END DO
700           END IF
701           IF (nsrf==is_lic) THEN
702              DO j = 1, knon
703                 i = ni(j)
704                 run_off_lic_0(i) = y_run_off_lic_0(j)
705              END DO
706           END IF
707           !$$$ PB ajout pour soil
708           ftsoil(:, :, nsrf) = 0.
709           DO k = 1, nsoilmx
710              DO j = 1, knon
711                 i = ni(j)
712                 ftsoil(i, k, nsrf) = ytsoil(j, k)
713              END DO
714           END DO
715    
716           DO j = 1, knon
717              i = ni(j)
718              DO k = 1, klev
719                 d_t(i, k) = d_t(i, k) + y_d_t(j, k)
720                 d_q(i, k) = d_q(i, k) + y_d_q(j, k)
721                 !$$$ PB        flux_t(i, k) = flux_t(i, k) + y_flux_t(j, k)
722                 !$$$         flux_q(i, k) = flux_q(i, k) + y_flux_q(j, k)
723                 d_u(i, k) = d_u(i, k) + y_d_u(j, k)
724                 d_v(i, k) = d_v(i, k) + y_d_v(j, k)
725                 !$$$  PB       flux_u(i, k) = flux_u(i, k) + y_flux_u(j, k)
726                 !$$$         flux_v(i, k) = flux_v(i, k) + y_flux_v(j, k)
727                 zcoefh(i, k) = zcoefh(i, k) + ycoefh(j, k)
728              END DO
729           END DO
730    
731           !cc diagnostic t, q a 2m et u, v a 10m
732    
733           DO j = 1, knon
734              i = ni(j)
735              uzon(j) = yu(j, 1) + y_d_u(j, 1)
736              vmer(j) = yv(j, 1) + y_d_v(j, 1)
737              tair1(j) = yt(j, 1) + y_d_t(j, 1)
738              qair1(j) = yq(j, 1) + y_d_q(j, 1)
739              zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &
740                   1)))*(ypaprs(j, 1)-ypplay(j, 1))
741              tairsol(j) = yts(j) + y_d_ts(j)
742              rugo1(j) = yrugos(j)
743              IF (nsrf==is_oce) THEN
744                 rugo1(j) = rugos(i, nsrf)
745              END IF
746              psfce(j) = ypaprs(j, 1)
747              patm(j) = ypplay(j, 1)
748    
749              qairsol(j) = yqsurf(j)
750           END DO
751    
752           CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, zgeo1, &
753                tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, yq10m, &
754                yu10m, yustar)
755           !IM 081204 END
756    
757           DO j = 1, knon
758              i = ni(j)
759              t2m(i, nsrf) = yt2m(j)
760              q2m(i, nsrf) = yq2m(j)
761    
762              ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
763              u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)
764              v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)
765    
766           END DO
767    
768           DO i = 1, knon
769              y_cd_h(i) = ycoefh(i, 1)
770              y_cd_m(i) = ycoefm(i, 1)
771           END DO
772           CALL hbtm(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, &
773                y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &
774                ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)
775    
776           DO j = 1, knon
777              i = ni(j)
778              pblh(i, nsrf) = ypblh(j)
779              plcl(i, nsrf) = ylcl(j)
780              capcl(i, nsrf) = ycapcl(j)
781              oliqcl(i, nsrf) = yoliqcl(j)
782              cteicl(i, nsrf) = ycteicl(j)
783              pblt(i, nsrf) = ypblt(j)
784              therm(i, nsrf) = ytherm(j)
785              trmb1(i, nsrf) = ytrmb1(j)
786              trmb2(i, nsrf) = ytrmb2(j)
787              trmb3(i, nsrf) = ytrmb3(j)
788           END DO
789    
790           DO j = 1, knon
791              DO k = 1, klev + 1
792                 i = ni(j)
793                 q2(i, k, nsrf) = yq2(j, k)
794              END DO
795           END DO
796           !IM "slab" ocean
797           IF (nsrf==is_oce) THEN
798              DO j = 1, knon
799                 ! on projette sur la grille globale
800                 i = ni(j)
801                 IF (pctsrf_new(i, is_oce)>epsfra) THEN
802                    flux_o(i) = y_flux_o(j)
803                 ELSE
804                    flux_o(i) = 0.
805                 END IF
806              END DO
807           END IF
808    
809           IF (nsrf==is_sic) THEN
810              DO j = 1, knon
811                 i = ni(j)
812                 ! On pondère lorsque l'on fait le bilan au sol :
813                 ! flux_g(i) = y_flux_g(j)*ypct(j)
814                 IF (pctsrf_new(i, is_sic)>epsfra) THEN
815                    flux_g(i) = y_flux_g(j)
816                 ELSE
817                    flux_g(i) = 0.
818                 END IF
819              END DO
820    
821           END IF
822           !nsrf.EQ.is_sic                                            
823           IF (ocean=='slab  ') THEN
824              IF (nsrf==is_oce) THEN
825                 tslab(1:klon) = ytslab(1:klon)
826                 seaice(1:klon) = y_seaice(1:klon)
827                 !nsrf                                                      
828              END IF
829              !OCEAN                                                      
830           END IF
831        END DO
832    
833        ! On utilise les nouvelles surfaces
834        ! A rajouter: conservation de l'albedo
835    
836    rugos(:, is_oce) = rugmer      rugos(:, is_oce) = rugmer
837    pctsrf = pctsrf_new      pctsrf = pctsrf_new
838    
839  END SUBROUTINE clmain    END SUBROUTINE clmain
840    
841    end module clmain_m

Legend:
Removed from v.30  
changed lines
  Added in v.40

  ViewVC Help
Powered by ViewVC 1.1.21