Index: /LMDZ5/branches/testing/DefLists/README
===================================================================
--- /LMDZ5/branches/testing/DefLists/README	(revision 2669)
+++ /LMDZ5/branches/testing/DefLists/README	(revision 2669)
@@ -0,0 +1,12 @@
+physiq.def                 Couples model version
+
+physiq.def_L39_AR40.0      IPSL-CM5A
+physiq.def_L39_NPv3.1      IPSL-CM5B
+physiq.def_L39_NPv3.2      Stabilized version
+physiq.def_L59_NPv4.12     First L59 version of the new new physics
+physiq.def_L79_NPv5.17h    (5.17h=5.3) IPSL-CM6.0.1
+physiq.def_NPv5.4          IPSL-CM6.0.2
+physiq.def_NPv5.5          IPSL-CM6.0.3
+physiq.def_NPv5.65         IPSL-CM6.0.4
+physiq.def_NPv5.70         IPSL-CM6.0.5
+physiq.def_NPv5.67         IPSL-CM6.0.5 / THC1 (test de sensibilité)
Index: /LMDZ5/branches/testing/DefLists/context_lmdz.xml
===================================================================
--- /LMDZ5/branches/testing/DefLists/context_lmdz.xml	(revision 2668)
+++ /LMDZ5/branches/testing/DefLists/context_lmdz.xml	(revision 2669)
@@ -31,5 +31,10 @@
     <grid_group id="vertical" axis_ref="presnivs" />
   </grid_> -->
-  
+
+<!-- Define Scalar grid for GHG, orbital parameters and solar constants -->
+  <grid_definition>
+    <grid id="grid_scalar" />
+    </grid_definition>
+ 
   <!-- Define groups of vertical axes -->
   <axis_definition>
Index: /LMDZ5/branches/testing/DefLists/field_def_lmdz.xml
===================================================================
--- /LMDZ5/branches/testing/DefLists/field_def_lmdz.xml	(revision 2668)
+++ /LMDZ5/branches/testing/DefLists/field_def_lmdz.xml	(revision 2669)
@@ -867,4 +867,24 @@
     </field_group>
 
+
+   <field_group id="field_scalar" operation="instant" freq_op="30d" grid_ref="grid_scalar">
+      <field id="R_ecc" long_name="R_ecc" unit="-" />
+      <field id="R_peri" long_name="R_peri" unit="-" />
+      <field id="R_incl" long_name="R_incl" unit="deg" />
+      <field id="solaire" long_name="solaire" unit="W/m2" />
+
+      <field id="rsun1" long_name="Fraction constante solaire bande 1" unit="W/m2" />
+      <field id="rsun2" long_name="Fraction constante solaire bande 2" unit="W/m2" />
+      <field id="rsun3" long_name="Fraction constante solaire bande 3" unit="W/m2" />
+      <field id="rsun4" long_name="Fraction constante solaire bande 4" unit="W/m2" />
+      <field id="rsun5" long_name="Fraction constante solaire bande 5" unit="W/m2" />
+      <field id="rsun6" long_name="Fraction constante solaire bande 6" unit="W/m2" />
+
+      <field id="co2_ppm" long_name="co2_ppm" unit="ppm" />
+      <field id="CH4_ppb" long_name="CH4_ppb" unit="ppb" />
+      <field id="N2O_ppb" long_name="N2O_ppb" unit="ppb" />
+      <field id="CFC11_ppt" long_name="CFC11_ppt" unit="ppt" />
+      <field id="CFC12_ppt" long_name="CFC12_ppt" unit="ppt" />
+   </field_group>
+
 </field_definition>
-
Index: /LMDZ5/branches/testing/DefLists/file_def_histmth_lmdz.xml
===================================================================
--- /LMDZ5/branches/testing/DefLists/file_def_histmth_lmdz.xml	(revision 2668)
+++ /LMDZ5/branches/testing/DefLists/file_def_histmth_lmdz.xml	(revision 2669)
@@ -10,4 +10,26 @@
 <!--            </field_group>  -->
         
+
+            <field_group operation="instant" freq_op="30d">
+                <field field_ref="R_ecc" level="1" name="R_ecc" />
+                <field field_ref="R_peri" level="1" name="R_peri" />
+                <field field_ref="R_incl" level="1" name="R_incl" />
+                <field field_ref="solaire" level="1" name="solaire" />
+
+                <field field_ref="rsun1" level="1" name="rsun1" />
+                <field field_ref="rsun2" level="1" name="rsun2" />
+                <field field_ref="rsun3" level="1" name="rsun3" />
+                <field field_ref="rsun4" level="1" name="rsun4" />
+                <field field_ref="rsun5" level="1" name="rsun5" />
+                <field field_ref="rsun6" level="1" name="rsun6" />
+
+                <field field_ref="co2_ppm" level="1" name="co2_ppm" />
+                <field field_ref="CH4_ppb" level="1" name="CH4_ppb" />
+                <field field_ref="N2O_ppb" level="1" name="N2O_ppb" />
+                <field field_ref="CFC11_ppt" level="1" name="CFC11_ppt" />
+                <field field_ref="CFC12_ppt" level="1" name="CFC12_ppt" />
+
+            </field_group>
+
             <!-- VARS 2D -->
             <field_group operation="average" freq_op="1ts">
Index: /LMDZ5/branches/testing/DefLists/gcm.def_144x142
===================================================================
--- /LMDZ5/branches/testing/DefLists/gcm.def_144x142	(revision 2668)
+++ /LMDZ5/branches/testing/DefLists/gcm.def_144x142	(revision 2669)
@@ -1,5 +1,5 @@
-## $Id: gcm.def 1403 2010-07-01 09:02:53Z fairhead $
+## $Id$
 ## nombre de pas par jour (multiple de iperiod) ( ici pour  dt = 1 min )      
-day_step=960
+day_step=720
 ## periode pour le pas Matsuno (en pas)
 iperiod=5
@@ -13,10 +13,9 @@
 niterh=2
 ## temps de dissipation des plus petites long.d ondes pour u,v (gradiv)  
-tetagdiv=1800.
+tetagdiv=3600.
 ## temps de dissipation des plus petites long.d ondes pour u,v(nxgradrot)
-#tetagrot=5400.
 tetagrot=7200.
 ## temps de dissipation des plus petites long.d ondes pour  h ( divgrad) 
-tetatemp=5400.
+tetatemp=7200.
 ## coefficient pour gamdissip                                            
 coefdis=0.
@@ -33,8 +32,6 @@
 read_start=y
 ## periode de la physique (en pas dynamiques, n'a de sens que si iflag_phys=1)                                       
+## iphysiq est modifie pour la nouvelle physique par le lmdz.driver pour les autres cas le valeur par defaut est utilise
 iphysiq=_AUTO_ : DEFAULT = 5
-#
-#
-### Parametres du zoom
 ## longitude en degres du centre du zoom                                 
 clon=0.
@@ -57,2 +54,4 @@
 ##  Fonction  f(y) avec y = Sin(latit.) si = .true. , sinon y = latit.         
 ysinus=y
+##
+ngroup=4
Index: /LMDZ5/branches/testing/DefLists/iodef_dev.xml
===================================================================
--- /LMDZ5/branches/testing/DefLists/iodef_dev.xml	(revision 2668)
+++ /LMDZ5/branches/testing/DefLists/iodef_dev.xml	(revision 2669)
@@ -1,5 +1,8 @@
 <?xml version="1.0"?>
 <simulation> 
-    <context id="LMDZ" calendar_type="D360" start_date="1980-01-01 00:00:00">
+    <context id="LMDZ">
+        
+        <!-- <calendar type="D360" start_date="1980-01-01 00:00:00" /> -->
+        
         <!-- Definition of model variables -->
         <field_definition prec="4" 
@@ -16,5 +19,5 @@
             <!-- 3D variables -->
             <field_group id="fields_3D"
-                         domain_ref="dom_glo">
+                         domain_ref="dom_glo" axis_ref="presnivs">
                 <field id="temperature"
                        long_name="Atmospheric temperature"
@@ -39,4 +42,5 @@
                       name="Xhistins"
                       output_freq="12ts" 
+                      type="one_file"
                       enabled=".TRUE."> <!-- NB: output_freq in physics ts-->
                     
@@ -46,6 +50,4 @@
                         <field field_ref="ps" />
                     </field_group>
-                    <!-- <field_group field_group_ref="fields_2D"
-                                 operation="instant" /> doesn't work ?!? -->
 
                     <!-- VARS 3D -->
@@ -58,6 +60,4 @@
                         <field field_ref="v" />
                     </field_group>
-                    <!-- <field_group field_group_ref="fields_3D"
-                                 operation="instant" /> doesn't work ?!? -->
                 </file>
 
@@ -65,4 +65,5 @@
                       name="Xdiurnalave"
                       output_freq="1d"
+                      type="one_file"
                       enabled=".true.">
                     <!-- VARS 2D -->
@@ -93,10 +94,4 @@
              wxios_add_vaxis -->
         <axis_definition>
-            <!--
-            <axis_group id="presnivs" 
-                        standard_name="Pseudo-pressure of model vertical levels" 
-                        unit="Pa">
-            </axis_group>
-            -->
             <axis id="presnivs" 
                         standard_name="Pseudo-pressure of model vertical levels" 
@@ -114,5 +109,5 @@
 
             <variable_group id="parameters" >
-                <variable id="using_server" type="boolean">true</variable>
+                <variable id="using_server" type="bool">false</variable>
                 <variable id="info_level" type="int">10</variable>
             </variable_group>
Index: DZ5/branches/testing/DefLists/physiq.def_L79_NPv5.4
===================================================================
--- /LMDZ5/branches/testing/DefLists/physiq.def_L79_NPv5.4	(revision 2668)
+++ 	(revision )
@@ -1,229 +1,0 @@
-##############################################################################RAMETRES ET CLES DE LA PHYSIQUE LMDZ
-###############################################################################
-## D: Valeur par default
-##
-##Parametres Orographiques et cdrags
-##**********************************
-###   Si=.T. ,  lecture du fichier limit avec la bonne annee 
-ok_limitvrai=n
-iflag_albedo=1
-#Cdrags 
-f_cdrag_stable=1.
-f_cdrag_ter=1.
-f_cdrag_oce=0.7
-f_ri_cd_min=0.6
-#cdmmax
-cdmmax=2.5E-3
-#cdhmax
-cdhmax=2.0E-3
-##  Orodr  ou  non   pour l orographie
-ok_orodr=y
-##  Orolf  ou  non   pour l orographie              
-ok_orolf=y
-#Rugoro
-f_rugoro=0.
-#
-### Random gravity waves
-#Activitate the stochastic parameterization of convective GWs (should be equal to ok_strato)
-ok_gwd_rando=y
-# Scale the momentum fluxes due to convective GWs
-gwd_rando_ruwmax=2.5
-#Scale the saturated momentum flux of  convective GWs 
-gwd_rando_sat=0.2
-
-ok_strato=y
-ok_hines=n
-sso_gkdrag=0.2
-sso_gfrcri=0.70
-sso_gkwake=0.40
-#
-# Rayonnement
-#*********************
-# activation du nouveau code de rayonnement RRTM
-# 0 : Ancien code et 1 : RRTM (D=0)
-iflag_rrtm=0
-#NSW=2
-#
-# Parametres nuages
-#******************
-#seuils de la precipitation des nuages strateformes (D: 2.6e-4 2.6e-4)
-cld_lc_lsc=0.0001
-cld_lc_con=0.0001
-#constante de temps pour eleminer eau lsc et convective (D: 3600. 3600.)
-cld_tau_lsc=900.
-cld_tau_con=900.
-#facteurs correctifs sur la vitesse de chute des cristaux de glace (D: 1 1)
-ffallv_lsc=1.
-ffallv_con=1.
-#coefficient sur la reevaporation de la pluie (D: 2.e-5 n)
-coef_eva=2e-5
-reevap_ice=y
-#calcul des proprietes des nuages convectifs (D:1 0.375 1.e-4) 
-iflag_cldcon=6
-fact_cldcon=1.
-facttemps=0.
-#calcul eau condensee et fraction nuageuse a partir des PDF (D:0, 0:version avec ratqs sinon nouvelles PDFs)  
-iflag_pdf=1
-#calcul epaisseur optique et emmissivite des nuages (D: y 1 0.01 0.3)
-ok_newmicro=y
-iflag_ratqs=4
-ratqsbas=0.002
-ratqshaut=0.312
-#rayon effectif pour les nuages de glace (D:35)
-rad_froid=35
-#rayons effectifs pour les nuages eau liq (tailles des gouttes d eau nuageuse, D: 13 9)
-rad_chau1=12
-rad_chau2=11
-##  Choix ou non  de  New oliq               
-new_oliq=y
-#Rayons min et max effectifs des particules des nuages
-rei_min=16.
-rei_max=61.29
-#
-# Flags Convection 
-#*****************
-#Convection scheme switch (D:2, 1:LMD, 2:Tiedtke, 3:KE New Physics, 30:KE AR4)
-iflag_con=3
-#output level of energy conservation diagnostics
-if_ebil=0
-#maximum efficiency of cld water->precipitation conversion (D: 0.993)
-epmax=0.9995
-#dry convective adjustment at beginning of Emanuel convective scheme (D: n 0)
-ok_adj_ema=n
-iflag_clw=0
-#Convective closure (D:1, 0:AR4, 2:ALE and ALP)
-iflag_clos=2
-#Convective entrainment mixing law (D:1, 0:AR4=flat PDF, 1=PDF)
-iflag_mix=1
-#weights of the bell shaped and flat PDF (used only if iflag_mix=1) (D: 1 0)
-qqa1=1.
-qqa2=0.
-#Parameters of the bell shaped PDF; Fmax=center (D= 0.65), Gammas=width (D=0.05)
-cvl_corr=1.0
-#Multiplicative factor of Emanuel convective scheme precipitation (D: 1.)
-Fmax=0.65
-#
-# Flags Couche limite
-#********************
-#shema de CL des thermiques (D:0, 0:ajustement sec,=>1:versions thermiques) 
-iflag_thermals=18
-#pas de temps du spliting pour les thermiques
-nsplit_thermals=1
-tau_thermals=0.
-iflag_thermals_ed=8
-fact_thermals_ed_dz=0.1
-iflag_thermals_optflux=0
-#shema de la couche de surface (D:1, 1:LMD, 8:Mellor-Yamada) 
-iflag_pbl=11
-#diffusion turbulente
-ksta_ter=1.e-7
-ksta=1.e-10
-#ok_kzmin : calcul Kzmin dans la CL de surface
-ok_kzmin=y
-#couplage ou non avec la convection
-iflag_coupl=5 (0:AR4 et 1:nouvelle physique)
-seuil_inversion=-0.08
-#rendement sur la puissance fournie par les thermiques a la convection
-alp_bl_k=0.5
-#
-## temperature dependency of qsat for iterative procedure
-## temporary flag to switch between the old (CMIP5, D=0) and new versions (=1) versions 
-iflag_fisrtilp_qsat=4 
-## Reactivation du calcul d'un zmax continu pour les thermiques
-## reactivation d'une fermeture humide des thermiques  iflag_thermals_closure=2 (D=1)
-iflag_thermals_closure=1
-#
-# Flags wakes
-#************
-#Wake scheme switch (D:0, 0:AR4, 1:New Physics)
-iflag_wake=1
-#Constant value to be added to Alp (D: 0.)
-alp_offset=0.
-
-### nrlmd le 10/04/2012: Declenchement ###
-#Déclelenchement par couche-limite (0: deterministe, 1: stochastique)
-iflag_trig_bl=1
-#Surface critique au LCL
-s_trig=1.2e7
-#Temps de decorrelation entre scenes de cumulus
-tau_trig_shallow=1200
-#Temps de decorrelation entre scenes de cumulonimbus
-tau_trig_deep=1200
-# Pour la stabilisation
-random_notrig_max=0.99
-ok_adjwk=y
-tmax_fonte_cv=288.15
-
-
-### nrlmd le 10/04/2012: Fermeture ###
-#Fermeture par couche-limite (0: deterministe, 1: statistique, 2: statistique+convergence grande échelle)
-iflag_clos_bl=1
-### fin nrlmd le 10/04/2012
-
-## Modification de la plage de temperature pour la phase mixte liquide/glace des nuages
-## Controle par les t_glace_min/max, exposant_glace, iflag_t_glace (D=0)
-t_glace_min=233.13
-t_glace_max=273.13
-exposant_glace=1.
-iflag_t_glace=1
-iflag_ice_thermo=1
-
-#
-## The source terms of the TKE prognostic equation are diagnosed from tendencies (du, dv, dT) 
-## associated with subrgrid scaleamotions and treated as an additional heat source.
-## Controled by a key, iflag_ener_conserv (D=-1):
-## 0 : no conservation, -1 : old adhoc correction for kinetic E only (used for CMIP5) 
-## 1 : conservation, 101 : conversion from kinetic to heat only, 110 conversion from potential to heat only
-iflag_ener_conserv=1
-
-#Activation of corrections for water conservation (D: n). 
-##The default is n for compatibility with earlier versions.
-##generally it should be y.
-ok_conserv_q=y 
-
-### distribution sous-maille sur la verticale 
-iflag_cloudth_vert=1
-
-### Convection
-#altitude, in hPa, above LCL at which buoybase is computed, def = original: -40
-dpbase=-40.
-#altitude, in hPa, above LCL below which Ep=0, def = original: 150, used only if flag_epKEorig=0
-pbcrit=150.0
-#altitude, in hPa, above LCL above which Ep=Epmax, def = original: 500, used only if flag_epKEorig=0
-ptcrit=500.0
-#reference fractional area of precipitating downdraughts, def = original: 0.01
-sigdz=0.01
-#unsaturated fraction of precipitating downdraughts, def = original: 0.15
-spfac=0.15
-#relaxation time of cloud base mass flux, def = original: 8000
-tau=8000.
-#flag for wb (= vert velocity at LFC); 0->wb=wbmax, 1->wb=f(plfc) bounded, 2->wb=f(plfc) linear, def = 1
-flag_wb=1
-#parameter used in the computation of wb, def = 6.
-wbmax=6.
-#
-ok_convstop=.False.
-#
-tau_stop=15000.
-#
-ok_intermittent=.False.
-#shedding coefficient (used when iflag_mix_adiab=1)
-coef_peel=0.25
-#
-### Clouds
-#Conversion cld water->precip: Ep=f(pressure) if 0; Ep=f(temperature) if 1, def=1
-flag_epKEorig=1
-#max condensed water concentration in kg/kg, def=0.0003, original: 0.0011, used only if flag_epKEorig=1
-elcrit=0.0003
-#temperature in degre Celsius, at which Ep=Epmax, def=-55., original: -55, used only if flag_epKEorig=1
-tlcrit=-55.0
-
-### Wakes
-#coefficient in Wape->C_star formula: C_star=stark*sqrt(2*Wape), def: 0.33
-stark=0.33
-#lifting efficiency of wakes; Alp=Alpk*Fip, Fip=power provided by wakes, def: 0.25
-alpk=0.25
-#wake density = number of wake centers per m2, def: 8.E-12
-wdens_ref=8.E-12
-#multiplicative factor of the damping by gravity waves, def: 4.
-coefgw=4.
Index: DZ5/branches/testing/DefLists/physiq.def_L79_NPv5.5
===================================================================
--- /LMDZ5/branches/testing/DefLists/physiq.def_L79_NPv5.5	(revision 2668)
+++ 	(revision )
@@ -1,236 +1,0 @@
-##############################################################################
-# PARAMETRES ET CLES DE LA PHYSIQUE LMDZ
-##############################################################################
-# D: Valeur par default
-#
-#Parametres Orographiques et cdrags
-#**********************************
-##   Si=.T. ,  lecture du fichier limit avec la bonne annee 
-ok_limitvrai=n
-iflag_albedo=1
-#Cdrags 
-f_cdrag_stable=1.
-f_cdrag_ter=1.
-f_cdrag_oce=0.7
-f_ri_cd_min=0.6
-#cdmmax
-cdmmax=2.5E-3
-#cdhmax
-cdhmax=2.0E-3
-##  Orodr  ou  non   pour l orographie
-ok_orodr=y
-##  Orolf  ou  non   pour l orographie              
-ok_orolf=y
-#Rugoro
-f_rugoro=0.
-#
-### Random gravity waves
-#Activitate the stochastic parameterization of convective GWs (should be equal to ok_strato)
-ok_gwd_rando=y
-# Scale the momentum fluxes due to convective GWs
-gwd_rando_ruwmax=2.5
-#Scale the saturated momentum flux of  convective GWs 
-gwd_rando_sat=0.2
-
-ok_strato=y
-ok_hines=n
-sso_gkdrag=0.2
-sso_gfrcri=0.70
-sso_gkwake=0.40
-#
-# Rayonnement
-#*********************
-# activation du nouveau code de rayonnement RRTM
-# 0 : Ancien code et 1 : RRTM (D=0)
-iflag_rrtm=1
-NSW=6
-#
-# Parametres nuages
-#******************
-#seuils de la precipitation des nuages strateformes (D: 2.6e-4 2.6e-4)
-cld_lc_lsc=0.00022
-cld_lc_con=0.00022
-#constante de temps pour eleminer eau lsc et convective (D: 3600. 3600.)
-cld_tau_lsc=1200.
-cld_tau_con=1200.
-#facteurs correctifs sur la vitesse de chute des cristaux de glace (D: 1 1)
-ffallv_lsc=0.67
-ffallv_con=0.67
-#coefficient sur la reevaporation de la pluie (D: 2.e-5 n)
-coef_eva=2e-5
-reevap_ice=y
-#calcul des proprietes des nuages convectifs (D:1 0.375 1.e-4) 
-iflag_cldcon=6
-fact_cldcon=1.
-facttemps=0.
-#calcul eau condensee et fraction nuageuse a partir des PDF (D:0, 0:version avec ratqs sinon nouvelles PDFs)  
-iflag_pdf=1
-#calcul epaisseur optique et emmissivite des nuages (D: y 1 0.01 0.3)
-ok_newmicro=y
-iflag_ratqs=4
-ratqsbas=0.002
-ratqshaut=0.312
-#rayon effectif pour les nuages de glace (D:35)
-rad_froid=35
-#rayons effectifs pour les nuages eau liq (tailles des gouttes d eau nuageuse, D: 13 9)
-rad_chau1=12
-rad_chau2=11
-##  Choix ou non  de  New oliq               
-new_oliq=y
-#Rayons min et max effectifs des particules des nuages
-rei_min=16.
-rei_max=61.29
-#
-# Flags Convection 
-#*****************
-#Convection scheme switch (D:2, 1:LMD, 2:Tiedtke, 3:KE New Physics, 30:KE AR4)
-iflag_con=3
-#output level of energy conservation diagnostics
-if_ebil=0
-#maximum efficiency of cld water->precipitation conversion (D: 0.993)
-epmax=0.999
-#dry convective adjustment at beginning of Emanuel convective scheme (D: n 0)
-ok_adj_ema=n
-iflag_clw=0
-#Convective closure (D:1, 0:AR4, 2:ALE and ALP)
-iflag_clos=2
-#Convective entrainment mixing law (D:1, 0:AR4=flat PDF, 1=PDF)
-iflag_mix=1
-#weights of the bell shaped and flat PDF (used only if iflag_mix=1) (D: 1 0)
-qqa1=1.
-qqa2=0.
-#Parameters of the bell shaped PDF; Fmax=center (D= 0.65), Gammas=width (D=0.05)
-cvl_corr=1.0
-#Multiplicative factor of Emanuel convective scheme precipitation (D: 1.)
-Fmax=0.65
-#
-# Flags Couche limite
-#********************
-#shema de CL des thermiques (D:0, 0:ajustement sec,=>1:versions thermiques) 
-iflag_thermals=18
-#pas de temps du spliting pour les thermiques
-nsplit_thermals=1
-tau_thermals=0.
-iflag_thermals_ed=8
-fact_thermals_ed_dz=0.1
-iflag_thermals_optflux=0
-#shema de la couche de surface (D:1, 1:LMD, 8:Mellor-Yamada) 
-iflag_pbl=11
-#diffusion turbulente
-ksta_ter=1.e-7
-ksta=1.e-10
-#ok_kzmin : calcul Kzmin dans la CL de surface
-ok_kzmin=y
-#couplage ou non avec la convection
-iflag_coupl=5 (0:AR4 et 1:nouvelle physique)
-seuil_inversion=-0.08
-#rendement sur la puissance fournie par les thermiques a la convection
-alp_bl_k=0.5
-#
-## temperature dependency of qsat for iterative procedure
-## temporary flag to switch between the old (CMIP5, D=0) and new versions (=1) versions 
-iflag_fisrtilp_qsat=4 
-## Reactivation du calcul d'un zmax continu pour les thermiques
-## reactivation d'une fermeture humide des thermiques  iflag_thermals_closure=2 (D=1)
-iflag_thermals_closure=1
-#
-# Flags wakes
-#************
-#Wake scheme switch (D:0, 0:AR4, 1:New Physics)
-iflag_wake=1
-#Constant value to be added to Alp (D: 0.)
-alp_offset=0.
-
-### nrlmd le 10/04/2012: Declenchement ###
-#Déclelenchement par couche-limite (0: deterministe, 1: stochastique)
-iflag_trig_bl=1
-#Surface critique au LCL
-s_trig=1.2e7
-#Temps de decorrelation entre scenes de cumulus
-tau_trig_shallow=1200
-#Temps de decorrelation entre scenes de cumulonimbus
-tau_trig_deep=1200
-# Pour la stabilisation
-random_notrig_max=0.99
-ok_adjwk=y
-tmax_fonte_cv=278.15
-
-
-### nrlmd le 10/04/2012: Fermeture ###
-#Fermeture par couche-limite (0: deterministe, 1: statistique, 2: statistique+convergence grande échelle)
-iflag_clos_bl=1
-### fin nrlmd le 10/04/2012
-
-## Modification de la plage de temperature pour la phase mixte liquide/glace des nuages
-## Controle par les t_glace_min/max, exposant_glace, iflag_t_glace (D=0)
-t_glace_min=233.13
-t_glace_max=273.13
-exposant_glace=1.
-iflag_t_glace=1
-iflag_ice_thermo=1
-
-#
-## The source terms of the TKE prognostic equation are diagnosed from tendencies (du, dv, dT) 
-## associated with subrgrid scaleamotions and treated as an additional heat source.
-## Controled by a key, iflag_ener_conserv (D=-1):
-## 0 : no conservation, -1 : old adhoc correction for kinetic E only (used for CMIP5) 
-## 1 : conservation, 101 : conversion from kinetic to heat only, 110 conversion from potential to heat only
-iflag_ener_conserv=1
-
-#Activation of corrections for water conservation (D: n). 
-##The default is n for compatibility with earlier versions.
-##generally it should be y.
-ok_conserv_q=y 
-
-iflag_cloudth_vert=1
-
-#
-lmixmin=10.
-
-#Bergeron
-iflag_bergeron=1
-
-
-### Convection
-#altitude, in hPa, above LCL at which buoybase is computed, def = original: -40
-dpbase=-40.
-#altitude, in hPa, above LCL below which Ep=0, def = original: 150, used only if flag_epKEorig=0
-pbcrit=150.0
-#altitude, in hPa, above LCL above which Ep=Epmax, def = original: 500, used only if flag_epKEorig=0
-ptcrit=500.0
-#reference fractional area of precipitating downdraughts, def = original: 0.01
-sigdz=0.01
-#unsaturated fraction of precipitating downdraughts, def = original: 0.15
-spfac=0.15
-#relaxation time of cloud base mass flux, def = original: 8000
-tau=8000.
-#flag for wb (= vert velocity at LFC); 0->wb=wbmax, 1->wb=f(plfc) bounded, 2->wb=f(plfc) linear, def = 1
-flag_wb=1
-#parameter used in the computation of wb, def = 6.
-wbmax=6.
-#
-ok_convstop=.False.
-#
-tau_stop=15000.
-#
-ok_intermittent=.False.
-#shedding coefficient (used when iflag_mix_adiab=1)
-coef_peel=0.25
-#
-### Clouds
-#Conversion cld water->precip: Ep=f(pressure) if 0; Ep=f(temperature) if 1, def=1
-flag_epKEorig=1
-#max condensed water concentration in kg/kg, def=0.0003, original: 0.0011, used only if flag_epKEorig=1
-elcrit=0.0003
-#temperature in degre Celsius, at which Ep=Epmax, def=-55., original: -55, used only if flag_epKEorig=1
-tlcrit=-55.0
-
-### Wakes
-#coefficient in Wape->C_star formula: C_star=stark*sqrt(2*Wape), def: 0.33
-stark=0.33
-#lifting efficiency of wakes; Alp=Alpk*Fip, Fip=power provided by wakes, def: 0.25
-alpk=0.25
-#wake density = number of wake centers per m2, def: 8.E-12
-wdens_ref=8.E-12
-#multiplicative factor of the damping by gravity waves, def: 4.
-coefgw=4.
Index: /LMDZ5/branches/testing/DefLists/physiq.def_NPv5.65
===================================================================
--- /LMDZ5/branches/testing/DefLists/physiq.def_NPv5.65	(revision 2668)
+++ /LMDZ5/branches/testing/DefLists/physiq.def_NPv5.65	(revision 2669)
@@ -36,6 +36,4 @@
 gwd_rando_sat=0.2
 
