Index: LMDZ6/branches/Amaury_dev/libf/phylmdiso/add_phys_tend_mod.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmdiso/add_phys_tend_mod.F90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmdiso/add_phys_tend_mod.F90	(revision 5144)
@@ -161,6 +161,8 @@
 #endif  
 #endif
+
+USE lmdz_yomcst
+
 IMPLICIT NONE
-  include "YOMCST.h"
 
 ! Arguments :
@@ -685,7 +687,7 @@
               , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_qbs_col, d_h_col
 USE lmdz_clesphys
+USE lmdz_yomcst
 
 IMPLICIT NONE
-  include "YOMCST.h"
 
 ! Arguments :
@@ -835,6 +837,6 @@
                     zh_qw_col, zh_ql_col, zh_qs_col, zh_qbs_col, zh_col)
 
+  USE lmdz_yomcst
 IMPLICIT NONE
-  include "YOMCST.h"
 
 INTEGER,                    INTENT(IN)    :: nlon,nlev
@@ -916,6 +918,7 @@
       , rain_lsc, snow_lsc
 USE climb_hq_mod, ONLY: d_h_col_vdf, f_h_bnd
+USE lmdz_yomcst
+
 IMPLICIT NONE
-include "YOMCST.h"
 
 ! Arguments :
Index: LMDZ6/branches/Amaury_dev/libf/phylmdiso/ajsec.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmdiso/ajsec.F90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmdiso/ajsec.F90	(revision 5144)
@@ -19,4 +19,5 @@
 #endif
 #endif
+USE lmdz_yomcst
   IMPLICIT NONE
   ! ======================================================================
@@ -29,5 +30,4 @@
   ! d_t-----output-R-Incrementation de la temperature
   ! ======================================================================
-  include "YOMCST.h"
   REAL paprs(klon, klev+1), pplay(klon, klev)
   REAL t(klon, klev), q(klon, klev)
@@ -313,4 +313,6 @@
 #endif
 #endif
+        USE lmdz_yomcst
+
   IMPLICIT NONE
   ! ======================================================================
@@ -323,5 +325,4 @@
   ! d_t-----output-R-Incrementation de la temperature
   ! ======================================================================
-  include "YOMCST.h"
   REAL paprs(klon, klev+1), pplay(klon, klev)
   REAL t(klon, klev), q(klon, klev)
@@ -579,4 +580,6 @@
 SUBROUTINE ajsec_old(paprs, pplay, t, d_t)
   USE dimphy
+        USE lmdz_yomcst
+
   IMPLICIT NONE
   ! ======================================================================
@@ -589,5 +592,4 @@
   ! d_t-----output-R-Incrementation de la temperature
   ! ======================================================================
-  include "YOMCST.h"
   REAL paprs(klon, klev+1), pplay(klon, klev)
   REAL t(klon, klev)
Index: LMDZ6/branches/Amaury_dev/libf/phylmdiso/calwake.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmdiso/calwake.F90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmdiso/calwake.F90	(revision 5144)
@@ -42,7 +42,8 @@
 #endif
 #endif
+USE lmdz_yomcst
+
   IMPLICIT NONE
   ! ======================================================================
-  include "YOMCST.h"
 
   ! Arguments
Index: LMDZ6/branches/Amaury_dev/libf/phylmdiso/change_srf_frac_mod.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmdiso/change_srf_frac_mod.F90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmdiso/change_srf_frac_mod.F90	(revision 5144)
@@ -38,7 +38,8 @@
 #endif
   USE lmdz_clesphys ! albedo SB
+  USE lmdz_yomcst
     
-    INCLUDE "YOMCST.h"
-
+
+  IMPLICIT NONE
 
 
Index: LMDZ6/branches/Amaury_dev/libf/phylmdiso/clift.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmdiso/clift.F90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmdiso/clift.F90	(revision 5144)
@@ -1,7 +1,8 @@
-
 ! $Header$
 
 SUBROUTINE clift(p, t, rr, rs, plcl, dplcldt, dplcldq)
-IMPLICIT NONE
+  USE lmdz_yomcst
+
+  IMPLICIT NONE
   ! ***************************************************************
   ! *                                                             *
@@ -41,7 +42,6 @@
   ! on utilise les constantes thermo du Centre Europeen: (sb)
 
-  include "YOMCST.h"
-  REAL :: p,t,rr,rs,plcl,dplcldt,dplcldq,cpd,cpv,cl,cpvmcl,eps,alv0,a,b
-  REAL :: rh,chi,alv
+  REAL :: p, t, rr, rs, plcl, dplcldt, dplcldq, cpd, cpv, cl, cpvmcl, eps, alv0, a, b
+  REAL :: rh, chi, alv
 
   cpd = rcpd
@@ -49,5 +49,5 @@
   cl = rcw
   cpvmcl = cl - cpv
-  eps = rd/rv
+  eps = rd / rv
   alv0 = rlvtt
 
@@ -57,18 +57,16 @@
   b = 122.0
 
