Index: /LMDZ5/branches/LMDZ_tree_FC/DefLists/CMIP6_ping_atmos.xml
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/DefLists/CMIP6_ping_atmos.xml	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/DefLists/CMIP6_ping_atmos.xml	(revision 2925)
@@ -16,5 +16,5 @@
    <field id="CMIP6_albsn"         field_ref="dummy_XY"         /> <!-- P1 (1.0) albsn : Albedo of the snow-covered surface, averaged over the grid cell. -->
    <field id="CMIP6_aod550volso4"  field_ref="dummy_XY"         /> <!-- P1 (1e-09) aod550volso4 : aerosol optical depth at 550 nm due to stratospheric volcanic aerosols  -->
-   <field id="CMIP6_areacella"     field_ref="dummy_XY"         /> <!-- P1 (m2) cell_area : For atmospheres with more than 1 mesh (e.g., staggered grids), report areas that apply to surface vertical fluxes of energy. -->
+   <field id="CMIP6_areacella"     field_ref="aire"         /> <!-- P1 (m2) cell_area : For atmospheres with more than 1 mesh (e.g., staggered grids), report areas that apply to surface vertical fluxes of energy. -->
    <field id="CMIP6_ares"          field_ref="dummy_XY"         /> <!-- P1 (s m-1) aerodynamic_resistance : Aerodynamic resistance -->
    <field id="CMIP6_cLand"         field_ref="dummy_XY"         /> <!-- P1 (kg m-2) cLand : as specified by C4MIP -->
@@ -22,5 +22,5 @@
    <field id="CMIP6_ccldncl"       field_ref="dummy_XY"         /> <!-- P1 (m-3) ccldncl : Droplets are liquid only.  Report concentration 'as seen from space' over convective liquid cloudy portion of grid cell.  This is the value from uppermost model layer with liquid cloud or, if available, it is better to sum over all liquid cloud tops, no matter where they occur, as long as they are seen from the top of the atmosphere. Weight by total liquid cloud top fraction of  (as seen from TOA) each time sample when computing monthly mean. -->
    <field id="CMIP6_cct"           field_ref="ptop"             /> <!-- P1 (Pa) air_pressure_at_convective_cloud_top : Where convective cloud is present in the grid cell, the instantaneous cloud top altitude should be that of the top of the highest level containing convective cloud. Missing data should be reported in the absence of convective cloud. The time mean should be calculated from these quantities averaging over occasions when convective cloud is present only, and should contain missing data for occasions when no convective cloud is present during the meaning period. -->
-   <field id="CMIP6_cfadDbze94"    field_ref="dummy_XYA"        /> <!-- P1 (1.0) histogram_of_equivalent_reflectivity_factor_over_height_above_reference_ellipsoid : CFAD (Cloud Frequency Altitude Diagrams) are frequency distributions of radar  reflectivity (or lidar scattering ratio) as a function of altitude. The variable cfadDbze94 is defined as the simulated relative frequency of occurrence of radar reflectivity in sampling volumes defined by altitude bins. The radar is observing at a frequency of 94GHz. -->
+   <field id="CMIP6_cfadDbze94"    field_ref="cfadDbze94"        /> <!-- P1 (1.0) histogram_of_equivalent_reflectivity_factor_over_height_above_reference_ellipsoid : CFAD (Cloud Frequency Altitude Diagrams) are frequency distributions of radar  reflectivity (or lidar scattering ratio) as a function of altitude. The variable cfadDbze94 is defined as the simulated relative frequency of occurrence of radar reflectivity in sampling volumes defined by altitude bins. The radar is observing at a frequency of 94GHz. -->
    <field id="CMIP6_cfadLidarsr532" field_ref="cfad_lidarsr532" /> <!-- P1 (1.0) histogram_of_backscattering_ratio_over_height_above_reference_ellipsoid : CFAD (Cloud Frequency Altitude Diagrams) are frequency distributions of radar  reflectivity (or lidar scattering ratio) as a function of altitude. The variable cfadLidarsr532 is defined as the simulated relative frequency of lidar scattering ratio in sampling volumes defined by altitude bins. The lidar is observing at a wavelength of 532nm. -->
    <field id="CMIP6_cfc113global"  field_ref="dummy_0d"         /> <!-- P1 (1e-12) mole_fraction_of_cfc113_in_air : unset -->
@@ -34,9 +34,9 @@
    <field id="CMIP6_cl"            field_ref="rneb"             /> <!-- P1 (%) cloud_area_fraction_in_atmosphere_layer : Percentage cloud cover, including both large-scale and convective cloud. -->
    <field id="CMIP6_clayFrac"      field_ref="dummy_XY"         /> <!-- P1 (1.0) clayFrac : Clay Fraction -->
-   <field id="CMIP6_clc"           field_ref="dummy_XYA"        /> <!-- P1 (%) convective_cloud_area_fraction_in_atmosphere_layer : Include only convective cloud. -->
+   <field id="CMIP6_clc"           field_ref="rnebcon"        /> <!-- P1 (%) convective_cloud_area_fraction_in_atmosphere_layer : Include only convective cloud. -->
    <field id="CMIP6_clcalipso"     field_ref="clcalipso"        /> <!-- P1 (%) cloud_area_fraction_in_atmosphere_layer : Percentage cloud cover at CALIPSO standard heights. -->
-   <field id="CMIP6_clcalipso2"    field_ref="dummy_XYA"        /> <!-- P1 (%) cloud_area_fraction_in_atmosphere_layer : Clouds detected by CALIPSO but below the detectability threshold of CloudSat -->
-   <field id="CMIP6_clcalipsoice"  field_ref="dummy_XYA"        /> <!-- P1 (%) ice_cloud_area_fraction_in_atmosphere_layer : CALIPSO ice cloud Fraction -->
-   <field id="CMIP6_clcalipsoliq"  field_ref="dummy_XYA"        /> <!-- P1 (%) clcalipsoliq : CALIPSO liquid cloud Fraction -->
+   <field id="CMIP6_clcalipso2"    field_ref="clcalipso2"        /> <!-- P1 (%) cloud_area_fraction_in_atmosphere_layer : Clouds detected by CALIPSO but below the detectability threshold of CloudSat -->
+   <field id="CMIP6_clcalipsoice"  field_ref="clcalipsoice"        /> <!-- P1 (%) ice_cloud_area_fraction_in_atmosphere_layer : CALIPSO ice cloud Fraction -->
+   <field id="CMIP6_clcalipsoliq"  field_ref="clcalipsoice"        /> <!-- P1 (%) clcalipsoliq : CALIPSO liquid cloud Fraction -->
    <field id="CMIP6_cldicemxrat27" field_ref="dummy_XYA"        /> <!-- P2 (1.0) cloud_ice_mixing_ratio : Cloud ice mixing ratio -->
    <field id="CMIP6_cldnci"        field_ref="dummy_XY"         /> <!-- P1 (m-3) number_concentration_of_ice_crystals_in_air_at_ice_cloud_top : Concentration 'as seen from space' over ice-cloud portion of grid cell.  This is the value from uppermost model layer with ice cloud or, if available, it is the sum over all ice cloud tops, no matter where they occur, as long as they are seen from the top of the atmosphere. Weight by total ice cloud top fraction (as seen from TOA) of each time sample when computing monthly mean. -->
@@ -47,20 +47,20 @@
    <field id="CMIP6_cli"           field_ref="iwcon"            /> <!-- P1 (kg kg-1) mass_fraction_of_cloud_ice_in_air : Includes both large-scale and convective cloud. This is calculated as the mass of cloud ice in the grid cell divided by the mass of air (including the water in all phases) in the grid cell. It includes precipitating hydrometeors ONLY if the precipitating hydrometeors affect the calculation of radiative transfer in model. -->
    <field id="CMIP6_clic"          field_ref="dummy_XYA"        /> <!-- P2 (1.0) mass_fraction_of_convective_cloud_ice_in_air : Calculated as the mass of convective cloud ice  in the grid cell divided by the mass of air (including the water in all phases) in the grid cell.  This includes precipitating hydrometeors ONLY if the precipitating hydrometeors affect the calculation of radiative transfer in model. -->
-   <field id="CMIP6_climodis"      field_ref="dummy_XY"         /> <!-- P1 (%) ice_cloud_area_fraction : MODIS Ice Cloud Fraction -->
+   <field id="CMIP6_climodis"      field_ref="climodis"         /> <!-- P1 (%) ice_cloud_area_fraction : MODIS Ice Cloud Fraction -->
    <field id="CMIP6_clis"          field_ref="dummy_XYA"        /> <!-- P2 (1.0) mass_fraction_of_stratiform_cloud_ice_in_air : Calculated as the mass of stratiform cloud ice  in the grid cell divided by the mass of air (including the water in all phases) in the grid cell.  This includes precipitating hydrometeors ONLY if the precipitating hydrometeors affect the calculation of radiative transfer in model. -->
-   <field id="CMIP6_clisccp"       field_ref="dummy_XYA"        /> <!-- P1 (%) isccp_cloud_area_fraction : Percentage cloud cover in optical depth categories. -->
+   <field id="CMIP6_clisccp"       field_ref="clisccp2"        /> <!-- P1 (%) isccp_cloud_area_fraction : Percentage cloud cover in optical depth categories. -->
    <field id="CMIP6_clivi"         field_ref="iwp"              /> <!-- P1 (kg m-2) atmosphere_mass_content_of_cloud_ice : calculate mass of ice water in the column divided by the area of the column (not just the area of the cloudy portion of the column). This includes precipitating frozen hydrometeors ONLY if the precipitating hydrometeors affect the calculation of radiative transfer in model. -->
    <field id="CMIP6_clivic"        field_ref="dummy_XY"         /> <!-- P1 (kg m-2 ) clivic : calculate mass of convective ice water in the column divided by the area of the column (not just the area of the cloudy portion of the column). This includes precipitating frozen hydrometeors ONLY if the precipitating hydrometeors affect the calculation of radiative transfer in model.   -->
    <field id="CMIP6_cllcalipso"    field_ref="cllcalipso"       /> <!-- P1 (%) cloud_area_fraction_in_atmosphere_layer : Percentage cloud cover in layer centred on 840hPa -->
    <field id="CMIP6_clmcalipso"    field_ref="clmcalipso"       /> <!-- P1 (%) cloud_area_fraction_in_atmosphere_layer : Percentage cloud cover in layer centred on 560hPa -->
-   <field id="CMIP6_clmisr"        field_ref="dummy_XYA"        /> <!-- P1 (%) cloud_area_fraction_in_atmosphere_layer : Cloud percentage in spectral bands and layers as observed by the Multi-angle Imaging SpectroRadiometer (MISR) instrument. -->
+   <field id="CMIP6_clmisr"        field_ref="clMISR"        /> <!-- P1 (%) cloud_area_fraction_in_atmosphere_layer : Cloud percentage in spectral bands and layers as observed by the Multi-angle Imaging SpectroRadiometer (MISR) instrument. -->
    <field id="CMIP6_cls"           field_ref="dummy_XYA"        /> <!-- P1 (%) stratiform_cloud_area_fraction_in_atmosphere_layer : unset -->
    <field id="CMIP6_clt"           field_ref="cldt"             /> <!-- P1 (1.0) cloud_area_fraction : Total cloud area fraction for the whole atmospheric column, as seen from the surface or the top of the atmosphere. Includes both large-scale and convective cloud. -->
    <field id="CMIP6_cltcalipso"    field_ref="cltcalipso"       /> <!-- P1 (%) cloud_area_fraction : unset -->
    <field id="CMIP6_cltisccp"      field_ref="tclisccp"         /> <!-- P1 (%) cloud_area_fraction : Percentage total cloud cover, simulating ISCCP observations. -->
-   <field id="CMIP6_cltmodis"      field_ref="dummy_XY"         /> <!-- P1 (%) cloud_area_fraction : MODIS Total Cloud Fraction -->
+   <field id="CMIP6_cltmodis"      field_ref="cltmodis"         /> <!-- P1 (%) cloud_area_fraction : MODIS Total Cloud Fraction -->
    <field id="CMIP6_clw"           field_ref="lwcon"            /> <!-- P1 (kg kg-1) mass_fraction_of_cloud_liquid_water_in_air : Includes both large-scale and convective cloud. Calculate as the mass of cloud liquid water in the grid cell divided by the mass of air (including the water in all phases) in the grid cells. Precipitating hydrometeors are included ONLY if the precipitating hydrometeors affect the calculation of radiative transfer in model. -->
-   <field id="CMIP6_clwc"          field_ref="dummy_XYA"        /> <!-- P2 (1.0) mass_fraction_of_convective_cloud_liquid_water_in_air : Calculated as the mass of convective cloud liquid water in the grid cell divided by the mass of air (including the water in all phases) in the grid cell.  This includes precipitating hydrometeors ONLY if the precipitating hydrometeors affect the calculation of radiative transfer in model. -->
-   <field id="CMIP6_clwmodis"      field_ref="dummy_XY"         /> <!-- P1 (%) clwmodis : MODIS Liquid Cloud Fraction -->
+   <field id="CMIP6_clwc"          field_ref="lcc3dcon"        /> <!-- P2 (1.0) mass_fraction_of_convective_cloud_liquid_water_in_air : Calculated as the mass of convective cloud liquid water in the grid cell divided by the mass of air (including the water in all phases) in the grid cell.  This includes precipitating hydrometeors ONLY if the precipitating hydrometeors affect the calculation of radiative transfer in model. -->
+   <field id="CMIP6_clwmodis"      field_ref="clwmodis"         /> <!-- P1 (%) clwmodis : MODIS Liquid Cloud Fraction -->
    <field id="CMIP6_clws"          field_ref="lcc3dstra"        /> <!-- P2 (1.0) mass_fraction_of_stratiform_cloud_liquid_water_in_air : Calculated as the mass of stratiform cloud liquid water in the grid cell divided by the mass of air (including the water in all phases) in the grid cell.  This includes precipitating hydrometeors ONLY if the precipitating hydrometeors affect the calculation of radiative transfer in model. -->
    <field id="CMIP6_clwvi"         field_ref="lwp"              /> <!-- P1 (kg m-2) atmosphere_cloud_condensed_water_content : Mass of condensed (liquid + ice) water in the column divided by the area of the column (not just the area of the cloudy portion of the column). Includes precipitating hydrometeors ONLY if the precipitating hydrometeors affect the calculation of radiative transfer in model. -->
@@ -82,5 +82,5 @@
    <field id="CMIP6_dgw"           field_ref="dummy_XY"         /> <!-- P1 (kg m-2) dgw : Change in Groundwater -->
    <field id="CMIP6_diabdrag"      field_ref="dummy_XYA"        /> <!-- P1 (m s-2) tendency_of_eastward_wind_due_to_numerical_artefacts : Other sub-grid scale/numerical zonal drag excluding that already provided for the parameterized orographic and non-ororgraphic gravity waves. This would be used to calculate the total 'diabatic drag'. Contributions to this additional drag such Rayleigh friction and diffusion that can be calculated from the monthly mean wind fields should not be included, but details (e.g. coefficients) of the friction and/or diffusion used in the model should be provided separately. -->
-   <field id="CMIP6_dmc"           field_ref="dummy_XYA"        /> <!-- P2 (kg m-2 s-1) atmosphere_net_upward_deep_convective_mass_flux : The net mass flux  represents the difference between the updraft and downdraft components.   This is calculated as the convective mass flux divided by the area of the whole grid cell (not just the area of the cloud). -->
+   <field id="CMIP6_dmc"           field_ref="upwd"        /> <!-- P2 (kg m-2 s-1) atmosphere_net_upward_deep_convective_mass_flux : The net mass flux  represents the difference between the updraft and downdraft components.   This is calculated as the convective mass flux divided by the area of the whole grid cell (not just the area of the cloud). -->
    <field id="CMIP6_dmlt"          field_ref="dummy_XY"         /> <!-- P1 (m) dmlt : Depth from surface to the zero degree isotherm. Above this isotherm T > 0o, and below this line T < 0o. -->
    <field id="CMIP6_drivw"         field_ref="dummy_XY"         /> <!-- P1 (kg m-2) drivw : Change in River Storage -->
@@ -135,6 +135,6 @@
    <field id="CMIP6_jo2"           field_ref="dummy_lat-P"      /> <!-- P1 (s-1) jo2 : rate of o2 -> o1d+o -->
    <field id="CMIP6_jo3"           field_ref="dummy_lat-P"      /> <!-- P1 (s-1) jo3 : sum of rates o3 -> o1d+o2 and o3 -> o+o2 -->
-   <field id="CMIP6_jpdftaureicemodis" field_ref="dummy_XYA"        /> <!-- P1 (%) cloud_area_fraction_in_atmosphere_layer : MODIS Optical Thickness-Particle Size joint  distribution, ice -->
-   <field id="CMIP6_jpdftaureliqmodis" field_ref="dummy_XYA"        /> <!-- P1 (%) cloud_area_fraction_in_atmosphere_layer : MODIS Optical Thickness-Particle Size joint  distribution, liquid -->
+   <field id="CMIP6_jpdftaureicemodis" field_ref="crimodis"        /> <!-- P1 (%) cloud_area_fraction_in_atmosphere_layer : MODIS Optical Thickness-Particle Size joint  distribution, ice -->
+   <field id="CMIP6_jpdftaureliqmodis" field_ref="crlmodis"        /> <!-- P1 (%) cloud_area_fraction_in_atmosphere_layer : MODIS Optical Thickness-Particle Size joint  distribution, liquid -->
    <field id="CMIP6_ksat"          field_ref="dummy_XY"         /> <!-- P1 (1e-6 m s-1) ksat : Saturated Hydraulic Conductivity -->
    <field id="CMIP6_latitude"      field_ref="dummy_COSPcurtain"/> <!-- P1 (degrees_north) latitude : latitude -->
