Index: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/CONFIG/SHARED/iodef_medusa_et_al.xml
===================================================================
--- branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/CONFIG/SHARED/iodef_medusa_et_al.xml	(revision 7720)
+++ branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/CONFIG/SHARED/iodef_medusa_et_al.xml	(revision 7766)
@@ -1,9 +1,9 @@
 <!-- GO6.0 version of the XIOS namelist -->
-<!-- XIOS2 version namelist -->
-
+<!-- XIOS version 2 -->
 <?xml version="1.0"?>
+
 <simulation> 
 
- <context id="nemo" time_origin="1950-01-01 00:00:00" >
+ <context id="nemo">
     
     <!-- $id$ -->
@@ -15,6 +15,6 @@
 ============================================================================================================
     -->
+    <field_definition src="./field_def_bgc.xml"/>
     <field_definition src="./field_def.xml"/>
-    <field_definition src="./field_def_bgc.xml"/>
     <!-- 
 ============================================================================================================
@@ -24,6 +24,6 @@
 ============================================================================================================
     -->
-    
-    <file_definition type="multiple_file" name="@expname@_@freq@_@startdate@_@enddate@" sync_freq="10d" min_digits="4">
+
+    <file_definition type="multiple_file" name="@expname@_@freq@_@startdate@_@enddate@" sync_freq="1mo" min_digits="4">
     
       <file_group id="1ts" output_freq="1ts"  output_level="10" enabled=".TRUE."/> <!-- 1 time step files -->
@@ -33,41 +33,128 @@
       <file_group id="3h" output_freq="3h"  output_level="10" enabled=".TRUE."/> <!-- 3h files -->     
       <file_group id="4h" output_freq="4h"  output_level="10" enabled=".TRUE."/> <!-- 4h files -->
-      <file_group id="6h" output_freq="6h"  output_level="10" enabled=".FALSE."> <!-- 6h files -->
-        <file id="file8" name_suffix="_UK_shelf" description="Forcings for the UK shelf seas model" >
-          <field field_ref="toce"         name="votemper"     grid_ref="grid_T_3D_zoom"    operation="instant" standard_name="sea_water_potential_temperature"/> 
-          <field field_ref="soce"         name="vosaline"     grid_ref="grid_T_3D_zoom"    operation="instant" standard_name="sea_water_salinity"/>   
-          <field field_ref="ssh"          name="sossheig"     grid_ref="grid_T_2D_zoom"    operation="instant" standard_name="sea_surface_height_above_geoid"/>
-          <field field_ref="uoce"         name="vozocrtx"     grid_ref="grid_U_3D_zoom"    operation="instant" standard_name="sea_water_x_velocity"/>
-          <field field_ref="voce"         name="vomecrty"     grid_ref="grid_V_3D_zoom"    operation="instant" standard_name="sea_water_y_velocity"/>
-          <field field_ref="ubar"         name="ubar"         grid_ref="grid_U_3D_zoom"    operation="instant" standard_name="ocean_barotropic_current_along_i-axis"/>
-          <field field_ref="vbar"         name="vbar"         grid_ref="grid_V_3D_zoom"    operation="instant" standard_name="ocean_barotropic_current_along_j-axis"/>
-        </file>
-      </file_group>
+      <file_group id="6h" output_freq="6h"  output_level="10" enabled=".TRUE."/> <!-- 6h files -->
      
       <file_group id="1d" output_freq="1d"  output_level="10" enabled=".TRUE."/> <!-- 1d files -->
       <file_group id="3d" output_freq="3d"  output_level="10" enabled=".TRUE."/> <!-- 3d files -->    
 
-      <file_group id="10d" output_freq="10d"  output_level="10" enabled=".TRUE.">  <!-- 10d files -->   
-    
+      <file_group id="10d" output_freq="10d" output_level="10" enabled=".TRUE."> <!-- 10d files -->
+
 	<file id="file1" name_suffix="_grid_T" description="ocean T grid variables" >
-          <field_group group_ref="SPIN-UP_T" />
+          <field field_ref="e3t"          name="e3t"        standard_name="cell_thickness" />
+          <field field_ref="toce"         name="votemper"        operation="average" freq_op="1mo" > @toce_e3t / @e3t </field>
+          <field field_ref="toce"         name="votemper2"       operation="average" freq_op="1mo" > @toce2_e3t / @e3t </field>
+          <field field_ref="soce"         name="vosaline"        operation="average" freq_op="1mo" > @soce_e3t / @e3t </field>
+          <field field_ref="soce"         name="vosaline2"       operation="average" freq_op="1mo" > @soce2_e3t / @e3t </field>
+          <field field_ref="ttrd_totad"   name="opottempadvect"  operation="average" freq_op="1mo" > @ttrd_totad_e3t * 1026.0 * 3991.86795711963 </field>
+          <field field_ref="ttrd_iso"     name="opottemppmdiff"  operation="average" freq_op="1mo" > @ttrd_iso_e3t * 1026.0 * 3991.86795711963 </field>
+          <field field_ref="ttrd_zdfp"    name="opottempdiff"    operation="average" freq_op="1mo" > @ttrd_zdfp_e3t * 1026.0 * 3991.86795711963 </field>
+          <field field_ref="ttrd_evd"     name="ttrd_evd"        operation="average" freq_op="1mo" > @ttrd_evd_e3t / @e3t </field>
+          <field field_ref="ttrd_qns"     name="ttrd_qns"        operation="average" freq_op="1mo" > @ttrd_qns_e3t / @e3t_surf </field>
+          <field field_ref="ttrd_qsr"     name="rsdoabsorb"      operation="average" freq_op="1mo" > @ttrd_qsr_e3t * 1026.0 * 3991.86795711963 </field>
+          <field field_ref="ttrd_bbl"     name="ttrd_bbl"        operation="average" freq_op="1mo" > @ttrd_bbl_e3t / @e3t </field>
+          <field field_ref="ttrd_tot"     name="opottemptend"    operation="average" freq_op="1mo" > @ttrd_tot_e3t * 1026.0 * 3991.86795711963 </field>
+          <field field_ref="strd_totad"   name="osaltadvect"     operation="average" freq_op="1mo" > @strd_totad_e3t * 1026.0 * 0.001 </field>
+          <field field_ref="strd_iso"     name="osaltpmdiff"     operation="average" freq_op="1mo" > @strd_iso_e3t * 1026.0 * 0.001  </field>
+          <field field_ref="strd_zdfp"    name="osaltdiff"       operation="average" freq_op="1mo" > @strd_zdfp_e3t * 1026.0 * 0.001  </field>
+          <field field_ref="strd_evd"     name="strd_evd"        operation="average" freq_op="1mo" > @strd_evd_e3t / @e3t </field>
+          <field field_ref="strd_tot"     name="osalttend"       operation="average" freq_op="1mo" > @strd_tot_e3t * 1026.0 * 0.001 </field>
+          <field field_ref="strd_bbl"     name="strd_bbl"        operation="average" freq_op="1mo" > @strd_bbl_e3t / @e3t </field>
+	  <field field_ref="ssh"          name="sossheig"   standard_name="sea_surface_height_above_geoid"                          />
+          <field field_ref="ssh2"         name="zossq"      standard_name="square_of_sea_surface_height_above_geoid"                />
+	  <field field_ref="empmr"        name="sowaflup"   standard_name="water_flux_out_of_sea_ice_and_sea_water"                 />
+	  <field field_ref="saltflx"      name="sosafldo"   standard_name="salt_flux_into_sea_water"                                />
+	  <field field_ref="qt"           name="sohefldo"   standard_name="surface_downward_heat_flux_in_sea_water"                 />
+	  <field field_ref="qsr"          name="soshfldo"   standard_name="downwelling_photosynthetic_radiative_flux_in_sea_water"               />
+          <field field_ref="mldkz5"       name="somixhgt"   standard_name="ocean_mixed_layer_thickness_defined_by_vertical_tracer_diffusivity"   />
+          <field field_ref="mldr10_1"     name="somxl010"   standard_name="ocean_mixed_layer_thickness_defined_by_sigma_theta"      />
+          <field field_ref="mldzint_1"    name="sokaraml"   standard_name="ocean_mixed_layer_thickness"                             />
+          <field field_ref="mldzint_2"    name="somlddbm"   standard_name="ocean_mixed_layer_thickness"                             />
+          <field field_ref="ice_cover"    name="soicecov"   standard_name="sea_ice_area_fraction"                                   />
+          <field field_ref="wspd"         name="sowindsp"   standard_name="wind_speed"                                              />          
+          <field field_ref="qhcisf"       name="sohflisf"   standard_name="ice_shelf_heat_flux"                                                        />
+          <field field_ref="qlatisf"      name="soqlatisf"  standard_name="ice_shelf_latent_heat_flux"                                                        />
+          <field field_ref="fwfisf"       name="sowflisf"   standard_name="ice_shelf_melting"                                                        />
+          <field field_ref="berg_total_melt"       name="berg_total_melt"          standard_name=""                                 />
+          <field field_ref="berg_total_heat_flux"  name="berg_total_heat_flux"     standard_name=""                                 />
+          <field field_ref="runoffs"       name="sorunoff"    standard_name="water_flux_into_sea_water_from_rivers"                 />
+          <field field_ref="hflx_rain_cea" name="hfrainds"    standard_name="temperature_flux_due_to_rainfall_expressed_as_heat_flux_into_sea_water"/>
+          <field field_ref="hflx_evap_cea" name="hfevapds"    standard_name="temperature_flux_due_to_evaporation_expressed_as_heat_flux_out_of_sea_water"/>
+          <field field_ref="evap_ao_cea"   name="evap_ao_cea" standard_name="water_evaporation_flux"                                />
+          <field field_ref="rain"          name="pr"          standard_name="rainfall_flux"                                         />
+          <field field_ref="snowpre"       name="snowpre"     standard_name="snowfall_flux"                                         />
+          <field field_ref="snow_ao_cea"   name="snow_ao_cea" standard_name="snowfall_flux_over_ocean"                              />
+          <field field_ref="snow_ai_cea"   name="snow_ai_cea" standard_name="snowfall_flux_over_sea_ice"                            />
 	</file>
 	
 	<file id="file2" name_suffix="_grid_U" description="ocean U grid variables" >
-          <field_group group_ref="SPIN-UP_U" />	
+          <field field_ref="e3u"          name="e3u"        standard_name="cell_thickness" />
+          <field field_ref="uoce"         name="vozocrtx"     operation="average" freq_op="1mo" > @uoce_e3u / @e3u </field>
+          <field field_ref="uoce"         name="vozocrtx2"    operation="average" freq_op="1mo" > @uoce2_e3u / @e3u </field>
+          <field field_ref="ut"           name="ut"           operation="average" freq_op="1mo" > @ut_e3u / @e3u </field>
+          <field field_ref="us"           name="us"           operation="average" freq_op="1mo" > @us_e3u / @e3u </field>
+	  <field field_ref="utau"         name="sozotaux"   standard_name="surface_downward_x_stress" />
 	</file>
 	
 	<file id="file3" name_suffix="_grid_V" description="ocean V grid variables" >
-          <field_group group_ref="SPIN-UP_V" />
-	</file>
-
-	<file id="file4" name_suffix="_grid_W" description="ocean W grid variables" >
-          <field_group group_ref="SPIN-UP_W" />
+          <field field_ref="e3v"          name="e3v"        standard_name="cell_thickness" />
+          <field field_ref="voce"         name="vomecrty"     operation="average" freq_op="1mo" > @voce_e3v / @e3v </field>
+          <field field_ref="voce"         name="vomecrty2"    operation="average" freq_op="1mo" > @voce2_e3v / @e3v </field>
+          <field field_ref="vt"           name="vt"           operation="average" freq_op="1mo" > @vt_e3v / @e3v </field>
+          <field field_ref="vs"           name="vs"           operation="average" freq_op="1mo" > @vs_e3v / @e3v </field>
+	  <field field_ref="vtau"         name="sometauy"   standard_name="surface_downward_y_stress" />
 	</file>
 	