-ok_strato=y
-ok_hines=n
 sso_gkdrag=0.2
 sso_gfrcri=0.70
Index: /LMDZ5/branches/testing/DefLists/physiq.def_NPv5.67
===================================================================
--- /LMDZ5/branches/testing/DefLists/physiq.def_NPv5.67	(revision 2668)
+++ /LMDZ5/branches/testing/DefLists/physiq.def_NPv5.67	(revision 2669)
@@ -36,6 +36,4 @@
 gwd_rando_sat=0.2
 
-ok_strato=y
-ok_hines=n
 sso_gkdrag=0.2
 sso_gfrcri=0.70
Index: /LMDZ5/branches/testing/DefLists/physiq.def_NPv5.70
===================================================================
--- /LMDZ5/branches/testing/DefLists/physiq.def_NPv5.70	(revision 2669)
+++ /LMDZ5/branches/testing/DefLists/physiq.def_NPv5.70	(revision 2669)
@@ -0,0 +1,240 @@
+#
+# $Id: physiq.def_NPv5.70 2929 2016-07-26 06:49:39Z jgipsl $
+#
+##############################################################################
+# PARAMETRES ET CLES DE LA PHYSIQUE LMDZ
+##############################################################################
+# D: Valeur par default
+#
+#Parametres Orographiques et cdrags
+#**********************************
+##   Si=.T. ,  lecture du fichier limit avec la bonne annee 
+ok_limitvrai=n
+iflag_albedo=1
+#Cdrags 
+f_cdrag_stable=1.
+f_cdrag_ter=1.
+f_cdrag_oce=0.6
+f_ri_cd_min=0.1
+#cdmmax
+cdmmax=2.5E-3
+#cdhmax
+cdhmax=2.0E-3
+##  Orodr  ou  non   pour l orographie
+ok_orodr=y
+##  Orolf  ou  non   pour l orographie              
+ok_orolf=y
+#Rugoro
+f_rugoro=0.
+#
+### Random gravity waves
+#Activitate the stochastic parameterization of convective GWs (should be equal to ok_strato)
+ok_gwd_rando=y
+# Scale the momentum fluxes due to convective GWs
+gwd_rando_ruwmax=2.5
+#Scale the saturated momentum flux of  convective GWs 
+gwd_rando_sat=0.2
+
+sso_gkdrag=0.2
+sso_gfrcri=0.70
+sso_gkwake=0.40
+#
+# Rayonnement
+#*********************
+# activation du nouveau code de rayonnement RRTM
+# 0 : Ancien code et 1 : RRTM (D=0)
+iflag_rrtm=1
+NSW=6
+#
+# Parametres nuages
+#******************
+#seuils de la precipitation des nuages strateformes (D: 2.6e-4 2.6e-4)
+cld_lc_lsc=0.0002
+cld_lc_con=0.0002
+#constante de temps pour eleminer eau lsc et convective (D: 3600. 3600.)
+cld_tau_lsc=900.
+cld_tau_con=900.
+#facteurs correctifs sur la vitesse de chute des cristaux de glace (D: 1 1)
+ffallv_lsc=0.5
+ffallv_con=0.5
+#coefficient sur la reevaporation de la pluie (D: 2.e-5 n)
+coef_eva=2e-4
+reevap_ice=y
+#calcul des proprietes des nuages convectifs (D:1 0.375 1.e-4) 
+iflag_cldcon=6
+fact_cldcon=1.
+facttemps=0.
+#calcul eau condensee et fraction nuageuse a partir des PDF (D:0, 0:version avec ratqs sinon nouvelles PDFs)  
+iflag_pdf=1
+#calcul epaisseur optique et emmissivite des nuages (D: y 1 0.01 0.3)
+ok_newmicro=y
+iflag_ratqs=4
+ratqsp0=37500.
+ratqsdp=20000.
+ratqsbas=0.002
+ratqshaut=0.4
+#rayon effectif pour les nuages de glace (D:35)
+rad_froid=35
+#rayons effectifs pour les nuages eau liq (tailles des gouttes d eau nuageuse, D: 13 9)
+rad_chau1=12
+rad_chau2=11
+##  Choix ou non  de  New oliq               
+new_oliq=y
+#Rayons min et max effectifs des particules des nuages
+rei_min=16.
+rei_max=61.29
+#
+# Flags Convection 
+#*****************
+#Convection scheme switch (D:2, 1:LMD, 2:Tiedtke, 3:KE New Physics, 30:KE AR4)
+iflag_con=3
+#output level of energy conservation diagnostics
+if_ebil=0
+#maximum efficiency of cld water->precipitation conversion (D: 0.993)
+epmax=0.999
+#dry convective adjustment at beginning of Emanuel convective scheme (D: n 0)
+ok_adj_ema=n
+iflag_clw=0
+#Convective closure (D:1, 0:AR4, 2:ALE and ALP)
+iflag_clos=2
+#Convective entrainment mixing law (D:1, 0:AR4=flat PDF, 1=PDF)
+iflag_mix=1
+#weights of the bell shaped and flat PDF (used only if iflag_mix=1) (D: 1 0)
+qqa1=1.
+qqa2=0.
+#Parameters of the bell shaped PDF; Fmax=center (D= 0.65), Gammas=width (D=0.05)
+cvl_corr=1.0
+#Multiplicative factor of Emanuel convective scheme precipitation (D: 1.)
+Fmax=0.65
+#
+# Flags Couche limite
+#********************
+#shema de CL des thermiques (D:0, 0:ajustement sec,=>1:versions thermiques) 
+iflag_thermals=18
+#pas de temps du spliting pour les thermiques
+nsplit_thermals=1
+tau_thermals=0.
+iflag_thermals_ed=8
+fact_thermals_ed_dz=0.1
+iflag_thermals_optflux=0
+#shema de la couche de surface (D:1, 1:LMD, 8:Mellor-Yamada) 
+iflag_pbl=11
+#diffusion turbulente
+ksta_ter=1.e-7
+ksta=1.e-10
+#ok_kzmin : calcul Kzmin dans la CL de surface
+ok_kzmin=y
+#couplage ou non avec la convection
+iflag_coupl=5 (0:AR4 et 1:nouvelle physique)
+seuil_inversion=-0.08
+#rendement sur la puissance fournie par les thermiques a la convection
+alp_bl_k=0.5
+#
+## temperature dependency of qsat for iterative procedure
+## temporary flag to switch between the old (CMIP5, D=0) and new versions (=1) versions 
+iflag_fisrtilp_qsat=4 
+## Reactivation du calcul d'un zmax continu pour les thermiques
+## reactivation d'une fermeture humide des thermiques  iflag_thermals_closure=2 (D=1)
+iflag_thermals_closure=1
+#
+# Flags wakes
+#************
+#Wake scheme switch (D:0, 0:AR4, 1:New Physics)
+iflag_wake=1
+#Constant value to be added to Alp (D: 0.)
+alp_offset=0.
+
+### nrlmd le 10/04/2012: Declenchement ###
+#Declenchement par couche-limite (0: deterministe, 1: stochastique)
+iflag_trig_bl=1
+#Surface critique au LCL
+s_trig=1.2e7
+#Temps de decorrelation entre scenes de cumulus
+tau_trig_shallow=1200
+#Temps de decorrelation entre scenes de cumulonimbus
+tau_trig_deep=1200
+# Pour la stabilisation
+random_notrig_max=0.99
+ok_adjwk=y
+tmax_fonte_cv=278.15
+
+
+### nrlmd le 10/04/2012: Fermeture ###
+#Fermeture par couche-limite (0: deterministe, 1: statistique, 2: statistique+convergence grande echelle)
+iflag_clos_bl=1
+### fin nrlmd le 10/04/2012
+
+## Modification de la plage de temperature pour la phase mixte liquide/glace des nuages
+## Controle par les t_glace_min/max, exposant_glace, iflag_t_glace (D=0)
+t_glace_min=233.13
+t_glace_max=273.13
+exposant_glace=1.
+iflag_t_glace=2
+iflag_ice_thermo=1
+
+#
+## The source terms of the TKE prognostic equation are diagnosed from tendencies (du, dv, dT) 
+## associated with subrgrid scaleamotions and treated as an additional heat source.
+## Controled by a key, iflag_ener_conserv (D=-1):
+## 0 : no conservation, -1 : old adhoc correction for kinetic E only (used for CMIP5) 
+## 1 : conservation, 101 : conversion from kinetic to heat only, 110 conversion from potential to heat only
+iflag_ener_conserv=1
+
+#Activation of corrections for water conservation (D: n). 
+##The default is n for compatibility with earlier versions.
+##generally it should be y.
+ok_conserv_q=y 
+
+iflag_cloudth_vert=2
+
+#
+lmixmin=10.
+#pbl_lmixmin_alpha=0.05
+
+#Bergeron
+iflag_bergeron=1
+
+
+### Convection
+#altitude, in hPa, above LCL at which buoybase is computed, def = original: -40
+dpbase=-40.
+#altitude, in hPa, above LCL below which Ep=0, def = original: 150, used only if flag_epKEorig=0
+pbcrit=150.0
+#altitude, in hPa, above LCL above which Ep=Epmax, def = original: 500, used only if flag_epKEorig=0
+ptcrit=500.0
+#reference fractional area of precipitating downdraughts, def = original: 0.01
+sigdz=0.01
+#unsaturated fraction of precipitating downdraughts, def = original: 0.15
+spfac=0.15
+#relaxation time of cloud base mass flux, def = original: 8000
+tau=8000.
+#flag for wb (= vert velocity at LFC); 0->wb=wbmax, 1->wb=f(plfc) bounded, 2->wb=f(plfc) linear, def = 1
+flag_wb=1
+#parameter used in the computation of wb, def = 6.
+wbmax=6.
+#
+ok_convstop=.False.
+#
+tau_stop=15000.
+#
+ok_intermittent=.False.
+#shedding coefficient (used when iflag_mix_adiab=1)
+coef_peel=0.25
+#
+### Clouds
+#Conversion cld water->precip: Ep=f(pressure) if 0; Ep=f(temperature) if 1, def=1
+flag_epKEorig=1
+#max condensed water concentration in kg/kg, def=0.0003, original: 0.0011, used only if flag_epKEorig=1
+elcrit=0.00025
+#temperature in degre Celsius, at which Ep=Epmax, def=-55., original: -55, used only if flag_epKEorig=1
+tlcrit=-40.
+
+### Wakes
+#coefficient in Wape->C_star formula: C_star=stark*sqrt(2*Wape), def: 0.33
+stark=0.33
+#lifting efficiency of wakes; Alp=Alpk*Fip, Fip=power provided by wakes, def: 0.25
+alpk=0.25
+#wake density = number of wake centers per m2, def: 8.E-12
+wdens_ref=8.E-12
+#multiplicative factor of the damping by gravity waves, def: 4.
+coefgw=4.
Index: /LMDZ5/branches/testing/libf/dynphy_lonlat/phydev/callphysiq_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/dynphy_lonlat/phydev/callphysiq_mod.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/dynphy_lonlat/phydev/callphysiq_mod.F90	(revision 2669)
@@ -12,5 +12,5 @@
                        jD_cur,jH_cur_split,zdt_split,                     &
                        zplev_omp,zplay_omp,                               &
-                       zphi_omp,zphis_omp,                                &
+                       zpk_omp,zphi_omp,zphis_omp,                        &
                        presnivs_omp,                                      &
                        zufi_omp,zvfi_omp,zrfi_omp,ztfi_omp,zqfi_omp,      &
@@ -34,4 +34,5 @@
   REAL,INTENT(IN) :: zplev_omp(klon,llm+1) ! interlayer pressure (Pa)
   REAL,INTENT(IN) :: zplay_omp(klon,llm) ! mid-layer pressure (Pa)
+  REAL,INTENT(IN) :: zpk_omp(klon,llm) ! Exner function
   REAL,INTENT(IN) :: zphi_omp(klon,llm) ! geopotential at midlayer
   REAL,INTENT(IN) :: zphis_omp(klon) ! surface geopotential
Index: /LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90
===================================================================
--- /LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90	(revision 2669)
@@ -111,4 +111,5 @@
   INTEGER :: flag_aerosol
   INTEGER :: flag_aerosol_strat
+  LOGICAL :: flag_bc_internal_mixture
   LOGICAL :: new_aod
   REAL    :: bl95_b0, bl95_b1
@@ -131,5 +132,5 @@
                    ok_ade, ok_aie, ok_cdnc, aerosol_couple,             &
                    flag_aerosol, flag_aerosol_strat, new_aod,           &
-                   bl95_b0, bl95_b1,                                    &
+                   flag_bc_internal_mixture, bl95_b0, bl95_b1,          &
                    read_climoz,                                         &
                    alp_offset)
Index: /LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90	(revision 2669)
@@ -9,5 +9,5 @@
                      nbp, communicator, &
                      punjours, pdayref,ptimestep, &
-                     rlatu,rlatv,rlonu,rlonv,aire,cu,cv,       &
+                     rlatudyn,rlatvdyn,rlonudyn,rlonvdyn,airedyn,cudyn,cvdyn, &
                      prad,pg,pr,pcpp,iflag_phys)
   USE dimphy, ONLY: init_dimphy
@@ -44,4 +44,6 @@
   USE mod_phys_lmdz_omp_data, ONLY: klon_omp 
 #endif
+  USE ioipsl_getin_p_mod, ONLY: getin_p
+  USE slab_heat_transp_mod, ONLY: ini_slab_transp_geom
   IMPLICIT NONE
 
@@ -52,6 +54,8 @@
 
   include "dimensions.h"
+  include "paramet.h"
   include "iniprint.h"
   include "tracstoke.h"
+  include "comgeom.h"
 
   REAL, INTENT (IN) :: prad ! radius of the planet (m)
@@ -65,11 +69,11 @@
   INTEGER, INTENT(IN) :: nbp ! number of physics columns for this MPI process
   INTEGER, INTENT(IN) :: communicator ! MPI communicator
-  REAL, INTENT (IN) :: rlatu(jj+1) ! latitudes of the physics grid
-  REAL, INTENT (IN) :: rlatv(jj) ! latitude boundaries of the physics grid
-  REAL, INTENT (IN) :: rlonv(ii+1) ! longitudes of the physics grid
-  REAL, INTENT (IN) :: rlonu(ii+1) ! longitude boundaries of the physics grid
-  REAL, INTENT (IN) :: aire(ii+1,jj+1) ! area of the dynamics grid (m2)
-  REAL, INTENT (IN) :: cu((ii+1)*(jj+1)) ! cu coeff. (u_covariant = cu * u)
-  REAL, INTENT (IN) :: cv((ii+1)*jj) ! cv coeff. (v_covariant = cv * v)
+  REAL, INTENT (IN) :: rlatudyn(jj+1) ! latitudes of the physics grid
+  REAL, INTENT (IN) :: rlatvdyn(jj) ! latitude boundaries of the physics grid
+  REAL, INTENT (IN) :: rlonvdyn(ii+1) ! longitudes of the physics grid
+  REAL, INTENT (IN) :: rlonudyn(ii+1) ! longitude boundaries of the physics grid
+  REAL, INTENT (IN) :: airedyn(ii+1,jj+1) ! area of the dynamics grid (m2)
+  REAL, INTENT (IN) :: cudyn((ii+1)*(jj+1)) ! cu coeff. (u_covariant = cu * u)
+  REAL, INTENT (IN) :: cvdyn((ii+1)*jj) ! cv coeff. (v_covariant = cv * v)
   INTEGER, INTENT (IN) :: pdayref ! reference day of for the simulation
   REAL, INTENT (IN) :: ptimestep !physics time step (s)
@@ -80,5 +84,8 @@
   CHARACTER (LEN=20) :: modname = 'iniphysiq'
   CHARACTER (LEN=80) :: abort_message
-
+  
+  LOGICAL :: slab_hdiff
+  INTEGER :: slab_ekman
+  CHARACTER (LEN = 6) :: type_ocean 
 
 #ifndef CPP_PARA
@@ -92,7 +99,7 @@
   CALL inigeomphy(ii,jj,nlayer, &
                nbp, communicator, &
-               rlatu,rlatv, &
-               rlonu,rlonv, &
-               aire,cu,cv)
+               rlatudyn,rlatvdyn, &
+               rlonudyn,rlonvdyn, &
+               airedyn,cudyn,cvdyn)
 
   ! --> now initialize things specific to the phylmd physics package
@@ -117,4 +124,18 @@
   ! Copy over "offline" settings
   CALL init_phystokenc(offline,istphy)
+
+  ! Initialization for slab heat transport
+  type_ocean="force"
+  CALL getin_p('type_ocean',type_ocean)
+  slab_hdiff=.FALSE.
+  CALL getin_p('slab_hdiff',slab_hdiff)
+  slab_ekman=0
+  CALL getin_p('slab_ekman',slab_ekman)
+  IF ((type_ocean=='slab').AND.(slab_hdiff.OR.(slab_ekman.GT.0))) THEN
+     CALL ini_slab_transp_geom(ip1jm,ip1jmp1,unsairez,fext,unsaire,&
+                                  cu,cuvsurcv,cv,cvusurcu, &
+                                  aire,apoln,apols, &
+                                  aireu,airev) 
+  END IF
 
   ! Initialize tracer names, numbers, etc. for physics
@@ -159,5 +180,5 @@
 #ifdef INCA
      CALL init_inca_dim(klon_omp,nbp_lev,nbp_lon,nbp_lat - 1, &
-          rlonu,rlatu,rlonv,rlatv)
+          rlonudyn,rlatudyn,rlonvdyn,rlatvdyn)
 #endif
   END IF
Index: DZ5/branches/testing/libf/phy_common/iophys.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phy_common/iophys.F90	(revision 2668)
+++ 	(revision )
@@ -1,132 +1,0 @@
-      subroutine iophys_ecrit(nom,lllm,titre,unite,px)
-      USE dimphy
-      USE mod_phys_lmdz_para
-      USE mod_grid_phy_lmdz
-      IMPLICIT NONE
-
-
-
-!  Ecriture de variables diagnostiques au choix dans la physique 
-!  dans un fichier NetCDF nomme  'diagfi'. Ces variables peuvent etre
-!  3d (ex : temperature), 2d (ex : temperature de surface), ou
-!  0d (pour un scalaire qui ne depend que du temps : ex : la longitude
-!  solaire)
-!  (ou encore 1d, dans le cas de testphys1d, pour sortir une colonne)
-!  La periode d'ecriture est donnee par 
-!  "ecritphy " regle dans le fichier de controle de run :  run.def
-!
-!    writediagfi peut etre appele de n'importe quelle subroutine
-!    de la physique, plusieurs fois. L'initialisation et la creation du
-!    fichier se fait au tout premier appel.
-!
-! WARNING : les variables dynamique (u,v,t,q,ps)
-!  sauvees par writediagfi avec une
-! date donnee sont legerement differentes que dans le fichier histoire car 
-! on ne leur a pas encore ajoute de la dissipation et de la physique !!!
-! IL est  RECOMMANDE d'ajouter les tendance physique a ces variables
-! avant l'ecriture dans diagfi (cf. physiq.F)
-!  
-! Modifs: Aug.2010 Ehouarn: enforce outputs to be real*4
-!
-!  parametres (input) :
-!  ----------
-!      unit : unite logique du fichier de sortie (toujours la meme)
-!      nom  : nom de la variable a sortir (chaine de caracteres)
-!      titre: titre de la variable (chaine de caracteres)
-!      unite : unite de la variable (chaine de caracteres)
-!      px : variable a sortir (real 0, 1, 2, ou 3d)
-!
-!=================================================================
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "netcdf.inc"
-#include "iotd.h"
-
-
-! Arguments on input:
-      integer lllm
-      character (len=*) :: nom,titre,unite
-      integer imjmax
-      parameter (imjmax=100000)
-      real px(klon_omp,lllm)
-      real xglo(klon_glo,lllm)
-      real zx(iim,jjp1,lllm)
-
-
-      CALL Gather(px,xglo)
-!$OMP MASTER
-      IF (is_mpi_root) THEN       
-        CALL Grid1Dto2D_glo(xglo,zx)
-        call iotd_ecrit(nom,lllm,titre,unite,zx)
-      ENDIF
-!$OMP END MASTER
-
-      return
-      end
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-      SUBROUTINE iophys_ini
-      USE mod_phys_lmdz_para
-      USE vertical_layers_mod, ONLY: presnivs
-      USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
-      USE dimphy, ONLY: klev
-
-      IMPLICIT NONE
-
-!=======================================================================
-!
-!   Auteur:  L. Fairhead  ,  P. Le Van, Y. Wanherdrick, F. Forget
-!   -------
-!
-!   Objet:
-!   ------
-!
-!   'Initialize' the diagfi.nc file: write down dimensions as well
-!   as time-independent fields (e.g: geopotential, mesh area, ...)
-!
-!=======================================================================
-!-----------------------------------------------------------------------
-!   Declarations:
-!   -------------
-
-real pi
-
-!   Arguments:
-!   ----------
-
-!$OMP MASTER
-    IF (is_mpi_root) THEN       
-pi=2.*asin(1.)
-call iotd_ini('phys.nc   ', &
-!iim,jjp1,llm,rlonv(1:iim)*180./pi,rlatu*180./pi,presnivs)
-size(lon_reg),size(lat_reg),klev,lon_reg(:)*180./pi,lat_reg*180./pi,presnivs)
-    ENDIF
-!$OMP END MASTER
-
-      END
-
-#ifdef und
-      SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
-      IMPLICIT none
-
-!=======================================================================
-      INTEGER nfield,nlon,iim,jjmp1, jjm
-      REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
-
-      INTEGER i, n, ig
-
-      jjm = jjmp1 - 1
-      DO n = 1, nfield
-         DO i=1,iim
-            ecrit(i,n) = fi(1,n)
-            ecrit(i+jjm*iim,n) = fi(nlon,n)
-         ENDDO
-         DO ig = 1, nlon - 2
-           ecrit(iim+ig,n) = fi(1+ig,n)
-         ENDDO
-      ENDDO
-      RETURN
-      END
-
-#endif
Index: /LMDZ5/branches/testing/libf/phydev/physiq_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phydev/physiq_mod.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phydev/physiq_mod.F90	(revision 2669)
@@ -25,6 +25,6 @@
 #ifdef CPP_XIOS
       USE xios, ONLY: xios_update_calendar
-      USE wxios, only: wxios_add_vaxis, wxios_set_timestep, wxios_closedef, &
-                       histwrite_phy
+      USE wxios, only: wxios_add_vaxis, wxios_set_cal, wxios_closedef
+      USE iophy, ONLY: histwrite_phy
 #endif
 
@@ -129,14 +129,14 @@
 !XIOS
     ! Declare available vertical axes to be used in output files:    
-    !CALL wxios_add_vaxis("presnivs", "dummy-not-used", klev, presnivs)
     CALL wxios_add_vaxis("presnivs", klev, presnivs)
 
-    ! Declare time step length (in s):
-    CALL wxios_set_timestep(dtime)
-
+    ! Declare calendar and time step
+    CALL wxios_set_cal(dtime,"earth_360d",1,1,1,0.0,1,1,1,0.0)
+    
     !Finalize the context:
     CALL wxios_closedef()
 #endif
 !$OMP END MASTER
+!$OMP BARRIER
 endif ! of if (debut)
 
Index: /LMDZ5/branches/testing/libf/phylmd/Dust/phytracr_spl_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/Dust/phytracr_spl_mod.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/Dust/phytracr_spl_mod.F90	(revision 2669)
@@ -1170,18 +1170,18 @@
        !print *,nbtr
        do it=1,nbtr
-        print *, it, tname(it+2)
-        if (tname(it+2) == 'PREC' ) then
+        print *, it, tname(it+nqo)
+        if (tname(it+nqo) == 'PREC' ) then
             id_prec=it
         endif
-        if (tname(it+2) == 'FINE' ) then
+        if (tname(it+nqo) == 'FINE' ) then
             id_fine=it
         endif
-        if (tname(it+2) == 'COSS' ) then
+        if (tname(it+nqo) == 'COSS' ) then
             id_coss=it
         endif
-        if (tname(it+2) == 'CODU' ) then
+        if (tname(it+nqo) == 'CODU' ) then
             id_codu=it
         endif
-        if (tname(it+2) == 'SCDU' ) then
+        if (tname(it+nqo) == 'SCDU' ) then
             id_scdu=it
         endif
@@ -2621,10 +2621,9 @@
          call iophys_ecrit('fav'//str2,1,'SOURCE','',source_tr(:,it))
       enddo
-#endif
-
       do it=1,nbtr
          write(str2,'(i2.2)') it
          call iophys_ecrit('TRB'//str2,klev,'SOURCE','',tr_seri(:,:,it))
       enddo
+#endif
 
 
@@ -2650,5 +2649,4 @@
 
 #ifdef IOPHYS_DUST
-
       do it=1,nbtr
          write(str2,'(i2.2)') it
@@ -2676,5 +2674,4 @@
 !
 #ifdef IOPHYS_DUST
-!
       print *,'INPUT TO PRECUREMISSION'
          call iophys_ecrit('ftsol',4,'ftsol','',ftsol)
Index: /LMDZ5/branches/testing/libf/phylmd/change_srf_frac_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/change_srf_frac_mod.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/change_srf_frac_mod.F90	(revision 2669)
@@ -87,6 +87,11 @@
        CALL limit_read_frac(itime, dtime, jour, pctsrf, is_modified)
     CASE ('slab')
+       IF (version_ocean == 'sicOBS'.OR. version_ocean == 'sicNO') THEN
+       ! Read fraction from limit.nc
+           CALL limit_read_frac(itime, dtime, jour, pctsrf, is_modified)
+       ELSE
        ! Get fraction from slab module
-       CALL ocean_slab_frac(itime, dtime, jour, pctsrf, is_modified)
+           CALL ocean_slab_frac(itime, dtime, jour, pctsrf, is_modified)
+       ENDIF
     CASE ('couple')
        ! Get fraction from the coupler
Index: /LMDZ5/branches/testing/libf/phylmd/cloudth.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/cloudth.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/cloudth.F90	(revision 2669)
@@ -9,5 +9,5 @@
 
 !===========================================================================
-! Auteur : Arnaud Octavio Jam (LMD/CNRS)
+! Author : Arnaud Octavio Jam (LMD/CNRS)
 ! Date : 25 Mai 2010
 ! Objet : calcule les valeurs de qc et rneb dans les thermiques
@@ -49,12 +49,11 @@
       REAL rneb(ngrid,klev)
       REAL t(ngrid,klev)
-      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
+      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
       REAL rdd,cppd,Lv
       REAL alth,alenv,ath,aenv
-      REAL sth,senv,sigma1s,sigma2s,xth,xenv
+      REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
       REAL Tbef,zdelta,qsatbef,zcor
       REAL qlbef  
-      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
-      
+      REAL ratqs(ngrid,klev) ! Determine the width of the vapour distribution
       REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
       REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) 
@@ -63,9 +62,4 @@
 
 
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! Gestion de deux versions de cloudth
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
       IF (iflag_cloudth_vert.GE.1) THEN
@@ -82,236 +76,4 @@
 ! Initialisation des variables r?elles
 !-------------------------------------------------------------------------------
-      sigma1(:,:)=0.
-      sigma2(:,:)=0.
-      qlth(:,:)=0.
-      qlenv(:,:)=0.  
-      qltot(:,:)=0.
-      rneb(:,:)=0.
-      qcloud(:)=0.
-      cth(:,:)=0.
-      cenv(:,:)=0.
-      ctot(:,:)=0.
-      qsatmmussig1=0.
-      qsatmmussig2=0.
-      rdd=287.04
-      cppd=1005.7
-      pi=3.14159 
-      Lv=2.5e6
-      sqrt2pi=sqrt(2.*pi)
-
-
-
-!-------------------------------------------------------------------------------
-! Calcul de la fraction du thermique et des ?cart-types des distributions
-!-------------------------------------------------------------------------------                 
-      do ind1=1,ngrid
-
-      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then 
-
-      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
-
-
-!      zqenv(ind1)=po(ind1)
-      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
-      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
-      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
-      qsatbef=MIN(0.5,qsatbef)
-      zcor=1./(1.-retv*qsatbef)
-      qsatbef=qsatbef*zcor
-      zqsatenv(ind1,ind2)=qsatbef
-
-
-
-
-      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)  
-      aenv=1./(1.+(alenv*Lv/cppd))
-      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) 
-
-
-
-
-      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
-      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
-      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
-      qsatbef=MIN(0.5,qsatbef)
-      zcor=1./(1.-retv*qsatbef)
-      qsatbef=qsatbef*zcor
-      zqsatth(ind1,ind2)=qsatbef
-            
-      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
-      ath=1./(1.+(alth*Lv/cppd))
-      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) 
-      
-      
-
-!------------------------------------------------------------------------------
-! Calcul des ?cart-types pour s
-!------------------------------------------------------------------------------
-
-!      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
-!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2) 
-!       if (paprs(ind1,ind2).gt.90000) then
-!       ratqs(ind1,ind2)=0.002
-!       else 
-!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
-!       endif
-       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
-       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) 
-!       sigma1s=ratqs(ind1,ind2)*po(ind1)
-!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003  
- 
-!------------------------------------------------------------------------------
-! Calcul de l'eau condens?e et de la couverture nuageuse
-!------------------------------------------------------------------------------
-      sqrt2pi=sqrt(2.*pi)
-      xth=sth/(sqrt(2.)*sigma2s)
-      xenv=senv/(sqrt(2.)*sigma1s)
-      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
-      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) 
-      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
-
-      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
-      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
-      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-      if (ctot(ind1,ind2).lt.1.e-10) then
-      ctot(ind1,ind2)=0.
-      qcloud(ind1)=zqsatenv(ind1,ind2) 
-
-      else   
-                
-      ctot(ind1,ind2)=ctot(ind1,ind2) 
-      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
-
-      endif                           
-      
-          
-!     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
-
-
-      else  ! gaussienne environnement seule
-      
-      zqenv(ind1)=po(ind1)
-      Tbef=t(ind1,ind2)
-      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
-      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
-      qsatbef=MIN(0.5,qsatbef)
-      zcor=1./(1.-retv*qsatbef)
-      qsatbef=qsatbef*zcor
-      zqsatenv(ind1,ind2)=qsatbef
-      
-
-!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
-      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
-      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)  
-      aenv=1./(1.+(alenv*Lv/cppd))
-      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) 
-      
-
-      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
-
-      sqrt2pi=sqrt(2.*pi)
-      xenv=senv/(sqrt(2.)*sigma1s)
-      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
-      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
-      
-      if (ctot(ind1,ind2).lt.1.e-3) then
-      ctot(ind1,ind2)=0.
-      qcloud(ind1)=zqsatenv(ind1,ind2) 
-
-      else   
-                
-      ctot(ind1,ind2)=ctot(ind1,ind2) 
-      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
-
-      endif    
- 
- 
- 
- 
- 
- 
-      endif   
-      enddo
-     
-      return
-      end
-
-
-
-!===========================================================================
-     SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
-     &           ztv,po,zqta,fraca, & 
-     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
-     &           ratqs,zqs,t)
-
-!===========================================================================
-! Auteur : Arnaud Octavio Jam (LMD/CNRS)
-! Date : 25 Mai 2010
-! Objet : calcule les valeurs de qc et rneb dans les thermiques
-!===========================================================================
-
-
-      USE ioipsl_getin_p_mod, ONLY : getin_p
-
-      IMPLICIT NONE
-
-#include "YOMCST.h"
-#include "YOETHF.h"
-#include "FCTTRE.h"
-#include "thermcell.h"
-#include "nuage.h"
-      
-      INTEGER itap,ind1,ind2
-      INTEGER ngrid,klev,klon,l,ig 
-      
-      REAL ztv(ngrid,klev)
-      REAL po(ngrid)
-      REAL zqenv(ngrid)   
-      REAL zqta(ngrid,klev)
-          
-      REAL fraca(ngrid,klev+1)
-      REAL zpspsk(ngrid,klev)
-      REAL paprs(ngrid,klev+1)
-      REAL ztla(ngrid,klev)
-      REAL zthl(ngrid,klev)
-
-      REAL zqsatth(ngrid,klev)
-      REAL zqsatenv(ngrid,klev)
-      
-      
-      REAL sigma1(ngrid,klev)                                                         
-      REAL sigma2(ngrid,klev)
-      REAL qlth(ngrid,klev)
-      REAL qlenv(ngrid,klev)
-      REAL qltot(ngrid,klev) 
-      REAL cth(ngrid,klev)  
-      REAL cenv(ngrid,klev)   
-      REAL ctot(ngrid,klev)
-      REAL rneb(ngrid,klev)
-      REAL t(ngrid,klev)                                                                  
-      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
-      REAL rdd,cppd,Lv,sqrt2,sqrtpi
-      REAL alth,alenv,ath,aenv
-      REAL sth,senv,sigma1s,sigma2s,xth,xenv
-      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
-      REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
-      REAL Tbef,zdelta,qsatbef,zcor
-      REAL qlbef  
-      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
-      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
-      ! (J Jouhaud, JL Dufresne, JB Madeleine)
-      REAL,SAVE :: vert_alpha
-      LOGICAL, SAVE :: firstcall = .TRUE.
-      
-      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
-      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) 
-      REAL zqs(ngrid), qcloud(ngrid)
-      REAL erf
-
-!------------------------------------------------------------------------------
-! Initialisation des variables r?elles
-!------------------------------------------------------------------------------
       sigma1(:,:)=0.
       sigma2(:,:)=0.