@@ -176,5 +176,5 @@
    <field id="CMIP6_oxloss"        field_ref="dummy_lat-P"      /> <!-- P1 (mol m-3 s-1) oxloss : total chemical loss rate for o+o1d+o3  -->
    <field id="CMIP6_oxprod"        field_ref="dummy_lat-P"      /> <!-- P1 (mol m-3 s-1) oxprod : total production rate of o+o1d+o3 including o2 photolysis and all o3 producing reactions -->
-   <field id="CMIP6_parasolRefl"   field_ref="dummy_XY"         /> <!-- P1 (1.0) toa_bidirectional_reflectance : Simulated reflectance from PARASOL as seen at the top of the atmosphere for 5 solar zenith angles. Valid only over ocean and for one viewing direction (viewing zenith angle of 30 degrees and relative azimuth angle 320 degrees). -->
+   <field id="CMIP6_parasolRefl"   field_ref="parasol_refl"         /> <!-- P1 (1.0) toa_bidirectional_reflectance : Simulated reflectance from PARASOL as seen at the top of the atmosphere for 5 solar zenith angles. Valid only over ocean and for one viewing direction (viewing zenith angle of 30 degrees and relative azimuth angle 320 degrees). -->
    <field id="CMIP6_parasolRefl_sea" field_ref="dummy_XY"         /> <!-- P1 (1.0) toa_bidirectional_reflectance : Simulated reflectance from PARASOL as seen at the top of the atmosphere for 5 solar zenith angles. Valid only over ocean and for one viewing direction (viewing zenith angle of 30 degrees and relative azimuth angle 320 degrees). -->
    <field id="CMIP6_pctisccp"      field_ref="ctpisccp"         /> <!-- P1 (Pa) air_pressure_at_cloud_top : ISCCP Mean Cloud Top Pressure. Time-means are weighted by the ISCCP Total Cloud Fraction {:cltisccp} - see  http://cfmip.metoffice.com/COSP.html -->
@@ -213,27 +213,27 @@
    <field id="CMIP6_rivo"          field_ref="dummy_XY"         /> <!-- P1 (m3 s-1) rivo : Outflow of River Water from Cell -->
    <field id="CMIP6_rld"           field_ref="rld"              /> <!-- P1 (W m-2) downwelling_longwave_flux_in_air : Downwelling Longwave Radiation (includes the fluxes at the surface and TOA) -->
-   <field id="CMIP6_rld4co2"       field_ref="dummy_XYA"        /> <!-- P1 (W m-2) downwelling_longwave_flux_in_air : Downwelling longwave radiation calculated using carbon dioxide concentrations increased fourfold (includes the fluxes at the surface and TOA) -->
+   <field id="CMIP6_rld4co2"       field_ref="rld4co2"        /> <!-- P1 (W m-2) downwelling_longwave_flux_in_air : Downwelling longwave radiation calculated using carbon dioxide concentrations increased fourfold (includes the fluxes at the surface and TOA) -->
    <field id="CMIP6_rldcs"         field_ref="rldcs"            /> <!-- P1 (W m-2) downwelling_longwave_flux_in_air_assuming_clear_sky : Downwelling clear-sky longwave radiation (includes the fluxes at the surface and TOA) -->
-   <field id="CMIP6_rldcs4co2"     field_ref="dummy_XYA"        /> <!-- P1 (W m-2) downwelling_longwave_flux_in_air_assuming_clear_sky : Downwelling clear-sky longwave radiation calculated using carbon dioxide concentrations increased fourfold (includes the fluxes at the surface and TOA) -->
+   <field id="CMIP6_rldcs4co2"     field_ref="rldcs4co2"        /> <!-- P1 (W m-2) downwelling_longwave_flux_in_air_assuming_clear_sky : Downwelling clear-sky longwave radiation calculated using carbon dioxide concentrations increased fourfold (includes the fluxes at the surface and TOA) -->
    <field id="CMIP6_rlds"          field_ref="LWdnSFC"          /> <!-- P1 (W m-2) surface_downwelling_longwave_flux_in_air : Surface Downwelling Longwave Radiation -->
    <field id="CMIP6_rlds_isf"      field_ref="dummy_??"         /> <!-- P1 (W m-2) surface_downwelling_longwave_flux_in_air : Surface Downwelling Longwave Radiation -->
    <field id="CMIP6_rldscs"        field_ref="LWdnSFCclr"       /> <!-- P1 (W m-2) surface_downwelling_longwave_flux_in_air_assuming_clear_sky : Surface downwelling clear-sky longwave radiation -->
-   <field id="CMIP6_rls"           field_ref="dummy_XY"         /> <!-- P1 (W m-2) surface_net_downward_longwave_flux : Net longwave surface radiation -->
+   <field id="CMIP6_rls"           field_ref="soll"         /> <!-- P1 (W m-2) surface_net_downward_longwave_flux : Net longwave surface radiation -->
    <field id="CMIP6_rls_land"      field_ref="dummy_XY"         /> <!-- P1 (W m-2) surface_net_downward_longwave_flux : Net longwave surface radiation -->
-   <field id="CMIP6_rlu"           field_ref="rli"              /> <!-- P1 (W m-2) upwelling_longwave_flux_in_air : Upwelling longwave radiation (includes the fluxes at the surface and TOA) -->
-   <field id="CMIP6_rlu4co2"       field_ref="dummy_XYA"        /> <!-- P1 (W m-2) upwelling_longwave_flux_in_air : Upwelling longwave radiation calculated using carbon dioxide concentrations increased fourfold (includes the fluxes at the surface and TOA) -->
+   <field id="CMIP6_rlu"           field_ref="rlu"              /> <!-- P1 (W m-2) upwelling_longwave_flux_in_air : Upwelling longwave radiation (includes the fluxes at the surface and TOA) -->
+   <field id="CMIP6_rlu4co2"       field_ref="rlu4co2"        /> <!-- P1 (W m-2) upwelling_longwave_flux_in_air : Upwelling longwave radiation calculated using carbon dioxide concentrations increased fourfold (includes the fluxes at the surface and TOA) -->
    <field id="CMIP6_rlucs"         field_ref="rlucs"            /> <!-- P1 (W m-2) upwelling_longwave_flux_in_air_assuming_clear_sky : Upwelling clear-sky longwave radiation  (includes the fluxes at the surface and TOA) -->
-   <field id="CMIP6_rlucs4co2"     field_ref="dummy_XYA"        /> <!-- P1 (W m-2) upwelling_longwave_flux_in_air_assuming_clear_sky : Upwelling clear-sky longwave radiation calculated using carbon dioxide concentrations increased fourfold (includes the fluxes at the surface and TOA) -->
+   <field id="CMIP6_rlucs4co2"     field_ref="rlucs4co2"        /> <!-- P1 (W m-2) upwelling_longwave_flux_in_air_assuming_clear_sky : Upwelling clear-sky longwave radiation calculated using carbon dioxide concentrations increased fourfold (includes the fluxes at the surface and TOA) -->
    <field id="CMIP6_rlus"          field_ref="LWupSFC"          /> <!-- P1 (W m-2) surface_upwelling_longwave_flux_in_air : Surface Upwelling Longwave Radiation -->
    <field id="CMIP6_rlus_isf"      field_ref="dummy_??"         /> <!-- P1 (W m-2) surface_upwelling_longwave_flux_in_air : Surface Upwelling Longwave Radiation -->
    <field id="CMIP6_rlut"          field_ref="topl"             /> <!-- P1 (W m-2) toa_outgoing_longwave_flux : at the top of the atmosphere (to be compared with satellite measurements) -->
-   <field id="CMIP6_rlut4co2"      field_ref="dummy_XY"         /> <!-- P1 (W m-2) toa_outgoing_longwave_flux : Top-of-atmosphere outgoing longwave radiation calculated using carbon dioxide concentrations increased fourfold -->
+   <field id="CMIP6_rlut4co2"      field_ref="rlut4co2"         /> <!-- P1 (W m-2) toa_outgoing_longwave_flux : Top-of-atmosphere outgoing longwave radiation calculated using carbon dioxide concentrations increased fourfold -->
    <field id="CMIP6_rlutcs"        field_ref="topl0"            /> <!-- P1 (W m-2) toa_outgoing_longwave_flux_assuming_clear_sky : TOA Outgoing Clear-sky Longwave Radiation -->
-   <field id="CMIP6_rlutcs4co2"    field_ref="dummy_XY"         /> <!-- P1 (W m-2) toa_outgoing_longwave_flux_assuming_clear_sky : Top-of-atmosphere outgoing clear-sky longwave radiation calculated using carbon dioxide concentrations increased fourfold -->
+   <field id="CMIP6_rlutcs4co2"    field_ref="rlutcs4co2"         /> <!-- P1 (W m-2) toa_outgoing_longwave_flux_assuming_clear_sky : Top-of-atmosphere outgoing clear-sky longwave radiation calculated using carbon dioxide concentrations increased fourfold -->
    <field id="CMIP6_rootdsl"       field_ref="dummy_XY"         /> <!-- P1 (kg m-3) rootdsl : Root Distribution -->
    <field id="CMIP6_rsd"           field_ref="rsd"              /> <!-- P1 (W m-2) downwelling_shortwave_flux_in_air : Downwelling shortwave radiation (includes the fluxes at the surface and top-of-atmosphere) -->
-   <field id="CMIP6_rsd4co2"       field_ref="dummy_XYA"        /> <!-- P1 (W m-2) downwelling_shortwave_flux_in_air : Downwelling shortwave radiation calculated using carbon dioxide concentrations increased fourfold -->
+   <field id="CMIP6_rsd4co2"       field_ref="rsd4co2"        /> <!-- P1 (W m-2) downwelling_shortwave_flux_in_air : Downwelling shortwave radiation calculated using carbon dioxide concentrations increased fourfold -->
    <field id="CMIP6_rsdcs"         field_ref="rsdcs"            /> <!-- P1 (W m-2) downwelling_shortwave_flux_in_air_assuming_clear_sky : Downwelling clear-sky shortwave radiation (includes the fluxes at the surface and top-of-atmosphere) -->
-   <field id="CMIP6_rsdcs4co2"     field_ref="dummy_XYA"        /> <!-- P1 (W m-2) downwelling_shortwave_flux_in_air_assuming_clear_sky : Downwelling clear-sky shortwave radiation calculated using carbon dioxide concentrations increased fourfold -->
+   <field id="CMIP6_rsdcs4co2"     field_ref="rsdcs4co2"        /> <!-- P1 (W m-2) downwelling_shortwave_flux_in_air_assuming_clear_sky : Downwelling clear-sky shortwave radiation calculated using carbon dioxide concentrations increased fourfold -->
    <field id="CMIP6_rsds"          field_ref="SWdnSFC"          /> <!-- P1 (W m-2) surface_downwelling_shortwave_flux_in_air : surface solar irradiance for UV calculations -->
    <field id="CMIP6_rsds_isf"      field_ref="dummy_??"         /> <!-- P1 (W m-2) surface_downwelling_shortwave_flux_in_air : surface solar irradiance for UV calculations -->
@@ -242,17 +242,17 @@
    <field id="CMIP6_rsdsdiff"      field_ref="dummy_XY"         /> <!-- P1 (W m-2) surface_diffuse_downwelling_shortwave_flux_in_air : unset -->
    <field id="CMIP6_rsdt"          field_ref="SWdnTOA"          /> <!-- P1 (W m-2) toa_incoming_shortwave_flux : Shortwave radiation incident at the top of the atmosphere -->
-   <field id="CMIP6_rss"           field_ref="dummy_XY"         /> <!-- P1 (W m-2) surface_net_downward_shortwave_flux : Net downward shortwave radiation at the surface -->
+   <field id="CMIP6_rss"           field_ref="sols"         /> <!-- P1 (W m-2) surface_net_downward_shortwave_flux : Net downward shortwave radiation at the surface -->
    <field id="CMIP6_rss_land"      field_ref="dummy_XY"         /> <!-- P1 (W m-2) surface_net_downward_shortwave_flux : Net downward shortwave radiation at the surface -->
    <field id="CMIP6_rsu"           field_ref="rsu"              /> <!-- P1 (W m-2) upwelling_shortwave_flux_in_air : Upwelling shortwave radiation  (includes also the fluxes at the surface and top of atmosphere) -->
-   <field id="CMIP6_rsu4co2"       field_ref="dummy_XYA"        /> <!-- P1 (W m-2) upwelling_shortwave_flux_in_air : Upwelling Shortwave Radiation calculated using carbon dioxide concentrations increased fourfold -->
+   <field id="CMIP6_rsu4co2"       field_ref="rsu4co2"        /> <!-- P1 (W m-2) upwelling_shortwave_flux_in_air : Upwelling Shortwave Radiation calculated using carbon dioxide concentrations increased fourfold -->
    <field id="CMIP6_rsucs"         field_ref="rsucs"            /> <!-- P1 (W m-2) upwelling_shortwave_flux_in_air_assuming_clear_sky : Upwelling clear-sky shortwave radiation  (includes the fluxes at the surface and TOA) -->
-   <field id="CMIP6_rsucs4co2"     field_ref="dummy_XYA"        /> <!-- P1 (W m-2) upwelling_shortwave_flux_in_air_assuming_clear_sky : Upwelling clear-sky shortwave radiation calculated using carbon dioxide concentrations increased fourfold -->
+   <field id="CMIP6_rsucs4co2"     field_ref="rsucs4co2"        /> <!-- P1 (W m-2) upwelling_shortwave_flux_in_air_assuming_clear_sky : Upwelling clear-sky shortwave radiation calculated using carbon dioxide concentrations increased fourfold -->
    <field id="CMIP6_rsus"          field_ref="SWupSFC"          /> <!-- P1 (W m-2) surface_upwelling_shortwave_flux_in_air : Surface Upwelling Shortwave Radiation -->
    <field id="CMIP6_rsus_isf"      field_ref="dummy_??"         /> <!-- P1 (W m-2) surface_upwelling_shortwave_flux_in_air : Surface Upwelling Shortwave Radiation -->
    <field id="CMIP6_rsuscs"        field_ref="SWupSFCclr"       /> <!-- P1 (W m-2) surface_upwelling_shortwave_flux_in_air_assuming_clear_sky : Surface upwelling clear-sky shortwave radiation -->
    <field id="CMIP6_rsut"          field_ref="SWupTOA"          /> <!-- P1 (W m-2) toa_outgoing_shortwave_flux : at the top of the atmosphere -->
-   <field id="CMIP6_rsut4co2"      field_ref="dummy_XY"         /> <!-- P1 (W m-2) toa_outgoing_shortwave_flux : TOA Outgoing Shortwave Radiation calculated using carbon dioxide concentrations increased fourfold -->
+   <field id="CMIP6_rsut4co2"      field_ref="rsut4co2"         /> <!-- P1 (W m-2) toa_outgoing_shortwave_flux : TOA Outgoing Shortwave Radiation calculated using carbon dioxide concentrations increased fourfold -->
    <field id="CMIP6_rsutcs"        field_ref="SWupTOAclr"       /> <!-- P1 (W m-2) toa_outgoing_shortwave_flux_assuming_clear_sky : Calculated in the absence of clouds. -->
-   <field id="CMIP6_rsutcs4co2"    field_ref="dummy_XY"         /> <!-- P1 (W m-2) toa_outgoing_shortwave_flux_assuming_clear_sky : TOA Outgoing Clear-Sky Shortwave Radiation calculated using carbon dioxide concentrations increased fourfold -->
+   <field id="CMIP6_rsutcs4co2"    field_ref="rsutcs4co2"         /> <!-- P1 (W m-2) toa_outgoing_shortwave_flux_assuming_clear_sky : TOA Outgoing Clear-Sky Shortwave Radiation calculated using carbon dioxide concentrations increased fourfold -->
    <field id="CMIP6_rsutna"        field_ref="dummy_XY"         /> <!-- P1 (W m-2) rsutna : Based on Ghan (2013, ACP) -->
    <field id="CMIP6_rsutnacs"      field_ref="dummy_XY"         /> <!-- P1 (W m-2) rsutnacs : Based on Ghan (2013, ACP) -->
@@ -263,11 +263,11 @@
    <field id="CMIP6_sci"           field_ref="ftime_th"         /> <!-- P1 (1.0) shallow_convection_time_fraction : Fraction of time that shallow convection occurs in the grid cell. -->
    <field id="CMIP6_scldncl"       field_ref="dummy_XY"         /> <!-- P1 (m-3) scldncl : Droplets are liquid only.  Report concentration "as seen from space" over stratiform liquid cloudy portion of grid cell.  This is the value from uppermost model layer with liquid cloud or, if available, it is better to sum over all liquid cloud tops, no matter where they occur, as long as they are seen from the top of the atmosphere. Weight by total liquid cloud top fraction of  (as seen from TOA) each time sample when computing monthly mean. -->