-	<file id="file5" name_suffix="_diaptr" description="process diagnostics" >
-          <field_group group_ref="SPIN-UP_diaptr" />
-	</file>	
-	
+	<file id="file4" name_suffix="_grid_W" description="ocean W grid variables" >
+	  <field field_ref="woce"         name="vovecrtz"   standard_name="upward_sea_water_velocity"                         />
+          <field field_ref="wt"           name="wt"         standard_name=""                                                  />
+          <field field_ref="ws"           name="ws"         standard_name=""                                                  />
+	  <field field_ref="avt"          name="votkeavt"   standard_name="ocean_vertical_heat_diffusivity"                   />
+          <field field_ref="avm"          name="votkeavm"   standard_name="ocean_vertical_momentum_diffusivity"               />
+          <field field_ref="aht2d"        name="soleahtw"   standard_name="vertical_component_of_ocean_xy_tracer_diffusivity" />
+	</file>
+        
+        <file id="file5" name_suffix="_diaptr" description="processed diagnostics" >
+	  <field field_ref="zomsfglo"     name="zomsfglo"   standard_name="meridional_streamfunction_global"           grid_ref="gznl_W_3D"       />
+          <field field_ref="zotemglo"     name="zotemglo"   standard_name="zonal_mean_temperature_global"              grid_ref="gznl_T_3D"       />
+          <field field_ref="zosrfglo"     name="zosrfglo"   standard_name="zonal_mean_surface_global"                  grid_ref="gznl_T_3D"       />
+	  <field field_ref="zosalglo"     name="zosalglo"   standard_name="zonal_mean_salinity_global"                 grid_ref="gznl_T_3D"       />          
+          <field field_ref="zomsfatl"     name="zomsfatl"   standard_name="meridional_streamfunction_atlantic"         grid_ref="gznl_W_3D"       />
+          <field field_ref="zotematl"     name="zotematl"   standard_name="zonal_mean_temperature_atlantic"            grid_ref="gznl_T_3D"       />
+          <field field_ref="zosrfatl"     name="zosrfatl"   standard_name="zonal_mean_surface_atlantic"                grid_ref="gznl_T_3D"       />
+	  <field field_ref="zosalatl"     name="zosalatl"   standard_name="zonal_mean_salinity_atlantic"               grid_ref="gznl_T_3D"       />          
+          <field field_ref="zomsfpac"     name="zomsfpac"   standard_name="meridional_streamfunction_pacific"          grid_ref="gznl_W_3D"       />
+          <field field_ref="zotempac"     name="zotempac"   standard_name="zonal_mean_temperature_pacific"             grid_ref="gznl_T_3D"       />
+          <field field_ref="zosrfpac"     name="zosrfpac"   standard_name="zonal_mean_surface_pacific"                 grid_ref="gznl_T_3D"       />
+	  <field field_ref="zosalpac"     name="zosalpac"   standard_name="zonal_mean_salinity_pacific"                grid_ref="gznl_T_3D"       />
+          <field field_ref="zomsfind"     name="zomsfind"   standard_name="meridional_streamfunction_indian"           grid_ref="gznl_W_3D"       />
+          <field field_ref="zotemind"     name="zotemind"   standard_name="zonal_mean_temperature_indian"              grid_ref="gznl_T_3D"       />
+          <field field_ref="zosrfind"     name="zosrfind"   standard_name="zonal_mean_surface_indian"                  grid_ref="gznl_T_3D"       />
+	  <field field_ref="zosalind"     name="zosalind"   standard_name="zonal_mean_salinity_indian"                 grid_ref="gznl_T_3D"       />          
+          <field field_ref="zomsfipc"     name="zomsfipc"   standard_name="meridional_streamfunction_indian_pacific"   grid_ref="gznl_W_3D"       />
+          <field field_ref="zotemipc"     name="zotemipc"   standard_name="zonal_mean_temperature_indian_pacific"      grid_ref="gznl_T_3D"       />
+          <field field_ref="zosrfipc"     name="zosrfipc"   standard_name="zonal_mean_surface_indian_pacific"          grid_ref="gznl_T_3D"       />
+	  <field field_ref="zosalipc"     name="zosalipc"   standard_name="zonal_mean_salinity_indian_pacific"         grid_ref="gznl_T_3D"       />          
+          <field field_ref="sophtadv"     name="sophtadv"   standard_name="advective_heat_transport"                   grid_ref="gznl_T_2D"       />
+          <field field_ref="sopstadv"     name="sopstadv"   standard_name="advective_salt_transport"                   grid_ref="gznl_T_2D"       />
+          <field field_ref="sophtldf"     name="sophtldf"   standard_name="diffusive_heat_transport"                   grid_ref="gznl_T_2D"       />
+          <field field_ref="sopstldf"     name="sopstldf"   standard_name="diffusive_salt_transport"                   grid_ref="gznl_T_2D"       />
+          <field field_ref="sophtadv_atl" name="sophtadv_atl"   standard_name="advective_heat_transport_atlantic"      grid_ref="gznl_T_2D"       />
+          <field field_ref="sopstadv_atl" name="sopstadv_atl"   standard_name="advective_salt_transport_atlantic"      grid_ref="gznl_T_2D"       />
+          <field field_ref="sophtldf_atl" name="sophtldf_atl"   standard_name="diffusive_heat_transport_atlantic"      grid_ref="gznl_T_2D"       />
+          <field field_ref="sopstldf_atl" name="sopstldf_atl"   standard_name="diffusive_salt_transport_atlantic"      grid_ref="gznl_T_2D"       />
+          <field field_ref="sophtadv_pac" name="sophtadv_pac"   standard_name="advective_heat_transport_pacific"       grid_ref="gznl_T_2D"       />
+          <field field_ref="sopstadv_pac" name="sopstadv_pac"   standard_name="advective_salt_transport_pacific"       grid_ref="gznl_T_2D"       />
+          <field field_ref="sophtldf_pac" name="sophtldf_pac"   standard_name="diffusive_heat_transport_pacific"       grid_ref="gznl_T_2D"       />
+          <field field_ref="sopstldf_pac" name="sopstldf_pac"   standard_name="diffusive_salt_transport_pacific"       grid_ref="gznl_T_2D"       />
+          <field field_ref="sophtadv_ind" name="sophtadv_ind"   standard_name="advective_heat_transport_indian"        grid_ref="gznl_T_2D"       />
+          <field field_ref="sopstadv_ind" name="sopstadv_ind"   standard_name="advective_salt_transport_indian"        grid_ref="gznl_T_2D"       />
+          <field field_ref="sophtldf_ind" name="sophtldf_ind"   standard_name="diffusive_heat_transport_indian"        grid_ref="gznl_T_2D"       />
+          <field field_ref="sopstldf_ind" name="sopstldf_ind"   standard_name="diffusive_salt_transport_indian"        grid_ref="gznl_T_2D"       />
+          <field field_ref="sophtadv_ipc" name="sophtadv_ipc"   standard_name="advective_heat_transport_indopacific"   grid_ref="gznl_T_2D"       />
+          <field field_ref="sopstadv_ipc" name="sopstadv_ipc"   standard_name="advective_salt_transport_indopacific"   grid_ref="gznl_T_2D"       />
+          <field field_ref="sophtldf_ipc" name="sophtldf_ipc"   standard_name="diffusive_heat_transport_indopacific"   grid_ref="gznl_T_2D"       />
+          <field field_ref="sopstldf_ipc" name="sopstldf_ipc"   standard_name="diffusive_salt_transport_indopacific"   grid_ref="gznl_T_2D"       />
+	</file>
+
         <file id="file6" name_suffix="_ptrc_T" description="Medusa sms variables" >
           <field_group group_ref="groupMEDUSA" />
@@ -89,27 +176,124 @@
       
   
-      <!-- file_group id="1m" output_freq="1mo" split_freq="1mo" output_level="10" enabled=".TRUE." -->  <!-- 1mo files -->   
-      <file_group id="1m" output_freq="1mo" output_level="10" enabled=".FALSE.">  <!-- 1mo files -->   
+      <file_group id="1m" output_freq="1mo" split_freq="1mo" output_level="10" enabled=".FALSE."> <!-- 1mo files -->   
+      <!-- <file_group id="1m" output_freq="1mo" output_level="10" enabled=".TRUE.">  - 1mo files -->   
 
         <file id="file9" name_suffix="_grid_T" description="ocean T grid variables" >
-          <field_group group_ref="SPIN-UP_T" />
-	</file>
-	
+          <field field_ref="e3t"          name="e3t"        standard_name="cell_thickness" />
+          <field field_ref="toce"         name="votemper"        operation="average" freq_op="1mo" > @toce_e3t / @e3t </field>
+          <field field_ref="toce"         name="votemper2"       operation="average" freq_op="1mo" > @toce2_e3t / @e3t </field>
+          <field field_ref="soce"         name="vosaline"        operation="average" freq_op="1mo" > @soce_e3t / @e3t </field>
+          <field field_ref="soce"         name="vosaline2"       operation="average" freq_op="1mo" > @soce2_e3t / @e3t </field>
+          <field field_ref="ttrd_totad"   name="opottempadvect"  operation="average" freq_op="1mo" > @ttrd_totad_e3t * 1026.0 * 3991.86795711963 </field>
+          <field field_ref="ttrd_iso"     name="opottemppmdiff"  operation="average" freq_op="1mo" > @ttrd_iso_e3t * 1026.0 * 3991.86795711963 </field>
+          <field field_ref="ttrd_zdfp"    name="opottempdiff"    operation="average" freq_op="1mo" > @ttrd_zdfp_e3t * 1026.0 * 3991.86795711963 </field>
+          <field field_ref="ttrd_evd"     name="ttrd_evd"        operation="average" freq_op="1mo" > @ttrd_evd_e3t / @e3t </field>
+          <field field_ref="ttrd_qns"     name="ttrd_qns"        operation="average" freq_op="1mo" > @ttrd_qns_e3t / @e3t_surf </field>
+          <field field_ref="ttrd_qsr"     name="rsdoabsorb"      operation="average" freq_op="1mo" > @ttrd_qsr_e3t * 1026.0 * 3991.86795711963 </field>
+          <field field_ref="ttrd_bbl"     name="ttrd_bbl"        operation="average" freq_op="1mo" > @ttrd_bbl_e3t / @e3t </field>
+          <field field_ref="ttrd_tot"     name="opottemptend"    operation="average" freq_op="1mo" > @ttrd_tot_e3t * 1026.0 * 3991.86795711963 </field>
+          <field field_ref="strd_totad"   name="osaltadvect"     operation="average" freq_op="1mo" > @strd_totad_e3t * 1026.0 * 0.001 </field>
+          <field field_ref="strd_iso"     name="osaltpmdiff"     operation="average" freq_op="1mo" > @strd_iso_e3t * 1026.0 * 0.001  </field>
+          <field field_ref="strd_zdfp"    name="osaltdiff"       operation="average" freq_op="1mo" > @strd_zdfp_e3t * 1026.0 * 0.001  </field>
+          <field field_ref="strd_evd"     name="strd_evd"        operation="average" freq_op="1mo" > @strd_evd_e3t / @e3t </field>
+          <field field_ref="strd_tot"     name="osalttend"       operation="average" freq_op="1mo" > @strd_tot_e3t * 1026.0 * 0.001 </field>
+          <field field_ref="strd_bbl"     name="strd_bbl"        operation="average" freq_op="1mo" > @strd_bbl_e3t / @e3t </field>
+          <field field_ref="ssh"          name="sossheig"   standard_name="sea_surface_height_above_geoid"                          />
+          <field field_ref="ssh2"         name="zossq"      standard_name="square_of_sea_surface_height_above_geoid"                />
+          <field field_ref="empmr"        name="sowaflup"   standard_name="water_flux_out_of_sea_ice_and_sea_water"                 />
+          <field field_ref="saltflx"      name="sosafldo"   standard_name="salt_flux_into_sea_water"                                />
+          <field field_ref="qt"           name="sohefldo"   standard_name="surface_downward_heat_flux_in_sea_water"                 />
+          <field field_ref="qsr"          name="soshfldo"   standard_name="downwelling_photosynthetic_radiative_flux_in_sea_water"               />
+          <field field_ref="mldkz5"       name="somixhgt"   standard_name="ocean_mixed_layer_thickness_defined_by_vertical_tracer_diffusivity"   />
+          <field field_ref="mldr10_1"     name="somxl010"   standard_name="ocean_mixed_layer_thickness_defined_by_sigma_theta"      />
+          <field field_ref="mldzint_1"    name="sokaraml"   standard_name="ocean_mixed_layer_thickness"                             />
+          <field field_ref="mldzint_2"    name="somlddbm"   standard_name="ocean_mixed_layer_thickness"                             />
+          <field field_ref="ice_cover"    name="soicecov"   standard_name="sea_ice_area_fraction"                                   />
+          <field field_ref="wspd"         name="sowindsp"   standard_name="wind_speed"                                              />
+          <field field_ref="qhcisf"       name="sohflisf"   standard_name="ice_shelf_heat_flux"                                                        />
+          <field field_ref="qlatisf"      name="soqlatisf"  standard_name="ice_shelf_latent_heat_flux"                                                        />
+          <field field_ref="fwfisf"       name="sowflisf"   standard_name="ice_shelf_melting"                                                        />
+          <field field_ref="berg_total_melt"       name="berg_total_melt"          standard_name=""                                 />
+          <field field_ref="berg_total_heat_flux"  name="berg_total_heat_flux"     standard_name=""                                 />
+          <field field_ref="runoffs"       name="sorunoff"    standard_name="water_flux_into_sea_water_from_rivers"                 />
+          <field field_ref="hflx_rain_cea" name="hfrainds"    standard_name="temperature_flux_due_to_rainfall_expressed_as_heat_flux_into_sea_water"/>
+          <field field_ref="hflx_evap_cea" name="hfevapds"    standard_name="temperature_flux_due_to_evaporation_expressed_as_heat_flux_out_of_sea_water"/>
+          <field field_ref="evap_ao_cea"   name="evap_ao_cea" standard_name="water_evaporation_flux"                                />
+          <field field_ref="rain"          name="pr"          standard_name="rainfall_flux"                                         />
+          <field field_ref="snowpre"       name="snowpre"     standard_name="snowfall_flux"                                         />
+          <field field_ref="snow_ao_cea"   name="snow_ao_cea" standard_name="snowfall_flux_over_ocean"                              />
+          <field field_ref="snow_ai_cea"   name="snow_ai_cea" standard_name="snowfall_flux_over_sea_ice"                            />
+        </file>
+
 	<file id="file10" name_suffix="_grid_U" description="ocean U grid variables" >
-          <field_group group_ref="SPIN-UP_U" />
-	</file>
-	
+          <field field_ref="e3u"          name="e3u"        standard_name="cell_thickness" />
+          <field field_ref="uoce"         name="vozocrtx"     operation="average" freq_op="1mo" > @uoce_e3u / @e3u </field>
+          <field field_ref="uoce"         name="vozocrtx2"    operation="average" freq_op="1mo" > @uoce2_e3u / @e3u </field>
+          <field field_ref="ut"           name="ut"           operation="average" freq_op="1mo" > @ut_e3u / @e3u </field>
+          <field field_ref="us"           name="us"           operation="average" freq_op="1mo" > @us_e3u / @e3u </field>
+          <field field_ref="utau"         name="sozotaux"   standard_name="surface_downward_x_stress" />
+        </file>
+
 	<file id="file11" name_suffix="_grid_V" description="ocean V grid variables" >
-          <field_group group_ref="SPIN-UP_V" />
-	</file>
-	
+          <field field_ref="e3v"          name="e3v"        standard_name="cell_thickness" />
+          <field field_ref="voce"         name="vomecrty"     operation="average" freq_op="1mo" > @voce_e3v / @e3v </field>
+          <field field_ref="voce"         name="vomecrty2"    operation="average" freq_op="1mo" > @voce2_e3v / @e3v </field>
+          <field field_ref="vt"           name="vt"           operation="average" freq_op="1mo" > @vt_e3v / @e3v </field>
+          <field field_ref="vs"           name="vs"           operation="average" freq_op="1mo" > @vs_e3v / @e3v </field>
+          <field field_ref="vtau"         name="sometauy"   standard_name="surface_downward_y_stress" />
+        </file>
+
 	<file id="file12" name_suffix="_grid_W" description="ocean W grid variables" >
-          <field_group group_ref="SPIN-UP_W" />
-	</file>
-	
+          <field field_ref="woce"         name="vovecrtz"   standard_name="upward_sea_water_velocity"                         />
+          <field field_ref="wt"           name="wt"         standard_name=""                                                  />
+          <field field_ref="ws"           name="ws"         standard_name=""                                                  />
+          <field field_ref="avt"          name="votkeavt"   standard_name="ocean_vertical_heat_diffusivity"                   />
+          <field field_ref="avm"          name="votkeavm"   standard_name="ocean_vertical_momentum_diffusivity"               />
+          <field field_ref="aht2d"        name="soleahtw"   standard_name="vertical_component_of_ocean_xy_tracer_diffusivity" />
+        </file>
+
 	<file id="file13" name_suffix="_grid_W" description="ocean W grid variables" >