-  rh = rr/rs
-  chi = t/(a-b*rh-t)
-  plcl = p*(rh**chi)
+  rh = rr / rs
+  chi = t / (a - b * rh - t)
+  plcl = p * (rh**chi)
 
-  alv = alv0 - cpvmcl*(t-273.15)
+  alv = alv0 - cpvmcl * (t - 273.15)
 
   ! -- sb: correction:
   ! DPLCLDQ = PLCL*CHI*( 1./RR - B*CHI/T/RS*ALOG(RH) )
-  dplcldq = plcl*chi*(1./rr+b*chi/t/rs*alog(rh))
+  dplcldq = plcl * chi * (1. / rr + b * chi / t / rs * alog(rh))
   ! sb --
 
-  dplcldt = plcl*chi*((a-b*rh*(1.+alv/rv/t))/t**2*chi*alog(rh)-alv/rv/t**2)
-
-
+  dplcldt = plcl * chi * ((a - b * rh * (1. + alv / rv / t)) / t**2 * chi * alog(rh) - alv / rv / t**2)
 
 END SUBROUTINE clift
Index: LMDZ6/branches/Amaury_dev/libf/phylmdiso/climb_hq_mod.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmdiso/climb_hq_mod.F90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmdiso/climb_hq_mod.F90	(revision 5144)
@@ -57,4 +57,5 @@
             )
   USE lmdz_compbl, ONLY: iflag_pbl, iflag_pbl_split, iflag_order2_sollw, ifl_pbltree
+  USE lmdz_yomcst
 #ifdef ISOVERIF
 USE isotopes_mod, ONLY: iso_eau,iso_HDO
@@ -126,5 +127,4 @@
 ! Include
 !****************************************************************************************
-    INCLUDE "YOMCST.h"
 
 #ifdef ISO
@@ -421,5 +421,8 @@
 ! where X is H or Q, and k the vertical level k=1,klev
 
-    INCLUDE "YOMCST.h"
+    USE lmdz_yomcst
+
+    IMPLICIT NONE
+
 ! Input arguments
 !****************************************************************************************
@@ -499,4 +502,5 @@
 
   USE lmdz_compbl, ONLY: iflag_pbl, iflag_pbl_split, iflag_order2_sollw, ifl_pbltree
+  USE lmdz_yomcst
 #ifdef ISOVERIF
 USE infotrac_phy, ONLY: nzone
@@ -552,9 +556,5 @@
 !#endif
 #endif
- 
-! Include
-!****************************************************************************************
-    INCLUDE "YOMCST.h"
-!****************************************************************************************
+
 ! 1) 
 ! Definition of some variables
Index: LMDZ6/branches/Amaury_dev/libf/phylmdiso/concvl.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmdiso/concvl.F90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmdiso/concvl.F90	(revision 5144)
@@ -43,5 +43,5 @@
   USE dimphy
   USE infotrac_phy, ONLY: nbtr
-  USE lmdz_YOETHF
+  USE lmdz_yoethf
 #ifdef ISO
   USE infotrac_phy, ONLY: ntraciso=>ntiso
@@ -67,4 +67,5 @@
   USE lmdz_conema3
   USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep
+  USE lmdz_yomcst
 
   IMPLICIT NONE
@@ -296,7 +297,5 @@
 !$OMP THREADPRIVATE(itap, igout)
 
-
-  include "YOMCST.h"
-  include "YOMCST2.h"
+ include "YOMCST2.h"
 
   IF (first) THEN
Index: LMDZ6/branches/Amaury_dev/libf/phylmdiso/cv3_enthalpmix.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmdiso/cv3_enthalpmix.F90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmdiso/cv3_enthalpmix.F90	(revision 5144)
@@ -18,6 +18,7 @@
   ! **************************************************************
 USE lmdz_cvthermo
-USE lmdz_YOETHF
+USE lmdz_yoethf
 USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep
+USE lmdz_yomcst
 
   IMPLICIT NONE
@@ -32,5 +33,4 @@
   ! ===============================================================
 
-  include "YOMCST.h"
 !inputs:
   INTEGER, INTENT (IN)                      :: nd, len
Index: LMDZ6/branches/Amaury_dev/libf/phylmdiso/cv3_estatmix.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmdiso/cv3_estatmix.F90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmdiso/cv3_estatmix.F90	(revision 5144)
@@ -19,6 +19,7 @@
   ! ****************************************************************
 USE lmdz_cvthermo
-USE lmdz_YOETHF
+USE lmdz_yoethf
 USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep
+USE lmdz_yomcst
 
   IMPLICIT NONE
@@ -33,5 +34,4 @@
   ! ===============================================================
 
-  include "YOMCST.h"
 !inputs:
   INTEGER, INTENT (IN)                      :: nd, len
Index: LMDZ6/branches/Amaury_dev/libf/phylmdiso/cv_driver.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmdiso/cv_driver.F90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmdiso/cv_driver.F90	(revision 5144)
@@ -1267,4 +1267,5 @@
 SUBROUTINE cv_thermo(iflag_con)
   USE lmdz_cvthermo