@@ -334,13 +96,7 @@
       sqrtpi=sqrt(pi)
 
-      IF (firstcall) THEN
-        vert_alpha=0.5
-        CALL getin_p('cloudth_vert_alpha',vert_alpha)
-        WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
-        firstcall=.FALSE.
-      ENDIF
-
-!-------------------------------------------------------------------------------
-! Calcul de la fraction du thermique et des ?cart-types des distributions
+
+!-------------------------------------------------------------------------------
+! Cloud fraction in the thermals and standard deviation of the PDFs
 !-------------------------------------------------------------------------------                 
       do ind1=1,ngrid
@@ -350,9 +106,7 @@
       zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
 
-
-!      zqenv(ind1)=po(ind1)
       Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
       zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
-      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
+      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
       qsatbef=MIN(0.5,qsatbef)
       zcor=1./(1.-retv*qsatbef)
@@ -361,11 +115,12 @@
 
 
-
-
-      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)  
-      aenv=1./(1.+(alenv*Lv/cppd))
-      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) 
-
-
+      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
+      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
+      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
+
+!po = qt de l'environnement ET des thermique
+!zqenv = qt environnement
+!zqsatenv = qsat environnement
+!zthl = Tl environnement
 
 
@@ -378,16 +133,238 @@
       zqsatth(ind1,ind2)=qsatbef
             
-      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
-      ath=1./(1.+(alth*Lv/cppd))
-      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) 
-      
-      
-
-!------------------------------------------------------------------------------
-! Calcul des ?cart-types pour s
-!------------------------------------------------------------------------------
-
-      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
-      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2) 
+      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
+      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
+      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
+      
+!zqta = qt thermals
+!zqsatth = qsat thermals
+!ztla = Tl thermals
+
+!------------------------------------------------------------------------------
+! s standard deviations
+!------------------------------------------------------------------------------
+
+!     tests
+!     sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
+!     sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+ratqs(ind1,ind2)*po(ind1)
+!     sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
+!     final option
+      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
+      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
+ 
+!------------------------------------------------------------------------------
+! Condensed water and cloud cover
+!------------------------------------------------------------------------------
+      xth=sth/(sqrt2*sigma2s)
+      xenv=senv/(sqrt2*sigma1s)
+      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))       !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
+      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))     !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
+      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
+
+      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2))
+      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
+      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
+
+      if (ctot(ind1,ind2).lt.1.e-10) then
+      ctot(ind1,ind2)=0.
+      qcloud(ind1)=zqsatenv(ind1,ind2) 
+      else
+      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
+      endif
+
+      else  ! Environnement only, follow the if l.110
+      
+      zqenv(ind1)=po(ind1)
+      Tbef=t(ind1,ind2)
+      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
+      qsatbef=MIN(0.5,qsatbef)
+      zcor=1./(1.-retv*qsatbef)
+      qsatbef=qsatbef*zcor
+      zqsatenv(ind1,ind2)=qsatbef
+
+!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
+      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
+      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)  
+      aenv=1./(1.+(alenv*Lv/cppd))
+      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
+      
+      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
+
+      xenv=senv/(sqrt2*sigma1s)
+      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
+      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
+
+      if (ctot(ind1,ind2).lt.1.e-3) then
+      ctot(ind1,ind2)=0.
+      qcloud(ind1)=zqsatenv(ind1,ind2)
+      else   
+      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
+      endif
+
+
+      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
+      enddo       ! from the loop on ngrid l.108
+      return
+      end
+
+
+
+!===========================================================================
+     SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
+     &           ztv,po,zqta,fraca, & 
+     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
+     &           ratqs,zqs,t)
+
+!===========================================================================
+! Auteur : Arnaud Octavio Jam (LMD/CNRS)
+! Date : 25 Mai 2010
+! Objet : calcule les valeurs de qc et rneb dans les thermiques
+!===========================================================================
+
+
+      USE ioipsl_getin_p_mod, ONLY : getin_p
+
+      IMPLICIT NONE
+
+#include "YOMCST.h"
+#include "YOETHF.h"
+#include "FCTTRE.h"
+#include "thermcell.h"
+#include "nuage.h"
+      
+      INTEGER itap,ind1,ind2
+      INTEGER ngrid,klev,klon,l,ig 
+      
+      REAL ztv(ngrid,klev)
+      REAL po(ngrid)
+      REAL zqenv(ngrid)   
+      REAL zqta(ngrid,klev)
+          
+      REAL fraca(ngrid,klev+1)
+      REAL zpspsk(ngrid,klev)
+      REAL paprs(ngrid,klev+1)
+      REAL ztla(ngrid,klev)
+      REAL zthl(ngrid,klev)
+
+      REAL zqsatth(ngrid,klev)
+      REAL zqsatenv(ngrid,klev)
+      
+      
+      REAL sigma1(ngrid,klev)                                                         
+      REAL sigma2(ngrid,klev)
+      REAL qlth(ngrid,klev)
+      REAL qlenv(ngrid,klev)
+      REAL qltot(ngrid,klev) 
+      REAL cth(ngrid,klev)  
+      REAL cenv(ngrid,klev)   
+      REAL ctot(ngrid,klev)
+      REAL rneb(ngrid,klev)
+      REAL t(ngrid,klev)                                                                  
+      REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
+      REAL rdd,cppd,Lv
+      REAL alth,alenv,ath,aenv
+      REAL sth,senv,sigma1s,sigma2s,xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
+      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
+      REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
+      REAL Tbef,zdelta,qsatbef,zcor
+      REAL qlbef  
+      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
+      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
+      ! (J Jouhaud, JL Dufresne, JB Madeleine)
+      REAL,SAVE :: vert_alpha
+      LOGICAL, SAVE :: firstcall = .TRUE.
+
+      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
+      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) 
+      REAL zqs(ngrid), qcloud(ngrid)
+      REAL erf
+
+!------------------------------------------------------------------------------
+! Initialize
+!------------------------------------------------------------------------------
+      sigma1(:,:)=0.
+      sigma2(:,:)=0.
+      qlth(:,:)=0.
+      qlenv(:,:)=0.  
+      qltot(:,:)=0.
+      rneb(:,:)=0.
+      qcloud(:)=0.
+      cth(:,:)=0.
+      cenv(:,:)=0.
+      ctot(:,:)=0.
+      qsatmmussig1=0.
+      qsatmmussig2=0.
+      rdd=287.04
+      cppd=1005.7
+      pi=3.14159 
+      Lv=2.5e6
+      sqrt2pi=sqrt(2.*pi)
+      sqrt2=sqrt(2.)
+      sqrtpi=sqrt(pi)
+
+      IF (firstcall) THEN
+        vert_alpha=0.5
+        CALL getin_p('cloudth_vert_alpha',vert_alpha)
+        WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
+        firstcall=.FALSE.
+      ENDIF
+
+!-------------------------------------------------------------------------------
+! Calcul de la fraction du thermique et des ecart-types des distributions
+!-------------------------------------------------------------------------------                 
+      do ind1=1,ngrid
+
+      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement
+
+      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv
+
+
+      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
+      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
+      qsatbef=MIN(0.5,qsatbef)
+      zcor=1./(1.-retv*qsatbef)
+      qsatbef=qsatbef*zcor
+      zqsatenv(ind1,ind2)=qsatbef
+
+
+      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
+      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
+      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
+
+!zqenv = qt environnement
+!zqsatenv = qsat environnement
+!zthl = Tl environnement
+
+
+      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
+      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
+      qsatbef=MIN(0.5,qsatbef)
+      zcor=1./(1.-retv*qsatbef)
+      qsatbef=qsatbef*zcor
+      zqsatth(ind1,ind2)=qsatbef
+            
+      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
+      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
+      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
+      
+      
+!zqta = qt thermals
+!zqsatth = qsat thermals
+!ztla = Tl thermals
+
+!------------------------------------------------------------------------------
+! s standard deviation
+!------------------------------------------------------------------------------
+
+      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
+      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
+!      tests
+!      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
+!      sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1)
+!      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2) 
+!      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2)
 !       if (paprs(ind1,ind2).gt.90000) then
 !       ratqs(ind1,ind2)=0.002
@@ -400,28 +377,12 @@
 !      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003  
  
-!------------------------------------------------------------------------------
-! Calcul de l'eau condens?e et de la couverture nuageuse
-!------------------------------------------------------------------------------
-      sqrt2pi=sqrt(2.*pi)
-      xth=sth/(sqrt(2.)*sigma2s)
-      xenv=senv/(sqrt(2.)*sigma1s)
-      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
-      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) 
-      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
-
-      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
-      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
-      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
-     
        IF (iflag_cloudth_vert == 1) THEN
 !-------------------------------------------------------------------------------
-!  Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
-!-------------------------------------------------------------------------------
-!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
-!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
+!  Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
+!-------------------------------------------------------------------------------
+
       deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
       deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
-!      deltasenv=aenv*0.01*po(ind1)
-!     deltasth=ath*0.01*zqta(ind1,ind2)    
+
       xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
       xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
@@ -435,4 +396,5 @@
       ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
 
+      ! Environment
       IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
       IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
@@ -441,17 +403,11 @@
 
       qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
-!      qlenv(ind1,ind2)=IntJ 
-!      print*, qlenv(ind1,ind2),'VERIF EAU'
-
-
+
+      ! Thermal
       IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
-!      IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2))
-!      IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1))
       IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
       IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
       IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
       qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
-!      qlth(ind1,ind2)=IntJ
-!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
       qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
 
@@ -459,5 +415,5 @@
 
 !-------------------------------------------------------------------------------
-!  Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
+!  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
 !-------------------------------------------------------------------------------
 !      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
@@ -470,8 +426,10 @@
       xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
       xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
+      exp_xenv1 = exp(-1.*xenv1**2)
+      exp_xenv2 = exp(-1.*xenv2**2)
       xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
       xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
-!     coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
-!     coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
+      exp_xth1 = exp(-1.*xth1**2)
+      exp_xth2 = exp(-1.*xth2**2)
       
       cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
@@ -479,33 +437,32 @@
       ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
 
-      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp(-1.*xenv2**2)
+      
+      !environnement
+      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
+      if (deltasenv .lt. 1.e-10) then 
+      qlenv(ind1,ind2)=IntJ
+      else
       IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
-      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
-      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp(-1.*xenv1**2)-exp(-1.*xenv2**2))
-
-!      IntI1=0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
-!      IntI2=xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
-!      IntI3=0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
-
+      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
+      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
       qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
-!      qlenv(ind1,ind2)=IntJ 
-!      print*, qlenv(ind1,ind2),'VERIF EAU'
-
-      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp(-1.*xth2**2)
+      endif
+
+
+      !thermique
+      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
+      if (deltasth .lt. 1.e-10) then ! Seuil a choisir !!!
+      qlth(ind1,ind2)=IntJ
+      else
       IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
-      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
-      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp(-1.*xth1**2)-exp(-1.*xth2**2))
-      
+      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
+      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
       qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
-!      qlth(ind1,ind2)=IntJ
-!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
+      endif
+
       qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
-      
-
 
 
       ENDIF ! of if (iflag_cloudth_vert==1 or 2)
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
       if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
@@ -515,5 +472,4 @@
       else 
                 
-      ctot(ind1,ind2)=ctot(ind1,ind2) 
       qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
 !      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
@@ -521,11 +477,6 @@
 
       endif  
-                       
-      
-          
-!     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
-
-
-      else  ! gaussienne environnement seule
+
+      else  ! Environment only
       
       zqenv(ind1)=po(ind1)
@@ -539,5 +490,5 @@
       
 
-!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
+!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
       zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
       alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)  
@@ -548,8 +499,7 @@
       sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
 
-      sqrt2pi=sqrt(2.*pi)
-      xenv=senv/(sqrt(2.)*sigma1s)
+      xenv=senv/(sqrt2*sigma1s)
       ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
-      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
+      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
       
       if (ctot(ind1,ind2).lt.1.e-3) then
@@ -559,16 +509,10 @@
       else   
                 
-      ctot(ind1,ind2)=ctot(ind1,ind2) 
       qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
 
       endif    
  
- 
- 
- 
- 
- 
-      endif   
-      enddo
+      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
+      enddo       ! from the loop on ngrid l.333
      
       return
Index: /LMDZ5/branches/testing/libf/phylmd/conf_phys_m.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/conf_phys_m.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/conf_phys_m.F90	(revision 2669)
@@ -19,5 +19,5 @@
        ok_ade, ok_aie, ok_cdnc, aerosol_couple, &
        flag_aerosol, flag_aerosol_strat, new_aod, &
-       bl95_b0, bl95_b1,&
+       flag_bc_internal_mixture, bl95_b0, bl95_b1,&
        read_climoz, &
        alp_offset)
@@ -64,4 +64,5 @@
     ! ok_cdnc, ok cloud droplet number concentration
     ! flag_aerosol_strat : flag pour les aerosols stratos
+    ! flag_bc_internal_mixture : use BC internal mixture if true
     ! bl95_b*: parameters in the formula to link CDNC to aerosol mass conc 
     !
@@ -77,4 +78,5 @@
     INTEGER              :: flag_aerosol
     INTEGER              :: flag_aerosol_strat
+    LOGICAL              :: flag_bc_internal_mixture
     LOGICAL              :: new_aod
     REAL                 :: bl95_b0, bl95_b1
@@ -95,4 +97,5 @@
     INTEGER, SAVE       :: flag_aerosol_omp
     INTEGER, SAVE       :: flag_aerosol_strat_omp
+    LOGICAL, SAVE       :: flag_bc_internal_mixture_omp
     LOGICAL, SAVE       :: new_aod_omp
     REAL,SAVE           :: bl95_b0_omp, bl95_b1_omp
@@ -398,4 +401,15 @@
     flag_aerosol_strat_omp = 0
     CALL getin('flag_aerosol_strat',flag_aerosol_strat_omp)
+
+    !
+    !Config Key  = flag_bc_internal_mixture
+    !Config Desc = state of mixture for BC aerosols
+    ! - n = external mixture
+    ! - y = internal mixture
+    !Config Def  = n
+    !Config Help = Used in physiq.F / aeropt
+    !
+    flag_bc_internal_mixture_omp = .false.
+    CALL getin('flag_bc_internal_mixture',flag_bc_internal_mixture_omp)
 
     ! Temporary variable for testing purpose!
@@ -2136,4 +2150,5 @@
     flag_aerosol=flag_aerosol_omp
     flag_aerosol_strat=flag_aerosol_strat_omp
+    flag_bc_internal_mixture=flag_bc_internal_mixture_omp
     new_aod=new_aod_omp
     aer_type = aer_type_omp
@@ -2311,4 +2326,15 @@
     IF (ok_aie .AND. .NOT. ok_cdnc) THEN
        CALL abort_physic('conf_phys', 'ok_cdnc must be set to y if ok_aie is activated',1)
+    ENDIF
+
+    ! BC internal mixture is only possible with RRTM & NSW=6 & flag_aerosol=6 or aerosol_couple
+    IF (flag_bc_internal_mixture .AND. NSW.NE.6) THEN 
+       CALL abort_physic('conf_phys', 'flag_bc_internal_mixture can only be activated with NSW=6',1)
+    ENDIF
+    IF (flag_bc_internal_mixture .AND. iflag_rrtm.NE.1) THEN 
+       CALL abort_physic('conf_phys', 'flag_bc_internal_mixture can only be activated with RRTM',1)
+    ENDIF
+    IF (flag_bc_internal_mixture .AND. flag_aerosol.NE.6) THEN 
+       CALL abort_physic('conf_phys', 'flag_bc_internal_mixture can only be activated with flag_aerosol=6',1)
     ENDIF
 
Index: /LMDZ5/branches/testing/libf/phylmd/cosp/cosp_output_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/cosp/cosp_output_mod.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/cosp/cosp_output_mod.F90	(revision 2669)
@@ -42,5 +42,5 @@
          "clmcalipso", "Lidar Mid-level Cloud Fraction", "%", (/ ('', i=1, 3) /))
   TYPE(ctrl_outcosp), SAVE :: o_clhcalipso = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
-         "clhcalipso", "Lidar Hight-level Cloud Fraction", "%", (/ ('', i=1, 3) /))
+         "clhcalipso", "Lidar High-level Cloud Fraction", "%", (/ ('', i=1, 3) /))
   TYPE(ctrl_outcosp), SAVE :: o_cltcalipso = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
          "cltcalipso", "Lidar Total Cloud Fraction", "%", (/ ('', i=1, 3) /)) 
@@ -145,5 +145,5 @@
          "clmmodis", "MODIS Mid-level Cloud Fraction", "1", (/ ('', i=1, 3) /))
   TYPE(ctrl_outcosp), SAVE :: o_clhmodis = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
-         "clhmodis", "MODIS Hight-level Cloud Fraction", "1", (/ ('', i=1, 3) /))
+         "clhmodis", "MODIS High-level Cloud Fraction", "1", (/ ('', i=1, 3) /))
   TYPE(ctrl_outcosp), SAVE :: o_cltmodis = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
          "cltmodis", "MODIS Total Cloud Fraction", "1", (/ ('', i=1, 3) /)) 
Index: /LMDZ5/branches/testing/libf/phylmd/cva_driver.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/cva_driver.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/cva_driver.F90	(revision 2669)
@@ -1077,11 +1077,11 @@
         DO k = 1,nd
         write (6, '(i4,5(1x,e13.6))'), &
-          k, mp(idcum(igout),k), water(idcum(igout),k), ice(idcum(igout),k), &
-           evap(idcum(igout),k), fondue(idcum(igout),k)
+          k, mp(igout,k), water(igout,k), ice(igout,k), &
+           evap(igout,k), fondue(igout,k)
         ENDDO
         Print *, 'cva_driver after cv3_unsat: wdtrainA, wdtrainM '
         DO k = 1,nd
         write (6, '(i4,2(1x,e13.6))'), &
-           k, wdtrainA(idcum(igout),k), wdtrainM(idcum(igout),k)
+           k, wdtrainA(igout,k), wdtrainM(igout,k)
         ENDDO
       ENDIF
@@ -1129,12 +1129,12 @@
 !
       IF (debut) THEN
-        PRINT *, ' cv3_yield -> fqd(1) = ', fqd(idcum(igout), 1)
+        PRINT *, ' cv3_yield -> fqd(1) = ', fqd(igout, 1)
       END IF !(debut) THEN
 !   
       IF (prt_level >= 10) THEN
         Print *, 'cva_driver after cv3_yield:ft(1) , ftd(1) ', &
-                    ft(idcum(igout),1), ftd(idcum(igout),1)
+                    ft(igout,1), ftd(igout,1)
         Print *, 'cva_driver after cv3_yield:fq(1) , fqd(1) ', &
-                    fq(idcum(igout),1), fqd(idcum(igout),1)
+                    fq(igout,1), fqd(igout,1)
       ENDIF
 !   
Index: /LMDZ5/branches/testing/libf/phylmd/dimphy.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/dimphy.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/dimphy.F90	(revision 2669)
@@ -9,7 +9,6 @@
   INTEGER,SAVE :: klevm1
   INTEGER,SAVE :: kflev
-  INTEGER,SAVE :: nslay
 
-!$OMP THREADPRIVATE(klon,kfdia,kidia,kdlon,nslay)
+!$OMP THREADPRIVATE(klon,kfdia,kidia,kdlon)
   REAL,save,allocatable,dimension(:) :: zmasq
 !$OMP THREADPRIVATE(zmasq)   
@@ -24,5 +23,4 @@
     
     klon=klon0
-    nslay=1 ! Slab, provisoire (F. Codron)
     kdlon=klon
     kidia=1
Index: /LMDZ5/branches/testing/libf/phylmd/fisrtilp.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/fisrtilp.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/fisrtilp.F90	(revision 2669)
@@ -624,5 +624,6 @@
                     zdelta = MAX(0.,SIGN(1.,t_glace_min_old-Tbef(i)))
                     else if (iflag_t_glace.ge.1) then
-                    zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
+                    zdelta = MAX(0.,SIGN(1.,t_glace_max-Tbef(i)))
+! BUG corrige par JYG   zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
                     endif
                  endif
Index: /LMDZ5/branches/testing/libf/phylmd/iophy.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/iophy.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/iophy.F90	(revision 2669)
@@ -17,5 +17,10 @@
 #ifdef CPP_XIOS
   INTERFACE histwrite_phy
-    MODULE PROCEDURE histwrite2d_phy,histwrite3d_phy,histwrite2d_phy_old,histwrite3d_phy_old,histwrite2d_xios,histwrite3d_xios
+!#ifdef CPP_XIOSnew
+    MODULE PROCEDURE histwrite2d_phy,histwrite3d_phy,histwrite2d_phy_old,histwrite3d_phy_old,histwrite2d_xios,histwrite3d_xios,histwrite0d_xios
+!#else
+!    MODULE PROCEDURE histwrite2d_phy,histwrite3d_phy,histwrite2d_phy_old,histwrite3d_phy_old,histwrite2d_xios,histwrite3d_xios
+!#endif
+
   END INTERFACE
 #else
@@ -1331,4 +1336,20 @@
   IF (prt_level >= 10) write(lunout,*)'End histrwrite3d_xios ',field_name
   END SUBROUTINE histwrite3d_xios
+
+#ifdef CPP_XIOS
+  SUBROUTINE histwrite0d_xios(field_name, field)
+  USE xios, only: xios_send_field
+  IMPLICIT NONE
+
+    CHARACTER(LEN=*), INTENT(IN) :: field_name
+    REAL, INTENT(IN) :: field ! --> scalar
+
+!$OMP MASTER
+   CALL xios_send_field(field_name, field)
+!$OMP END MASTER
+
+  END SUBROUTINE histwrite0d_xios
+#endif
+
 #endif
 end module iophy
Index: /LMDZ5/branches/testing/libf/phylmd/iophys.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/iophys.F90	(revision 2669)
+++ /LMDZ5/branches/testing/libf/phylmd/iophys.F90	(revision 2669)
@@ -0,0 +1,132 @@
+      subroutine iophys_ecrit(nom,lllm,titre,unite,px)
+      USE dimphy
+      USE mod_phys_lmdz_para
+      USE mod_grid_phy_lmdz
+      IMPLICIT NONE
+
+
+
+!  Ecriture de variables diagnostiques au choix dans la physique 
+!  dans un fichier NetCDF nomme  'diagfi'. Ces variables peuvent etre
+!  3d (ex : temperature), 2d (ex : temperature de surface), ou
+!  0d (pour un scalaire qui ne depend que du temps : ex : la longitude
+!  solaire)
+!  (ou encore 1d, dans le cas de testphys1d, pour sortir une colonne)
+!  La periode d'ecriture est donnee par 
+!  "ecritphy " regle dans le fichier de controle de run :  run.def
+!
+!    writediagfi peut etre appele de n'importe quelle subroutine
+!    de la physique, plusieurs fois. L'initialisation et la creation du
+!    fichier se fait au tout premier appel.
+!
+! WARNING : les variables dynamique (u,v,t,q,ps)
+!  sauvees par writediagfi avec une
+! date donnee sont legerement differentes que dans le fichier histoire car 
+! on ne leur a pas encore ajoute de la dissipation et de la physique !!!
+! IL est  RECOMMANDE d'ajouter les tendance physique a ces variables
+! avant l'ecriture dans diagfi (cf. physiq.F)
+!  
+! Modifs: Aug.2010 Ehouarn: enforce outputs to be real*4
+!
+!  parametres (input) :
+!  ----------
+!      unit : unite logique du fichier de sortie (toujours la meme)
+!      nom  : nom de la variable a sortir (chaine de caracteres)
+!      titre: titre de la variable (chaine de caracteres)
+!      unite : unite de la variable (chaine de caracteres)
+!      px : variable a sortir (real 0, 1, 2, ou 3d)
+!
+!=================================================================
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "netcdf.inc"
+#include "iotd.h"
+
+
+! Arguments on input:
+      integer lllm
+      character (len=*) :: nom,titre,unite
+      integer imjmax
+      parameter (imjmax=100000)
+      real px(klon_omp,lllm)
+      real xglo(klon_glo,lllm)
+      real zx(iim,jjp1,lllm)
+
+
+      CALL Gather(px,xglo)
+!$OMP MASTER
+      IF (is_mpi_root) THEN       
+        CALL Grid1Dto2D_glo(xglo,zx)
+        call iotd_ecrit(nom,lllm,titre,unite,zx)
+      ENDIF
+!$OMP END MASTER
+
+      return
+      end
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+      SUBROUTINE iophys_ini
+      USE mod_phys_lmdz_para
+      USE vertical_layers_mod, ONLY: presnivs
+      USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
+      USE dimphy, ONLY: klev
+
+      IMPLICIT NONE
+
+!=======================================================================
+!
+!   Auteur:  L. Fairhead  ,  P. Le Van, Y. Wanherdrick, F. Forget
+!   -------
+!
+!   Objet:
+!   ------
+!
+!   'Initialize' the diagfi.nc file: write down dimensions as well
+!   as time-independent fields (e.g: geopotential, mesh area, ...)
+!
+!=======================================================================
+!-----------------------------------------------------------------------
+!   Declarations:
+!   -------------
+
+real pi
+
+!   Arguments:
+!   ----------
+
+!$OMP MASTER
+    IF (is_mpi_root) THEN       
+pi=2.*asin(1.)
+call iotd_ini('phys.nc   ', &
+!iim,jjp1,llm,rlonv(1:iim)*180./pi,rlatu*180./pi,presnivs)
+size(lon_reg),size(lat_reg),klev,lon_reg(:)*180./pi,lat_reg*180./pi,presnivs)
+    ENDIF
+!$OMP END MASTER
+
+      END
+
+#ifdef und
+      SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
+      IMPLICIT none
+
+!=======================================================================
+      INTEGER nfield,nlon,iim,jjmp1, jjm
+      REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
+
+      INTEGER i, n, ig
+
+      jjm = jjmp1 - 1
+      DO n = 1, nfield
+         DO i=1,iim
+            ecrit(i,n) = fi(1,n)
+            ecrit(i+jjm*iim,n) = fi(nlon,n)
+         ENDDO
+         DO ig = 1, nlon - 2
+           ecrit(iim+ig,n) = fi(1+ig,n)
+         ENDDO
+      ENDDO
+      RETURN
+      END
+
+#endif
Index: /LMDZ5/branches/testing/libf/phylmd/limit_slab.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/limit_slab.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/limit_slab.F90	(revision 2669)
@@ -164,8 +164,4 @@
   diff_sst(:) = diff_sst_save(:)
   diff_siv(:) = diff_siv_save(:)