-          <field_group group_ref="SPIN-UP_diaptr" />
-	</file>
-	
+          <field field_ref="zomsfglo"     name="zomsfglo"   standard_name="meridional_streamfunction_global"           grid_ref="gznl_W_3D"       />
+          <field field_ref="zotemglo"     name="zotemglo"   standard_name="zonal_mean_temperature_global"              grid_ref="gznl_T_3D"       />
+          <field field_ref="zosrfglo"     name="zosrfglo"   standard_name="zonal_mean_surface_global"                  grid_ref="gznl_T_3D"       />
+          <field field_ref="zosalglo"     name="zosalglo"   standard_name="zonal_mean_salinity_global"                 grid_ref="gznl_T_3D"       />
+          <field field_ref="zomsfatl"     name="zomsfatl"   standard_name="meridional_streamfunction_atlantic"         grid_ref="gznl_W_3D"       />
+          <field field_ref="zotematl"     name="zotematl"   standard_name="zonal_mean_temperature_atlantic"            grid_ref="gznl_T_3D"       />
+          <field field_ref="zosrfatl"     name="zosrfatl"   standard_name="zonal_mean_surface_atlantic"                grid_ref="gznl_T_3D"       />
+          <field field_ref="zosalatl"     name="zosalatl"   standard_name="zonal_mean_salinity_atlantic"               grid_ref="gznl_T_3D"       />
+          <field field_ref="zomsfpac"     name="zomsfpac"   standard_name="meridional_streamfunction_pacific"          grid_ref="gznl_W_3D"       />
+          <field field_ref="zotempac"     name="zotempac"   standard_name="zonal_mean_temperature_pacific"             grid_ref="gznl_T_3D"       />
+          <field field_ref="zosrfpac"     name="zosrfpac"   standard_name="zonal_mean_surface_pacific"                 grid_ref="gznl_T_3D"       />
+          <field field_ref="zosalpac"     name="zosalpac"   standard_name="zonal_mean_salinity_pacific"                grid_ref="gznl_T_3D"       />
+          <field field_ref="zomsfind"     name="zomsfind"   standard_name="meridional_streamfunction_indian"           grid_ref="gznl_W_3D"       />
+          <field field_ref="zotemind"     name="zotemind"   standard_name="zonal_mean_temperature_indian"              grid_ref="gznl_T_3D"       />
+          <field field_ref="zosrfind"     name="zosrfind"   standard_name="zonal_mean_surface_indian"                  grid_ref="gznl_T_3D"       />
+          <field field_ref="zosalind"     name="zosalind"   standard_name="zonal_mean_salinity_indian"                 grid_ref="gznl_T_3D"       />
+          <field field_ref="zomsfipc"     name="zomsfipc"   standard_name="meridional_streamfunction_indian_pacific"   grid_ref="gznl_W_3D"       />
+          <field field_ref="zotemipc"     name="zotemipc"   standard_name="zonal_mean_temperature_indian_pacific"      grid_ref="gznl_T_3D"       />
+          <field field_ref="zosrfipc"     name="zosrfipc"   standard_name="zonal_mean_surface_indian_pacific"          grid_ref="gznl_T_3D"       />
+          <field field_ref="zosalipc"     name="zosalipc"   standard_name="zonal_mean_salinity_indian_pacific"         grid_ref="gznl_T_3D"       />
+          <field field_ref="sophtadv"     name="sophtadv"   standard_name="advective_heat_transport"                   grid_ref="gznl_T_2D"       />
+          <field field_ref="sopstadv"     name="sopstadv"   standard_name="advective_salt_transport"                   grid_ref="gznl_T_2D"       />
+          <field field_ref="sophtldf"     name="sophtldf"   standard_name="diffusive_heat_transport"                   grid_ref="gznl_T_2D"       />
+          <field field_ref="sopstldf"     name="sopstldf"   standard_name="diffusive_salt_transport"                   grid_ref="gznl_T_2D"       />
+          <field field_ref="sophtadv_atl" name="sophtadv_atl"   standard_name="advective_heat_transport_atlantic"      grid_ref="gznl_T_2D"       />
+          <field field_ref="sopstadv_atl" name="sopstadv_atl"   standard_name="advective_salt_transport_atlantic"      grid_ref="gznl_T_2D"       />
+          <field field_ref="sophtldf_atl" name="sophtldf_atl"   standard_name="diffusive_heat_transport_atlantic"      grid_ref="gznl_T_2D"       />
+          <field field_ref="sopstldf_atl" name="sopstldf_atl"   standard_name="diffusive_salt_transport_atlantic"      grid_ref="gznl_T_2D"       />
+          <field field_ref="sophtadv_pac" name="sophtadv_pac"   standard_name="advective_heat_transport_pacific"       grid_ref="gznl_T_2D"       />
+          <field field_ref="sopstadv_pac" name="sopstadv_pac"   standard_name="advective_salt_transport_pacific"       grid_ref="gznl_T_2D"       />
+          <field field_ref="sophtldf_pac" name="sophtldf_pac"   standard_name="diffusive_heat_transport_pacific"       grid_ref="gznl_T_2D"       />
+          <field field_ref="sopstldf_pac" name="sopstldf_pac"   standard_name="diffusive_salt_transport_pacific"       grid_ref="gznl_T_2D"       />
+          <field field_ref="sophtadv_ind" name="sophtadv_ind"   standard_name="advective_heat_transport_indian"        grid_ref="gznl_T_2D"       />
+          <field field_ref="sopstadv_ind" name="sopstadv_ind"   standard_name="advective_salt_transport_indian"        grid_ref="gznl_T_2D"       />
+          <field field_ref="sophtldf_ind" name="sophtldf_ind"   standard_name="diffusive_heat_transport_indian"        grid_ref="gznl_T_2D"       />
+          <field field_ref="sopstldf_ind" name="sopstldf_ind"   standard_name="diffusive_salt_transport_indian"        grid_ref="gznl_T_2D"       />
+          <field field_ref="sophtadv_ipc" name="sophtadv_ipc"   standard_name="advective_heat_transport_indopacific"   grid_ref="gznl_T_2D"       />
+          <field field_ref="sopstadv_ipc" name="sopstadv_ipc"   standard_name="advective_salt_transport_indopacific"   grid_ref="gznl_T_2D"       />
+          <field field_ref="sophtldf_ipc" name="sophtldf_ipc"   standard_name="diffusive_heat_transport_indopacific"   grid_ref="gznl_T_2D"       />
+          <field field_ref="sopstldf_ipc" name="sopstldf_ipc"   standard_name="diffusive_salt_transport_indopacific"   grid_ref="gznl_T_2D"       />
+        </file>
+
         <file id="file14" name_suffix="_ptrc_T" description="Medusa sms variables" >
           <field_group group_ref="groupMEDUSA" />
@@ -147,5 +331,5 @@
 ============================================================================================================
     -->
-
+    
    <axis_definition>
       <axis id="deptht" standard_name="depth" long_name="Vertical T levels" unit="m" positive="down" />
@@ -155,25 +339,95 @@
       <axis id="nfloat" long_name="Float number"      unit="-"  />
       <axis id="icbcla" long_name="Iceberg class"     unit="-"  />
+      <axis id="ncatice" long_name="Ice category"       unit="1"               />
+      <axis id="iax_20C" long_name="20 degC isotherm"   unit="degC"            />
+      <axis id="iax_28C" long_name="28 degC isotherm"   unit="degC"            />
+      <axis id="deptht_surface" axis_ref="deptht" >
+         <zoom_axis begin=" 1 " n=" 1 " />
+      </axis>
+
    </axis_definition>
     
    <domain_definition src="./domain_def.xml"/>
    
-   <grid_definition>    
-     <grid id="grid_T_2D" domain_ref="grid_T"/>
-     <grid id="grid_T_3D" domain_ref="grid_T" axis_ref="deptht"/>
-     <grid id="grid_U_2D" domain_ref="grid_U"/>
-     <grid id="grid_U_3D" domain_ref="grid_U" axis_ref="depthu"/>
-     <grid id="grid_V_2D" domain_ref="grid_V"/>
-     <grid id="grid_V_3D" domain_ref="grid_V" axis_ref="depthv"/>
-     <grid id="grid_W_2D" domain_ref="grid_W"/>
-     <grid id="grid_W_3D" domain_ref="grid_W" axis_ref="depthw"/>
-     <grid id="grid_T_2D_zoom" domain_ref="zoomT"/>
-     <grid id="grid_T_3D_zoom" domain_ref="zoomT" axis_ref="deptht"/>
-     <grid id="grid_U_3D_zoom" domain_ref="zoomU" axis_ref="depthu"/>
-     <grid id="grid_V_3D_zoom" domain_ref="zoomV" axis_ref="depthv"/>
-     <grid id="gznl_T_2D" domain_ref="gznl"/>
-     <grid id="gznl_T_3D" domain_ref="gznl" axis_ref="deptht"/>
-     <grid id="gznl_W_3D" domain_ref="gznl" axis_ref="depthw"/>
-    </grid_definition>   
+    <grid_definition>
+        <!--  -->
+       <grid id="grid_T_2D" >
+         <domain id="grid_T" />
+       </grid>
+        <!--  -->
+       <grid id="grid_T_3D_ncatice" >
+         <domain id="grid_T" />
+         <axis id="ncatice" />
+       </grid>
+        <!--  -->
+       <grid id="grid_T_3D" >
+         <domain id="grid_T" />
+         <axis id="deptht" />
+       </grid>
+        <!--  -->
+       <grid id="grid_U_2D" >
+         <domain id="grid_U" />
+       </grid>
+        <!--  -->
+       <grid id="grid_U_3D" >
+         <domain id="grid_U" />
+         <axis id="depthu" />
+       </grid>
+        <!--  -->
+       <grid id="grid_V_2D" >
+         <domain id="grid_V" />
+       </grid>
+        <!--  -->
+       <grid id="grid_V_3D" >
+         <domain id="grid_V" />
+         <axis id="depthv" />
+       </grid>
+        <!--  -->
+       <grid id="grid_W_2D" >
+         <domain id="grid_W" />
+       </grid>
+        <!--  -->
+       <grid id="grid_W_3D" >
+         <domain id="grid_W" />
+         <axis id="depthw" />
+       </grid>
+        <!--  -->
+       <grid id="grid_1point" >
+         <domain id="1point"/>
+       </grid>
+        <!--  -->
+       <grid id="grid_T_nfloat" >
+         <domain id="grid_T" />
+         <axis id="nfloat" />
+       </grid>
+        <!--  -->
+       <grid id="grid_EqT" >
+         <domain id="EqT" />
+       </grid>
+        <!--  -->
+       <grid id="gznl_T_2D">
+         <domain id="ptr" />
+       </grid>
+        <!--  -->
+       <grid id="gznl_T_3D">
+         <domain id="ptr" />
+         <axis id="deptht" />
+       </grid>
+        <!--  -->
+       <grid id="gznl_W_2D">
+         <domain id="ptr" />
+       </grid>
+        <!--  -->
+       <grid id="gznl_W_3D">
+         <domain id="ptr" />
+         <axis id="depthw" />
+       </grid>
+        <!--  -->
+       <grid id="grid_T_surface_extract">
+         <domain domain_ref="grid_T" />
+         <axis   axis_ref="deptht_surface" />
+       </grid>
+
+    </grid_definition>
   </context>
   
@@ -183,13 +437,8 @@
       <variable_definition>
 	
-     <!-- 
-        We must have buffer_size > jpi*jpj*jpk*8 (with jpi and jpj the subdomain size)
--->
-	  <variable id="buffer_size"               type="integer">30000000</variable>
-	  <variable id="buffer_server_factor_size" type="integer">2</variable>
-	  <variable id="info_level"                type="integer">0</variable>
-	  <variable id="using_server"              type="boolean">true</variable>
-	  <variable id="using_oasis"               type="boolean">true</variable>
-	  <variable id="oasis_codes_id"            type="string" >toyoce</variable>
+	  <variable id="info_level"            type="int">0</variable>
+	  <variable id="using_server"          type="bool">false</variable>
+	  <variable id="using_oasis"           type="bool">true</variable>
+	  <variable id="oasis_codes_id"        type="string" >toyoce</variable>
 	
       </variable_definition>
Index: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90
===================================================================
--- branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90	(revision 7720)
+++ branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90	(revision 7766)
@@ -445,5 +445,5 @@
       !!----------------------------------------------------------------------
       USE lib_mpp , ONLY: ctl_warn
-      INTEGER ::   ierr(7)        ! Local variables
+      INTEGER ::   ierr(8)        ! Local variables
       !!----------------------------------------------------------------------
       ierr(:) = 0
Index: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
===================================================================
--- branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 7720)
+++ branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 7766)
@@ -533,4 +533,5 @@
       foxy_cons(:,:) = 0.0      !! oxygen consumption
       foxy_anox(:,:) = 0.0      !! unrealised oxygen consumption
+      !!
 # endif
       ftot_n(:,:)   = 0.0       !! N inventory 
@@ -571,4 +572,13 @@
       f_riv_c(:,:)   = 0.0      !! riverine C   input 
       f_riv_alk(:,:) = 0.0      !! riverine alk input 
+      !! 
+      !! Jpalm -- 06-03-2017 -- Forgotten var to init
+      f_omarg(:,:) = 0.0        !!
+      f_omcal(:,:) = 0.0 
+      xFree(:,:) = 0.0          !! state variables for iron-ligand system
+      fcomm_resp(:,:) = 0.0 
+      fprn_ml(:,:) = 0.0        !! mixed layer PP diagnostics
+      fprd_ml(:,:) = 0.0        !! mixed layer PP diagnostics
+
       !!
       !! allocate and initiate 2D diag
@@ -724,5 +734,5 @@
          ENDIF     
 !!
-!! skip FDT_XXXX, RG_XXXXF, FDS_XXXX, RGSXXXXF diagnostics here
+!! skip FDT_XXXX, RG_XXXXF, FDS_XXXX, RGS_XXXXF diagnostics here
 !!
          IF( med_diag%REMINN%dgsave ) THEN
Index: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcdms_medusa.F90
===================================================================
--- branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcdms_medusa.F90	(revision 7720)
+++ branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcdms_medusa.F90	(revision 7766)
@@ -5,5 +5,5 @@
    !!======================================================================
    !! History :
-   !!  -   !  2014-08  (J. Palmiéri - A. Yool)    added for UKESM1 project
+   !!  -   !  2014-08  (J. Palmieri - A. Yool)    added for UKESM1 project
    !!----------------------------------------------------------------------
 #if defined key_medusa && defined key_roam
@@ -44,5 +44,5 @@
       !!
       !! Title  : Calculates DMS ocean surface concentration
-      !! Author : Julien Palmiéri and Andrew Yool
+      !! Author : Julien Palmieri and Andrew Yool
       !! Date   : 08/08/14 
       !!
@@ -118,5 +118,5 @@
 !
 ! AXY (13/03/15): Simo & Dachs (2002)
-        fq1 = (-1 * log(mld)) + 5.7
+        fq1 = (-1.0 * log(mld)) + 5.7
         fq2 = (55.8 * cmr) + 0.6
         if (cmr < 0.02) then
@@ -136,5 +136,5 @@
 !        
 ! AXY (13/03/15): Halloran et al. (2010)
-        fq1 = (-1 * log(mld)) + 5.7
+        fq1 = (-1.0 * log(mld)) + 5.7
         fq2 = (55.8 * cmr) + 0.6
         fq3 = (90.0 / mld)