-   <field id="CMIP6_sconcdust"     field_ref="dummy_XY"         /> <!-- P1 (kg m-3) mass_concentration_of_dust_dry_aerosol_in_air : mass concentration of dust dry aerosol in air in model lowest layer -->
-   <field id="CMIP6_sconcso4"      field_ref="dummy_XY"         /> <!-- P1 (kg m-3) mass_concentration_of_sulfate_dry_aerosol_in_air : mass concentration of sulfate dry aerosol in air in model lowest layer. -->
-   <field id="CMIP6_sconcss"       field_ref="dummy_XY"         /> <!-- P1 (kg m-3) mass_concentration_of_seasalt_dry_aerosol_in_air : mass concentration of seasalt dry aerosol in air in model lowest layer -->
+   <field id="CMIP6_sconcdust"     field_ref="sconcdust"         /> <!-- P1 (kg m-3) mass_concentration_of_dust_dry_aerosol_in_air : mass concentration of dust dry aerosol in air in model lowest layer -->
+   <field id="CMIP6_sconcso4"      field_ref="sconcso4"         /> <!-- P1 (kg m-3) mass_concentration_of_sulfate_dry_aerosol_in_air : mass concentration of sulfate dry aerosol in air in model lowest layer. -->
+   <field id="CMIP6_sconcss"       field_ref="sconcss"         /> <!-- P1 (kg m-3) mass_concentration_of_seasalt_dry_aerosol_in_air : mass concentration of seasalt dry aerosol in air in model lowest layer -->
    <field id="CMIP6_sfcWind"       field_ref="wind10m"          /> <!-- P1 (m s-1) wind_speed : near-surface (usually, 10 meters) wind speed. -->
-   <field id="CMIP6_sfcWindmax"    field_ref="dummy_XY"         /> <!-- P1 (m s-1) wind_speed : Daily maximum near-surface (usually, 10 meters) wind speed. -->
-   <field id="CMIP6_sftlf"         field_ref="dummy_XY"         /> <!-- P1 (1) land_area_fraction : Please express "X_area_fraction" as the fraction of horizontal area occupied by X. -->
-   <field id="CMIP6_sic"           field_ref="dummy_XY"         /> <!-- P1 (1.0) sea_ice_area_fraction : fraction of grid cell covered by sea ice. -->
+   <field id="CMIP6_sfcWindmax"    field_ref="wind10max"         /> <!-- P1 (m s-1) wind_speed : Daily maximum near-surface (usually, 10 meters) wind speed. -->
+   <field id="CMIP6_sftlf"         field_ref="fract_ter"         /> <!-- P1 (1) land_area_fraction : Please express "X_area_fraction" as the fraction of horizontal area occupied by X. -->
+   <field id="CMIP6_sic"           field_ref="fract_sic"         /> <!-- P1 (1.0) sea_ice_area_fraction : fraction of grid cell covered by sea ice. -->
    <field id="CMIP6_siltFrac"      field_ref="dummy_XY"         /> <!-- P1 (1.0) siltFrac : Silt Fraction -->
    <field id="CMIP6_slbnosn"       field_ref="dummy_XY"         /> <!-- P1 (kg m-2 s-1) slbnosn : Sublimation of the snow free area -->
@@ -276,5 +276,5 @@
    <field id="CMIP6_sltnorth"      field_ref="dummy_basin_zonal_mean"/> <!-- P2 (kg s-1) northward_ocean_salt_transport : unset -->
    <field id="CMIP6_sltnortha"     field_ref="dummy_basin_zonal_mean"/> <!-- P1 (kg s-1) northward_ocean_salt_transport : unset -->
-   <field id="CMIP6_smc"           field_ref="dummy_XYA"        /> <!-- P2 (kg m-2 s-1) atmosphere_net_upward_shallow_convective_mass_flux : The net mass flux represents the difference between the updraft and downdraft components.  For models with a distinct shallow convection scheme, this is calculated as convective mass flux divided by the area of the whole grid cell (not just the area of the cloud). -->
+   <field id="CMIP6_smc"           field_ref="f_th"        /> <!-- P2 (kg m-2 s-1) atmosphere_net_upward_shallow_convective_mass_flux : The net mass flux represents the difference between the updraft and downdraft components.  For models with a distinct shallow convection scheme, this is calculated as convective mass flux divided by the area of the whole grid cell (not just the area of the cloud). -->
    <field id="CMIP6_snmsl"         field_ref="dummy_XY"         /> <!-- P1 (kg m-2 s-1) snmsl : Water flowing out of snowpack -->
    <field id="CMIP6_snowmxrat27"   field_ref="dummy_XYA"        /> <!-- P2 (1.0) mass_fraction_of_snow_in_air : Snow mixing ratio -->
@@ -288,8 +288,8 @@
    <field id="CMIP6_ta"            field_ref="ta"               /> <!-- P1 (K) air_temperature : Air Temperature -->
    <field id="CMIP6_ta27"          field_ref="dummy_XYA"        /> <!-- P3 (K) air_temperature : Air Temperature -->
-   <field id="CMIP6_ta500"         field_ref="dummy_XY"         /> <!-- P1 (K) air_temperature : Temperature on the 500 hPa surface -->
+   <field id="CMIP6_ta500"         field_ref="t500"             /> <!-- P1 (K) air_temperature : Temperature on the 500 hPa surface -->
    <field id="CMIP6_ta700"         field_ref="t700"             /> <!-- P1 (K) air_temperature : Air temperature at 700hPa -->
    <field id="CMIP6_ta7h"          field_ref="dummy_XYA"        /> <!-- P2 (K) air_temperature : Air Temperature -->
-   <field id="CMIP6_ta850"         field_ref="dummy_XY"         /> <!-- P1 (K) air_temperature : Air temperature at 850hPa -->
+   <field id="CMIP6_ta850"         field_ref="t850"             /> <!-- P1 (K) air_temperature : Air temperature at 850hPa -->
    <field id="CMIP6_tas"           field_ref="t2m"              /> <!-- P1 (K) air_temperature : near-surface (usually, 2 meter) air temperature -->
    <field id="CMIP6_tas_isf"       field_ref="dummy_XY"         /> <!-- P1 (K) air_temperature : near-surface (usually, 2 meter) air temperature -->
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90	(revision 2925)
@@ -43,5 +43,6 @@
     sig1, ftsol, clwcon, fm_therm, wake_Cstar,  pctsrf,  entr_therm,radpas, f0,&
     zmax0,fevap, rnebcon,falb_dir, wake_fip,    agesno,  detr_therm, pbl_tke,  &
-    phys_state_var_init
+    phys_state_var_init, ql_ancien, qs_ancien, prlw_ancien, prsw_ancien, &
+    prw_ancien
   USE comconst_mod, ONLY: pi, dtvr
 
@@ -201,4 +202,9 @@
   t_ancien   = 273.15
   q_ancien   = 0.
+  ql_ancien = 0.
+  qs_ancien = 0.
+  prlw_ancien = 0.
+  prsw_ancien = 0.
+  prw_ancien = 0.
   agesno     = 0.
 
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cloudth_mod.F90
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cloudth_mod.F90	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cloudth_mod.F90	(revision 2925)
@@ -311,5 +311,7 @@
       ! (J Jouhaud, JL Dufresne, JB Madeleine)
       REAL,SAVE :: vert_alpha
+      !$OMP THREADPRIVATE(vert_alpha)
       LOGICAL, SAVE :: firstcall = .TRUE.
+      !$OMP THREADPRIVATE(firstcall)
       
       REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
@@ -858,6 +860,8 @@
       ! Change the width of the PDF used for vertical subgrid scale heterogeneity
       ! (J Jouhaud, JL Dufresne, JB Madeleine)
-      REAL,SAVE :: vert_alpha
+      REAL,SAVE :: vert_alpha, vert_alpha_th
+      !$OMP THREADPRIVATE(vert_alpha, vert_alpha_th)
       LOGICAL, SAVE :: firstcall = .TRUE.
+      !$OMP THREADPRIVATE(firstcall)
 
       REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
@@ -893,4 +897,9 @@
         CALL getin_p('cloudth_vert_alpha',vert_alpha)
         WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
+        ! The factor used for the thermal is equal to that of the environment
+        !   if nothing is explicitly specified in the def file
+        vert_alpha_th=vert_alpha
+        CALL getin_p('cloudth_vert_alpha_th',vert_alpha_th)
+        WRITE(*,*) 'cloudth_vert_alpha_th = ', vert_alpha_th
         firstcall=.FALSE.
       ENDIF
@@ -997,5 +1006,5 @@
       qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
 
-      ELSE IF (iflag_cloudth_vert == 3) THEN
+      ELSE IF (iflag_cloudth_vert >= 3) THEN
 
 !-------------------------------------------------------------------------------
@@ -1006,6 +1015,11 @@
 !      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
 !      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
-      deltasenv=aenv*vert_alpha*sigma1s
-      deltasth=ath*vert_alpha*sigma2s
+      IF (iflag_cloudth_vert == 3) THEN
+        deltasenv=aenv*vert_alpha*sigma1s
+        deltasth=ath*vert_alpha_th*sigma2s
+      ELSE IF (iflag_cloudth_vert == 4) THEN
+        deltasenv=vert_alpha*sigma1s
+        deltasth=vert_alpha_th*sigma2s
+      ENDIF
       
       xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
@@ -1049,5 +1063,5 @@
 
 
-      ENDIF ! of if (iflag_cloudth_vert==1 or 3)
+      ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
 
       if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/comsoil.h
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/comsoil.h	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/comsoil.h	(revision 2925)
@@ -3,5 +3,7 @@
 !
 
-      common /comsoil/inertie_sol,inertie_sno,inertie_ice
-      real inertie_sol,inertie_sno,inertie_ice
+      common /comsoil/inertie_sol,inertie_sno,inertie_sic,inertie_lic,  &
+     &                iflag_sic
+      real inertie_sol,inertie_sno,inertie_sic,inertie_lic
+      integer iflag_sic
 !$OMP THREADPRIVATE(/comsoil/)
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/conf_phys_m.F90
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/conf_phys_m.F90	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/conf_phys_m.F90	(revision 2925)
@@ -173,5 +173,7 @@
     REAL,SAVE :: exposant_glace_omp
     REAL,SAVE :: rei_min_omp, rei_max_omp
-    REAL,SAVE :: inertie_sol_omp,inertie_sno_omp,inertie_ice_omp
+    INTEGER,SAVE :: iflag_sic_omp
+    REAL,SAVE :: inertie_sol_omp,inertie_sno_omp,inertie_sic_omp
+    REAL,SAVE :: inertie_lic_omp
     REAL,SAVE :: qsol0_omp
     REAL,SAVE :: evap0_omp
@@ -1137,11 +1139,27 @@
     !-----------------------------------------------------------------------
     !
-    !Config Key  = inertie_ice
+    !Config Key  = iflag_sic
+    !Config Desc =  
+    !Config Def  = 0
+    !Config Help = 
+    !
+    iflag_sic_omp = 0
+    CALL getin('iflag_sic',iflag_sic_omp)
+    !
+    !Config Key  = inertie_sic
     !Config Desc =  
     !Config Def  = 2000.
     !Config Help = 
     !
-    inertie_ice_omp = 2000.
-    CALL getin('inertie_ice',inertie_ice_omp)
+    inertie_sic_omp = 2000.
+    CALL getin('inertie_sic',inertie_sic_omp)
+    !
+    !Config Key  = inertie_lic
+    !Config Desc =  
+    !Config Def  = 2000.
+    !Config Help = 
+    !
+    inertie_lic_omp = 2000.
+    CALL getin('inertie_lic',inertie_lic_omp)
     !
     !Config Key  = inertie_sno
@@ -2154,6 +2172,8 @@
     evap0 = evap0_omp
     albsno0 = albsno0_omp
+    iflag_sic = iflag_sic_omp
     inertie_sol = inertie_sol_omp
-    inertie_ice = inertie_ice_omp
+    inertie_sic = inertie_sic_omp
+    inertie_lic = inertie_lic_omp
     inertie_sno = inertie_sno_omp
     rad_froid = rad_froid_omp
@@ -2572,6 +2592,8 @@
     write(lunout,*)' evap0 = ', evap0
     write(lunout,*)' albsno0 = ', albsno0
+    write(lunout,*)' iflag_sic = ', iflag_sic
     write(lunout,*)' inertie_sol = ', inertie_sol