-! For Debug purpose
-!  PRINT *,'limit_slab sst',MINVAL(diff_sst(:)),MAXVAL(diff_sst(:))
-!  PRINT *,'limit_slab siv',MINVAL(diff_siv(:)),MAXVAL(diff_siv(:))
-!  PRINT *,'limit_slab bils',MINVAL(lmt_bils(:)),MAXVAL(lmt_bils(:))
 
 END SUBROUTINE limit_slab
Index: /LMDZ5/branches/testing/libf/phylmd/ocean_slab_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/ocean_slab_mod.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/ocean_slab_mod.F90	(revision 2669)
@@ -9,4 +9,6 @@
   USE indice_sol_mod
   USE surface_data
+  USE mod_grid_phy_lmdz, ONLY: klon_glo
+  USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
 
   IMPLICIT NONE
@@ -14,43 +16,68 @@
   PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice, ocean_slab_ice
 
-!****************************************************************************************
+!***********************************************************************************
 ! Global saved variables
-!****************************************************************************************
-
+!***********************************************************************************
+  ! number of slab vertical layers
+  INTEGER, PUBLIC, SAVE :: nslay
+  !$OMP THREADPRIVATE(nslay)
+  ! timestep for coupling (update slab temperature) in timesteps
   INTEGER, PRIVATE, SAVE                           :: cpl_pas
   !$OMP THREADPRIVATE(cpl_pas)
+  ! cyang = 1/heat capacity of top layer (rho.c.H)
   REAL, PRIVATE, SAVE                              :: cyang
   !$OMP THREADPRIVATE(cyang)
+  ! depth of slab layers (1 or 2)
   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: slabh
   !$OMP THREADPRIVATE(slabh)
+  ! slab temperature
   REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: tslab
   !$OMP THREADPRIVATE(tslab)
+  ! heat flux convergence due to Ekman
+  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_ekman
+  !$OMP THREADPRIVATE(dt_ekman)
+  ! heat flux convergence due to horiz diffusion
+  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_hdiff
+  !$OMP THREADPRIVATE(dt_hdiff)
+  ! fraction of ocean covered by sea ice (sic / (oce+sic))
   REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: fsic
   !$OMP THREADPRIVATE(fsic)
+  ! temperature of the sea ice
   REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: tice
   !$OMP THREADPRIVATE(tice)
+  ! sea ice thickness, in kg/m2
   REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: seaice
   !$OMP THREADPRIVATE(seaice)
+  ! net surface heat flux, weighted by open ocean fraction
   REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bils
   !$OMP THREADPRIVATE(slab_bils)
+  ! slab_bils accumulated over cpl_pas timesteps
   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bils_cum
   !$OMP THREADPRIVATE(bils_cum)
+  ! net heat flux into the ocean below the ice : conduction + solar radiation
   REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bilg
   !$OMP THREADPRIVATE(slab_bilg)
+  ! slab_bilg over cpl_pas timesteps
   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bilg_cum
   !$OMP THREADPRIVATE(bilg_cum)
-
-!****************************************************************************************
+  ! wind stress saved over cpl_pas timesteps
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: taux_cum
+  !$OMP THREADPRIVATE(taux_cum)
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: tauy_cum
+  !$OMP THREADPRIVATE(tauy_cum)
+
+!***********************************************************************************
 ! Parameters (could be read in def file: move to slab_init)
-!****************************************************************************************
+!***********************************************************************************
 ! snow and ice physical characteristics:
     REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp
     REAL, PARAMETER :: t_melt=273.15   ! melting ice temp
     REAL, PARAMETER :: sno_den=300. !mean snow density, kg/m3
-    REAL, PARAMETER :: ice_den=917.
-    REAL, PARAMETER :: sea_den=1025.
-    REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity
-    REAL, PARAMETER :: sno_cond=0.31*sno_den
+    REAL, PARAMETER :: ice_den=917. ! ice density
+    REAL, PARAMETER :: sea_den=1025. ! sea water density
+    REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity of ice
+    REAL, PARAMETER :: sno_cond=0.31*sno_den ! conductivity of snow
     REAL, PARAMETER :: ice_cap=2067.   ! specific heat capacity, snow and ice
+    REAL, PARAMETER :: sea_cap=3995.   ! specific heat capacity, snow and ice
     REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice
 
@@ -74,36 +101,90 @@
     REAL, PARAMETER :: alb_ice_dry=0.75 !dry thick ice
     REAL, PARAMETER :: alb_ice_wet=0.66 !melting thick ice
-    REAL, PARAMETER :: pen_frac=0.3 !fraction of penetrating shortwave (no snow)
+    REAL, PARAMETER :: pen_frac=0.3 !fraction of shortwave penetrating into the
+    ! ice (no snow)
     REAL, PARAMETER :: pen_ext=1.5 !extinction of penetrating shortwave (m-1)
 
-!****************************************************************************************
+! horizontal transport
+   LOGICAL, PUBLIC, SAVE :: slab_hdiff
+   !$OMP THREADPRIVATE(slab_hdiff)
+   REAL, PRIVATE, SAVE    :: coef_hdiff ! coefficient for horizontal diffusion
+   !$OMP THREADPRIVATE(coef_hdiff)
+   INTEGER, PUBLIC, SAVE :: slab_ekman, slab_cadj ! Ekman, conv adjustment
+   !$OMP THREADPRIVATE(slab_ekman)
+   !$OMP THREADPRIVATE(slab_cadj)
+
+!***********************************************************************************
 
 CONTAINS
 !
-!****************************************************************************************
+!***********************************************************************************
 !
   SUBROUTINE ocean_slab_init(dtime, pctsrf_rst)
   !, seaice_rst etc
 
-    use IOIPSL
-
-    ! For ok_xxx vars (Ekman...)
-    INCLUDE "clesphys.h"
+    USE ioipsl_getin_p_mod, ONLY : getin_p
+    USE mod_phys_lmdz_transfert_para, ONLY : gather
+    USE slab_heat_transp_mod, ONLY : ini_slab_transp
 
     ! Input variables
-!****************************************************************************************
+!***********************************************************************************
     REAL, INTENT(IN)                         :: dtime
 ! Variables read from restart file
-    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst
+    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst 
+    ! surface fractions from start file
 
 ! Local variables
-!****************************************************************************************
+!************************************************************************************
     INTEGER                :: error
+    REAL, DIMENSION(klon_glo) :: zmasq_glo
     CHARACTER (len = 80)   :: abort_message
     CHARACTER (len = 20)   :: modname = 'ocean_slab_intit'
 
-!****************************************************************************************
+!***********************************************************************************
+! Define some parameters
+!***********************************************************************************
+! Number of slab layers
+    nslay=2
+    CALL getin_p('slab_layers',nslay)
+    print *,'number of slab layers : ',nslay
+! Layer thickness
+    ALLOCATE(slabh(nslay), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation slabh'
+       CALL abort_physic(modname,abort_message,1)
+    ENDIF
+    slabh(1)=50.
+    IF (nslay.GT.1) THEN 
+        slabh(2)=150.
+    END IF
+
+! cyang = 1/heat capacity of top layer (rho.c.H)
+    cyang=1/(slabh(1)*sea_den*sea_cap)
+
+! cpl_pas  coupling period (update of tslab and ice fraction)
+! pour un calcul a chaque pas de temps, cpl_pas=1
+    cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
+    CALL getin_p('cpl_pas',cpl_pas)
+    print *,'cpl_pas',cpl_pas
+
+! Horizontal diffusion
+    slab_hdiff=.FALSE.
+    CALL getin_p('slab_hdiff',slab_hdiff)
+    coef_hdiff=25000.
+    CALL getin_p('coef_hdiff',coef_hdiff)
+! Ekman transport
+    slab_ekman=0
+    CALL getin_p('slab_ekman',slab_ekman)
+! Convective adjustment
+    IF (nslay.EQ.1) THEN
+        slab_cadj=0
+    ELSE
+        slab_cadj=1
+    END IF
+    CALL getin_p('slab_cadj',slab_cadj)
+
+!************************************************************************************
 ! Allocate surface fraction read from restart file
-!****************************************************************************************
+!************************************************************************************
     ALLOCATE(fsic(klon), stat = error)
     IF (error /= 0) THEN
@@ -112,11 +193,12 @@
     ENDIF
     fsic(:)=0.
+    !zmasq = continent fraction
     WHERE (1.-zmasq(:)>EPSFRA)
         fsic(:) = pctsrf_rst(:,is_sic)/(1.-zmasq(:))
     END WHERE
 
-!****************************************************************************************
-! Allocate saved variables
-!****************************************************************************************
+!************************************************************************************
+! Allocate saved fields
+!************************************************************************************
     ALLOCATE(tslab(klon,nslay), stat=error)
        IF (error /= 0) CALL abort_physic &
@@ -136,5 +218,5 @@
     bils_cum(:) = 0.0   
 
-    IF (version_ocean=='sicINT') THEN
+    IF (version_ocean=='sicINT') THEN ! interactive sea ice
         ALLOCATE(slab_bilg(klon), stat = error)
         IF (error /= 0) THEN
@@ -161,55 +243,69 @@
     END IF
 
-!****************************************************************************************
-! Define some parameters
-!****************************************************************************************
-! Layer thickness
-    ALLOCATE(slabh(nslay), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation slabh'
-       CALL abort_physic(modname,abort_message,1)
+    IF (slab_hdiff) THEN !horizontal diffusion
+        ALLOCATE(dt_hdiff(klon,nslay), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation dt_hdiff'
+           CALL abort_physic(modname,abort_message,1)
+        ENDIF
+        dt_hdiff(:,:) = 0.0   
     ENDIF
-    slabh(1)=50.
-! cyang = 1/heat capacity of top layer (rho.c.H)
-    cyang=1/(slabh(1)*4.228e+06)
-
-! cpl_pas periode de couplage avec slab (update tslab, pctsrf)
-! pour un calcul à chaque pas de temps, cpl_pas=1
-    cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
-    CALL getin('cpl_pas',cpl_pas)
-    print *,'cpl_pas',cpl_pas
+
+    IF (slab_ekman.GT.0) THEN ! ekman transport
+        ALLOCATE(dt_ekman(klon,nslay), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation dt_ekman'
+           CALL abort_physic(modname,abort_message,1)
+        ENDIF
+        dt_ekman(:,:) = 0.0   
+        ALLOCATE(taux_cum(klon), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation taux_cum'
+           CALL abort_physic(modname,abort_message,1)
+        ENDIF
+        taux_cum(:) = 0.0   
+        ALLOCATE(tauy_cum(klon), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation tauy_cum'
+           CALL abort_physic(modname,abort_message,1)
+        ENDIF
+        tauy_cum(:) = 0.0   
+    ENDIF
+
+! Initialize transport
+    IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN
+      CALL gather(zmasq,zmasq_glo)
+!$OMP MASTER  ! Only master thread
+      IF (is_mpi_root) THEN 
+          CALL ini_slab_transp(zmasq_glo)
+      END IF
+!$OMP END MASTER
+    END IF
 
  END SUBROUTINE ocean_slab_init
 !
-!****************************************************************************************
+!***********************************************************************************
 !
   SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf_chg, is_modified)
 
-    USE limit_read_mod
-
-!    INCLUDE "clesphys.h"
+! this routine sends back the sea ice and ocean fraction to the main physics
+! routine. Called only with interactive sea ice
 
 ! Arguments
-!****************************************************************************************
-    INTEGER, INTENT(IN)                        :: itime   ! numero du pas de temps courant
-    INTEGER, INTENT(IN)                        :: jour    ! jour a lire dans l'annee
-    REAL   , INTENT(IN)                        :: dtime   ! pas de temps de la physique (en s)
+!************************************************************************************
+    INTEGER, INTENT(IN)                        :: itime   ! current timestep 
+    INTEGER, INTENT(IN)                        :: jour    !  day in year (not 
+    REAL   , INTENT(IN)                        :: dtime   ! physics timestep (s)
     REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf_chg  ! sub-surface fraction
-    LOGICAL, INTENT(OUT)                       :: is_modified ! true if pctsrf is modified at this time step
-
-! Local variables
-!****************************************************************************************
-
-    IF (version_ocean == 'sicOBS'.OR. version_ocean == 'sicNO') THEN   
-       CALL limit_read_frac(itime, dtime, jour, pctsrf_chg, is_modified)
-    ELSE
+    LOGICAL, INTENT(OUT)                       :: is_modified ! true if pctsrf is 
+                                                         ! modified at this time step
+
        pctsrf_chg(:,is_oce)=(1.-fsic(:))*(1.-zmasq(:))
        pctsrf_chg(:,is_sic)=fsic(:)*(1.-zmasq(:))
        is_modified=.TRUE.
-    END IF
 
   END SUBROUTINE ocean_slab_frac
 !
-!****************************************************************************************
+!************************************************************************************
 !
   SUBROUTINE ocean_slab_noice( & 
@@ -224,58 +320,83 @@
     
     USE calcul_fluxs_mod
+    USE slab_heat_transp_mod, ONLY: divgrad_phy, slab_ekman1, slab_ekman2
+    USE mod_phys_lmdz_para
 
     INCLUDE "clesphys.h"
 
+! This routine 
+! (1) computes surface turbulent fluxes over points with some open ocean    
+! (2) reads additional Q-flux (everywhere)
+! (3) computes horizontal transport (diffusion & Ekman)
+! (4) updates slab temperature every cpl_pas ; creates new ice if needed.
+
+! Note : 
+! klon total number of points
+! knon number of points with open ocean (varies with time)
+! knindex gives position of the knon points within klon.
+! In general, local saved variables have klon values
+! variables exchanged with PBL module have knon.
+
 ! Input arguments
-!****************************************************************************************
-    INTEGER, INTENT(IN)                  :: itime
-    INTEGER, INTENT(IN)                  :: jour
-    INTEGER, INTENT(IN)                  :: knon
-    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
-    REAL, INTENT(IN)                     :: dtime
-    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
-    REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragq, cdragm
-    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
-    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
-    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ
-    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
-    REAL, DIMENSION(klon), INTENT(IN)    :: ps
-    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness
-    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
-    REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
+!***********************************************************************************
+    INTEGER, INTENT(IN)                  :: itime ! current timestep INTEGER,
+    INTEGER, INTENT(IN)                  :: jour  ! day in year (for Q-Flux) 
+    INTEGER, INTENT(IN)                  :: knon  ! number of points 
+    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex 
+    REAL, INTENT(IN) :: dtime  ! timestep (s)
+    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay 
+    REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragq, cdragm 
+    ! drag coefficients
+    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow 
+    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum ! near surface T, q
+    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ 
+    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV 
+    ! exchange coefficients for boundary layer scheme
+    REAL, DIMENSION(klon), INTENT(IN)    :: ps  ! surface pressure
+    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness ! surface wind
+    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in ! surface temperature
+    REAL, DIMENSION(klon), INTENT(INOUT) :: radsol ! net surface radiative flux
 
 ! In/Output arguments
-!****************************************************************************************
-    REAL, DIMENSION(klon), INTENT(INOUT) :: snow
+!************************************************************************************
+    REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! in kg/m2
     
 ! Output arguments
-!****************************************************************************************
+!************************************************************************************
     REAL, DIMENSION(klon), INTENT(OUT)   :: qsurf
     REAL, DIMENSION(klon), INTENT(OUT)   :: evap, fluxsens, fluxlat
     REAL, DIMENSION(klon), INTENT(OUT)   :: flux_u1, flux_v1
-    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new
+    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new ! new surface tempearture
     REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l      
     REAL, DIMENSION(klon), INTENT(OUT)   :: qflux
 
 ! Local variables
-!****************************************************************************************
-    INTEGER               :: i,ki
-    ! surface fluxes
+!************************************************************************************
+    INTEGER               :: i,ki,k
+    REAL                  :: t_cadj
+    !  for surface heat fluxes
     REAL, DIMENSION(klon) :: cal, beta, dif_grnd
-    ! flux correction
+    ! for Q-Flux computation: d/dt SST, d/dt ice volume (kg/m2), surf fluxes 
     REAL, DIMENSION(klon) :: diff_sst, diff_siv, lmt_bils
-    ! surface wind stress
+    ! for surface wind stress
     REAL, DIMENSION(klon) :: u0, v0
     REAL, DIMENSION(klon) :: u1_lay, v1_lay
-    ! ice creation
+    ! for new ice creation
     REAL                  :: e_freeze, h_new, dfsic
-
-!****************************************************************************************
-! 1) Flux calculation
-!
-!****************************************************************************************
-    cal(:)      = 0.
-    beta(:)     = 1.
-    dif_grnd(:) = 0.
+    ! horizontal diffusion and Ekman local vars
+    ! dimension = global domain (klon_glo) instead of // subdomains 
+    REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_glo, dt_ekman_glo
+    ! dt_ekman_glo saved for diagnostic, dt_ekman_tmp used for time loop
+    REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_tmp, dt_ekman_tmp
+    REAL, DIMENSION(klon_glo,nslay) :: tslab_glo
+    REAL, DIMENSION(klon_glo) :: taux_glo,tauy_glo
+
+!****************************************************************************************
+! 1) Surface fluxes calculation
+!    
+!****************************************************************************************
+    cal(:)      = 0. ! infinite thermal inertia
+    beta(:)     = 1. ! wet surface
+    dif_grnd(:) = 0. ! no diffusion into ground
     
 ! Suppose zero surface speed
@@ -285,4 +406,5 @@
     v1_lay(:) = v1(:) - v0(:)
 
+! Compute latent & sensible fluxes
     CALL calcul_fluxs(knon, is_oce, dtime, &
          tsurf_in, p1lay, cal, beta, cdragh, cdragq, ps, &
@@ -292,12 +414,6 @@
          tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
 
-! - Flux calculation at first modele level for U and V
-    CALL calcul_flux_wind(knon, dtime, &
-         u0, v0, u1, v1, gustiness, cdragm, &
-         AcoefU, AcoefV, BcoefU, BcoefV, &
-         p1lay, temp_air, &
-         flux_u1, flux_v1)  
-
-! Accumulate total fluxes locally
+! save total cumulated heat fluxes locally
+! radiative + turbulent + melt of falling snow 
     slab_bils(:)=0.
     DO i=1,knon
@@ -306,13 +422,28 @@
                       -precip_snow(i)*ice_lat*(1.+snow_wfact*fsic(ki)))
         bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
-! Also taux, tauy, saved vars...
     END DO
 
-!****************************************************************************************
-! 2) Get global variables lmt_bils and diff_sst from file limit_slab.nc
+!  Compute surface wind stress
+    CALL calcul_flux_wind(knon, dtime, &
+         u0, v0, u1, v1, gustiness, cdragm, &
+         AcoefU, AcoefV, BcoefU, BcoefV, &
+         p1lay, temp_air, &
+         flux_u1, flux_v1)  
+
+! save cumulated wind stress
+    IF (slab_ekman.GT.0) THEN 
+      DO i=1,knon
+          ki=knindex(i)
+          taux_cum(ki)=taux_cum(ki)+flux_u1(i)*(1.-fsic(ki))/cpl_pas
+          tauy_cum(ki)=tauy_cum(ki)+flux_v1(i)*(1.-fsic(ki))/cpl_pas
+      END DO
+    ENDIF
+
+!****************************************************************************************
+! 2) Q-Flux : get global variables lmt_bils, diff_sst and diff_siv from file limit_slab.nc
 !
 !****************************************************************************************
     lmt_bils(:)=0.
-    CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) ! global pour un processus
+    CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) 
     ! lmt_bils and diff_sst,siv saved by limit_slab
     qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.-diff_siv(:)*ice_den*ice_lat/86400.
@@ -320,71 +451,160 @@
 
 !****************************************************************************************
-! 3) Recalculate new temperature
-!
+! 3) Recalculate new temperature (add Surf fluxes, Q-Flux, Ocean transport)
+!    Bring to freezing temp and make sea ice if necessary
+!  
 !***********************************************o*****************************************
     tsurf_new=tsurf_in
     IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
-        ! Compute transport
-        ! Add read QFlux and SST tendency
-        tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas 
-        ! Add cumulated surface fluxes
-        tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime
-        ! Update surface temperature
-        SELECT CASE(version_ocean)
-        CASE('sicNO')
-            DO i=1,knon
-                ki=knindex(i)
-                tsurf_new(i)=tslab(ki,1)
+! ***********************************
+!  Horizontal transport 
+! ***********************************
+      IF (slab_ekman.GT.0) THEN 
+          ! copy wind stress to global var
+          CALL gather(taux_cum,taux_glo)
+          CALL gather(tauy_cum,tauy_glo)
+      END IF
+
+      IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN 
+        CALL gather(tslab,tslab_glo)
+      ! Compute horiz transport on one process only
+!$OMP MASTER  ! Only master thread
+        IF (is_mpi_root) THEN ! Only master processus          
+          IF (slab_hdiff) THEN 
+              dt_hdiff_glo(:,:)=0.
+          END IF
+          IF (slab_ekman.GT.0) THEN 
+              dt_ekman_glo(:,:)=0.
+          END IF
+          DO i=1,cpl_pas ! time splitting for numerical stability
+            IF (slab_ekman.GT.0) THEN
+              SELECT CASE (slab_ekman)
+                CASE (1)
+                  CALL slab_ekman1(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp)
+                CASE (2)
+                  CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp)
+                CASE DEFAULT
+                  dt_ekman_tmp(:,:)=0.
+              END SELECT
+              dt_ekman_glo(:,:)=dt_ekman_glo(:,:)+dt_ekman_tmp(:,:)
+              ! convert dt_ekman from K.s-1.(kg.m-2) to K.s-1   
+              DO k=1,nslay
+                dt_ekman_tmp(:,k)=dt_ekman_tmp(:,k)/(slabh(k)*sea_den)
+              ENDDO
+              tslab_glo=tslab_glo+dt_ekman_tmp*dtime
+            ENDIF
+            IF (slab_hdiff) THEN ! horizontal diffusion
+              ! laplacian of slab T
+              CALL divgrad_phy(nslay,tslab_glo,dt_hdiff_tmp)
+              ! multiply by diff coef and normalize to 50m slab equivalent
+              dt_hdiff_tmp=dt_hdiff_tmp*coef_hdiff*50./SUM(slabh)
+              dt_hdiff_glo(:,:)=dt_hdiff_glo(:,:)+ dt_hdiff_tmp(:,:)
+              tslab_glo=tslab_glo+dt_hdiff_tmp*dtime
+            END IF
+          END DO ! time splitting
+          IF (slab_hdiff) THEN
+            !dt_hdiff_glo saved in W/m2
+            DO k=1,nslay
+              dt_hdiff_glo(:,k)=dt_hdiff_glo(:,k)*slabh(k)*sea_den*sea_cap/cpl_pas
             END DO
-        CASE('sicOBS') ! check for sea ice or tslab below freezing
-            DO i=1,knon
-                ki=knindex(i)
-                IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN
-                    tslab(ki,1)=t_freeze
-                END IF
-                tsurf_new(i)=tslab(ki,1)
-            END DO
-        CASE('sicINT')
-            DO i=1,knon
-                ki=knindex(i)
-                IF (fsic(ki).LT.epsfra) THEN ! Free of ice
-                    IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
-                        ! quantity of new ice formed
-                        e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
-                        ! new ice
-                        tice(ki)=t_freeze
-                        fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
-                        IF (fsic(ki).GT.ice_frac_min) THEN
-                            seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
-                            tslab(ki,1)=t_freeze
-                        ELSE
-                            fsic(ki)=0.
-                        END IF
-                        tsurf_new(i)=t_freeze
-                    ELSE
-                        tsurf_new(i)=tslab(ki,1)
-                    END IF 
-                ELSE ! ice present
-                    tsurf_new(i)=t_freeze
-                    IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
-                        ! quantity of new ice formed over open ocean
-                        e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
-                                 /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
-                        ! new ice height and fraction
-                        h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new
-                        dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new)
-                        h_new=MIN(e_freeze/dfsic,h_ice_max)
-                        ! update tslab to freezing over open ocean only
-                        tslab(ki,1)=tslab(ki,1)*fsic(ki)+t_freeze*(1.-fsic(ki))
-                        ! update sea ice
-                        seaice(ki)=(h_new*dfsic+seaice(ki)*fsic(ki)) &
-                                   /(dfsic+fsic(ki))
-                        fsic(ki)=fsic(ki)+dfsic
-                        ! update snow?
-                    END IF !freezing
-                END IF ! sea ice present
-            END DO
-        END SELECT
-        bils_cum(:)=0.0! clear cumulated fluxes
+          END IF
+          IF (slab_ekman.GT.0) THEN 
+            ! dt_ekman_glo saved in W/m2
+            dt_ekman_glo(:,:)=dt_ekman_glo(:,:)*sea_cap/cpl_pas
+          END IF
+        END IF ! is_mpi_root
+!$OMP END MASTER
+!$OMP BARRIER
+        ! Send new fields back to all processes
+        CALL Scatter(tslab_glo,tslab)
+        IF (slab_hdiff) THEN
+          CALL Scatter(dt_hdiff_glo,dt_hdiff)
+        END IF
+        IF (slab_ekman.GT.0) THEN 
+          CALL Scatter(dt_ekman_glo,dt_ekman)
+          ! clear wind stress
+          taux_cum(:)=0.
+          tauy_cum(:)=0.
+        END IF
+      ENDIF ! transport
+
+! ***********************************
+! Other heat fluxes
+! ***********************************
+      ! Add read QFlux
+      tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas 
+      ! Add cumulated surface fluxes
+      tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime
+      ! Convective adjustment if 2 layers 
+      IF ((nslay.GT.1).AND.(slab_cadj.GT.0)) THEN
+        DO i=1,klon
+          IF (tslab(i,2).GT.tslab(i,1)) THEN
+            ! mean (mass-weighted) temperature
+            t_cadj=SUM(tslab(i,:)*slabh(:))/SUM(slabh(:))
+            tslab(i,1)=t_cadj
+            tslab(i,2)=t_cadj
+          END IF
+        END DO
+      END IF
+! ***********************************
+! Update surface temperature and ice 
+! ***********************************
+      SELECT CASE(version_ocean)
+      CASE('sicNO') ! no sea ice even below freezing !
+          DO i=1,knon
+              ki=knindex(i)
+              tsurf_new(i)=tslab(ki,1)
+          END DO
+      CASE('sicOBS') ! "realistic" case, for prescribed sea ice
+        ! tslab cannot be below freezing, or above it if there is sea ice
+          DO i=1,knon
+              ki=knindex(i)
+              IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN
+                  tslab(ki,1)=t_freeze
+              END IF
+              tsurf_new(i)=tslab(ki,1)
+          END DO
+      CASE('sicINT') ! interactive sea ice
+          DO i=1,knon
+              ki=knindex(i)
+              IF (fsic(ki).LT.epsfra) THEN ! Free of ice
+                  IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
+                      ! quantity of new ice formed
+                      e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
+                      ! new ice
+                      tice(ki)=t_freeze
+                      fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
+                      IF (fsic(ki).GT.ice_frac_min) THEN
+                          seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
+                          tslab(ki,1)=t_freeze
+                      ELSE
+                          fsic(ki)=0.
+                      END IF
+                      tsurf_new(i)=t_freeze
+                  ELSE
+                      tsurf_new(i)=tslab(ki,1)
+                  END IF 
+              ELSE ! ice present
+                  tsurf_new(i)=t_freeze
+                  IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
+                      ! quantity of new ice formed over open ocean
+                      e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
+                               /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
+                      ! new ice height and fraction
+                      h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new
+                      dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new)
+                      h_new=MIN(e_freeze/dfsic,h_ice_max)
+                      ! update tslab to freezing over open ocean only
+                      tslab(ki,1)=tslab(ki,1)*fsic(ki)+t_freeze*(1.-fsic(ki))
+                      ! update sea ice
+                      seaice(ki)=(h_new*dfsic+seaice(ki)*fsic(ki)) &
+                                 /(dfsic+fsic(ki))
+                      fsic(ki)=fsic(ki)+dfsic
+                      ! update snow?
+                  END IF ! tslab below freezing
+              END IF ! sea ice present
+          END DO
+      END SELECT
+      bils_cum(:)=0.0! clear cumulated fluxes
     END IF ! coupling time
   END SUBROUTINE ocean_slab_noice