@@ -179,2 +179,3 @@
 
 
+
Index: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcdms_medusa.F90_COLIN4
===================================================================
--- branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcdms_medusa.F90_COLIN4	(revision 7720)
+++ 	(revision )
@@ -1,219 +1,0 @@
-MODULE trcdms_medusa
-   !!======================================================================
-   !!                         ***  MODULE trcdms_medusa  ***
-   !! TOP :   MEDUSA
-   !!======================================================================
-   !! History :
-   !!  -   !  2014-08  (J. Palmiéri - A. Yool)    added for UKESM1 project
-   !!----------------------------------------------------------------------
-#if defined key_medusa && defined key_roam
-   !!----------------------------------------------------------------------
-   !!                                        MEDUSA DMS surface concentration
-   !!----------------------------------------------------------------------
-   !!   trc_dms_medusa        :  
-   !!----------------------------------------------------------------------
-      USE oce_trc
-      USE trc
-      USE sms_medusa
-      USE lbclnk
-      USE prtctl_trc      ! Print control for debugging
-      USE in_out_manager  ! I/O manager
-
-      IMPLICIT NONE
-      PRIVATE
-
-      PUBLIC   trc_dms_medusa    ! called in trc_bio_medusa
-
-   !!* Substitution
-#  include "domzgr_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
-   !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
-
-CONTAINS
-
-!=======================================================================
-!
-   SUBROUTINE trc_dms_medusa( chn, chd, mld, xqsr, xdin,   &  !! inputs
-     &  dms_andr, dms_simo, dms_aran, dms_hall )              !! outputs
-!      
-!=======================================================================
-      !!
-      !! Title  : Calculates DMS ocean surface concentration
-      !! Author : Julien Palmiéri and Andrew Yool
-      !! Date   : 08/08/14 
-      !!
-      !! DMS module is called in trc_bio's huge jk,jj,ji loop
-      !! --> DMS concentration is calculated in a specific cell 
-      !! (no need of ji,jj,jk)
-      !!
-      !! AXY (13/03/15): amend to include all four schemes tested
-      !!                 during winter/spring 2015; these are:
-      !!
-      !!                 1. Anderson et al. (2001); this uses fields
-      !!                    of surface chl, irradiance and nutrients
-      !!                    to empirically estimate DMS via a broken
-      !!                    stick approach
-      !!
-      !!                 2. Simo & Dachs (2002); this uses fields of
-      !!                    surface chl and mixed layer depth
-      !!
-      !!                 3. Aranami & Tsunogai (2004); this is an
-      !!                    embellishment of Simo & Dachs
-      !!
-      !!                 4. Halloran et al. (2010); this is an
-      !!                    alternative embellishment of Sim & Dachs
-      !!                    and is included because it is formally
-      !!                    published (and different from the above)
-      !!
-      !! AXY (08/07/15): amend to remove Julien's original calculation
-      !!                 as this is now superfluous; the four schemes 
-      !!                 are calculated and one is chosen to be passed
-      !!                 to the atmosphere in trc_bio_medusa
-      !!
-!=======================================================================
-
-      IMPLICIT NONE
-!
-      REAL(wp), INTENT( in )    :: chn                  !! non-diatom chlorophyll    (mg/m3)
-      REAL(wp), INTENT( in )    :: chd                  !! diatom chlorophyll        (mg/m3)
-      REAL(wp), INTENT( in )    :: mld                  !! mix layer depth           (m)
-      REAL(wp), INTENT( in )    :: xqsr                 !! surface irradiance        (W/m2)
-      REAL(wp), INTENT( in )    :: xdin                 !! surface DIN               (mmol N/m3)
-      REAL(wp), INTENT( inout ) :: dms_andr             !! DMS surface concentration (nmol/L) 
-      REAL(wp), INTENT( inout ) :: dms_simo             !! DMS surface concentration (nmol/L) 
-      REAL(wp), INTENT( inout ) :: dms_aran             !! DMS surface concentration (nmol/L) 
-      REAL(wp), INTENT( inout ) :: dms_hall             !! DMS surface concentration (nmol/L) 
-!
-      REAL(wp) :: CHL, cmr, sw_dms
-      REAL(wp) :: Jterm, Qterm
-      !! temporary variables
-      REAL(wp) :: fq1,fq2,fq3
-! AXY (22/02/17): Colin Jones hack v4
-      REAL(wp) :: fracmin, mldeep1, mldeep2, mlda
-      REAL(wp) :: mldcntl1, mldcntl2, mldcntl3, mlduse
-! 
-!=======================================================================
-!
-! AXY (13/03/15): per remarks above, the following calculations estimate
-!                 DMS using all of the schemes examined for UKESM1
-!
-      CHL = 0.0
-      CHL = chn+chd                                 !! mg/m3 
-      cmr = CHL / mld
-!
-! AXY (13/03/15): Anderson et al. (2001)
-        Jterm = xqsr + 1.0e-6
-        !! this next line makes a hard-coded assumption about the 
-        !! half-saturation constant of MEDUSA (which should be
-        !! done properly; perhaps even scaled with the proportion
-        !! of diatoms and non-diatoms)
-        Qterm = xdin / (xdin + 0.5)
-        fq1 = log10(CHL * Jterm * Qterm)
-        if (fq1 > 1.72) then
-           dms_andr = (8.24 * (fq1 - 1.72)) + 2.29
-        else
-           dms_andr = 2.29
-        endif
-!
-! AXY (13/03/15): Simo & Dachs (2002)
-        fq1 = (-1.0 * log(mld)) + 5.7
-        fq2 = (55.8 * cmr) + 0.6
-        if (cmr < 0.02) then
-           dms_simo = fq1
-        else
-           dms_simo = fq2
-        endif
-!           
-! AXY (13/03/15): Aranami & Tsunogai (2004)
-        fq1 = 60.0 / mld
-        fq2 = (55.8 * cmr) + 0.6
-        if (cmr < 0.02) then
-           dms_aran = fq1
-        else
-           dms_aran = fq2
-        endif
-!        
-! AXY (13/03/15): Halloran et al. (2010)
-        fq1 = (-1.0 * log(mld)) + 5.7
-        fq2 = (55.8 * cmr) + 0.6
-        fq3 = (90.0 / mld)
-        if (cmr < 0.02) then
-           dms_hall = fq1
-        else
-           dms_hall = fq2
-        endif
-        if (mld > 182.5) then
-           dms_hall = fq3
-        endif
-!
-! AXY (22/02/17): Aranami & Tsunogai (2004); Colin Jones hack v4
-!
-! This version modifies Aranami & Tsunogai (2004) to increase DMS
-! especially in the Southern Ocean; to avoid adding extra outputs,
-! it currently hijacks and overwrites Simo & Dachs (2002)
-!
-        fracmin  = 0.1   ! minimum value by which MLD is multiplied 
-        mldeep1  = 12.0  ! MLD depth where modification begins to be applied
-        mldeep2  = 100.0 ! MLD depth where we maximise the effect to fracmin
-        mlda     = 3.0   ! power term of 1 - Q**mlda
-        mldcntl1 = 1.0   !
-        mldcntl2 = 1.0   ! initialised; just in case ...
-        mldcntl3 = 1.0   !
-!
-        if (mld .ge. mldeep1) then
-           ! mldcntl1 lies between 0.15 and 1.0 depending on the value of Q;
-           ! as Q increases towards 1, mldcntl1 tends to a minimum of 0.1
-           mldcntl1 = min(1.0, max(fracmin, (1.0 - Qterm**mlda)))
-           ! a second control lies between = 5 (mld = mldeep1) and = 1
-           ! (mld = mldeep2); this is used to ensure mldcntl1 primarily
-           ! operates on deep MLD with high Q values
-           mldcntl2 = min(5.0, max(1.0, (mldeep2 / mld)))
-           mldcntl3 = min(mldcntl1 * mldcntl2, 1.0)
-           mlduse   = mld * mldcntl3
-        else
-           mlduse   = mld
-        endif
-        cmr = CHL / mlduse
-        fq1 = 60.0 / mlduse
-        fq2 = (55.8 * cmr) + 0.6
-        if (cmr < 0.02) then
-           dms_simo = fq1
-        else
-           dms_simo = fq2
-        endif
-
-  END SUBROUTINE trc_dms_medusa
-
-
-!=======================================================================
-!=======================================================================
-!=======================================================================
-
-#else
-   !!======================================================================
-   !!  Dummy module :                                   No MEDUSA bio-model
-   !!======================================================================
-
-CONTAINS
-
-!=======================================================================
-!
-   SUBROUTINE trc_dms_medusa( kt )                                        !! EMPTY Routine
-!      
-!
-      INTEGER, INTENT( in ) ::   kt
-!
-
-      WRITE(*,*) 'trc_dms_medusa: You should not have seen this print! error?'
-
-   END SUBROUTINE trc_dms_medusa
-#endif
-
-   !!======================================================================
-END MODULE  trcdms_medusa
-
-
Index: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90
===================================================================
--- branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 7720)
+++ branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 7766)
@@ -84,13 +84,9 @@
       &    xburial,                                           &
       &  jriver_n,jriver_si,jriver_c,jriver_alk,jriver_dep,   &
-      &    friver_dep,                                        &
       &  xsdiss,                                              &
-      &  vsed,xhr,                                            &
       &  sedlam,sedlostpoc,jpkb,jdms,jdms_input,jdms_model
 #if defined key_roam
       NAMELIST/natroam/ xthetaphy,xthetazoo,xthetanit,        &
-      &    xthetarem,xo2min,                                  &
-      &    f3_pH,f3_h2co3,f3_hco3,f3_co3,f3_omcal,f3_omarg,   &
-      &    f2_ccd_cal,f2_ccd_arg
+      &    xthetarem,xo2min 
 #endif
       NAMELIST/natopt/xkg0,xkr0,xkgp,xkrp,xlg,xlr,rpig