-    write(lunout,*)' inertie_ice = ', inertie_ice
+    write(lunout,*)' inertie_sic = ', inertie_sic
+    write(lunout,*)' inertie_lic = ', inertie_lic
     write(lunout,*)' inertie_sno = ', inertie_sno
     write(lunout,*)' f_cdrag_ter = ',f_cdrag_ter
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cv3_enthalpmix.F90
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cv3_enthalpmix.F90	(revision 2925)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cv3_enthalpmix.F90	(revision 2925)
@@ -0,0 +1,215 @@
+SUBROUTINE cv3_enthalpmix(len, nd, iflag, plim1, plim2, p, ph, &
+                       t, q, u, v, w, &
+                       wi, nk, tmix, thmix, qmix, qsmix, umix, vmix, plcl)
+  ! **************************************************************
+  ! *
+  ! CV3_ENTHALPMIX   Brassage adiabatique d'une couche d'epaisseur *
+  ! arbitraire.                                   *
+  ! *
+  ! written by   : Grandpeix Jean-Yves, 28/12/2001, 13.14.24    *
+  ! modified by :  Filiberti M-A 06/2005 vectorisation          *
+  ! **************************************************************
+
+  IMPLICIT NONE
+  ! ==============================================================
+
+  ! vertmix : determines theta, t, q, qs, u and v of the mixture generated by
+  ! adiabatic mixing of air between plim1 and plim2 with weighting w.
+  ! If plim1 and plim2 fall within the same model layer, then theta, ... v
+  ! are those of that layer.
+  ! A minimum value (dpmin) is imposed upon plim1-plim2
+
+  ! ===============================================================
+
+  include "cvthermo.h"
+  include "YOETHF.h"
+  include "YOMCST.h"
+  include "FCTTRE.h"
+!inputs:
+  INTEGER, INTENT (IN)                      :: nd, len
+  INTEGER, DIMENSION (len), INTENT (IN)     :: nk
+  REAL, DIMENSION (len), INTENT (IN)        :: plim1, plim2
+  REAL, DIMENSION (len,nd), INTENT (IN)     :: t, q
+  REAL, DIMENSION (len,nd), INTENT (IN)     :: u, v
+  REAL, DIMENSION (nd), INTENT (IN)         :: w
+  REAL, DIMENSION (len,nd), INTENT (IN)     :: p
+  REAL, DIMENSION (len,nd+1), INTENT (IN)   :: ph
+!input/output:
+  INTEGER, DIMENSION (len), INTENT (INOUT)  ::  iflag
+!outputs:
+  REAL, DIMENSION (len), INTENT (OUT)       :: tmix, thmix, qmix
+  REAL, DIMENSION (len), INTENT (OUT)       :: umix, vmix
+  REAL, DIMENSION (len), INTENT (OUT)       :: qsmix
+  REAL, DIMENSION (len), INTENT (OUT)       :: plcl
+  REAL, DIMENSION (len,nd), INTENT (OUT)    :: wi
+!internal variables :
+  INTEGER i, j
+  INTEGER niflag7
+  INTEGER, DIMENSION(len)                   :: j1, j2
+  REAL                                      :: a, b
+  REAL                                      :: cpn
+  REAL                                      :: x, y, p0, p0m1, zdelta, zcor
+  REAL                                      :: dpmin=1.
+!$OMP THREADPRIVATE(dpmin)
+  REAL, DIMENSION(len)                      :: plim2p  ! = min(plim2(:),plim1(:)-dpmin)
+  REAL, DIMENSION(len)                      :: akm     ! mixture enthalpy
+  REAL, DIMENSION(len)                      :: dpw, coef
+  REAL, DIMENSION(len)                      :: rdcp, a2, b2, pnk
+  REAL, DIMENSION(len)                      :: rh, chi
+  REAL, DIMENSION(len)                      :: eqwght
+  REAL, DIMENSION(len,nd)                   :: p1, p2
+
+
+!!  print *,' ->cv3_vertmix, plim1,plim2 ', plim1,plim2   !jyg
+  plim2p(:) = min(plim2(:),plim1(:)-dpmin)
+  j1(:)=nd
+  j2(:) = 0
+  DO j = 1, nd
+    DO i = 1, len
+      IF (plim1(i)<=ph(i,j)) j1(i) = j
+!!!      IF (plim2p(i)>=ph(i,j+1) .AND. plim2p(i)<ph(i,j)) j2(i) = j
+      IF (plim2p(i)< ph(i,j)) j2(i) = j
+    END DO
+  END DO
+
+  DO j = 1, nd
+    DO i = 1, len
+      wi(i, j) = 0.
+    END DO
+  END DO
+  DO i = 1, len
+    akm(i) = 0.
+    qmix(i) = 0.
+    umix(i) = 0.
+    vmix(i) = 0.
+    dpw(i) = 0.
+    a2(i) = 0.0
+    b2(i) = 0.
+    pnk(i) = p(i, nk(i))
+  END DO
+  eqwght(:) = 0.
+
+  p0 = 1000.
+  p0m1 = 1./p0
+
+  DO i = 1, len
+    IF (j2(i) < j1(i)) THEN
+      coef(i) = 1.
+      eqwght(i) = 1.
+    ELSE
+      coef(i) = 1./(plim1(i)-plim2p(i))
+    ENDIF
+  END DO
+
+!!  print *,'cv3_vertmix, j1,j2,coef ', j1,j2,coef  !jyg
+
+  DO j = 1, nd
+    DO i = 1, len
+      IF (j>=j1(i) .AND. j<=j2(i)) THEN
+        p1(i, j) = min(ph(i,j), plim1(i))
+        p2(i, j) = max(ph(i,j+1), plim2p(i))
+        ! CRtest:couplage thermiques: deja normalise
+        ! wi(i,j) = w(j)
+        ! print*,'wi',wi(i,j)
+        wi(i, j) = w(j)*(p1(i,j)-p2(i,j))*coef(i)+eqwght(i)
+        dpw(i) = dpw(i) + wi(i, j)
+
+!!  print *,'cv3_vertmix, j, wi(1,j),dpw ', j, wi(1,j),dpw  !jyg
+
+      END IF
+    END DO
+  END DO
+
+  ! CR:print
+  ! do i=1,len
+  ! print*,'plim',plim1(i),plim2p(i)
+  ! enddo
+  DO j = 1, nd
+    DO i = 1, len
+      IF (j>=j1(i) .AND. j<=j2(i)) THEN
+        wi(i, j) = wi(i, j)/dpw(i)
+        akm(i) = akm(i) + (cpd*(1.-q(i,j))+q(i,j)*cpv)*t(i, j)*wi(i, j)
+        qmix(i) = qmix(i) + q(i, j)*wi(i, j)
+        umix(i) = umix(i) + u(i, j)*wi(i, j)
+        vmix(i) = vmix(i) + v(i, j)*wi(i, j)
+      END IF
+    END DO
+  END DO
+
+  DO i = 1, len
+    rdcp(i) = (rrd*(1.-qmix(i))+qmix(i)*rrv)/(cpd*(1.-qmix(i))+qmix(i)*cpv)
+  END DO
+
+
+!!  print *,'cv3_vertmix, rdcp ', rdcp  !jyg
+
+
+
+  DO j = 1, nd
+    DO i = 1, len
+      IF (j>=j1(i) .AND. j<=j2(i)) THEN
+        ! c            x=(.5*(p1(i,j)+p2(i,j))*p0m1)**rdcp(i)
+        y = (.5*(p1(i,j)+p2(i,j))/pnk(i))**rdcp(i)
+        ! c            a2(i)=a2(i)+(cpd*(1.-qmix(i))+qmix(i)*cpv)*x*wi(i,j)
+        b2(i) = b2(i) + (cpd*(1.-qmix(i))+qmix(i)*cpv)*y*wi(i, j)
+      END IF
+    END DO
+  END DO
+
+  DO i = 1, len
+    tmix(i) = akm(i)/b2(i)
+    thmix(i) = tmix(i)*(p0/pnk(i))**rdcp(i)
+    ! print*,'thmix akm',akm(i),b2(i)
+    ! print*,'thmix t',tmix(i),p0
+    ! print*,'thmix p',pnk(i),rdcp(i)
+    ! print*,'thmix',thmix(i)
+    ! c         thmix(i) = akm(i)/a2(i)
+    ! c         tmix(i)= thmix(i)*(pnk(i)*p0m1)**rdcp(i)
+    zdelta = max(0., sign(1.,rtt-tmix(i)))
+    qsmix(i) = r2es*foeew(tmix(i), zdelta)/(pnk(i)*100.)
+    qsmix(i) = min(0.5, qsmix(i))
+    zcor = 1./(1.-retv*qsmix(i))
+    qsmix(i) = qsmix(i)*zcor
+  END DO
+
+  ! -------------------------------------------------------------------
+  ! --- Calculate lifted condensation level of air at parcel origin level
+  ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
+  ! -------------------------------------------------------------------
+
+  a = 1669.0 ! convect3
+  b = 122.0 ! convect3
+
+
+  niflag7 = 0
+  DO i = 1, len
+
+    IF (iflag(i)/=7) THEN ! modif sb Jun7th 2002
+
+      rh(i) = qmix(i)/qsmix(i)
+      chi(i) = tmix(i)/(a-b*rh(i)-tmix(i)) ! convect3
+      ! ATTENTION, la LIGNE DESSOUS A ETE RAJOUTEE ARBITRAIREMENT ET
+      ! MASQUE UN PB POTENTIEL
+      chi(i) = max(chi(i), 0.)
+      rh(i) = max(rh(i), 0.)
+      plcl(i) = pnk(i)*(rh(i)**chi(i))
+      IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) &
+          iflag(i) = 8
+
+    ELSE
+
+      niflag7 = niflag7 + 1
+      plcl(i) = plim2p(i)
+
+    END IF ! iflag=7
+
+    ! print*,'NIFLAG7  =',niflag7
+
+  END DO
+
+!!  print *,' cv3_vertmix->'  !jyg
+
+
+  RETURN
+END SUBROUTINE cv3_enthalpmix
+
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cv3_estatmix.F90
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cv3_estatmix.F90	(revision 2925)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cv3_estatmix.F90	(revision 2925)
@@ -0,0 +1,198 @@
+SUBROUTINE cv3_estatmix(len, nd, iflag, plim1, plim2, p, ph, &
+                       t, q, u, v, h, gz, w, &
+                       wi, nk, tmix, thmix, qmix, qsmix, umix, vmix, plcl)
+  ! **************************************************************
+  ! *
+  ! CV3_ESTATMIX  Determine the properties of an adiabatic updraft  *
+  !                made of air coming from several layers by        *
+  !                mixing static energy                             *
+  !                                                                 *
+  ! written by   : Grandpeix Jean-Yves, 28/12/2001, 13.14.24        *
+  ! modified by :  Filiberti M-A 06/2005 vectorisation              *
+  ! ****************************************************************
+
+  IMPLICIT NONE
+  ! ==============================================================
+
+  ! estatmix : determines theta, t, q, qs, u and v of the lifted mixture
+  ! made of air between plim1 and plim2 with weighting w.
+  ! If plim1 and plim2 fall within the same model layer, then theta, ... v
+  ! are those of that layer.
+  ! A minimum value (dpmin) is imposed upon plim1-plim2
+
+  ! ===============================================================
+
+  include "cvthermo.h"
+  include "YOETHF.h"
+  include "YOMCST.h"
+  include "FCTTRE.h"
+!inputs:
+  INTEGER, INTENT (IN)                      :: nd, len
+  INTEGER, DIMENSION (len), INTENT (IN)     :: nk
+  REAL, DIMENSION (len), INTENT (IN)        :: plim1, plim2
+  REAL, DIMENSION (len,nd), INTENT (IN)     :: t, q
+  REAL, DIMENSION (len,nd), INTENT (IN)     :: u, v
+  REAL, DIMENSION (len,nd), INTENT (IN)     :: h ! static energy of the layers
+  REAL, DIMENSION (len,nd), INTENT (IN)     :: gz
+  REAL, DIMENSION (nd), INTENT (IN)         :: w
+  REAL, DIMENSION (len,nd), INTENT (IN)     :: p
+  REAL, DIMENSION (len,nd+1), INTENT (IN)   :: ph
+!input/output:
+  INTEGER, DIMENSION (len), INTENT (INOUT)  ::  iflag
+!outputs:
+  REAL, DIMENSION (len), INTENT (OUT)       :: tmix, thmix, qmix
+  REAL, DIMENSION (len), INTENT (OUT)       :: umix, vmix
+  REAL, DIMENSION (len), INTENT (OUT)       :: qsmix
+  REAL, DIMENSION (len), INTENT (OUT)       :: plcl
+  REAL, DIMENSION (len,nd), INTENT (OUT)    :: wi
+!internal variables :
+  INTEGER i, j
+  INTEGER niflag7
+  INTEGER, DIMENSION(len)                   :: j1, j2
+  REAL                                      :: a, b
+  REAL                                      :: cpn
+  REAL                                      :: x, y, p0, zdelta, zcor
+  REAL                                      :: dpmin=1.
+!$OMP THREADPRIVATE(dpmin)
+  REAL, DIMENSION(len)                      :: plim2p  ! = min(plim2(:),plim1(:)-dpmin)
+  REAL, DIMENSION(len)                      :: dpw, coef
+  REAL, DIMENSION(len)                      :: hmix ! static energy of the updraft
+  REAL, DIMENSION(len)                      :: rdcp, pnk
+  REAL, DIMENSION(len)                      :: rh, chi
+  REAL, DIMENSION(len)                      :: eqwght
+  REAL, DIMENSION(len,nd)                   :: p1, p2
+
+
+!!  print *,' ->cv3_vertmix, plim1,plim2 ', plim1,plim2   !jyg
+  plim2p(:) = min(plim2(:),plim1(:)-dpmin)
+  j1(:)=nd
+  j2(:) = 0
+  DO j = 1, nd
+    DO i = 1, len
+      IF (plim1(i)<=ph(i,j)) j1(i) = j
+!!!      IF (plim2p(i)>=ph(i,j+1) .AND. plim2p(i)<ph(i,j)) j2(i) = j
+      IF (plim2p(i)< ph(i,j)) j2(i) = j
+    END DO
+  END DO
+
+  DO j = 1, nd
+    DO i = 1, len
+      wi(i, j) = 0.
+    END DO
+  END DO
+  DO i = 1, len
+    hmix(i) = 0.
+    qmix(i) = 0.
+    umix(i) = 0.
+    vmix(i) = 0.
+    dpw(i) = 0.
+    pnk(i) = p(i, nk(i))
+  END DO
+  eqwght(:) = 0.
+
+  p0 = 1000.
+
+  DO i = 1, len
+    IF (j2(i) < j1(i)) THEN
+      coef(i) = 1.
+      eqwght(i) = 1.
+    ELSE
+      coef(i) = 1./(plim1(i)-plim2p(i))
+    ENDIF
+  END DO
+
+!!  print *,'cv3_vertmix, j1,j2,coef ', j1,j2,coef  !jyg
+
+  DO j = 1, nd
+    DO i = 1, len
+      IF (j>=j1(i) .AND. j<=j2(i)) THEN
+        p1(i, j) = min(ph(i,j), plim1(i))
+        p2(i, j) = max(ph(i,j+1), plim2p(i))
+        ! CRtest:couplage thermiques: deja normalise
+        ! wi(i,j) = w(j)
+        ! print*,'wi',wi(i,j)
+        wi(i, j) = w(j)*(p1(i,j)-p2(i,j))*coef(i)+eqwght(i)
+        dpw(i) = dpw(i) + wi(i, j)
+
+!!  print *,'cv3_vertmix, j, wi(1,j),dpw ', j, wi(1,j),dpw  !jyg
+
+      END IF
+    END DO
+  END DO
+
+  ! CR:print
+  ! do i=1,len
+  ! print*,'plim',plim1(i),plim2p(i)
+  ! enddo
+  DO j = 1, nd
+    DO i = 1, len
+      IF (j>=j1(i) .AND. j<=j2(i)) THEN
+        wi(i, j) = wi(i, j)/dpw(i)
+        hmix(i) = hmix(i) + h(i, j)*wi(i, j)
+        qmix(i) = qmix(i) +  q(i, j)*wi(i, j)
+        umix(i) = umix(i) +  u(i, j)*wi(i, j)
+        vmix(i) = vmix(i) +  v(i, j)*wi(i, j)
+      END IF
+    END DO
+  END DO
+
+  DO i = 1, len
+    rdcp(i) = (rrd*(1.-qmix(i))+qmix(i)*rrv)/(cpd*(1.-qmix(i))+qmix(i)*cpv)
+  END DO
+
+
+!!  print *,'cv3_vertmix, rdcp ', rdcp  !jyg
+
+  DO i = 1, len
+    tmix(i) = (hmix(i) - gz(i,1))/(cpd*(1.-qmix(i)) + qmix(i)*cpv)
+    !      (Use of Cpv since we are dealing with dry static energy)
+    thmix(i) = tmix(i)*(p0/pnk(i))**rdcp(i)
+    ! print*,'tmix thmix hmix ',tmix(i),thmix(i),hmix(i)
+    zdelta = max(0., sign(1.,rtt-tmix(i)))
+    qsmix(i) = r2es*foeew(tmix(i), zdelta)/(pnk(i)*100.)
+    qsmix(i) = min(0.5, qsmix(i))
+    zcor = 1./(1.-retv*qsmix(i))
+    qsmix(i) = qsmix(i)*zcor
+  END DO
+
+  ! -------------------------------------------------------------------
+  ! --- Calculate lifted condensation level of air at parcel origin level
+  ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
+  ! -------------------------------------------------------------------
+
+  a = 1669.0 ! convect3
+  b = 122.0 ! convect3
+
+
+  niflag7 = 0
+  DO i = 1, len
+
+    IF (iflag(i)/=7) THEN ! modif sb Jun7th 2002
+
+      rh(i) = qmix(i)/qsmix(i)
+      chi(i) = tmix(i)/(a-b*rh(i)-tmix(i)) ! convect3
+      ! ATTENTION, la LIGNE DESSOUS A ETE RAJOUTEE ARBITRAIREMENT ET
+      ! MASQUE UN PB POTENTIEL
+      chi(i) = max(chi(i), 0.)
+      rh(i) = max(rh(i), 0.)
+      plcl(i) = pnk(i)*(rh(i)**chi(i))
+      IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) &
+          iflag(i) = 8
+
+    ELSE
+
+      niflag7 = niflag7 + 1
+      plcl(i) = plim2p(i)
+
+    END IF ! iflag=7
+
+    ! print*,'NIFLAG7  =',niflag7
+
+  END DO
+
+!!  print *,' cv3_vertmix->'  !jyg
+
+
+  RETURN
+END SUBROUTINE cv3_estatmix
+
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cv3_routines.F90
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cv3_routines.F90	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cv3_routines.F90	(revision 2925)
@@ -111,4 +111,9 @@
      ok_optim_yield=.False.
      CALL getin_p('ok_optim_yield',ok_optim_yield)
+     ok_homo_tend=.TRUE.
+     CALL getin_p('ok_homo_tend',ok_homo_tend)
+     ok_entrain=.TRUE.
+     CALL getin_p('ok_entrain',ok_entrain)
+
      coef_peel=0.25
      CALL getin_p('coef_peel',coef_peel)
@@ -277,5 +282,5 @@
 
 SUBROUTINE cv3_feed(len, nd, ok_conserv_q, &
-                    t, q, u, v, p, ph, hm, gz, &
+                    t, q, u, v, p, ph, h, gz, &
                     p1feed, p2feed, wght, &
                     wghti, tnk, thnk, qnk, qsnk, unk, vnk, &
@@ -283,4 +288,5 @@
 
   USE mod_phys_lmdz_transfert_para, ONLY : bcast
+  USE add_phys_tend_mod, ONLY: fl_cor_ebil
   IMPLICIT NONE
 
@@ -292,4 +298,5 @@
 ! - here, nk(i)=minorig
 ! - icb defined differently (plcl compared with ph instead of p)
+! - dry static energy as argument instead of moist static energy
 
 ! Main differences with convect3:
@@ -307,5 +314,5 @@
   REAL, DIMENSION (len, nd), INTENT (IN)             :: t, q, p
   REAL, DIMENSION (len, nd), INTENT (IN)             :: u, v
-  REAL, DIMENSION (len, nd), INTENT (IN)             :: hm, gz
+  REAL, DIMENSION (len, nd), INTENT (IN)             :: h, gz
   REAL, DIMENSION (len, nd+1), INTENT (IN)           :: ph
   REAL, DIMENSION (len), INTENT (IN)                 :: p1feed
@@ -378,14 +385,26 @@
     pup(i) = p2feed(i)
   END DO
-  CALL cv3_vertmix(len, nd, iflag, p1feed, pup, p, ph, &
-                   t, q, u, v, wght, &
-                   wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclup)
+  IF (fl_cor_ebil >=2 ) THEN
+    CALL cv3_estatmix(len, nd, iflag, p1feed, pup, p, ph, &
+                     t, q, u, v, h, gz, wght, &
+                     wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclup)
+  ELSE
+    CALL cv3_enthalpmix(len, nd, iflag, p1feed, pup, p, ph, &
+                       t, q, u, v, wght, &
+                       wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclup)
+  ENDIF  ! (fl_cor_ebil >=2 ) 
 ! 1.b- LCL associated with ph(nk+1)
   DO i = 1, len
     plo(i) = ph(i, nk(i)+1)
   END DO
-  CALL cv3_vertmix(len, nd, iflag, p1feed, plo, p, ph, &
-                   t, q, u, v, wght, &
-                   wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plcllo)
+  IF (fl_cor_ebil >=2 ) THEN
+    CALL cv3_estatmix(len, nd, iflag, p1feed, plo, p, ph, &
+                     t, q, u, v, h, gz, wght, &
+                     wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plcllo)
+  ELSE
+    CALL cv3_enthalpmix(len, nd, iflag, p1feed, plo, p, ph, &
+                       t, q, u, v, wght, &
+                       wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plcllo)
+  ENDIF  ! (fl_cor_ebil >=2 ) 
 ! 2- Iterations
   niter = 5
@@ -434,7 +453,13 @@
 !jyg>
 
-    CALL cv3_vertmix(len, nd, iflag, p1feed, pfeed, p, ph, &
-                   t, q, u, v, wght, &
-                   wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclfeed)
+    IF (fl_cor_ebil >=2 ) THEN
+      CALL cv3_estatmix(len, nd, iflag, p1feed, pfeed, p, ph, &
+                       t, q, u, v, h, gz, wght, &
+                       wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclfeed)
+    ELSE
+      CALL cv3_enthalpmix(len, nd, iflag, p1feed, pfeed, p, ph, &
+                         t, q, u, v, wght, &
+                         wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclfeed)
+    ENDIF  ! (fl_cor_ebil >=2 ) 
 !jyg20140217<
     IF (ok_new_feed) THEN
@@ -1647,4 +1672,10 @@
       DO i = 1, ncum
         IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
+!jyg<   (energy conservation tests)
+!!          hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)*tp(i,k))*ep(i, k)*clw(i, k)
+!!          hp(i, k) = ( hnk(i) + (lv(i,k)+(cpd-cpv)*t(i,k))*ep(i, k)*clw(i, k) ) / &
+!!                     (1. - ep(i,k)*clw(i,k))
+!!          hp(i, k) = ( hnk(i) + (lv(i,k)+(cpd-cl)*t(i,k))*ep(i, k)*clw(i, k) ) / &
+!!                     (1. - ep(i,k)*clw(i,k))
           hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)*t(i,k))*ep(i, k)*clw(i, k)
         END IF