@@ -530,7 +750,7 @@
            seaice(ki) = MAX(0.0,seaice(ki)-(evap(i)-snow_evap)*dtime)
         ENDIF
-! Melt / Freeze from above if Tsurf>0
+! Melt / Freeze snow from above if Tsurf>0
         IF (tsurf_new(i).GT.t_melt) THEN
-            ! energy available for melting snow (in kg/m2 of snow)
+            ! energy available for melting snow (in kg of melted snow /m2)
             e_melt = MIN(MAX(snow(i)*(tsurf_new(i)-t_melt)*ice_cap/2. &
                /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i))
@@ -672,5 +892,5 @@
     END IF ! coupling time
     
-    !tests fraction glace
+    !tests ice fraction 
     WHERE (fsic.LT.ice_frac_min)
         tslab(:,1)=tslab(:,1)-fsic*seaice*ice_lat*cyang
@@ -691,9 +911,14 @@
     IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
     IF (ALLOCATED(fsic)) DEALLOCATE(fsic)
+    IF (ALLOCATED(tice)) DEALLOCATE(tice)
+    IF (ALLOCATED(seaice)) DEALLOCATE(seaice)
     IF (ALLOCATED(slab_bils)) DEALLOCATE(slab_bils)
     IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg)
     IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum)
     IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum)
-    IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
+    IF (ALLOCATED(taux_cum)) DEALLOCATE(taux_cum)
+    IF (ALLOCATED(tauy_cum)) DEALLOCATE(tauy_cum)
+    IF (ALLOCATED(dt_ekman)) DEALLOCATE(dt_ekman)
+    IF (ALLOCATED(dt_hdiff)) DEALLOCATE(dt_hdiff)
 
   END SUBROUTINE ocean_slab_final
Index: /LMDZ5/branches/testing/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/pbl_surface_mod.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/pbl_surface_mod.F90	(revision 2669)
@@ -3095,4 +3095,19 @@
                 endif
                 mfois(nsrf) = mfois(nsrf) + 1
+                ! F. Codron sensible default values for ocean and sea ice
+                IF (nsrf.EQ.is_oce) THEN
+                    tsurf(i,nsrf) = 271.35 ! Freezing sea water
+                    DO k=1,nsw
+                      alb_dir(i,k,nsrf) = 0.06 ! typical Ocean albedo
+                      alb_dif(i,k,nsrf) = 0.06
+                    END DO
+                ELSE IF (nsrf.EQ.is_sic) THEN
+                    tsurf(i,nsrf) = 273.15 ! Melting ice 
+                    DO k=1,nsw
+                      alb_dir(i,k,nsrf) = 0.3 ! thin ice
+                      alb_dif(i,k,nsrf) = 0.3
+                    END DO
+                END IF
+                ! F. Codron
              ELSE
                 ! The continents have changed. The new fraction receives the mean sum of the existent fractions
Index: /LMDZ5/branches/testing/libf/phylmd/phyetat0.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/phyetat0.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/phyetat0.F90	(revision 2669)
@@ -3,5 +3,5 @@
 SUBROUTINE phyetat0 (fichnom, clesphy0, tabcntr0)
 
-  USE dimphy, only: klon, zmasq, klev, nslay
+  USE dimphy, only: klon, zmasq, klev
   USE iophy, ONLY : init_iophy_new
   USE ocean_cpl_mod,    ONLY : ocean_cpl_init
@@ -25,5 +25,5 @@
   USE carbon_cycle_mod, ONLY : carbon_cycle_tr, carbon_cycle_cpl, co2_send
   USE indice_sol_mod, only: nbsrf, is_ter, epsfra, is_lic, is_oce, is_sic
-  USE ocean_slab_mod, ONLY: tslab, seaice, tice, ocean_slab_init
+  USE ocean_slab_mod, ONLY: nslay, tslab, seaice, tice, ocean_slab_init
   USE time_phylmdz_mod, ONLY: init_iteration, pdtphys, itau_phy
 
@@ -443,5 +443,15 @@
   IF ( type_ocean == 'slab' ) THEN
       CALL ocean_slab_init(dtime, pctsrf)
-      found=phyetat0_get(nslay,tslab,"tslab","tslab",0.)
+      IF (nslay.EQ.1) THEN
+        found=phyetat0_get(1,tslab,"tslab01","tslab",0.)
+        IF (.NOT. found) THEN
+            found=phyetat0_get(1,tslab,"tslab","tslab",0.)
+        END IF
+      ELSE
+          DO i=1,nslay
+            WRITE(str2,'(i2.2)') i
+            found=phyetat0_get(1,tslab(:,i),"tslab"//str2,"tslab",0.)  
+          END DO
+      END IF
       IF (.NOT. found) THEN 
           PRINT*, "phyetat0: Le champ <tslab> est absent"
@@ -453,6 +463,6 @@
 
       ! Sea ice variables
-      found=phyetat0_get(1,tice,"slab_tice","slab_tice",0.)
       IF (version_ocean == 'sicINT') THEN
+          found=phyetat0_get(1,tice,"slab_tice","slab_tice",0.)
           IF (.NOT. found) THEN 
               PRINT*, "phyetat0: Le champ <tice> est absent"
@@ -460,4 +470,5 @@
                   tice(:)=ftsol(:,is_sic)
           END IF 
+          found=phyetat0_get(1,seaice,"seaice","seaice",0.)
           IF (.NOT. found) THEN
               PRINT*, "phyetat0: Le champ <seaice> est absent"
Index: /LMDZ5/branches/testing/libf/phylmd/phyredem.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/phyredem.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/phyredem.F90	(revision 2669)
@@ -30,5 +30,5 @@
   USE indice_sol_mod, ONLY: nbsrf, is_oce, is_sic, is_ter, is_lic, epsfra
   USE surface_data, ONLY: type_ocean, version_ocean
-  USE ocean_slab_mod, ONLY : tslab, seaice, tice, fsic
+  USE ocean_slab_mod, ONLY : nslay, tslab, seaice, tice, fsic
   USE time_phylmdz_mod, ONLY: annee_ref, day_end, itau_phy, pdtphys
 
@@ -58,5 +58,5 @@
 
   INTEGER isoil, nsrf,isw
-  CHARACTER (len=7) :: str7
+  CHARACTER (len=2) :: str2
   CHARACTER (len=256) :: nam, lnam
   INTEGER           :: it, iiq
@@ -304,5 +304,12 @@
   ! Restart variables for Slab ocean
   IF (type_ocean == 'slab') THEN
-      CALL put_field("tslab", "Slab ocean temperature", tslab)
+      IF (nslay.EQ.1) THEN
+        CALL put_field("tslab", "Slab ocean temperature", tslab)
+      ELSE
+        DO it=1,nslay
+          WRITE(str2,'(i2.2)') it
+          CALL put_field("tslab"//str2, "Slab ocean temperature", tslab(:,it))
+        END DO
+      END IF
       IF (version_ocean == 'sicINT') THEN
           CALL put_field("seaice", "Slab seaice (kg/m2)", seaice)
Index: /LMDZ5/branches/testing/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 2669)
@@ -677,4 +677,8 @@
   TYPE(ctrl_out), SAVE :: o_slab_sic = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
     'seaice', 'Epaisseur banquise slab', 'kg/m2', (/ ('', i=1, 9) /))
+  TYPE(ctrl_out), SAVE :: o_slab_hdiff = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
+    'slab_hdiff', 'Horizontal diffusion', 'W/m2', (/ ('', i=1, 9) /))
+  TYPE(ctrl_out), SAVE :: o_slab_ekman = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
+    'slab_ekman', 'Ekman heat transport', 'W/m2', (/ ('', i=1, 9) /))
   TYPE(ctrl_out), SAVE :: o_ale_bl = ctrl_out((/ 1, 1, 1, 10, 10, 10, 11, 11, 11 /), &
     'ale_bl', 'ALE BL', 'm2/s2', (/ ('', i=1, 9) /))
Index: /LMDZ5/branches/testing/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/phys_output_write_mod.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/phys_output_write_mod.F90	(revision 2669)
@@ -24,5 +24,5 @@
     ! defined and initialised in phys_output_mod.F90
 
-    USE dimphy, only: klon, klev, klevp1, nslay
+    USE dimphy, only: klon, klev, klevp1
     USE mod_phys_lmdz_para, ONLY: is_north_pole_phy,is_south_pole_phy
     USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
@@ -89,4 +89,5 @@
          o_slab_qflux, o_tslab, o_slab_bils, &
          o_slab_bilg, o_slab_sic, o_slab_tice, &
+         o_slab_hdiff, o_slab_ekman, &
          o_weakinv, o_dthmin, o_cldtau, &
          o_cldemi, o_pr_con_l, o_pr_con_i, &
@@ -291,5 +292,6 @@
   
 
-    USE ocean_slab_mod, only: tslab, slab_bils, slab_bilg, tice, seaice
+    USE ocean_slab_mod, only: nslay, tslab, slab_bils, slab_bilg, tice, &
+        seaice, slab_ekman,slab_hdiff, dt_ekman, dt_hdiff
     USE pbl_surface_mod, only: snow
     USE indice_sol_mod, only: nbsrf
@@ -312,4 +314,7 @@
     USE phys_cal_mod, only : mth_len
 
+#ifdef CPP_RRTM
+    USE YOESW, ONLY : RSUN
+#endif
 
     IMPLICIT NONE
@@ -355,4 +360,6 @@
 #endif
     REAL, PARAMETER :: un_jour=86400.
+    INTEGER ISW
+    CHARACTER*1 ch1
 
     ! On calcul le nouveau tau:
@@ -403,4 +410,30 @@
        CALL histwrite_phy(o_contfracOR, pctsrf(:,is_ter))
        CALL histwrite_phy(o_aireTER, paire_ter)
+!
+#ifdef CPP_XIOS
+       CALL histwrite_phy("R_ecc",R_ecc)
+       CALL histwrite_phy("R_peri",R_peri)
+       CALL histwrite_phy("R_incl",R_incl)
+       CALL histwrite_phy("solaire",solaire)
+!
+#ifdef CPP_RRTM
+      IF (iflag_rrtm.EQ.1) THEN
+        DO ISW=1, NSW
+          WRITE(ch1,'(i1)') ISW
+!         zx_tmp_0d=RSUN(ISW)
+!         CALL histwrite_phy("rsun"//ch1,zx_tmp_0d)
+          CALL histwrite_phy("rsun"//ch1,RSUN(ISW))
+        ENDDO
+      ENDIF
+#endif
+!
+       CALL histwrite_phy("co2_ppm",co2_ppm)
+       CALL histwrite_phy("CH4_ppb",CH4_ppb)
+       CALL histwrite_phy("N2O_ppb",N2O_ppb)
+       CALL histwrite_phy("CFC11_ppt",CFC11_ppt)
+       CALL histwrite_phy("CFC12_ppt",CFC12_ppt)
+!
+#endif
+
 !!! Champs 2D !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 ! Simulateur AIRS
@@ -951,5 +984,5 @@
               CALL histwrite_phy(o_tslab, zx_tmp_fi2d)
           ELSE
-              CALL histwrite_phy(o_tslab, tslab)
+              CALL histwrite_phy(o_tslab, tslab(:,1:nslay))
           END IF
           IF (version_ocean=='sicINT') THEN
@@ -957,4 +990,20 @@
               CALL histwrite_phy(o_slab_tice, tice)
               CALL histwrite_phy(o_slab_sic, seaice)
+          END IF
+          IF (slab_hdiff) THEN
+            IF (nslay.EQ.1) THEN
+                zx_tmp_fi2d(:)=dt_hdiff(:,1)
+                CALL histwrite_phy(o_slab_hdiff, zx_tmp_fi2d)
+            ELSE
+                CALL histwrite_phy(o_slab_hdiff, dt_hdiff(:,1:nslay))
+            END IF
+          END IF
+          IF (slab_ekman.GT.0) THEN
+            IF (nslay.EQ.1) THEN
+                zx_tmp_fi2d(:)=dt_ekman(:,1)
+                CALL histwrite_phy(o_slab_ekman, zx_tmp_fi2d)
+            ELSE
+                CALL histwrite_phy(o_slab_ekman, dt_ekman(:,1:nslay))
+            END IF
           END IF
        ENDIF !type_ocean == force/slab
Index: /LMDZ5/branches/testing/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/physiq_mod.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/physiq_mod.F90	(revision 2669)
@@ -233,5 +233,7 @@
     USE ioipsl_getin_p_mod, ONLY : getin_p
 
+#ifndef CPP_XIOS
     USE paramLMDZ_phy_mod
+#endif
 
     USE cmp_seri_mod
@@ -582,6 +584,6 @@
     LOGICAL,SAVE :: ok_adjwk=.FALSE.
     !$OMP THREADPRIVATE(ok_adjwk) 
-    REAL,SAVE :: oliqmax=999.
-    !$OMP THREADPRIVATE(oliqmax) 
+    REAL,SAVE :: oliqmax=999.,oicemax=999.
+    !$OMP THREADPRIVATE(oliqmax,oicemax) 
     REAL, SAVE :: alp_offset
     !$OMP THREADPRIVATE(alp_offset)
@@ -983,4 +985,6 @@
     INTEGER, SAVE :: flag_aerosol 
     !$OMP THREADPRIVATE(flag_aerosol) 
+    LOGICAL, SAVE :: flag_bc_internal_mixture
+    !$OMP THREADPRIVATE(flag_bc_internal_mixture)
     LOGICAL, SAVE :: new_aod
     !$OMP THREADPRIVATE(new_aod) 
@@ -1011,8 +1015,4 @@
     ! edges of pressure intervals for ozone climatologies, in Pa, in strictly
     ! ascending order
-
-    integer, save:: co3i = 0
-    !     time index in NetCDF file of current ozone fields
-    !$OMP THREADPRIVATE(co3i) 
 
     integer ro3i
@@ -1140,5 +1140,5 @@
             ok_ade, ok_aie, ok_cdnc, aerosol_couple,  &
             flag_aerosol, flag_aerosol_strat, new_aod, &
-            bl95_b0, bl95_b1, &
+            flag_bc_internal_mixture, bl95_b0, bl95_b1, &
                                 ! nv flags pour la convection et les
                                 ! poches froides
@@ -1197,4 +1197,5 @@
        CALL getin_p('ok_adjwk',ok_adjwk)
        CALL getin_p('oliqmax',oliqmax)
+       CALL getin_p('oicemax',oicemax)
        CALL getin_p('ratqsp0',ratqsp0)
        CALL getin_p('ratqsdp',ratqsdp)
@@ -1551,5 +1552,7 @@
        WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3)
 
+#ifndef CPP_XIOS
        CALL ini_paramLMDZ_phy(dtime,nid_ctesGCM)
+#endif
 
 #endif
@@ -1906,51 +1909,47 @@
     ! Prescrire l'ozone et calculer l'albedo sur l'ocean.
     !
-    if (read_climoz >= 1) then
-       ! Ozone from a file
-       ! Update required ozone index:
-       ro3i = int((days_elapsed + jh_cur - jh_1jan) / year_len * 360.) + 1
-       if (ro3i == 361) ro3i = 360
-       ! (This should never occur, except perhaps because of roundup
-       ! error. See documentation.)
-       if (ro3i /= co3i) then
-          ! Update ozone field:
-          if (read_climoz == 1) then
-             call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
-                  press_in_edg=press_climoz, paprs=paprs, v3=wo)
-          else
-             ! read_climoz == 2
-             call regr_pr_av(ncid_climoz, (/"tro3         ", &
-                  "tro3_daylight"/), julien=ro3i, press_in_edg=press_climoz, &
-                  paprs=paprs, v3=wo)
-          end if
-          ! Convert from mole fraction of ozone to column density of ozone in a
-          ! cell, in kDU:
-          forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
-               * zmasse / dobson_u / 1e3
-          ! (By regridding ozone values for LMDZ only once every 360th of
-          ! year, we have already neglected the variation of pressure in one
-          ! 360th of year. So do not recompute "wo" at each time step even if
-          ! "zmasse" changes a little.)
-          co3i = ro3i
-       end if
-    ELSEIF (MOD(itap-1,lmt_pas) == 0) THEN
+    ! Update ozone if day change
+    IF (MOD(itap-1,lmt_pas) == 0) THEN
+      IF (read_climoz == 0) THEN
        ! Once per day, update ozone from Royer:
-
-       IF (solarlong0<-999.) then
-          ! Generic case with evolvoing season
-          zzz=real(days_elapsed+1)
-       ELSE IF (abs(solarlong0-1000.)<1.e-4) then
-          ! Particular case with annual mean insolation
-          zzz=real(90) ! could be revisited
-          IF (read_climoz/=-1) THEN
-             abort_message ='read_climoz=-1 is recommended when ' &
-                  // 'solarlong0=1000.'
-             CALL abort_physic (modname,abort_message,1)
-          ENDIF
-       ELSE
+        IF (solarlong0<-999.) then
+           ! Generic case with evolvoing season
+           zzz=real(days_elapsed+1)
+        ELSE IF (abs(solarlong0-1000.)<1.e-4) then
+           ! Particular case with annual mean insolation
+           zzz=real(90) ! could be revisited
+           IF (read_climoz/=-1) THEN
+              abort_message ='read_climoz=-1 is recommended when ' &
+                   // 'solarlong0=1000.'
+              CALL abort_physic (modname,abort_message,1)
+           ENDIF
+        ELSE
           ! Case where the season is imposed with solarlong0
           zzz=real(90) ! could be revisited
-       ENDIF
-       wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
+        ENDIF
+
+        wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
+      ELSE
+        ro3i = int((days_elapsed + jh_cur - jh_1jan) / year_len * 360.) + 1   
+        if (ro3i == 361) ro3i = 360
+        if (read_climoz == 1) then
+           call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
+                press_in_edg=press_climoz, paprs=paprs, v3=wo)
+        else
+           ! read_climoz == 2
+           call regr_pr_av(ncid_climoz, (/"tro3         ", &
+                "tro3_daylight"/), julien=ro3i, press_in_edg=press_climoz, &
+                paprs=paprs, v3=wo)
+        end if
+        ! Convert from mole fraction of ozone to column density of ozone in a
+        ! cell, in kDU:
+        forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
+             * zmasse / dobson_u / 1e3
+        ! (By regridding ozone values for LMDZ only once every 360th of
+        ! year, we have already neglected the variation of pressure in one
+        ! 360th of year. So do not recompute "wo" at each time step even if
+        ! "zmasse" changes a little.)
+     
+      ENDIF
     ENDIF
     !
@@ -3034,5 +3033,5 @@
          'lsc',abortphy,flag_inhib_tend)
     rain_num(:)=0.
-    DO k = 1, klev  ! re-evaporation de l'eau liquide nuageuse
+    DO k = 1, klev
        DO i = 1, klon
           IF (ql_seri(i,k)>oliqmax) THEN
@@ -3042,4 +3041,14 @@
        ENDDO
     ENDDO
+    IF (nqo==3) THEN
+    DO k = 1, klev
+       DO i = 1, klon
+          IF (qs_seri(i,k)>oicemax) THEN
+             rain_num(i)=rain_num(i)+(qs_seri(i,k)-oicemax)*zmasse(i,k)/pdtphys
+             qs_seri(i,k)=oicemax
+          ENDIF
+       ENDDO
+    ENDDO
+    ENDIF
 
     !---------------------------------------------------------------------------