+  USE lmdz_yomcst
 
   IMPLICIT NONE
@@ -1273,6 +1274,4 @@
   ! Set thermodynamical constants for convectL
   ! -------------------------------------------------------------
-
-  include "YOMCST.h"
 
   INTEGER iflag_con
Index: LMDZ6/branches/Amaury_dev/libf/phylmdiso/isotopes_routines_mod.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmdiso/isotopes_routines_mod.F90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmdiso/isotopes_routines_mod.F90	(revision 5144)
@@ -12843,5 +12843,5 @@
     USE isotrac_mod, ONLY: izone_cont,index_zone,index_iso
 #endif
-USE lmdz_YOETHF
+USE lmdz_yoethf
 USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep
 
@@ -13407,5 +13407,5 @@
 &       bassin_map
 #endif
-USE lmdz_YOETHF
+USE lmdz_yoethf
 USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep
 
@@ -14318,5 +14318,5 @@
 &       bassin_map
 #endif
-USE lmdz_YOETHF
+USE lmdz_yoethf
 USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep
 
Index: LMDZ6/branches/Amaury_dev/libf/phylmdiso/lmdz_lscp_old.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmdiso/lmdz_lscp_old.F90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmdiso/lmdz_lscp_old.F90	(revision 5144)
@@ -68,6 +68,7 @@
 #endif
 
-USE lmdz_YOETHF
+USE lmdz_yoethf
   USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep
+  USE lmdz_yomcst
 
   IMPLICIT NONE
@@ -97,5 +98,4 @@
   !   fl_cor_ebil= 0 pour reproduire anciens bugs
   !======================================================================
-  include "YOMCST.h"
 
   ! Principaux inputs:
Index: LMDZ6/branches/Amaury_dev/libf/phylmdiso/lmdz_wake.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmdiso/lmdz_wake.F90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmdiso/lmdz_wake.F90	(revision 5144)
@@ -42,4 +42,5 @@
 #endif
 USE lmdz_cvthermo
+USE lmdz_yomcst
 
   IMPLICIT NONE
@@ -136,6 +137,4 @@
   ! Déclaration de variables
   ! -------------------------------------------------------------------------
-
-  include "YOMCST.h"
 
   ! Arguments en entree
Index: LMDZ6/branches/Amaury_dev/libf/phylmdiso/phyaqua_mod.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmdiso/phyaqua_mod.F90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmdiso/phyaqua_mod.F90	(revision 5144)
@@ -1,3 +1,2 @@
-
 ! $Id: phyaqua_mod.F90 3579 2019-10-09 13:11:07Z fairhead $
 
@@ -9,5 +8,5 @@
 CONTAINS
 
-  SUBROUTINE iniaqua(nlon,year_len,iflag_phys)
+  SUBROUTINE iniaqua(nlon, year_len, iflag_phys)
 
     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -34,26 +33,26 @@
     USE phys_state_var_mod
     USE time_phylmdz_mod, ONLY: day_ref, ndays, pdtphys, &
-                                day_ini,day_end
+            day_ini, day_end
     USE indice_sol_mod
     USE lmdz_physical_constants, ONLY: pi
-!    USE ioipsl
+    !    USE ioipsl
     USE lmdz_phys_para, ONLY: is_master
     USE lmdz_phys_transfert_para, ONLY: bcast
     USE lmdz_grid_phy
     USE lmdz_ioipsl_getin_p, ONLY: getin_p
-    USE phys_cal_mod , ONLY: calend, year_len_phy => year_len
+    USE phys_cal_mod, ONLY: calend, year_len_phy => year_len
     USE lmdz_clesphys
 #ifdef ISO
   USE infotrac_phy, ONLY: niso
 #endif
+    USE lmdz_yomcst
 
     IMPLICIT NONE
 
-    include "YOMCST.h"
     include "dimsoil.h"
 
     INTEGER, INTENT (IN) :: nlon, year_len, iflag_phys
     ! IM ajout latfi, lonfi
-!    REAL, INTENT (IN) :: lonfi(nlon), latfi(nlon)
+    !    REAL, INTENT (IN) :: lonfi(nlon), latfi(nlon)
 
     INTEGER type_profil, type_aqua
@@ -103,5 +102,5 @@
 
     INTEGER, SAVE :: read_climoz = 0 ! read ozone climatology
-!$OMP THREADPRIVATE(read_climoz)
+    !$OMP THREADPRIVATE(read_climoz)
 
     ! -------------------------------------------------------------------------
@@ -119,5 +118,5 @@
 
     INTEGER longcles
-    PARAMETER (longcles=20)
+    PARAMETER (longcles = 20)
     REAL clesphy0(longcles)
 
@@ -132,6 +131,6 @@
 
     ! Local
-    CHARACTER (LEN=20) :: modname='phyaqua'
-    CHARACTER (LEN=80) :: abort_message
+    CHARACTER (LEN = 20) :: modname = 'phyaqua'
+    CHARACTER (LEN = 80) :: abort_message
 
 
@@ -145,24 +144,23 @@
 
     !IF (calend .EQ. "earth_360d") Then