Index: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90_save
===================================================================
--- branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90_save	(revision 7720)
+++ 	(revision )
@@ -1,1096 +1,0 @@
-MODULE trcnam_medusa
-   !!======================================================================
-   !!                      ***  MODULE trcnam_medusa  ***
-   !! TOP :   initialisation of some run parameters for MEDUSA bio-model
-   !!======================================================================
-   !! History :   2.0  !  2007-12  (C. Ethe, G. Madec) Original code
-   !!              -   !  2008-08  (K. Popova) adaptation for MEDUSA
-   !!              -   !  2008-11  (A. Yool) continuing adaptation for MEDUSA
-   !!              -   !  2010-03  (A. Yool) updated for branch inclusion
-   !!              -   !  2011-04  (A. Yool) updated for ROAM project
-   !!              -   !  2013-05  (A. Yool) renamed (from trclsm) for v3.5
-   !!----------------------------------------------------------------------
-#if defined key_medusa
-   !!----------------------------------------------------------------------
-   !!   'key_medusa'   :                                       MEDUSA model
-   !!----------------------------------------------------------------------
-   !! trc_nam_medusa      : MEDUSA model initialisation
-   !!----------------------------------------------------------------------
-   USE oce_trc         ! Ocean variables
-   USE par_trc         ! TOP parameters
-   USE trc             ! TOP variables
-   USE sms_medusa      ! sms trends
-   USE iom             ! I/O manager
-
-   !! AXY (04/02/14): necessary to find NaNs on HECTOR
-   USE, INTRINSIC :: ieee_arithmetic 
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   trc_nam_medusa   ! called by trcnam.F90 module
-
-   !!* Substitution
-#  include "domzgr_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
-   !! $Id: trcnam_medusa.F90 1162 2008-07-02 09:33:33Z cetlod $ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
-CONTAINS
-
-   SUBROUTINE trc_nam_medusa
-      !!----------------------------------------------------------------------
-      !!                     ***  trc_nam_medusa  ***  
-      !!
-      !! ** Purpose :   read MEDUSA namelist
-      !!
-      !! ** input   :   file 'namelist.trc.sms' containing the following
-      !!             namelist: natbio, natopt, and natdbi ("key_trc_diabio")
-      !!
-      !! ekp: namelist nabio contains ALL parameters of the ecosystem 
-      !!      point sourses and sinks PLUS sediment exchange
-      !!      dia_bio - used by Lobster to output all point terms 
-      !!                (sourses and sinks of bio)
-      !!      dia_add - additional diagnostics for biology such as 
-      !!                primary production (2d depth integrated field or 3d)
-      !!----------------------------------------------------------------------
-      !!
-      INTEGER            :: ji,jj,jk
-      REAL(wp)           :: fthk, fdep, fdep1
-      REAL(wp)           :: q1, q2, q3
-      !
-      NAMELIST/natbio/ xxi,xaln,xald,jphy,xvpn,xvpd,          &
-      &    xsin0,xnsi0,xuif,jliebig,                          &
-      &    xthetam,xthetamd,xnln,xnld,xsld,xfln,xfld,         &
-      &  xgmi,xgme,xkmi,xkme,xphi,xbetan,xbetac,xkc,          &
-      &    xpmipn,xpmid,xpmepn,xpmepd,xpmezmi,xpmed,          &
-      &  xmetapn,xmetapd,xmetazmi,xmetazme,                   &
-      &  jmpn,xmpn,xkphn,jmpd,xmpd,xkphd,jmzmi,xmzmi,xkzmi,   &
-      &    jmzme,xmzme,xkzme,jmd,jsfd,xmd,xmdc,               &
-      &  xthetapn,xthetapd,xthetazmi,xthetazme,xthetad,       &
-      &    xrfn,xrsn,vsed,xhr,                                &
-      &  jiron,xfe_mass,xfe_sol,xfe_sed,xLgT,xk_FeL,xk_sc_Fe, &
-      &  jexport,jfdfate,jrratio,jocalccd,xridg_r0,           &
-      &    xfdfrac1,xfdfrac2,xfdfrac3,                        &
-      &    xcaco3a,xcaco3b,xmassc,xmassca,xmasssi,xprotca,    &
-      &    xprotsi,xfastc,xfastca,xfastsi,                    &
-      &  jorgben,jinorgben,xsedn,xsedfe,xsedsi,xsedc,xsedca,  &
-      &    xburial,                                           &
-      &  jriver_n,jriver_si,jriver_c,jriver_alk,jriver_dep,   &
-      &    friver_dep,                                        &
-      &  xsdiss,                                              &
-      &  vsed,xhr,                                            &
-      &  sedlam,sedlostpoc,jpkb,jdms
-#if defined key_roam
-      NAMELIST/natroam/ xthetaphy,xthetazoo,xthetanit,        &
-      &    xthetarem,xo2min,                                  &
-      &    f3_pH,f3_h2co3,f3_hco3,f3_co3,f3_omcal,f3_omarg,   &
-      &    f2_ccd_cal,f2_ccd_arg
-#endif
-      NAMELIST/natopt/xkg0,xkr0,xkgp,xkrp,xlg,xlr,rpig
-      INTEGER :: jl, jn
-      TYPE(DIAG), DIMENSION(jp_medusa_2d)  :: meddia2d
-      TYPE(DIAG), DIMENSION(jp_medusa_3d)  :: meddia3d
-      TYPE(DIAG), DIMENSION(jp_medusa_trd) :: meddiabio
-      CHARACTER(LEN=32)   ::   clname
-      !!
-      NAMELIST/nammeddia/ meddia3d, meddia2d     ! additional diagnostics
-
-      !!----------------------------------------------------------------------
-
-      IF(lwp) WRITE(numout,*)
-      clname = 'namelist_medusa'
-      IF(lwp) WRITE(numout,*) ' trc_nam_medusa: read MEDUSA namelist'
-      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
-# if defined key_debug_medusa
-      CALL flush(numout)
-# endif
-
-
-      CALL ctl_opn( numnatm, TRIM( clname ), 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
-
-# if defined key_debug_medusa
-      CALL flush(numout)
-      IF (lwp) write (numout,*) '------------------------------'
-      IF (lwp) write (numout,*) 'Jpalm - debug'
-      IF (lwp) write (numout,*) 'open namelist_medusa -- OK'
-      IF (lwp) write (numout,*) 'Now, read namilists inside :'
-      IF (lwp) write (numout,*) ' '
-# endif
-      !
-!     IF( .NOT.lk_iomput .AND. ln_diatrc ) THEN
-! should be done via iom_put exclusively but there are a lot of them; this is the easy first attempt !
-         !
-         ! Namelist nammeddia
-         ! -------------------
-         DO jl = 1, jp_medusa_2d
-            WRITE(meddia2d(jl)%sname,'("2D_",I1)') jl                      ! short name
-            WRITE(meddia2d(jl)%lname,'("2D DIAGNOSTIC NUMBER ",I2)') jl    ! long name
-            meddia2d(jl)%units = ' '                                       ! units
-         END DO
-# if defined key_debug_medusa
-      CALL flush(numout)
-# endif
-         !                                 ! 3D output arrays
-         DO jl = 1, jp_medusa_3d
-            WRITE(meddia3d(jl)%sname,'("3D_",I1)') jl                      ! short name
-            WRITE(meddia3d(jl)%lname,'("3D DIAGNOSTIC NUMBER ",I2)') jl    ! long name
-            meddia3d(jl)%units = ' '                                       ! units
-         END DO
-# if defined key_debug_medusa
-      IF (lwp) write (numout,*) '------------------------------'
-      IF (lwp) write (numout,*) 'Jpalm - debug'
-      IF (lwp) write (numout,*) 'Just before reading namelist_medusa :: nammeddia'
-      IF (lwp) write (numout,*) ' '
-      CALL flush(numout)
-# endif
-
-         REWIND( numnatm )               !
-         READ  ( numnatm, nammeddia )
-
-# if defined key_debug_medusa
-      IF (lwp) write (numout,*) '------------------------------'
-      IF (lwp) write (numout,*) 'Jpalm - debug'
-      IF (lwp) write (numout,*) 'reading namelist_medusa :: nammeddia OK'
-      IF (lwp) write (numout,*) 'Check number of variable in nammeddia:'
-      IF (lwp) write (numout,*) 'jp_medusa_2d: ',jp_medusa_2d ,'jp_medusa_3d: ',jp_medusa_3d
-      IF (lwp) write (numout,*) ' '
-      CALL flush(numout)
-# endif
-         DO jl = 1, jp_medusa_2d
-            jn = jp_msa0_2d + jl - 1
-            ctrc2d(jn) = meddia2d(jl)%sname
-            ctrc2l(jn) = meddia2d(jl)%lname
-            ctrc2u(jn) = meddia2d(jl)%units
-         END DO
-
-         DO jl = 1, jp_medusa_3d
-            jn = jp_msa0_3d + jl - 1
-            ctrc3d(jn) = meddia3d(jl)%sname
-            ctrc3l(jn) = meddia3d(jl)%lname
-            ctrc3u(jn) = meddia3d(jl)%units
-         END DO
-
-         IF(lwp) THEN                   ! control print
-# if defined key_debug_medusa
-      IF (lwp) write (numout,*) '------------------------------'
-      IF (lwp) write (numout,*) 'Jpalm - debug'
-      IF (lwp) write (numout,*) 'Var name assignation OK'
-      IF (lwp) write (numout,*) 'next check var names'
-      IF (lwp) write (numout,*) ' '
-            CALL flush(numout)
-# endif
-            WRITE(numout,*)
-            WRITE(numout,*) ' Namelist : natadd'
-            DO jl = 1, jp_medusa_3d
-               jn = jp_msa0_3d + jl - 1
-               WRITE(numout,*) '  3d diag nb : ', jn, '    short name : ', ctrc3d(jn), &
-                 &             '  long name  : ', ctrc3l(jn), '   unit : ', ctrc3u(jn)
-            END DO
-            WRITE(numout,*) ' '
-
-            DO jl = 1, jp_medusa_2d
-               jn = jp_msa0_2d + jl - 1
-               WRITE(numout,*) '  2d diag nb : ', jn, '    short name : ', ctrc2d(jn), &
-                 &             '  long name  : ', ctrc2l(jn), '   unit : ', ctrc2u(jn)
-            END DO
-            WRITE(numout,*) ' '
-         ENDIF
-         !
-# if defined key_debug_medusa
-            CALL flush(numout)
-# endif
-
-!     ENDIF
-
-      ! 1.4 namelist natbio : biological parameters
-      ! -------------------------------------------
-      
-      xxi         = 0.
-      xaln        = 0.
-      xald        = 0.
-      jphy        = 0
-      xvpn        = 0.
-      xvpd        = 0.
-      xthetam     = 0.
-      xthetamd    = 0.
-!!
-      xsin0       = 0.
-      xnsi0       = 0.
-      xuif        = 0.
-!!
-      jliebig     = 0
-      xnln        = 0.
-      xnld        = 0.
-      xsld        = 0.
-      xfln        = 0.
-      xfld        = 0.
-!!
-      xgmi        = 0.
-      xgme        = 0.
-      xkmi        = 0.
-      xkme        = 0.
-      xphi	  = 0.
-      xbetan      = 0.
-      xbetac      = 0.
-      xkc         = 0.
-      xpmipn      = 0.
-      xpmid       = 0.
-      xpmepn      = 0.
-      xpmepd      = 0.
-      xpmezmi     = 0.
-      xpmed       = 0.
-!!
-      xmetapn     = 0.
-      xmetapd     = 0.
-      xmetazmi    = 0.
-      xmetazme    = 0.
-!!
-      jmpn        = 0
-      xmpn        = 0.
-      xkphn       = 0.
-      jmpd        = 0
-      xmpd        = 0.
-      xkphd       = 0.
-      jmzmi       = 0
-      xmzmi       = 0.
-      xkzmi       = 0.
-      jmzme       = 0
-      xmzme       = 0.
-      xkzme       = 0.
-!!
-      jmd         = 0
-      jsfd        = 0
-      xmd         = 0.
-      xmdc        = 0.
-!!
-      xthetapn    = 0.
-      xthetapd    = 0.
-      xthetazmi   = 0.
-      xthetazme   = 0.
-      xthetad     = 0.
-      xrfn        = 0.
-      xrsn        = 0.  !: (NOT USED HERE; RETAINED FOR LOBSTER)
-!!
-      jiron       = 0
-      xfe_mass    = 0.
-      xfe_sol     = 0.
-      xfe_sed     = 0.
-      xLgT        = 0.
-      xk_FeL	  = 0.
-      xk_sc_Fe    = 0.
-!!
-      jexport     = 0
-      jfdfate     = 0
-      jrratio     = 0
-      jocalccd    = 0
-      xridg_r0    = 0.
-      xfdfrac1	  = 0.
-      xfdfrac2	  = 0.
-      xfdfrac3	  = 0.
-      xcaco3a	  = 0.
-      xcaco3b	  = 0.
-      xmassc	  = 0.
-      xmassca	  = 0.
-      xmasssi	  = 0.
-      xprotca	  = 0.
-      xprotsi	  = 0.
-      xfastc	  = 0.
-      xfastca	  = 0.
-      xfastsi	  = 0.
-!!
-      jorgben     = 0
-      jinorgben   = 0
-      xsedn       = 0.
-      xsedfe      = 0.
-      xsedsi      = 0.
-      xsedc       = 0.
-      xsedca      = 0.
-      xburial     = 0.
-!!
-      jriver_n    = 0
-      jriver_si   = 0
-      jriver_c    = 0
-      jriver_alk  = 0
-      jriver_dep  = 1
-!!
-      xsdiss	  = 0.
-!!
-      vsed        = 0.
-      xhr         = 0.
-!!
-      sedlam	  = 0.
-      sedlostpoc  = 0.
-      jpkb	  = 0.
-      jdms        = 0
-            
-      REWIND(numnatm)
-      READ(numnatm,natbio)
-
-!! Primary production and chl related quantities
-!!       xxi         :  conversion factor from gC to mmolN 
-!!       xaln        :  Chl-a specific initial slope of P-I curve for non-diatoms 
-!!       xald        :  Chl-a specific initial slope of P-I curve for diatoms
-!!       jphy        :  phytoplankton T-dependent growth switch
-!!       xvpn        :  maximum growth rate for non-diatoms
-!!       xvpd        :  maximum growth rate for diatoms
-!!       xthetam     :  maximum Chl to C ratio for non-diatoms      
-!!       xthetamd    :  maximum Chl to C ratio for diatoms      
-!!
-!! Diatom silicon parameters
-!!       xsin0       :  minimum diatom Si:N ratio
-!!       xnsi0       :  minimum diatom N:Si ratio
-!!       xuif        :  hypothetical growth ratio at infinite Si:N ratio
-!!
-!! Nutrient limitation
-!!       jliebig     :  Liebig nutrient uptake switch
-!!       xnln        :  half-sat constant for DIN uptake by non-diatoms 
-!!       xnld        :  half-sat constant for DIN uptake by diatoms 
-!!       xsl         :  half-sat constant for Si uptake by diatoms 
-!!       xfld        :  half-sat constant for Fe uptake by diatoms  
-!!       xfln        :  half-sat constant for Fe uptake by non-datoms 
-!!
-!! Grazing
-!!       xgmi        :  microzoo maximum growth rate 
-!!       xgme        :  mesozoo maximum growth rate 
-!!       xkmi        :  microzoo grazing half-sat parameter
-!!       xkme        :  mesozoo grazing half-sat parameter
-!!       xphi        :  micro/mesozoo grazing inefficiency
-!!       xbetan      :  micro/mesozoo N assimilation efficiency
-!!       xbetac      :  micro/mesozoo C assimilation efficiency
-!!       xkc         :  micro/mesozoo net C growth efficiency
-!!       xpmipn      :  grazing preference of microzoo for non-diatoms
-!!       xpmid       :  grazing preference of microzoo for diatoms
-!!       xpmepn      :  grazing preference of mesozoo for non-diatoms 
-!!       xpmepd      :  grazing preference of mesozoo for diatoms
-!!       xpmezmi     :  grazing preference of mesozoo for microzoo
-!!       xpmed       :  grazing preference of mesozoo for detritus
-!!
-!! Metabolic losses
-!!       xmetapn     :  non-diatom metabolic loss rate
-!!       xmetapd     :  diatom     metabolic loss rate
-!!       xmetazmi    :  microzoo   metabolic loss rate
-!!       xmetazme    :  mesozoo    metabolic loss rate
-!!
-!! Mortality/Remineralisation
-!!       jmpn        :  non-diatom mortality functional form
-!!       xmpn        :  non-diatom mortality rate
-!!       xkphn       :  non-diatom mortality half-sat constant
-!!       jmpd        :  diatom     mortality functional form
-!!       xmpd        :  diatom mortality rate
-!!       xkphd       :  diatom mortality half-sat constant
-!!       jmzmi       :  microzoo   mortality functional form
-!!       xmzmi       :  microzoo mortality rate
-!!       xkzmi       :  microzoo mortality half-sat constant
-!!       jmzme       :  mesozoo    mortality functional form
-!!       xmzme       :  mesozoo mortality rate
-!!       xkzme       :  mesozoo mortality half-sat constant
-!!
-!! Remineralisation
-!!       jmd         :  detritus T-dependent remineralisation switch
-!!       jsfd        :  accelerate seafloor detritus remin. switch
-!!       xmd         :  detrital nitrogen remineralisation rate 
-!!       xmdc        :  detrital carbon remineralisation rate 
-!!
-!! Stochiometric ratios
-!!       xthetapn    :  non-diatom C:N ratio
-!!       xthetapd    :  diatom C:N ratio
-!!       xthetazmi   :  microzoo C:N ratio
-!!       xthetazme   :  mesozoo C:N ratio
-!!       xthetad     :  detritus C:N ratio
-!!       xrfn        :  phytoplankton Fe:N ratio
-!!	 xrsn        :  diatom Si:N ratio (*NOT* used)
-!!
-!! Iron parameters
-!!       jiron       :  iron scavenging submodel switch
-!!       xfe_mass    :  iron atomic mass
-!!	 xfe_sol     :  aeolian iron solubility
-!!	 xfe_sed     :  sediment iron input
-!!	 xLgT	     :  total ligand concentration (umol/m3)
-!!	 xk_FeL	     :  dissociation constant for (Fe + L)
-!!	 xk_sc_Fe    :  scavenging rate of "free" iron
-!!	 
-!! Fast-sinking detritus parameters
-!!       jexport     :  fast detritus remineralisation switch
-!!       jfdfate     :  fate of fast detritus at seafloor switch
-!!       jrratio     :  rain ratio switch
-!!       jocalccd    :  CCD switch
-!!       xridg_r0    :  Ridgwell rain ratio coefficient
-!!       xfdfrac1    :  fast-sinking fraction of diatom nat. mort. losses
-!!       xfdfrac2    :  fast-sinking fraction of meszooplankton mort. losses
-!!       xfdfrac3    :  fast-sinking fraction of diatom silicon grazing losses
-!!       xcaco3a     :  polar (high latitude) CaCO3 fraction
-!!       xcaco3b     :  equatorial (low latitude) CaCO3 fraction
-!!       xmassc      :  organic C mass:mole ratio, C106 H175 O40 N16 P1
-!!       xmassca     :  calcium carbonate mass:mole ratio, CaCO3
-!!       xmasssi     :  biogenic silicon mass:mole ratio, (H2SiO3)n
-!!       xprotca     :  calcium carbonate protection ratio
-!!       xprotsi     :  biogenic silicon protection ratio
-!!       xfastc      :  organic C remineralisation length scale
-!!       xfastca     :  calcium carbonate dissolution length scale
-!!       xfastsi     :  biogenic silicon dissolution length scale
-!!
-!! Benthic 
-!!       jorgben     :  does   organic detritus go to the benthos?
-!!       jinorgben   :  does inorganic detritus go to the benthos?
-!!       xsedn       :  organic   nitrogen sediment remineralisation rate 
-!!       xsedfe      :  organic   iron     sediment remineralisation rate 
-!!       xsedsi      :  inorganic silicon  sediment dissolution      rate 
-!!       xsedc       :  organic   carbon   sediment remineralisation rate 
-!!       xsedca      :  inorganic carbon   sediment dissolution      rate 
-!!       xburial     :  burial rate of seafloor detritus
-!!
-!! Riverine inputs
-!!       jriver_n    :  riverine N          input?
-!!       jriver_si   :  riverine Si         input?
-!!       jriver_c    :  riverine C          input?
-!!       jriver_alk  :  riverine alkalinity input?
-!!       jriver_dep  :  depth of riverine   input?
-!!
-!! Miscellaneous
-!!       xsdiss      :  diatom frustule dissolution rate
-!!
-!! Gravitational sinking      
-!!       vsed        :  detritus gravitational sinking rate 
-!!       xhr         :  coeff for Martin's remineralisation profile
-!!
-!! Additional parameters
-!!       sedlam      :  time coeff of POC in sediments
-!!      sedlostpoc   :  sediment geol loss for POC
-!!       jpkb        :  vertical layer for diagnostic of the vertical flux 
-!!                      NOTE that in LOBSTER it is a first vertical layers where 
-!!                      biology is active  
-!!
-!! UKESM1 - new diagnostics  !! Jpalm
-!!       jdms        :  include dms diagnostics
-!!
-!!
-!!
-
-      IF(lwp) THEN
-!!
-!! AXY (08/11/13): compilation key notification
-         WRITE(numout,*) '=== Compilation keys'
-#if defined key_roam
-         WRITE(numout,*)     &
-         &   ' key_roam                                                               = ACTIVE'
-#else
-         WRITE(numout,*)     &
-         &   ' key_roam                                                               = INACTIVE'
-#endif        
-#if defined key_axy_carbchem
-         WRITE(numout,*)     &
-         &   ' key_axy_carbchem                                                       = ACTIVE'
-#else
-         WRITE(numout,*)     &
-         &   ' key_axy_carbchem                                                       = INACTIVE'
-#endif        
-#if defined key_bs_axy_zforce
-         WRITE(numout,*)     &
-         &   ' key_bs_axy_zforce                                                      = ACTIVE'
-#else
-         WRITE(numout,*)     &
-         &   ' key_bs_axy_zforce                                                      = INACTIVE'
-#endif        
-#if defined key_bs_axy_yrlen
-         WRITE(numout,*)     &
-         &   ' key_bs_axy_yrlen                                                       = ACTIVE'
-#else
-         WRITE(numout,*)     &
-         &   ' key_bs_axy_yrlen                                                       = INACTIVE'
-#endif        
-#if defined key_deep_fe_fix
-         WRITE(numout,*)     &
-         &   ' key_deep_fe_fix                                                        = ACTIVE'
-#else
-         WRITE(numout,*)     &
-         &   ' key_deep_fe_fix                                                        = INACTIVE'
-#endif        
-#if defined key_axy_nancheck
-         WRITE(numout,*)     &
-         &   ' key_axy_nancheck                                                       = ACTIVE'
-#else
-         WRITE(numout,*)     &
-         &   ' key_axy_nancheck                                                       = INACTIVE'
-#endif        
-# if defined key_axy_pi_co2
-         WRITE(numout,*)     &
-         &   ' key_axy_pi_co2                                                         = ACTIVE'
-#else
-         WRITE(numout,*)     &
-         &   ' key_axy_pi_co2                                                         = INACTIVE'
-# endif
-         WRITE(numout,*) ' '
-
-         WRITE(numout,*) 'natbio'
-         WRITE(numout,*) ' '
-!!
-!! Primary production and chl related quantities
-         WRITE(numout,*) '=== Primary production'
-         WRITE(numout,*)     &
-         &   ' conversion factor from gC to mmolN,                        xxi         =', xxi
-         WRITE(numout,*)     &
-         &   ' Chl-a specific initial slope of P-I curve for non-diatoms, xaln        = ', xaln
-         WRITE(numout,*)     &
-         &   ' Chl-a specific initial slope of P-I curve for diatoms,     xald        = ', xald
-         if (jphy.eq.1) then
-            WRITE(numout,*) &
-            &   ' phytoplankton growth is *temperature-dependent*            jphy        = ', jphy
-         elseif (jphy.eq.0) then
-            WRITE(numout,*) &
-            &   ' phytoplankton growth is *temperature-independent*          jphy        = ', jphy
-         endif
-         WRITE(numout,*)     &
-         &   ' maximum growth rate for non-diatoms,                       xvpn        = ', xvpn
-         WRITE(numout,*)     &
-         &   ' maximum growth rate for diatoms,                           xvpn        = ', xvpd
-         WRITE(numout,*)     &
-         &   ' maximum Chl to C ratio for non-diatoms,                    xthetam     = ', xthetam
-         WRITE(numout,*)     &
-         &   ' maximum Chl to C ratio for diatoms,                        xthetamd    = ', xthetamd
-!!
-!! Diatom silicon parameters
-         WRITE(numout,*) '=== Diatom silicon parameters'
-         WRITE(numout,*)     &
-         &   ' minimum diatom Si:N ratio,                                 xsin0       = ', xsin0
-         WRITE(numout,*)     &
-         &   ' minimum diatom N:Si ratio,                                 xnsi0       = ', xnsi0
-         WRITE(numout,*)     &
-         &   ' hypothetical growth ratio at infinite Si:N ratio,          xuif        = ', xuif
-!!
-!! Nutrient limitation
-         WRITE(numout,*) '=== Nutrient limitation'
-         if (jliebig.eq.1) then
-            WRITE(numout,*) &
-            &   ' nutrient uptake is a Liebig Law (= most limiting) function jliebig     = ', jliebig
-         elseif (jliebig.eq.0) then
-            WRITE(numout,*) &
-            &   ' nutrient uptake is a multiplicative function               jliebig     = ', jliebig
-         endif
-         WRITE(numout,*)     &
-         &   ' half-sat constant for DIN uptake by non-diatoms,           xnln        = ', xnln
-         WRITE(numout,*)     &
-         &   ' half-sat constant for DIN uptake by diatoms,               xnld        = ', xnld
-         WRITE(numout,*)     &
-         &   ' half-sat constant for Si uptake by diatoms,                xsld        = ', xsld
-         WRITE(numout,*)     &
-         &   ' half-sat constant for Fe uptake by non-diatoms,            xfln        = ', xfln
-         WRITE(numout,*)     &
-         &   ' half-sat constant for Fe uptake by diatoms,                xfld        = ', xfld
-!!
-!! Grazing
-         WRITE(numout,*) '=== Zooplankton grazing'
-         WRITE(numout,*)     &
-         &   ' microzoo maximum growth rate,                              xgmi        = ', xgmi
-         WRITE(numout,*)     &
-         &   ' mesozoo maximum growth rate,                               xgme        = ', xgme
-         WRITE(numout,*)     &
-         &   ' microzoo grazing half-sat parameter,                       xkmi        = ', xkmi
-         WRITE(numout,*)     &
-         &   ' mesozoo grazing half-sat parameter,                        xkme        = ', xkme
-         WRITE(numout,*)     &
-         &   ' micro/mesozoo grazing inefficiency,                        xphi        = ', xphi
-         WRITE(numout,*)     &
-         &   ' micro/mesozoo N assimilation efficiency,                   xbetan      = ', xbetan
-         WRITE(numout,*)     &
-         &   ' micro/mesozoo C assimilation efficiency,                   xbetac      = ', xbetan
-         WRITE(numout,*)     &
-         &   ' micro/mesozoo net C growth efficiency,                     xkc         = ', xkc
-         WRITE(numout,*)     &
-         &   ' grazing preference of microzoo for non-diatoms,            xpmipn      = ', xpmipn
-         WRITE(numout,*)     &
-         &   ' grazing preference of microzoo for detritus,               xpmid       = ', xpmid
-         WRITE(numout,*)     &
-         &   ' grazing preference of mesozoo for non-diatoms,             xpmepn      = ', xpmepn
-         WRITE(numout,*)     &
-         &   ' grazing preference of mesozoo for diatoms,                 xpmepd      = ', xpmepd
-         WRITE(numout,*)     &
-         &   ' grazing preference of mesozoo for microzoo,                xpmezmi     = ', xpmezmi
-         WRITE(numout,*)     &
-         &   ' grazing preference of mesozoo for detritus,                xpmed       = ', xpmed
-!!
-!! Metabolic losses
-         WRITE(numout,*) '=== Metabolic losses'
-         WRITE(numout,*)     &
-         &   ' non-diatom metabolic loss rate,                            xmetapn     = ', xmetapn
-         WRITE(numout,*)     &
-         &   ' diatom     metabolic loss rate,                            xmetapd     = ', xmetapd
-         WRITE(numout,*)     &
-         &   ' microzoo   metabolic loss rate,                            xmetazmi    = ', xmetazmi
-         WRITE(numout,*)     &
-         &   ' mesozoo    metabolic loss rate,                            xmetazme    = ', xmetazme
-!!
-!! Mortality losses
-         WRITE(numout,*) '=== Mortality losses'
-         if (jmpn.eq.1) then
-            WRITE(numout,*)     &
-            &   ' non-diatom mortality functional form,            LINEAR    jmpn        = ', jmpn
-         elseif (jmpn.eq.2) then
-            WRITE(numout,*)     &
-            &   ' non-diatom mortality functional form,         QUADRATIC    jmpn        = ', jmpn
-         elseif (jmpn.eq.3) then
-            WRITE(numout,*)     &
-            &   ' non-diatom mortality functional form,        HYPERBOLIC    jmpn        = ', jmpn
-         elseif (jmpn.eq.4) then
-            WRITE(numout,*)     &
-            &   ' non-diatom mortality functional form,           SIGMOID    jmpn        = ', jmpn
-         endif
-         WRITE(numout,*)     &
-         &   ' non-diatom mortality rate,                                 xmpn        = ', xmpn
-         WRITE(numout,*)     &
-         &   ' non-diatom mortality half-sat constant                     xkphn       = ', xkphn
-         if (jmpd.eq.1) then
-            WRITE(numout,*)     &
-            &   ' diatom mortality functional form,                LINEAR    jmpd        = ', jmpd
-         elseif (jmpd.eq.2) then
-            WRITE(numout,*)     &
-            &   ' diatom mortality functional form,             QUADRATIC    jmpd        = ', jmpd
-         elseif (jmpd.eq.3) then
-            WRITE(numout,*)     &
-            &   ' diatom mortality functional form,            HYPERBOLIC    jmpd        = ', jmpd
-         elseif (jmpd.eq.4) then
-            WRITE(numout,*)     &
-            &   ' diatom mortality functional form,               SIGMOID    jmpd        = ', jmpd
-         endif
-         WRITE(numout,*)     &
-         &   ' diatom mortality rate,                                     xmpd        = ', xmpd
-         WRITE(numout,*)     &
-         &   ' diatom mortality half-sat constant                         xkphd       = ', xkphd
-         if (jmzmi.eq.1) then
-            WRITE(numout,*)     &
-            &   ' microzoo mortality functional form,              LINEAR    jmzmi       = ', jmzmi
-         elseif (jmzmi.eq.2) then
-            WRITE(numout,*)     &
-            &   ' microzoo mortality functional form,           QUADRATIC    jmzmi       = ', jmzmi
-         elseif (jmzmi.eq.3) then
-            WRITE(numout,*)     &
-            &   ' microzoo mortality functional form,          HYPERBOLIC    jmzmi       = ', jmzmi
-         elseif (jmzmi.eq.4) then
-            WRITE(numout,*)     &
-            &   ' microzoo mortality functional form,             SIGMOID    jmzmi       = ', jmzmi
-         endif
-         WRITE(numout,*)     &
-         &   ' microzoo mortality rate,                                   xmzmi       = ', xmzmi
-         WRITE(numout,*)     &
-         &   ' mesozoo mortality half-sat constant,                       xkzmi       = ', xkzmi
-         if (jmzme.eq.1) then
-            WRITE(numout,*)     &
-            &   ' mesozoo mortality functional form,               LINEAR    jmzme       = ', jmzme
-         elseif (jmzme.eq.2) then
-            WRITE(numout,*)     &
-            &   ' mesozoo mortality functional form,            QUADRATIC    jmzme       = ', jmzme
-         elseif (jmzme.eq.3) then
-            WRITE(numout,*)     &
-            &   ' mesozoo mortality functional form,           HYPERBOLIC    jmzme       = ', jmzme
-         elseif (jmzme.eq.4) then
-            WRITE(numout,*)     &
-            &   ' mesozoo mortality functional form,              SIGMOID    jmzme       = ', jmzme
-         endif
-         WRITE(numout,*)     &
-         &   ' mesozoo mortality rate,                                    xmzme       = ', xmzme
-         WRITE(numout,*)     &
-         &   ' mesozoo mortality half-sat constant,                       xkzme       = ', xkzme
-!!
-!! Remineralisation
-         WRITE(numout,*) '=== Remineralisation'
-         if (jmd.eq.1) then
-            WRITE(numout,*) &
-            &   ' detritus remineralisation is *temperature-dependent*       jmd         = ', jmd
-         elseif (jmd.eq.0) then
-            WRITE(numout,*) &
-            &   ' detritus remineralisation is *temperature-independent*     jmd         = ', jmd
-         endif
-         if (jsfd.eq.1) then
-            WRITE(numout,*) &
-            &   ' detritus seafloor remineralisation is *accelerated*        jsfd        = ', jsfd
-         else
-            WRITE(numout,*) &
-            &   ' detritus seafloor remineralisation occurs at same rate     jsfd        = ', jsfd
-         endif
-         WRITE(numout,*)     &
-         &   ' detrital nitrogen remineralisation rate,                   xmd         = ', xmd
-         WRITE(numout,*)     &
-         &   ' detrital carbon remineralisation rate,                     xmdc        = ', xmdc
-!!
-!! Stochiometric ratios
-         WRITE(numout,*) '=== Stoichiometric ratios'
-         WRITE(numout,*)     &
-         &   ' non-diatom C:N ratio,                                      xthetapn    = ', xthetapn
-         WRITE(numout,*)     &
-         &   ' diatom C:N ratio,                                          xthetapd    = ', xthetapd
-         WRITE(numout,*)     &
-         &   ' microzoo C:N ratio,                                        xthetazmi   = ', xthetazmi
-         WRITE(numout,*)     &
-         &   ' mesozoo C:N ratio,                                         xthetazme   = ', xthetazme
-         WRITE(numout,*)     &
-         &   ' detritus C:N ratio,                                        xthetad     = ', xthetad
-         WRITE(numout,*)     &
-         &   ' phytoplankton Fe:N ratio,                                  xrfn        = ', xrfn
-         WRITE(numout,*)     &
-         &   ' diatom Si:N ratio,                                         xrsn        = ', xrsn
-!!	  
-!! Iron parameters
-         WRITE(numout,*) '=== Iron parameters'
-         if (jiron.eq.1) then
-            WRITE(numout,*)     &
-            &   ' Dutkiewicz et al. (2005) iron scavenging                   jiron       = ', jiron
-         elseif (jiron.eq.2) then
-            WRITE(numout,*)     &
-            &   ' Moore et al. (2004) iron scavenging                        jiron       = ', jiron
-         elseif (jiron.eq.3) then
-            WRITE(numout,*)     &
-            &   ' Moore et al. (2008) iron scavenging                        jiron       = ', jiron
-         elseif (jiron.eq.4) then
-            WRITE(numout,*)     &
-            &   ' Galbraith et al. (2010) iron scavenging                    jiron       = ', jiron
-         else
-            WRITE(numout,*)     &
-            &   ' There is **no** iron scavenging                            jiron       = ', jiron
-         endif
-         WRITE(numout,*)     &
-         &   ' iron atomic mass,                                          xfe_mass    = ', xfe_mass
-         WRITE(numout,*)     &
-         &   ' aeolian iron solubility,                                   xfe_sol     = ', xfe_sol
-         WRITE(numout,*)     &
-         &   ' sediment iron input,                                       xfe_sed     = ', xfe_sed
-         WRITE(numout,*)     &
-         &   ' total ligand concentration (umol/m3),                      xLgT        = ', xLgT
-         WRITE(numout,*)     &
-         &   ' dissociation constant for (Fe + L),                        xk_FeL      = ', xk_FeL
-         WRITE(numout,*)     &
-         &   ' scavenging rate for free iron,                             xk_sc_Fe    = ', xk_sc_Fe
-!!
-!! Fast-sinking detritus parameters
-         WRITE(numout,*) '=== Fast-sinking detritus'
-         if (jexport.eq.1) then
-            WRITE(numout,*) &
-            &   ' fast-detritus remin. uses Dunne et al. (2007; ballast)     jexport     = ', jexport
-         elseif (jexport.eq.2) then
-            WRITE(numout,*) &
-            &   ' fast-detritus remin. uses Martin et al. (1987)             jexport     = ', jexport
-         elseif (jexport.eq.2) then
-            WRITE(numout,*) &
-            &   ' fast-detritus remin. uses Henson et al. (2011)             jexport     = ', jexport
-         endif
-         if (jfdfate.eq.1) then
-            WRITE(numout,*) &
-            &   ' fast-detritus reaching seafloor becomes slow-detritus      jfdfate     = ', jfdfate
-         elseif (jfdfate.eq.0) then
-            WRITE(numout,*) &
-            &   ' fast-detritus reaching seafloor instantly remineralised    jfdfate     = ', jfdfate
-         endif
-#if defined key_roam
-         if (jrratio.eq.0) then
-            WRITE(numout,*) &
-            &   ' Dunne et al. (2005) rain ratio submodel                    jrratio     = ', jrratio
-         elseif (jrratio.eq.1) then
-            WRITE(numout,*) &
-            &   ' Ridgwell et al. (2007) rain ratio submodel (surface omega) jrratio     = ', jrratio
-         elseif (jrratio.eq.2) then
-            WRITE(numout,*) &
-            &   ' Ridgwell et al. (2007) rain ratio submodel (3D omega)      jrratio     = ', jrratio
-         endif
-#else          
-         jrratio = 0
-         WRITE(numout,*) &
-         &   ' Dunne et al. (2005) rain ratio submodel                    jrratio     = ', jrratio
-#endif          
-#if defined key_roam
-         if (jocalccd.eq.0) then
-            WRITE(numout,*) &
-            &   ' Default, fixed CCD used                                    jocalccd    = ', jocalccd
-         elseif (jocalccd.eq.1) then
-            WRITE(numout,*) &
-            &   ' Calculated, dynamic CCD used                               jocalccd    = ', jocalccd
-         endif
-#else          
-         jocalccd = 0
-         WRITE(numout,*) &
-         &   ' Default, fixed CCD used                                    jocalccd    = ', jocalccd
-#endif
-         WRITE(numout,*)     &
-         &   ' Ridgwell rain ratio coefficient,                           xridg_r0    = ', xridg_r0
-         WRITE(numout,*)     &
-         &   ' fast-sinking fraction of diatom nat. mort. losses,         xfdfrac1    = ', xfdfrac1
-         WRITE(numout,*)     &
-         &   ' fast-sinking fraction of mesozooplankton mort. losses,     xfdfrac2    = ', xfdfrac2
-         WRITE(numout,*)     &
-         &   ' fast-sinking fraction of diatom silicon grazing losses,    xfdfrac3    = ', xfdfrac3
-         WRITE(numout,*)     &
-         &   ' polar (high latitude) CaCO3 fraction,                      xcaco3a     = ', xcaco3a
-         WRITE(numout,*)     &
-         &   ' equatorial (low latitude) CaCO3 fraction,                  xcaco3b     = ', xcaco3b
-         WRITE(numout,*)     &
-         &   ' organic C mass:mole ratio, C106 H175 O40 N16 P1,           xmassc      = ', xmassc
-         WRITE(numout,*)     &
-         &   ' calcium carbonate mass:mole ratio, CaCO3,                  xmassca     = ', xmassca
-         WRITE(numout,*)     &
-         &   ' biogenic silicon mass:mole ratio, (H2SiO3)n,               xmasssi     = ', xmasssi
-         WRITE(numout,*)     &
-         &   ' calcium carbonate protection ratio,                        xprotca     = ', xprotca
-         WRITE(numout,*)     &
-         &   ' biogenic silicon protection ratio,                         xprotsi     = ', xprotsi
-         WRITE(numout,*)     &
-         &   ' organic C remineralisation length scale,                   xfastc      = ', xfastc
-         WRITE(numout,*)     &
-         &   ' calcium carbonate dissolution length scale,                xfastca     = ', xfastca
-         WRITE(numout,*)     &
-         &   ' biogenic silicon dissolution length scale,                 xfastsi     = ', xfastsi
-!!
-!! Benthos parameters
-         WRITE(numout,*) '=== Benthos parameters'
-         WRITE(numout,*)     &
-         &   ' does   organic detritus go to the benthos?,                jorgben     = ', jorgben
-         WRITE(numout,*)     &
-         &   ' does inorganic detritus go to the benthos?,                jinorgben   = ', jinorgben
-!!
-!! Some checks on parameters related to benthos parameters
-         if (jorgben.eq.1 .and. jsfd.eq.1) then
-            !! slow detritus going to benthos at an accelerated rate
-            WRITE(numout,*) '  === WARNING! ==='
-            WRITE(numout,*) '  jsfd *and* jorgben are active - please check that you wish this'
-            WRITE(numout,*) '  === WARNING! ==='
-         endif
-         if (jorgben.eq.1 .and. jfdfate.eq.1) then
-            !! fast detritus going to benthos but via slow detritus
-            WRITE(numout,*) '  === WARNING! ==='
-            WRITE(numout,*) '  jfdfate *and* jorgben are active - please check that you wish this'
-            WRITE(numout,*) '  === WARNING! ==='
-         endif
-         if (jorgben.eq.0 .and. jinorgben.eq.1) then
-            !! inorganic fast detritus going to benthos but organic fast detritus is not
-            WRITE(numout,*) '  === WARNING! ==='
-            WRITE(numout,*) '  jinorgben is active but jorgben is not - please check that you wish this'
-            WRITE(numout,*) '  === WARNING! ==='
-         endif
-         WRITE(numout,*)     &
-         &   ' organic   nitrogen sediment remineralisation rate,         xsedn       = ', xsedn
-         WRITE(numout,*)     &
-         &   ' organic   iron     sediment remineralisation rate,         xsedfe      = ', xsedfe
-         WRITE(numout,*)     &
-         &   ' inorganic silicon  sediment remineralisation rate,         xsedsi      = ', xsedsi
-         WRITE(numout,*)     &
-         &   ' organic   carbon   sediment remineralisation rate,         xsedc       = ', xsedc
-         WRITE(numout,*)     &
-         &   ' inorganic carbon   sediment remineralisation rate,         xsedca      = ', xsedca
-         WRITE(numout,*)     &
-         &   ' burial rate of seafloor detritus,                          xburial     = ', xburial
-!!
-!! Riverine inputs
-         WRITE(numout,*) '=== Riverine inputs'
-         if (jriver_n.eq.0) then
-            WRITE(numout,*)     &
-            &   ' *no* riverine N input,                                     jriver_n    = ', jriver_n
-         elseif (jriver_n.eq.1) then
-            WRITE(numout,*)     &
-            &   ' riverine N concentrations specified,                       jriver_n    = ', jriver_n
-         elseif (jriver_n.eq.2) then
-            WRITE(numout,*)     &
-            &   ' riverine N inputs specified,                               jriver_n    = ', jriver_n
-         endif
-         if (jriver_si.eq.0) then
-            WRITE(numout,*)     &
-            &   ' *no* riverine Si input,                                    jriver_si   = ', jriver_si
-         elseif (jriver_si.eq.1) then
-            WRITE(numout,*)     &
-            &   ' riverine Si concentrations specified,                      jriver_si   = ', jriver_si
-         elseif (jriver_si.eq.2) then
-            WRITE(numout,*)     &
-            &   ' riverine Si inputs specified,                              jriver_si   = ', jriver_si
-         endif
-         if (jriver_c.eq.0) then
-            WRITE(numout,*)     &
-            &   ' *no* riverine C input,                                     jriver_c    = ', jriver_c
-         elseif (jriver_c.eq.1) then
-            WRITE(numout,*)     &
-            &   ' riverine C concentrations specified,                       jriver_c    = ', jriver_c
-         elseif (jriver_c.eq.2) then
-            WRITE(numout,*)     &
-            &   ' riverine C inputs specified,                               jriver_c    = ', jriver_c
-         endif
-         if (jriver_alk.eq.0) then
-            WRITE(numout,*)     &
-            &   ' *no* riverine alkalinity input,                            jriver_alk  = ', jriver_alk
-         elseif (jriver_alk.eq.1) then
-            WRITE(numout,*)     &
-            &   ' riverine alkalinity concentrations specified,              jriver_alk  = ', jriver_alk
-         elseif (jriver_alk.eq.2) then
-            WRITE(numout,*)     &
-            &   ' riverine alkalinity inputs specified,                      jriver_alk  = ', jriver_alk
-         endif
-         !! AXY (19/07/12): prevent (gross) stupidity on part of user
-         if (jriver_dep.lt.1.or.jriver_dep.ge.jpk) then
-            jriver_dep = 1
-         endif
-         WRITE(numout,*)     &
-         &   ' riverine input applied to down to depth k = ...            jriver_dep  = ', jriver_dep
-!!
-!! Miscellaneous
-         WRITE(numout,*) '=== Miscellaneous'
-         WRITE(numout,*)     &
-         &   ' diatom frustule dissolution rate,                          xsdiss      = ', xsdiss
-!!
-!! Gravitational sinking      
-         WRITE(numout,*) '=== Gravitational sinking'
-         WRITE(numout,*)     &
-         &   ' detritus gravitational sinking rate,                       vsed        = ', vsed
-         WRITE(numout,*)     & 
-         &   ' coefficient for Martin et al. (1987) remineralisation,     xhr         = ', xhr
-!!
-!! Non-Medusa parameters
-         WRITE(numout,*) '=== Non-Medusa parameters'
-         WRITE(numout,*)     & 
-         &   ' time coeff of POC in sediments,                            sedlam      = ', sedlam
-         WRITE(numout,*)     &
-         &   ' Sediment geol loss for POC,                                sedlostpoc  = ', sedlostpoc
-         WRITE(numout,*)     &
-         &   ' Vert layer for diagnostic of vertical flux,                jpkp        = ', jpkb
-!!
-!! UKESM1 - new diagnostics  !! Jpalm
-         WRITE(numout,*) '=== UKESM1-related parameters'
-         WRITE(numout,*)     &
-         &   ' include DMS diagnostic?,                                   jdms        = ', jdms
-!!
-      ENDIF
-!!
-!! Key depth positions (with thanks to Andrew Coward for bug-fixing this bit)
-      DO jk = 1,jpk
-         !! level thickness
-         fthk  = e3t_1d(jk)
-         !! level depth (top of level)
-         fdep  = gdepw_1d(jk)
-         !! level depth (bottom of level)
-         fdep1 = fdep + fthk
-         !!
-         if (fdep.lt.100.0.AND.fdep1.gt.100.0) then        !  100 m
-            i0100 = jk
-         elseif (fdep.lt.150.0.AND.fdep1.gt.150.0) then    !  150 m (for BASIN)
-            i0150 = jk
-         elseif (fdep.lt.200.0.AND.fdep1.gt.200.0) then    !  200 m
-            i0200 = jk
-         elseif (fdep.lt.500.0.AND.fdep1.gt.500.0) then    !  500 m
-            i0500 = jk
-         elseif (fdep.lt.1000.0.AND.fdep1.gt.1000.0) then  ! 1000 m
-            i1000 = jk
-         elseif (fdep1.lt.1100.0) then                     ! 1100 m (for Moore et al. sedimentary iron)
-            i1100 = jk
-         endif
-      enddo
-      !!
-      IF(lwp) THEN
-          WRITE(numout,*) '=== Position of key depths'
-          WRITE(numout,*)     & 
-          &   ' jk position of  100 m horizon                              i0100       = ', i0100
-          WRITE(numout,*)     &
-          &   ' jk position of  150 m horizon                              i0150       = ', i0150
-          WRITE(numout,*)     & 
-          &   ' jk position of  200 m horizon                              i0200       = ', i0200
-          WRITE(numout,*)     & 
-          &   ' jk position of  500 m horizon                              i0500       = ', i0500
-          WRITE(numout,*)     & 
-          &   ' jk position of 1000 m horizon                              i1000       = ', i1000
-          WRITE(numout,*)     & 
-          &   ' jk position of 1100 m horizon [*]                          i1100       = ', i1100
-          WRITE(numout,*) 'Got here ' , SIZE(friver_dep)
-          CALL flush(numout)
-      ENDIF
-
-#if defined key_roam
-
-      ! 1.4b namelist natroam : ROAM parameters
-      ! ---------------------------------------
-      
-      xthetaphy = 0.
-      xthetazoo = 0.
-      xthetanit = 0.
-      xthetarem = 0.
-      xo2min    = 0.
-
-      READ(numnatm,natroam)
-
-!! ROAM carbon, alkalinity and oxygen cycle parameters
-!!       xthetaphy :  oxygen evolution/consumption by phytoplankton
-!!       xthetazoo :  oxygen consumption by zooplankton
-!!       xthetanit :  oxygen consumption by nitrogen remineralisation
-!!       xthetarem :  oxygen consumption by carbon remineralisation
-!!       xo2min    :  oxygen minimum concentration
-
-      IF(lwp) THEN
-          WRITE(numout,*) 'natroam'
-          WRITE(numout,*) ' '
-!!
-!! ROAM carbon, alkalinity and oxygen cycle parameters
-          WRITE(numout,*) '=== ROAM carbon, alkalinity and oxygen cycle parameters'
-          WRITE(numout,*)     &
-          &   ' oxygen evolution/consumption by phytoplankton              xthetaphy   = ', xthetaphy
-          WRITE(numout,*)     &
-          &   ' oxygen consumption by zooplankton                          xthetazoo   = ', xthetazoo
-          WRITE(numout,*)     &
-          &   ' oxygen consumption by nitrogen remineralisation            xthetanit   = ', xthetanit
-          WRITE(numout,*)     &
-          &   ' oxygen consumption by carbon remineralisation              xthetarem   = ', xthetarem
-          WRITE(numout,*)     &
-          &   ' oxygen minimum concentration                               xo2min      = ', xo2min
-       ENDIF
-
-#endif
-
-      CALL flush(numout)
-
-      ! 1.5 namelist natopt : parameters for optic
-      ! ------------------------------------------
-
-      xkg0  = 0.
-      xkr0  = 0.
-      xkgp  = 0.
-      xkrp  = 0.
-      xlg   = 0.
-      xlr   = 0.
-      rpig  = 0.
-
-      READ(numnatm,natopt)
-
-      IF(lwp) THEN
-         WRITE(numout,*) 'natopt'
-         WRITE(numout,*) ' '
-         WRITE(numout,*) ' green   water absorption coeff  xkg0  = ',xkg0
-         WRITE(numout,*) ' red water absorption coeff      xkr0  = ',xkr0
-         WRITE(numout,*) ' pigment red absorption coeff    xkrp  = ',xkrp
-         WRITE(numout,*) ' pigment green absorption coeff  xkgp  = ',xkgp
-         WRITE(numout,*) ' green chl exposant              xlg   = ',xlg
-         WRITE(numout,*) ' red   chl exposant              xlr   = ',xlr
-         WRITE(numout,*) ' chla/chla+phea ratio            rpig  = ',rpig
-         WRITE(numout,*) ' '
-
-      ENDIF
-
-      IF(lwp) THEN
-         WRITE(numout,*) 'NaN check'
-         WRITE(numout,*) ' '
-         q1 = -1.
-         q2 = 0.
-         q3 = log(q1)
-         write (numout,*) 'q3 = ', q3
-         if ( ieee_is_nan( q3 ) ) then
-            write (numout,*) 'NaN detected'
-         else
-            write (numout,*) 'NaN not detected'
-         endif
-         WRITE(numout,*) ' '
-       ENDIF
-
-   END SUBROUTINE trc_nam_medusa
-   
-#else
-   !!----------------------------------------------------------------------
-   !!  Dummy module :                                             No MEDUSA
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE trc_nam_medusa                      ! Empty routine
-   END  SUBROUTINE  trc_nam_medusa
-#endif  
-
-   !!======================================================================
-END MODULE trcnam_medusa
Index: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcopt_medusa.F90
===================================================================
--- branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcopt_medusa.F90	(revision 7720)
+++ branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcopt_medusa.F90	(revision 7766)
@@ -156,4 +156,7 @@
 
       ! Euphotic layer depth
+      !! Jpalm -- 06-03-2017 -- add init xze, to avoid halo problems within the
+      !!                        writing process
+      xze(:,:) = 0.0
       DO jj = 1, jpj
          DO ji = 1, jpi
Index: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsed_medusa.F90
===================================================================
--- branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsed_medusa.F90	(revision 7720)
+++ branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsed_medusa.F90	(revision 7766)
@@ -41,5 +41,6 @@
 
    PUBLIC   trc_sed_medusa     ! called in ???
-   PUBLIC   trc_sed_medusa_sbc      
+   PUBLIC   trc_sed_medusa_sbc     
+   PUBLIC   trc_sed_medusa_dust 
 
    !! * Module variables
@@ -126,5 +127,4 @@
 
       !! AXY (20/11/14): alter this to report on first MEDUSA call
-      !! IF( kt == nit000 ) THEN
       IF( kt == nittrc000 ) THEN
          IF(lwp) WRITE(numout,*)
@@ -133,59 +133,4 @@
 	 IF(lwp) WRITE(numout,*) ' kt =',kt
       ENDIF
-
-      !! AXY (04/11/13): replace this with a call in trc_ini_medusa
-      !! AXY (25/02/10)
-      !! call routine for populating CCD array if this is the first time-step
-      !! IF( kt == nittrc000 ) CALL medusa_ccd( kt )
-
-      !! AXY (04/11/13): replace this with a call in trc_ini_medusa
-      !! AXY (26/01/12)
-      !! call routine for populating river arrays if this is the first time-step
-      !! IF( kt == nittrc000 ) CALL medusa_river( kt )
-
-      !! AXY (10/02/09)
-      !!IF( (jnt == 1) .and. (bdustfer) )  CALL trc_sed_medusa_sbc( kt )
-
-      !! JPALM -- 31-03-2016 -- rewrite trc_sed_medusa_sbc.
-      !! IF (kt == nittrc000 ) CALL trc_sed_medusa_sbc 
-
-      !! JPALM -- 20-07-2016 -- adapt dust forcing fields reading and conversion
-      !!                     To read dust dep in kg-dust/m2/s instead of g-Fe/m2/month 
-      !!                     So all forcings and coupling dust dep are in the same SI units
-      !!                     and then convert in mmol-Fe/m2/day
-
-      IF( bdustfer ) THEN
-        !! IF( kt == nittrc000 .OR. ( kt /= nittrc000 .AND. ntimes_dust > 1 ) ) THEN
-            CALL fld_read( kt, 1, sf_dust )
-            dust(:,:) = sf_dust(1)%fnow(:,:,1)
-         !! ENDIF
-      ELSEIF (lk_oasis) THEN
-         dust = Dust_in_cpl
-      ELSE
-         dust(:,:) = 0.0
-      ENDIF
-      !!
-      zirondep(:,:) = 0.e0     !! Initialisation of deposition variables
-      zirondep(:,:) = dust(:,:) * Fe_dust_mratio / xfe_mass * 1.e6 * 86400.  !! mmol-Fe/m2/d
-      
-      !! JPALM -- 20-07-2016 -- Zirondep and zsidep are not used.
-      !!                     So comment out the following lines. but keep them
-      !!                     as we may want to used them later on
-      !!================================================     
-      !!
-      !! zirondep(:,:,:) = 0.e0     !! Initialisation of deposition variables
-      !! zsidep  (:,:)   = 0.e0
-      !!
-      !! Iron and Si deposition at the surface
-      !! -------------------------------------
-      !!
-      !! DO jj = 1, jpj
-      !!    DO ji = 1, jpi
-      !!       zirondep(ji,jj,1) = (dustsolub * dust(ji,jj) / (55.85 * rmtss) + 3.e-10 / ryyss) &
-      !!       & * rfact2 / fse3t(ji,jj,1)
-      !!       zsidep  (ji,jj)   = 8.8 * 0.075 * dust(ji,jj) * rfact2 / &
-      !!       & (fse3t(ji,jj,1) * 28.1 * rmtss)
-      !!    END DO
-      !! END DO
 
       ! sedimentation of detrital nitrogen : upstream scheme
@@ -412,4 +357,83 @@
    END SUBROUTINE trc_sed_medusa_sbc
 
+   !! ======================================================================
+   !! ======================================================================
+   !! ======================================================================
+
+   !! AXY & JPALM (28/02/17)
+
+   SUBROUTINE trc_sed_medusa_dust( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_sed_medusa_dust  ***
+      !!
+      !! ** Purpose : compute current dust *before* trc_bio_medusa call
+      !!
+      !! ** Method  : does what it says on the tin
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index      
+
+      !! AXY (20/11/14): alter this to report on first MEDUSA call
+      IF( kt == nittrc000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) ' trc_sed_medusa_dust: MEDUSA dust timestep'
+         IF(lwp) WRITE(numout,*) ' ~~~~~~~'
+	 IF(lwp) WRITE(numout,*) ' kt =',kt
+      ENDIF
+
+      !! AXY (04/11/13): replace this with a call in trc_ini_medusa
+      !! AXY (25/02/10)
+      !! call routine for populating CCD array if this is the first time-step
+      !! IF( kt == nittrc000 ) CALL medusa_ccd( kt )
+
+      !! AXY (04/11/13): replace this with a call in trc_ini_medusa
+      !! AXY (26/01/12)
+      !! call routine for populating river arrays if this is the first time-step
+      !! IF( kt == nittrc000 ) CALL medusa_river( kt )
+
+      !! AXY (10/02/09)
+      !! IF( (jnt == 1) .and. (bdustfer) )  CALL trc_sed_medusa_sbc( kt )
+
+      !! JPALM -- 31-03-2016 -- rewrite trc_sed_medusa_sbc.
+      !! IF (kt == nittrc000 ) CALL trc_sed_medusa_sbc 
+
+      !! JPALM -- 20-07-2016 -- adapt dust forcing fields reading and conversion
+      !!                     To read dust dep in kg-dust/m2/s instead of g-Fe/m2/month 
+      !!                     So all forcings and coupling dust dep are in the same SI units
+      !!                     and then convert in mmol-Fe/m2/day
+
+      IF( bdustfer ) THEN
+            CALL fld_read( kt, 1, sf_dust )
+            dust(:,:) = sf_dust(1)%fnow(:,:,1)
+      ELSEIF (lk_oasis) THEN
+         dust = Dust_in_cpl
+      ELSE
+         dust(:,:) = 0.0
+      ENDIF
+      !!
+      zirondep(:,:) = 0.e0     !! Initialisation of deposition variables
+      zirondep(:,:) = dust(:,:) * Fe_dust_mratio / xfe_mass * 1.e6 * 86400.  !! mmol-Fe/m2/d
+      
+      !! JPALM -- 20-07-2016 -- Zirondep and zsidep are not used.
+      !!                     So comment out the following lines. but keep them
+      !!                     as we may want to used them later on
+      !!================================================     
+      !!
+      !! zirondep(:,:,:) = 0.e0     !! Initialisation of deposition variables
+      !! zsidep  (:,:)   = 0.e0
+      !!
+      !! Iron and Si deposition at the surface
+      !! -------------------------------------
+      !!
+      !! DO jj = 1, jpj
+      !!    DO ji = 1, jpi
+      !!       zirondep(ji,jj,1) = (dustsolub * dust(ji,jj) / (55.85 * rmtss) + 3.e-10 / ryyss) &
+      !!       & * rfact2 / fse3t(ji,jj,1)
+      !!       zsidep  (ji,jj)   = 8.8 * 0.075 * dust(ji,jj) * rfact2 / &
+      !!       & (fse3t(ji,jj,1) * 28.1 * rmtss)
+      !!    END DO
+      !! END DO
+
+   END SUBROUTINE trc_sed_medusa_dust
+
 #else
    !!======================================================================
@@ -425,2 +449,3 @@
    !!======================================================================
 END MODULE  trcsed_medusa
+
Index: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsms_medusa.F90
===================================================================
--- branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsms_medusa.F90	(revision 7720)
+++ branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsms_medusa.F90	(revision 7766)
@@ -58,14 +58,21 @@
       ENDIF
 
-      CALL trc_avg_medusa( kt )   ! rolling average module
+      CALL trc_avg_medusa( kt ) ! rolling average module
 # if defined key_debug_medusa
-         IF(lwp) WRITE(numout,*) ' MEDUSA done trc_avg_medusa'
-         CALL flush(numout)
+      IF(lwp) WRITE(numout,*) ' MEDUSA done trc_avg_medusa'
+      CALL flush(numout)
+# endif
+      
+      CALL trc_opt_medusa( kt ) ! optical model
+# if defined key_debug_medusa
+      IF(lwp) WRITE(numout,*) ' MEDUSA done trc_opt_medusa'
+      CALL flush(numout)
 # endif
 
-      CALL trc_opt_medusa( kt )   ! optical model
+      !! AXY & JPALM (28/02/17): call dust before trc_bio_medusa (because of coupling)
+      CALL trc_sed_medusa_dust( kt ) ! dust submodel
 # if defined key_debug_medusa
-         IF(lwp) WRITE(numout,*) ' MEDUSA done trc_opt_medusa'
-         CALL flush(numout)
+      IF(lwp) WRITE(numout,*) ' MEDUSA done trc_sed_medusa_dust'
+      CALL flush(numout)
 # endif
 
@@ -75,15 +82,15 @@
       CALL flush(numout)
 # else
-      CALL trc_bio_medusa( kt )   ! biological model
-# if defined key_debug_medusa
-         IF(lwp) WRITE(numout,*) ' MEDUSA done trc_bio_medusa'
-         CALL flush(numout)
-# endif
-
-      CALL trc_sed_medusa( kt )   ! sedimentation model
-# if defined key_debug_medusa
-         IF(lwp) WRITE(numout,*) ' MEDUSA done trc_sed_medusa'
-         CALL flush(numout)
-# endif
+      CALL trc_bio_medusa( kt ) ! biological model
+#  if defined key_debug_medusa
+      IF(lwp) WRITE(numout,*) ' MEDUSA done trc_bio_medusa'
+      CALL flush(numout)
+#  endif
+      
+      CALL trc_sed_medusa( kt ) ! sedimentation model
+#  if defined key_debug_medusa
+      IF(lwp) WRITE(numout,*) ' MEDUSA done trc_sed_medusa'
+      CALL flush(numout)
+#  endif
 # endif
 
@@ -103,2 +110,3 @@
    !!======================================================================
 END MODULE trcsms_medusa
+
Index: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 7720)
+++ branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 7766)
@@ -686,5 +686,5 @@
       !!----------------------------------------------------------------------
       INTEGER  :: jk, jn