@@ -3388,5 +3397,5 @@
                    !
                    CALL readaerosol_optic_rrtm( debut, aerosol_couple, &
-                        new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
+                        new_aod, flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
                         pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
                         tr_seri, mass_solu_aero, mass_solu_aero_pi,  &
@@ -4517,5 +4526,7 @@
 
 
+#ifndef CPP_XIOS
     CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
+#endif
 
 #endif
Index: /LMDZ5/branches/testing/libf/phylmd/readaerosol_optic.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/readaerosol_optic.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/readaerosol_optic.F90	(revision 2669)
@@ -166,7 +166,7 @@
   m_allaer_pi(:,:,id_CSSO4M_phy)  = 0.                   ! CSSO4M pre-ind
   m_allaer_pi(:,:,id_SSSSM_phy)   = sssupco_pi(:,:)      ! SSSSM pre-ind
-  m_allaer_pi(:,:,id_ASSSM_phy)   = sscoarse_pi(:,:)     ! CSSSM pre-ind
-  m_allaer_pi(:,:,id_CIDUSTM_phy) = ssacu_pi(:,:)        ! ASSSM pre-ind
-  m_allaer_pi(:,:,id_AIBCM_phy)   = cidust_pi(:,:)       ! CIDUSTM pre-ind
+  m_allaer_pi(:,:,id_CSSSM_phy)   = sscoarse_pi(:,:)     ! CSSSM pre-ind
+  m_allaer_pi(:,:,id_ASSSM_phy)   = ssacu_pi(:,:)        ! ASSSM pre-ind
+  m_allaer_pi(:,:,id_CIDUSTM_phy) = cidust_pi(:,:)       ! CIDUSTM pre-ind
   m_allaer_pi(:,:,id_AIBCM_phy)   = bcins_pi(:,:)        ! AIBCM pre-ind
   m_allaer_pi(:,:,id_AIPOMM_phy)  = pomins_pi(:,:)       ! AIPOMM pre-ind
Index: /LMDZ5/branches/testing/libf/phylmd/rrtm/aeropt_5wv_rrtm.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/rrtm/aeropt_5wv_rrtm.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/rrtm/aeropt_5wv_rrtm.F90	(revision 2669)
@@ -3,8 +3,9 @@
 !
 
-SUBROUTINE AEROPT_5WV_RRTM(&
-   pdel, m_allaer, delt, &
-   RHcl, ai, flag_aerosol, &
-   pplay, t_seri, &
+SUBROUTINE AEROPT_5WV_RRTM(  &
+   pdel, m_allaer,           &
+   RHcl, ai, flag_aerosol,   &
+   flag_bc_internal_mixture, & 
+   pplay, t_seri,            &
    tausum, tau )
 
@@ -54,8 +55,8 @@
   !
   REAL, DIMENSION(klon,klev), INTENT(in)   :: pdel
-  REAL, INTENT(in)                         :: delt
   REAL, DIMENSION(klon,klev,naero_tot), INTENT(in) :: m_allaer
   REAL, DIMENSION(klon,klev), INTENT(in)   :: RHcl     ! humidite relative ciel clair
   INTEGER,INTENT(in)                       :: flag_aerosol
+  LOGICAL,INTENT(in)                       :: flag_bc_internal_mixture
   REAL, DIMENSION(klon,klev), INTENT(in)   :: pplay
   REAL, DIMENSION(klon,klev), INTENT(in)   :: t_seri
@@ -106,8 +107,20 @@
 
   !
+  ! BC internal mixture
+  !
+  INTEGER, PARAMETER ::  nbclassbc = 5  ! Added by Rong Wang/OB for the 5 fractions
+                                       ! of BC in the soluble mode:
+                                       ! bc_content/0.001, 0.01, 0.02, 0.05, ! 0.1/
+  ! for Maxwell-Garnet internal mixture
+  ! Detailed theory can be found in R. Wang Estimation of global black carbon ! direct
+  ! radiative forcing and its uncertainty constrained by observations. J.
+  ! Geophys. Res. Atmos. Added by R. Wang and OB
+  REAL :: alpha_MG_5wv(nbre_RH,las,nbclassbc)
+
+  !
   ! Proprietes optiques
   !
-  REAL :: fact_RH(nbre_RH)
-  INTEGER :: n
+  REAL :: fact_RH(nbre_RH), BC_massfra
+  INTEGER :: n, classbc
 
 ! From here on we look at the optical parameters at 5 wavelengths:  
@@ -117,80 +130,45 @@
  DATA alpha_aers_5wv/ & 
                                 ! bc soluble 
-       7.930,7.930,7.930,7.930,7.930,7.930,     & 
-       7.930,7.930,10.893,12.618,14.550,16.613, & 
-       7.658,7.658,7.658,7.658,7.658,7.658,     & 
-       7.658,7.658,10.351,11.879,13.642,15.510, & 
-       7.195,7.195,7.195,7.195,7.195,7.195,     & 
-       7.195,7.195,9.551,10.847,12.381,13.994,  & 
-       6.736,6.736,6.736,6.736,6.736,6.736,     & 
-       6.736,6.736,8.818,9.938,11.283,12.687,   & 
-       6.277,6.277,6.277,6.277,6.277,6.277,     & 
-       6.277,6.277,8.123,9.094,10.275,11.501,   & 
+       7.930,7.930,7.930,7.930,7.930,7.930,7.930,7.930,10.893,12.618,14.550,16.613, & 
+       7.658,7.658,7.658,7.658,7.658,7.658,7.658,7.658,10.351,11.879,13.642,15.510, & 
+       7.195,7.195,7.195,7.195,7.195,7.195,7.195,7.195,9.551,10.847,12.381,13.994,  & 
+       6.736,6.736,6.736,6.736,6.736,6.736,6.736,6.736,8.818,9.938,11.283,12.687,   & 
+       6.277,6.277,6.277,6.277,6.277,6.277,6.277,6.277,8.123,9.094,10.275,11.501,   & 
                                 ! pom soluble 
-       6.676,6.676,6.676,6.676,6.710,6.934,   & 
-       7.141,7.569,8.034,8.529,9.456,10.511,  & 
-       5.109,5.109,5.109,5.109,5.189,5.535,   & 
-       5.960,6.852,8.008,9.712,12.897,19.676, & 
-       3.718,3.718,3.718,3.718,3.779,4.042,   & 
-       4.364,5.052,5.956,7.314,9.896,15.688,  & 
-       2.849,2.849,2.849,2.849,2.897,3.107,   & 
-       3.365,3.916,4.649,5.760,7.900,12.863,  & 
-       2.229,2.229,2.229,2.229,2.268,2.437,   & 
-       2.645,3.095,3.692,4.608,6.391,10.633,  & 
+       6.676,6.676,6.676,6.676,6.710,6.934,7.141,7.569,8.034,8.529,9.456,10.511,  & 
+       5.109,5.109,5.109,5.109,5.189,5.535,5.960,6.852,8.008,9.712,12.897,19.676, & 
+       3.718,3.718,3.718,3.718,3.779,4.042,4.364,5.052,5.956,7.314,9.896,15.688,  & 
+       2.849,2.849,2.849,2.849,2.897,3.107,3.365,3.916,4.649,5.760,7.900,12.863,  & 
+       2.229,2.229,2.229,2.229,2.268,2.437,2.645,3.095,3.692,4.608,6.391,10.633,  & 
                                 ! Sulfate (Accumulation) 
-       5.751,6.215,6.690,7.024,7.599,8.195,      & 
-       9.156,10.355,12.660,14.823,18.908,24.508, & 
-       4.320,4.675,5.052,5.375,5.787,6.274,      & 
-       7.066,8.083,10.088,12.003,15.697,21.133,  & 
-       3.079,3.351,3.639,3.886,4.205,4.584,      & 
-       5.206,6.019,7.648,9.234,12.391,17.220,    & 
-       2.336,2.552,2.781,2.979,3.236,3.540,      & 
-       4.046,4.711,6.056,7.388,10.093,14.313,    & 
-       1.777,1.949,2.134,2.292,2.503,2.751,      & 
-       3.166,3.712,4.828,5.949,8.264,11.922,     & 
+       5.751,6.215,6.690,7.024,7.599,8.195,9.156,10.355,12.660,14.823,18.908,24.508, & 
+       4.320,4.675,5.052,5.375,5.787,6.274,7.066,8.083,10.088,12.003,15.697,21.133,  & 
+       3.079,3.351,3.639,3.886,4.205,4.584,5.206,6.019,7.648,9.234,12.391,17.220,    & 
+       2.336,2.552,2.781,2.979,3.236,3.540,4.046,4.711,6.056,7.388,10.093,14.313,    & 
+       1.777,1.949,2.134,2.292,2.503,2.751,3.166,3.712,4.828,5.949,8.264,11.922,     & 
                                 ! Sulfate (Coarse) 
-       5.751,6.215,6.690,7.024,7.599,8.195,      & 
-       9.156,10.355,12.660,14.823,18.908,24.508, & 
-       4.320,4.675,5.052,5.375,5.787,6.274,      & 
-       7.066,8.083,10.088,12.003,15.697,21.133,  & 
-       3.079,3.351,3.639,3.886,4.205,4.584,      & 
-       5.206,6.019,7.648,9.234,12.391,17.220,    & 
-       2.336,2.552,2.781,2.979,3.236,3.540,      & 
-       4.046,4.711,6.056,7.388,10.093,14.313,    & 
-       1.777,1.949,2.134,2.292,2.503,2.751,      & 
-       3.166,3.712,4.828,5.949,8.264,11.922,     & 
+       5.751,6.215,6.690,7.024,7.599,8.195,9.156,10.355,12.660,14.823,18.908,24.508, & 
+       4.320,4.675,5.052,5.375,5.787,6.274,7.066,8.083,10.088,12.003,15.697,21.133,  & 
+       3.079,3.351,3.639,3.886,4.205,4.584,5.206,6.019,7.648,9.234,12.391,17.220,    & 
+       2.336,2.552,2.781,2.979,3.236,3.540,4.046,4.711,6.056,7.388,10.093,14.313,    & 
+       1.777,1.949,2.134,2.292,2.503,2.751,3.166,3.712,4.828,5.949,8.264,11.922,     & 
                            ! seasalt seasalt Super Coarse Soluble (SS)
-        0.218, 0.272, 0.293, 0.316, 0.343, 0.380, &
-        0.429, 0.501, 0.636, 0.755, 0.967, 1.495, & 
-        0.221, 0.275, 0.297, 0.320, 0.348, 0.383, &
-        0.432, 0.509, 0.640, 0.759, 0.972, 1.510, &
-        0.224, 0.279, 0.301, 0.324, 0.352, 0.388, &
-        0.438, 0.514, 0.647, 0.768, 0.985, 1.514, &
-        0.227, 0.282, 0.303, 0.327, 0.356, 0.392, &
-        0.441, 0.518, 0.652, 0.770, 0.987, 1.529, &
-        0.230, 0.285, 0.306, 0.330, 0.359, 0.396, &
-        0.446, 0.522, 0.656, 0.777, 0.993, 1.539, &
+        0.218, 0.272, 0.293, 0.316, 0.343, 0.380, 0.429, 0.501, 0.636, 0.755, 0.967, 1.495, & 
+        0.221, 0.275, 0.297, 0.320, 0.348, 0.383, 0.432, 0.509, 0.640, 0.759, 0.972, 1.510, &
+        0.224, 0.279, 0.301, 0.324, 0.352, 0.388, 0.438, 0.514, 0.647, 0.768, 0.985, 1.514, &
+        0.227, 0.282, 0.303, 0.327, 0.356, 0.392, 0.441, 0.518, 0.652, 0.770, 0.987, 1.529, &
+        0.230, 0.285, 0.306, 0.330, 0.359, 0.396, 0.446, 0.522, 0.656, 0.777, 0.993, 1.539, &
                            ! seasalt seasalt Coarse Soluble (CS)      
-        0.578, 0.706, 0.756, 0.809, 0.876, 0.964, &
-        1.081, 1.256, 1.577, 1.858, 2.366, 3.613, &
-        0.598, 0.725, 0.779, 0.833, 0.898, 0.990, &
-        1.109, 1.290, 1.609, 1.889, 2.398, 3.682, &
-        0.619, 0.750, 0.802, 0.857, 0.927, 1.022, &
-        1.141, 1.328, 1.648, 1.939, 2.455, 3.729, &
-        0.633, 0.767, 0.820, 0.879, 0.948, 1.044, &
-        1.167, 1.353, 1.683, 1.969, 2.491, 3.785, & 
-        0.648, 0.785, 0.838, 0.896, 0.967, 1.066, &
-        1.192, 1.381, 1.714, 2.006, 2.531, 3.836, &
+        0.578, 0.706, 0.756, 0.809, 0.876, 0.964, 1.081, 1.256, 1.577, 1.858, 2.366, 3.613, &
+        0.598, 0.725, 0.779, 0.833, 0.898, 0.990, 1.109, 1.290, 1.609, 1.889, 2.398, 3.682, &
+        0.619, 0.750, 0.802, 0.857, 0.927, 1.022, 1.141, 1.328, 1.648, 1.939, 2.455, 3.729, &
+        0.633, 0.767, 0.820, 0.879, 0.948, 1.044, 1.167, 1.353, 1.683, 1.969, 2.491, 3.785, & 
+        0.648, 0.785, 0.838, 0.896, 0.967, 1.066, 1.192, 1.381, 1.714, 2.006, 2.531, 3.836, &
                            ! seasalt seasalt Accumulation Soluble (AS)
-        4.432, 5.899, 6.505, 7.166, 7.964, 7.962, &
-        9.232,11.257,14.979,18.337,24.223,37.811, &
-        3.298, 4.569, 5.110, 5.709, 6.446, 6.268, &
-        7.396, 9.246,12.787,16.113,22.197,37.136, &
-        2.340, 3.358, 3.803, 4.303, 4.928, 4.696, &
-        5.629, 7.198,10.308,13.342,19.120,34.296, &
-        1.789, 2.626, 2.999, 3.422, 3.955, 3.730, &
-        4.519, 5.864, 8.593,11.319,16.653,31.331, &
-        1.359, 2.037, 2.343, 2.693, 3.139, 2.940, &
-        3.596, 4.729, 7.076, 9.469,14.266,28.043 /
+        4.432, 5.899, 6.505, 7.166, 7.964, 7.962, 9.232,11.257,14.979,18.337,24.223,37.811, &
+        3.298, 4.569, 5.110, 5.709, 6.446, 6.268, 7.396, 9.246,12.787,16.113,22.197,37.136, &
+        2.340, 3.358, 3.803, 4.303, 4.928, 4.696, 5.629, 7.198,10.308,13.342,19.120,34.296, &
+        1.789, 2.626, 2.999, 3.422, 3.955, 3.730, 4.519, 5.864, 8.593,11.319,16.653,31.331, &
+        1.359, 2.037, 2.343, 2.693, 3.139, 2.940, 3.596, 4.729, 7.076, 9.469,14.266,28.043 /
 
   DATA alpha_aeri_5wv/ &
@@ -201,4 +179,38 @@
                                  ! pom insoluble 
         5.042, 3.101, 1.890, 1.294, 0.934/ 
+
+  ! internal mixture
+  DATA alpha_MG_5wv/ &
+     ! bc content = 0.001
+       7.930,7.930,7.930,7.930,7.930,7.930,7.930,7.930,10.893,12.618,14.550,16.613, & 
+       7.658,7.658,7.658,7.658,7.658,7.658,7.658,7.658,10.351,11.879,13.642,15.510, & 
+       7.195,7.195,7.195,7.195,7.195,7.195,7.195,7.195,9.551,10.847,12.381,13.994,  & 
+       6.736,6.736,6.736,6.736,6.736,6.736,6.736,6.736,8.818,9.938,11.283,12.687,   & 
+       6.277,6.277,6.277,6.277,6.277,6.277,6.277,6.277,8.123,9.094,10.275,11.501,   & 
+     ! bc content = 0.01
+       7.930,7.930,7.930,7.930,7.930,7.930,7.930,7.930,10.893,12.618,14.550,16.613, & 
+       7.658,7.658,7.658,7.658,7.658,7.658,7.658,7.658,10.351,11.879,13.642,15.510, & 
+       7.195,7.195,7.195,7.195,7.195,7.195,7.195,7.195,9.551,10.847,12.381,13.994,  & 
+       6.736,6.736,6.736,6.736,6.736,6.736,6.736,6.736,8.818,9.938,11.283,12.687,   & 
+       6.277,6.277,6.277,6.277,6.277,6.277,6.277,6.277,8.123,9.094,10.275,11.501,   & 
+     ! bc content = 0.02
+       7.930,7.930,7.930,7.930,7.930,7.930,7.930,7.930,10.893,12.618,14.550,16.613, & 
+       7.658,7.658,7.658,7.658,7.658,7.658,7.658,7.658,10.351,11.879,13.642,15.510, & 
+       7.195,7.195,7.195,7.195,7.195,7.195,7.195,7.195,9.551,10.847,12.381,13.994,  & 
+       6.736,6.736,6.736,6.736,6.736,6.736,6.736,6.736,8.818,9.938,11.283,12.687,   & 
+       6.277,6.277,6.277,6.277,6.277,6.277,6.277,6.277,8.123,9.094,10.275,11.501,   & 
+     ! bc content = 0.05
+       7.930,7.930,7.930,7.930,7.930,7.930,7.930,7.930,10.893,12.618,14.550,16.613, & 
+       7.658,7.658,7.658,7.658,7.658,7.658,7.658,7.658,10.351,11.879,13.642,15.510, & 
+       7.195,7.195,7.195,7.195,7.195,7.195,7.195,7.195,9.551,10.847,12.381,13.994,  & 
+       6.736,6.736,6.736,6.736,6.736,6.736,6.736,6.736,8.818,9.938,11.283,12.687,   & 
+       6.277,6.277,6.277,6.277,6.277,6.277,6.277,6.277,8.123,9.094,10.275,11.501,   & 
+     ! bc content = 0.10
+       7.930,7.930,7.930,7.930,7.930,7.930,7.930,7.930,10.893,12.618,14.550,16.613, & 
+       7.658,7.658,7.658,7.658,7.658,7.658,7.658,7.658,10.351,11.879,13.642,15.510, & 
+       7.195,7.195,7.195,7.195,7.195,7.195,7.195,7.195,9.551,10.847,12.381,13.994,  & 
+       6.736,6.736,6.736,6.736,6.736,6.736,6.736,6.736,8.818,9.938,11.283,12.687,   & 
+       6.277,6.277,6.277,6.277,6.277,6.277,6.277,6.277,8.123,9.094,10.275,11.501 /
+
   ! 
   ! Initialisations
@@ -251,5 +263,5 @@
      aerosol_name(8) = id_ASSSM_phy
      aerosol_name(9) = id_CIDUSTM_phy
-     aerosol_name(10) = id_CSSO4M_phy
+     aerosol_name(10)= id_CSSO4M_phy
   ENDIF
 
@@ -289,9 +301,9 @@
         soluble=.TRUE.
         spsol=3
-        fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD
+        fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for AOD
     ELSEIF (aerosol_name(m).EQ.id_CSSO4M_phy) THEN
         soluble=.TRUE.
         spsol=4
-        fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD
+        fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for AOD
     ELSEIF (aerosol_name(m).EQ.id_SSSSM_phy) THEN 
         soluble=.TRUE.
@@ -325,4 +337,35 @@
       IF (soluble) THEN            ! For soluble aerosol
 
+        !--treat special case of soluble BC internal mixture
+        IF (spsol.EQ.1 .AND. flag_bc_internal_mixture) THEN
+
+          DO k=1, klev
+            DO i=1, klon
+
+             BC_massfra = m_allaer(i,k,id_ASBCM_phy)/(m_allaer(i,k,id_ASBCM_phy)+m_allaer(i,k,id_ASSO4M_phy))
+
+             IF (BC_massfra.GE.0.10) THEN
+               classbc = 5
+             ELSEIF  (BC_massfra.GE.0.05) THEN
+               classbc = 4
+             ELSEIF  (BC_massfra.GE.0.02) THEN
+               classbc = 3
+             ELSEIF  (BC_massfra.GE.0.01) THEN
+               classbc = 2
+             ELSE
+               classbc = 1
+             ENDIF
+
+              tau_ae5wv_int = alpha_MG_5wv(RH_num(i,k),la,classbc)+DELTA(i,k)* &
+                             (alpha_MG_5wv(RH_num(i,k)+1,la,classbc) - & 
+                              alpha_MG_5wv(RH_num(i,k),la,classbc))
+              tau(i,k,la,aerindex) = m_allaer(i,k,aerindex)/1.e6*zdh(i,k)*tau_ae5wv_int*fac
+              tausum(i,la,aerindex)=tausum(i,la,aerindex)+tau(i,k,la,aerindex)
+            ENDDO
+          ENDDO
+
+        !--other cases of soluble aerosols
+        ELSE
+
           DO k=1, klev
             DO i=1, klon
@@ -334,6 +377,9 @@
             ENDDO
           ENDDO
+
+        ENDIF
  
-      ELSE                         ! For insoluble aerosol
+      ! cases of insoluble aerosol
+      ELSE                         
 
         DO k=1, klev
Index: /LMDZ5/branches/testing/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90	(revision 2669)
@@ -3,8 +3,8 @@
 !
 SUBROUTINE AEROPT_6BANDS_RRTM ( &
-     pdel, m_allaer, delt, RHcl, &
+     pdel, m_allaer, RHcl, &
      tau_allaer, piz_allaer, &
      cg_allaer, m_allaer_pi, &
-     flag_aerosol, zrho )
+     flag_aerosol, flag_bc_internal_mixture, zrho )
 
   USE dimphy
@@ -28,9 +28,9 @@
   !
   REAL, DIMENSION(klon,klev),     INTENT(in)  :: pdel
-  REAL,                           INTENT(in)  :: delt
   REAL, DIMENSION(klon,klev,naero_tot),   INTENT(in)  :: m_allaer
   REAL, DIMENSION(klon,klev,naero_tot),   INTENT(in)  :: m_allaer_pi
   REAL, DIMENSION(klon,klev),     INTENT(in)  :: RHcl       ! humidite relative ciel clair
   INTEGER,                        INTENT(in)  :: flag_aerosol
+  LOGICAL,                        INTENT(in)  :: flag_bc_internal_mixture
   REAL, DIMENSION(klon,klev),     INTENT(in)  :: zrho
   !
@@ -71,7 +71,10 @@
 
   REAL, DIMENSION(klon,klev,naero_tot,nbands_sw_rrtm) ::  tau_ae
-  REAL, DIMENSION(klon,klev,naero_tot,nbands_sw_rrtm) ::  tau_ae_pi
   REAL, DIMENSION(klon,klev,naero_tot,nbands_sw_rrtm) ::  piz_ae
   REAL, DIMENSION(klon,klev,naero_tot,nbands_sw_rrtm) ::  cg_ae
+
+  REAL, DIMENSION(klon,klev,naero_tot,nbands_sw_rrtm) ::  tau_ae_pi
+  REAL, DIMENSION(klon,klev,id_ASBCM_phy:id_ASBCM_phy,nbands_sw_rrtm) :: piz_ae_pi
+  REAL, DIMENSION(klon,klev,id_ASBCM_phy:id_ASBCM_phy,nbands_sw_rrtm) :: cg_ae_pi
   !
   ! Proprietes optiques
@@ -83,7 +86,24 @@
   REAL:: piz_aers_6bands(nbre_RH,nbands_sw_rrtm,naero_soluble)     !-- unit
   REAL:: piz_aeri_6bands(nbands_sw_rrtm,naero_insoluble)        !-- unit
-
-  INTEGER :: id
-  REAL :: tmp_var, tmp_var_pi
+  ! 
+  ! BC internal mixture
+  !
+  INTEGER, PARAMETER ::  nbclassbc = 5  ! Added by Rong Wang/OB for the 5 fractions
+                                       ! of BC in the soluble mode:
+                                       ! bc_content/0.001, 0.01, 0.02, 0.05, 0.1/
+  ! for Maxwell-Garnet internal mixture
+  ! Detailed theory can be found in R. Wang Estimation of global black carbon direct
+  ! radiative forcing and its uncertainty constrained by observations. J.
+  ! Geophys. Res. Atmos. Added by R. Wang and OB
+  REAL :: alpha_MG_6bands(nbre_RH,nbands_sw_rrtm,nbclassbc)
+  REAL :: cg_MG_6bands(nbre_RH,nbands_sw_rrtm,nbclassbc)
+  REAL :: piz_MG_6bands(nbre_RH,nbands_sw_rrtm,nbclassbc)
+  !
+  INTEGER :: id, classbc, classbc_pi
+  REAL :: tmp_var, tmp_var_pi, BC_massfra, BC_massfra_pi
+
+  !
+  REAL, PARAMETER :: tau_min = 1.e-15
+!  REAL, PARAMETER :: tau_min = 1.e-7
 
 !***************************************************************************
@@ -268,4 +288,118 @@
   0.973, 0.973, 0.972, 0.940, 0.816, 0.663 /
 
+! Added by R. Wang (July 31 2016)
+! properties for BC assuming Maxwell-Garnett rule and internal mixture
+
+  DATA alpha_MG_6bands/ &
+     ! bc content = 0.001
+  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
+  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
+  5.900, 5.900, 5.900, 5.900, 5.900, 6.502, 7.151, 8.496, 9.818, 9.965,10.124,11.564, &
+  4.284, 4.284, 4.284, 4.284, 4.284, 4.721, 5.193, 6.169, 7.129, 7.236, 7.352, 8.397, &
+  2.163, 2.163, 2.163, 2.163, 2.163, 2.384, 2.622, 3.115, 3.600, 3.654, 3.712, 4.240, &
+  0.966, 0.966, 0.966, 0.966, 0.966, 1.065, 1.171, 1.392, 1.608, 1.632, 1.658, 1.894, &
+     ! bc content = 0.01
+  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
+  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
+  5.900, 5.900, 5.900, 5.900, 5.900, 6.502, 7.151, 8.496, 9.818, 9.965,10.124,11.564, &
+  4.284, 4.284, 4.284, 4.284, 4.284, 4.721, 5.193, 6.169, 7.129, 7.236, 7.352, 8.397, &
+  2.163, 2.163, 2.163, 2.163, 2.163, 2.384, 2.622, 3.115, 3.600, 3.654, 3.712, 4.240, &
+  0.966, 0.966, 0.966, 0.966, 0.966, 1.065, 1.171, 1.392, 1.608, 1.632, 1.658, 1.894, &
+     ! bc content = 0.02
+  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
+  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
+  5.900, 5.900, 5.900, 5.900, 5.900, 6.502, 7.151, 8.496, 9.818, 9.965,10.124,11.564, &
+  4.284, 4.284, 4.284, 4.284, 4.284, 4.721, 5.193, 6.169, 7.129, 7.236, 7.352, 8.397, &
+  2.163, 2.163, 2.163, 2.163, 2.163, 2.384, 2.622, 3.115, 3.600, 3.654, 3.712, 4.240, &
+  0.966, 0.966, 0.966, 0.966, 0.966, 1.065, 1.171, 1.392, 1.608, 1.632, 1.658, 1.894, &
+     ! bc content = 0.05
+  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
+  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
+  5.900, 5.900, 5.900, 5.900, 5.900, 6.502, 7.151, 8.496, 9.818, 9.965,10.124,11.564, &
+  4.284, 4.284, 4.284, 4.284, 4.284, 4.721, 5.193, 6.169, 7.129, 7.236, 7.352, 8.397, &
+  2.163, 2.163, 2.163, 2.163, 2.163, 2.384, 2.622, 3.115, 3.600, 3.654, 3.712, 4.240, &
+  0.966, 0.966, 0.966, 0.966, 0.966, 1.065, 1.171, 1.392, 1.608, 1.632, 1.658, 1.894, &
+     ! bc content = 0.10
+  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
+  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
+  5.900, 5.900, 5.900, 5.900, 5.900, 6.502, 7.151, 8.496, 9.818, 9.965,10.124,11.564, &
+  4.284, 4.284, 4.284, 4.284, 4.284, 4.721, 5.193, 6.169, 7.129, 7.236, 7.352, 8.397, &
+  2.163, 2.163, 2.163, 2.163, 2.163, 2.384, 2.622, 3.115, 3.600, 3.654, 3.712, 4.240, &
+  0.966, 0.966, 0.966, 0.966, 0.966, 1.065, 1.171, 1.392, 1.608, 1.632, 1.658, 1.894 /
+
+  DATA cg_MG_6bands/ &
+     ! bc content = 0.001
+  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
+  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
+  0.643, 0.643, 0.643, 0.654, 0.662, 0.670, 0.677, 0.688, 0.698, 0.704, 0.707, 0.715, &
+  0.513, 0.513, 0.513, 0.522, 0.530, 0.536, 0.542, 0.552, 0.560, 0.565, 0.568, 0.575, &
+  0.321, 0.321, 0.321, 0.323, 0.325, 0.327, 0.328, 0.331, 0.333, 0.334, 0.335, 0.337, &
+  0.153, 0.153, 0.153, 0.149, 0.145, 0.142, 0.139, 0.135, 0.130, 0.128, 0.127, 0.123, &
+     ! bc content = 0.01
+  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
+  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
+  0.643, 0.643, 0.643, 0.654, 0.662, 0.670, 0.677, 0.688, 0.698, 0.704, 0.707, 0.715, &
+  0.513, 0.513, 0.513, 0.522, 0.530, 0.536, 0.542, 0.552, 0.560, 0.565, 0.568, 0.575, &
+  0.321, 0.321, 0.321, 0.323, 0.325, 0.327, 0.328, 0.331, 0.333, 0.334, 0.335, 0.337, &
+  0.153, 0.153, 0.153, 0.149, 0.145, 0.142, 0.139, 0.135, 0.130, 0.128, 0.127, 0.123, &
+     ! bc content = 0.02
+  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
+  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
+  0.643, 0.643, 0.643, 0.654, 0.662, 0.670, 0.677, 0.688, 0.698, 0.704, 0.707, 0.715, &
+  0.513, 0.513, 0.513, 0.522, 0.530, 0.536, 0.542, 0.552, 0.560, 0.565, 0.568, 0.575, &
+  0.321, 0.321, 0.321, 0.323, 0.325, 0.327, 0.328, 0.331, 0.333, 0.334, 0.335, 0.337, &
+  0.153, 0.153, 0.153, 0.149, 0.145, 0.142, 0.139, 0.135, 0.130, 0.128, 0.127, 0.123, &
+     ! bc content = 0.05
+  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
+  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
+  0.643, 0.643, 0.643, 0.654, 0.662, 0.670, 0.677, 0.688, 0.698, 0.704, 0.707, 0.715, &
+  0.513, 0.513, 0.513, 0.522, 0.530, 0.536, 0.542, 0.552, 0.560, 0.565, 0.568, 0.575, &
+  0.321, 0.321, 0.321, 0.323, 0.325, 0.327, 0.328, 0.331, 0.333, 0.334, 0.335, 0.337, &
+  0.153, 0.153, 0.153, 0.149, 0.145, 0.142, 0.139, 0.135, 0.130, 0.128, 0.127, 0.123, &
+     ! bc content = 0.10
+  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
+  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
+  0.643, 0.643, 0.643, 0.654, 0.662, 0.670, 0.677, 0.688, 0.698, 0.704, 0.707, 0.715, &
+  0.513, 0.513, 0.513, 0.522, 0.530, 0.536, 0.542, 0.552, 0.560, 0.565, 0.568, 0.575, &
+  0.321, 0.321, 0.321, 0.323, 0.325, 0.327, 0.328, 0.331, 0.333, 0.334, 0.335, 0.337, &
+  0.153, 0.153, 0.153, 0.149, 0.145, 0.142, 0.139, 0.135, 0.130, 0.128, 0.127, 0.123 /
+
+  DATA piz_MG_6bands/ &
+     ! bc content = 0.001
+  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
+  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
+  0.445, 0.445, 0.445, 0.445, 0.445, 0.521, 0.583, 0.679, 0.741, 0.747, 0.753, 0.798, &
+  0.394, 0.394, 0.394, 0.394, 0.394, 0.477, 0.545, 0.649, 0.718, 0.724, 0.730, 0.779, &
+  0.267, 0.267, 0.267, 0.267, 0.267, 0.365, 0.446, 0.571, 0.652, 0.660, 0.667, 0.725, &
+  0.121, 0.121, 0.121, 0.121, 0.121, 0.139, 0.155, 0.178, 0.193, 0.195, 0.196, 0.207, &
+     ! bc content = 0.01
+  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
+  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
+  0.445, 0.445, 0.445, 0.445, 0.445, 0.521, 0.583, 0.679, 0.741, 0.747, 0.753, 0.798, &
+  0.394, 0.394, 0.394, 0.394, 0.394, 0.477, 0.545, 0.649, 0.718, 0.724, 0.730, 0.779, &
+  0.267, 0.267, 0.267, 0.267, 0.267, 0.365, 0.446, 0.571, 0.652, 0.660, 0.667, 0.725, &
+  0.121, 0.121, 0.121, 0.121, 0.121, 0.139, 0.155, 0.178, 0.193, 0.195, 0.196, 0.207, &
+     ! bc content = 0.02
+  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
+  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
+  0.445, 0.445, 0.445, 0.445, 0.445, 0.521, 0.583, 0.679, 0.741, 0.747, 0.753, 0.798, &
+  0.394, 0.394, 0.394, 0.394, 0.394, 0.477, 0.545, 0.649, 0.718, 0.724, 0.730, 0.779, &
+  0.267, 0.267, 0.267, 0.267, 0.267, 0.365, 0.446, 0.571, 0.652, 0.660, 0.667, 0.725, &
+  0.121, 0.121, 0.121, 0.121, 0.121, 0.139, 0.155, 0.178, 0.193, 0.195, 0.196, 0.207, &
+     ! bc content = 0.05
+  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
+  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
+  0.445, 0.445, 0.445, 0.445, 0.445, 0.521, 0.583, 0.679, 0.741, 0.747, 0.753, 0.798, &
+  0.394, 0.394, 0.394, 0.394, 0.394, 0.477, 0.545, 0.649, 0.718, 0.724, 0.730, 0.779, &
+  0.267, 0.267, 0.267, 0.267, 0.267, 0.365, 0.446, 0.571, 0.652, 0.660, 0.667, 0.725, &
+  0.121, 0.121, 0.121, 0.121, 0.121, 0.139, 0.155, 0.178, 0.193, 0.195, 0.196, 0.207, &
+     ! bc content = 0.10
+  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
+  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
+  0.445, 0.445, 0.445, 0.445, 0.445, 0.521, 0.583, 0.679, 0.741, 0.747, 0.753, 0.798, &
+  0.394, 0.394, 0.394, 0.394, 0.394, 0.477, 0.545, 0.649, 0.718, 0.724, 0.730, 0.779, &
+  0.267, 0.267, 0.267, 0.267, 0.267, 0.365, 0.446, 0.571, 0.652, 0.660, 0.667, 0.725, &
+  0.121, 0.121, 0.121, 0.121, 0.121, 0.139, 0.155, 0.178, 0.193, 0.195, 0.196, 0.207 /
+
 !----BEGINNING OF CALCULATIONS
 
@@ -342,7 +476,10 @@
 
   tau_ae(:,:,:,:)=0.
-  tau_ae_pi(:,:,:,:)=0.
   piz_ae(:,:,:,:)=0.
   cg_ae(:,:,:,:)=0.
+
+  tau_ae_pi(:,:,:,:)=0.
+  piz_ae_pi(:,:,:,:)=0.
+  cg_ae_pi(:,:,:,:)=0.
     
   DO m=1,nb_aer   ! tau is only computed for each mass
@@ -357,9 +494,9 @@
         soluble=.TRUE.
         spsol=3
-        fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD
+        fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for AOD
      ELSEIF  (aerosol_name(m).EQ.id_CSSO4M_phy) THEN
         soluble=.TRUE.
         spsol=4
-        fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD
+        fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for AOD
      ELSEIF (aerosol_name(m).EQ.id_SSSSM_phy) THEN 
          soluble=.TRUE.
@@ -384,35 +521,126 @@
      ENDIF
 
+    !--shortname for aerosol index
     id=aerosol_name(m)
 
     IF (soluble) THEN
 
-       DO k=1, klev
-         DO i=1, klon
-           tmp_var=m_allaer(i,k,spsol)/1.e6*zdh(i,k)*fac
-           tmp_var_pi=m_allaer_pi(i,k,spsol)/1.e6*zdh(i,k)*fac
-
-           DO inu=1,NSW
-
-             tau_ae2b_int= alpha_aers_6bands(RH_num(i,k),inu,spsol)+ & 
-                           delta(i,k)* (alpha_aers_6bands(RH_num(i,k)+1,inu,spsol) - & 
-                           alpha_aers_6bands(RH_num(i,k),inu,spsol))
+       !--here we treat the special case of soluble BC internal mixture with Maxwell-Garnett rule
+       IF (spsol.EQ.1 .AND. flag_bc_internal_mixture) THEN
+
+         DO k=1, klev
+           DO i=1, klon
+
+             tmp_var=m_allaer(i,k,spsol)/1.e6*zdh(i,k)*fac
+             tmp_var_pi=m_allaer_pi(i,k,spsol)/1.e6*zdh(i,k)*fac
+
+             ! Calculate the dry BC/(BC+SUL) mass ratio for all (natural+anthropogenic) aerosols
+             BC_massfra = m_allaer(i,k,id_ASBCM_phy)/(m_allaer(i,k,id_ASBCM_phy)+m_allaer(i,k,id_ASSO4M_phy))
+
+             IF (BC_massfra.GE.0.10) THEN
+               classbc = 5
+             ELSEIF  (BC_massfra.GE.0.05) THEN
+               classbc = 4
+             ELSEIF  (BC_massfra.GE.0.02) THEN
+               classbc = 3
+             ELSEIF  (BC_massfra.GE.0.01) THEN
+               classbc = 2
+             ELSE
+               classbc = 1
+             ENDIF
+
+             ! Calculate the dry BC/(BC+SUL) mass ratio for natural aerosols
+             BC_massfra_pi = m_allaer_pi(i,k,id_ASBCM_phy)/(m_allaer_pi(i,k,id_ASBCM_phy)+m_allaer_pi(i,k,id_ASSO4M_phy))
+
+             IF (BC_massfra_pi.GE.0.10) THEN
+               classbc_pi = 5
+             ELSEIF  (BC_massfra_pi.GE.0.05) THEN
+               classbc_pi = 4
+             ELSEIF  (BC_massfra_pi.GE.0.02) THEN
+               classbc_pi = 3
+             ELSEIF  (BC_massfra_pi.GE.0.01) THEN
+               classbc_pi = 2
+             ELSE
+               classbc_pi = 1
+             ENDIF
+
+             DO inu=1,NSW
+
+               !--all aerosols
+               tau_ae2b_int= alpha_MG_6bands(RH_num(i,k),inu,classbc)+                 &
+                             delta(i,k)* (alpha_MG_6bands(RH_num(i,k)+1,inu,classbc) - & 
+                             alpha_MG_6bands(RH_num(i,k),inu,classbc))
                    
-             piz_ae2b_int = piz_aers_6bands(RH_num(i,k),inu,spsol) + & 
-                            delta(i,k)* (piz_aers_6bands(RH_num(i,k)+1,inu,spsol) - & 
-                            piz_aers_6bands(RH_num(i,k),inu,spsol))
+               piz_ae2b_int = piz_MG_6bands(RH_num(i,k),inu,classbc) +                 &
+                              delta(i,k)* (piz_MG_6bands(RH_num(i,k)+1,inu,classbc) -  &
+                              piz_MG_6bands(RH_num(i,k),inu,classbc))
                    
-             cg_ae2b_int = cg_aers_6bands(RH_num(i,k),inu,spsol) + & 
-                           delta(i,k)* (cg_aers_6bands(RH_num(i,k)+1,inu,spsol) - & 
-                           cg_aers_6bands(RH_num(i,k),inu,spsol))
-
-             tau_ae(i,k,id,inu)    = tmp_var*tau_ae2b_int
-             tau_ae_pi(i,k,id,inu) = tmp_var_pi* tau_ae2b_int
-             piz_ae(i,k,id,inu)    = piz_ae2b_int
-             cg_ae(i,k,id,inu)     = cg_ae2b_int
-                      
+               cg_ae2b_int = cg_MG_6bands(RH_num(i,k),inu,classbc) +                   & 
+                             delta(i,k)* (cg_MG_6bands(RH_num(i,k)+1,inu,classbc) -    &
+                             cg_MG_6bands(RH_num(i,k),inu,classbc))
+
+               tau_ae(i,k,id,inu)    = tmp_var*tau_ae2b_int
+               piz_ae(i,k,id,inu)    = piz_ae2b_int
+               cg_ae(i,k,id,inu)     = cg_ae2b_int
+
+               !--natural aerosols
+               tau_ae2b_int= alpha_MG_6bands(RH_num(i,k),inu,classbc_pi)+                 &
+                             delta(i,k)* (alpha_MG_6bands(RH_num(i,k)+1,inu,classbc_pi) - & 
+                             alpha_MG_6bands(RH_num(i,k),inu,classbc_pi))
+                   
+               piz_ae2b_int = piz_MG_6bands(RH_num(i,k),inu,classbc_pi) +                 &
+                              delta(i,k)* (piz_MG_6bands(RH_num(i,k)+1,inu,classbc_pi) -  &
+                              piz_MG_6bands(RH_num(i,k),inu,classbc_pi))
+                   
+               cg_ae2b_int = cg_MG_6bands(RH_num(i,k),inu,classbc_pi) +                   &
+                             delta(i,k)* (cg_MG_6bands(RH_num(i,k)+1,inu,classbc_pi) -    &
+                             cg_MG_6bands(RH_num(i,k),inu,classbc_pi))
+
+               tau_ae_pi(i,k,id,inu) = tmp_var_pi* tau_ae2b_int
+               piz_ae_pi(i,k,id,inu) = piz_ae2b_int
+               cg_ae_pi(i,k,id,inu)  = cg_ae2b_int
+                        
+             ENDDO
            ENDDO
          ENDDO
-       ENDDO
+
+       !--else treat all other cases of soluble aerosols
+       ELSE
+
+         DO k=1, klev
+           DO i=1, klon
+             tmp_var=m_allaer(i,k,spsol)/1.e6*zdh(i,k)*fac
+             tmp_var_pi=m_allaer_pi(i,k,spsol)/1.e6*zdh(i,k)*fac
+
+             DO inu=1,NSW
+
+               tau_ae2b_int= alpha_aers_6bands(RH_num(i,k),inu,spsol)+ & 
+                             delta(i,k)* (alpha_aers_6bands(RH_num(i,k)+1,inu,spsol) - & 
+                             alpha_aers_6bands(RH_num(i,k),inu,spsol))
+                   
+               piz_ae2b_int = piz_aers_6bands(RH_num(i,k),inu,spsol) + & 
+                            delta(i,k)* (piz_aers_6bands(RH_num(i,k)+1,inu,spsol) - & 
+                              piz_aers_6bands(RH_num(i,k),inu,spsol))
+                   
+               cg_ae2b_int = cg_aers_6bands(RH_num(i,k),inu,spsol) + & 
+                             delta(i,k)* (cg_aers_6bands(RH_num(i,k)+1,inu,spsol) - & 
+                             cg_aers_6bands(RH_num(i,k),inu,spsol))
+
+               tau_ae(i,k,id,inu)    = tmp_var*tau_ae2b_int
+               tau_ae_pi(i,k,id,inu) = tmp_var_pi* tau_ae2b_int
+               piz_ae(i,k,id,inu)    = piz_ae2b_int
+               cg_ae(i,k,id,inu)     = cg_ae2b_int
+                        
+             ENDDO
+           ENDDO
+         ENDDO
+
+         !--external mixture case for soluble BC
+         IF (spsol.EQ.1) THEN
+           piz_ae_pi(:,:,id,:) = piz_ae(:,:,id,:)
+           cg_ae_pi(:,:,id,:)  = cg_ae(:,:,id,:)
+         ENDIF
+
+       ENDIF
         
      ELSE    ! For all aerosol insoluble components
@@ -443,5 +671,5 @@
      DO k=1, klev
        DO i=1, klon
-!--anthropogenic aerosol
+!--all (natural + anthropogenic) aerosol
          tau_allaer(i,k,2,inu)=tau_ae(i,k,id_ASSO4M_phy,inu)+tau_ae(i,k,id_CSSO4M_phy,inu)+ &
                                tau_ae(i,k,id_ASBCM_phy,inu)+tau_ae(i,k,id_AIBCM_phy,inu)+   &
@@ -449,32 +677,35 @@
                                tau_ae(i,k,id_ASSSM_phy,inu)+tau_ae(i,k,id_CSSSM_phy,inu)+   &
                                tau_ae(i,k,id_SSSSM_phy,inu)+ tau_ae(i,k,id_CIDUSTM_phy,inu)
-         tau_allaer(i,k,2,inu)=MAX(tau_allaer(i,k,2,inu),1e-15)
-
-         piz_allaer(i,k,2,inu)=(tau_ae(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)+ &
-                                tau_ae(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)+ &
-                                tau_ae(i,k,id_ASBCM_phy,inu)*piz_ae(i,k,id_ASBCM_phy,inu)+ &
-                                tau_ae(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)+ &
-                                tau_ae(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)+ &
-                                tau_ae(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)+ &
-                                tau_ae(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)+ &
-                                tau_ae(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)+ &
-                                tau_ae(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)+ &
+         tau_allaer(i,k,2,inu)=MAX(tau_allaer(i,k,2,inu),tau_min)
+
+         piz_allaer(i,k,2,inu)=(tau_ae(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)+   &
+                                tau_ae(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)+   &
+                                tau_ae(i,k,id_ASBCM_phy,inu)*piz_ae(i,k,id_ASBCM_phy,inu)+     &
+                                tau_ae(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)+     &
+                                tau_ae(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)+   &
+                                tau_ae(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)+   &
+                                tau_ae(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)+     &
+                                tau_ae(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)+     &
+                                tau_ae(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)+     &
                                 tau_ae(i,k,id_CIDUSTM_phy,inu)*piz_ae(i,k,id_CIDUSTM_phy,inu)) &
                                 /tau_allaer(i,k,2,inu)
-         piz_allaer(i,k,2,inu)=MAX(piz_allaer(i,k,2,inu),0.01)
+         piz_allaer(i,k,2,inu)=MIN(MAX(piz_allaer(i,k,2,inu),0.01),1.0)
+         IF (tau_allaer(i,k,2,inu).LE.tau_min) piz_allaer(i,k,2,inu)=1.0
 
          cg_allaer(i,k,2,inu)=(tau_ae(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)*cg_ae(i,k,id_ASSO4M_phy,inu)+ &
                                tau_ae(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)*cg_ae(i,k,id_CSSO4M_phy,inu)+ &
-                               tau_ae(i,k,id_ASBCM_phy,inu)*piz_ae(i,k,id_ASBCM_phy,inu)*cg_ae(i,k,id_ASBCM_phy,inu)+ &
-                               tau_ae(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)*cg_ae(i,k,id_AIBCM_phy,inu)+ &
+                               tau_ae(i,k,id_ASBCM_phy,inu)*piz_ae(i,k,id_ASBCM_phy,inu)*cg_ae(i,k,id_ASBCM_phy,inu)+    &
+                               tau_ae(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)*cg_ae(i,k,id_AIBCM_phy,inu)+    &
                                tau_ae(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)*cg_ae(i,k,id_ASPOMM_phy,inu)+ &
                                tau_ae(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)*cg_ae(i,k,id_AIPOMM_phy,inu)+ &
-                               tau_ae(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)*cg_ae(i,k,id_ASSSM_phy,inu)+ &
-                               tau_ae(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)*cg_ae(i,k,id_CSSSM_phy,inu)+ &
-                               tau_ae(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)*cg_ae(i,k,id_SSSSM_phy,inu)+ &
+                               tau_ae(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)*cg_ae(i,k,id_ASSSM_phy,inu)+    &
+                               tau_ae(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)*cg_ae(i,k,id_CSSSM_phy,inu)+    &
+                               tau_ae(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)*cg_ae(i,k,id_SSSSM_phy,inu)+    &
                                tau_ae(i,k,id_CIDUSTM_phy,inu)*piz_ae(i,k,id_CIDUSTM_phy,inu)*cg_ae(i,k,id_CIDUSTM_phy,inu))/ &
                                (tau_allaer(i,k,2,inu)*piz_allaer(i,k,2,inu))
+         cg_allaer(i,k,2,inu)=MIN(MAX(cg_allaer(i,k,2,inu),0.0),1.0)
 
 !--natural aerosol
+!--ASBCM aerosols take _pi value because of internal mixture option
          tau_allaer(i,k,1,inu)=tau_ae_pi(i,k,id_ASSO4M_phy,inu)+tau_ae_pi(i,k,id_CSSO4M_phy,inu)+ &
                                tau_ae_pi(i,k,id_ASBCM_phy,inu)+tau_ae_pi(i,k,id_AIBCM_phy,inu)+   &
@@ -482,30 +713,32 @@
                                tau_ae_pi(i,k,id_ASSSM_phy,inu)+tau_ae_pi(i,k,id_CSSSM_phy,inu)+   &
                                tau_ae_pi(i,k,id_SSSSM_phy,inu)+ tau_ae_pi(i,k,id_CIDUSTM_phy,inu)
-         tau_allaer(i,k,1,inu)=MAX(tau_allaer(i,k,1,inu),1e-15)
-
-         piz_allaer(i,k,1,inu)=(tau_ae_pi(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)+ &
-                                tau_ae_pi(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)+ &
-                                tau_ae_pi(i,k,id_ASBCM_phy,inu)*piz_ae(i,k,id_ASBCM_phy,inu)+ &
-                                tau_ae_pi(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)+ &
-                                tau_ae_pi(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)+ &
-                                tau_ae_pi(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)+ &
-                                tau_ae_pi(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)+ &
-                                tau_ae_pi(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)+ &
-                                tau_ae_pi(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)+ &
+         tau_allaer(i,k,1,inu)=MAX(tau_allaer(i,k,1,inu),tau_min)
+
+         piz_allaer(i,k,1,inu)=(tau_ae_pi(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)+   &
+                                tau_ae_pi(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)+   &
+                                tau_ae_pi(i,k,id_ASBCM_phy,inu)*piz_ae_pi(i,k,id_ASBCM_phy,inu)+  &
+                                tau_ae_pi(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)+     &
+                                tau_ae_pi(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)+   &
+                                tau_ae_pi(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)+   &
+                                tau_ae_pi(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)+     &
+                                tau_ae_pi(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)+     &
+                                tau_ae_pi(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)+     &
                                 tau_ae_pi(i,k,id_CIDUSTM_phy,inu)*piz_ae(i,k,id_CIDUSTM_phy,inu)) &
                                 /tau_allaer(i,k,1,inu)
-         piz_allaer(i,k,1,inu)=MAX(piz_allaer(i,k,1,inu),0.01)
-
-         cg_allaer(i,k,1,inu)=(tau_ae_pi(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)*cg_ae(i,k,id_ASSO4M_phy,inu)+ &
-                               tau_ae_pi(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)*cg_ae(i,k,id_CSSO4M_phy,inu)+ &
-                               tau_ae_pi(i,k,id_ASBCM_phy,inu)*piz_ae(i,k,id_ASBCM_phy,inu)*cg_ae(i,k,id_ASBCM_phy,inu)+ &
-                               tau_ae_pi(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)*cg_ae(i,k,id_AIBCM_phy,inu)+ &
-                               tau_ae_pi(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)*cg_ae(i,k,id_ASPOMM_phy,inu)+ &
-                               tau_ae_pi(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)*cg_ae(i,k,id_AIPOMM_phy,inu)+ &
-                               tau_ae_pi(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)*cg_ae(i,k,id_ASSSM_phy,inu)+ &
-                               tau_ae_pi(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)*cg_ae(i,k,id_CSSSM_phy,inu)+ &
-                               tau_ae_pi(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)*cg_ae(i,k,id_SSSSM_phy,inu)+ &
+         piz_allaer(i,k,1,inu)=MIN(MAX(piz_allaer(i,k,1,inu),0.01),1.0)
+         IF (tau_allaer(i,k,1,inu).LE.tau_min) piz_allaer(i,k,1,inu)=1.0
+
+         cg_allaer(i,k,1,inu)=(tau_ae_pi(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)*cg_ae(i,k,id_ASSO4M_phy,inu)+    &
+                               tau_ae_pi(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)*cg_ae(i,k,id_CSSO4M_phy,inu)+    &
+                               tau_ae_pi(i,k,id_ASBCM_phy,inu)*piz_ae_pi(i,k,id_ASBCM_phy,inu)*cg_ae_pi(i,k,id_ASBCM_phy,inu)+ &
+                               tau_ae_pi(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)*cg_ae(i,k,id_AIBCM_phy,inu)+       &
+                               tau_ae_pi(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)*cg_ae(i,k,id_ASPOMM_phy,inu)+    &
+                               tau_ae_pi(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)*cg_ae(i,k,id_AIPOMM_phy,inu)+    &
+                               tau_ae_pi(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)*cg_ae(i,k,id_ASSSM_phy,inu)+       &
+                               tau_ae_pi(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)*cg_ae(i,k,id_CSSSM_phy,inu)+       &
+                               tau_ae_pi(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)*cg_ae(i,k,id_SSSSM_phy,inu)+       &
                                tau_ae_pi(i,k,id_CIDUSTM_phy,inu)*piz_ae(i,k,id_CIDUSTM_phy,inu)*cg_ae(i,k,id_CIDUSTM_phy,inu))/ &
                                (tau_allaer(i,k,1,inu)*piz_allaer(i,k,1,inu))
+         cg_allaer(i,k,1,inu)=MIN(MAX(cg_allaer(i,k,1,inu),0.0),1.0)
 
         ENDDO
@@ -513,5 +746,5 @@
     ENDDO
    
-!--waveband 2 and all aerosol
+!--waveband 2 and all aerosol (third index = 2)
   inu=2
   DO i=1, klon
Index: /LMDZ5/branches/testing/libf/phylmd/rrtm/readaerosol_optic_rrtm.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/rrtm/readaerosol_optic_rrtm.F90	(revision 2668)
+++ /LMDZ5/branches/testing/libf/phylmd/rrtm/readaerosol_optic_rrtm.F90	(revision 2669)
@@ -2,5 +2,5 @@
 !
 SUBROUTINE readaerosol_optic_rrtm(debut, aerosol_couple,  &
-     new_aod, flag_aerosol, itap, rjourvrai, &
+     new_aod, flag_aerosol, flag_bc_internal_mixture, itap, rjourvrai, &
      pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
      tr_seri, mass_solu_aero, mass_solu_aero_pi, &
@@ -32,4 +32,5 @@
   LOGICAL, INTENT(IN)                      :: new_aod
   INTEGER, INTENT(IN)                      :: flag_aerosol
+  LOGICAL, INTENT(IN)                      :: flag_bc_internal_mixture
   INTEGER, INTENT(IN)                      :: itap
   REAL, INTENT(IN)                         :: rjourvrai
@@ -303,15 +304,17 @@
 !--new aerosol properties
   ! aeropt_6bands for rrtm
-  CALL aeropt_6bands_rrtm( &
-       pdel, m_allaer, pdtphys, rhcl, & 
-       tau_aero, piz_aero, cg_aero,   &
-       m_allaer_pi, flag_aerosol, &
-       zrho ) 
+  CALL aeropt_6bands_rrtm(          &
+       pdel, m_allaer, rhcl,        & 
+       tau_aero, piz_aero, cg_aero, &
+       m_allaer_pi, flag_aerosol,   &
+       flag_bc_internal_mixture, zrho ) 
 
   ! aeropt_5wv only for validation and diagnostics
-  CALL aeropt_5wv_rrtm(                    &
-       pdel, m_allaer,                &
-       pdtphys, rhcl, aerindex,       & 
-       flag_aerosol, pplay, t_seri,   &
+  CALL aeropt_5wv_rrtm(              &
+       pdel, m_allaer,               &
+       rhcl, aerindex,               & 
+       flag_aerosol,                 & 
+       flag_bc_internal_mixture,     & 
+       pplay, t_seri,                &
        tausum_aero, tau3d_aero )
 
Index: /LMDZ5/branches/testing/libf/phylmd/slab_heat_transp_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/slab_heat_transp_mod.F90	(revision 2669)
+++ /LMDZ5/branches/testing/libf/phylmd/slab_heat_transp_mod.F90	(revision 2669)
@@ -0,0 +1,757 @@
+!
+MODULE slab_heat_transp_mod
+!
+! Slab ocean : temperature tendencies due to horizontal diffusion 
+! and / or Ekman transport
+
+USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo
+IMPLICIT NONE
+
+  ! Variables copied over from dyn3d dynamics:
+  REAL,SAVE,ALLOCATABLE :: fext(:) ! Coriolis f times cell area
+  !$OMP THREADPRIVATE(fext)
+  REAL,SAVE,ALLOCATABLE :: unsairez(:) ! 1/cell area
+  !$OMP THREADPRIVATE(unsairez)
+  REAL,SAVE,ALLOCATABLE :: unsaire(:)
+  !$OMP THREADPRIVATE(unsaire)
+  REAL,SAVE,ALLOCATABLE :: cu(:) ! cell longitude dim (m)
+  !$OMP THREADPRIVATE(cu)
+  REAL,SAVE,ALLOCATABLE :: cv(:) ! cell latitude dim (m)
+  !$OMP THREADPRIVATE(cv)
+  REAL,SAVE,ALLOCATABLE :: cuvsurcv(:) ! cu/cv (v points)
+  !$OMP THREADPRIVATE(cuvsurcv)
+  REAL,SAVE,ALLOCATABLE :: cvusurcu(:) ! cv/cu (u points)
+  !$OMP THREADPRIVATE(cvusurcu)
+  REAL,SAVE,ALLOCATABLE :: aire(:) ! cell area
+  !$OMP THREADPRIVATE(aire)
+  REAL,SAVE :: apoln ! area of north pole points
+  !$OMP THREADPRIVATE(apoln)
+  REAL,SAVE :: apols ! area of south pole points
+  !$OMP THREADPRIVATE(apols)
+  REAL,SAVE,ALLOCATABLE :: aireu(:) ! area of u cells 
+  !$OMP THREADPRIVATE(aireu)
+  REAL,SAVE,ALLOCATABLE :: airev(:) ! area of v cells
+  !$OMP THREADPRIVATE(airev)
+
+  ! Local variables for horiz mass flux in slab
+  LOGICAL,SAVE :: alpha_var
+  !$OMP THREADPRIVATE(alpha_var)
+  LOGICAL,SAVE :: slab_upstream
+  !$OMP THREADPRIVATE(slab_upstream)
+
+  REAL,SAVE,ALLOCATABLE :: zmasqu(:) ! continent mask for zonal mass flux
+  !$OMP THREADPRIVATE(zmasqu)
+  REAL,SAVE,ALLOCATABLE :: zmasqv(:) ! continent mask for merid mass flux
+  !$OMP THREADPRIVATE(zmasqv)
+  REAL,SAVE,ALLOCATABLE :: unsfv(:) ! 1/f, v points
+  !$OMP THREADPRIVATE(unsfv)
+  REAL,SAVE,ALLOCATABLE :: unsev(:) ! 1/epsilon (drag)
+  !$OMP THREADPRIVATE(unsev)
+  REAL,SAVE,ALLOCATABLE :: unsfu(:) ! 1/F, u points
+  !$OMP THREADPRIVATE(unsfu)
+  REAL,SAVE,ALLOCATABLE :: unseu(:)
+  !$OMP THREADPRIVATE(unseu)
+
+  ! Routines from dyn3d, valid on global dynamics grid only:
+  PRIVATE :: gr_fi_dyn, gr_dyn_fi ! to go between 1D nd 2D horiz grid
+  PRIVATE :: gr_scal_v,gr_v_scal,gr_scal_u ! change on 2D grid U,V, T points
+  PRIVATE :: grad,diverg
+
+CONTAINS
+ 
+  SUBROUTINE ini_slab_transp_geom(ip1jm,ip1jmp1,unsairez_,fext_,unsaire_,&
+                                  cu_,cuvsurcv_,cv_,cvusurcu_, &
+                                  aire_,apoln_,apols_, &
+                                  aireu_,airev_)
+    ! number of points in lon, lat
+    IMPLICIT NONE
+    ! Routine copies some geometry variables from the dynamical core
+    ! see global vars for meaning
+    INTEGER,INTENT(IN) :: ip1jm
+    INTEGER,INTENT(IN) :: ip1jmp1
+    REAL,INTENT(IN) :: unsairez_(ip1jm)
+    REAL,INTENT(IN) :: fext_(ip1jm)
+    REAL,INTENT(IN) :: unsaire_(ip1jmp1)
+    REAL,INTENT(IN) :: cu_(ip1jmp1)
+    REAL,INTENT(IN) :: cuvsurcv_(ip1jm)
+    REAL,INTENT(IN) :: cv_(ip1jm)
+    REAL,INTENT(IN) :: cvusurcu_(ip1jmp1)
+    REAL,INTENT(IN) :: aire_(ip1jmp1)
+    REAL,INTENT(IN) :: apoln_
+    REAL,INTENT(IN) :: apols_
+    REAL,INTENT(IN) :: aireu_(ip1jmp1)
+    REAL,INTENT(IN) :: airev_(ip1jm)
+
+
+    ! Sanity check on dimensions
+    if ((ip1jm.ne.((nbp_lon+1)*(nbp_lat-1))).or. &
+        (ip1jmp1.ne.((nbp_lon+1)*nbp_lat))) then
+      write(*,*) "ini_slab_transp_geom Error: wrong array sizes"
+      stop
+    endif
+! Allocations could be done only on master process/thread...
+    allocate(unsairez(ip1jm))
+    unsairez(:)=unsairez_(:)
+    allocate(fext(ip1jm))
+    fext(:)=fext_(:)
+    allocate(unsaire(ip1jmp1))
+    unsaire(:)=unsaire_(:)
+    allocate(cu(ip1jmp1))
+    cu(:)=cu_(:)
+    allocate(cuvsurcv(ip1jm))
+    cuvsurcv(:)=cuvsurcv_(:)
+    allocate(cv(ip1jm))
+    cv(:)=cv_(:)
+    allocate(cvusurcu(ip1jmp1))
+    cvusurcu(:)=cvusurcu_(:)
+    allocate(aire(ip1jmp1))
+    aire(:)=aire_(:)
+    apoln=apoln_
+    apols=apols_
+    allocate(aireu(ip1jmp1))
+    aireu(:)=aireu_(:)
+    allocate(airev(ip1jm))
+    airev(:)=airev_(:)
+
+  END SUBROUTINE ini_slab_transp_geom
+
+  SUBROUTINE ini_slab_transp(zmasq)
+
+!    USE ioipsl_getin_p_mod, only: getin_p
+    USE IOIPSL, ONLY : getin 
+    IMPLICIT NONE
+
+    REAL zmasq(klon_glo) ! ocean / continent mask, 1=continent
+    REAL zmasq_2d((nbp_lon+1)*nbp_lat)
+    REAL ff((nbp_lon+1)*(nbp_lat-1)) ! Coriolis parameter
+    REAL eps ! epsilon friction timescale (s-1)
+    INTEGER :: slab_ekman 
+    INTEGER i
+    INTEGER :: iim,iip1,jjp1,ip1jm,ip1jmp1
+
+! Some definition for grid size
+    ip1jm=(nbp_lon+1)*(nbp_lat-1)
+    ip1jmp1=(nbp_lon+1)*nbp_lat
+    iim=nbp_lon
+    iip1=nbp_lon+1
+    jjp1=nbp_lat
+    ip1jm=(nbp_lon+1)*(nbp_lat-1)
+    ip1jmp1=(nbp_lon+1)*nbp_lat
+
+! Define ocean / continent mask (no flux into continent cell)
+    allocate(zmasqu(ip1jmp1))
+    allocate(zmasqv(ip1jm))
+    zmasqu=1.
+    zmasqv=1.
+
+    ! convert mask to 2D grid
+    CALL gr_fi_dyn(1,iip1,jjp1,zmasq,zmasq_2d)
+    ! put flux mask to 0 at boundaries of continent cells
+    DO i=1,ip1jmp1-1
+      IF (zmasq_2d(i).GT.1e-5 .OR. zmasq_2d(i+1).GT.1e-5) THEN
+              zmasqu(i)=0.
+      ENDIF
+    END DO
+    DO i=iip1,ip1jmp1,iip1
+      zmasqu(i)=zmasqu(i-iim)
+    END DO
+    DO i=1,ip1jm
+      IF (zmasq_2d(i).GT.1e-5 .OR. zmasq_2d(i+iip1).GT.1e-5) THEN
+              zmasqv(i)=0.
+      END IF
+    END DO
+    slab_ekman=2
+    CALL getin("slab_ekman",slab_ekman)
+! Coriolis and friction for Ekman transport
+    IF (slab_ekman.GT.0) THEN
+      allocate(unsfv(ip1jm))
+      allocate(unsev(ip1jm))
+      allocate(unsfu(ip1jmp1))
+      allocate(unseu(ip1jmp1))
+
+      eps=1e-5 ! Drag
+      CALL getin('slab_eps',eps)
+      print *,'epsilon=',eps
+      ff=fext*unsairez ! Coriolis
+      ! coefs to convert tau_x, tau_y to Ekman mass fluxes
+      ! on 2D grid v points
+      DO i=1,ip1jm
+        unsev(i)=eps/(ff(i)*ff(i)+eps**2)
+        unsfv(i)=ff(i)/(ff(i)*ff(i)+eps**2)
+      ENDDO
+      ! compute values on 2D u grid
+      CALL gr_v_scal(1,unsfv,unsfu)
+      CALL gr_v_scal(1,unsev,unseu)
+    END IF
+
+! Other options for Ekman transport
+    ! Alpha variable?
+      alpha_var=.FALSE.
+      CALL getin('slab_alphav',alpha_var)
+      print *,'alpha variable',alpha_var
+!  centered ou upstream scheme for meridional transport
+      slab_upstream=.FALSE.
+      CALL getin('slab_upstream',slab_upstream)
+      print *,'upstream slab scheme',slab_upstream
+  
+  END SUBROUTINE ini_slab_transp
+
+  SUBROUTINE divgrad_phy(nlevs,temp,delta)
+! Computes temperature tendency due to horizontal diffusion :
+! T Laplacian, later multiplied by diffusion coef and time-step
+
+    IMPLICIT NONE
+ 
+    INTEGER, INTENT(IN) :: nlevs ! nlevs : slab layers
+    REAL, INTENT(IN)   :: temp(klon_glo,nlevs) ! slab temperature
+    REAL , INTENT(OUT) :: delta(klon_glo,nlevs) ! temp laplacian (heat flux div.)
+    REAL :: delta_2d((nbp_lon+1)*nbp_lat,nlevs)
+    REAL ghx((nbp_lon+1)*nbp_lat,nlevs), ghy((nbp_lon+1)*(nbp_lat-1),nlevs)
+    INTEGER :: ll,iip1,jjp1
+ 
+    iip1=nbp_lon+1
+    jjp1=nbp_lat
+    
+    ! transpose temp to 2D horiz. grid
+    CALL gr_fi_dyn(nlevs,iip1,jjp1,temp,delta_2d)
+    ! computes gradient (proportional to heat flx)
+    CALL grad(nlevs,delta_2d,ghx,ghy)
+    ! put flux to 0 at ocean / continent boundary
+    DO ll=1,nlevs
+        ghx(:,ll)=ghx(:,ll)*zmasqu
+        ghy(:,ll)=ghy(:,ll)*zmasqv
+    END DO
+    ! flux divergence
+    CALL diverg(nlevs,ghx,ghy,delta_2d)
+    ! laplacian back to 1D grid
+    CALL gr_dyn_fi(nlevs,iip1,jjp1,delta_2d,delta)
+ 
+    RETURN
+  END SUBROUTINE divgrad_phy
+
+  SUBROUTINE slab_ekman1(tx_phy,ty_phy,ts_phy,dt_phy)
+! 1.5 Layer Ekman transport temperature tendency
+! note : tendency dt later multiplied by (delta t)/(rho.H)
+! to convert from divergence of heat fluxes to T
+
+      IMPLICIT NONE
+      
+      ! tx, ty : wind stress (different grids)
+      ! fluxm, fluz : mass *or heat* fluxes
+      ! dt : temperature tendency
+      INTEGER ij
+
+      ! ts surface temp, td deep temp (diagnosed)
+      REAL ts_phy(klon_glo)
+      REAL ts((nbp_lon+1)*nbp_lat),td((nbp_lon+1)*nbp_lat)
+      ! zonal and meridional wind stress 
+      REAL tx_phy(klon_glo),ty_phy(klon_glo)
+      REAL tyu((nbp_lon+1)*nbp_lat),txu((nbp_lon+1)*nbp_lat)
+      REAL txv((nbp_lon+1)*(nbp_lat-1)),tyv((nbp_lon+1)*(nbp_lat-1))
+      ! zonal and meridional Ekman mass fluxes at u, v points (2D grid)
+      REAL fluxz((nbp_lon+1)*nbp_lat),fluxm((nbp_lon+1)*(nbp_lat-1))
+      ! vertical  and absolute mass fluxes (to estimate alpha)
+      REAL fluxv((nbp_lon+1)*nbp_lat),fluxt((nbp_lon+1)*(nbp_lat-1))
+      ! temperature tendency
+      REAL dt((nbp_lon+1)*nbp_lat),dt_phy(klon_glo)
+      REAL alpha((nbp_lon+1)*nbp_lat) ! deep temperature coef
+
+      INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1
+
+! Grid definitions
+      iim=nbp_lon
+      iip1=nbp_lon+1
+      iip2=nbp_lon+2
+      jjp1=nbp_lat
+      ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
+      ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
+      ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
+
+! Convert taux,y to 2D  scalar grid
+! Note: 2D grid size = iim*jjm. iip1=iim+1
+! First and last points in zonal direction are the same 
+! we use 1 index ij from 1 to (iim+1)*(jjm+1)
+      CALL gr_fi_dyn(1,iip1,jjp1,tx_phy,txu)
+      CALL gr_fi_dyn(1,iip1,jjp1,ty_phy,tyu)
+! convert to u,v grid (Arakawa C)
+! Multiply by f or eps to get mass flux
+      ! Meridional mass flux
+      CALL gr_scal_v(1,txu,txv) ! wind stress at v points
+      CALL gr_scal_v(1,tyu,tyv)
+      fluxm=tyv*unsev-txv*unsfv ! in kg.s-1.m-1 (zonal distance)
+      ! Zonal mass flux
+      CALL gr_scal_u(1,txu,txu) ! wind stress at u points
+      CALL gr_scal_u(1,tyu,tyu)
+      fluxz=tyu*unsfu+txu*unseu
+            
+! Correct flux: continent mask and horiz grid size
+      ! multiply m-flux by mask and dx: flux in kg.s-1
+      fluxm=fluxm*cv*cuvsurcv*zmasqv 
+      ! multiply z-flux by mask and dy: flux in kg.s-1
+      fluxz=fluxz*cu*cvusurcu*zmasqu
+
+! Compute vertical  and absolute mass flux (for variable alpha)
+      IF (alpha_var) THEN
+        DO ij=iip2,ip1jm
+        fluxv(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1)
+        fluxt(ij)=ABS(fluxz(ij))+ABS(fluxz(ij-1)) &
+               +ABS(fluxm(ij))+ABS(fluxm(ij-iip1))
+        ENDDO
+        DO ij=iip1,ip1jmi1,iip1
+            fluxt(ij+1)=fluxt(ij+iip1)
+            fluxv(ij+1)=fluxv(ij+iip1)
+        END DO
+        fluxt(1)=SUM(ABS(fluxm(1:iim)))
+        fluxt(ip1jmp1)=SUM(ABS(fluxm(ip1jm-iim:ip1jm-1)))
+        fluxv(1)=-SUM(fluxm(1:iim))
+        fluxv(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1))
+        fluxt=MAX(fluxt,1.e-10)
+      ENDIF
+
+! Compute alpha coefficient. 
+! Tdeep = Tsurf * alpha + 271.35 * (1-alpha)
+      IF (alpha_var) THEN
+          ! increase alpha (and Tdeep) in downwelling regions
+          ! and decrease in upwelling regions
+          ! to avoid "hot spots" where there is surface convergence
+          DO ij=iip2,ip1jm
+              alpha(ij)=2./3.-fluxv(ij)/fluxt(ij)/3.
+          ENDDO
+          alpha(1:iip1)=2./3.-fluxv(1)/fluxt(1)/3.
+          alpha(ip1jm+1:ip1jmp1)=2./3.-fluxv(ip1jmp1)/fluxt(ip1jmp1)/3.
+      ELSE
+          alpha(:)=2./3.
+          ! Tsurf-Tdeep ~ 10° in the Tropics
+      ENDIF
+
+! Estimate deep temperature
+      CALL gr_fi_dyn(1,iip1,jjp1,ts_phy,ts)
+      DO ij=1,ip1jmp1
+         td(ij)=271.35+(ts(ij)-271.35)*alpha(ij)
+         td(ij)=MIN(td(ij),ts(ij))
+      END DO
+       
+! Meridional heat flux: multiply mass flux by (ts-td)
+! flux in K.kg.s-1
+      IF (slab_upstream) THEN
+        ! upstream scheme to avoid hot spots
+        DO ij=1,ip1jm
+        IF (fluxm(ij).GE.0.) THEN
+           fluxm(ij)=fluxm(ij)*(ts(ij+iip1)-td(ij))
+        ELSE
+           fluxm(ij)=fluxm(ij)*(ts(ij)-td(ij+iip1))
+        END IF
+        END DO
+      ELSE
+        ! centered scheme better in mid-latitudes
+        DO ij=1,ip1jm
+          fluxm(ij)=fluxm(ij)*(ts(ij+iip1)+ts(ij)-td(ij)-td(ij+iip1))/2.
+        END DO
+      ENDIF
+
+! Zonal heat flux
+! upstream scheme
+      DO ij=iip2,ip1jm
+          fluxz(ij)=fluxz(ij)*(ts(ij)+ts(ij+1)-td(ij+1)-td(ij))/2.
+      END DO
+      DO ij=iip1*2,ip1jmp1,iip1
+           fluxz(ij)=fluxz(ij-iim)
+      END DO
+
+! temperature tendency = divergence of heat fluxes
+! dt in K.s-1.kg.m-2 (T trend times mass/horiz surface)
+      DO ij=iip2,ip1jm
+         dt(ij)=(fluxz(ij-1)-fluxz(ij)+fluxm(ij)-fluxm(ij-iip1)) &
+               /aire(ij) ! aire : grid area 
+      END DO
+      DO ij=iip1,ip1jmi1,iip1
+         dt(ij+1)=dt(ij+iip1) 
+      END DO
+! special treatment at the Poles
+      dt(1)=SUM(fluxm(1:iim))/apoln
+      dt(ip1jmp1)=-SUM(fluxm(ip1jm-iim:ip1jm-1))/apols
+      dt(2:iip1)=dt(1)
+      dt(ip1jm+1:ip1jmp1)=dt(ip1jmp1)
+
+! tendencies back to 1D grid
+      CALL gr_dyn_fi(1,iip1,jjp1,dt,dt_phy)
+
+      RETURN
+  END SUBROUTINE slab_ekman1
+
+  SUBROUTINE slab_ekman2(tx_phy,ty_phy,ts_phy,dt_phy)
+! Temperature tendency for 2-layers slab ocean
+! note : tendency dt later multiplied by (delta time)/(rho.H)
+! to convert from divergence of heat fluxes to T
+      
+      IMPLICIT NONE
+      
+      ! Here, temperature and flux variables are on 2 layers
+      INTEGER ij
+
+      REAL tx_phy(klon_glo),ty_phy(klon_glo)
+      REAL txv((nbp_lon+1)*(nbp_lat-1)), tyv((nbp_lon+1)*(nbp_lat-1))
+      REAL tyu((nbp_lon+1)*nbp_lat),txu((nbp_lon+1)*nbp_lat)
+      ! slab temperature on  1D, 2D grid
+      REAL ts_phy(klon_glo,2), ts((nbp_lon+1)*nbp_lat,2)
+      ! zonal and meridional Ekman mass flux at u, v points (2D grid)
+      REAL fluxz((nbp_lon+1)*nbp_lat), fluxm((nbp_lon+1)*(nbp_lat-1))
+      ! vertical mass flux between the 2 layers
+      REAL fluxv((nbp_lon+1)*nbp_lat)
+      ! zonal and meridional heat fluxes
+      REAL fluxtz((nbp_lon+1)*nbp_lat,2)
+      REAL fluxtm((nbp_lon+1)*(nbp_lat-1),2)
+      ! temperature tendency (in K.s-1.kg.m-2)
+      REAL dt((nbp_lon+1)*nbp_lat,2), dt_phy(klon_glo,2)
+
+      INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1
+
+! Grid definitions
+      iim=nbp_lon
+      iip1=nbp_lon+1
+      iip2=nbp_lon+2
+      jjp1=nbp_lat
+      ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
+      ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
+      ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
+
+! Convert taux,y to 2D  scalar grid
+      CALL gr_fi_dyn(1,iip1,jjp1,tx_phy,txu)
+      CALL gr_fi_dyn(1,iip1,jjp1,ty_phy,tyu)
+! Multiply taux,y by f or eps, and convert to 2D u,v grids
+! (Arakawa C grid)
+      ! Meridional flux
+      CALL gr_scal_v(1,txu,txv) ! wind stress at v points
+      CALL gr_scal_v(1,tyu,tyv)
+      fluxm=tyv*unsev-txv*unsfv
+      ! Zonal flux
+      CALL gr_scal_u(1,txu,txu) ! wind stress at u points
+      CALL gr_scal_u(1,tyu,tyu)
+      fluxz=tyu*unsfu+txu*unseu
+            
+! Convert temperature to 2D grid 
+      CALL gr_fi_dyn(2,iip1,jjp1,ts_phy,ts)
+       
+!  Mass fluxes (apply continent mask, multiply by horiz grid dimension) 
+      fluxm=fluxm*cv*cuvsurcv*zmasqv
+      fluxz=fluxz*cu*cvusurcu*zmasqu
+!  Vertical mass flux from mass budget (divergence of horiz fluxes)
+      DO ij=iip2,ip1jm
+        fluxv(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1)
+      ENDDO
+      DO ij=iip1,ip1jmi1,iip1
+         fluxv(ij+1)=fluxv(ij+iip1)
+      END DO
+!  vertical mass flux at Poles
+      fluxv(1)=-SUM(fluxm(1:iim))     
+      fluxv(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1))
+      fluxv=fluxv
+
+! Meridional heat fluxes 
+      DO ij=1,ip1jm
+          ! centered scheme
+          fluxtm(ij,1)=fluxm(ij)*(ts(ij+iip1,1)+ts(ij,1))/2.
+          fluxtm(ij,2)=-fluxm(ij)*(ts(ij+iip1,2)+ts(ij,2))/2.
+      END DO
+
+! Zonal heat fluxes
+! Schema upstream      
+      DO ij=iip2,ip1jm
+      IF (fluxz(ij).GE.0.) THEN
+             fluxtz(ij,1)=fluxz(ij)*ts(ij,1)
+             fluxtz(ij,2)=-fluxz(ij)*ts(ij+1,2)
+      ELSE
+             fluxtz(ij,1)=fluxz(ij)*ts(ij+1,1)
+             fluxtz(ij,2)=-fluxz(ij)*ts(ij,2)
+      ENDIF
+      ENDDO
+      DO ij=iip1*2,ip1jmp1,iip1
+             fluxtz(ij,:)=fluxtz(ij-iim,:)
+      END DO
+                   
+! Temperature tendency :
+      DO ij=iip2,ip1jm
+! divergence of horizontal heat fluxes
+         dt(ij,:)=fluxtz(ij-1,:)-fluxtz(ij,:) &
+                 +fluxtm(ij,:)-fluxtm(ij-iip1,:)
+! + vertical heat flux (mass flux * T, upstream scheme)
+         IF (fluxv(ij).GT.0.) THEN
+           dt(ij,1)=dt(ij,1)+fluxv(ij)*ts(ij,2)
+           dt(ij,2)=dt(ij,2)-fluxv(ij)*ts(ij,2)
+         ELSE
+           dt(ij,1)=dt(ij,1)+fluxv(ij)*ts(ij,1)
+           dt(ij,2)=dt(ij,2)-fluxv(ij)*ts(ij,1)
+         ENDIF
+         ! divide by cell area
+         dt(ij,:)=dt(ij,:)/aire(ij)
+      END DO
+      DO ij=iip1,ip1jmi1,iip1
+         dt(ij+1,:)=dt(ij+iip1,:) 
+      END DO
+! Pôles
+      dt(1,:)=SUM(fluxtm(1:iim,:),dim=1)
+        IF (fluxv(1).GT.0.) THEN
+          dt(1,1)=dt(1,1)+fluxv(1)*ts(1,2)
+          dt(1,2)=dt(1,2)-fluxv(1)*ts(1,2)
+        ELSE
+          dt(1,1)=dt(1,1)+fluxv(1)*ts(1,1)
+          dt(1,2)=dt(1,2)-fluxv(1)*ts(1,1)
+        ENDIF
+      dt(1,:)=dt(1,:)/apoln
+      dt(ip1jmp1,:)=-SUM(fluxtm(ip1jm-iim:ip1jm-1,:),dim=1)
+       IF (fluxv(ip1jmp1).GT.0.) THEN
+         dt(ip1jmp1,1)=dt(ip1jmp1,1)+fluxv(ip1jmp1)*ts(ip1jmp1,2)
+         dt(ip1jmp1,2)=dt(ip1jmp1,2)-fluxv(ip1jmp1)*ts(ip1jmp1,2)
+       ELSE
+         dt(ip1jmp1,1)=dt(ip1jmp1,1)+fluxv(ip1jmp1)*ts(ip1jmp1,1)
+         dt(ip1jmp1,2)=dt(ip1jmp1,2)-fluxv(ip1jmp1)*ts(ip1jmp1,1)
+       ENDIF
+      dt(ip1jmp1,:)=dt(ip1jmp1,:)/apols
+      dt(2:iip1,1)=dt(1,1)
+      dt(2:iip1,2)=dt(1,2)
+      dt(ip1jm+1:ip1jmp1,1)=dt(ip1jmp1,1)
+      dt(ip1jm+1:ip1jmp1,2)=dt(ip1jmp1,2)
+
+! T tendency back to 1D grid...
+      CALL gr_dyn_fi(2,iip1,jjp1,dt,dt_phy)
+
+      RETURN
+  END SUBROUTINE slab_ekman2
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  SUBROUTINE gr_fi_dyn(nfield,im,jm,pfi,pdyn)
+  ! Transfer a variable from 1D "physics" grid to 2D "dynamics" grid
+  IMPLICIT NONE
+
+  INTEGER,INTENT(IN) :: im,jm,nfield
+  REAL,INTENT(IN) :: pfi(klon_glo,nfield) ! on 1D grid
+  REAL,INTENT(OUT) :: pdyn(im,jm,nfield) ! on 2D grid
+
+  INTEGER :: i,j,ifield,ig
+
+  DO ifield=1,nfield
+    ! Handle poles
+    DO i=1,im
+      pdyn(i,1,ifield)=pfi(1,ifield)
+      pdyn(i,jm,ifield)=pfi(klon_glo,ifield)
+    ENDDO
+    ! Other points
+    DO j=2,jm-1
+      ig=2+(j-2)*(im-1)
+      CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
+      pdyn(im,j,ifield)=pdyn(1,j,ifield)
+    ENDDO
+  ENDDO ! of DO ifield=1,nfield
+
+  END SUBROUTINE gr_fi_dyn
+
+  SUBROUTINE gr_dyn_fi(nfield,im,jm,pdyn,pfi)
+  ! Transfer a variable from 2D "dynamics" grid to 1D "physics" grid
+  IMPLICIT NONE
+
+  INTEGER,INTENT(IN) :: im,jm,nfield
+  REAL,INTENT(IN) :: pdyn(im,jm,nfield) ! on 2D grid
+  REAL,INTENT(OUT) :: pfi(klon_glo,nfield) ! on 1D grid
+
+  INTEGER j,ifield,ig
+
+  ! Sanity check:
+  IF(klon_glo.NE.2+(jm-2)*(im-1)) THEN
+    WRITE(*,*) "gr_dyn_fi error, wrong sizes"
+    STOP
+  ENDIF
+
+  ! Handle poles
+  CALL SCOPY(nfield,pdyn,im*jm,pfi,klon_glo)
+  CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(klon_glo,1),klon_glo)
+  ! Other points
+  DO ifield=1,nfield
+    DO j=2,jm-1
+      ig=2+(j-2)*(im-1)
+      CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
+    ENDDO
+  ENDDO
+
+  END SUBROUTINE gr_dyn_fi
+
+  SUBROUTINE  grad(klevel,pg,pgx,pgy)
+  ! compute the covariant components pgx,pgy of the gradient of pg
+  ! pgx = d(pg)/dx * delta(x) = delta(pg)
+  IMPLICIT NONE
+
+  INTEGER,INTENT(IN) :: klevel
+  REAL,INTENT(IN) :: pg((nbp_lon+1)*nbp_lat,klevel)
+  REAL,INTENT(OUT) :: pgx((nbp_lon+1)*nbp_lat,klevel)
+  REAL,INTENT(OUT) :: pgy((nbp_lon+1)*(nbp_lat-1),klevel)
+
+  INTEGER :: l,ij
+  INTEGER :: iim,iip1,ip1jm,ip1jmp1
+
+  iim=nbp_lon
+  iip1=nbp_lon+1
+  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
+  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
+
+  DO l=1,klevel
+    DO ij=1,ip1jmp1-1
+      pgx(ij,l)=pg(ij+1,l)-pg(ij,l)
+    ENDDO
+    ! correction for pgx(ip1,j,l) ...
+    ! ... pgx(iip1,j,l)=pgx(1,j,l) ...
+    DO ij=iip1,ip1jmp1,iip1
+      pgx(ij,l)=pgx(ij-iim,l)
+    ENDDO
+    DO ij=1,ip1jm
+      pgy(ij,l)=pg(ij,l)-pg(ij+iip1,l)
+    ENDDO
+  ENDDO
+
+  END SUBROUTINE grad
+
+  SUBROUTINE diverg(klevel,x,y,div)
+  ! computes the divergence of a vector field of components
+  ! x,y. x and y being covariant components
+  IMPLICIT NONE
+
+  INTEGER,INTENT(IN) :: klevel
+  REAL,INTENT(IN) :: x((nbp_lon+1)*nbp_lat,klevel)
+  REAL,INTENT(IN) :: y((nbp_lon+1)*(nbp_lat-1),klevel)
+  REAL,INTENT(OUT) :: div((nbp_lon+1)*nbp_lat,klevel)
+
+  INTEGER :: l,ij
+  INTEGER :: iim,iip1,iip2,ip1jm,ip1jmp1,ip1jmi1
+
+  REAL :: aiy1(nbp_lon+1),aiy2(nbp_lon+1)
+  REAL :: sumypn,sumyps
+  REAL,EXTERNAL :: SSUM
+
+  iim=nbp_lon
+  iip1=nbp_lon+1
+  iip2=nbp_lon+2
+  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
+  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
+  ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
+
+  DO l=1,klevel
+    DO ij=iip2,ip1jm-1
+      div(ij+1,l)= &
+        cvusurcu(ij+1)*x(ij+1,l)-cvusurcu(ij)*x(ij,l)+ &
+        cuvsurcv(ij-iim)*y(ij-iim,l)-cuvsurcv(ij+1)*y(ij+1,l)
+    ENDDO
+    ! correction for div(1,j,l) ...
+    ! ... div(1,j,l)= div(iip1,j,l) ...
+    DO ij=iip2,ip1jm,iip1
+      div(ij,l)=div(ij+iim,l)
+    ENDDO
+    ! at the poles
+    DO ij=1,iim
+      aiy1(ij)=cuvsurcv(ij)*y(ij,l)
+      aiy2(ij)=cuvsurcv(ij+ip1jmi1)*y(ij+ip1jmi1,l)
+    ENDDO
+    sumypn=SSUM(iim,aiy1,1)/apoln
+    sumyps=SSUM(iim,aiy2,1)/apols
+    DO ij=1,iip1
+      div(ij,l)=-sumypn
+      div(ij+ip1jm,l)=sumyps
+    ENDDO
+    ! End (poles)
+  ENDDO ! of DO l=1,klevel
+
+  !!! CALL filtreg( div, jjp1, klevel, 2, 2, .TRUE., 1 )
+  DO l=1,klevel
+    DO ij=iip2,ip1jm
+      div(ij,l)=div(ij,l)*unsaire(ij)
+    ENDDO
+  ENDDO
+
+  END SUBROUTINE diverg
+
+  SUBROUTINE gr_v_scal(nx,x_v,x_scal)
+  ! convert values from v points to scalar points on C-grid
+  ! used to  compute unsfu, unseu (u points, but depends only on latitude)
+  IMPLICIT NONE
+
+  INTEGER,INTENT(IN) :: nx ! number of levels or fields
+  REAL,INTENT(IN) :: x_v((nbp_lon+1)*(nbp_lat-1),nx)
+  REAL,INTENT(OUT) :: x_scal((nbp_lon+1)*nbp_lat,nx)
+
+  INTEGER :: l,ij
+  INTEGER :: iip1,iip2,ip1jm,ip1jmp1
+
+  iip1=nbp_lon+1
+  iip2=nbp_lon+2
+  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
+  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
+
+  DO l=1,nx
+    DO ij=iip2,ip1jm
+      x_scal(ij,l)= &
+                   (airev(ij-iip1)*x_v(ij-iip1,l)+airev(ij)*x_v(ij,l)) &
+                  /(airev(ij-iip1)+airev(ij))
+    ENDDO
+    DO ij=1,iip1
+      x_scal(ij,l)=0.
+    ENDDO
+    DO ij=ip1jm+1,ip1jmp1
+      x_scal(ij,l)=0.
+    ENDDO
+  ENDDO
+
+  END SUBROUTINE gr_v_scal
+ 
+  SUBROUTINE gr_scal_v(nx,x_scal,x_v)
+  ! convert values from scalar points to v points on C-grid
+  ! used to compute wind stress at V points
+  IMPLICIT NONE
+
+  INTEGER,INTENT(IN) :: nx ! number of levels or fields
+  REAL,INTENT(OUT) :: x_v((nbp_lon+1)*(nbp_lat-1),nx)
+  REAL,INTENT(IN) :: x_scal((nbp_lon+1)*nbp_lat,nx)
+
+  INTEGER :: l,ij
+  INTEGER :: iip1,ip1jm
+
+  iip1=nbp_lon+1
+  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
+
+      DO l=1,nx
+        DO ij=1,ip1jm
+          x_v(ij,l)= &
+            (cu(ij)*cvusurcu(ij)*x_scal(ij,l)+ &
+            cu(ij+iip1)*cvusurcu(ij+iip1)*x_scal(ij+iip1,l)) &
+            /(cu(ij)*cvusurcu(ij)+cu(ij+iip1)*cvusurcu(ij+iip1))
+        ENDDO
+      ENDDO
+
+  END SUBROUTINE gr_scal_v
+
+  SUBROUTINE gr_scal_u(nx,x_scal,x_u)
+  ! convert values from scalar points to U points on C-grid
+  ! used to compute wind stress at U points
+  IMPLICIT NONE
+
+  INTEGER,INTENT(IN) :: nx
+  REAL,INTENT(OUT) :: x_u((nbp_lon+1)*nbp_lat,nx)
+  REAL,INTENT(IN) :: x_scal((nbp_lon+1)*nbp_lat,nx)
+
+  INTEGER :: l,ij
+  INTEGER :: iip1,jjp1,ip1jmp1
+
+  iip1=nbp_lon+1
+  jjp1=nbp_lat
+  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
+
+  DO l=1,nx
+     DO ij=1,ip1jmp1-1
+        x_u(ij,l)= &
+         (aire(ij)*x_scal(ij,l)+aire(ij+1)*x_scal(ij+1,l)) &
+         /(aire(ij)+aire(ij+1))
+     ENDDO
+  ENDDO
+
+  CALL SCOPY(nx*jjp1,x_u(1,1),iip1,x_u(iip1,1),iip1)
+
+  END SUBROUTINE gr_scal_u
+
+END MODULE slab_heat_transp_mod
Index: /LMDZ5/branches/testing/tools/job_ada.sh
===================================================================
--- /LMDZ5/branches/testing/tools/job_ada.sh	(revision 2669)
+++ /LMDZ5/branches/testing/tools/job_ada.sh	(revision 2669)
@@ -0,0 +1,57 @@
+#!/bin/bash
+
+
+##############################
+# Gestion des chemins d'acces
+##############################
+localdir=`pwd`
+dir_tools=`echo $0 | sed -e s:\`basename $0\`::`
+modipsl=`cd $dir_tools ; ( pwd | sed -e 's:modipsl.*.$:modipsl:' )`
+if [ ! -d $modipsl/modeles/LMDZ5 ]  ; then
+    echo Dans $0
+    echo Probleme pour indentifier l arborescence modipsl
+    exit
+fi
+
+
+##############################
+# Lien local pour l'executable
+##############################
+if [ -f gcm.e ] ; then cp gcm.e gcm.$$ ; fi
+ln -sf $3 gcm.e
+
+
+###################
+# Creation d'un job
+###################
+cat <<eod>| tmp
+# @ job_type = parallel
+# @ job_name = NOMSIMU
+# @ output   = \$(job_name).\$(jobid)
+# @ error    = \$(job_name).\$(jobid)
+# @ total_tasks = $1
+# @ parallel_threads = $2
+# @ wall_clock_limit = 00:30:00
+# @ queue
+
+ulimit -s unlimited
+export OMP_STACKSIZE=800M
+export LD_LIBRARY_PATH=\$LD_LIBRARY_PATH:/smplocal/pub/NetCDF/4.1.3/lib:/smplocal/pub/HDF5/1.8.9/seq/lib
+
+cd $localdir
+\rm -f hist*
+
+poe  ./gcm.e
+
+for type in day hf mth ; do
+   \rm -f \$type.nc
+   $modipsl/bin/rebuild -o \$type.nc hist\${type}_*nc
+done
+
+eod
+
+###################
+# soumission du job
+###################
+llsubmit tmp
+
Index: /LMDZ5/branches/testing/tools/run_para.sh
===================================================================
--- /LMDZ5/branches/testing/tools/run_para.sh	(revision 2669)
+++ /LMDZ5/branches/testing/tools/run_para.sh	(revision 2669)
@@ -0,0 +1,37 @@
+#!/bin/bash
+
+export OMP_STACKSIZE=800M  # Max stack size for open MP
+export OMP_NUM_THREADS=$2  # Number of openMP processes (threads)
+ulimit -s unlimited        # Allowing illimted memory for ...
+mkdir -p WORK$$ ; mv *_0*nc WORK$$ # saving old files
+
+dir_tools=`echo $0 | sed -e s:\`basename $0\`::`
+modipsl=`cd $dir_tools ; ( pwd | sed -e 's:modipsl.*.$:modipsl:' )`
+if [ ! -d $modipsl/modeles/LMDZ5 ]  ; then
+   echo Dans $0
+   echo Probleme pour indentifier l arborescence modipsl
+   exit
+fi
+
+# On recupere le chemin pour mpirun dans arch.fcm
+mpirun=`head -1 $modipsl/modeles/LMDZ5/arch.fcm | awk ' { print $2 } ' | sed -e 's/mpif90/mpirun/'`
+$mpirun  -np $1 $3
+
+#############################################################################
+# Reconstructing netcdf files from fragments comming from parallel execution
+# file_0000.nc file_0001.nc ... > file.nc
+# Using the rebuild commend available on ..../modipsl/bin/rebuild
+#############################################################################
+
+for f in `ls *_0000.nc` ; do
+   file=`echo $f | sed -e 's/_0000.nc//'` ; echo Rebuild for $file
+   if [ -f $file.nc ] ; then
+       echo The $file.nc file already exists
+       echo remove it and run 
+       echo $modipsl/bin/rebuild -o $file.nc ${file}_0*.nc
+       echo manually
+   else
+      $modipsl/bin/rebuild -o $file.nc ${file}_0*.nc
+      if [ -f $file.nc ] ; then \rm -f ${file}_0*.nc ; fi
+   fi
+done
