Index: /LMDZ5/trunk/DefLists/field_def_lmdz.xml
===================================================================
--- /LMDZ5/trunk/DefLists/field_def_lmdz.xml	(revision 2579)
+++ /LMDZ5/trunk/DefLists/field_def_lmdz.xml	(revision 2580)
@@ -837,4 +837,34 @@
     </field_group>
 
+    <field_group id="fields_AIRS" domain_ref="dom_glo" operation="instant" freq_op="3h">
+      <field id="alt_tropo" long_name="Tropopause pressure" unit="hPa" />
+      <field id="map_prop_hc" long_name="Proportion of high clouds" unit=" "/>
+      <field id="map_prop_hist" long_name="Proportion of high ice semi-transp clouds" unit=" "/>
+      <field id="map_emis_hc" long_name="Emissivity of high clouds" unit=" "/>
+      <field id="map_iwp_hc" long_name="Ice water path of high clouds" unit="g/m2"/>
+      <field id="map_deltaz_hc" long_name="geom thickness of high clouds" unit="m"/>
+      <field id="map_pcld_hc" long_name="cloud pressure of high clouds" unit="hPa"/>
+      <field id="map_tcld_hc" long_name="cloud temperature of high clouds" unit="K"/>
+      <field id="map_emis_hist" long_name="Emissivity of high ice st clouds" unit=" "/>
+      <field id="map_iwp_hist" long_name="Ice water path of high ice st clouds" unit="g/m2"/>
+      <field id="map_deltaz_hist" long_name="geom thickness of high ice st clouds" unit="m"/>
+      <field id="map_rad_hist" long_name="ice crystals radius in high ice st clouds" unit="m"/>
+      <field id="map_emis_Cb" long_name="Emissivity of high Cb clouds" unit=" "/>
+      <field id="map_pcld_Cb" long_name="cloud pressure of high Cb clouds" unit="hPa"/>
+      <field id="map_tcld_Cb" long_name="cloud temperature of high Cb clouds" unit="K"/>
+      <field id="map_emis_Anv" long_name="Emissivity of high Anv clouds" unit=" "/>
+      <field id="map_pcld_Anv" long_name="cloud pressure of high Anv clouds" unit="hPa"/>
+      <field id="map_tcld_Anv" long_name="cloud temperature of high Anv clouds" unit="K"/>
+      <field id="map_emis_ThCi" long_name="Emissivity of high ThCi clouds" unit=" "/>
+      <field id="map_pcld_ThCi" long_name="cloud pressure of high ThCi clouds" unit="hPa"/>
+      <field id="map_tcld_ThCi" long_name="cloud temperature of high ThCi clouds" unit="K"/>
+      <field id="map_ntot" long_name="total AIRS cloud fraction" unit=" "/>
+      <field id="map_hc" long_name="high clouds AIRS cloud fraction" unit=" "/>
+      <field id="map_hist" long_name="high clouds ice st AIRS cloud fraction" unit=" "/>
+      <field id="map_Cb" long_name="high clouds Cb AIRS cloud fraction" unit=" "/>
+      <field id="map_ThCi" long_name="high clouds ThCi AIRS cloud fraction" unit=" "/>
+      <field id="map_Anv" long_name="high clouds Anv AIRS cloud fraction" unit=" "/>
+    </field_group>
+
 </field_definition>
 
Index: /LMDZ5/trunk/DefLists/file_def_histhf_lmdz.xml
===================================================================
--- /LMDZ5/trunk/DefLists/file_def_histhf_lmdz.xml	(revision 2579)
+++ /LMDZ5/trunk/DefLists/file_def_histhf_lmdz.xml	(revision 2580)
@@ -411,4 +411,34 @@
             </field_group>
 
+            <field_group operation="instant" freq_op="1ts" >
+                 <field field_ref="alt_tropo" level="10" />
+                 <field field_ref="map_prop_hc" level="10" />
+                 <field field_ref="map_prop_hist" level="10" />
+                 <field field_ref="map_emis_hc" level="10" />
+                 <field field_ref="map_iwp_hc" level="10" />
+                 <field field_ref="map_deltaz_hc" level="10" />
+                 <field field_ref="map_pcld_hc" level="10" />
+                 <field field_ref="map_tcld_hchc" level="10" />
+                 <field field_ref="map_emis_hist" level="10" />
+                 <field field_ref="map_iwp_hist" level="10" />
+                 <field field_ref="map_deltaz_hist" level="10" />
+                 <field field_ref="map_rad_hist" level="10" />
+                 <field field_ref="map_emis_Cb" level="10" />
+                 <field field_ref="map_pcld_Cb" level="10" />
+                 <field field_ref="map_tcld_Cb" level="10" />
+                 <field field_ref="map_emis_Anv" level="10" />
+                 <field field_ref="map_pcld_Anv" level="10" />
+                 <field field_ref="map_tcld_Anv" level="10" />
+                 <field field_ref="map_emis_ThCi" level="10" />
+                 <field field_ref="map_pcld_ThCi" level="10" />
+                 <field field_ref="map_tcld_ThCi" level="10" />
+                 <field field_ref="map_ntot" level="10" />
+                 <field field_ref="map_hc" level="10" />
+                 <field field_ref="map_hist" level="10" />
+                 <field field_ref="map_Cb" level="10" />
+                 <field field_ref="map_ThCi" level="10" />
+                 <field field_ref="map_Anv" level="10" />
+            </field_group>
+
             <!-- VARS 3D -->
             <field_group operation="average" freq_op="1ts" axis_ref="presnivs">
Index: /LMDZ5/trunk/libf/phylmd/clesphys.h
===================================================================
--- /LMDZ5/trunk/libf/phylmd/clesphys.h	(revision 2579)
+++ /LMDZ5/trunk/libf/phylmd/clesphys.h	(revision 2580)
@@ -78,6 +78,7 @@
        REAL ecrit_LES
        REAL freq_ISCCP, ecrit_ISCCP
-       REAL freq_COSP
+       REAL freq_COSP, freq_AIRS
        LOGICAL :: ok_cosp,ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP
+       LOGICAL :: ok_airs
        INTEGER :: ip_ebil_phy, iflag_rrtm, iflag_ice_thermo, NSW, iflag_albedo
        LOGICAL :: ok_chlorophyll
@@ -102,5 +103,5 @@
      &     , pasphys            , freq_outNMC, freq_calNMC              &
      &     , lonmin_ins, lonmax_ins, latmin_ins, latmax_ins             &
-     &     , freq_ISCCP, ecrit_ISCCP, freq_COSP                         &
+     &     , freq_ISCCP, ecrit_ISCCP, freq_COSP, freq_AIRS              &
      &     , cvl_corr                                                   &
      &     , qsol0,albsno0,evap0                                        &
@@ -121,5 +122,5 @@
      &     , lev_histins, lev_histLES, lev_histdayNMC, levout_histNMC   &
      &     , ok_histNMC                                                 &
-     &     , type_run, ok_regdyn, ok_cosp                               &
+     &     , type_run, ok_regdyn, ok_cosp, ok_airs                      &
      &     , ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP                     &
      &     , ip_ebil_phy                                                &
Index: /LMDZ5/trunk/libf/phylmd/conf_phys_m.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/conf_phys_m.F90	(revision 2579)
+++ /LMDZ5/trunk/libf/phylmd/conf_phys_m.F90	(revision 2580)
@@ -98,5 +98,5 @@
     REAL,SAVE           :: bl95_b0_omp, bl95_b1_omp
     REAL,SAVE           :: freq_ISCCP_omp, ecrit_ISCCP_omp
-    REAL,SAVE           :: freq_COSP_omp
+    REAL,SAVE           :: freq_COSP_omp, freq_AIRS_omp 
     real,SAVE           :: fact_cldcon_omp, facttemps_omp,ratqsbas_omp
     real,SAVE           :: tau_cld_cv_omp, coefw_cld_cv_omp
@@ -190,5 +190,5 @@
     REAL, SAVE :: freq_outNMC_omp(3), freq_calNMC_omp(3)
     CHARACTER*4, SAVE :: type_run_omp
-    LOGICAL,SAVE :: ok_cosp_omp
+    LOGICAL,SAVE :: ok_cosp_omp, ok_airs_omp
     LOGICAL,SAVE :: ok_mensuelCOSP_omp,ok_journeCOSP_omp,ok_hfCOSP_omp
     REAL,SAVE :: lonmin_ins_omp, lonmax_ins_omp, latmin_ins_omp, latmax_ins_omp
@@ -462,4 +462,13 @@
     call getin('freq_COSP', freq_COSP_omp)
 
+    !Config Key  = freq_AIRS
+    !Config Desc = Frequence d'appel du simulateur AIRS en secondes;
+    !              par defaut 10800, i.e. 3 heures
+    !Config Def  = 10800.
+    !Config Help = Used in ini_histdayAIRS.h
+    !
+    freq_AIRS_omp = 10800.
+    call getin('freq_AIRS', freq_AIRS_omp)
+
     !
     !Config Key  = ip_ebil_phy
@@ -1609,4 +1618,13 @@
     ok_cosp_omp = .false.
     call getin('ok_cosp',ok_cosp_omp)
+
+    !
+    !Config Key  = ok_airs
+    !Config Desc =
+    !Config Def  = .false.
+    !Config Help =
+    !
+    ok_airs_omp = .false.
+    call getin('ok_airs',ok_airs_omp)
 
     !