-      year_len_phy = year_len
+    year_len_phy = year_len
     !END IF
-    
+
     IF (year_len/=360) THEN
-      write (*,*) year_len
-      write (*,*) 'iniaqua: 360 day calendar is required !'
+      write (*, *) year_len
+      write (*, *) 'iniaqua: 360 day calendar is required !'
       stop
     endif
 
-    type_aqua = iflag_phys/100
-    type_profil = iflag_phys - type_aqua*100
+    type_aqua = iflag_phys / 100
+    type_profil = iflag_phys - type_aqua * 100
     PRINT *, 'iniaqua:type_aqua, type_profil', type_aqua, type_profil
 
     IF (klon/=nlon) THEN
       WRITE (*, *) 'iniaqua: klon=', klon, ' nlon=', nlon
-      abort_message= 'probleme de dimensions dans iniaqua'
-      CALL abort_physic(modname,abort_message,1)
+      abort_message = 'probleme de dimensions dans iniaqua'
+      CALL abort_physic(modname, abort_message, 1)
     END IF
     CALL phys_state_var_init(read_climoz)
-
 
     read_climoz = 0
@@ -192,5 +190,5 @@
     solaire = 1365.
     CALL getin_p('solaire', solaire)
- 
+
     ! CALL getin('albedo',albedo) ! albedo is set below, depending on
     ! type_aqua
@@ -269,9 +267,9 @@
     IF (grid_type==unstructured) THEN
       CALL writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
-                             phy_fter, phy_foce, phy_flic, phy_fsic)
+              phy_fter, phy_foce, phy_flic, phy_fsic)
     ELSE
-      
-       CALL writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
-                     phy_fter, phy_foce, phy_flic, phy_fsic)
+
+      CALL writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
+              phy_fter, phy_foce, phy_flic, phy_fsic)
     ENDIF
 
@@ -284,5 +282,5 @@
 
     timestep = pdtphys
-    radpas = nint(rday/timestep/float(nbapp_rad))
+    radpas = nint(rday / timestep / float(nbapp_rad))
 
     DO i = 1, longcles
@@ -348,6 +346,5 @@
     snsrf(:, :) = 0. ! couverture de neige des sous surface
     z0m(:, :) = rugos ! couverture de neige des sous surface
-    z0h=z0m
-
+    z0h = z0m
 
     CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil)
@@ -358,5 +355,5 @@
     PRINT *, 'iniaqua: before phyredem'
 
-    pbl_tke(:,:,:) = 1.e-8 
+    pbl_tke(:, :, :) = 1.e-8
     falb1 = albedo
     falb2 = albedo
@@ -376,33 +373,32 @@
     detr_therm = 0.
     ale_bl = 0.
-    ale_bl_trig =0. 
-    alp_bl =0.
-    treedrg(:,:,:)=0.
+    ale_bl_trig = 0.
+    alp_bl = 0.
+    treedrg(:, :, :) = 0.
 
     u10m = 0.
     v10m = 0.
 
-    ql_ancien   = 0.
-    qs_ancien   = 0.
-    u_ancien    = 0.
-    v_ancien    = 0.
-    prw_ancien  = 0.
+    ql_ancien = 0.
+    qs_ancien = 0.
+    u_ancien = 0.
+    v_ancien = 0.
+    prw_ancien = 0.
     prlw_ancien = 0.
-    prsw_ancien = 0.  
-
-    ale_wake    = 0.
-    ale_bl_stat = 0.  
-
-
-!ym error : the sub surface dimension is the third not second : forgotten for iniaqua
-!    falb_dir(:,is_ter,:)=0.08; falb_dir(:,is_lic,:)=0.6
-!    falb_dir(:,is_oce,:)=0.5;  falb_dir(:,is_sic,:)=0.6
-    falb_dir(:,:,is_ter)=0.08; falb_dir(:,:,is_lic)=0.6
-    falb_dir(:,:,is_oce)=0.5;  falb_dir(:,:,is_sic)=0.6
-
-!ym falb_dif has been forgotten, initialize with defaukt value found in phyetat0 or 0 ?
-!ym probably the uninitialized value was 0 for standard (regular grid) case
-    falb_dif(:,:,:)=0
-
+    prsw_ancien = 0.
+
+    ale_wake = 0.
+    ale_bl_stat = 0.
+
+
+    !ym error : the sub surface dimension is the third not second : forgotten for iniaqua
+    !    falb_dir(:,is_ter,:)=0.08; falb_dir(:,is_lic,:)=0.6
+    !    falb_dir(:,is_oce,:)=0.5;  falb_dir(:,is_sic,:)=0.6
+    falb_dir(:, :, is_ter) = 0.08; falb_dir(:, :, is_lic) = 0.6
+    falb_dir(:, :, is_oce) = 0.5;  falb_dir(:, :, is_sic) = 0.6
+
+    !ym falb_dif has been forgotten, initialize with defaukt value found in phyetat0 or 0 ?
+    !ym probably the uninitialized value was 0 for standard (regular grid) case
+    falb_dif(:, :, :) = 0
 
     CALL phyredem('startphy.nc')