@@ -2984,4 +3015,7 @@
                      qcondc, wd, &
                      ftd, fqd, qnk, qtc, sigt, tau_cld_cv, coefw_cld_cv)
+
+    USE print_control_mod, ONLY: lunout, prt_level
+    USE add_phys_tend_mod, only : fl_cor_ebil
 
   IMPLICIT NONE
@@ -3220,6 +3254,13 @@
 
       IF ((0.01*grav*work(il)*am(il))>=delti) iflag(il) = 1 !consist vect
-      ft(il, 1) = ft(il, 1) + 0.01*grav*work(il)*am(il) * &
-                                   (t(il,2)-t(il,1)+(gz(il,2)-gz(il,1))/cpn(il,1))
+!jyg<
+        IF (fl_cor_ebil >= 2) THEN
+          ft(il, 1) = ft(il, 1) + 0.01*grav*work(il)*am(il) * &
+                    ((t(il,2)-t(il,1))*cpn(il,2)+gz(il,2)-gz(il,1))/cpn(il,1)
+        ELSE
+          ft(il, 1) = ft(il, 1) + 0.01*grav*work(il)*am(il) * &
+                    (t(il,2)-t(il,1)+(gz(il,2)-gz(il,1))/cpn(il,1))
+        ENDIF
+!>jyg
     END IF ! iflag
   END DO
@@ -3529,7 +3570,15 @@
 
 ! sature
-        ft(il, i) = ft(il, i) + 0.01*grav*dpinv * &
+!jyg<
+        IF (fl_cor_ebil >= 2) THEN
+          ft(il, i) = ft(il, i) + 0.01*grav*dpinv * &
+              ( amp1(il)*( (t(il,i+1)-t(il,i))*cpn(il,i+1) + gz(il,i+1)-gz(il,i))*cpinv - &
+                ad(il)*( (t(il,i)-t(il,i-1))*cpn(il,i-1) + gz(il,i)-gz(il,i-1))*cpinv)
+        ELSE
+          ft(il, i) = ft(il, i) + 0.01*grav*dpinv * &
                      (amp1(il)*(t(il,i+1)-t(il,i) + (gz(il,i+1)-gz(il,i))*cpinv) - &
                       ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv))
+        ENDIF
+!>jyg
 
 
@@ -3538,7 +3587,5 @@
                                     t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
         END IF
-
-
-
+!
 ! sb: on ne fait pas encore la correction permettant de mieux
 ! conserver l'eau:
@@ -3873,14 +3920,16 @@
 
 !!!!      do 700 i=1,icb(il)-1
-  DO i = 1, nl
-    DO il = 1, ncum
-      IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
-        ftd(il, i) = esum(il)*t_wake(il, i)/(th_wake(il,i)*hsum(il))
-        fqd(il, i) = fsum(il)/gsum(il)
-        ft(il, i) = ftd(il, i) + asum(il)*t(il, i)/(th(il,i)*dsum(il))
-        fr(il, i) = fqd(il, i) + bsum(il)/csum(il)
-      END IF
-    END DO
-  END DO
+  IF (ok_homo_tend) THEN
+    DO i = 1, nl
+      DO il = 1, ncum
+        IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
+          ftd(il, i) = esum(il)*t_wake(il, i)/(th_wake(il,i)*hsum(il))
+          fqd(il, i) = fsum(il)/gsum(il)
+          ft(il, i) = ftd(il, i) + asum(il)*t(il, i)/(th(il,i)*dsum(il))
+          fr(il, i) = fqd(il, i) + bsum(il)/csum(il)
+        END IF
+      END DO
+    END DO
+  ENDIF
 
 !jyg<
@@ -3920,5 +3969,7 @@
   END DO
 !
-!    print *,' YIELD : alpha_qpos ',alpha_qpos(1)
+    IF (prt_level .GE. 5) THEN
+      print *,' CV3_YIELD : alpha_qpos ',alpha_qpos(1)
+    ENDIF
 !
   DO il = 1, ncum
Index: DZ5/branches/LMDZ_tree_FC/libf/phylmd/cv3_vertmix.F90
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cv3_vertmix.F90	(revision 2924)
+++ 	(revision )
@@ -1,214 +1,0 @@
-SUBROUTINE cv3_vertmix(len, nd, iflag, plim1, plim2, p, ph, &
-                       t, q, u, v, w, &
-                       wi, nk, tmix, thmix, qmix, qsmix, umix, vmix, plcl)
-  ! **************************************************************
-  ! *
-  ! CV3_VERTMIX   Brassage adiabatique d'une couche d'epaisseur *
-  ! arbitraire.                                   *
-  ! *
-  ! written by   : Grandpeix Jean-Yves, 28/12/2001, 13.14.24    *
-  ! modified by :  Filiberti M-A 06/2005 vectorisation          *
-  ! **************************************************************
-
-  IMPLICIT NONE
-  ! ==============================================================
-
-  ! vertmix : determines theta, t, q, qs, u and v of the mixture generated by
-  ! adiabatic mixing of air between plim1 and plim2 with weighting w.
-  ! If plim1 and plim2 fall within the same model layer, then theta, ... v
-  ! are those of that layer.
-  ! A minimum value (dpmin) is imposed upon plim1-plim2
-
-  ! ===============================================================
-
-  include "cvthermo.h"
-  include "YOETHF.h"
-  include "YOMCST.h"
-  include "FCTTRE.h"
-!inputs:
-  INTEGER, INTENT (IN)                      ::  nd, len
-  INTEGER, DIMENSION (len), INTENT (IN)     ::  nk
-  REAL, DIMENSION (nd), INTENT (IN)        ::  w
-  REAL, DIMENSION (len), INTENT (IN)        :: plim1, plim2
-  REAL, DIMENSION (len,nd), INTENT (IN)     :: t, q
-  REAL, DIMENSION (len,nd), INTENT (IN)     :: u, v
-  REAL, DIMENSION (len,nd), INTENT (IN)     :: p
-  REAL, DIMENSION (len,nd+1), INTENT (IN)   :: ph
-!input/output:
-  INTEGER, DIMENSION (len), INTENT (INOUT)  ::  iflag
-!outputs:
-  REAL, DIMENSION (len), INTENT (OUT)       :: tmix, thmix, qmix
-  REAL, DIMENSION (len), INTENT (OUT)       :: umix, vmix
-  REAL, DIMENSION (len), INTENT (OUT)       :: qsmix
-  REAL, DIMENSION (len), INTENT (OUT)       :: plcl
-  REAL, DIMENSION (len,nd), INTENT (OUT)    :: wi
-!internal variables :
-  INTEGER i, j
-  INTEGER niflag7
-  INTEGER, DIMENSION(len)                   :: j1, j2
-  REAL                                      :: a, b
-  REAL                                      :: cpn
-  REAL                                      :: x, y, p0, p0m1, zdelta, zcor
-  REAL                                      :: dpmin=1.
-!$OMP THREADPRIVATE(dpmin)
-  REAL, DIMENSION(len)                      :: plim2p  ! = min(plim2(:),plim1(:)-dpmin)
-  REAL, DIMENSION(len)                      :: ahm, dpw, coef
-  REAL, DIMENSION(len)                      :: rdcp, a2, b2, pnk
-  REAL, DIMENSION(len)                      :: rh, chi
-  REAL, DIMENSION(len)                      :: eqwght
-  REAL, DIMENSION(len,nd)                   :: p1, p2
-
-
-!!  print *,' ->cv3_vertmix, plim1,plim2 ', plim1,plim2   !jyg
-  plim2p(:) = min(plim2(:),plim1(:)-dpmin)
-  j1(:)=nd
-  j2(:) = 0
-  DO j = 1, nd
-    DO i = 1, len
-      IF (plim1(i)<=ph(i,j)) j1(i) = j
-!!!      IF (plim2p(i)>=ph(i,j+1) .AND. plim2p(i)<ph(i,j)) j2(i) = j
-      IF (plim2p(i)< ph(i,j)) j2(i) = j
-    END DO
-  END DO
-
-  DO j = 1, nd
-    DO i = 1, len
-      wi(i, j) = 0.
-    END DO
-  END DO
-  DO i = 1, len
-    ahm(i) = 0.
-    qmix(i) = 0.
-    umix(i) = 0.
-    vmix(i) = 0.
-    dpw(i) = 0.
-    a2(i) = 0.0
-    b2(i) = 0.
-    pnk(i) = p(i, nk(i))
-  END DO
-  eqwght(:) = 0.
-
-  p0 = 1000.
-  p0m1 = 1./p0
-
-  DO i = 1, len
-    IF (j2(i) < j1(i)) THEN
-      coef(i) = 1.
-      eqwght(i) = 1.
-    ELSE
-      coef(i) = 1./(plim1(i)-plim2p(i))
-    ENDIF
-  END DO
-
-!!  print *,'cv3_vertmix, j1,j2,coef ', j1,j2,coef  !jyg
-
-  DO j = 1, nd
-    DO i = 1, len
-      IF (j>=j1(i) .AND. j<=j2(i)) THEN
-        p1(i, j) = min(ph(i,j), plim1(i))
-        p2(i, j) = max(ph(i,j+1), plim2p(i))
-        ! CRtest:couplage thermiques: deja normalise
-        ! wi(i,j) = w(j)
-        ! print*,'wi',wi(i,j)
-        wi(i, j) = w(j)*(p1(i,j)-p2(i,j))*coef(i)+eqwght(i)
-        dpw(i) = dpw(i) + wi(i, j)
-
-!!  print *,'cv3_vertmix, j, wi(1,j),dpw ', j, wi(1,j),dpw  !jyg
-
-      END IF
-    END DO
-  END DO
-
-  ! CR:print
-  ! do i=1,len
-  ! print*,'plim',plim1(i),plim2p(i)
-  ! enddo
-  DO j = 1, nd
-    DO i = 1, len
-      IF (j>=j1(i) .AND. j<=j2(i)) THEN
-        wi(i, j) = wi(i, j)/dpw(i)
-        ahm(i) = ahm(i) + (cpd*(1.-q(i,j))+q(i,j)*cpv)*t(i, j)*wi(i, j)
-        qmix(i) = qmix(i) + q(i, j)*wi(i, j)
-        umix(i) = umix(i) + u(i, j)*wi(i, j)
-        vmix(i) = vmix(i) + v(i, j)*wi(i, j)
-      END IF
-    END DO
-  END DO
-
-  DO i = 1, len
-    rdcp(i) = (rrd*(1.-qmix(i))+qmix(i)*rrv)/(cpd*(1.-qmix(i))+qmix(i)*cpv)
-  END DO
-
-
-!!  print *,'cv3_vertmix, rdcp ', rdcp  !jyg
-
-
-
-  DO j = 1, nd
-    DO i = 1, len
-      IF (j>=j1(i) .AND. j<=j2(i)) THEN
-        ! c            x=(.5*(p1(i,j)+p2(i,j))*p0m1)**rdcp(i)
-        y = (.5*(p1(i,j)+p2(i,j))/pnk(i))**rdcp(i)
-        ! c            a2(i)=a2(i)+(cpd*(1.-qmix(i))+qmix(i)*cpv)*x*wi(i,j)
-        b2(i) = b2(i) + (cpd*(1.-qmix(i))+qmix(i)*cpv)*y*wi(i, j)
-      END IF
-    END DO
-  END DO
-
-  DO i = 1, len
-    tmix(i) = ahm(i)/b2(i)
-    thmix(i) = tmix(i)*(p0/pnk(i))**rdcp(i)
-    ! print*,'thmix ahm',ahm(i),b2(i)
-    ! print*,'thmix t',tmix(i),p0
-    ! print*,'thmix p',pnk(i),rdcp(i)
-    ! print*,'thmix',thmix(i)
-    ! c         thmix(i) = ahm(i)/a2(i)
-    ! c         tmix(i)= thmix(i)*(pnk(i)*p0m1)**rdcp(i)
-    zdelta = max(0., sign(1.,rtt-tmix(i)))
-    qsmix(i) = r2es*foeew(tmix(i), zdelta)/(pnk(i)*100.)
-    qsmix(i) = min(0.5, qsmix(i))
-    zcor = 1./(1.-retv*qsmix(i))
-    qsmix(i) = qsmix(i)*zcor
-  END DO
-
-  ! -------------------------------------------------------------------
-  ! --- Calculate lifted condensation level of air at parcel origin level
-  ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
-  ! -------------------------------------------------------------------
-
-  a = 1669.0 ! convect3
-  b = 122.0 ! convect3
-
-
-  niflag7 = 0
-  DO i = 1, len
-
-    IF (iflag(i)/=7) THEN ! modif sb Jun7th 2002
-
-      rh(i) = qmix(i)/qsmix(i)
-      chi(i) = tmix(i)/(a-b*rh(i)-tmix(i)) ! convect3
-      ! ATTENTION, la LIGNE DESSOUS A ETE RAJOUTEE ARBITRAIREMENT ET
-      ! MASQUE UN PB POTENTIEL
-      chi(i) = max(chi(i), 0.)
-      rh(i) = max(rh(i), 0.)
-      plcl(i) = pnk(i)*(rh(i)**chi(i))
-      IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) &
-          iflag(i) = 8
-
-    ELSE
-
-      niflag7 = niflag7 + 1
-      plcl(i) = plim2p(i)
-
-    END IF ! iflag=7
-
-    ! print*,'NIFLAG7  =',niflag7
-
-  END DO
-
-!!  print *,' cv3_vertmix->'  !jyg
-
-
-  RETURN
-END SUBROUTINE cv3_vertmix
-
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cv3p_mixing.F90
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cv3p_mixing.F90	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cv3p_mixing.F90	(revision 2925)
@@ -14,4 +14,6 @@
 
   USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
+  USE ioipsl_getin_p_mod, ONLY: getin_p
+  USE add_phys_tend_mod, ONLY: fl_cor_ebil
 
   IMPLICIT NONE
@@ -110,5 +112,5 @@
             fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max
 !>jyg
-
+!
   END IF
 
@@ -165,58 +167,64 @@
   DO i = minorig + 1, nl
 
-    DO j = minorig, nl
-      DO il = 1, ncum
-        IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) &
-                         .AND. (j<=inb(il))) THEN
-
-          rti = qnk(il) - ep(il, i)*clw(il, i)
-          bf2 = 1. + lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
-!jyg(from aj)<
-          IF (cvflag_ice) THEN
-! print*,cvflag_ice,'cvflag_ice dans do 700'
-            IF (t(il,j)<=263.15) THEN
-              bf2 = 1. + (lf(il,j)+lv(il,j))*(lv(il,j)+frac(il,j)* &
-                   lf(il,j))*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
-            END IF
-          END IF
-!>jyg
-          anum = h(il, j) - hp(il, i) + (cpv-cpd)*t(il, j)*(rti-rr(il,j))
-          denom = h(il, i) - hp(il, i) + (cpd-cpv)*(rr(il,i)-rti)*t(il, j)
-          dei = denom
-          IF (abs(dei)<0.01) dei = 0.01
-          Sij(il, i, j) = anum/dei
-          Sij(il, i, i) = 1.0
-          altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il,i,j))*rti - rs(il, j)
-          altem = altem/bf2
-          cwat = clw(il, j)*(1.-ep(il,j))
-          stemp = Sij(il, i, j)
-          IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
+    IF (ok_entrain) THEN
+      DO j = minorig, nl
+        DO il = 1, ncum
+          IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) &
+                           .AND. (j<=inb(il))) THEN
+
+            rti = qnk(il) - ep(il, i)*clw(il, i)
+            bf2 = 1. + lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
 !jyg(from aj)<
             IF (cvflag_ice) THEN
-              anum = anum - (lv(il,j)+frac(il,j)*lf(il,j))*(rti-rs(il,j)-cwat*bf2)
-              denom = denom + (lv(il,j)+frac(il,j)*lf(il,j))*(rr(il,i)-rti)
-            ELSE
-              anum = anum - lv(il, j)*(rti-rs(il,j)-cwat*bf2)
-              denom = denom + lv(il, j)*(rr(il,i)-rti)
+! print*,cvflag_ice,'cvflag_ice dans do 700'
+              IF (t(il,j)<=263.15) THEN
+                bf2 = 1. + (lf(il,j)+lv(il,j))*(lv(il,j)+frac(il,j)* &
+                     lf(il,j))*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
+              END IF
             END IF
 !>jyg