@@ -2111,4 +2129,5 @@
     ecrit_ISCCP = ecrit_ISCCP_omp
     freq_COSP = freq_COSP_omp
+    freq_AIRS = freq_AIRS_omp
     ok_ade = ok_ade_omp
     ok_aie = ok_aie_omp
@@ -2161,4 +2180,6 @@
     type_run = type_run_omp
     ok_cosp = ok_cosp_omp
+    ok_airs = ok_airs_omp
+
     ok_mensuelCOSP = ok_mensuelCOSP_omp
     ok_journeCOSP = ok_journeCOSP_omp
@@ -2308,4 +2329,5 @@
     write(lunout,*)' Frequence appel simulateur ISCCP, ecrit_ISCCP =', ecrit_ISCCP
     write(lunout,*)' Frequence appel simulateur COSP, freq_COSP =', freq_COSP
+    write(lunout,*)' Frequence appel simulateur AIRS, freq_AIRS =', freq_AIRS
     write(lunout,*)' Sortie bilan d''energie, ip_ebil_phy =', ip_ebil_phy
     write(lunout,*)' Excentricite = ',R_ecc
@@ -2423,4 +2445,6 @@
     write(lunout,*)' type_run = ',type_run 
     write(lunout,*)' ok_cosp = ',ok_cosp
+    write(lunout,*)' ok_airs = ',ok_airs
+
     write(lunout,*)' ok_mensuelCOSP = ',ok_mensuelCOSP
     write(lunout,*)' ok_journeCOSP = ',ok_journeCOSP
Index: /LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 2579)
+++ /LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 2580)
@@ -41,4 +41,173 @@
 
 !!! 2D
+
+! Marine
+
+  TYPE(ctrl_out), SAVE :: o_alt_tropo = ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'alt_tropo','Tropopause pressure','hPa',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+  TYPE(ctrl_out), SAVE :: o_map_prop_hc = ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_prop_hc','Proportion of high clouds',' ',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+  TYPE(ctrl_out), SAVE :: o_map_prop_hist = &
+  ctrl_out((/1,1,1,1,1,1,10,10,10/),&
+  'map_prop_hist','Proportion of high ice semi-transp clouds',' ',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+  TYPE(ctrl_out), SAVE :: o_map_emis_hc = &
+  ctrl_out((/1,1,1,1,1,1,10,10,10/),&
+  'map_emis_hc','Emissivity of high clouds',' ',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+  TYPE(ctrl_out), SAVE :: o_map_iwp_hc = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_iwp_hc','Ice water path of high clouds','g/m2',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+  TYPE(ctrl_out), SAVE :: o_map_deltaz_hc = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_deltaz_hc','geom thickness of high clouds','m',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+  TYPE(ctrl_out), SAVE :: o_map_pcld_hc = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_pcld_hc','cloud pressure of high clouds','hPa',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+   TYPE(ctrl_out), SAVE :: o_map_tcld_hc = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_tcld_hc','cloud temperature of high clouds','K',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+
+  TYPE(ctrl_out), SAVE :: o_map_emis_hist = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_emis_hist','Emissivity of high ice st clouds',' ',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+  TYPE(ctrl_out), SAVE :: o_map_iwp_hist = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_iwp_hist','Ice water path of high ice st clouds','g/m2',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+  TYPE(ctrl_out), SAVE :: o_map_deltaz_hist = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_deltaz_hist','geom thickness of high ice st clouds','m',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+  TYPE(ctrl_out), SAVE :: o_map_rad_hist = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_rad_hist','ice crystals radius in high ice st clouds','µm',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+
+ TYPE(ctrl_out), SAVE :: o_map_emis_Cb = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_emis_Cb','Emissivity of high Cb clouds',' ',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+ TYPE(ctrl_out), SAVE :: o_map_pcld_Cb = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_pcld_Cb','cloud pressure of high Cb clouds','hPa',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+ TYPE(ctrl_out), SAVE :: o_map_tcld_Cb = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_tcld_Cb','cloud temperature of high Cb clouds','K',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+
+ TYPE(ctrl_out), SAVE :: o_map_emis_Anv = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_emis_Anv','Emissivity of high Anv clouds',' ',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+ TYPE(ctrl_out), SAVE :: o_map_pcld_Anv = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_pcld_Anv','cloud pressure of high Anv clouds','hPa',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+  TYPE(ctrl_out), SAVE :: o_map_tcld_Anv = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_tcld_Anv','cloud temperature of high Anv clouds','K',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+  TYPE(ctrl_out), SAVE :: o_map_emis_ThCi = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_emis_ThCi','Emissivity of high ThCi clouds',' ',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+  TYPE(ctrl_out), SAVE :: o_map_pcld_ThCi = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_pcld_ThCi','cloud pressure of high ThCi clouds','hPa',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+  TYPE(ctrl_out), SAVE :: o_map_tcld_ThCi = &
+  ctrl_out((/10,10,1,10,10,10,10,10,10/),&
+  'map_tcld_ThCi','cloud temperature of high ThCi clouds','K',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+   TYPE(ctrl_out), SAVE :: o_map_ntot = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_ntot','total AIRS cloud fraction',' ',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+  TYPE(ctrl_out), SAVE :: o_map_hc = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_hc','high clouds AIRS cloud fraction',' ',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+  TYPE(ctrl_out), SAVE :: o_map_hist = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_hist','high clouds ice st AIRS cloud fraction',' ',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+  TYPE(ctrl_out), SAVE :: o_map_Cb = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_Cb','high clouds Cb AIRS cloud fraction',' ',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+ TYPE(ctrl_out), SAVE :: o_map_ThCi = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_ThCi','high clouds ThCi AIRS cloud fraction',' ',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+ TYPE(ctrl_out), SAVE :: o_map_Anv = &
+  ctrl_out((/1,1,1,1,1,10,10,10,10/),&
+  'map_Anv','high clouds Anv AIRS cloud fraction',' ',&
+   (/ "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)", "inst(X)",&
+    "inst(X)", "inst(X)","inst(X)" /))
+
+
+! Fin Marine
+
   TYPE(ctrl_out), SAVE :: o_flat = ctrl_out((/ 5, 1, 10, 10, 5, 10, 11, 11, 11 /), &
     'flat', 'Latent heat flux', 'W/m2', (/ ('', i=1, 9) /))
Index: /LMDZ5/trunk/libf/phylmd/phys_output_var_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phys_output_var_mod.F90	(revision 2579)
+++ /LMDZ5/trunk/libf/phylmd/phys_output_var_mod.F90	(revision 2580)
@@ -24,4 +24,20 @@
   REAL, ALLOCATABLE :: bils_latent(:) ! bilan de chaleur au sol
   !$OMP THREADPRIVATE(bils_ec,bils_ech,bils_tke,bils_diss,bils_kinetic,bils_enthalp,bils_latent)
+
+! Marine
+! Variables de sortie du simulateur AIRS
+
+  REAL, ALLOCATABLE :: map_prop_hc(:),map_prop_hist(:)
+  REAL, ALLOCATABLE :: alt_tropo(:)
+  REAL, ALLOCATABLE :: &
+     & map_emis_hc(:),map_iwp_hc(:),map_deltaz_hc(:),&
+     & map_pcld_hc(:),map_tcld_hc(:)
+  REAL, ALLOCATABLE :: map_emis_hist(:),map_iwp_hist(:),map_deltaz_hist(:),&
+     & map_rad_hist(:)
+  REAL, ALLOCATABLE :: map_ntot(:),map_hc(:),map_hist(:)
+  REAL, ALLOCATABLE :: map_Cb(:),map_ThCi(:),map_Anv(:)
+  REAL, ALLOCATABLE :: map_emis_Cb(:),map_pcld_Cb(:),map_tcld_Cb(:)
+  REAL, ALLOCATABLE :: map_emis_ThCi(:),map_pcld_ThCi(:),map_tcld_ThCi(:)
+  REAL, ALLOCATABLE :: map_emis_Anv(:),map_pcld_Anv(:),map_tcld_Anv(:)
 
 
@@ -103,4 +119,20 @@
     allocate (bils_ec(klon),bils_ech(klon),bils_tke(klon),bils_diss(klon),bils_kinetic(klon),bils_enthalp(klon),bils_latent(klon))
 
+! Marine
+! Variables de sortie simulateur AIRS
+
+      allocate (map_prop_hc(klon),map_prop_hist(klon))
+      allocate (alt_tropo(klon))
+      allocate &
+     & (map_emis_hc(klon),map_iwp_hc(klon),map_deltaz_hc(klon))
+      allocate (map_pcld_hc(klon),map_tcld_hc(klon))
+      allocate (map_emis_hist(klon),map_iwp_hist(klon),map_deltaz_hist(klon))
+      allocate (map_rad_hist(klon))
+      allocate (map_ntot(klon),map_hc(klon),map_hist(klon))
+      allocate (map_Cb(klon),map_ThCi(klon),map_Anv(klon))
+      allocate (map_emis_Cb(klon),map_pcld_Cb(klon),map_tcld_Cb(klon))
+      allocate (map_emis_ThCi(klon),map_pcld_ThCi(klon),map_tcld_ThCi(klon))
+      allocate (map_emis_Anv(klon),map_pcld_Anv(klon),map_tcld_Anv(klon))
+
     IF (ok_hines) allocate(zustr_gwd_hines(klon), zvstr_gwd_hines(klon))
     IF (.not.ok_hines.and.ok_gwd_rando) &
@@ -118,4 +150,20 @@
     deallocate (bils_ec,bils_ech,bils_tke,bils_diss,bils_kinetic,bils_enthalp,bils_latent)
 
+! Marine
+! Variables de sortie simulateur AIRS
+
+      deallocate (map_prop_hc,map_prop_hist)
+      deallocate (alt_tropo)
+      deallocate &
+     & (map_emis_hc,map_iwp_hc,map_deltaz_hc)
+      deallocate (map_pcld_hc,map_tcld_hc)
+      deallocate (map_emis_hist,map_iwp_hist,map_deltaz_hist)
+      deallocate (map_rad_hist)
+      deallocate (map_ntot,map_hc,map_hist)
+      deallocate (map_Cb,map_ThCi,map_Anv)
+      deallocate (map_emis_Cb,map_pcld_Cb,map_tcld_Cb)
+      deallocate (map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi)
+      deallocate (map_emis_Anv,map_pcld_Anv,map_tcld_Anv)
+
   END SUBROUTINE phys_output_var_end
 
Index: /LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 2579)
+++ /LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 2580)
@@ -168,5 +168,16 @@
          o_sens_prec_sol_oce, o_sens_prec_sol_sic, &
          o_lat_prec_liq_oce, o_lat_prec_liq_sic, &
-         o_lat_prec_sol_oce, o_lat_prec_sol_sic
+         o_lat_prec_sol_oce, o_lat_prec_sol_sic, &
+! Marine
+         o_map_prop_hc, o_map_prop_hist, o_map_emis_hc, o_map_iwp_hc, &
+         o_map_deltaz_hc, o_map_pcld_hc, o_map_tcld_hc, &
+         o_map_emis_hist, o_map_iwp_hist, o_map_deltaz_hist, &
+         o_map_rad_hist, &
+         o_map_emis_Cb, o_map_pcld_Cb, o_map_tcld_Cb, &
+         o_map_emis_ThCi, o_map_pcld_ThCi, o_map_tcld_ThCi, &
+         o_map_emis_Anv, o_map_pcld_Anv, o_map_tcld_Anv, &
+         o_map_ntot, o_map_hc,o_map_hist,o_map_Cb,o_map_ThCi,o_map_Anv, &
+         o_alt_tropo
+
 
     USE phys_state_var_mod, only: pctsrf, paire_ter, rain_fall, snow_fall, &
@@ -264,5 +275,18 @@
          zustr_gwd_hines, zvstr_gwd_hines,zustr_gwd_rando, zvstr_gwd_rando, &
          zustr_gwd_front, zvstr_gwd_front,     &
-         sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o
+         sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o, &
+! Marine
+         map_prop_hc, map_prop_hist, &
+         map_emis_hc,map_iwp_hc,map_deltaz_hc,&
+         map_pcld_hc,map_tcld_hc,&
+         map_emis_hist,map_iwp_hist,map_deltaz_hist,&
+         map_rad_hist,&
+         map_ntot,map_hc,map_hist,&
+         map_Cb,map_ThCi,map_Anv,&
+         map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
+         map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
+         map_emis_Anv,map_pcld_Anv,map_tcld_Anv, &
+         alt_tropo
+
   
 
@@ -380,4 +404,42 @@
        CALL histwrite_phy(o_aireTER, paire_ter)
 !!! Champs 2D !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Simulateur AIRS
+     IF (ok_airs) then
+      CALL histwrite_phy(o_alt_tropo,alt_tropo)
+ 
+      CALL histwrite_phy(o_map_prop_hc,map_prop_hc)
+      CALL histwrite_phy(o_map_prop_hist,map_prop_hist)
+
+      CALL histwrite_phy(o_map_emis_hc,map_emis_hc)
+      CALL histwrite_phy(o_map_iwp_hc,map_iwp_hc)
+      CALL histwrite_phy(o_map_deltaz_hc,map_deltaz_hc)
+      CALL histwrite_phy(o_map_pcld_hc,map_pcld_hc)
+      CALL histwrite_phy(o_map_tcld_hc,map_tcld_hc)
+
+      CALL histwrite_phy(o_map_emis_hist,map_emis_hist)
+      CALL histwrite_phy(o_map_iwp_hist,map_iwp_hist)
+      CALL histwrite_phy(o_map_deltaz_hist,map_deltaz_hist)
+
+      CALL histwrite_phy(o_map_ntot,map_ntot)
+      CALL histwrite_phy(o_map_hc,map_hc)
+      CALL histwrite_phy(o_map_hist,map_hist)
+
+      CALL histwrite_phy(o_map_Cb,map_Cb)
+      CALL histwrite_phy(o_map_ThCi,map_ThCi)
+      CALL histwrite_phy(o_map_Anv,map_Anv)
+
+      CALL histwrite_phy(o_map_emis_Cb,map_emis_Cb)
+      CALL histwrite_phy(o_map_pcld_Cb,map_pcld_Cb)
+      CALL histwrite_phy(o_map_tcld_Cb,map_tcld_Cb)
+
+      CALL histwrite_phy(o_map_emis_ThCi,map_emis_ThCi)
+      CALL histwrite_phy(o_map_pcld_ThCi,map_pcld_ThCi)
+      CALL histwrite_phy(o_map_tcld_ThCi,map_tcld_ThCi)
+
+      CALL histwrite_phy(o_map_emis_Anv,map_emis_Anv)
+      CALL histwrite_phy(o_map_pcld_Anv,map_pcld_Anv)
+      CALL histwrite_phy(o_map_tcld_Anv,map_tcld_Anv)
+     ENDIF
+
        CALL histwrite_phy(o_flat, zxfluxlat)
        CALL histwrite_phy(o_ptstar, ptstar)
Index: /LMDZ5/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/physiq_mod.F90	(revision 2579)
+++ /LMDZ5/trunk/libf/phylmd/physiq_mod.F90	(revision 2580)
@@ -3948,4 +3948,28 @@
 #endif
     ENDIF  !ok_cosp
+
+
+! Marine
+
+  IF (ok_airs) then
+
+  IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/dtime)).EQ.0) THEN
+  write(*,*) 'je vais appeler simu_airs, ok_airs, freq_airs=', &
+     & ok_airs, freq_airs
+  call simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
+     & map_prop_hc,map_prop_hist,&
+     & map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
+     & map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
+     & map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
+     & map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
+     & map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
+     & map_ntot,map_hc,map_hist,&
+     & map_Cb,map_ThCi,map_Anv,&
+     & alt_tropo )
+  ENDIF
+
+  ENDIF  ! ok_airs
+
+
     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !AA