@@ -410,5 +406,4 @@
     PRINT *, 'iniaqua: after phyredem'
     CALL phys_state_var_end
-
 
   END SUBROUTINE iniaqua
@@ -419,4 +414,6 @@
   SUBROUTINE zenang_an(cycle_diurne, gmtime, rlat, rlon, rmu0, fract)
     USE dimphy
+    USE lmdz_yomcst
+
     IMPLICIT NONE
     ! ====================================================================
@@ -450,7 +447,5 @@
     ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
     ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
-    ! ================================================================
-    include "YOMCST.h"
-    ! ================================================================
+
     LOGICAL cycle_diurne
     REAL gmtime
@@ -461,9 +456,7 @@
     REAL pi_local
 
-
     REAL rmu0m(klon), rmu0a(klon)
 
-
-    pi_local = 4.0*atan(1.0)
+    pi_local = 4.0 * atan(1.0)
 
     ! ================================================================
@@ -475,11 +468,11 @@
       ! Calcule du flux moyen
       IF (abs(rlat(i))<=28.75) THEN
-        rmu0m(i) = (210.1924+206.6059*cos(0.0174533*rlat(i))**2)/1365.
+        rmu0m(i) = (210.1924 + 206.6059 * cos(0.0174533 * rlat(i))**2) / 1365.
       ELSE IF (abs(rlat(i))<=43.75) THEN
-        rmu0m(i) = (187.4562+236.1853*cos(0.0174533*rlat(i))**2)/1365.
+        rmu0m(i) = (187.4562 + 236.1853 * cos(0.0174533 * rlat(i))**2) / 1365.
       ELSE IF (abs(rlat(i))<=71.25) THEN
-        rmu0m(i) = (162.4439+284.1192*cos(0.0174533*rlat(i))**2)/1365.
+        rmu0m(i) = (162.4439 + 284.1192 * cos(0.0174533 * rlat(i))**2) / 1365.
       ELSE
-        rmu0m(i) = (172.8125+183.7673*cos(0.0174533*rlat(i))**2)/1365.
+        rmu0m(i) = (172.8125 + 183.7673 * cos(0.0174533 * rlat(i))**2) / 1365.
       END IF
     END DO
@@ -495,7 +488,7 @@
 
       DO i = 1, klon
-        rmu0a(i) = 2.*rmu0m(i)*sqrt(2.)*pi_local/(4.-pi_local)
-        rmu0(i) = rmu0a(i)*abs(sin(pi_local*gmtime+pi_local*rlon(i)/360.)) - &
-          rmu0a(i)/sqrt(2.)
+        rmu0a(i) = 2. * rmu0m(i) * sqrt(2.) * pi_local / (4. - pi_local)
+        rmu0(i) = rmu0a(i) * abs(sin(pi_local * gmtime + pi_local * rlon(i) / 360.)) - &
+                rmu0a(i) / sqrt(2.)
       END DO
 
@@ -522,10 +515,9 @@
       DO i = 1, klon
         fract(i) = 0.5
-        rmu0(i) = rmu0m(i)*2.
+        rmu0(i) = rmu0m(i) * 2.
       END DO
 
     END IF
 
-
   END SUBROUTINE zenang_an
 
@@ -533,5 +525,5 @@
 
   SUBROUTINE writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
-      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
+          phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
 
     USE lmdz_phys_para, ONLY: is_omp_master, klon_mpi
@@ -553,38 +545,37 @@
 
     REAL :: phy_mpi(klon_mpi, 360) ! temporary variable, to store phy_***(:)
-      ! on the whole physics grid
- 
+    ! on the whole physics grid
+
     PRINT *, 'writelim: Ecriture du fichier limit'
 
     CALL gather_omp(phy_foce, phy_mpi)
-    IF (is_omp_master) CALL xios_send_field('foce_limout',phy_mpi)
+    IF (is_omp_master) CALL xios_send_field('foce_limout', phy_mpi)
 
     CALL gather_omp(phy_fsic, phy_mpi)
-    IF (is_omp_master) CALL xios_send_field('fsic_limout',phy_mpi)
-      
+    IF (is_omp_master) CALL xios_send_field('fsic_limout', phy_mpi)
+
     CALL gather_omp(phy_fter, phy_mpi)
-    IF (is_omp_master) CALL xios_send_field('fter_limout',phy_mpi)
-      
+    IF (is_omp_master) CALL xios_send_field('fter_limout', phy_mpi)
+
     CALL gather_omp(phy_flic, phy_mpi)
-    IF (is_omp_master) CALL xios_send_field('flic_limout',phy_mpi)
+    IF (is_omp_master) CALL xios_send_field('flic_limout', phy_mpi)
 
     CALL gather_omp(phy_sst, phy_mpi)