-            IF (abs(denom)<0.01) denom = 0.01
-            Sij(il, i, j) = anum/denom
+            anum = h(il, j) - hp(il, i) + (cpv-cpd)*t(il, j)*(rti-rr(il,j))
+            denom = h(il, i) - hp(il, i) + (cpd-cpv)*(rr(il,i)-rti)*t(il, j)
+            dei = denom
+            IF (abs(dei)<0.01) dei = 0.01
+            Sij(il, i, j) = anum/dei
+            Sij(il, i, i) = 1.0
             altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il,i,j))*rti - rs(il, j)
-            altem = altem - (bf2-1.)*cwat
-          END IF
-          IF (Sij(il,i,j)>0.0) THEN
+            altem = altem/bf2
+            cwat = clw(il, j)*(1.-ep(il,j))
+            stemp = Sij(il, i, j)
+            IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
+!jyg(from aj)<
+              IF (cvflag_ice) THEN
+                anum = anum - (lv(il,j)+frac(il,j)*lf(il,j))*(rti-rs(il,j)-cwat*bf2)
+                denom = denom + (lv(il,j)+frac(il,j)*lf(il,j))*(rr(il,i)-rti)
+              ELSE
+                anum = anum - lv(il, j)*(rti-rs(il,j)-cwat*bf2)
+                denom = denom + lv(il, j)*(rr(il,i)-rti)
+              END IF
+!>jyg
+              IF (abs(denom)<0.01) denom = 0.01
+              Sij(il, i, j) = anum/denom
+              altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il,i,j))*rti - rs(il, j)
+              altem = altem - (bf2-1.)*cwat
+            END IF
+            IF (Sij(il,i,j)>0.0) THEN
 !!!                 Ment(il,i,j)=m(il,i)
-            Ment(il, i, j) = 1.
-            elij(il, i, j) = altem
-            elij(il, i, j) = amax1(0.0, elij(il,i,j))
-            nent(il, i) = nent(il, i) + 1
-          END IF
-
-          Sij(il, i, j) = amax1(0.0, Sij(il,i,j))
-          Sij(il, i, j) = amin1(1.0, Sij(il,i,j))
-        END IF ! new
-      END DO
-    END DO
+              Ment(il, i, j) = 1.
+              elij(il, i, j) = altem
+              elij(il, i, j) = amax1(0.0, elij(il,i,j))
+              nent(il, i) = nent(il, i) + 1
+            END IF
+
+            Sij(il, i, j) = amax1(0.0, Sij(il,i,j))
+            Sij(il, i, j) = amin1(1.0, Sij(il,i,j))
+          END IF ! new
+        END DO
+      END DO
+    ELSE  ! (ok_entrain)
+      DO il = 1,ncum
+        nent(il,i) = 0
+      ENDDO
+    ENDIF ! (ok_entrain)
 
 !jygdebug<
@@ -243,9 +251,12 @@
         uent(il, i, i) = unk(il)
         vent(il, i, i) = vnk(il)
+        IF (fl_cor_ebil .GE. 2) THEN
+          hent(il, i, i) = hp(il,i)
+        ENDIF
         elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
         Sij(il, i, i) = 0.0
       END IF
     END DO
-  END DO
+  END DO ! i = minorig + 1, nl
 
 !jyg!  DO j = 1, ntra
@@ -652,5 +663,10 @@
         vent(il, i, i) = vnk(il)
         elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
-        Sij(il, i, i) = 0.0
+        IF (fl_cor_ebil .GE. 2) THEN
+          hent(il, i, i) = hp(il,i)
+          Sigij(il, i, i) = 0.0
+        ELSE
+          Sij(il, i, i) = 0.0
+        ENDIF
       END IF
     END DO ! il
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cv3param.h
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cv3param.h	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cv3param.h	(revision 2925)
@@ -7,5 +7,7 @@
 !------------------------------------------------------------
 
+      logical ok_homo_tend
       logical ok_optim_yield
+      logical ok_entrain
       logical ok_convstop
       logical ok_intermittent
@@ -41,5 +43,7 @@
                       ,noff, minorig, nl, nlp, nlm  &
                       ,ok_convstop, ok_intermittent &
-                      ,ok_optim_yield
+                      ,ok_optim_yield &
+                      ,ok_entrain &
+                      ,ok_homo_tend
 !$OMP THREADPRIVATE(/cv3param/)
 
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cva_driver.F90
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cva_driver.F90	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/cva_driver.F90	(revision 2925)
@@ -40,4 +40,5 @@
 
   USE print_control_mod, ONLY: prt_level, lunout
+  USE add_phys_tend_mod, ONLY: fl_cor_ebil
   IMPLICIT NONE
 
@@ -730,5 +731,5 @@
              PRINT *, 'cva_driver -> cv3_feed'
     CALL cv3_feed(len, nd, ok_conserv_q, &                 ! nd->na
-                  t1, q1, u1, v1, p1, ph1, hm1, gz1, & 
+                  t1, q1, u1, v1, p1, ph1, h1, gz1, & 
                   p1feed1, p2feed1, wght1, &
                   wghti1, tnk1, thnk1, qnk1, qsnk1, unk1, vnk1, &
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/dyn1d/1DUTILS.h
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/dyn1d/1DUTILS.h	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/dyn1d/1DUTILS.h	(revision 2925)
@@ -482,4 +482,18 @@
        forc_omega =0
        CALL getin('forc_omega',forc_omega)
+
+!Config  Key  = forc_u
+!Config  Desc = forcage ou non par u
+!Config  Def  = false
+!Config  Help = forcage ou non par u
+       forc_u =0
+       CALL getin('forc_u',forc_u)
+
+!Config  Key  = forc_v
+!Config  Desc = forcage ou non par v
+!Config  Def  = false
+!Config  Help = forcage ou non par v
+       forc_v =0
+       CALL getin('forc_v',forc_v)
 
 !Config  Key  = forc_w
@@ -653,12 +667,9 @@
 
       modname = 'dyn1deta0 : '
-!!      nmq(1)="vap"
-!!      nmq(2)="cond"
-!!      do iq=3,nqtot
-!!        write(nmq(iq),'("tra",i1)') iq-2
-!!      enddo
-      DO iq = 1,nqtot
-        nmq(iq) = trim(tname(iq))
-      ENDDO
+      nmq(1)="vap"
+      nmq(2)="cond"
+      do iq=3,nqtot
+        write(nmq(iq),'("tra",i1)') iq-2
+      enddo
       print*,'in dyn1deta0 ',fichnom,klon,klev,nqtot
       CALL open_startphy(fichnom)
@@ -807,11 +818,8 @@
       CALL open_restartphy(fichnom)
       print*,'redm1 ',fichnom,klon,klev,nqtot
-!!      nmq(1)="vap"
-!!      nmq(2)="cond"
-!!      nmq(3)="tra1"
-!!      nmq(4)="tra2"
-      DO iq = 1,nqtot
-        nmq(iq) = trim(tname(iq))
-      ENDDO
+      nmq(1)="vap"
+      nmq(2)="cond"
+      nmq(3)="tra1"
+      nmq(4)="tra2"
 
       modname = 'dyn1dredem'
@@ -1396,5 +1404,4 @@
       cor(:) = rkappa*temp*(1.+q(:,1)*rv/rd)/(play*(1.+q(:,1)))
 
-
       do k=2,llm-1
 
@@ -2789,5 +2796,4 @@
          hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
          vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
-         dtrad_mod_cas(l)= dtrad_prof_cas(k2) - frac*(dtrad_prof_cas(k2)-dtrad_prof_cas(k1))
      
          else !play>plev_prof_cas(1)
@@ -2816,5 +2822,4 @@
          hq_mod_cas(l)= frac1*hq_prof_cas(k1) - frac2*hq_prof_cas(k2)
          vq_mod_cas(l)= frac1*vq_prof_cas(k1) - frac2*vq_prof_cas(k2)
-         dtrad_mod_cas(l)= frac1*dtrad_prof_cas(k1) - frac2*dtrad_prof_cas(k2)
 
          endif ! play.le.plev_prof_cas(1)
@@ -2845,5 +2850,4 @@
          hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact                            !jyg
          vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact                            !jyg
-         dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact                      !jyg
  
         endif ! play
@@ -5171,5 +5175,4 @@
          hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
          vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
-         dtrad_mod_cas(l)= dtrad_prof_cas(k2) - frac*(dtrad_prof_cas(k2)-dtrad_prof_cas(k1))
      
          else !play>plev_prof_cas(1)
@@ -5208,5 +5211,4 @@
          hq_mod_cas(l)= frac1*hq_prof_cas(k1) - frac2*hq_prof_cas(k2)
          vq_mod_cas(l)= frac1*vq_prof_cas(k1) - frac2*vq_prof_cas(k2)
-         dtrad_mod_cas(l)= frac1*dtrad_prof_cas(k1) - frac2*dtrad_prof_cas(k2)
 
          endif ! play.le.plev_prof_cas(1)
@@ -5245,5 +5247,4 @@
          hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact                     !jyg
          vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact                     !jyg
-         dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact               !jyg
  
         endif ! play
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/dyn1d/1D_interp_cases.h
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/dyn1d/1D_interp_cases.h	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/dyn1d/1D_interp_cases.h	(revision 2925)
@@ -37,5 +37,5 @@
 
        alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
-       d_th_adv(l) = ht_gcssold(l)
+       d_t_adv(l) = ht_gcssold(l)
        d_q_adv(l,1) = hq_gcssold(l)
        dt_cooling(l) = 0.0
@@ -84,5 +84,5 @@
 
        alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
-       d_th_adv(l) = alpha*omega(l)/rcpd-(ht_mod(l)+vt_mod(l))
+       d_t_adv(l) = alpha*omega(l)/rcpd-(ht_mod(l)+vt_mod(l))
        d_q_adv(l,1) = -(hq_mod(l)+vq_mod(l))
        dt_cooling(l) = 0.0
@@ -183,5 +183,5 @@
 
        alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
-       d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-d_t_dyn_z(l)
+       d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-d_t_dyn_z(l)
        d_q_adv(l,1) = hq_mod(l)-d_q_dyn_z(l)
        d_u_adv(l) = hu_mod(l)-d_u_dyn_z(l)
@@ -224,5 +224,5 @@
        ug(l)= ug_mod(l)
        vg(l)= vg_mod(l)
-       d_th_adv(l)=ht_mod(l)
+       d_t_adv(l)=ht_mod(l)
        d_q_adv(l,1)=hq_mod(l)
       enddo
@@ -315,10 +315,10 @@
 !calcul de l'advection totale
         if (cptadvw) then
-        d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-d_t_dyn_z(l)
+        d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-d_t_dyn_z(l)
 !        print*,'temp vert adv',l,ht_mod(l),vt_mod(l),-d_t_dyn_z(l)
         d_q_adv(l,1) = hq_mod(l)-d_q_dyn_z(l)
 !        print*,'q vert adv',l,hq_mod(l),vq_mod(l),-d_q_dyn_z(l)
         else
-        d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod(l)+vt_mod(l))
+        d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod(l)+vt_mod(l))
         d_q_adv(l,1) = (hq_mod(l)+vq_mod(l))
         endif
@@ -392,9 +392,9 @@
        alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
 !calcul de l'advection totale
-!        d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-omega(l)*d_t_z(l)
+!        d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-omega(l)*d_t_z(l)
 !attention: on impose dth
-        d_th_adv(l) = alpha*omega(l)/rcpd+                                  &
+        d_t_adv(l) = alpha*omega(l)/rcpd+                                  &
      &         ht_mod(l)*(play(l)/pzero)**rkappa-omega(l)*d_t_z(l)
-!        d_th_adv(l) = 0.
+!        d_t_adv(l) = 0.
 !        print*,'temp vert adv',l,ht_mod(l),vt_mod(l),-d_t_dyn_z(l)
         d_q_adv(l,1) = hq_mod(l)-omega(l)*d_q_z(l)
@@ -417,10 +417,9 @@
 !---------------------------------------------------------------------
       if (forcing_rico) then
-!       call lstendH(llm,omega,dt_dyn,dq_dyn,du_dyn, dv_dyn,
-!     :  q,temp,u,v,play)
+!      call lstendH(llm,omega,dt_dyn,dq_dyn,du_dyn, dv_dyn,q,temp,u,v,play)
        call lstendH(llm,nqtot,omega,dt_dyn,dq_dyn,q,temp,u,v,play)
 
         do l=1,llm
-       d_th_adv(l) =  (dth_rico(l) +  dt_dyn(l))
+       d_t_adv(l) =  (dth_rico(l) +  dt_dyn(l))
        d_q_adv(l,1) = (dqh_rico(l) +  dq_dyn(l,1))
        d_q_adv(l,2) = 0.
@@ -458,27 +457,27 @@
        vg(l)= v_mod(l)
        IF((phi(l)/RG).LT.1000) THEN
-         d_th_adv(l) = (adv_theta_prof + rad_theta_prof)/3600.
+         d_t_adv(l) = (adv_theta_prof + rad_theta_prof)/3600.
          d_q_adv(l,1) = adv_qt_prof/1000./3600.
          d_q_adv(l,2) = 0.0
 !        print *,'INF1000: phi dth dq1 dq2',
-!    :  phi(l)/RG,d_th_adv(l),d_q_adv(l,1),d_q_adv(l,2)
+!    :  phi(l)/RG,d_t_adv(l),d_q_adv(l,1),d_q_adv(l,2)
        ELSEIF ((phi(l)/RG).GE.1000.AND.(phi(l)/RG).lt.3000) THEN
          fact=((phi(l)/RG)-1000.)/2000.
          fact=1-fact
-         d_th_adv(l) = (adv_theta_prof + rad_theta_prof)*fact/3600.
+         d_t_adv(l) = (adv_theta_prof + rad_theta_prof)*fact/3600.
          d_q_adv(l,1) = adv_qt_prof*fact/1000./3600.
          d_q_adv(l,2) = 0.0
 !        print *,'SUP1000: phi fact dth dq1 dq2',
-!    :  phi(l)/RG,fact,d_th_adv(l),d_q_adv(l,1),d_q_adv(l,2)
+!    :  phi(l)/RG,fact,d_t_adv(l),d_q_adv(l,1),d_q_adv(l,2)
        ELSE
-         d_th_adv(l) = 0.0
+         d_t_adv(l) = 0.0
          d_q_adv(l,1) = 0.0
          d_q_adv(l,2) = 0.0
 !        print *,'SUP3000: phi dth dq1 dq2',
-!    :  phi(l)/RG,d_th_adv(l),d_q_adv(l,1),d_q_adv(l,2)
+!    :  phi(l)/RG,d_t_adv(l),d_q_adv(l,1),d_q_adv(l,2)
        ENDIF
       dt_cooling(l) = 0.0  
 !     print *,'Interp armcu: phi dth dq1 dq2',
-!    :  l,phi(l),d_th_adv(l),d_q_adv(l,1),d_q_adv(l,2)
+!    :  l,phi(l),d_t_adv(l),d_q_adv(l,1),d_q_adv(l,2)
       enddo
       endif ! forcing_armcu
@@ -554,9 +553,9 @@
        alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
 !
-!      d_th_adv(l) = 0.0
+!      d_t_adv(l) = 0.0
 !      d_q_adv(l,1) = 0.0
 !CR:test advection=0
 !calcul de l'advection verticale
-        d_th_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
+        d_t_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
 !        print*,'temp adv',l,-d_t_dyn_z(l)
         d_q_adv(l,1) = -d_q_dyn_z(l)
@@ -637,9 +636,9 @@
        alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
 !
-!      d_th_adv(l) = 0.0
+!      d_t_adv(l) = 0.0
 !      d_q_adv(l,1) = 0.0
 !CR:test advection=0
 !calcul de l'advection verticale
-        d_th_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
+        d_t_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
 !        print*,'temp adv',l,-d_t_dyn_z(l)
         d_q_adv(l,1) = -d_q_dyn_z(l)
@@ -715,5 +714,7 @@
 
 !Calcul de l advection verticale
+
       d_t_dyn_z(:)=w_mod_cas(:)*d_t_z(:)
+
       d_q_dyn_z(:)=w_mod_cas(:)*d_q_z(:)
       d_u_dyn_z(:)=w_mod_cas(:)*d_u_z(:)
@@ -764,5 +765,5 @@
 
       do l = 1, llm
-       omega(l) = w_mod_cas(l)
+       omega(l) = w_mod_cas(l)  ! juste car w_mod_cas en Pa/s (MPL 20170310)
        omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
        alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
@@ -782,9 +783,9 @@
 
         if ((tend_t.eq.1).and.(tend_w.eq.0)) then
-!           d_th_adv(l)=alpha*omega(l)/rcpd+dt_mod_cas(l)
-           d_th_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l)
+!           d_t_adv(l)=alpha*omega(l)/rcpd+dt_mod_cas(l)
+           d_t_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l)
         else if ((tend_t.eq.1).and.(tend_w.eq.1)) then
-!           d_th_adv(l)=alpha*omega(l)/rcpd+ht_mod_cas(l)-d_t_dyn_z(l)
-           d_th_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l)
+!           d_t_adv(l)=alpha*omega(l)/rcpd+ht_mod_cas(l)-d_t_dyn_z(l)
+           d_t_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l)
         endif
 