-      REAL(wp) :: ztraf, zmin, zmax, zmean, zdrift
+      REAL(wp) :: ztraf, zmin, zmax, zmean, zdrift, areasf
       REAL(wp), DIMENSION(jpi,jpj) :: zvol
       !!----------------------------------------------------------------------
@@ -696,5 +696,6 @@
       ENDIF
       !
-         zvol(:,:) = e1e2t(:,:) * fse3t_a(:,:,1) * tmask(:,:,1)
+      zvol(:,:) = e1e2t(:,:) * fse3t_a(:,:,1) * tmask(:,:,1)
+      areasf = glob_sum(zvol(:,:))
       DO jn = 1, jptra
          ztraf = glob_sum( tra(:,:,1,jn) * zvol(:,:) )
@@ -705,5 +706,5 @@
             CALL mpp_max( zmax )      ! max over the global domain
          END IF
-         zmean  = ztraf / areatot
+         zmean  = ztraf / areasf
          IF(lwp) WRITE(numout,9001) jn, TRIM( ctrcnm(jn) ), zmean, zmin, zmax
       END DO
@@ -734,5 +735,6 @@
       zvol(:,:) = e1e2t(:,:) * fse3t_a(:,:,1) * tmask(:,:,1)
       ztraf = glob_sum( dgtr(:,:) * zvol(:,:) )
-      areasf = glob_sum(e1e2t(:,:) * tmask(:,:,1) )
+      !! areasf = glob_sum(e1e2t(:,:) * tmask(:,:,1) )
+      areasf = glob_sum(zvol(:,:))
       zmin  = MINVAL( dgtr(:,:), mask= ((tmask(:,:,1).NE.0.)) )
       zmax  = MAXVAL( dgtr(:,:), mask= ((tmask(:,:,1).NE.0.)) )
@@ -741,5 +743,5 @@
          CALL mpp_max( zmax )      ! max over the global domain
       END IF
-      zmean  = ztraf / areatot
+      zmean  = ztraf / areasf
       IF(lwp) WRITE(numout,9002) TRIM( names ), zmean, zmin, zmax
       !
Index: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/trcstp.F90
===================================================================
--- branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 7720)
+++ branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 7766)
@@ -146,5 +146,5 @@
             ztrai = ztrai + glob_sum( trn(:,:,:,jn) * cvol(:,:,:)   )
          END DO
-         IF( lwp ) WRITE(numstr,9300) kt,  ztrai / areatot
+         IF( numstr /= -1 ) WRITE(numstr,9300) kt,  ztrai / areatot
 9300     FORMAT(i10,e18.10)
       ENDIF