-    IF (is_omp_master) CALL xios_send_field('sst_limout',phy_mpi)
+    IF (is_omp_master) CALL xios_send_field('sst_limout', phy_mpi)
 
     CALL gather_omp(phy_bil, phy_mpi)
-    IF (is_omp_master) CALL xios_send_field('bils_limout',phy_mpi)
+    IF (is_omp_master) CALL xios_send_field('bils_limout', phy_mpi)
 
     CALL gather_omp(phy_alb, phy_mpi)
-    IF (is_omp_master) CALL xios_send_field('alb_limout',phy_mpi)
+    IF (is_omp_master) CALL xios_send_field('alb_limout', phy_mpi)
 
     CALL gather_omp(phy_rug, phy_mpi)
-    IF (is_omp_master) CALL xios_send_field('rug_limout',phy_mpi)
+    IF (is_omp_master) CALL xios_send_field('rug_limout', phy_mpi)
 
   END SUBROUTINE writelim_unstruct
 
 
-
   SUBROUTINE writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
-      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
+          phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
 
     USE lmdz_phys_para, ONLY: is_master
@@ -592,6 +583,6 @@
     USE lmdz_phys_transfert_para, ONLY: gather
     USE phys_cal_mod, ONLY: year_len
-    USE netcdf, ONLY:nf90_clobber,nf90_close,nf90_noerr,nf90_strerror,nf90_put_att,nf90_def_var,&
-            nf90_def_dim,nf90_create,nf90_put_var,nf90_unlimited,nf90_global,nf90_64bit_offset,&
+    USE netcdf, ONLY: nf90_clobber, nf90_close, nf90_noerr, nf90_strerror, nf90_put_att, nf90_def_var, &
+            nf90_def_dim, nf90_create, nf90_put_var, nf90_unlimited, nf90_global, nf90_64bit_offset, &
             nf90_enddef
     USE lmdz_cppkeys_wrapper, ONLY: nf90_format
@@ -611,5 +602,5 @@
 
     REAL :: phy_glo(klon_glo, year_len) ! temporary variable, to store phy_***(:)
-      ! on the whole physics grid
+    ! on the whole physics grid
     INTEGER :: k
     INTEGER ierr
@@ -627,8 +618,8 @@
       PRINT *, 'writelim: Ecriture du fichier limit'
 
-      ierr = nf90_create('limit.nc', IOR(nf90_clobber,nf90_64bit_offset), nid)
+      ierr = nf90_create('limit.nc', IOR(nf90_clobber, nf90_64bit_offset), nid)
 
       ierr = nf90_put_att(nid, nf90_global, 'title', &
-        'Fichier conditions aux limites')
+              'Fichier conditions aux limites')
       !        ierr = nf90_def_dim (nid, "points_physiques", klon, ndim)
       ierr = nf90_def_dim(nid, 'points_physiques', klon_glo, ndim)
@@ -643,13 +634,13 @@
       ierr = nf90_def_var(nid, 'NAT', nf90_format, dims, id_nat)
       ierr = nf90_put_att(nid, id_nat, 'title', &
-        'Nature du sol (0,1,2,3)')
+              'Nature du sol (0,1,2,3)')
 
       ierr = nf90_def_var(nid, 'SST', nf90_format, dims, id_sst)
       ierr = nf90_put_att(nid, id_sst, 'title', &
-        'Temperature superficielle de la mer')
+              'Temperature superficielle de la mer')
 
       ierr = nf90_def_var(nid, 'BILS', nf90_format, dims, id_bils)
       ierr = nf90_put_att(nid, id_bils, 'title', &
-        'Reference flux de chaleur au sol')
+              'Reference flux de chaleur au sol')
 
       ierr = nf90_def_var(nid, 'ALB', nf90_format, dims, id_alb)
@@ -660,11 +651,11 @@
 
       ierr = nf90_def_var(nid, 'FTER', nf90_format, dims, id_fter)
-      ierr = nf90_put_att(nid, id_fter, 'title','Frac. Land')
+      ierr = nf90_put_att(nid, id_fter, 'title', 'Frac. Land')
       ierr = nf90_def_var(nid, 'FOCE', nf90_format, dims, id_foce)
-      ierr = nf90_put_att(nid, id_foce, 'title','Frac. Ocean')
+      ierr = nf90_put_att(nid, id_foce, 'title', 'Frac. Ocean')
       ierr = nf90_def_var(nid, 'FSIC', nf90_format, dims, id_fsic)
-      ierr = nf90_put_att(nid, id_fsic, 'title','Frac. Sea Ice')
+      ierr = nf90_put_att(nid, id_fsic, 'title', 'Frac. Sea Ice')
       ierr = nf90_def_var(nid, 'FLIC', nf90_format, dims, id_flic)
-      ierr = nf90_put_att(nid, id_flic, 'title','Frac. Land Ice')
+      ierr = nf90_put_att(nid, id_flic, 'title', 'Frac. Land Ice')
 
       ierr = nf90_enddef(nid)