@@ -816,5 +817,4 @@
       ENDIF
       endif ! forcing_case
-
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -877,24 +877,44 @@
       d_th_z(:)=0.
       d_q_z(:)=0.
+      d_u_z(:)=0.
+      d_v_z(:)=0.
       d_t_dyn_z(:)=0.
       d_th_dyn_z(:)=0.
       d_q_dyn_z(:)=0.
+      d_u_dyn_z(:)=0.
+      d_v_dyn_z(:)=0.
       DO l=2,llm-1
        d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
        d_th_z(l)=(teta(l+1)-teta(l-1))/(play(l+1)-play(l-1))
        d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
+       d_u_z(l)=(u(l+1)-u(l-1))/(play(l+1)-play(l-1))
+       d_v_z(l)=(v(l+1)-v(l-1))/(play(l+1)-play(l-1))
       ENDDO
       d_t_z(1)=d_t_z(2)
       d_th_z(1)=d_th_z(2)
       d_q_z(1)=d_q_z(2)
+      d_u_z(1)=d_u_z(2)
+      d_v_z(1)=d_v_z(2)
       d_t_z(llm)=d_t_z(llm-1)
       d_th_z(llm)=d_th_z(llm-1)
       d_q_z(llm)=d_q_z(llm-1)
+      d_u_z(llm)=d_u_z(llm-1)
+      d_v_z(llm)=d_v_z(llm-1)
 
 !Calcul de l advection verticale
-      d_t_dyn_z(:)=w_mod_cas(:)*d_t_z(:)
-      d_th_dyn_z(:)=w_mod_cas(:)*d_th_z(:)
-      d_q_dyn_z(:)=w_mod_cas(:)*d_q_z(:)
-
+! Modif w_mod_cas -> omega_mod_cas (MM+MPL 20170310)
+      d_t_dyn_z(:)=omega_mod_cas(:)*d_t_z(:)
+      d_th_dyn_z(:)=omega_mod_cas(:)*d_th_z(:)
+      d_q_dyn_z(:)=omega_mod_cas(:)*d_q_z(:)
+      d_u_dyn_z(:)=omega_mod_cas(:)*d_u_z(:)
+      d_v_dyn_z(:)=omega_mod_cas(:)*d_v_z(:)
+
+!geostrophic wind 
+      if (forc_geo.eq.1) then
+        do l=1,llm
+        ug(l) = ug_mod_cas(l) 
+        vg(l) = vg_mod_cas(l) 
+        enddo
+      endif
 !wind nudging 
       if (nudging_u.gt.0.) then
@@ -902,8 +922,8 @@
            u(l)=u(l)+timestep*(u_mod_cas(l)-u(l))/(nudge_u)
         enddo
-      else
-        do l=1,llm
-        ug(l) = u_mod_cas(l) 
-        enddo
+!     else
+!       do l=1,llm
+!          u(l) = u_mod_cas(l) 
+!       enddo
       endif
 
@@ -912,8 +932,8 @@
            v(l)=v(l)+timestep*(v_mod_cas(l)-v(l))/(nudge_v)
         enddo
-      else
-        do l=1,llm
-        vg(l) = v_mod_cas(l) 
-        enddo
+!     else
+!       do l=1,llm
+!          v(l) = v_mod_cas(l) 
+!       enddo
       endif
 
@@ -922,8 +942,8 @@
            w(l)=w(l)+timestep*(w_mod_cas(l)-w(l))/(nudge_w)
         enddo
-      else
-        do l=1,llm
-        w(l) = w_mod_cas(l) 
-        enddo
+ !    else
+ !      do l=1,llm
+ !         w(l) = w_mod_cas(l) 
+ !      enddo
       endif
 
@@ -941,43 +961,46 @@
 
       do l = 1, llm
-       omega(l) = w_mod_cas(l)
+! Modif w_mod_cas -> omega_mod_cas (MM+MPL 20170309)
+       omega(l) = omega_mod_cas(l)
        omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
        alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
 
-!calcul advection
-!       if ((tend_u.eq.1).and.(tend_w.eq.0)) then
-!          d_u_adv(l)=du_mod_cas(l)
-!       else if ((tend_u.eq.1).and.(tend_w.eq.1)) then
-!          d_u_adv(l)=hu_mod_cas(l)-d_u_dyn_z(l)
+!calcul advections
+        if ((forc_u.eq.1).and.(forc_w.eq.0)) then
+           d_u_adv(l)=du_mod_cas(l)
+        else if ((forc_u.eq.1).and.(forc_w.eq.1)) then
+           d_u_adv(l)=hu_mod_cas(l)-d_u_dyn_z(l)
+        endif
+
+        if ((forc_v.eq.1).and.(forc_w.eq.0)) then
+           d_v_adv(l)=dv_mod_cas(l)
+        else if ((forc_v.eq.1).and.(forc_w.eq.1)) then
+           d_v_adv(l)=hv_mod_cas(l)-d_v_dyn_z(l)
+        endif
+
+! Puisque dth a ete converti en dt, on traite de la meme facon 
+! les flags tadv et thadv
+        if ((tadv.eq.1.or.thadv.eq.1) .and. (forc_w.eq.0)) then
+!          d_t_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l)
+           d_t_adv(l)=alpha*omega(l)/rcpd+dt_mod_cas(l)
+        else if ((tadv.eq.1.or.thadv.eq.1) .and. (forc_w.eq.1)) then
+!          d_t_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l)
+           d_t_adv(l)=alpha*omega(l)/rcpd+ht_mod_cas(l)-d_t_dyn_z(l)
+        endif
+
+!       if ((thadv.eq.1) .and. (forc_w.eq.0)) then
+!          d_t_adv(l)=alpha*omega(l)/rcpd-dth_mod_cas(l)
+!          d_t_adv(l)=alpha*omega(l)/rcpd+dth_mod_cas(l)
+!       else if ((thadv.eq.1) .and. (forc_w.eq.1)) then
+!          d_t_adv(l)=alpha*omega(l)/rcpd-hth_mod_cas(l)-d_t_dyn_z(l)
+!          d_t_adv(l)=alpha*omega(l)/rcpd+hth_mod_cas(l)-d_t_dyn_z(l)
 !       endif
-!
-!       if ((tend_v.eq.1).and.(tend_w.eq.0)) then
-!          d_v_adv(l)=dv_mod_cas(l)
-!       else if ((tend_v.eq.1).and.(tend_w.eq.1)) then
-!          d_v_adv(l)=hv_mod_cas(l)-d_v_dyn_z(l)
-!       endif
-!
-!-----------------------------------------------------
-        if (tadv.eq.1 .or. tadvh.eq.1) then
-           d_t_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l)
-        else if (tadvv.eq.1) then
-! ATTENTION d_t_dyn_z pas calcule (voir twpice)
-           d_t_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l)
-        endif
-        print *,'interp_case d_t_dyn_z=',d_t_dyn_z(l),d_q_dyn_z(l)
-
-! Verifier le signe !!
-        if (thadv.eq.1 .or. thadvh.eq.1) then
-           d_th_adv(l)=dth_mod_cas(l)
-           print *,'dthadv=',d_th_adv(l)*86400.
-        else if (thadvv.eq.1) then
-           d_th_adv(l)=hth_mod_cas(l)-d_th_dyn_z(l)
-        endif
- 
-! Verifier le signe !!
-        if ((qadv.eq.1).and.(forc_w.eq.0)) then
+
+        if ((qadv.eq.1) .and. (forc_w.eq.0)) then
            d_q_adv(l,1)=dq_mod_cas(l)
-        else if ((qadvh.eq.1).and.(forc_w.eq.1)) then
+!          d_q_adv(l,1)=-1*dq_mod_cas(l)
+        else if ((qadv.eq.1) .and. (forc_w.eq.1)) then
            d_q_adv(l,1)=hq_mod_cas(l)-d_q_dyn_z(l)
+!          d_q_adv(l,1)=-1*hq_mod_cas(l)-d_q_dyn_z(l)
         endif
          
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/dyn1d/1D_nudge_sandu_astex.h
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/dyn1d/1D_nudge_sandu_astex.h	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/dyn1d/1D_nudge_sandu_astex.h	(revision 2925)
@@ -18,5 +18,5 @@
 !     print *,'l dq1 relax dqadv',l,dq(l,1),relax_q(l,1),d_q_adv(l,1)
 !     print *,'l dq2 relax dqadv',l,dq(l,2),relax_q(l,2),d_q_adv(l,2)
-!     print *,'l dt  relax dtadv',l,dt_phys(l),relax_thl(l),d_th_adv(l)
+!     print *,'l dt  relax dtadv',l,dt_phys(l),relax_thl(l),d_t_adv(l)
       enddo
         u(1:mxcalc)=u(1:mxcalc) + timestep*(                                &          
@@ -32,3 +32,4 @@
      &      dq(1:mxcalc,2) - relax_q(1:mxcalc,2)+d_q_adv(1:mxcalc,2))
         temp(1:mxcalc)=temp(1:mxcalc)+timestep*(                            &
-     &      dt_phys(1:mxcalc)-relax_thl(1:mxcalc)+d_th_adv(1:mxcalc))
+     &      dt_phys(1:mxcalc)-relax_thl(1:mxcalc)+d_t_adv(1:mxcalc))
+
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/dyn1d/1D_read_forc_cases.h
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/dyn1d/1D_read_forc_cases.h	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/dyn1d/1D_read_forc_cases.h	(revision 2925)
@@ -76,4 +76,5 @@
         dt_cooling(l)=thlpcar(kmax)-frac*(thlpcar(kmax)-thlpcar(kmax-1))
         do k=2,kmax
+          print *,'k l height(k) height(k-1) zlay(l) frac=',k,l,height(k),height(k-1),zlay(l),frac
           frac = (height(k)-zlay(l))/(height(k)-height(k-1))
           if(l==1) print*,'k, height, tttprof',k,height(k),tttprof(k)
@@ -227,5 +228,5 @@
 !?       omega2(l)=-rho(l)*omega(l)
        alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
-       d_th_adv(l) = alpha*omega(l)/rcpd-(ht_mod(l)+vt_mod(l))
+       d_t_adv(l) = alpha*omega(l)/rcpd-(ht_mod(l)+vt_mod(l))
        d_q_adv(l,1) = -(hq_mod(l)+vq_mod(l))
        d_q_adv(l,2) = 0.0
@@ -280,5 +281,5 @@
 !on applique le forcage total au premier pas de temps
 !attention: signe different de toga
-       d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod(l)+vt_mod(l))
+       d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod(l)+vt_mod(l))
        d_q_adv(l,1) = (hq_mod(l)+vq_mod(l))
        d_q_adv(l,2) = 0.0
@@ -341,7 +342,7 @@
 !on applique le forcage total au premier pas de temps
 !attention: signe different de toga
-       d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
+       d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
 !forcage en th
-!       d_th_adv(l) = ht_mod(l)
+!       d_t_adv(l) = ht_mod(l)
        d_q_adv(l,1) = hq_mod(l)
        d_q_adv(l,2) = 0.0
@@ -442,7 +443,7 @@
 !on applique le forcage total au premier pas de temps
 !attention: signe different de toga
-       d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
+       d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
 !forcage en th
-!       d_th_adv(l) = ht_mod(l)
+!       d_t_adv(l) = ht_mod(l)
        d_q_adv(l,1) = hq_mod(l)
        d_q_adv(l,2) = 0.0
@@ -557,7 +558,7 @@
 !on applique le forcage total au premier pas de temps
 !attention: signe different de toga
-!      d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
+!      d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
 !forcage en th
-       d_th_adv(l) = ht_mod(l)
+       d_t_adv(l) = ht_mod(l)
        d_q_adv(l,1) = hq_mod(l)
        d_q_adv(l,2) = 0.0
@@ -657,16 +658,16 @@
 ! Advective forcings are given in K or g/kg ... per HOUR
 !      IF(height(l).LT.1000) THEN
-!        d_th_adv(l) = (adv_theta_prof + rad_theta_prof)/3600.
+!        d_t_adv(l) = (adv_theta_prof + rad_theta_prof)/3600.
 !        d_q_adv(l,1) = adv_qt_prof/1000./3600.
 !        d_q_adv(l,2) = 0.0
 !      ELSEIF (height(l).GE.1000.AND.height(l).LT.3000) THEN
-!        d_th_adv(l) = (adv_theta_prof + rad_theta_prof)*
+!        d_t_adv(l) = (adv_theta_prof + rad_theta_prof)*
 !    :               (1-(height(l)-1000.)/2000.)
-!        d_th_adv(l) = d_th_adv(l)/3600.
+!        d_t_adv(l) = d_t_adv(l)/3600.
 !        d_q_adv(l,1) = adv_qt_prof*(1-(height(l)-1000.)/2000.)
 !        d_q_adv(l,1) = d_q_adv(l,1)/1000./3600.
 !        d_q_adv(l,2) = 0.0
 !      ELSE
-!        d_th_adv(l) = 0.0
+!        d_t_adv(l) = 0.0
 !        d_q_adv(l,1) = 0.0
 !        d_q_adv(l,2) = 0.0
@@ -751,7 +752,7 @@
 !?       omega2(l)=-rho(l)*omega(l)
        alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
-!      d_th_adv(l) = alpha*omega(l)/rcpd+vt_mod(l)
+!      d_t_adv(l) = alpha*omega(l)/rcpd+vt_mod(l)
 !      d_q_adv(l,1) = vq_mod(l)
-       d_th_adv(l) = alpha*omega(l)/rcpd
+       d_t_adv(l) = alpha*omega(l)/rcpd
        d_q_adv(l,1) = 0.0
        d_q_adv(l,2) = 0.0
@@ -827,7 +828,7 @@
 !      omega2(l)=-rho(l)*omega(l)
        alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
-!      d_th_adv(l) = alpha*omega(l)/rcpd+vt_mod(l)
+!      d_t_adv(l) = alpha*omega(l)/rcpd+vt_mod(l)
 !      d_q_adv(l,1) = vq_mod(l)
-       d_th_adv(l) = alpha*omega(l)/rcpd
+       d_t_adv(l) = alpha*omega(l)/rcpd
        d_q_adv(l,1) = 0.0
        d_q_adv(l,2) = 0.0
@@ -890,9 +891,10 @@
 !on applique le forcage total au premier pas de temps
 !attention: signe different de toga
-       d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l))
+       d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l))
        d_q_adv(l,1) = (hq_mod_cas(l)+vq_mod_cas(l))
        d_q_adv(l,2) = 0.0
        d_u_adv(l) = (hu_mod_cas(l)+vu_mod_cas(l))
-       d_u_adv(l) = (hv_mod_cas(l)+vv_mod_cas(l))
+! correction bug d_u -> d_v (MM+MPL 20170310)
+       d_v_adv(l) = (hv_mod_cas(l)+vv_mod_cas(l))
       enddo     
 
@@ -971,8 +973,9 @@
        q(l,2) = ql_mod_cas(l)
        u(l) = u_mod_cas(l)
-       ug(l)= u_mod_cas(l)
+       ug(l)= ug_mod_cas(l)
        v(l) = v_mod_cas(l)
-       vg(l)= v_mod_cas(l)
-       omega(l) = w_mod_cas(l)
+       vg(l)= vg_mod_cas(l)
+! Modif w_mod_cas -> omega_mod_cas (MM+MPL 20170309)
+       omega(l) = omega_mod_cas(l)
        omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
 
@@ -980,5 +983,5 @@
 !on applique le forcage total au premier pas de temps
 !attention: signe different de toga
-       d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l))
+       d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l))
        d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l))
 !      d_q_adv(l,1) = (hq_mod_cas(l)+vq_mod_cas(l))
@@ -987,6 +990,7 @@
 !      d_u_adv(l) = (hu_mod_cas(l)+vu_mod_cas(l))
        d_u_adv(l) = du_mod_cas(l)
-!      d_u_adv(l) = (hv_mod_cas(l)+vv_mod_cas(l))
-       d_u_adv(l) = dv_mod_cas(l)
+!      d_v_adv(l) = (hv_mod_cas(l)+vv_mod_cas(l))
+! correction bug d_u -> d_v (MM+MPL 20170310)
+       d_v_adv(l) = dv_mod_cas(l)
       enddo     
 
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/dyn1d/compar1d.h
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/dyn1d/compar1d.h	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/dyn1d/compar1d.h	(revision 2925)
@@ -41,5 +41,5 @@
 
       integer :: tadv, tadvv, tadvh, qadv, qadvv, qadvh, thadv, thadvv, thadvh, trad
-      integer :: forc_omega, forc_w, forc_geo, forc_ustar
+      integer :: forc_omega, forc_u, forc_v, forc_w, forc_geo, forc_ustar
       real    :: nudging_u, nudging_v, nudging_w, nudging_t, nudging_q
       common/com_par1d/                                                 &
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/dyn1d/lmdz1d.F90
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/dyn1d/lmdz1d.F90	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/dyn1d/lmdz1d.F90	(revision 2925)
@@ -194,5 +194,5 @@
       real :: du_phys(llm),dv_phys(llm),dt_phys(llm)
       real :: dt_dyn(llm)
-      real :: dt_cooling(llm),d_t_adv(llm),d_th_adv(llm),d_t_nudge(llm)
+      real :: dt_cooling(llm),d_t_adv(llm),d_t_nudge(llm)
       real :: d_u_nudge(llm),d_v_nudge(llm)
       real :: du_adv(llm),dv_adv(llm)