Index: /LMDZ5/trunk/libf/phylmd/simu_airs.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/simu_airs.F90	(revision 2580)
+++ /LMDZ5/trunk/libf/phylmd/simu_airs.F90	(revision 2580)
@@ -0,0 +1,1300 @@
+        
+        module m_simu_airs
+
+        implicit none
+
+        real, parameter :: tau_thresh = 0.05 ! seuil nuages detectables
+        real, parameter :: p_thresh = 445.   ! seuil nuages hauts
+        real, parameter :: em_min = 0.2      ! seuils nuages semi-transp
+        real, parameter :: em_max = 0.85
+        real, parameter :: undef = 999.
+
+        contains
+
+        real function search_tropopause(P,T,alt,N) result(P_tropo)
+! this function searches for the tropopause pressure in [hPa].
+! The search is based on ideology described in
+! Reichler et al., Determining the tropopause height from gridded data,
+! GRL, 30(20) 2042, doi:10.1029/2003GL018240, 2003
+
+        integer N,i,i_lev,first_point,exit_flag,i_dir
+        real P(N),T(N),alt(N),slope(N)
+        real P_min, P_max, slope_limit,slope_2km, &
+     & delta_alt_limit,tmp,delta_alt
+        parameter(P_min=75.0, P_max=470.0)   ! hPa
+        parameter(slope_limit=0.002)         ! 2 K/km converted to K/m
+        parameter(delta_alt_limit=2000.0)    ! 2000 meters
+
+        do i=1,N-1
+        slope(i)=-(T(i+1)-T(i))/(alt(i+1)-alt(i))
+        end do
+        slope(N)=slope(N-1)
+
+        if (P(1).gt.P(N)) then
+        i_dir= 1
+        i=1
+        else
+        i_dir=-1
+        i=N
+        end if
+
+        first_point=0
+        exit_flag=0
+        do while(exit_flag.eq.0)
+        if (P(i).ge.P_min.and.P(i).le.P_max) then
+        if (first_point.gt.0) then
+        delta_alt=alt(i)-alt(first_point)
+        if (delta_alt.ge.delta_alt_limit) then
+        slope_2km=(T(first_point)-T(i))/delta_alt
+        if (slope_2km.lt.slope_limit) then
+        exit_flag=1
+        else
+        first_point=0
+        end if
+        end if
+        end if
+        if (first_point.eq.0.and.slope(i).lt.slope_limit) first_point=i
+        end if
+        i=i+i_dir
+        if (i.le.1.or.i.ge.N) exit_flag=1
+        end do
+
+        if (first_point.le.0) P_tropo=65.4321
+        if (first_point.eq.1) P_tropo=64.5432
+        if (first_point.gt.1) then
+        tmp=(slope_limit-slope(first_point))/(slope(first_point+1)- &
+     & slope(first_point))*(P(first_point+1)-P(first_point))
+        P_tropo=P(first_point)+tmp
+        ! print*, 'P_tropo= ', tmp, P(first_point), P_tropo
+        end if
+
+! Ajout Marine
+        if (P_tropo .lt. 60. .or. P_tropo .gt. 470.) then
+        P_tropo = 999.
+        endif
+
+        return
+        end function search_tropopause
+
+
+
+        subroutine cloud_structure(len_cs, rneb_cs, temp_cs, &
+     & emis_cs, iwco_cs, &
+     & pres_cs, dz_cs, rhodz_cs, rad_cs, &
+     & cc_tot_cs, cc_hc_cs, cc_hist_cs, &
+     & cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs, &
+     & pcld_hc_cs, tcld_hc_cs, &
+     & em_hc_cs, iwp_hc_cs, deltaz_hc_cs, &
+     & pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs, &
+     & pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs, &
+     & pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs, &
+     & em_hist_cs, iwp_hist_cs, deltaz_hist_cs, rad_hist_cs)
+
+
+     
+        integer :: i, n, nss
+
+        integer, intent(in) :: len_cs
+        real, dimension(:), intent(in) :: rneb_cs, temp_cs
+        real, dimension(:), intent(in) :: emis_cs, iwco_cs, rad_cs
+        real, dimension(:), intent(in) :: pres_cs, dz_cs, rhodz_cs
+
+        real, intent(out) :: cc_tot_cs, cc_hc_cs, cc_hist_cs, &
+     & cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs, &
+     & pcld_hc_cs, tcld_hc_cs, em_hc_cs, iwp_hc_cs, &
+     & pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs, &
+     & pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs, &
+     & pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs, &
+     & em_hist_cs, iwp_hist_cs, &
+     & deltaz_hc_cs, deltaz_hist_cs, rad_hist_cs
+
+        real, dimension(len_cs) :: rneb_ord
+        real :: rneb_min
+
+        real, dimension(:), allocatable :: s, s_hc, s_hist, rneb_max
+        real, dimension(:), allocatable :: sCb, sThCi, sAnv
+        real, dimension(:), allocatable :: iwp_ss, pcld_ss, tcld_ss,&
+     & emis_ss
+        real, dimension(:), allocatable :: deltaz_ss, rad_ss
+
+
+        write(*,*) 'dans cloud_structure'
+
+        call ordonne(len_cs, rneb_cs, rneb_ord)
+       
+
+! Definition des sous_sections
+
+        rneb_min = rneb_ord(1)
+        nss = 1
+
+        do i = 2, size(rneb_cs)
+        if (rneb_ord(i) .gt. rneb_min) then
+        nss = nss+1
+        rneb_min = rneb_ord(i)
+        endif
+        enddo
+
+        allocate (s(nss))
+        allocate (s_hc(nss))
+        allocate (s_hist(nss))
+        allocate (rneb_max(nss))
+        allocate (emis_ss(nss))
+        allocate (pcld_ss(nss))
+        allocate (tcld_ss(nss))
+        allocate (iwp_ss(nss))
+        allocate (deltaz_ss(nss))
+        allocate (rad_ss(nss))
+        allocate (sCb(nss))
+        allocate (sThCi(nss))
+        allocate (sAnv(nss))
+
+        rneb_min = rneb_ord(1)
+        n = 1
+        s(1) = rneb_ord(1)
+        s_hc(1) = rneb_ord(1)
+        s_hist(1) = rneb_ord(1)
+        sCb(1) = rneb_ord(1)
+        sThCi(1) = rneb_ord(1)
+        sAnv(1) = rneb_ord(1)
+
+        rneb_max(1) = rneb_ord(1)
+
+        do i = 2, size(rneb_cs)
+        if (rneb_ord(i) .gt. rneb_min) then
+        n = n+1
+        s(n) = rneb_ord(i)-rneb_min
+        s_hc(n) = rneb_ord(i)-rneb_min
+        s_hist(n) = rneb_ord(i)-rneb_min
+        sCb(n) = rneb_ord(i)-rneb_min
+        sThCi(n) = rneb_ord(i)-rneb_min
+        sAnv(n) = rneb_ord(i)-rneb_min
+
+        rneb_max(n) = rneb_ord(i)
+        rneb_min = rneb_ord(i)
+        endif
+        enddo
+
+! Appel de sous_section
+
+        do i = 1, nss
+         call sous_section(len_cs, rneb_cs, temp_cs, &
+     &  emis_cs, iwco_cs, &
+     &  pres_cs, dz_cs, rhodz_cs, rad_cs, rneb_ord, &
+     &  rneb_max(i),s(i),s_hc(i),s_hist(i), &
+     &  sCb(i), sThCi(i), sAnv(i), &
+     &  emis_ss(i), &
+     &  pcld_ss(i), tcld_ss(i), iwp_ss(i), deltaz_ss(i), rad_ss(i))
+        enddo
+
+! Caracteristiques de la structure nuageuse
+
+        cc_tot_cs = 0.
+        cc_hc_cs = 0.
+        cc_hist_cs = 0.
+
+        cc_Cb_cs = 0.
+        cc_ThCi_cs = 0.
+        cc_Anv_cs = 0.
+
+        em_hc_cs = 0.
+        iwp_hc_cs = 0.
+        deltaz_hc_cs = 0.
+
+        em_hist_cs = 0.
+        iwp_hist_cs = 0.
+        deltaz_hist_cs = 0.
+        rad_hist_cs = 0.
+
+        pcld_hc_cs = 0.
+        tcld_hc_cs = 0.
+
+        pcld_Cb_cs = 0.
+        tcld_Cb_cs = 0.
+        em_Cb_cs = 0.
+
+        pcld_ThCi_cs = 0.
+        tcld_ThCi_cs = 0.
+        em_ThCi_cs = 0.
+
+        pcld_Anv_cs = 0.
+        tcld_Anv_cs = 0.
+        em_Anv_cs = 0.
+
+         do i = 1, nss
+
+        cc_tot_cs = cc_tot_cs + s(i)
+        cc_hc_cs = cc_hc_cs + s_hc(i)
+        cc_hist_cs = cc_hist_cs + s_hist(i)
+
+        cc_Cb_cs = cc_Cb_cs + sCb(i)
+        cc_ThCi_cs = cc_ThCi_cs + sThCi(i)
+        cc_Anv_cs = cc_Anv_cs + sAnv(i)
+
+        iwp_hc_cs = iwp_hc_cs + s_hc(i)*iwp_ss(i)
+        em_hc_cs = em_hc_cs + s_hc(i)*emis_ss(i)
+        deltaz_hc_cs = deltaz_hc_cs + s_hc(i)*deltaz_ss(i)
+
+        iwp_hist_cs = iwp_hist_cs + s_hist(i)*iwp_ss(i)
+        em_hist_cs = em_hist_cs + s_hist(i)*emis_ss(i)
+        deltaz_hist_cs = deltaz_hist_cs + s_hist(i)*deltaz_ss(i)
+        rad_hist_cs = rad_hist_cs + s_hist(i)*rad_ss(i)
+
+        pcld_hc_cs = pcld_hc_cs + s_hc(i)*pcld_ss(i)
+        tcld_hc_cs = tcld_hc_cs + s_hc(i)*tcld_ss(i)
+
+        pcld_Cb_cs = pcld_Cb_cs + sCb(i)*pcld_ss(i)
+        tcld_Cb_cs = tcld_Cb_cs + sCb(i)*tcld_ss(i)
+        em_Cb_cs = em_Cb_cs + sCb(i)*emis_ss(i)
+
+        pcld_ThCi_cs = pcld_ThCi_cs + sThCi(i)*pcld_ss(i)
+        tcld_ThCi_cs = tcld_ThCi_cs + sThCi(i)*tcld_ss(i)
+        em_ThCi_cs = em_ThCi_cs + sThCi(i)*emis_ss(i)
+
+        pcld_Anv_cs = pcld_Anv_cs + sAnv(i)*pcld_ss(i)
+        tcld_Anv_cs = tcld_Anv_cs + sAnv(i)*tcld_ss(i)
+        em_Anv_cs = em_Anv_cs + sAnv(i)*emis_ss(i)
+
+        enddo
+
+        deallocate(s)
+        deallocate (s_hc)
+        deallocate (s_hist)
+        deallocate (rneb_max)
+        deallocate (emis_ss)
+        deallocate (pcld_ss)
+        deallocate (tcld_ss)
+        deallocate (iwp_ss)
+        deallocate (deltaz_ss)
+        deallocate (rad_ss)
+        deallocate (sCb)
+        deallocate (sThCi)
+        deallocate (sAnv)
+
+       call normal_undef(pcld_hc_cs,cc_hc_cs)
+       call normal_undef(tcld_hc_cs,cc_hc_cs)
+       call normal_undef(iwp_hc_cs,cc_hc_cs)
+       call normal_undef(em_hc_cs,cc_hc_cs)
+       call normal_undef(deltaz_hc_cs,cc_hc_cs)
+       
+       call normal_undef(iwp_hist_cs,cc_hist_cs)
+       call normal_undef(em_hist_cs,cc_hist_cs)
+       call normal_undef(deltaz_hist_cs,cc_hist_cs)
+       call normal_undef(rad_hist_cs,cc_hist_cs)
+
+       call normal_undef(pcld_Cb_cs,cc_Cb_cs)
+       call normal_undef(tcld_Cb_cs,cc_Cb_cs)
+       call normal_undef(em_Cb_cs,cc_Cb_cs)
+
+       call normal_undef(pcld_ThCi_cs,cc_ThCi_cs)
+       call normal_undef(tcld_ThCi_cs,cc_ThCi_cs)
+       call normal_undef(em_ThCi_cs,cc_ThCi_cs)
+    
+       call normal_undef(pcld_Anv_cs,cc_Anv_cs)
+       call normal_undef(tcld_Anv_cs,cc_Anv_cs)
+       call normal_undef(em_Anv_cs,cc_Anv_cs)
+
+
+! Tests
+
+        if (cc_tot_cs .gt. maxval(rneb_cs) .and. &
+     & abs(cc_tot_cs-maxval(rneb_cs)) .gt. 1.e-4 )  then
+        write(*,*) 'cc_tot_cs > max rneb_cs'
+        write(*,*) cc_tot_cs, maxval(rneb_cs)
+        STOP
+        endif
+
+        if (iwp_hc_cs .lt. 0.) then
+        write(*,*) 'cloud_structure:'
+        write(*,*) 'iwp_hc_cs < 0'
+        STOP
+        endif
+ 
+        end subroutine cloud_structure
+
+
+        subroutine normal_undef(num, den)
+
+        real, intent(in) :: den
+        real, intent(inout) :: num
+
+        if (den .ne. 0) then
+        num = num/den
+        else
+        num = undef
+        endif
+
+        end subroutine normal_undef
+
+
+        subroutine normal2_undef(res,num,den)
+
+        real, intent(in) :: den
+        real, intent(in) :: num
+        real, intent(out) :: res
+
+        if (den .ne. 0.) then
+        res = num/den
+        else
+        res = undef
+        endif
+
+        end subroutine normal2_undef
+
+
+        subroutine sous_section(len_cs, rneb_cs, temp_cs, &
+     & emis_cs, iwco_cs, &
+     & pres_cs, dz_cs, rhodz_cs, rad_cs, rneb_ord, &
+     & rnebmax, stot, shc, shist, &
+     & sCb, sThCi, sAnv, &
+     & emis, pcld, tcld, iwp, deltaz, rad)
+
+        integer, intent(in) :: len_cs
+        real, dimension(len_cs), intent(in) :: rneb_cs, temp_cs
+        real, dimension(len_cs), intent(in) :: emis_cs, iwco_cs, &
+     & rneb_ord
+        real, dimension(len_cs), intent(in) :: pres_cs, dz_cs, rad_cs
+        real, dimension(len_cs), intent(in) :: rhodz_cs
+        real, dimension(len_cs) :: tau_cs, w
+        real, intent(in) :: rnebmax
+        real, intent(inout) :: stot, shc, shist
+        real, intent(inout) :: sCb, sThCi, sAnv
+        real, intent(out) :: emis, pcld, tcld, iwp, deltaz, rad
+
+        integer :: i, ideb, ibeg, iend, nuage, visible
+        real :: som, som_tau, som_iwc, som_dz, som_rad, tau
+
+
+! Ponderation: 1 pour les nuages, 0 pour les trous
+
+        do i = 1, len_cs
+        if (rneb_cs(i) .ge. rnebmax) then
+        w(i) = 1.
+        else
+        w(i) = 0.
+        endif
+        enddo
+
+! Calcul des epaisseurs optiques a partir des emissivites
+
+        som = 0.
+        do i = 1, len_cs
+        if (emis_cs(i) .eq. 1.) then
+        tau_cs(i) = 10.
+        else
+        tau_cs(i) = -log(1.-emis_cs(i))
+        endif
+        som = som+tau_cs(i)
+        enddo
+
+
+        ideb = 1
+        nuage = 0
+        visible = 0
+
+
+! Boucle sur les nuages
+        do while (ideb .ne. 0 .and. ideb .le. len_cs)   
+
+
+! Definition d'un nuage de la sous-section
+
+        call topbot(ideb, w, ibeg, iend)
+        ideb = iend+1
+
+        if (ibeg .gt. 0) then
+
+        nuage = nuage + 1
+
+! On determine les caracteristiques du nuage
+! (ep. optique, ice water path, pression, temperature)
+
+        call caract(ibeg, iend, temp_cs, tau_cs, iwco_cs, &
+     & pres_cs, dz_cs, rhodz_cs, rad_cs, pcld, tcld, &
+     & som_tau, som_iwc, som_dz, som_rad)
+
+! On masque le nuage s'il n'est pas detectable
+
+        call masque(ibeg, iend, som_tau, visible, w)
+
+        endif
+
+! Fin boucle nuages
+        enddo
+
+
+! Analyse du nuage detecte
+
+        call topbot(1, w, ibeg, iend)
+
+        if (ibeg .gt. 0) then
+
+        call caract(ibeg, iend, temp_cs, tau_cs, iwco_cs, &
+     & pres_cs, dz_cs, rhodz_cs, rad_cs, pcld, tcld, &
+     & som_tau, som_iwc, som_dz, som_rad)
+
+        tau = som_tau
+        emis = 1. - exp(-tau)
+        iwp = som_iwc
+        deltaz = som_dz
+        rad = som_rad
+
+        if (pcld .gt. p_thresh) then
+
+        shc = 0.
+        shist = 0.
+        sCb = 0.
+        sThCi = 0.
+        sAnv = 0.
+
+        else
+
+        if (emis .lt. em_min .or. emis .gt. em_max  &
+     & .or. tcld .gt. 230.) then
+        shist = 0.
+        endif
+
+        if (emis .lt. 0.98) then
+        sCb = 0.
+        endif
+
+        if (emis .gt. 0.5 .or. emis .lt. 0.1) then
+        sThCi = 0.
+        endif
+
+        if (emis .le. 0.5 .or. emis .ge. 0.98) then
+        sAnv = 0.
+        endif
+
+        endif
+
+        else
+
+        tau = 0.
+        emis = 0.
+        iwp = 0.
+        deltaz = 0.
+        pcld = 0.
+        tcld = 0.
+        stot = 0.
+        shc = 0.
+        shist = 0.
+        rad = 0.
+        sCb = 0.
+        sThCi = 0.
+        sAnv = 0.
+
+        endif
+
+
+! Tests
+
+        if (iwp .lt. 0.) then
+        write(*,*) 'ideb iwp =', ideb, iwp
+        STOP
+        endif
+
+        if (deltaz .lt. 0.) then
+        write(*,*) 'ideb deltaz =', ideb, deltaz
+        STOP
+        endif
+
+        if (emis .lt. 0.048 .and. emis .ne. 0.) then
+        write(*,*) 'ideb emis =', ideb, emis
+        STOP
+        endif
+
+        end subroutine sous_section
+
+
+        subroutine masque (ibeg, iend, som_tau, &
+     & visible, w)
+
+        integer, intent(in) :: ibeg, iend
+        real, intent(in) :: som_tau
+
+        integer, intent(inout) :: visible
+        real, dimension(:), intent(inout) :: w
+
+        integer :: i
+
+
+
+! Masque
+
+! Cas ou il n'y a pas de nuage visible au-dessus
+
+        if (visible .eq. 0) then
+
+        if (som_tau .lt. tau_thresh) then
+        do i = ibeg, iend
+        w(i) = 0.
+        enddo
+        else
+        visible = 1
+        endif
+
+! Cas ou il y a un nuage visible au-dessus
+
+        else
+
+        do i = ibeg, iend
+        w(i) = 0.
+        enddo
+
+        endif
+
+
+        end subroutine masque
+
+
+         subroutine caract (ibeg, iend, temp_cs, tau_cs, iwco_cs, &
+     & pres_cs, dz_cs, rhodz_cs, rad_cs, pcld, tcld, &
+     & som_tau, som_iwc, som_dz, som_rad)
+
+        integer, intent(in) :: ibeg, iend
+        real, dimension(:), intent(in) :: tau_cs, iwco_cs, temp_cs
+        real, dimension(:), intent(in) :: pres_cs, dz_cs, rad_cs
+        real, dimension(:), intent(in) :: rhodz_cs
+        real, intent(out) :: som_tau, som_iwc, som_dz, som_rad
+        real , intent(out) :: pcld, tcld
+
+        integer :: i, ibase, imid
+
+! Somme des epaisseurs optiques et des contenus en glace sur le nuage
+
+        som_tau = 0.
+        som_iwc = 0.
+        som_dz = 0.
+        som_rad = 0.
+        ibase = -100
+
+        do i = ibeg, iend
+
+        som_tau = som_tau + tau_cs(i)
+
+        som_dz = som_dz + dz_cs(i)
+        som_iwc = som_iwc + iwco_cs(i)*1000*rhodz_cs(i)  ! en g/m2
+        som_rad = som_rad + rad_cs(i)*dz_cs(i)
+
+        if (som_tau .gt. 3. .and. ibase .eq. -100) then
+        ibase = i-1
+        endif
+
+        enddo
+
+        if (som_dz .ne. 0.) then
+        som_rad = som_rad/som_dz
+        else
+        write(*,*) 'som_dez = 0 STOP'
+        write(*,*) 'ibeg, iend =', ibeg, iend
+        do i = ibeg, iend
+        write(*,*) dz_cs(i), rhodz_cs(i)
+        enddo
+        STOP
+        endif
+
+! Determination de Pcld et Tcld
+
+       if (ibase .lt. ibeg) then
+       ibase = ibeg
+       endif
+
+       imid = (ibeg+ibase)/2
+
+       pcld = pres_cs(imid)/100.        ! pcld en hPa
+       tcld = temp_cs(imid)
+
+
+       end subroutine caract
+ 
+        subroutine topbot(ideb,w,ibeg,iend)
+
+        integer, intent(in) :: ideb
+        real, dimension(:), intent(in) :: w
+        integer, intent(out) :: ibeg, iend
+
+        integer :: i, itest
+
+        itest = 0
+        ibeg = 0
+        iend = 0
+
+        do i = ideb, size(w)
+
+        if (w(i) .eq. 1. .and. itest .eq. 0) then
+        ibeg = i
+        itest = 1
+        endif
+
+        enddo
+
+
+        i = ibeg
+        do while (w(i) .eq. 1. .and. i .le. size(w))
+        i = i+1
+        enddo
+        iend = i-1
+
+
+        end subroutine topbot
+
+        subroutine ordonne(len_cs, rneb_cs, rneb_ord)
+
+        integer, intent(in) :: len_cs
+        real, dimension(:), intent(in) :: rneb_cs
+        real, dimension(:), intent(out) :: rneb_ord
+
+        integer :: i, j, ind_min
+
+        real, dimension(len_cs) :: rneb
+        real :: rneb_max
+
+
+        do i = 1, size(rneb_cs)
+        rneb(i) = rneb_cs(i)
+        enddo
+
+        do j = 1, size(rneb_cs)
+
+        rneb_max = 100.
+
+        do i = 1, size(rneb_cs)
+        if (rneb(i) .lt. rneb_max) then
+        rneb_max = rneb(i)
+        ind_min = i
+        endif
+        enddo
+
+        rneb_ord(j) = rneb_max
+        rneb(ind_min) = 100.
+
+        enddo
+        
+        end subroutine ordonne
+
+ 
+        subroutine sim_mesh(rneb_1D, temp_1D, emis_1D, &
+     & iwcon_1D, rad_1D, &
+     & pres, dz, &
+     & rhodz_1D, cc_tot_mesh, cc_hc_mesh, cc_hist_mesh, pcld_hc_mesh,&
+     & tcld_hc_mesh, &
+     & em_hc_mesh, iwp_hc_mesh, deltaz_hc_mesh, &
+     & cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh, &
+     & pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh, &
+     & pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh, &
+     & pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh, &
+     & em_hist_mesh, iwp_hist_mesh, deltaz_hist_mesh, rad_hist_mesh)
+
+       USE dimphy
+
+       real, dimension(klev), intent(in) :: rneb_1D, temp_1D, emis_1D, &
+     & iwcon_1D, rad_1D
+        real, dimension(klev), intent(in) :: pres, dz, rhodz_1D
+        real, intent(out) :: cc_tot_mesh, cc_hc_mesh, cc_hist_mesh
+        real, intent(out) :: cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh
+        real, intent(out) :: em_hc_mesh, pcld_hc_mesh, tcld_hc_mesh, &
+     & iwp_hc_mesh
+
+        real, intent(out) :: pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh
+        real, intent(out) :: pcld_ThCi_mesh, tcld_ThCi_mesh, &
+     & em_ThCi_mesh
+        real, intent(out) :: pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh
+
+        real, intent(out) :: em_hist_mesh, iwp_hist_mesh, rad_hist_mesh
+        real, intent(out) :: deltaz_hc_mesh, deltaz_hist_mesh
+
+        real, dimension(:), allocatable :: rneb_cs, temp_cs, emis_cs, &
+     & iwco_cs
+        real, dimension(:), allocatable :: pres_cs, dz_cs, rad_cs, &
+     & rhodz_cs
+
+        integer :: i,j,l
+        integer :: ltop, itop, ibot, num_cs, N_cs, len_cs, ics
+
+        real :: som_emi_hc,som_pcl_hc,som_tcl_hc,som_iwp_hc,som_hc,&
+     & som_hist
+        real :: som_emi_hist, som_iwp_hist, som_deltaz_hc, &
+     & som_deltaz_hist
+        real :: som_rad_hist
+        real :: som_Cb, som_ThCi, som_Anv
+        real :: som_emi_Cb, som_tcld_Cb, som_pcld_Cb
+        real :: som_emi_Anv, som_tcld_Anv, som_pcld_Anv
+        real :: som_emi_ThCi, som_tcld_ThCi, som_pcld_ThCi
+        real :: tsom_tot, tsom_hc, tsom_hist
+        real :: prod, prod_hh
+
+        real :: cc_tot_cs, cc_hc_cs, cc_hist_cs
+        real :: cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs
+        real :: pcld_hc_cs, tcld_hc_cs
+        real :: em_hc_cs, iwp_hc_cs, deltaz_hc_cs
+        real :: pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs
+        real :: pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs
+        real :: pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs
+        real :: em_hist_cs, iwp_hist_cs, deltaz_hist_cs, rad_hist_cs
+
+        real, dimension(klev) :: test_tot, test_hc, test_hist
+        real, dimension(klev) :: test_pcld, test_tcld, test_em, test_iwp
+
+
+        do j = 1, klev
+        write(*,*) 'simu_airs, j, rneb_1D =', rneb_1D(j)
+        enddo
+
+! Definition des structures nuageuses, de la plus haute a la plus basse
+
+        num_cs = 0
+        ltop = klev-1
+
+        prod = 1.
+
+        som_emi_hc = 0.
+        som_emi_hist = 0.
+        som_pcl_hc = 0.
+        som_tcl_hc = 0.
+        som_iwp_hc = 0.
+        som_iwp_hist = 0.
+        som_deltaz_hc = 0.
+        som_deltaz_hist = 0.
+        som_rad_hist = 0.
+        som_hc = 0.
+        som_hist = 0.
+
+        som_Cb = 0.
+        som_ThCi = 0.
+        som_Anv = 0.
+
+        som_emi_Cb = 0.
+        som_pcld_Cb = 0.
+        som_tcld_Cb = 0.
+
+        som_emi_ThCi = 0.
+        som_pcld_ThCi = 0.
+        som_tcld_ThCi = 0.
+
+        som_emi_Anv = 0.
+        som_pcld_Anv = 0.
+        som_tcld_Anv = 0.
+
+        tsom_tot = 0.
+        tsom_hc = 0.
+        tsom_hist = 0.
+
+        do while (ltop .ge. 1)   ! Boucle sur toute la colonne
+
+        itop = 0
+
+        do l = ltop,1,-1
+
+        if (itop .eq. 0 .and. rneb_1D(l) .gt. 0.001 ) then
+        itop = l
+        endif
+
+        enddo
+
+        ibot = itop
+
+        do while (rneb_1D(ibot) .gt. 0.001 .and. ibot .ge. 1)
+        ibot = ibot -1
+        enddo
+
+
+        ibot = ibot+1
+
+        if (itop .gt. 0) then    ! itop > 0
+
+        num_cs = num_cs +1
+        len_cs = itop-ibot+1
+
+! Allocation et definition des variables de la structure nuageuse
+! le premier indice denote ce qui est le plus haut
+
+        allocate (rneb_cs(len_cs))
+        allocate (temp_cs(len_cs))
+        allocate (emis_cs(len_cs))
+        allocate (iwco_cs(len_cs))
+        allocate (pres_cs(len_cs))
+        allocate (dz_cs(len_cs))
+        allocate (rad_cs(len_cs))
+        allocate (rhodz_cs(len_cs))
+
+        ics = 0
+
+        do i = itop, ibot, -1
+        ics = ics + 1
+        rneb_cs(ics) = rneb_1D(i)
+        temp_cs(ics) = temp_1D(i)
+        emis_cs(ics) = emis_1D(i)
+        iwco_cs(ics) = iwcon_1D(i)
+        rad_cs(ics) = rad_1D(i)
+        pres_cs(ics) = pres(i)
+        dz_cs(ics) = dz(i)
+        rhodz_cs(ics) = rhodz_1D(i)
+        enddo
+
+! Appel du sous_programme cloud_structure
+
+        call cloud_structure(len_cs,rneb_cs,temp_cs,emis_cs,iwco_cs,&
+     & pres_cs, dz_cs, rhodz_cs, rad_cs, &
+     & cc_tot_cs, cc_hc_cs, cc_hist_cs, &
+     & cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs, &
+     & pcld_hc_cs, tcld_hc_cs, &
+     & em_hc_cs, iwp_hc_cs, deltaz_hc_cs, &
+     & pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs, &
+     & pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs, &
+     & pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs, &
+     & em_hist_cs, iwp_hist_cs, deltaz_hist_cs, rad_hist_cs)
+
+
+        deallocate (rneb_cs)
+        deallocate (temp_cs)
+        deallocate (emis_cs)
+        deallocate (iwco_cs)
+        deallocate (pres_cs)
+        deallocate (dz_cs)
+        deallocate (rad_cs)
+        deallocate (rhodz_cs)
+
+
+! Pour la couverture nuageuse sur la maille
+
+        prod_hh = prod
+
+        prod = prod*(1.-cc_tot_cs)
+
+! Pour les autres variables definies sur la maille
+
+        som_emi_hc = som_emi_hc + em_hc_cs*cc_hc_cs*prod_hh
+        som_iwp_hc = som_iwp_hc + iwp_hc_cs*cc_hc_cs*prod_hh
+        som_deltaz_hc = som_deltaz_hc + deltaz_hc_cs*cc_hc_cs*prod_hh
+
+        som_emi_Cb = som_emi_Cb + em_Cb_cs*cc_Cb_cs*prod_hh
+        som_tcld_Cb = som_tcld_Cb + tcld_Cb_cs*cc_Cb_cs*prod_hh
+        som_pcld_Cb = som_pcld_Cb + pcld_Cb_cs*cc_Cb_cs*prod_hh
+
+        som_emi_ThCi = som_emi_ThCi + em_ThCi_cs*cc_ThCi_cs*prod_hh
+        som_tcld_ThCi = som_tcld_ThCi + tcld_ThCi_cs*cc_ThCi_cs*prod_hh
+        som_pcld_ThCi = som_pcld_ThCi + pcld_ThCi_cs*cc_ThCi_cs*prod_hh
+
+        som_emi_Anv = som_emi_Anv + em_Anv_cs*cc_Anv_cs*prod_hh
+        som_tcld_Anv = som_tcld_Anv + tcld_Anv_cs*cc_Anv_cs*prod_hh
+        som_pcld_Anv = som_pcld_Anv + pcld_Anv_cs*cc_Anv_cs*prod_hh
+
+        som_emi_hist = som_emi_hist + em_hist_cs*cc_hist_cs*prod_hh
+        som_iwp_hist = som_iwp_hist + iwp_hist_cs*cc_hist_cs*prod_hh
+        som_deltaz_hist = som_deltaz_hist + &
+     & deltaz_hist_cs*cc_hist_cs*prod_hh
+        som_rad_hist = som_rad_hist + rad_hist_cs*cc_hist_cs*prod_hh
+
+        som_pcl_hc = som_pcl_hc + pcld_hc_cs*cc_hc_cs*prod_hh
+        som_tcl_hc = som_tcl_hc + tcld_hc_cs*cc_hc_cs*prod_hh
+
+        som_hc = som_hc + cc_hc_cs*prod_hh
+        som_hist = som_hist + cc_hist_cs*prod_hh
+
+        som_Cb = som_Cb + cc_Cb_cs*prod_hh
+        som_ThCi = som_ThCi + cc_ThCi_cs*prod_hh
+        som_Anv = som_Anv + cc_Anv_cs*prod_hh
+
+
+! Pour test
+        
+        call test_bornes('cc_tot_cs     ',cc_tot_cs,1.,0.)
+        call test_bornes('cc_hc_cs      ',cc_hc_cs,1.,0.)
+        call test_bornes('cc_hist_cs    ',cc_hist_cs,1.,0.)
+        call test_bornes('pcld_hc_cs    ',pcld_hc_cs,1200.,0.)
+        call test_bornes('tcld_hc_cs    ',tcld_hc_cs,1000.,100.)
+        call test_bornes('em_hc_cs      ',em_hc_cs,1000.,0.048)
+
+        test_tot(num_cs) = cc_tot_cs
+        test_hc(num_cs) = cc_hc_cs
+        test_hist(num_cs) = cc_hist_cs
+        test_pcld(num_cs) = pcld_hc_cs
+        test_tcld(num_cs) = tcld_hc_cs
+        test_em(num_cs) = em_hc_cs
+        test_iwp(num_cs) = iwp_hc_cs
+
+        tsom_tot = tsom_tot + cc_tot_cs
+        tsom_hc = tsom_hc + cc_hc_cs
+        tsom_hist = tsom_hist + cc_hist_cs
+
+
+        endif                   ! itop > 0
+
+        ltop = ibot -1
+
+        enddo                   ! fin de la boucle sur la colonne
+
+        N_CS = num_cs
+
+
+! Determination des variables de sortie
+
+        if (N_CS .gt. 0) then   ! if N_CS>0
+
+        cc_tot_mesh = 1. - prod
+
+        cc_hc_mesh = som_hc
+        cc_hist_mesh = som_hist
+
+        cc_Cb_mesh = som_Cb
+        cc_ThCi_mesh = som_ThCi
+        cc_Anv_mesh = som_Anv
+
+        call normal2_undef(pcld_hc_mesh,som_pcl_hc, &
+     & cc_hc_mesh)
+        call normal2_undef(tcld_hc_mesh,som_tcl_hc, &
+     & cc_hc_mesh)
+        call normal2_undef(em_hc_mesh,som_emi_hc, &
+     & cc_hc_mesh)
+        call normal2_undef(iwp_hc_mesh,som_iwp_hc, &
+     & cc_hc_mesh)
+        call normal2_undef(deltaz_hc_mesh,som_deltaz_hc, &
+     & cc_hc_mesh)
+
+        call normal2_undef(em_Cb_mesh,som_emi_Cb, &
+     & cc_Cb_mesh)
+        call normal2_undef(tcld_Cb_mesh,som_tcld_Cb, &
+     & cc_Cb_mesh)
+        call normal2_undef(pcld_Cb_mesh,som_pcld_Cb, &
+     & cc_Cb_mesh)
+
+        call normal2_undef(em_ThCi_mesh,som_emi_ThCi, &
+     & cc_ThCi_mesh)
+        call normal2_undef(tcld_ThCi_mesh,som_tcld_ThCi, &
+     & cc_ThCi_mesh)
+        call normal2_undef(pcld_ThCi_mesh,som_pcld_ThCi, &
+     & cc_ThCi_mesh)
+
+       call normal2_undef(em_Anv_mesh,som_emi_Anv, &
+     & cc_Anv_mesh)
+        call normal2_undef(tcld_Anv_mesh,som_tcld_Anv, &
+     & cc_Anv_mesh)
+        call normal2_undef(pcld_Anv_mesh,som_pcld_Anv, &
+     & cc_Anv_mesh)
+
+
+        call normal2_undef(em_hist_mesh,som_emi_hist, &
+     & cc_hist_mesh)
+        call normal2_undef(iwp_hist_mesh,som_iwp_hist, &
+     & cc_hist_mesh)
+        call normal2_undef(deltaz_hist_mesh,som_deltaz_hist, &
+     & cc_hist_mesh)
+        call normal2_undef(rad_hist_mesh,som_rad_hist, &
+     & cc_hist_mesh)
+
+
+! Tests 
+
+        ! Tests
+
+       if (cc_tot_mesh .gt. tsom_tot .and. &
+     & abs(cc_tot_mesh-tsom_tot) .gt. 1.e-4) then
+        write(*,*) 'cc_tot_mesh > tsom_tot'
+        write(*,*) cc_tot_mesh, tsom_tot
+        STOP
+        endif
+
+        if (cc_tot_mesh .lt. maxval(test_tot(1:N_CS)) .and. &
+     & abs(cc_tot_mesh-maxval(test_tot(1:N_CS))) .gt. 1.e-4) then
+        write(*,*) 'cc_tot_mesh < max'
+        write(*,*) cc_tot_mesh, maxval(test_tot(1:N_CS))
+        STOP
+        endif
+
+        if (cc_hc_mesh .gt. tsom_hc .and. &
+     & abs(cc_hc_mesh-tsom_hc) .gt. 1.e-4) then
+        write(*,*) 'cc_hc_mesh > tsom_hc'
+        write(*,*) cc_hc_mesh, tsom_hc
+        STOP
+        endif
+
+        if (cc_hc_mesh .lt. maxval(test_hc(1:N_CS)) .and. &
+     & abs(cc_hc_mesh-maxval(test_hc(1:N_CS))) .gt. 1.e-4) then
+        write(*,*) 'cc_hc_mesh < max'
+        write(*,*) cc_hc_mesh, maxval(test_hc(1:N_CS))
+        STOP
+        endif
+
+        if (cc_hist_mesh .gt. tsom_hist .and. &
+     & abs(cc_hist_mesh-tsom_hist) .gt. 1.e-4) then
+        write(*,*) 'cc_hist_mesh > tsom_hist'
+        write(*,*) cc_hist_mesh, tsom_hist
+        STOP
+        endif
+
+        if (cc_hist_mesh .lt. 0.) then
+        write(*,*) 'cc_hist_mesh < 0'
+        write(*,*) cc_hist_mesh
+        STOP
+        endif
+
+        if ((pcld_hc_mesh .gt. maxval(test_pcld(1:N_CS)) .or. &
+     & pcld_hc_mesh .lt. minval(test_pcld(1:N_CS))) .and. &
+     & abs(pcld_hc_mesh-maxval(test_pcld(1:N_CS))) .gt. 1. .and. &
+     & maxval(test_pcld(1:N_CS)) .ne. 999. &
+     & .and. minval(test_pcld(1:N_CS)) .ne. 999.) then
+        write(*,*) 'pcld_hc_mesh est faux'
+        write(*,*) pcld_hc_mesh, maxval(test_pcld(1:N_CS)), &
+     & minval(test_pcld(1:N_CS))
+        STOP
+        endif
+
+       if ((tcld_hc_mesh .gt. maxval(test_tcld(1:N_CS)) .or. &
+     & tcld_hc_mesh .lt. minval(test_tcld(1:N_CS))) .and. &
+     & abs(tcld_hc_mesh-maxval(test_tcld(1:N_CS))) .gt. 0.1 .and. &
+     & maxval(test_tcld(1:N_CS)) .ne. 999. &
+     & .and. minval(test_tcld(1:N_CS)) .ne. 999.) then
+        write(*,*) 'tcld_hc_mesh est faux'
+        write(*,*) tcld_hc_mesh, maxval(test_tcld(1:N_CS)), &
+     & minval(test_tcld(1:N_CS))
+        endif
+
+        if ((em_hc_mesh .gt. maxval(test_em(1:N_CS)) .or. &
+     & em_hc_mesh .lt. minval(test_em(1:N_CS))) .and. &
+     & abs(em_hc_mesh-maxval(test_em(1:N_CS))) .gt. 1.e-4 .and. &
+     & minval(test_em(1:N_CS)) .ne. 999. .and. &
+     & maxval(test_em(1:N_CS)) .ne. 999. ) then
+        write(*,*) 'em_hc_mesh est faux'
+        write(*,*) em_hc_mesh, maxval(test_em(1:N_CS)), &
+     & minval(test_em(1:N_CS))
+        STOP
+        endif
+
+        else               ! if N_CS>0
+
+        cc_tot_mesh = 0.
+        cc_hc_mesh = 0.
+        cc_hist_mesh = 0.
+
+        cc_Cb_mesh = 0.
+        cc_ThCi_mesh = 0.
+        cc_Anv_mesh = 0.
+
+        iwp_hc_mesh = undef 
+        deltaz_hc_mesh = undef 
+        em_hc_mesh = undef 
+        iwp_hist_mesh = undef 
+        deltaz_hist_mesh = undef 
+        rad_hist_mesh = undef 
+        em_hist_mesh = undef 
+        pcld_hc_mesh = undef 
+        tcld_hc_mesh = undef 
+
+        pcld_Cb_mesh = undef 
+        tcld_Cb_mesh = undef 
+        em_Cb_mesh = undef 
+
+        pcld_ThCi_mesh = undef 
+        tcld_ThCi_mesh = undef 
+        em_ThCi_mesh = undef 
+
+        pcld_Anv_mesh = undef 
+        tcld_Anv_mesh = undef 
+        em_Anv_mesh = undef 
+
+
+        endif                  ! if N_CS>0
+
+        end subroutine sim_mesh
+
+        subroutine test_bornes(sx,x,bsup,binf)
+
+        real, intent(in) :: x, bsup, binf
+        character*14, intent(in) :: sx
+
+        if (x .gt. bsup .or. x .lt. binf) then
+        write(*,*) sx, 'est faux'
+        write(*,*) sx, x
+        STOP
+        endif
+ 
+        end subroutine test_bornes
+
+        end module m_simu_airs
+
+
+        subroutine simu_airs &
+     & (itap, rneb_airs, temp_airs, cldemi_airs, iwcon0_airs, rad_airs, &
+     & geop_airs, pplay_airs, paprs_airs, &
+     & map_prop_hc,map_prop_hist,&
+     & map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
+     & map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
+     & map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
+     & map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
+     & map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
+     & map_ntot,map_hc,map_hist,&
+     & map_Cb,map_ThCi,map_Anv,alt_tropo )
+
+        USE dimphy
+        USE m_simu_airs
+
+        IMPLICIT NONE
+
+        include "YOMCST.h"
+
+        integer,intent(in) :: itap
+
+        real, dimension(klon,klev), intent(in) :: &
+     & rneb_airs, temp_airs, cldemi_airs, iwcon0_airs, &
+     & rad_airs, geop_airs, pplay_airs, paprs_airs
+
+       real, dimension(klon,klev) :: &
+     & rhodz_airs, rho_airs, iwcon_airs
+
+        real, dimension(klon),intent(out) :: alt_tropo
+
+        real, dimension(klev) :: rneb_1D, temp_1D, &
+     & emis_1D, rad_1D, pres_1D, alt_1D, &
+     & rhodz_1D, dz_1D, iwcon_1D
+
+        integer :: i, j
+
+        real :: cc_tot_mesh, cc_hc_mesh, cc_hist_mesh
+        real :: cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh
+        real :: pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh
+        real :: em_hist_mesh, iwp_hist_mesh
+        real :: deltaz_hc_mesh, deltaz_hist_mesh, rad_hist_mesh
+        real :: pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh
+        real :: pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh
+        real :: pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh
+
+        real, dimension(klon),intent(out) :: map_prop_hc, map_prop_hist
+        real, dimension(klon),intent(out) :: map_emis_hc, map_iwp_hc
+        real, dimension(klon),intent(out) :: map_deltaz_hc, map_pcld_hc
+        real, dimension(klon),intent(out) :: map_tcld_hc
+        real, dimension(klon),intent(out) :: map_emis_Cb,map_pcld_Cb,map_tcld_Cb
+        real, dimension(klon),intent(out) :: &
+     & map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi
+        real, dimension(klon),intent(out) :: &
+     & map_emis_Anv,map_pcld_Anv,map_tcld_Anv
+        real, dimension(klon),intent(out) :: &
+     & map_emis_hist,map_iwp_hist,map_deltaz_hist,&
+     & map_rad_hist
+        real, dimension(klon),intent(out) :: map_ntot,map_hc,map_hist
+        real, dimension(klon),intent(out) :: map_Cb,map_ThCi,map_Anv
+ 
+ 
+        write(*,*) 'simu_airs'
+        write(*,*) 'itap, klon, klev', itap, klon, klev
+        write(*,*) 'RG, RD =', RG, RD
+
+
+! Definition des variables 1D
+
+        do i = 1, klon
+        do j = 1, klev-1
+        rhodz_airs(i,j) = &
+     & (paprs_airs(i,j)-paprs_airs(i,j+1))/RG
+        enddo
+        rhodz_airs(i,klev) = 0.
+        enddo
+
+        do i = 1, klon
+        do j = 1,klev
+        rho_airs(i,j) = &
+     & pplay_airs(i,j)/(temp_airs(i,j)*RD)
+
+        if (rneb_airs(i,j) .gt. 0.001) then
+        iwcon_airs(i,j) = iwcon0_airs(i,j)/rneb_airs(i,j)
+        else
+        iwcon_airs(i,j) = 0.
+        endif
+ 
+        enddo
+        enddo
+
+!=============================================================================
+
+        do i = 1, klon  ! boucle sur les points de grille
+
+!=============================================================================
+        
+        do j = 1,klev
+
+        rneb_1D(j) = rneb_airs(i,j)
+        temp_1D(j) = temp_airs(i,j)
+        emis_1D(j) = cldemi_airs(i,j)
+        iwcon_1D(j) = iwcon_airs(i,j)
+        rad_1D(j) = rad_airs(i,j)
+        pres_1D(j) = pplay_airs(i,j)
+        alt_1D(j) = geop_airs(i,j)/RG
+        rhodz_1D(j) = rhodz_airs(i,j)
+        dz_1D(j) = rhodz_airs(i,j)/rho_airs(i,j)
+
+        enddo
+
+        alt_tropo(i) = &
+     & search_tropopause(pres_1D/100.,temp_1D,alt_1D,klev)
+
+
+! Appel du ss-programme sim_mesh
+
+!        if (itap .eq. 1 ) then
+
+        call sim_mesh(rneb_1D, temp_1D, emis_1D, iwcon_1D, rad_1D, &
+     & pres_1D, dz_1D, rhodz_1D, &
+     & cc_tot_mesh, cc_hc_mesh, cc_hist_mesh, &
+     & pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh, &
+     & deltaz_hc_mesh,&
+     & cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh, &
+     & pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh, &
+     & pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh, &
+     & pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh, &
+     & em_hist_mesh, iwp_hist_mesh, deltaz_hist_mesh, rad_hist_mesh)
+
+         write(*,*) '===================================='
+         write(*,*) 'itap, i:', itap, i 
+         write(*,*) 'cc_tot, cc_hc, cc_hist, pcld_hc, tcld_hc, em_hc, &
+     & iwp_hc, em_hist, iwp_hist ='
+         write(*,*) cc_tot_mesh, cc_hc_mesh, cc_hist_mesh
+         write(*,*) pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh
+         write(*,*)  em_hist_mesh, iwp_hist_mesh
+
+!        endif
+
+! Definition des variables a ecrire dans le fichier de sortie
+
+        call normal2_undef(map_prop_hc(i),cc_hc_mesh, &
+     & cc_tot_mesh)
+        call normal2_undef(map_prop_hist(i),cc_hist_mesh, &
+     & cc_tot_mesh)
+
+       map_emis_hc(i) = em_hc_mesh
+       map_iwp_hc(i) = iwp_hc_mesh
+       map_deltaz_hc(i) = deltaz_hc_mesh
+       map_pcld_hc(i) = pcld_hc_mesh
+       map_tcld_hc(i) = tcld_hc_mesh
+
+       map_emis_Cb(i) = em_Cb_mesh
+       map_pcld_Cb(i) = pcld_Cb_mesh
+       map_tcld_Cb(i) = tcld_Cb_mesh
+
+       map_emis_ThCi(i) = em_ThCi_mesh
+       map_pcld_ThCi(i) = pcld_ThCi_mesh
+       map_tcld_ThCi(i) = tcld_ThCi_mesh
+
+       map_emis_Anv(i) = em_Anv_mesh
+       map_pcld_Anv(i) = pcld_Anv_mesh
+       map_tcld_Anv(i) = tcld_Anv_mesh
+
+       map_emis_hist(i) = em_hist_mesh
+       map_iwp_hist(i) = iwp_hist_mesh
+       map_deltaz_hist(i) = deltaz_hist_mesh
+       map_rad_hist(i) = rad_hist_mesh
+
+       map_ntot(i) = cc_tot_mesh
+       map_hc(i) = cc_hc_mesh
+       map_hist(i) = cc_hist_mesh
+
+       map_Cb(i) = cc_Cb_mesh
+       map_ThCi(i) = cc_ThCi_mesh
+       map_Anv(i) = cc_Anv_mesh
+
+
+        enddo         ! fin boucle sur les points de grille
+
+        
+
+        end subroutine simu_airs
+