@@ -780,5 +771,5 @@
   SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
     USE dimphy
-    USE phys_cal_mod , ONLY: year_len
+    USE phys_cal_mod, ONLY: year_len
     IMPLICIT NONE
 
@@ -787,12 +778,12 @@
     INTEGER imn, imx, amn, amx, kmn, kmx
     INTEGER p, pplus, nlat_max
-    PARAMETER (nlat_max=72)
+    PARAMETER (nlat_max = 72)
     REAL x_anom_sst(nlat_max)
-    CHARACTER (LEN=20) :: modname='profil_sst'
-    CHARACTER (LEN=80) :: abort_message
+    CHARACTER (LEN = 20) :: modname = 'profil_sst'
+    CHARACTER (LEN = 80) :: abort_message
 
     IF (klon/=nlon) THEN
-       abort_message='probleme de dimensions dans profil_sst'
-       CALL abort_physic(modname,abort_message,1)
+      abort_message = 'probleme de dimensions dans profil_sst'
+      CALL abort_physic(modname, abort_message, 1)
     ENDIF
     WRITE (*, *) ' profil_sst: type_profil=', type_profil
@@ -805,5 +796,5 @@
         ! Méthode 1 "Control" faible plateau à l'Equateur
         DO j = 1, klon
-          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**2)
+          phy_sst(j, i) = 273. + 27. * (1 - sin(1.5 * rlatd(j))**2)
           ! PI/3=1.047197551
           IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
@@ -815,5 +806,5 @@
         ! Méthode 2 "Flat" fort plateau à l'Equateur
         DO j = 1, klon
-          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**4)
+          phy_sst(j, i) = 273. + 27. * (1 - sin(1.5 * rlatd(j))**4)
           IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
             phy_sst(j, i) = 273.
@@ -821,11 +812,10 @@
         END DO
       END IF
-
 
       IF (type_profil==3) THEN
         ! Méthode 3 "Qobs" plateau réel à l'Equateur
         DO j = 1, klon
-          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
-            rlatd(j))**4)
+          phy_sst(j, i) = 273. + 0.5 * 27. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 * &
+                  rlatd(j))**4)
           IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
             phy_sst(j, i) = 273.
@@ -837,6 +827,6 @@
         ! Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur
         DO j = 1, klon
-          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
-            rlatd(j))**4)
+          phy_sst(j, i) = 273. + 0.5 * 29. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 * &
+                  rlatd(j))**4)
           IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
             phy_sst(j, i) = 273.
@@ -848,6 +838,6 @@
         ! Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur
         DO j = 1, klon
-          phy_sst(j, i) = 273. + 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
-            *rlatd(j))**4)
+          phy_sst(j, i) = 273. + 2. + 0.5 * 27. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 &
+                  * rlatd(j))**4)
           IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
             phy_sst(j, i) = 273. + 2.
@@ -863,5 +853,4 @@
         END DO
       END IF
-
 
       IF (type_profil==7) THEN
@@ -876,7 +865,7 @@
         ! Méthode 8 profil anomalies SST du modèle couplé AR4
         DO j = 1, klon
-          IF (rlatd(j)==rlatd(j-1)) THEN
+          IF (rlatd(j)==rlatd(j - 1)) THEN
             phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
-              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
+                    0.5 * 27. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 * rlatd(j))**4)
             IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
               phy_sst(j, i) = 273. + x_anom_sst(pplus)
@@ -886,5 +875,5 @@
             pplus = 73 - p
             phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
-              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
+                    0.5 * 27. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 * rlatd(j))**4)
             IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
               phy_sst(j, i) = 273. + x_anom_sst(pplus)
@@ -898,6 +887,6 @@
         ! Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur
         DO j = 1, klon
-          phy_sst(j, i) = 273. - 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
-            *rlatd(j))**4)
+          phy_sst(j, i) = 273. - 2. + 0.5 * 27. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 &
+                  * rlatd(j))**4)
           IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
             phy_sst(j, i) = 273. - 2.
@@ -905,11 +894,10 @@
         END DO
       END IF
-
 
       IF (type_profil==10) THEN
         ! Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur
         DO j = 1, klon
-          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
-            *rlatd(j))**4)
+          phy_sst(j, i) = 273. + 4. + 0.5 * 27. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 &
+                  * rlatd(j))**4)
           IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
             phy_sst(j, i) = 273. + 4.
@@ -921,6 +909,6 @@
         ! Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur
         DO j = 1, klon
-          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
-            rlatd(j))**4)
+          phy_sst(j, i) = 273. + 0.5 * 27. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 * &
+                  rlatd(j))**4)
           IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
             phy_sst(j, i) = 273.
@@ -932,6 +920,6 @@
         ! Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur
         DO j = 1, klon
-          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
-            *rlatd(j))**4)
+          phy_sst(j, i) = 273. + 4. + 0.5 * 27. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 &
+                  * rlatd(j))**4)
           IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
             phy_sst(j, i) = 273. + 4.