@@ -206,5 +206,5 @@
       REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_adv
       REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_nudge
-!      REAL, ALLOCATABLE, DIMENSION(:):: d_th_adv
+!      REAL, ALLOCATABLE, DIMENSION(:):: d_t_adv
 
 !---------------------------------------------------------------------
@@ -271,5 +271,5 @@
       dt_dyn(:)=0.
       dt_cooling(:)=0.
-      d_th_adv(:)=0.
+      d_t_adv(:)=0.
       d_t_nudge(:)=0.
       d_u_nudge(:)=0.
@@ -404,4 +404,11 @@
        heure_ini_cas=0.
        pdt_cas=1800.         ! forcing frequency
+      elseif (forcing_type .eq.105) THEN ! bomex starts 16-12-2004 0h
+       forcing_case2 = .true.
+       year_ini_cas=1969
+       mth_ini_cas=6
+       day_deb=24
+       heure_ini_cas=0.
+       pdt_cas=1800.         ! forcing frequency
       elseif (forcing_type .eq.40) THEN
        forcing_GCSSold = .true.
@@ -573,5 +580,5 @@
       allocate(d_q_adv(llm,nqtot)) 
       allocate(d_q_nudge(llm,nqtot)) 
-!      allocate(d_th_adv(llm)) 
+!      allocate(d_t_adv(llm)) 
 
       q(:,:) = 0.
@@ -1057,5 +1064,5 @@
          fcoriolis=0.0
          dt_cooling=0.0
-         d_th_adv=0.0
+         d_t_adv=0.0
          d_q_adv=0.0
        endif
@@ -1143,6 +1150,6 @@
         if (prt_level.ge.3) then
           print *,                                                          &
-     &    'physiq-> temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1) ',         &
-     &              temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1)
+     &    'physiq-> temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1) ',         &
+     &              temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1)
            print* ,'dv_phys=',dv_phys
            print* ,'dv_age=',dv_age
@@ -1155,5 +1162,5 @@
         temp(1:mxcalc)=temp(1:mxcalc)+timestep*(                            &
      &              dt_phys(1:mxcalc)                                       &
-     &             +d_th_adv(1:mxcalc)                                      &
+     &             +d_t_adv(1:mxcalc)                                      &
      &             +d_t_nudge(1:mxcalc)                                      &
      &             +dt_cooling(1:mxcalc))  ! Taux de chauffage ou refroid.
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/ener_conserv.F90
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/ener_conserv.F90	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/ener_conserv.F90	(revision 2925)
@@ -27,4 +27,6 @@
 USE phys_state_var_mod, ONLY : du_gwd_front,du_gwd_rando
 USE phys_output_var_mod, ONLY : bils_ec,bils_ech,bils_tke,bils_kinetic,bils_enthalp,bils_latent,bils_diss
+USE add_phys_tend_mod, ONLY : fl_cor_ebil 
+
 
 IMPLICIT none
@@ -62,9 +64,13 @@
    DO k = 1, klev
    DO i = 1, klon
-      ZRCPD = RCPD*(1.0+RVTMP2*(pqn(i,k)+pqln(i,k)+pqsn(i,k)))
-      d_t_ec(i,k)=0.5/ZRCPD &
- &      *(puo(i,k)**2+pvo(i,k)**2-pun(i,k)**2-pvn(i,k)**2)
-      ENDDO
-      ENDDO
+     IF (fl_cor_ebil .GT. 0) then
+       ZRCPD = RCPD*(1.0+RVTMP2*(pqn(i,k)+pqln(i,k)+pqsn(i,k)))
+     ELSE
+       ZRCPD = RCPD*(1.0+RVTMP2*pqn(i,k))
+     ENDIF
+     d_t_ec(i,k)=0.5/ZRCPD &
+ &     *(puo(i,k)**2+pvo(i,k)**2-pun(i,k)**2-pvn(i,k)**2)
+   ENDDO
+   ENDDO
 !-jld ec_conser
 
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/fisrtilp.F90
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/fisrtilp.F90	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/fisrtilp.F90	(revision 2925)
@@ -103,11 +103,12 @@
   ! Options du programme:
   !
-  REAL seuil_neb ! un nuage existe vraiment au-dela
-  PARAMETER (seuil_neb=0.001)
+  REAL, SAVE :: seuil_neb=0.001 ! un nuage existe vraiment au-dela
+  !$OMP THREADPRIVATE(seuil_neb)
+
 
   INTEGER ninter ! sous-intervals pour la precipitation
   PARAMETER (ninter=5)
-  LOGICAL evap_prec ! evaporation de la pluie
-  PARAMETER (evap_prec=.TRUE.)
+  INTEGER,SAVE :: iflag_evap_prec=1 ! evaporation de la pluie
+  !$OMP THREADPRIVATE(iflag_evap_prec)
   !
   LOGICAL cpartiel ! condensation partielle
@@ -124,4 +125,7 @@
   !
   INTEGER i, k, n, kk
+  INTEGER,save::itap=0
+  !$OMP THREADPRIVATE(itap)
+
   REAL qsl, qsi
   real zct      ,zcl
@@ -154,5 +158,5 @@
   REAL t_glace_min_old
   REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
-  REAL zchau      ,zfroi      ,zfice(klon),zneb(klon)
+  REAL zchau      ,zfroi      ,zfice(klon),zneb(klon),znebprecip(klon)
   REAL zmelt, zpluie, zice
   REAL dzfice(klon)
@@ -212,4 +216,8 @@
 !CR: pour iflag_ice_thermo=2, on active que la convection
 !  ice_thermo = iflag_ice_thermo .GE. 1
+
+itap=itap+1
+znebprecip(:)=0.
+
    ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
   zdelq=0.0
@@ -218,9 +226,12 @@
   IF (appel1er) THEN
      CALL getin_p('iflag_oldbug_fisrtilp',iflag_oldbug_fisrtilp)
+     CALL getin_p('iflag_evap_prec',iflag_evap_prec)
+     CALL getin_p('seuil_neb',seuil_neb)
      write(lunout,*)' iflag_oldbug_fisrtilp =',iflag_oldbug_fisrtilp 
      !
      WRITE(lunout,*) 'fisrtilp, ninter:', ninter
-     WRITE(lunout,*) 'fisrtilp, evap_prec:', evap_prec
+     WRITE(lunout,*) 'fisrtilp, iflag_evap_prec:', iflag_evap_prec
      WRITE(lunout,*) 'fisrtilp, cpartiel:', cpartiel
+     
      IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
         WRITE(lunout,*) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
@@ -387,5 +398,5 @@
      !   - zmqc: masse de precip qui doit etre thermalisee
      !
-     IF (evap_prec) THEN
+     IF (iflag_evap_prec>=1) THEN
         DO i = 1, klon
 !          S'il y a des precipitations
@@ -480,7 +491,13 @@
 !        S'il y a des precipitations
          IF (zrfl(i)+zifl(i).GT.0.) THEN
+
+         IF (iflag_evap_prec==1) THEN
+            znebprecip(i)=zneb(i)
+         ELSE
+            znebprecip(i)=MAX(zneb(i),znebprecip(i))
+         ENDIF
      
         ! Evap max pour ne pas saturer la fraction sous le nuage
-         zqev0 = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
+         zqev0 = MAX (0.0, (zqs(i)-zq(i))*znebprecip(i) )
 
          !JAM
@@ -567,4 +584,5 @@
          zrfl(i) = zrfln(i)
          zifl(i) = zifln(i)
+!        print*,'REEVAP ',itap,k,znebprecip(1),zqev0,zqev,zqevi,zrfl(1)
 
 !CR ATTENTION: deplacement de la fonte de la glace
@@ -595,4 +613,8 @@
            end if
 
+           ELSE
+              ! Si on n'a plus de pluies, on reinitialise a 0 la farcion
+              ! sous nuageuse utilisee pour la pluie.
+              znebprecip(i)=0.
            ENDIF ! (zrfl(i)+zifl(i).GT.0.)
         ENDDO
@@ -603,5 +625,5 @@
      ! Fin evaporation de la precipitation
      ! ----------------------------------------------------------------
-     ENDIF ! (evap_prec)
+     ENDIF ! (iflag_evap_prec>=1)
      !
      ! Calculer Qs et L/Cp*dQs/dT:
@@ -702,5 +724,5 @@
                    qcloud,ctot,zpspsk,paprs,ztla,zthl, &
                    ratqs,zqs,t)
-              elseif (iflag_cloudth_vert==3) then
+              elseif (iflag_cloudth_vert>=3) then
                call cloudth_v3(klon,klev,k,ztv, &
                    zq,zqta,fraca, &
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/oasis.F90
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/oasis.F90	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/oasis.F90	(revision 2925)
@@ -80,5 +80,7 @@
 
   TYPE(FLD_CPL), DIMENSION(maxsend), SAVE, PUBLIC :: infosend   ! Information for sending coupling fields
+!$OMP THREADPRIVATE(infosend)
   TYPE(FLD_CPL), DIMENSION(maxrecv), SAVE, PUBLIC :: inforecv   ! Information for receiving coupling fields
+!$OMP THREADPRIVATE(inforecv)
   
   LOGICAL,SAVE :: cpl_current
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/physiq_mod.F90	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/physiq_mod.F90	(revision 2925)
@@ -245,5 +245,5 @@
 
     USE cmp_seri_mod
-    USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, prt_enerbil, &
+    USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
   &      fl_ebil, fl_cor_ebil
 
@@ -1682,4 +1682,5 @@
                pdtphys, &
                annee_ref, &
+               year_cur, &
                day_ref,  &
                day_ini, &
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/sisvat/sisvat.F
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/sisvat/sisvat.F	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/sisvat/sisvat.F	(revision 2925)
@@ -7070,5 +7070,5 @@
 
       DO ig = 1,knonv
-        ztherm_i(ig)   = inertie_ice
+        ztherm_i(ig)   = inertie_lic
         IF (isnoSV(ig) > 0) ztherm_i(ig)   = inertie_sno
       ENDDO
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/soil.F90
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/soil.F90	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/soil.F90	(revision 2925)
@@ -152,19 +152,33 @@
 !-----------------------------------------------------------------------
 !   Calcul de l'inertie thermique a partir de la variable rnat.
-!   on initialise a inertie_ice meme au-dessus d'un point de mer au cas 
+!   on initialise a inertie_sic meme au-dessus d'un point de mer au cas 
 !   ou le point de mer devienne point de glace au pas suivant
 !   on corrige si on a un point de terre avec ou sans glace
 !
-!-----------------------------------------------------------------------
+!   iophys can be used to write the ztherm_i variable in a phys.nc file
+!   and check the results; to do so, add "CALL iophys_ini" in physiq_mod
+!   and add knindex to the list of inputs in all the calls to soil.F90
+!   (and to soil.F90 itself !)
+!-----------------------------------------------------------------------
+
   IF (indice == is_sic) THEN
      DO ig = 1, knon
-        ztherm_i(ig)   = inertie_ice
+        ztherm_i(ig)   = inertie_sic
+     ENDDO
+     IF (iflag_sic == 0) THEN
+       DO ig = 1, knon
+         IF (snow(ig) > 0.0) ztherm_i(ig)   = inertie_sno
+       ENDDO
+!      Otherwise sea-ice keeps the same inertia, even when covered by snow
+     ENDIF
+!    CALL iophys_ecrit_index('ztherm_sic', 1, 'ztherm_sic', 'USI', &
+!      knon, knindex, ztherm_i)
+  ELSE IF (indice == is_lic) THEN
+     DO ig = 1, knon
+        ztherm_i(ig)   = inertie_lic
         IF (snow(ig) > 0.0) ztherm_i(ig)   = inertie_sno
      ENDDO
-  ELSE IF (indice == is_lic) THEN
-     DO ig = 1, knon
-        ztherm_i(ig)   = inertie_ice
-        IF (snow(ig) > 0.0) ztherm_i(ig)   = inertie_sno
-     ENDDO
+!    CALL iophys_ecrit_index('ztherm_lic', 1, 'ztherm_lic', 'USI', &
+!      knon, knindex, ztherm_i)
   ELSE IF (indice == is_ter) THEN
      DO ig = 1, knon
@@ -172,8 +186,13 @@
         IF (snow(ig) > 0.0) ztherm_i(ig)   = inertie_sno
      ENDDO
+!    CALL iophys_ecrit_index('ztherm_ter', 1, 'ztherm_ter', 'USI', &
+!      knon, knindex, ztherm_i)
   ELSE IF (indice == is_oce) THEN
      DO ig = 1, knon
-        ztherm_i(ig)   = inertie_ice
-     ENDDO
+!       This is just in case, but SST should be used by the model anyway
+        ztherm_i(ig)   = inertie_sic
+     ENDDO
+!    CALL iophys_ecrit_index('ztherm_oce', 1, 'ztherm_oce', 'USI', &
+!      knon, knindex, ztherm_i)
   ELSE
      WRITE(lunout,*) "valeur d indice non prevue", indice
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/surf_land_orchidee_mod.F90
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/surf_land_orchidee_mod.F90	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/surf_land_orchidee_mod.F90	(revision 2925)
@@ -32,5 +32,4 @@
   PUBLIC  :: surf_land_orchidee
 
-  LOGICAL, ALLOCATABLE, SAVE :: flag_omp(:)
 CONTAINS
 !
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/surf_land_orchidee_noopenmp_mod.F90
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/surf_land_orchidee_noopenmp_mod.F90	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/surf_land_orchidee_noopenmp_mod.F90	(revision 2925)
@@ -144,5 +144,7 @@
     INTEGER                                   :: error
     INTEGER, SAVE                             :: nb_fields_cpl ! number of fields for the climate-carbon coupling (between ATM and ORCHIDEE). 
+    !$OMP THREADPRIVATE(nb_fields_cpl)
     REAL, SAVE, ALLOCATABLE, DIMENSION(:,:)   :: fields_cpl    ! Fluxes for the climate-carbon coupling
+    !$OMP THREADPRIVATE(fields_cpl)
     REAL, DIMENSION(klon)                     :: swdown_vrai
     CHARACTER (len = 20)                      :: modname = 'surf_land_orchidee'
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/surf_land_orchidee_noz0h_mod.F90
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/surf_land_orchidee_noz0h_mod.F90	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/surf_land_orchidee_noz0h_mod.F90	(revision 2925)
@@ -33,5 +33,4 @@
   PUBLIC  :: surf_land_orchidee
 
-  LOGICAL, ALLOCATABLE, SAVE :: flag_omp(:)
 CONTAINS
 !
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/surface_data.F90
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/surface_data.F90	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/surface_data.F90	(revision 2925)
@@ -9,10 +9,11 @@
   
   LOGICAL, SAVE          :: ok_veget      ! true for use of vegetation model ORCHIDEE
-  CHARACTER(len=10), SAVE :: type_veget    ! orchidee/y/bucket/n/betaclim
   !$OMP THREADPRIVATE(ok_veget)
-  ! Martin
+
+  CHARACTER(len=10), SAVE :: type_veget   ! orchidee/y/bucket/n/betaclim
+  !$OMP THREADPRIVATE(type_veget)
+
   LOGICAL, SAVE          :: ok_snow       ! true for coupling to snow model SISVAT
   !$OMP THREADPRIVATE(ok_snow)
-  ! Martin
 
   CHARACTER(len=6), SAVE :: type_ocean    ! force/slab/couple
Index: /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/wake.F90
===================================================================
--- /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/wake.F90	(revision 2924)
+++ /LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/wake.F90	(revision 2925)
@@ -178,4 +178,6 @@
   LOGICAL, SAVE                                         :: flag_wk_check_trgl
   !$OMP THREADPRIVATE(flag_wk_check_trgl)
+  INTEGER, SAVE                                         :: iflag_wk_check_trgl
+  !$OMP THREADPRIVATE(iflag_wk_check_trgl)
 
   REAL                                                  :: delta_t_min
@@ -347,4 +349,8 @@
   CALL getin_p('flag_wk_check_trgl ', flag_wk_check_trgl)
   WRITE(*,*) 'flag_wk_check_trgl=', flag_wk_check_trgl
+  WRITE(*,*) 'flag_wk_check_trgl OBSOLETE. Utilisr iflag_wk_check_trgl plutot'
+  iflag_wk_check_trgl=0 ; IF (flag_wk_check_trgl) iflag_wk_check_trgl=1
+  CALL getin_p('iflag_wk_check_trgl ', iflag_wk_check_trgl)
+  WRITE(*,*) 'iflag_wk_check_trgl=', iflag_wk_check_trgl
 
   first=.false.
@@ -1793,5 +1799,5 @@
   ! Filter out bad wakes
 
-  IF (flag_wk_check_trgl) THEN
+  IF (iflag_wk_check_trgl>=1) THEN
     ! Check triangular shape of dth profile
     DO i = 1, klon
@@ -1805,5 +1811,5 @@
           wape2(i) = -1.
           !! print *,'wake, rej 1'
-        ELSE IF (abs(2.*sum_dth(i)/(hw0(i)*dthmin(i)) - 1.) > 0.5) THEN
+        ELSE IF (iflag_wk_check_trgl==1.AND.abs(2.*sum_dth(i)/(hw0(i)*dthmin(i)) - 1.) > 0.5) THEN
           wape2(i) = -1.
           !! print *,'wake, rej 2'