@@ -943,6 +931,6 @@
         ! Méthode 13 "Qmax" plateau réel à l'Equateur augmenté !
         DO j = 1, klon
-          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
-            rlatd(j))**4)
+          phy_sst(j, i) = 273. + 0.5 * 29. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 * &
+                  rlatd(j))**4)
           IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
             phy_sst(j, i) = 273.
@@ -954,6 +942,6 @@
         ! Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K !
         DO j = 1, klon
-          phy_sst(j, i) = 273. + 2. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
-            *rlatd(j))**4)
+          phy_sst(j, i) = 273. + 2. + 0.5 * 29. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 &
+                  * rlatd(j))**4)
           IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
             phy_sst(j, i) = 273.
@@ -963,27 +951,25 @@
 
       IF (type_profil==20) THEN
-      PRINT*,'Profile SST 20'
-!     Méthode 13 "Qmax2K" plateau réel é|  l'Equateur augmenté +2K
-
-      do j=1,klon
-        phy_sst(j,i)=248.+55.*(1-sin(rlatd(j))**2)
-      enddo
+        PRINT*, 'Profile SST 20'
+        !     Méthode 13 "Qmax2K" plateau réel é|  l'Equateur augmenté +2K
+
+        do j = 1, klon
+          phy_sst(j, i) = 248. + 55. * (1 - sin(rlatd(j))**2)
+        enddo
       endif
 
       IF (type_profil==21) THEN
-      PRINT*,'Profile SST 21'
-!     Méthode 13 "Qmax2K" plateau réel é|  l'Equateur augmenté +2K
-      do j=1,klon
-        phy_sst(j,i)=252.+55.*(1-sin(rlatd(j))**2)
-      enddo
+        PRINT*, 'Profile SST 21'
+        !     Méthode 13 "Qmax2K" plateau réel é|  l'Equateur augmenté +2K
+        do j = 1, klon
+          phy_sst(j, i) = 252. + 55. * (1 - sin(rlatd(j))**2)
+        enddo
       endif
 
-
-
     END DO
 
     ! IM beg : verif profil SST: phy_sst
-    amn = min(phy_sst(1,1), 1000.)
-    amx = max(phy_sst(1,1), -1000.)
+    amn = min(phy_sst(1, 1), 1000.)
+    amx = max(phy_sst(1, 1), -1000.)
     imn = 1
     kmn = 1
@@ -992,10 +978,10 @@
     DO k = 1, year_len
       DO i = 2, nlon
-        IF (phy_sst(i,k)<amn) THEN
+        IF (phy_sst(i, k)<amn) THEN
           amn = phy_sst(i, k)
           imn = i
           kmn = k
         END IF
-        IF (phy_sst(i,k)>amx) THEN
+        IF (phy_sst(i, k)>amx) THEN
           amx = phy_sst(i, k)
           imx = i
@@ -1009,5 +995,4 @@
     ! IM end : verif profil SST: phy_sst
 
-
   END SUBROUTINE profil_sst
 
Index: LMDZ6/branches/Amaury_dev/libf/phylmdiso/phyetat0_mod.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmdiso/phyetat0_mod.F90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmdiso/phyetat0_mod.F90	(revision 5144)
@@ -62,4 +62,5 @@
 #endif
 #endif
+USE lmdz_yomcst
 
   IMPLICIT NONE
@@ -69,5 +70,4 @@
   !======================================================================
   include "dimsoil.h"
-  include "YOMCST.h"
   !======================================================================
   CHARACTER*(*) fichnom
Index: LMDZ6/branches/Amaury_dev/libf/phylmdiso/physiq_mod.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmdiso/physiq_mod.F90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmdiso/physiq_mod.F90	(revision 5144)
@@ -424,6 +424,7 @@
     USE lmdz_conema3
     USE lmdz_dimpft, ONLY: nvm_lmdz
-    USE lmdz_YOETHF
+    USE lmdz_yoethf
     USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep
+    USE lmdz_yomcst
 
     IMPLICIT NONE
@@ -1268,6 +1269,4 @@
     REAL :: ro3i ! 0<=ro3i<=360 ; required time index in NetCDF file for
                  ! the ozone fields, old method.
-
-    include "YOMCST.h"
 
     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Index: LMDZ6/branches/Amaury_dev/libf/phylmdiso/reevap.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmdiso/reevap.F90	(revision 5143)
+++ LMDZ6/branches/Amaury_dev/libf/phylmdiso/reevap.F90	(revision 5144)
@@ -16,6 +16,7 @@
 #endif
 
-    USE lmdz_YOETHF
+    USE lmdz_yoethf
     USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep
+    USE lmdz_yomcst
 
     IMPLICIT NONE
@@ -32,7 +33,4 @@
 
     !--------Stochastic Boundary Layer Triggering: ALE_BL--------
-    !---Propri\'et\'es du thermiques au LCL
-    include "YOMCST.h"
-    !IM 100106 BEG : pouvoir sortir les ctes de la physique
 
 DO ixt=1,1+ntiso
