Index: /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/AMM12/EXP00/iodef.xml
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/AMM12/EXP00/iodef.xml	(revision 6070)
+++ /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/AMM12/EXP00/iodef.xml	(revision 6071)
@@ -25,9 +25,66 @@
       <file_group id="1ts" output_freq="1ts"  output_level="10" enabled=".TRUE."/> <!-- 1 time step files -->
 
-      <file_group id="1h" output_freq="1h"  output_level="10" enabled=".TRUE."/> <!-- 1h files -->
+<!--old      <file_group id="1h" output_freq="1h"  output_level="10" enabled=".TRUE."/> old --> <!-- 1h files --> 
+   <!-- TMB files -->
+      <file_group id="1h" output_freq="1h"  output_level="10" enabled=".TRUE."  >                      
+
+   <file id="file9" name_suffix="_shelftmb_grid_T" description="ocean T grid variables" enabled=".TRUE." >
+     <field field_ref="top_temp"           name="votemper_top"   operation="instant" enabled=".TRUE." />
+     <field field_ref="mid_temp"           name="votemper_mid"   operation="instant" enabled=".TRUE." />
+     <field field_ref="bot_temp"           name="votemper_bot"   operation="instant" enabled=".TRUE." />
+     <field field_ref="sshnmasked"         name="sossheig"       operation="instant" enabled=".TRUE." /> 
+     <field field_ref="top_sal"            name="vosaline_top"   operation="instant" enabled=".TRUE." />
+     <field field_ref="mid_sal"            name="vosaline_mid"   operation="instant" enabled=".TRUE." />
+     <field field_ref="bot_sal"            name="vosaline_bot"   operation="instant" enabled=".TRUE." />
+   </file>
+
+   <file id="file10" name_suffix="_shelftmb_grid_U" description="TMB ocean U grid variables" enabled=".TRUE." >
+     <field field_ref="top_u"         name="vozocrtx_top"   operation="instant" />
+     <field field_ref="mid_u"         name="vozocrtx_mid"   operation="instant" />
+     <field field_ref="bot_u"         name="vozocrtx_bot"   operation="instant" />
+     <field field_ref="baro_u"        name="vobtcrtx"       operation="instant" />
+   </file>
+
+   <file id="file11" name_suffix="_shelftmb_grid_V" description="TMB ocean V grid variables" enabled=".TRUE." >
+     <field field_ref="top_v"         name="vomecrty_top"   operation="instant" />
+     <field field_ref="mid_v"         name="vomecrty_mid"   operation="instant" />
+     <field field_ref="bot_v"         name="vomecrty_bot"   operation="instant" />
+     <field field_ref="baro_v"        name="vobtcrty"       operation="instant" />
+   </file>
+
+      </file_group>
+
       <file_group id="2h" output_freq="2h"  output_level="10" enabled=".TRUE."/> <!-- 2h files -->
       <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=".TRUE."/> <!-- 6h files -->
+
+      <file_group id="25h_mean" output_freq="1d" output_level="10" enabled=".TRUE.">
+
+
+     <file id="file12" name_suffix="_25hourm_grid_T" description="ocean T grid variables, 25h meaned" enabled=".TRUE." >
+     <field field_ref="temper25h"   name="votemper"  operation="instant" enabled=".TRUE."/>
+     <field field_ref="tempis25h"   name="votempis"  operation="instant" enabled=".TRUE."/>
+     <field field_ref="salin25h"    name="vosaline"  operation="instant" enabled=".TRUE."/>
+     <field field_ref="ssh25h"      name="sossheig"  operation="instant" enabled=".TRUE."/>
+    </file>
+
+     <file id="file13" name_suffix="_25hourm_grid_U" description="ocean U grid variables, 25h meaned" enabled=".TRUE." >
+     <field field_ref="vozocrtx25h" name="vozocrtx"  operation="instant" enabled=".TRUE."/>
+    </file>
+
+     <file id="file14" name_suffix="_25hourm_grid_V" description="ocean V grid variables, 25h meaned" enabled=".TRUE." >
+     <field field_ref="vomecrty25h" name="vomecrty"  operation="instant" enabled=".TRUE."/>
+    </file>
+
+     <file id="file15" name_suffix="_25hourm_grid_W" description="ocean W grid variables, 25h meaned" enabled=".TRUE." >
+     <field field_ref="vomecrtz25h"          name="vomerctz"  operation="instant" enabled=".TRUE."/>
+     <field field_ref="woce"                 name="vovecrtz"  operation="instant" enabled=".TRUE."/>
+     <field field_ref="avt25h"               name="votkeavt"  operation="instant" enabled=".TRUE."/>
+     <field field_ref="avm25h"               name="votkeavm"  operation="instant" enabled=".TRUE."/>
+     <field field_ref="tke25h"               name="votke"  operation="instant" enabled=".TRUE."/>
+     <field field_ref="mxln25h"              name="mxln"   operation="instant" enabled=".TRUE."/>
+    </file>
+      </file_group>
      
       <file_group id="1d" output_freq="1d"  output_level="10" enabled=".TRUE."> <!-- 1d files -->
@@ -98,4 +155,5 @@
       <file_group id="5y"  output_freq="5y" output_level="10" enabled=".TRUE."/> <!-- real 5y files -->
       <file_group id="10y" output_freq="10y" output_level="10" enabled=".TRUE."/> <!-- real 10y files -->
+
 
    </file_definition>
Index: /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg	(revision 6070)
+++ /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg	(revision 6071)
@@ -446,2 +446,13 @@
 !-----------------------------------------------------------------------
 /
+!-----------------------------------------------------------------------
+&nam_diatmb  !  Top Middle Bottom Output
+!-----------------------------------------------------------------------
+   ln_diatmb  = .true.    !  Choose Top Middle and Bottom output or not
+/
+!-----------------------------------------------------------------------
+&nam_dia25h  !  25h Mean Output
+!-----------------------------------------------------------------------
+   ln_dia25h  = .true.    ! Choose 25h mean output or not
+/
+
Index: /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist_cfg
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist_cfg	(revision 6070)
+++ /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist_cfg	(revision 6071)
@@ -259,5 +259,5 @@
 /
 !-----------------------------------------------------------------------
-&namobs       !  observation usage                                      ('key_diaobs')
+&namobs       !  observation usage
 !-----------------------------------------------------------------------
 /
Index: /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/ORCA2_LIM3/cpp_ORCA2_LIM3.fcm
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/ORCA2_LIM3/cpp_ORCA2_LIM3.fcm	(revision 6070)
+++ /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/ORCA2_LIM3/cpp_ORCA2_LIM3.fcm	(revision 6071)
@@ -1,1 +1,1 @@
- bld::tool::fppkeys key_trabbl key_lim3 key_zdftke key_zdfddm key_zdftmx  key_mpp_mpi key_diaobs key_asminc key_iomput key_nosignedzero
+bld::tool::fppkeys key_trabbl key_lim3 key_zdftke key_zdfddm key_zdftmx key_iomput key_mpp_mpi key_asminc key_nosignedzero
Index: /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/ORCA2_LIM_CFC_C14b/EXP00/1_namelist_cfg
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/ORCA2_LIM_CFC_C14b/EXP00/1_namelist_cfg	(revision 6070)
+++ /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/ORCA2_LIM_CFC_C14b/EXP00/1_namelist_cfg	(revision 6071)
@@ -868,5 +868,4 @@
    ln_ignmis  = .true.     ! Logical switch for ignoring missing files   
                            !     endailyavtypes   ENACT daily average types                    
-   ln_grid_global = .true.
    ln_grid_search_lookup = .false.
 / 
Index: /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 6070)
+++ /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 6071)
@@ -40,5 +40,8 @@
          <field id="sstgrad2"     long_name="square of module of sst gradient"                                                        unit="degC2/m2" />
          <field id="sbt"          long_name="sea bottom temperature"                                                                  unit="degC"     />
-	 
+         <field id="sst_wl"       long_name="Delta SST of warm layer"                                                                 unit="degC"     />
+         <field id="sst_cs"       long_name="Delta SST of cool skin"                                                                  unit="degC"     />
+	 <field id="temp_3m"      long_name="temperature at 3m"                                                                       unit="degC"     />
+         
          <field id="sss"          long_name="sea surface salinity"             standard_name="sea_surface_salinity"   unit="1e-3" />
          <field id="sss2"         long_name="square of sea surface salinity"                                          unit="1e-6"  > sss * sss </field >
@@ -197,10 +200,16 @@
 
          <!-- * variable related to ice shelf forcing * -->
-         <field id="fwfisf"       long_name="Ice shelf melting"                             unit="kg/m2/s"  />
-         <field id="qisf"         long_name="Ice Shelf Heat Flux"                           unit="W/m2"     />
-         <field id="isfgammat"    long_name="transfert coefficient for isf (temperature)"   unit="m/s"      />
-         <field id="isfgammas"    long_name="transfert coefficient for isf (salinity)"      unit="m/s"      />
-         <field id="stbl"         long_name="salinity in the Losh tbl"                      unit="1e-3"     />
-         <field id="ttbl"         long_name="temperature in the Losh tbl"                   unit="degC"     />
+         <field id="fwfisf"       long_name="Ice shelf melting"                                            unit="Kg/m2/s"  />
+         <field id="qisf"         long_name="Ice Shelf Heat Flux"                                          unit="W/m2"     />
+         <field id="isfgammat"    long_name="transfert coefficient for isf (temperature) "                 unit="m/s"      />
+         <field id="isfgammas"    long_name="transfert coefficient for isf (salinity)    "                 unit="m/s"      />
+ 	 <field id="stbl"         long_name="salinity in the Losh tbl                    "                 unit="PSU"      />
+	 <field id="ttbl"         long_name="temperature in the Losh tbl                 "                 unit="C"        />
+	 <field id="utbl"         long_name="zonal current in the Losh tbl at T point    "                 unit="m/s"      />
+	 <field id="vtbl"         long_name="merid current in the Losh tbl at T point    "                 unit="m/s"      />
+	 <field id="thermald"     long_name="thermal driving of ice shelf melting        "                 unit="C"        />
+	 <field id="tfrz"         long_name="top freezing point (used to compute melt)   "                 unit="C"        />
+	 <field id="tinsitu"      long_name="top insitu temperature (used to cmpt melt)  "                 unit="C"        />
+	 <field id="ustar"        long_name="ustar at T point used in ice shelf melting  "                 unit="m/s"      />
 
          <!-- *_oce variables available with ln_blk_clio or ln_blk_core -->
@@ -1039,5 +1048,52 @@
       <field field_ref="woce"         name="wo"       long_name="ocean vertical velocity"  />
     </field_group>
-    
+
+       <!-- TMB diagnostic output -->
+    <field_group  id="1h_grid_T_tmb" grid_ref="grid_T_2D" operation="instant">
+     <field id="top_temp"           name="votemper_top"  unit="degC"  />
+     <field id="mid_temp"           name="votemper_mid"  unit="degC"  />
+     <field id="bot_temp"           name="votemper_bot"  unit="degC"  />
+     <field id="top_sal"            name="vosaline_top"  unit="psu"   />
+     <field id="mid_sal"            name="vosaline_mid"  unit="psu"   />
+     <field id="bot_sal"            name="vosaline_bot"  unit="psu"   />
+     <field id="sshnmasked"         name="sossheig"      unit="m"     /> 
+   </field_group>
+    <field_group  id="1h_grid_U_tmb" grid_ref="grid_U_2D" operation="instant">
+     <field id="top_u"           name="vozocrtx_top"  unit="m/s"  />
+     <field id="mid_u"           name="vozocrtx_mid"  unit="m/s"  />
+     <field id="bot_u"           name="vozocrtx_bot"  unit="m/s"  />
+     <field id="baro_u"          name="vobtcrtx"      unit="m/s"  />
+   </field_group>
+    <field_group  id="1h_grid_V_tmb" grid_ref="grid_V_2D" operation="instant">
+     <field id="top_v"           name="vomecrty_top"  unit="m/s"  />
+     <field id="mid_v"           name="vomecrty_mid"  unit="m/s"  />
+     <field id="bot_v"           name="vomecrty_bot"  unit="m/s"  />
+     <field id="baro_v"          name="vobtcrty"      unit="m/s"  />
+   </field_group>
+       <!-- 25h diagnostic output -->
+      <field_group id="25h_grid_T" grid_ref="grid_T_3D" operation="instant">
+   <field id="temper25h"         name="potential temperature 25h mean"    unit="degC" />
+   <field id="tempis25h"         name="insitu temperature 25h mean"    unit="degC" />
+   <field id="salin25h"          name="salinity 25h mean"                 unit="psu"  />
+   <field id="ssh25h"            name="sea surface height 25h mean"  grid_ref="grid_T_2D"      unit="m"    />
+      </field_group>
+
+      <field_group id="25h_grid_U" grid_ref="grid_U_3D" operation="instant" >
+   <field id="vozocrtx25h"         name="i current 25h mean"    unit="m/s"   />
+      </field_group>
+
+      <field_group id="25h_grid_V" grid_ref="grid_V_3D" operation="instant">
+   <field id="vomecrty25h"         name="j current 25h mean"    unit="m/s"    />
+      </field_group>
+
+      <field_group id="25h_grid_W" grid_ref="grid_W_3D" operation="instant">
+   <field id="vomecrtz25h"         name="k current 25h mean"                 unit="m/s"      />
+   <field id="avt25h"             name="vertical diffusivity25h mean"       unit="m2/s" />
+   <field id="avm25h"              name="vertical viscosity 25h mean"        unit="m2/s" />
+   <field id="tke25h"              name="turbulent kinetic energy 25h mean" />
+   <field id="mxln25h"             name="mixing length 25h mean"             unit="m" />
+      </field_group>
+
+
 
    </field_definition>
Index: /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 6070)
+++ /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 6071)
@@ -31,4 +31,5 @@
    nn_itend    =    5475   !  last  time step (std 5475)
    nn_date0    =  010101   !  date at nit_0000 (format yyyymmdd) used if ln_rstart=F or (ln_rstart=T and nn_rstctl=0 or 1)
+   nn_time0    =       0   !  initial time of day in hhmm
    nn_leapy    =       0   !  Leap year calendar (1) or not (0)
    ln_rstart   = .false.   !  start from rest (F) or from a restart file (T)
@@ -42,4 +43,5 @@
       cn_ocerst_out = "restart"  !  suffix of ocean restart name (output)
       cn_ocerst_outdir = "."     !  directory in which to write output ocean restarts
+   ln_iscpl    = .false.   !  cavity evolution forcing or coupling to ice sheet model
    nn_istate   =       0   !  output the initial state (1) or not (0)
    ln_rst_list = .false.   !  output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F)
@@ -132,4 +134,5 @@
    nn_msh      =    1      !  create (=1) a mesh file or not (=0)
    rn_hmin     =   -3.     !  min depth of the ocean (>0) or min number of ocean level (<0)
+   rn_isfhmin  =    1.00   !  treshold (m) to discriminate grounding ice to floating ice
    rn_e3zps_min=   20.     !  partial step thickness is set larger than the minimum of
    rn_e3zps_rat=    0.1    !  rn_e3zps_min and rn_e3zps_rat*e3t, with 0<rn_e3zps_rat<1
@@ -230,9 +233,10 @@
 !!   namsbc_cpl      CouPLed            formulation                     ("key_oasis3" )
 !!   namsbc_sas      StAndalone Surface module
-!!   namtra_qsr      penetrative solar radiation                        (ln_traqsr  =T)
-!!   namsbc_rnf      river runoffs                                      (ln_rnf     =T)
-!!   namsbc_isf      ice shelf melting/freezing                         (nn_isf     >0)
-!!   namsbc_apr      Atmospheric Pressure                               (ln_apr_dyn =T)
-!!   namsbc_ssr      sea surface restoring term (for T and/or S)        (ln_ssr     =T)
+!!   namtra_qsr      penetrative solar radiation
+!!   namsbc_rnf      river runoffs
+!!   namsbc_isf      ice shelf melting/freezing
+!!   namsbc_iscpl    coupling option between land ice model and ocean
+!!   namsbc_apr      Atmospheric Pressure
+!!   namsbc_ssr      sea surface restoring term (for T and/or S)
 !!   namsbc_alb      albedo parameters
 !!   namsbc_wave     external fields from wave model                    (ln_wave    =T)
@@ -279,9 +283,5 @@
                            !     =2 annual global mean of e-p-r set to zero
    ln_apr_dyn  = .false.   !  Patm gradient added in ocean & ice Eqs.   (T => fill namsbc_apr )
-   nn_isf      = 0         !  ice shelf melting/freezing                (/=0 => fill namsbc_isf)
-                           !  0 =no isf                  1 = presence of ISF
-                           !  2 = bg03 parametrisation   3 = rnf file for isf
-                           !  4 = ISF fwf specified
-                           !  option 1 and 4 need ln_isfcav = .true. (domzgr)
+   ln_isf      = .false.   !  ice shelf                                 (T   => fill namsbc_isf)
    ln_wave = .false.       !  coupling with surface wave                    (T => fill namsbc_wave)
    nn_lsm  = 0             !  =0 land/sea mask for input fields is not applied (keep empty land/sea mask filename field) ,
@@ -453,34 +453,42 @@
 &namsbc_isf    !  Top boundary layer (ISF)                              (nn_isf >0)
 !-----------------------------------------------------------------------
-!              ! file name ! frequency (hours) ! variable ! time interp. ! clim   ! 'yearly'/ ! weights  ! rotation !
-!              !           !  (if <0  months)  !   name   !   (logical)  ! (T/F)  ! 'monthly' ! filename ! pairing  !
+!              ! file name ! frequency (hours) ! variable ! time interpol. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
+!              !           !  (if <0  months)  !   name   !    (logical)   !  (T/F)  ! 'monthly' ! filename ! pairing  ! filename      !
 ! nn_isf == 4
-   sn_qisf      = 'rnfisf' ,         -12       ,'sohflisf',    .false.   , .true. , 'yearly'  ,  ''      ,   ''
-   sn_fwfisf    = 'rnfisf' ,         -12       ,'sowflisf',    .false.   , .true. , 'yearly'  ,  ''      ,   ''
+   sn_fwfisf   = 'rnfisf'  ,         -12       ,'sowflisf',     .false.    , .true.  , 'yearly'  ,  ''      ,   ''     , ''
 ! nn_isf == 3
-   sn_rnfisf    = 'runoffs',         -12       ,'sofwfisf',    .false.   , .true. , 'yearly'  ,  ''      ,   ''
+   sn_rnfisf   = 'rnfisf'  ,         -12       ,'sofwfisf',     .false.    , .true.  , 'yearly'  ,  ''      ,   ''     , ''
 ! nn_isf == 2 and 3
-   sn_depmax_isf = 'runoffs',        -12       ,'sozisfmax',   .false.   , .true. , 'yearly'  ,  ''      ,   ''
-   sn_depmin_isf = 'runoffs',        -12       ,'sozisfmin',   .false.   , .true. , 'yearly'  ,  ''      ,   ''
+   sn_depmax_isf = 'rnfisf' ,        -12       ,'sozisfmax' ,   .false.    , .true.  , 'yearly'  ,  ''      ,   ''     , ''
+   sn_depmin_isf = 'rnfisf' ,        -12       ,'sozisfmin' ,   .false.    , .true.  , 'yearly'  ,  ''      ,   ''     , ''
 ! nn_isf == 2
-   sn_Leff_isf = 'rnfisf'  ,           0       ,'Leff'    ,    .false.   , .true. , 'yearly'  ,  ''      ,   ''
-
+   sn_Leff_isf = 'rnfisf'  ,         -12       ,'Leff'    ,     .false.    , .true.  , 'yearly'  ,  ''      ,   ''     , ''
 ! for all case
-   ln_divisf   = .true.   ! apply isf melting as a mass flux or in the salinity trend. (maybe I should remove this option as for runoff?)
+   nn_isf      = 1         !  ice shelf melting/freezing
+                           !  1 = presence of ISF    2 = bg03 parametrisation 
+                           !  3 = rnf file for isf   4 = ISF fwf specified
+                           !  option 1 and 4 need ln_isfcav = .true. (domzgr)
 ! only for nn_isf = 1 or 2
    rn_gammat0  = 1.e-4    ! gammat coefficient used in blk formula
    rn_gammas0  = 1.e-4    ! gammas coefficient used in blk formula
+! only for nn_isf = 1 or 4
+   rn_hisf_tbl =  30.      ! thickness of the top boundary layer    (Losh et al. 2008)
+                          ! 0 => thickness of the tbl = thickness of the first wet cell
 ! only for nn_isf = 1
-   nn_isfblk   =  1       ! 1 ISOMIP ; 2 conservative (3 equation formulation, Jenkins et al. 1991 ??)
-   rn_hisf_tbl =  30.     ! thickness of the top boundary layer           (Losh et al. 2008)
-                          ! 0 => thickness of the tbl = thickness of the first wet cell
-   ln_conserve = .true.   ! conservative case (take into account meltwater advection)
+   nn_isfblk   = 1        ! 1 ISOMIP  like: 2 equations formulation (Hunter et al., 2006)
+                          ! 2 ISOMIP+ like: 3 equations formulation (Asay-Davis et al., 2015)
    nn_gammablk = 1        ! 0 = cst Gammat (= gammat/s)
                           ! 1 = velocity dependend Gamma (u* * gammat/s)  (Jenkins et al. 2010)
-                          !     if you want to keep the cd as in global config, adjust rn_gammat0 to compensate
-                          ! 2 = velocity and stability dependent Gamma    Holland et al. 1999
-/
-!-----------------------------------------------------------------------
-&namsbc_apr    !   Atmospheric pressure forcing (in ocean or bulk)      (ln_apr_dyn=T)
+                          ! 2 = velocity and stability dependent Gamma    (Holland et al. 1999)
+/
+!-----------------------------------------------------------------------
+&namsbc_iscpl  !   land ice / ocean coupling option                     
+!-----------------------------------------------------------------------
+   nn_drown  = 10       ! number of iteration of the extrapolation loop (fill the new wet cells)
+   ln_hsb    = .false.  ! activate conservation module (conservation exact after a time of rn_fiscpl)
+   nn_fiscpl = 43800    ! (number of time step) conservation period (maybe should be fix to the coupling frequencey of restart frequency)
+/
+!-----------------------------------------------------------------------
+&namsbc_apr    !   Atmospheric pressure used as ocean forcing or in bulk
 !-----------------------------------------------------------------------
 !              !  file name ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
@@ -1067,4 +1075,5 @@
                            !     (no physical validity of the results)
    nn_timing   =    0      !  timing by routine activated (=1) creates timing.output file, or not (=0)
+   nn_diacfl   =    0      !  Write out cfl diagnostics (=1) in cfl_diagnostics.ascii, or not (=0)
 /
 !-----------------------------------------------------------------------
@@ -1176,62 +1185,44 @@
 !!               ***  Observation & Assimilation  ***
 !!======================================================================
-!!   namobs       observation and model comparison                      ('key_diaobs')
+!!   namobs       observation and model comparison
 !!   nam_asminc   assimilation increments                               ('key_asminc')
 !!======================================================================
 !
 !-----------------------------------------------------------------------
-&namobs       !  observation usage switch                               ('key_diaobs')
-!-----------------------------------------------------------------------
-   ln_t3d     = .false.    ! Logical switch for T profile observations
-   ln_s3d     = .false.    ! Logical switch for S profile observations
-   ln_ena     = .false.    ! Logical switch for ENACT insitu data set
-   ln_cor     = .false.    ! Logical switch for Coriolis insitu data set
-   ln_profb   = .false.    ! Logical switch for feedback insitu data set
-   ln_sla     = .false.    ! Logical switch for SLA observations
-   ln_sladt   = .false.    ! Logical switch for AVISO SLA data
-   ln_slafb   = .false.    ! Logical switch for feedback SLA data
-   ln_ssh     = .false.    ! Logical switch for SSH observations
-   ln_sst     = .false.    ! Logical switch for SST observations
-   ln_reysst  = .false.    ! Logical switch for Reynolds observations
-   ln_ghrsst  = .false.    ! Logical switch for GHRSST observations
-   ln_sstfb   = .false.    ! Logical switch for feedback SST data
-   ln_sss     = .false.    ! Logical switch for SSS observations
-   ln_seaice  = .false.    ! Logical switch for Sea Ice observations
-   ln_vel3d   = .false.    ! Logical switch for velocity observations
-   ln_velavcur= .false     ! Logical switch for velocity daily av. cur.
-   ln_velhrcur= .false     ! Logical switch for velocity high freq. cur.
-   ln_velavadcp=.false.    ! Logical switch for velocity daily av. ADCP
-   ln_velhradcp=.false.    ! Logical switch for velocity high freq. ADCP
-   ln_velfb   = .false.    ! Logical switch for feedback velocity data
-   ln_grid_global=.false.  ! Global distribtion of observations
-   ln_grid_search_lookup = .false. !  Logical switch for obs grid search w/lookup table
-   grid_search_file = 'grid_search'  !  Grid search lookup file header
+&namobs       !  observation usage switch
+!-----------------------------------------------------------------------
+   ln_diaobs  = .false.             ! Logical switch for the observation operator
+   ln_t3d     = .false.             ! Logical switch for T profile observations
+   ln_s3d     = .false.             ! Logical switch for S profile observations
+   ln_sla     = .false.             ! Logical switch for SLA observations
+   ln_sst     = .false.             ! Logical switch for SST observations
+   ln_sic     = .false.             ! Logical switch for Sea Ice observations
+   ln_vel3d   = .false.             ! Logical switch for velocity observations
+   ln_altbias = .false.             ! Logical switch for altimeter bias correction
+   ln_nea     = .false.             ! Logical switch for rejection of observations near land
+   ln_grid_global = .true.          ! Logical switch for global distribution of observations
+   ln_grid_search_lookup = .false.  ! Logical switch for obs grid search w/lookup table
+   ln_ignmis  = .true.              ! Logical switch for ignoring missing files
+   ln_s_at_t  = .false.             ! Logical switch for computing model S at T obs if not there
+   ln_sstnight = .false.            ! Logical switch for calculating night-time average for SST obs
 ! All of the *files* variables below are arrays. Use namelist_cfg to add more files
-   enactfiles = 'enact.nc' !  ENACT input observation file names (specify full array in namelist_cfg)
-   coriofiles = 'corio.nc' !  Coriolis input observation file name
-   profbfiles = 'profiles_01.nc' ! Profile feedback input observation file name
-   ln_profb_enatim = .false !        Enact feedback input time setting switch
-   slafilesact = 'sla_act.nc' !  Active SLA input observation file names
-   slafilespas = 'sla_pass.nc' ! Passive SLA input observation file names
-   slafbfiles = 'sla_01.nc' ! slafbfiles: Feedback SLA input observation file names
-   sstfiles = 'ghrsst.nc'   ! GHRSST input observation file names
-   sstfbfiles = 'sst_01.nc' ! Feedback SST input observation file names
-   seaicefiles = 'seaice_01.nc' ! Sea Ice input observation file names
-   velavcurfiles  = 'velavcurfile.nc'  ! Vel. cur. daily av. input file name
-   velhrcurfiles  = 'velhrcurfile.nc'  ! Vel. cur. high freq. input file name
-   velavadcpfiles = 'velavadcpfile.nc' ! Vel. ADCP daily av. input file name
-   velhradcpfiles = 'velhradcpfile.nc' ! Vel. ADCP high freq. input file name
-   velfbfiles = 'velfbfile.nc' ! Vel. feedback input observation file name
-   dobsini = 20000101.000000  !  Initial date in window YYYYMMDD.HHMMSS
-   dobsend = 20010101.000000  !  Final date in window YYYYMMDD.HHMMSS
-   n1dint = 0  !               Type of vertical interpolation method
-   n2dint = 0  !               Type of horizontal interpolation method
-   ln_nea = .false.   !        Rejection of observations near land switch
-   nmsshc     = 0     !        MSSH correction scheme
-   mdtcorr = 1.61     !        MDT  correction
-   mdtcutoff = 65.0   !        MDT cutoff for computed correction
-   ln_altbias = .false.    ! Logical switch for alt bias
-   ln_ignmis  = .true.     ! Logical switch for ignoring missing files
-   endailyavtypes = 820    ! ENACT daily average types - array (use namelist_cfg to set more values)
+   cn_profbfiles = 'profiles_01.nc'    ! Profile feedback input observation file names
+   cn_slafbfiles = 'sla_01.nc'         ! SLA feedback input observation file names
+   cn_sstfbfiles = 'sst_01.nc'         ! SST feedback input observation file names
+   cn_sicfbfiles = 'sic_01.nc'         ! SIC feedback input observation file names
+   cn_velfbfiles = 'vel_01.nc'         ! Velocity feedback input observation file names
+   cn_altbiasfile = 'altbias.nc'       ! Altimeter bias input file name
+   cn_gridsearchfile = 'gridsearch.nc' ! Grid search file name
+   rn_gridsearchres = 0.5              ! Grid search resolution
+   rn_dobsini = 00010101.000000        ! Initial date in window YYYYMMDD.HHMMSS
+   rn_dobsend = 00010102.000000        ! Final date in window YYYYMMDD.HHMMSS
+   nn_1dint = 0                        ! Type of vertical interpolation method
+   nn_2dint = 0                        ! Type of horizontal interpolation method
+   nn_msshc = 0                        ! MSSH correction scheme
+   rn_mdtcorr = 1.61                   ! MDT  correction
+   rn_mdtcutoff = 65.0                 ! MDT cutoff for computed correction
+   nn_profdavtypes = -1                ! Profile daily average types - array
+   ln_sstbias = .false. 
+   cn_sstbias_files = 'sstbias.nc'
 /
 !-----------------------------------------------------------------------
@@ -1253,2 +1244,45 @@
     nn_divdmp = 0          !  Number of iterations of divergence damping operator
 /
+!-----------------------------------------------------------------------
+&namsbc_wave   ! External fields from wave model
+!-----------------------------------------------------------------------
+!              !  file name  ! frequency (hours) ! variable     ! time interp. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
+!              !             !  (if <0  months)  !   name       !   (logical)  !  (T/F)  ! 'monthly' ! filename ! pairing  ! filename      !
+   sn_cdg      =  'cdg_wave' ,        1          , 'drag_coeff' ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
+   sn_usd      =  'sdw_wave' ,        1          , 'u_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
+   sn_vsd      =  'sdw_wave' ,        1          , 'v_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
+   sn_wn       =  'sdw_wave' ,        1          , 'wave_num'   ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
+!
+   cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files
+/
+!-----------------------------------------------------------------------
+&namdyn_nept  !   Neptune effect (simplified: lateral and vertical diffusions removed)
+!-----------------------------------------------------------------------
+   ! Suggested lengthscale values are those of Eby & Holloway (1994) for a coarse model
+   ln_neptsimp       = .false.  ! yes/no use simplified neptune
+
+   ln_smooth_neptvel = .false.  ! yes/no smooth zunep, zvnep
+   rn_tslse          =  1.2e4   ! value of lengthscale L at the equator
+   rn_tslsp          =  3.0e3   ! value of lengthscale L at the pole
+   ! Specify whether to ramp down the Neptune velocity in shallow
+   ! water, and if so the depth range controlling such ramping down
+   ln_neptramp       = .true.   ! ramp down Neptune velocity in shallow water
+   rn_htrmin         =  100.0   ! min. depth of transition range
+   rn_htrmax         =  200.0   ! max. depth of transition range
+/
+!-----------------------------------------------------------------------
+&namdiu !   Cool skin and warm layer modesl
+!-----------------------------------------------------------------------
+   ln_diurnal = .false.
+   ln_diurnal_only = .false.
+/
+!-----------------------------------------------------------------------
+&nam_diatmb  !  Top Middle Bottom Output
+!-----------------------------------------------------------------------
+   ln_diatmb  = .false.    !  Choose Top Middle and Bottom output or not
+/
+!-----------------------------------------------------------------------
+&nam_dia25h  !  25h Mean Output
+!-----------------------------------------------------------------------
+   ln_dia25h  = .false.    ! Choose 25h mean output or not
+/
Index: /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/cfg.txt
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/cfg.txt	(revision 6070)
+++ /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/cfg.txt	(revision 6071)
@@ -7,6 +7,6 @@
 AMM12 OPA_SRC
 ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC
-GYRE OPA_SRC
 ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC
 ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC
 ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC
+GYRE OPA_SRC
Index: /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/uspcfg.txt
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/uspcfg.txt	(revision 6070)
+++ /branches/2015/dev_merge_2015/NEMOGCM/CONFIG/uspcfg.txt	(revision 6071)
@@ -1,2 +1,3 @@
 ORCA1_CICE # ORCA2_LIM # OPA_SRC TOP_SRC  # http://gws-access.ceda.ac.uk/public/nemo/uspconfigs/ORCA1_CICE/v3.6.0/ORCA1_CICE_ctl.txt
 ISOMIP     # GYRE      # OPA_SRC          # http://gws-access.ceda.ac.uk/public/nemo/uspconfigs/ISOMIP/v3.6.0/ISOMIP_ctl.txt
+ISOMIP_r5151_UKMO_ISF     # GYRE      # OPA_SRC          # http://gws-access.ceda.ac.uk/public/nemo/uspconfigs/ISOMIP_r5151_UKMO_ISF/v3.6.0/ISOMIP_ctl.txt
Index: /branches/2015/dev_merge_2015/NEMOGCM/SETTE/BATCH_TEMPLATE/batch-XC40_METO
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/SETTE/BATCH_TEMPLATE/batch-XC40_METO	(revision 6070)
+++ /branches/2015/dev_merge_2015/NEMOGCM/SETTE/BATCH_TEMPLATE/batch-XC40_METO	(revision 6071)
@@ -2,7 +2,7 @@
 #!
 #PBS -N nemo_sette
-#PBS -l walltime=00:30:00
+#PBS -l walltime=00:08:00
 #PBS -j oe
-#PBS -q QUEUE
+#PBS -q normal
 #PBS -l SELECT
 
@@ -77,12 +77,6 @@
 # Run SPMD case
 #
-    if [ $OCORES -gt 32 ] ; then
-        #Still more than one node so run in parallel queue
-        echo time aprun -b  -n $OCORES -N $O_PER_NODE ./opa
-             time aprun -b  -n $OCORES -N $O_PER_NODE ./opa
-    else
-       echo time mpiexec -n $OCORES ./opa
-            time mpiexec -n $OCORES ./opa
-    fi
+    echo time aprun -b  -n $OCORES -N $O_PER_NODE ./opa
+         time aprun -b  -n $OCORES -N $O_PER_NODE ./opa
   fi
 #
Index: /branches/2015/dev_merge_2015/NEMOGCM/SETTE/input_ORCA2_LIM.cfg
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/SETTE/input_ORCA2_LIM.cfg	(revision 6070)
+++ /branches/2015/dev_merge_2015/NEMOGCM/SETTE/input_ORCA2_LIM.cfg	(revision 6071)
@@ -1,1 +1,1 @@
-ORCA2_LIM_nemo_v3.6.tar ORCA2_LIM_nemo_v3.6
+ORCA2_LIM_nemo_v3.6st.tar ORCA2_LIM_nemo_v3.6
Index: /branches/2015/dev_merge_2015/NEMOGCM/SETTE/prepare_job.sh
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/SETTE/prepare_job.sh	(revision 6070)
+++ /branches/2015/dev_merge_2015/NEMOGCM/SETTE/prepare_job.sh	(revision 6071)
@@ -266,17 +266,7 @@
                                 echo NB_XNODES=${NB_XNODES} 
                                 echo Total NB_NODES=${NB_NODES}
-                                if [ ${NB_NODES} -eq 1 ] ; then
-                                   QUEUE=shared
-                                   #Not using XIOS in detatched mode and using less than one node so should be ok on shared node
-                                   #Load snplauncher module to allow use of mpiexec
-                                   SELECT="select=1:ncpus=$((NXIO_PROC + NB_PROC))":mem=15GB
-                                   module load cray-snplauncher
-                                   echo 'Shared Queue'
-                                else
-                                   QUEUE=normal
-                                   SELECT="select=$NB_NODES"
-                                   module unload cray-snplauncher #Make sure snplauncher module is not loaded
-                                   echo 'Normal Queue'
-                                fi
+                                QUEUE=normal
+                                SELECT="select=$NB_NODES"
+                                module unload cray-snplauncher #Make sure snplauncher module is not loaded
                                 ;;
                         openmpi_NAVITI_MERCATOR)
@@ -330,6 +320,5 @@
                     ;;
               XC40_METO*)
-                    cat run_sette_test.job | sed -e"s/QUEUE/${QUEUE}/" \
-                                                 -e"s/SELECT/${SELECT}/" > run_sette_test1.job
+                    cat run_sette_test.job | sed -e"s/SELECT/${SELECT}/" > run_sette_test1.job
                     mv run_sette_test1.job run_sette_test.job
                     ;;
Index: /branches/2015/dev_merge_2015/NEMOGCM/SETTE/sette.sh
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/SETTE/sette.sh	(revision 6070)
+++ /branches/2015/dev_merge_2015/NEMOGCM/SETTE/sette.sh	(revision 6071)
@@ -1137,5 +1137,5 @@
     export TEST_NAME="REPRO_4_4"
     cd ${CONFIG_DIR}
-    . ./makenemo -m ${CMP_NAM} -n ORCA2_LIM_OBS -r ORCA2_LIM -j 8 add_key "key_mpp_rep key_diaobs key_asminc" del_key ${DEL_KEYS}
+    . ./makenemo -m ${CMP_NAM} -n ORCA2_LIM_OBS -r ORCA2_LIM -j 8 add_key "key_mpp_rep key_asminc" del_key ${DEL_KEYS}
     cd ${SETTE_DIR}
     . ./param.cfg
@@ -1154,13 +1154,11 @@
     set_namelist namelist_cfg jpnj 4
     set_namelist namelist_cfg jpnij 16
-#   set_namelist namelist_cfg nn_solv 2   
+    set_namelist namelist_cfg ln_diaobs .true.
     set_namelist namelist_cfg ln_t3d .true.
     set_namelist namelist_cfg ln_s3d .true.
-    set_namelist namelist_cfg ln_profb .true.
     set_namelist namelist_cfg ln_sst .true.
-    set_namelist namelist_cfg ln_sstfb .true.
     set_namelist namelist_cfg ln_sla .true.
-    set_namelist namelist_cfg ln_slafb .true.
-    set_namelist namelist_cfg ln_seaice .false.
+    set_namelist namelist_cfg ln_sic .true.
+    set_namelist namelist_cfg ln_vel3d .true.
     set_namelist namelist_cfg ln_bkgwri .true.
     set_namelist namelist_cfg ln_trainc .true.
@@ -1193,13 +1191,11 @@
     set_namelist namelist_cfg jpnj 8
     set_namelist namelist_cfg jpnij 16
-#   set_namelist namelist_cfg nn_solv 2
+    set_namelist namelist_cfg ln_diaobs .true.
     set_namelist namelist_cfg ln_t3d .true.
     set_namelist namelist_cfg ln_s3d .true.
-    set_namelist namelist_cfg ln_profb .true.
     set_namelist namelist_cfg ln_sst .true.
-    set_namelist namelist_cfg ln_sstfb .true.
     set_namelist namelist_cfg ln_sla .true.
-    set_namelist namelist_cfg ln_slafb .true.
-    set_namelist namelist_cfg ln_seaice .false.
+    set_namelist namelist_cfg ln_sic .true.
+    set_namelist namelist_cfg ln_vel3d .true.
     set_namelist namelist_cfg ln_bkgwri .true.
     set_namelist namelist_cfg ln_trainc .true.
Index: /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obs_prof_io.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obs_prof_io.F90	(revision 6071)
+++ /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obs_prof_io.F90	(revision 6071)
@@ -0,0 +1,32 @@
+MODULE obs_prof_io
+   !!======================================================================
+   !!                       ***  MODULE obs_prof_io  ***
+   !! Observation operators : I/O for ENACT and Coriolis files
+   !!======================================================================
+   !! History : 
+   !!             !  09-01  (K. Mogensen) Initial version
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   read_enactfile    :  Read a obfbdata structure from an ENACT file.
+   !!   read_coriofile    :  Read a obfbdata structure from an Coriolis file.
+   !!----------------------------------------------------------------------
+   USE par_kind
+   USE obs_utils
+   USE obs_fbm
+   USE obs_conv
+   USE netcdf
+   IMPLICIT NONE
+
+   INTEGER, PARAMETER :: imaxlev = 10000
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: obs_prof_io.F90 2287 2010-10-18 07:53:52Z smasson $
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+#include "obsprof_io.h90"
+
+END MODULE obs_prof_io
Index: /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obs_sla_io.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obs_sla_io.F90	(revision 6071)
+++ /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obs_sla_io.F90	(revision 6071)
@@ -0,0 +1,30 @@
+MODULE obs_sla_io
+   !!======================================================================
+   !!                       ***  MODULE obs_sla_io  ***
+   !! Observation operators : I/O for AVISO SLA files
+   !!======================================================================
+   !! History : 
+   !!             !  09-01  (K. Mogensen) Initial version
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   read_enactfile    :  Read a obfbdata structure from an AVISO SLA file
+   !!----------------------------------------------------------------------
+   USE par_kind
+   USE obs_utils
+   USE obs_fbm
+   use obs_sla_types
+   USE netcdf
+
+   IMPLICIT NONE
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: obs_sla_io.F90 2287 2010-10-18 07:53:52Z smasson $
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+#include "obssla_io.h90"
+
+END MODULE obs_sla_io
Index: /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obs_sla_types.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obs_sla_types.F90	(revision 6071)
+++ /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obs_sla_types.F90	(revision 6071)
@@ -0,0 +1,19 @@
+MODULE obs_sla_types
+   !!======================================================================
+   !!                       ***  MODULE obs_sla_type  ***
+   !! Observation operators : Satellite identifiers for SLA data
+   !!======================================================================
+   !! History : 
+   !!             !  09-01  (K. Mogensen) Initial version
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: obs_sla_types.F90 2287 2010-10-18 07:53:52Z smasson $
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+   IMPLICIT NONE
+
+#include "obssla_types.h90"
+
+END MODULE obs_sla_types
Index: /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obs_vel_io.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obs_vel_io.F90	(revision 6071)
+++ /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obs_vel_io.F90	(revision 6071)
@@ -0,0 +1,32 @@
+MODULE obs_vel_io
+   !!======================================================================
+   !!                       ***  MODULE obs_vel_io  ***
+   !! Observation operators : I/O for TAO velocity data
+   !!======================================================================
+   !! History : 
+   !!             !  09-01  (K. Mogensen) Initial version based on obs_read_taovel
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   read_taofile    :  Read a obfbdata structure from an TAO velocity file.
+   !!----------------------------------------------------------------------
+   USE par_kind
+   USE obs_utils
+   USE obs_fbm
+   USE obs_conv
+   USE in_out_manager
+   USE netcdf
+   IMPLICIT NONE
+
+   INTEGER, PARAMETER :: imaxlev = 10000
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: obs_vel_io.F90 2287 2010-10-18 07:53:52Z smasson $
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+#include "obsvel_io.h90"
+
+END MODULE obs_vel_io
Index: /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obsprof_io.h90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obsprof_io.h90	(revision 6071)
+++ /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obsprof_io.h90	(revision 6071)
@@ -0,0 +1,834 @@
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: obsprof_io.h90 2287 2010-10-18 07:53:52Z smasson $
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+   SUBROUTINE read_enactfile( cdfilename, inpfile, kunit, ldwp, ldgrid )
+      !!---------------------------------------------------------------------
+      !!
+      !!                     ** ROUTINE read_enactfile **
+      !!
+      !! ** Purpose : Read from file the profile ENACT observations.
+      !!
+      !! ** Method  : The data file is a NetCDF file. 
+      !!
+      !! ** Action  :
+      !!
+      !! History : 
+      !!          ! 09-01 (K. Mogensen) Original based on old versions
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      CHARACTER(LEN=*) :: cdfilename ! Input filename
+      TYPE(obfbdata)   :: inpfile    ! Output obfbdata structure
+      INTEGER          :: kunit      ! Unit for output
+      LOGICAL          :: ldwp       ! Print info
+      LOGICAL          :: ldgrid     ! Save grid info in data structure
+      !! * Local declarations
+      INTEGER :: iobs                ! Number of observations
+      INTEGER :: ilev                      ! Number of levels
+      INTEGER :: i_file_id
+      INTEGER :: i_obs_id
+      INTEGER :: i_lev_id
+      INTEGER :: i_phi_id
+      INTEGER :: i_lam_id
+      INTEGER :: i_depth_id
+      INTEGER :: i_var_id
+      INTEGER :: i_pl_num_id
+      INTEGER :: i_reference_date_time_id
+      INTEGER :: i_format_version_id
+      INTEGER :: i_juld_id
+      INTEGER :: i_data_type_id
+      INTEGER :: i_wmo_inst_type_id
+      INTEGER :: i_qc_var_id
+      INTEGER :: i_dc_ref_id
+      INTEGER :: i_qc_flag_id
+      CHARACTER(LEN=40) :: cl_fld_lam
+      CHARACTER(LEN=40) :: cl_fld_phi
+      CHARACTER(LEN=40) :: cl_fld_depth 
+      CHARACTER(LEN=40) :: cl_fld_var_tp
+      CHARACTER(LEN=40) :: cl_fld_var_s
+      CHARACTER(LEN=40) :: cl_fld_var_ti
+      CHARACTER(LEN=40) :: cl_fld_var_juld_qc
+      CHARACTER(LEN=40) :: cl_fld_var_pos_qc
+      CHARACTER(LEN=40) :: cl_fld_var_depth_qc 
+      CHARACTER(LEN=40) :: cl_fld_var_qc_t
+      CHARACTER(LEN=40) :: cl_fld_var_qc_s
+      CHARACTER(LEN=40) :: cl_fld_var_prof_qc_t
+      CHARACTER(LEN=40) :: cl_fld_var_prof_qc_s
+      CHARACTER(LEN=40) :: cl_fld_reference_date_time
+      CHARACTER(LEN=40) :: cl_fld_juld
+      CHARACTER(LEN=40) :: cl_fld_data_type
+      CHARACTER(LEN=40) :: cl_fld_pl_num
+      CHARACTER(LEN=40) :: cl_fld_format_version
+      CHARACTER(LEN=40) :: cl_fld_wmo_inst_type
+      CHARACTER(LEN=40) :: cl_fld_qc_flags_profiles
+      CHARACTER(LEN=40) :: cl_fld_qc_flags_levels
+
+      CHARACTER(LEN=14), PARAMETER :: cl_name = 'read_enactfile'
+      CHARACTER(LEN=16)            :: cl_data_type = ''
+      CHARACTER(LEN=4 )            :: cl_format_version = ''
+      INTEGER, DIMENSION(1) :: istart1, icount1
+      INTEGER, DIMENSION(2) :: istart2, icount2
+      CHARACTER(len=imaxlev) :: clqc
+      CHARACTER(len=1) :: cqc
+      INTEGER :: ji, jk
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: iqc1
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: iqc2
+
+      !-----------------------------------------------------------------------
+      ! Initialization
+      !-----------------------------------------------------------------------
+
+      cl_fld_lam                 = 'LONGITUDE'
+      cl_fld_phi                 = 'LATITUDE'
+      cl_fld_depth               = 'DEPH_CORRECTED'
+      cl_fld_reference_date_time = 'REFERENCE_DATE_TIME'
+      cl_fld_juld                = 'JULD'
+      cl_fld_data_type           = 'DATA_TYPE'
+      cl_fld_format_version      = 'FORMAT_VERSION'
+      cl_fld_wmo_inst_type       = 'WMO_INST_TYPE'
+      cl_fld_pl_num              = 'PLATFORM_NUMBER'
+
+      cl_fld_var_qc_t            = 'POTM_CORRECTED_QC'
+      cl_fld_var_prof_qc_t       = 'PROFILE_POTM_QC'
+      cl_fld_var_tp              = 'POTM_CORRECTED'
+      cl_fld_var_qc_s            = 'PSAL_CORRECTED_QC'
+      cl_fld_var_prof_qc_s       = 'PROFILE_PSAL_QC'
+      cl_fld_var_s               = 'PSAL_CORRECTED'
+      cl_fld_var_depth_qc        = 'DEPH_CORRECTED_QC'
+      cl_fld_var_juld_qc         = 'JULD_QC'
+      cl_fld_var_pos_qc          = 'POSITION_QC'
+      cl_fld_var_ti              = 'TEMP'
+      cl_fld_qc_flags_profiles   = 'QC_FLAGS_PROFILES'
+      cl_fld_qc_flags_levels     = 'QC_FLAGS_LEVELS'
+
+      icount1(1) = 1 
+
+      !-----------------------------------------------------------------------
+      ! Open file
+      !-----------------------------------------------------------------------
+
+      CALL chkerr( nf90_open( TRIM( cdfilename ), nf90_nowrite, &
+            &      i_file_id ),           cl_name, __LINE__ )
+
+      !-----------------------------------------------------------------------
+      ! Read the heading of the file
+      !-----------------------------------------------------------------------
+
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_data_type,      &
+         &         i_data_type_id ),      cl_name, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_data_type_id,        &
+         &         cl_data_type ),        cl_name, __LINE__ )
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_format_version, &
+         &         i_format_version_id ), cl_name, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_format_version_id,   &
+         &         cl_format_version ),   cl_name, __LINE__ )
+      
+      CALL str_c_to_for( cl_data_type )
+      CALL str_c_to_for( cl_format_version )
+      
+      IF(ldwp)WRITE(kunit,*)
+      IF(ldwp)WRITE(kunit,*) ' read_enactfile :' 
+      IF(ldwp)WRITE(kunit,*) ' ~~~~~~~~~~~~~~~~'
+      IF(ldwp)WRITE(kunit,*) '               Data type           = ', &
+         &                TRIM( ADJUSTL( cl_data_type ) )
+      IF(ldwp)WRITE(kunit,*) '               Format version      = ',  &
+         &                TRIM( ADJUSTL( cl_format_version ) )
+      
+      IF ( ( ( INDEX( cl_data_type,"ENACT v1.0" ) == 1 ) .OR.   &
+         &   ( INDEX( cl_data_type,"ENACT v1.4" ) == 1 ) .OR.   &
+         &   ( INDEX( cl_data_type,"ENACT v1.5" ) == 1 ) .OR.   &
+         &   ( INDEX( cl_data_type,"ENSEMBLES EN3 v1" ) == 1 )      ) &
+         &   .AND.                                              &
+         &   ( INDEX( cl_format_version,"2.0"   ) == 1 ) ) THEN
+         IF(ldwp)WRITE(kunit,*)'               Valid input file'
+      ELSE
+         CALL fatal_error( 'Invalid input file', __LINE__ )
+      ENDIF
+
+      !---------------------------------------------------------------------
+      ! Read the number of observations and levels to allocate arrays
+      !---------------------------------------------------------------------
+
+      CALL chkerr( nf90_inq_dimid        ( i_file_id, 'N_PROF', i_obs_id ),         &
+         &         cl_name, __LINE__ )
+      CALL chkerr( nf90_inquire_dimension( i_file_id, i_obs_id, len = iobs ),     &
+         &         cl_name, __LINE__ )
+      CALL chkerr( nf90_inq_dimid        ( i_file_id, 'N_LEVELS', i_lev_id ),     &
+         &         cl_name, __LINE__ )
+      CALL chkerr( nf90_inquire_dimension( i_file_id, i_lev_id, len = ilev ), &
+         &         cl_name, __LINE__ )
+      IF(ldwp)WRITE(kunit,*) '               No. of data records = ', iobs
+      IF(ldwp)WRITE(kunit,*) '               No. of levels       = ', ilev
+      IF(ldwp)WRITE(kunit,*) 
+      IF (ilev > imaxlev) THEN
+         CALL fatal_error( 'Increase imaxlev in obs_prof_io.F90', __LINE__ )
+      ENDIF
+
+      !---------------------------------------------------------------------
+      ! Allocate arrays
+      !---------------------------------------------------------------------
+
+      CALL init_obfbdata( inpfile )
+      CALL alloc_obfbdata( inpfile, 2, iobs, ilev, 0, 1, ldgrid )
+      inpfile%cname(1) = 'POTM'
+      inpfile%cname(2) = 'PSAL'
+      inpfile%coblong(1) = 'Potential temperature'
+      inpfile%coblong(2) = 'Practical salinity'
+      inpfile%cobunit(1) = 'Degrees Celsius'
+      inpfile%cobunit(2) = 'PSU'
+      inpfile%cextname(1) = 'TEMP'
+      inpfile%cextlong(1) = 'Insitu temperature'
+      inpfile%cextunit(1) = 'Degrees Celsius'
+
+      !---------------------------------------------------------------------
+      ! Read the QC atributes
+      !---------------------------------------------------------------------
+
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_var_qc_t, i_qc_var_id ),                         &
+         &         cl_name, __LINE__ )        
+      istart2(1) = 1
+      icount2(2) = 1
+      icount2(1) = ilev
+      DO ji = 1, iobs 
+         istart2(2) = ji
+         CALL chkerr( nf90_get_var  ( i_file_id, i_qc_var_id, clqc,                                   &
+            &                         start = istart2, count = icount2),                              &
+            &         cl_name, __LINE__ )
+         DO jk = 1, ilev
+            inpfile%ivlqc(jk,ji,1) = IACHAR( clqc(jk:jk) ) - IACHAR( '0' )
+         END DO
+      END DO
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_var_qc_s, i_qc_var_id ),                         &
+         &         cl_name, __LINE__ )
+      DO ji = 1, iobs
+         istart2(2) = ji
+         CALL chkerr( nf90_get_var  ( i_file_id, i_qc_var_id, clqc,                                   &
+            &                         start = istart2, count = icount2),                              &
+            &         cl_name, __LINE__ )
+         DO jk = 1, ilev
+            inpfile%ivlqc(jk,ji,2) = IACHAR( clqc(jk:jk) ) - IACHAR( '0' )
+         END DO
+      END DO
+      ! No depth QC in files
+      DO ji = 1, iobs
+         DO jk = 1, ilev
+            inpfile%idqc(jk,ji)  = 1
+            inpfile%idqcf(:,jk,ji) = 0
+         END DO
+      END DO
+
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_var_prof_qc_t, i_qc_var_id ),                    &
+         &         cl_name,  __LINE__ )
+      DO ji = 1,iobs
+         istart1(1) = ji
+         CALL chkerr( nf90_get_var  ( i_file_id, i_qc_var_id, cqc,                                    &
+            &                         start = istart1, count = icount1),                              &
+            &         cl_name, __LINE__ )
+         inpfile%ivqc(ji,1) = IACHAR( cqc ) - IACHAR( '0' )
+      END DO
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_var_prof_qc_s, i_qc_var_id ),                    & 
+         &         cl_name, __LINE__ ) 
+      DO ji = 1,iobs
+         istart1(1) = ji
+         CALL chkerr( nf90_get_var  ( i_file_id, i_qc_var_id, cqc,                                    &
+            &                         start = istart1, count = icount1),                              &
+            &         cl_name, __LINE__ )
+         inpfile%ivqc(ji,2) = IACHAR( cqc ) - IACHAR( '0' )
+      END DO
+!!      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_var_juld_qc, i_qc_var_id ),                       &
+!!         &         cl_name, __LINE__ ) 
+!!      !DO ji = 1,iobs
+!!         istart1(1) = ji
+!!         CALL chkerr( nf90_get_var  ( i_file_id, i_qc_var_id, cqc,                                    &
+!!            &                         start = istart1, count = icount1),                              &
+!!            &         cl_name, __LINE__ )
+!!         inpfile%itqc(ji)    = IACHAR( cqc ) - IACHAR( '0' )
+!!         inpfile%itqcf(:,ji) = 0
+!!      END DO
+      ! Since the flags are not set in the ENACT files we reset them to 0
+      inpfile%itqc(:)    = 1      
+      inpfile%itqcf(:,:) = 0
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_var_pos_qc, i_qc_var_id ),                       &
+         &         cl_name, __LINE__ ) 
+      DO ji = 1,iobs
+         istart1(1) = ji
+         CALL chkerr( nf90_get_var  ( i_file_id, i_qc_var_id, cqc,                                    &
+            &                         start = istart1, count = icount1),                              &
+            &         cl_name, __LINE__ )
+         inpfile%ipqc(ji)    = IACHAR( cqc ) - IACHAR( '0' )
+         inpfile%ipqcf(:,ji) = 0
+      END DO
+      DO ji = 1,iobs
+         inpfile%ioqc(ji) = MIN( inpfile%ivqc(ji,1), inpfile%ivqc(ji,2) )
+      END DO
+      IF ( nf90_inq_varid( i_file_id, cl_fld_qc_flags_profiles, i_qc_flag_id ) == nf90_noerr ) THEN
+         ALLOCATE( &
+            & iqc1(iobs) &
+            & )
+         CALL chkerr( nf90_get_var  ( i_file_id, i_qc_flag_id, iqc1 ),                                &
+            &         cl_name, __LINE__ )
+         DO ji = 1,iobs
+            inpfile%ioqcf(1,ji)   = iqc1(ji)
+            inpfile%ivqcf(1,ji,:) = iqc1(ji)
+            inpfile%ioqcf(2,ji)   = 0
+            inpfile%ivqcf(2,ji,:) = 0
+         END DO
+         DEALLOCATE( &
+            & iqc1 &
+            & )
+      ELSE
+         IF(ldwp) WRITE(kunit,*)'No QC profile flags in file'
+         inpfile%ioqcf(:,:)   = 0
+         inpfile%ivqcf(:,:,:) = 0
+      ENDIF
+      IF ( nf90_inq_varid( i_file_id, cl_fld_qc_flags_levels, i_qc_flag_id ) == nf90_noerr ) THEN
+         ALLOCATE( &
+            & iqc2(ilev,iobs) &
+            & )
+         CALL chkerr( nf90_get_var  ( i_file_id, i_qc_flag_id, iqc2 ),                                &
+            &         cl_name, __LINE__ )
+         DO ji = 1,iobs
+            DO jk = 1,ilev
+               inpfile%ivlqcf(1,jk,ji,:) = iqc2(jk,ji)
+               inpfile%ivlqcf(2,jk,ji,:) = 0
+            END DO
+         END DO
+         DEALLOCATE( &
+            & iqc2 &
+            & )
+      ELSE
+         IF(ldwp) WRITE(kunit,*)'No QC level flags in file'
+         inpfile%ivlqcf(:,:,:,:) = 0
+      ENDIF
+
+      !---------------------------------------------------------------------
+      ! Read the time/position variables 
+      !---------------------------------------------------------------------
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_juld, i_juld_id ),                               &
+         &         cl_name, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_juld_id, inpfile%ptim ),                              &
+         &         cl_name, __LINE__ )
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_depth, i_depth_id ),                             &
+            &         cl_name, __LINE__ )         
+      CALL chkerr( nf90_get_var  ( i_file_id, i_depth_id, inpfile%pdep ),                             &
+         &         cl_name, __LINE__ )
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_phi, i_phi_id ),                                 &
+         &         cl_name, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_phi_id, inpfile%pphi ),                               &
+         &         cl_name, __LINE__ )
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_lam, i_lam_id ),                                 &
+         &         cl_name, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_lam_id, inpfile%plam ),                               &
+         &         cl_name, __LINE__ )
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_reference_date_time, i_reference_date_time_id ), &
+         &         cl_name, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_reference_date_time_id, inpfile%cdjuldref ),          &
+         &         cl_name, __LINE__ )
+      
+      !---------------------------------------------------------------------
+      ! Read the platform information
+      !---------------------------------------------------------------------
+
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_wmo_inst_type, i_wmo_inst_type_id ),             &
+         &         cl_name, __LINE__ )          
+      CALL chkerr( nf90_get_var  ( i_file_id, i_wmo_inst_type_id, inpfile%cdtyp ),                    &
+         &         cl_name, __LINE__ )
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_pl_num, i_pl_num_id ),                           &
+         &         cl_name, __LINE__ )         
+      CALL chkerr( nf90_get_var  ( i_file_id, i_pl_num_id, inpfile%cdwmo ),                           &
+         &         cl_name, __LINE__ )         
+
+      !---------------------------------------------------------------------
+      ! Read the variables
+      !---------------------------------------------------------------------
+
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_var_tp, i_var_id ),                              &
+         &         cl_name, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, inpfile%pob(:,:,1) ),                         &
+         &         cl_name, __LINE__ )
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_var_s, i_var_id ),                               &
+         &         cl_name, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, inpfile%pob(:,:,2) ),                         &
+         &         cl_name, __LINE__ )
+ 
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_var_ti, i_var_id ),                              &
+         &         cl_name, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, inpfile%pext(:,:,1) ),                        &
+         &         cl_name, __LINE__ )
+
+      !---------------------------------------------------------------------
+      ! Close file
+      !---------------------------------------------------------------------
+
+      CALL chkerr( nf90_close( i_file_id ),           cl_name, __LINE__ )
+
+      !---------------------------------------------------------------------
+      ! Set file indexes
+      !---------------------------------------------------------------------
+      DO ji = 1, inpfile%nobs
+         inpfile%kindex(ji) = ji
+      END DO
+
+   END SUBROUTINE read_enactfile
+
+   SUBROUTINE read_coriofile( cdfilename, inpfile, kunit, ldwp, ldgrid )
+      !!---------------------------------------------------------------------
+      !!
+      !!                     ** ROUTINE read_coriofile **
+      !!
+      !! ** Purpose : Read from file the profile CORIO observations.
+      !!
+      !! ** Method  : The data file is a NetCDF file. 
+      !!
+      !! ** Action  :
+      !!
+      !! History : 
+      !!          ! 09-01 (K. Mogensen) Original based on old versions
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      CHARACTER(LEN=*) :: cdfilename ! Input filename
+      TYPE(obfbdata)   :: inpfile    ! Output enactfile structure
+      INTEGER          :: kunit      ! Unit for output
+      LOGICAL          :: ldwp       ! Print info
+      LOGICAL          :: ldgrid     ! Save grid info in data structure
+      INTEGER  :: &
+         & iobs, &
+         & ilev
+      !! * Local declarations
+      INTEGER :: &
+         & i_file_id,                &
+         & i_obs_id,                 &
+         & i_lev_id,                 &
+         & i_phi_id,                 & 
+         & i_lam_id,                 &
+         & i_depth_id,               &
+         & i_pres_id,                &
+         & i_var_id,                 &
+         & i_pl_num_id,              &
+         & i_format_version_id,      &
+         & i_juld_id,                &
+         & i_data_type_id,           &
+         & i_wmo_inst_type_id,       &
+         & i_qc_var_id,              &
+         & i_dc_ref_id
+      CHARACTER(LEN=40) :: & 
+         & cl_fld_lam,                 &
+         & cl_fld_phi,                 &
+         & cl_fld_depth,               & 
+         & cl_fld_depth_qc,            & 
+         & cl_fld_pres,                & 
+         & cl_fld_pres_qc,             & 
+         & cl_fld_var_t,               &
+         & cl_fld_var_s,               &
+         & cl_fld_var_ti,              &
+         & cl_fld_var_pos_qc,          &
+         & cl_fld_var_qc_t,            &
+         & cl_fld_var_qc_s,            &
+         & cl_fld_var_prof_qc_t,       &
+         & cl_fld_var_prof_qc_s,       &
+         & cl_fld_dc_ref,              &
+         & cl_fld_juld,                &
+         & cl_fld_pl_num,              &
+         & cl_fld_wmo_inst_type
+      CHARACTER(LEN=14), PARAMETER :: &
+         & cl_name = 'read_coriofile'
+      CHARACTER(LEN=4 )            :: &
+         & cl_format_version = ''
+      INTEGER, DIMENSION(1) :: &
+         & istart1, icount1
+      INTEGER, DIMENSION(2) :: &
+         & istart2, icount2
+      CHARACTER(len=imaxlev) :: &
+         & clqc
+      CHARACTER(len=1) :: &
+         & cqc
+      CHARACTER(len=256) :: &
+         & cdjulref
+      INTEGER :: &
+         & ji, jk
+      INTEGER :: &
+         & iformat
+      LOGICAL :: & 
+         & lsal
+      REAL(fbdp), DIMENSION(:,:), ALLOCATABLE :: &
+         & zpres
+      INTEGER, DIMENSION(:,:), ALLOCATABLE :: &
+         & ipresqc
+      CHARACTER(len=256) :: &
+         & cerr
+      !-----------------------------------------------------------------------
+      ! Initialization
+      !-----------------------------------------------------------------------
+
+      icount1(1) = 1 
+      lsal = .TRUE.
+
+      !-----------------------------------------------------------------------
+      ! Open file
+      !-----------------------------------------------------------------------
+
+      CALL chkerr( nf90_open( TRIM( cdfilename ), nf90_nowrite, &
+            &      i_file_id ),           cl_name, __LINE__ )
+
+      !-----------------------------------------------------------------------
+      ! Check format and set variables accordingly
+      !-----------------------------------------------------------------------
+
+      IF ( ( nf90_inq_dimid( i_file_id, 'N_PROF', i_obs_id ) == nf90_noerr ) .AND. &
+         & ( nf90_inq_dimid( i_file_id, 'N_LEVELS', i_lev_id ) == nf90_noerr ) ) THEN
+         iformat = 1
+      ELSEIF ( ( nf90_inq_dimid( i_file_id, 'mN_PROF', i_obs_id ) == nf90_noerr ) .AND. &
+         & ( nf90_inq_dimid( i_file_id, 'mN_ZLEV', i_lev_id ) == nf90_noerr ) ) THEN
+         iformat = 2
+      ELSE
+         WRITE(cerr,'(2A)')'Invalid data format in ',cl_name
+         CALL fatal_error( cerr, __LINE__ )
+      ENDIF
+      IF ( iformat == 1 ) THEN
+         cl_fld_lam                 = 'LONGITUDE'
+         cl_fld_phi                 = 'LATITUDE'  
+         cl_fld_depth               = 'DEPH'
+         cl_fld_depth_qc            = 'DEPH_QC'
+         cl_fld_pres                = 'PRES'
+         cl_fld_pres_qc             = 'PRES_QC'
+         cl_fld_juld                = 'JULD'
+         cl_fld_wmo_inst_type       = 'WMO_INST_TYPE'
+         cl_fld_dc_ref              = 'DC_REFERENCE'
+         cl_fld_pl_num              = 'PLATFORM_NUMBER'
+         cl_fld_var_qc_t            = 'TEMP_QC'
+         cl_fld_var_prof_qc_t       = 'PROFILE_TEMP_QC'
+         cl_fld_var_t               = 'TEMP'
+         cl_fld_var_qc_s            = 'PSAL_QC'
+         cl_fld_var_prof_qc_s       = 'PROFILE_PSAL_QC'
+         cl_fld_var_s               = 'PSAL'
+         cl_fld_var_pos_qc          = 'POSITION_QC'
+      ELSEIF ( iformat==2 ) THEN
+         cl_fld_lam                 = 'LONGITUDE'
+         cl_fld_phi                 = 'LATITUDE'  
+         cl_fld_depth               = 'DEPH'
+         cl_fld_depth_qc            = 'QC_DEPH'
+         cl_fld_pres                = 'PRES'
+         cl_fld_pres_qc             = 'QC_PRES'
+         cl_fld_juld                = 'JULD'
+         cl_fld_wmo_inst_type       = 'INST_TYPE'
+         cl_fld_dc_ref              = 'REFERENCE'
+         cl_fld_pl_num              = 'PLATFORM_NUMBER'
+         cl_fld_var_qc_t            = 'QC_TEMP'
+         cl_fld_var_prof_qc_t       = 'Q_PROFILE_TEMP'
+         cl_fld_var_t               = 'TEMP'
+         cl_fld_var_qc_s            = 'QC_PSAL'
+         cl_fld_var_prof_qc_s       = 'Q_PROFILE_PSAL'
+         cl_fld_var_s               = 'PSAL'
+         cl_fld_var_pos_qc          = 'Q_POSITION'
+      ENDIF
+
+      !-----------------------------------------------------------------------
+      ! Read the heading of the file
+      !-----------------------------------------------------------------------
+
+      IF(ldwp)WRITE(kunit,*)
+      IF(ldwp)WRITE(kunit,*) ' read_coriofile :' 
+      IF(ldwp)WRITE(kunit,*) ' ~~~~~~~~~~~~~~~~'
+      IF(ldwp)WRITE(kunit,*) '               Format version      = ', iformat
+
+      !---------------------------------------------------------------------
+      ! Read the number of observations and levels to allocate arrays
+      !---------------------------------------------------------------------
+
+      CALL chkerr( nf90_inquire_dimension( i_file_id, i_obs_id, len = iobs ), &
+         &         cl_name, __LINE__ )
+      CALL chkerr( nf90_inquire_dimension( i_file_id, i_lev_id, len = ilev ), &
+         &         cl_name, __LINE__ )
+      IF(ldwp)WRITE(kunit,*) '               No. of data records = ', iobs
+      IF(ldwp)WRITE(kunit,*) '               No. of levels       = ', ilev
+      IF(ldwp)WRITE(kunit,*) 
+      IF (ilev > imaxlev) THEN
+         CALL fatal_error( 'Increase imaxlev in obs_prof_io.F90', __LINE__ )
+      ENDIF
+
+      !---------------------------------------------------------------------
+      ! Allocate arrays
+      !---------------------------------------------------------------------
+
+      CALL init_obfbdata( inpfile )
+      CALL alloc_obfbdata( inpfile, 2, iobs, ilev, 0, 1, ldgrid )
+      inpfile%cname(1) = 'POTM'
+      inpfile%cname(2) = 'PSAL'
+      inpfile%coblong(1) = 'Potential temperature'
+      inpfile%coblong(2) = 'Practical salinity'
+      inpfile%cobunit(1) = 'Degrees Celsius'
+      inpfile%cobunit(2) = 'PSU'
+      inpfile%cextname(1) = 'TEMP'
+      inpfile%cextlong(1) = 'Insitu temperature'
+      inpfile%cextunit(1) = 'Degrees Celsius'
+      ALLOCATE( &
+         & zpres(ilev,iobs),  &
+         & ipresqc(ilev,iobs) &
+         & )
+      !---------------------------------------------------------------------
+      ! Get julian data reference (iformat==2)
+      !---------------------------------------------------------------------
+
+      IF (iformat==2) THEN
+         CALL chkerr ( nf90_get_att( i_file_id, nf90_global, &
+            &                        "Reference_date_time", cdjulref ), &
+            &         cl_name, __LINE__ )
+         inpfile%cdjuldref = cdjulref(7:10)//cdjulref(4:5)// &
+            & cdjulref(1:2)//cdjulref(12:13)//cdjulref(15:16)//cdjulref(18:19)
+      ENDIF
+
+      !---------------------------------------------------------------------
+      ! Read the QC attributes
+      !---------------------------------------------------------------------
+
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_var_qc_t, i_qc_var_id ),                                &
+         &         cl_name, __LINE__ )
+      istart2(1) = 1
+      icount2(2) = 1
+      icount2(1) = ilev
+      DO ji = 1, iobs 
+         istart2(2) = ji
+         CALL chkerr( nf90_get_var  ( i_file_id, i_qc_var_id, clqc,                                   &
+            &                         start = istart2, count = icount2),                              &
+            &         cl_name, __LINE__ )
+         DO jk = 1, ilev
+            inpfile%ivlqc(jk,ji,1) = IACHAR( clqc(jk:jk) ) - IACHAR( '0' )
+         END DO
+      END DO
+      IF ( nf90_inq_varid( i_file_id, cl_fld_var_qc_s, i_qc_var_id ) == nf90_noerr ) THEN
+         DO ji = 1, iobs
+            istart2(2) = ji
+            CALL chkerr( nf90_get_var  ( i_file_id, i_qc_var_id, clqc,                                &
+               &                         start = istart2, count = icount2),                           &
+               &         cl_name, __LINE__ )
+            DO jk = 1, ilev
+               inpfile%ivlqc(jk,ji,2) = IACHAR( clqc(jk:jk) ) - IACHAR( '0' )
+            END DO
+         END DO
+      ELSE
+         inpfile%ivlqc(:,:,2) = 4
+         inpfile%pob(:,:,2) = fbrmdi
+         lsal = .FALSE.
+      ENDIF
+
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_var_prof_qc_t, i_qc_var_id ),                    &
+         &         cl_name,  __LINE__ )
+      DO ji = 1,iobs
+         istart1(1) = ji
+         CALL chkerr( nf90_get_var  ( i_file_id, i_qc_var_id, cqc,                                    &
+            &                         start = istart1, count = icount1),                              &
+            &         cl_name, __LINE__ )
+         inpfile%ivqc(ji,1) = IACHAR( cqc ) - IACHAR( '0' )
+      END DO
+      IF (lsal) THEN
+         CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_var_prof_qc_s, i_qc_var_id ),                 & 
+            &         cl_name, __LINE__ ) 
+         DO ji = 1,iobs
+            istart1(1) = ji
+            CALL chkerr( nf90_get_var  ( i_file_id, i_qc_var_id, cqc,                                 &
+               &                         start = istart1, count = icount1),                           &
+               &         cl_name, __LINE__ )
+            inpfile%ivqc(ji,2) = IACHAR( cqc ) - IACHAR( '0' )
+         END DO
+      ELSE
+         inpfile%ivqc(:,2) = 4
+      ENDIF
+      DO ji = 1,iobs
+         inpfile%ioqc(ji) = MIN( inpfile%ivqc(ji,1), inpfile%ivqc(ji,2) )
+      END DO
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_var_pos_qc, i_qc_var_id ),                       &
+         &         cl_name, __LINE__ ) 
+      DO ji = 1, iobs
+         istart1(1) = ji
+         CALL chkerr( nf90_get_var  ( i_file_id, i_qc_var_id, cqc,                                    &
+            &                         start = istart1, count = icount1),                              &
+            &         cl_name, __LINE__ )
+         inpfile%ipqc(ji)  = IACHAR( cqc ) - IACHAR( '0' )
+      END DO
+      
+      !---------------------------------------------------------------------
+      ! Read the time/position variables 
+      !---------------------------------------------------------------------
+
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_juld, i_juld_id ),                               &
+         &         cl_name, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_juld_id, inpfile%ptim ),                              &
+         &         cl_name, __LINE__ )
+      IF (iformat==1) THEN
+         CALL chkerr ( nf90_get_att( i_file_id, i_juld_id, &
+            &                        "units", cdjulref ), &
+            &         cl_name, __LINE__ )
+         inpfile%cdjuldref = cdjulref(12:15)//cdjulref(17:18)// &
+            & cdjulref(20:21)//cdjulref(23:24)//cdjulref(26:27)//&
+            & cdjulref(29:30)
+      ENDIF
+      
+      IF ( nf90_inq_varid( i_file_id, cl_fld_depth, i_depth_id ) == nf90_noerr ) THEN
+         CALL chkerr( nf90_get_var  ( i_file_id, i_depth_id, inpfile%pdep ),                          &
+            &         cl_name, __LINE__ )
+         CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_depth_qc, i_qc_var_id ),                      & 
+            &         cl_name, __LINE__ ) 
+         DO ji = 1, iobs
+            istart2(2) = ji
+            CALL chkerr( nf90_get_var  ( i_file_id, i_qc_var_id, clqc,                                &
+               &                         start = istart2, count = icount2),                           &
+               &         cl_name, __LINE__ )
+            DO jk = 1, ilev
+               inpfile%idqc(jk,ji) = IACHAR( clqc(jk:jk) ) - IACHAR( '0' )
+            END DO
+         END DO
+      ELSE
+         inpfile%pdep(:,:) = fbrmdi
+         inpfile%idqc(:,:) = 4
+      ENDIF
+
+      IF ( nf90_inq_varid( i_file_id, cl_fld_pres, i_pres_id ) == nf90_noerr ) THEN
+         CALL chkerr( nf90_get_var  ( i_file_id, i_pres_id, zpres ),                                  &
+            &         cl_name, __LINE__ )
+         CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_pres_qc, i_qc_var_id ),                       & 
+            &         cl_name, __LINE__ ) 
+         DO ji = 1, iobs
+            istart2(2) = ji
+            CALL chkerr( nf90_get_var  ( i_file_id, i_qc_var_id, clqc,                                &
+               &                         start = istart2, count = icount2),                           &
+               &         cl_name, __LINE__ )
+            DO jk = 1, ilev
+               ipresqc(jk,ji) = IACHAR( clqc(jk:jk) ) - IACHAR( '0' )
+            END DO
+         END DO
+      ELSE
+         zpres(:,:) = fbrmdi
+         ipresqc(:,:) = 4
+      ENDIF
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_phi, i_phi_id ),                                 &
+         &         cl_name, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_phi_id, inpfile%pphi ),                               &
+         &         cl_name, __LINE__ )
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_lam, i_lam_id ),                                 &
+         &         cl_name, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_lam_id, inpfile%plam ),                               &
+         &         cl_name, __LINE__ )
+      
+      !---------------------------------------------------------------------
+      ! Read the platform information
+      !---------------------------------------------------------------------
+
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_wmo_inst_type, i_wmo_inst_type_id ),             &
+         &         cl_name, __LINE__ )          
+      CALL chkerr( nf90_get_var  ( i_file_id, i_wmo_inst_type_id, inpfile%cdtyp ),                    &
+         &         cl_name, __LINE__ )
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_pl_num, i_pl_num_id ),                           &
+         &         cl_name, __LINE__ )         
+      CALL chkerr( nf90_get_var  ( i_file_id, i_pl_num_id, inpfile%cdwmo ),                           &
+         &         cl_name, __LINE__ )         
+
+      
+      !---------------------------------------------------------------------
+      ! Read the variables
+      !---------------------------------------------------------------------
+
+      CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_var_t, i_var_id ),                               &
+         &         cl_name, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, inpfile%pext(:,:,1) ),                        &
+         &         cl_name, __LINE__ )
+
+      IF (lsal) THEN      
+         CALL chkerr( nf90_inq_varid( i_file_id, cl_fld_var_s, i_var_id ),                            &
+            &         cl_name, __LINE__ )
+         CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, inpfile%pob(:,:,2) ),                      &
+            &         cl_name, __LINE__ )
+      ENDIF
+
+      !---------------------------------------------------------------------
+      ! Close file
+      !---------------------------------------------------------------------
+         
+      CALL chkerr( nf90_close( i_file_id ),           cl_name, __LINE__ )
+
+      !---------------------------------------------------------------------
+      ! Set file indexes
+      !---------------------------------------------------------------------
+      DO ji = 1, inpfile%nobs
+         inpfile%kindex(ji) = ji
+      END DO
+      
+      !---------------------------------------------------------------------
+      !  Coriolis data conversion from insitu to potential temperature
+      !---------------------------------------------------------------------
+      !---------------------------------------------------------------------
+      !  Convert pressure to depth if depth not present
+      !---------------------------------------------------------------------
+      DO ji = 1, inpfile%nobs
+         IF ( inpfile%pphi(ji) < 9999.0 ) THEN
+            DO jk = 1, inpfile%nlev
+               IF ( inpfile%pdep(jk,ji) >= 9999.0 ) THEN
+                  IF ( zpres(jk,ji) < 9999.0 ) THEN
+                     inpfile%pdep(jk,ji) = &
+                        & p_to_dep( REAL(zpres(jk,ji),wp), REAL(inpfile%pphi(ji),wp) )
+                     inpfile%idqc(jk,ji) = ipresqc(jk,ji)
+                  ENDIF
+               ENDIF
+            END DO
+         ENDIF
+      END DO
+      
+      !---------------------------------------------------------------------
+      !  Convert depth to pressure if pressure not present
+      !---------------------------------------------------------------------
+      DO ji = 1, inpfile%nobs
+         IF ( inpfile%pphi(ji) < 9999.0 ) THEN
+            DO jk = 1, inpfile%nlev
+               IF ( zpres(jk,ji) >= 9999.0 ) THEN
+                  IF ( inpfile%pdep(jk,ji) < 9999.0 ) THEN
+                     zpres(jk,ji) = dep_to_p( REAL(inpfile%pdep(jk,ji),wp), &
+                        &                            REAL(inpfile%pphi(ji),wp) )
+                     ipresqc(jk,ji) = inpfile%idqc(jk,ji)
+                  ENDIF
+               ENDIF
+            END DO
+         ENDIF
+      END DO
+      
+      !---------------------------------------------------------------------
+      !  Convert insitu temperature to potential temperature if
+      !  salinity, insitu temperature and pressure are present
+      !---------------------------------------------------------------------
+      DO ji = 1, inpfile%nobs
+         DO jk = 1, inpfile%nlev
+            IF (( inpfile%pob(jk,ji,2) < 9999.0 ) .AND. &
+               &( inpfile%pext(jk,ji,1) < 9999.0 ) .AND. &
+               &( zpres(jk,ji) < 9999.0 ) ) THEN
+               inpfile%pob(jk,ji,1) = potemp( REAL(inpfile%pob(jk,ji,2), wp),  &
+                  &                           REAL(inpfile%pext(jk,ji,1), wp), &
+                  &                           REAL(zpres(jk,ji),wp),  &
+                  &                           0.0_wp )
+            ELSE
+               inpfile%pob(jk,ji,1) = fbrmdi
+            ENDIF
+         END DO
+      END DO
+
+      !---------------------------------------------------------------------
+      !  Initialize flags since they are not in the CORIOLIS input files
+      !---------------------------------------------------------------------
+
+      inpfile%ioqcf(:,:)      = 0
+      inpfile%ipqcf(:,:)      = 0
+      inpfile%itqcf(:,:)      = 0
+      inpfile%idqcf(:,:,:)    = 0
+      inpfile%ivqcf(:,:,:)    = 0
+      inpfile%ivlqcf(:,:,:,:) = 0
+
+   END SUBROUTINE read_coriofile
Index: /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obssla_io.h90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obssla_io.h90	(revision 6071)
+++ /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obssla_io.h90	(revision 6071)
@@ -0,0 +1,378 @@
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: obssla_io.h90 2287 2010-10-18 07:53:52Z smasson $
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+   SUBROUTINE read_avisofile( cdfilename, inpfile, kunit, ldwp, ldgrid )
+      !!---------------------------------------------------------------------
+      !!
+      !!                     ** ROUTINE read_avisofile **
+      !!
+      !! ** Purpose : Read from file the AVISO SLA observations.
+      !!
+      !! ** Method  : The data file is a NetCDF file. 
+      !!
+      !! ** Action  :
+      !!
+      !! References : http://www.aviso.oceanobs.com
+      !!
+      !! History : 
+      !!          ! 09-01 (K. Mogensen) Original based on old versions
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      CHARACTER(LEN=*) :: cdfilename ! Input filename
+      TYPE(obfbdata)   :: inpfile    ! Output obfbdata structure
+      INTEGER          :: kunit      ! Unit for output
+      LOGICAL          :: ldwp       ! Print info
+      LOGICAL          :: ldgrid     ! Save grid info in data structure
+      !! * Local declarations
+      CHARACTER(LEN=14),PARAMETER :: cpname = 'read_avisofile'
+      INTEGER :: i_file_id     ! netcdf IDS
+      INTEGER :: i_tracks_id
+      INTEGER :: i_cycles_id
+      INTEGER :: i_data_id
+      INTEGER :: i_var_id
+      INTEGER, PARAMETER :: imaxdim = 2    ! Assumed maximum for no. dims. in file
+      INTEGER, DIMENSION(2) :: idims       ! Dimensions in file
+      INTEGER :: iilen         ! Length of netCDF attributes
+      INTEGER :: itype         ! Typeof netCDF attributes
+      REAL(fbdp) :: zsca       ! Scale factor
+      REAL(fbdp) :: zfill      ! Fill value
+      CHARACTER(len=3) ::  cmission      ! Mission global attribute
+      INTEGER :: itracks       ! Maximum number of passes in file
+      INTEGER :: icycles       ! Maximum number of cycles for each pass
+      INTEGER :: idata         ! Number of data per parameter in current file
+      REAL(fbdp) :: zdeltat    ! Time gap getween two measurements in seconds
+      INTEGER, DIMENSION(:), POINTER :: & 
+         & iptracks,     &  ! List of passes contained in current file
+         & ipnbpoints,   &  ! Number of points per pass
+         & ipdataindexes    ! Index of data in theoretical profile
+      INTEGER, DIMENSION(:,:), POINTER :: & 
+         & ipcycles         ! List of cycles per pass
+      REAL(fbdp), DIMENSION(:), POINTER :: &
+         & zphi,   &        ! Latitudes
+         & zlam             ! Longitudes
+      REAL(fbdp), DIMENSION(:,:), POINTER :: &
+         & zbegindates      ! Date of point with index 0
+      REAL(fbdp) :: zbeginmiss    ! Missing data for dates
+      REAL(fbsp), DIMENSION(:,:), POINTER :: &
+         & zsla             ! SLA data
+      REAL(fbdp), DIMENSION(:), POINTER :: &
+         & zjuld            ! Julian date
+      LOGICAL, DIMENSION(:), POINTER :: &
+         & llskip           ! Skip observation
+      CHARACTER(len=14) :: cdjuldref    ! Julian data reference
+      INTEGER :: imission   ! Mission number converted from Mission global 
+                            ! netCDF atttribute.
+      CHARACTER(len=255) :: ctmp
+      INTEGER :: iobs
+      INTEGER :: jl
+      INTEGER :: jm
+      INTEGER :: jj
+      INTEGER :: ji
+      INTEGER :: jk
+      INTEGER :: jobs
+      INTEGER :: jcycle
+
+      ! Open the file
+
+      CALL chkerr( nf90_open( TRIM( cdfilename ), nf90_nowrite, i_file_id ), &
+         &         cpname, __LINE__ )
+
+      ! Get the netCDF dimensions
+      
+      CALL chkerr( nf90_inq_dimid( i_file_id, 'Tracks', i_tracks_id ),  &
+         &         cpname, __LINE__ )
+      CALL chkerr( nf90_inquire_dimension( i_file_id, i_tracks_id, &
+         &                                 len = itracks ),  &
+         &         cpname, __LINE__ )
+      
+      CALL chkerr( nf90_inq_dimid( i_file_id, 'Cycles', i_cycles_id ),  &
+         &         cpname, __LINE__ )
+      CALL chkerr( nf90_inquire_dimension( i_file_id, i_cycles_id, &
+         &                                 len = icycles ),  &
+         &      cpname, __LINE__ )
+      
+      CALL chkerr( nf90_inq_dimid( i_file_id, 'Data', i_data_id ),  &
+         &         cpname, __LINE__ )
+      CALL chkerr( nf90_inquire_dimension( i_file_id, i_data_id, &
+         &                                 len = idata ),  &
+         &         cpname, __LINE__ )
+
+      ! Allocate memory for input data
+
+      ALLOCATE( &
+         & iptracks     ( itracks          ), &  
+         & ipnbpoints   ( itracks          ), &  
+         & ipcycles     ( icycles, itracks ), &
+         & zphi         ( idata            ), &   
+         & zlam         ( idata            ), &  
+         & zbegindates  ( icycles, itracks ), &
+         & ipdataindexes( idata            ), & 
+         & zsla         ( icycles, idata   ), &
+         & zjuld        ( idata*icycles    ), &
+         & llskip       ( idata*icycles    )  &
+         & )
+
+      ! Get time gap getween two measurements in seconds
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, 'DeltaT', i_var_id ), & 
+         &         cpname, __LINE__ )
+      CALL chkdim( i_file_id, i_var_id, 0, idims, cpname, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, zdeltat), &
+         &         cpname, __LINE__ )
+      zsca = 1.0
+      IF (nf90_inquire_attribute( i_file_id, i_var_id, 'scale_factor') &
+         &                      == nf90_noerr) THEN
+         CALL chkerr( nf90_get_att( i_file_id, i_var_id, &
+            &                     'scale_factor',zsca), cpname,  __LINE__)
+      ENDIF
+      zdeltat = zsca * zdeltat
+         
+      ! Get List of passes contained in current file      
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, 'Tracks', i_var_id ), & 
+         &         cpname, __LINE__ )
+      idims(1) = itracks
+      CALL chkdim( i_file_id, i_var_id, 1, idims, cpname, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, iptracks), &
+         &         cpname, __LINE__ )
+      
+      ! Get Number of points per pass
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, 'NbPoints', i_var_id ), & 
+         &         cpname, __LINE__ )
+      idims(1) = itracks
+      CALL chkdim( i_file_id, i_var_id, 1, idims, cpname, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, ipnbpoints),&
+         &         cpname, __LINE__ )
+      
+      ! Get list of cycles per pass
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, 'Cycles', i_var_id ), & 
+         &         cpname, __LINE__ )
+      idims(1) = icycles
+      idims(2) = itracks
+      CALL chkdim( i_file_id, i_var_id, 2, idims, cpname, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, ipcycles),   &
+         &         cpname, __LINE__ )
+
+      ! Get longitudes
+
+      CALL chkerr( nf90_inq_varid( i_file_id, 'Longitudes', i_var_id ), & 
+         &         cpname, __LINE__ )
+      idims(1) = idata
+      CALL chkdim( i_file_id, i_var_id, 1, idims, cpname, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, zlam),   &
+         &         cpname, __LINE__ )
+      zsca = 1.0
+      IF (nf90_inquire_attribute( i_file_id, i_var_id, 'scale_factor') &
+         &                         == nf90_noerr) THEN
+         CALL chkerr( nf90_get_att( i_file_id, i_var_id, &
+            &                     'scale_factor',zsca), cpname,  __LINE__)
+      ENDIF
+      zlam(:) = zsca * zlam(:)
+      
+      ! Get latitudes
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, 'Latitudes', i_var_id ), & 
+         &         cpname, __LINE__ )
+      idims(1) = idata
+      CALL chkdim( i_file_id, i_var_id, 1, idims, cpname, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, zphi),   &
+         &         cpname, __LINE__ )
+      zsca = 1.0
+      IF (nf90_inquire_attribute( i_file_id, i_var_id, 'scale_factor') &
+         &                       == nf90_noerr) THEN
+         CALL chkerr( nf90_get_att( i_file_id, i_var_id, &
+            &                     'scale_factor',zsca), cpname,  __LINE__)
+      ENDIF
+      zphi(:) = zsca * zphi(:)
+      
+      ! Get date of point with index 0
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, 'BeginDates', i_var_id ), & 
+         &        cpname, __LINE__ )
+      idims(1) = icycles
+      idims(2) = itracks
+      CALL chkdim( i_file_id, i_var_id, 2, idims, cpname, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, zbegindates), &
+         &         cpname, __LINE__ )
+      CALL chkerr( nf90_inquire_attribute( i_file_id, i_var_id, &
+         &                                 'units', len = iilen,  &
+         &                                  xtype = itype), cpname, __LINE__ )
+      IF (( itype /= NF90_CHAR ) .OR. ( iilen > 255 )) THEN
+         CALL fatal_error('Error decoding BeginDates unit', __LINE__ )
+      ENDIF
+      CALL chkerr( nf90_get_att( i_file_id, i_var_id, 'units', &
+         &                       ctmp ), cpname, __LINE__ )
+      jl=1
+      DO jm = 1, LEN(TRIM(ctmp))
+         IF ((ctmp(jm:jm)>='0').AND.(ctmp(jm:jm)<='9')) THEN
+            cdjuldref(jl:jl) = ctmp(jm:jm)
+            jl=jl+1
+         ENDIF
+         IF (jl>14) EXIT
+      END DO
+      CALL chkerr( nf90_inquire_attribute( i_file_id, i_var_id, '_FillValue', &
+         &                                  xtype = itype), cpname, __LINE__ )
+      IF ( itype /= NF90_DOUBLE ) THEN
+         CALL fatal_error('Error decoding BeginDates missing data', __LINE__ )
+      ENDIF
+      CALL chkerr( nf90_get_att( i_file_id, i_var_id, '_FillValue', &
+         &                       zbeginmiss ), cpname, __LINE__ )
+
+      ! Get indices of data in theoretical profile
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, 'DataIndexes', i_var_id ), & 
+         &         cpname, __LINE__ )
+      idims(1) = idata
+      CALL chkdim( i_file_id, i_var_id, 1, idims, cpname, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, ipdataindexes),   &
+         &       cpname, __LINE__ )
+      
+      ! Get SLA data
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, 'SLA', i_var_id ), & 
+         &         cpname, __LINE__ )
+      idims(1) = icycles
+      idims(2) = idata
+      CALL chkdim( i_file_id, i_var_id, 2, idims, cpname, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, zsla),   &
+         &         cpname, __LINE__ )
+      zsca = 1.0
+      IF (nf90_inquire_attribute( i_file_id, i_var_id, 'scale_factor') &
+         &                        == nf90_noerr) THEN
+         CALL chkerr( nf90_get_att( i_file_id, i_var_id, &
+            &                       'scale_factor',zsca), cpname, __LINE__ )
+      ENDIF
+      zfill = 0.0
+      IF (nf90_inquire_attribute( i_file_id, i_var_id, '_FillValue') &
+         &                      == nf90_noerr) THEN
+         CALL chkerr( nf90_get_att( i_file_id, i_var_id, &
+            &                     '_FillValue',zfill), cpname,  __LINE__ )
+      ENDIF
+      WHERE(zsla(:,:) /=  zfill)
+         zsla(:,:) = zsca * zsla(:,:)
+      ELSEWHERE
+         zsla(:,:) = fbrmdi
+      END WHERE
+      
+      ! Get the global Mission netCDF attribute
+      
+      cmission='   '
+      CALL chkerr( nf90_inquire_attribute( i_file_id, nf90_global, &
+         &                                 'Mission', len = iilen,  &
+         &                                 xtype = itype), cpname, __LINE__ )
+      IF (( itype /= NF90_CHAR ) .OR. ( iilen > 3 )) THEN
+         CALL fatal_error('Error decoding Mission global attribute', __LINE__ )
+      ENDIF
+      CALL chkerr( nf90_get_att( i_file_id, nf90_global, 'Mission', &
+         &                       cmission ), cpname, __LINE__ )
+      
+      ! Convert it to an integer
+      imission = 0
+      DO jm = 1, imaxmissions
+         IF (cmission==cmissions(jm)) THEN
+            imission = jm
+            EXIT
+         ENDIF
+      END DO
+      
+      ! Close the file
+      
+      CALL chkerr( nf90_close( i_file_id ), cpname, __LINE__ )
+
+      ! Compute Julian dates for all observations
+      
+      jm = 0
+      jl = 0
+      DO jj = 1, itracks
+         DO ji = 1, ipnbpoints(jj)
+            jl = jl + 1
+            DO jk = 1, icycles
+               jm = jm + 1
+               IF (zbegindates(jk,jj)==zbeginmiss) THEN
+                  llskip(jm) = .TRUE.
+                  zjuld(jm)  = fbrmdi
+               ELSE
+                  llskip(jm) = .FALSE.
+                  zjuld(jm)  = zbegindates(jk,jj) + &
+                     &         (ipdataindexes(jl) * zdeltat / 86400._wp )
+               ENDIF
+            END DO
+         END DO
+      END DO
+
+      ! Get rid of missing data
+
+      jm = 0
+      DO jobs = 1, idata
+         DO jcycle = 1, icycles
+            jm = jm + 1
+            IF (zsla(jcycle,jobs) == fbrmdi) llskip(jm) = .TRUE.
+         END DO
+      END DO
+      
+      ! Allocate obfbdata
+
+      iobs = COUNT( .NOT.llskip(:) )
+      CALL init_obfbdata( inpfile )
+      CALL alloc_obfbdata( inpfile, 1, iobs, 1, 0, 0, ldgrid )
+      inpfile%cname(1) = 'SLA'
+
+      ! Fill the obfbdata structure from input data
+
+      inpfile%cdjuldref = cdjuldref
+      iobs = 0
+      jm = 0
+      DO jobs = 1, idata
+         DO jcycle = 1, icycles
+            jm = jm + 1
+            IF (llskip(jm)) CYCLE
+            iobs = iobs + 1
+            ! Characters
+            WRITE(inpfile%cdwmo(iobs),'(A3,A5)') cmissions(imission), '     '
+            WRITE(inpfile%cdtyp(iobs),'(I4)') imission
+            ! Real values
+            inpfile%plam(iobs)         = zlam(jobs)
+            inpfile%pphi(iobs)         = zphi(jobs)
+            inpfile%pob(1,iobs,1)      = zsla(jcycle,jobs)
+            inpfile%ptim(iobs)         = zjuld(jm)
+            inpfile%pdep(1,iobs)       = 0.0
+            ! Integers
+            inpfile%kindex(iobs)       = iobs
+            inpfile%ioqc(iobs)      = 1
+            inpfile%ivqc(iobs,1)    = 1
+            inpfile%ivlqc(1,iobs,1) = 1
+            inpfile%ipqc(iobs)         = 0 
+            inpfile%ipqcf(:,iobs)      = 0
+            inpfile%itqc(iobs)         = 0
+            inpfile%itqcf(:,iobs)      = 0
+            inpfile%ivqcf(:,iobs,1)    = 0
+            inpfile%ioqcf(:,iobs)      = 0
+            inpfile%idqc(1,iobs)       = 0
+            inpfile%idqcf(1,1,iobs)    = 0
+            inpfile%ivlqcf(:,1,iobs,1) = 0
+         END DO
+      END DO
+
+
+      ! Deallocate memory for input data
+
+      DEALLOCATE( &
+         & iptracks,      &
+         & ipnbpoints,    &
+         & ipcycles,      &
+         & zphi,          &
+         & zlam,          &  
+         & zbegindates,   &
+         & ipdataindexes, &
+         & zsla,          &
+         & zjuld,         &
+         & llskip         &
+         & )
+
+   END SUBROUTINE read_avisofile
+
Index: /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obssla_types.h90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obssla_types.h90	(revision 6071)
+++ /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/OBSTOOLS/src/obssla_types.h90	(revision 6071)
@@ -0,0 +1,9 @@
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: obssla_types.h90 2287 2010-10-18 07:53:52Z smasson $
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+   INTEGER, PARAMETER :: imaxmissions=8
+   CHARACTER(len=3) :: cmissions(0:imaxmissions) = &
+      & (/ 'XXX', 'E1 ', 'E2 ', 'TP ', 'TPM', 'G2 ', 'J1 ', 'EN ', 'J2 ' /)
Index: /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/REBUILD_NEMO/icb_combrest.py
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/REBUILD_NEMO/icb_combrest.py	(revision 6071)
+++ /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/REBUILD_NEMO/icb_combrest.py	(revision 6071)
@@ -0,0 +1,255 @@
+from netCDF4 import Dataset
+from argparse import ArgumentParser
+import numpy as np
+import sys
+
+#
+# Basic iceberg trajectory restart post-processing python script.
+# This script collects iceberg information from the distributed restarts written
+# out by each processing region and writes the information into a global restart file.
+# The global restart file must already exist and contain the collated 2D spatial fields
+# as prepared by utilities such as rebuild_nemo. This python script simply adds the
+# iceberg position and state data that is held using the unlimited dimension 'n' which
+# has been ignored by rebuild_nemo. Each processing region that contains icebergs will
+# (probably) have a different number of icebergs (indicated by differing values for the
+# current size of the unlimited dimension). This script collects all icebergs into a 
+# single unordered list.
+#
+
+parser = ArgumentParser(description='produce a global trajectory restart file from distributed output\
+                                     files, e.g. \n python ./icb_pp.py \
+                                     -f  icebergs_00692992_restart_ -n 480 -o icebergs_restart.nc [-O]')
+
+parser.add_argument('-f',dest='froot',help='fileroot_of_distrbuted_data; root name of \
+                     distributed trajectory restart file (usually completed with XXXX.nc, where \
+                     XXXX is the 4 digit processor number)', 
+                     default='icebergs_00692992_restart_')
+
+parser.add_argument('-n',dest='fnum',help='number of distributed files to process', 
+                     type=int, default=None)
+
+parser.add_argument('-o',dest='fout',help='global_iceberg_restart; file name to append the \
+                     global iceberg restart data to.', default='icebergs_restart.nc')
+
+parser.add_argument('-O',dest='fcre',help='Create the output file from scratch rather than \
+                     append to an existing file.', \
+                     action='store_true', default=False)
+
+args = parser.parse_args()
+
+default_used = 0
+if args.froot is None:
+    pathstart = 'icebergs_00692992_restart_'
+    default_used = 1
+else:
+    pathstart = args.froot
+
+if args.fnum is None:
+    procnum = 0
+    default_used = 1
+else:
+    procnum = args.fnum
+
+if args.fout is None:
+    pathout = 'icebergs_restart.nc'
+    default_used = 1
+else:
+    pathout = args.fout
+
+if default_used == 1:
+   print('At least one default value will be used; command executing is:')
+   print('icb_combrest.py -f ',pathstart,' -n ',procnum,' -o ',pathout)
+
+if procnum < 1:
+   print('Need some files to collate! procnum = ',procnum)
+   sys.exit()
+
+icu = []
+times = []
+ntraj = 0
+nk = 0
+#
+# Loop through all distributed datasets to obtain the total number 
+# of icebergs to transfer
+#
+for n in range(procnum):
+ nn = '%4.4d' % n
+ try:
+   fw = Dataset(pathstart+nn+'.nc')
+ except:
+   print 'Error: unable to open input file: ' + pathstart+nn+'.nc'
+   sys.exit()
+ for d in fw.dimensions :
+  if d == 'n' :
+   if len(fw.dimensions['n']) > 0:
+#    print 'icebergs found in: ' + pathstart+nn+'.nc' 
+     if len(fw.dimensions['k']) > nk :
+       nk = len(fw.dimensions['k'])
+     ntraj = ntraj + len(fw.dimensions['n'])
+ fw.close()
+#
+print(ntraj, ' icebergs found across all datasets')
+#
+# Declare 2-D arrays to receive the data from all files
+#
+lons = np.zeros(ntraj)
+lats = np.zeros(ntraj)
+xis  = np.zeros(ntraj)
+yjs  = np.zeros(ntraj)
+uvs  = np.zeros(ntraj)
+vvs  = np.zeros(ntraj)
+mas  = np.zeros(ntraj)
+ths  = np.zeros(ntraj)
+wis  = np.zeros(ntraj)
+les  = np.zeros(ntraj)
+dys  = np.zeros(ntraj)
+mss  = np.zeros(ntraj)
+msb  = np.zeros(ntraj)
+hds  = np.zeros(ntraj)
+yrs  = np.zeros(ntraj      , dtype=int)
+num  = np.zeros((ntraj, nk), dtype=int)
+#
+# loop through distributed datasets again, this time
+# collecting all trajectory data
+#
+nt = 0
+for n in range(procnum):
+ nn = '%4.4d' % n
+ fw = Dataset(pathstart+nn+'.nc')
+ for d in fw.dimensions :
+  if d == 'n' :
+# Note many distributed datafiles will contain no iceberg data
+# so skip quickly over these
+    m  = len(fw.dimensions['n'])
+    if m > 0:
+#     print pathstart+nn+'.nc'
+      lons[nt:nt+m]  = fw.variables['lon'][:]
+      lats[nt:nt+m]  = fw.variables['lat'][:]
+      xis[nt:nt+m]   = fw.variables['xi'][:]
+      yjs[nt:nt+m]   = fw.variables['yj'][:]
+      uvs[nt:nt+m]   = fw.variables['uvel'][:]
+      vvs[nt:nt+m]   = fw.variables['vvel'][:]
+      mas[nt:nt+m]   = fw.variables['mass'][:]
+      ths[nt:nt+m]   = fw.variables['thickness'][:]
+      wis[nt:nt+m]   = fw.variables['width'][:]
+      les[nt:nt+m]   = fw.variables['length'][:]
+      dys[nt:nt+m]   = fw.variables['day'][:]
+      mss[nt:nt+m]   = fw.variables['mass_scaling'][:]
+      msb[nt:nt+m]   = fw.variables['mass_of_bits'][:]
+      hds[nt:nt+m]   = fw.variables['heat_density'][:]
+      yrs[nt:nt+m]   = fw.variables['year'][:]
+      num[nt:nt+m,:] = fw.variables['number'][:,:]
+      nt = nt + m
+ fw.close()
+
+# Finally create the output file and write out the collated sets
+#
+if args.fcre :
+  try:
+    fo = Dataset(pathout, 'w', format='NETCDF4')
+  except:
+    print 'Error accessing output file: ' + pathout
+    print 'Check it is a writable location.'
+    sys.exit()
+else :
+  try:
+    fo = Dataset(pathout, 'r+', format='NETCDF4')
+  except:
+    print 'Error accessing output file: ' + pathout
+    print 'Check it exists and is writable.'
+    print 'Or run adding the -O option to create an output file which will'
+    print 'contain the iceberg state data only.'
+    sys.exit()
+#
+add_k = 1
+for d in fo.dimensions :
+  if d == 'n' :
+    print 'Error: dimension n already exists in output file'
+    sys.exit()
+  if d == 'k' :
+    add_k = 0
+onn  = fo.createDimension('n', None)
+if add_k > 0 :
+ onk = fo.createDimension('k', nk)
+olon = fo.createVariable('lon', 'f8',('n'))
+olat = fo.createVariable('lat', 'f8',('n'))
+oxis = fo.createVariable('xi', 'f8',('n'))
+oyjs = fo.createVariable('yj', 'f8',('n'))
+ouvs = fo.createVariable('uvel', 'f8',('n'))
+ovvs = fo.createVariable('vvel', 'f8',('n'))
+omas = fo.createVariable('mass', 'f8',('n'))
+oths = fo.createVariable('thickness', 'f8',('n'))
+owis = fo.createVariable('width', 'f8',('n'))
+oles = fo.createVariable('length', 'f8',('n'))
+odys = fo.createVariable('day', 'f8',('n'))
+omss = fo.createVariable('mass_scaling', 'f8',('n'))
+omsb = fo.createVariable('mass_of_bits', 'f8',('n'))
+ohds = fo.createVariable('heat_density', 'f8',('n'))
+oyrs = fo.createVariable('year', 'i4',('n'))
+onum = fo.createVariable('number', 'i4',('n','k'))
+#
+olon[:] = lons
+olon.long_name = "longitude"
+olon.units = "degrees_E" 
+#
+olat[:] = lats
+olat.long_name = "latitude"
+olat.units = "degrees_N"
+#
+oxis[:] = xis
+oxis.long_name = "x grid box position"
+oxis.units = "fractional"
+#
+oyjs[:] = yjs
+oyjs.long_name = "y grid box position"
+oyjs.units = "fractional"
+#
+ouvs[:] = uvs
+ouvs.long_name = "zonal velocity"
+ouvs.units = "m/s"
+#
+ovvs[:] = vvs
+ovvs.long_name = "meridional velocity"
+ovvs.units = "m/s"
+#
+omas[:] = mas
+omas.long_name = "mass"
+omas.units = "kg"
+#
+oths[:] = ths
+oths.long_name = "thickness"
+oths.units = "m"
+#
+owis[:] = wis
+owis.long_name = "width"
+owis.units = "m"
+#
+oles[:] = les
+oles.long_name = "length"
+oles.units = "m"
+#
+odys[:] = dys
+odys.long_name = "year day of calving event"
+odys.units = "days"
+#
+omss[:] = mss
+omss.long_name = "scaling factor for mass of calving berg"
+omss.units = "none"
+#
+omsb[:] = msb
+omsb.long_name = "mass of bergy bits"
+omsb.units = "kg"
+#
+ohds[:] = hds
+ohds.long_name = "heat density"
+ohds.units = "J/kg"
+#
+oyrs[:] = yrs
+oyrs.long_name = "calendar year of calving event"
+oyrs.units = "years"
+#
+onum[:,:] = num
+onum.long_name = "iceberg number on this processor"
+onum.units = "count"
+#
+fo.close()
Index: /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/SCOORD_GEN/namelist
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/SCOORD_GEN/namelist	(revision 6071)
+++ /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/SCOORD_GEN/namelist	(revision 6071)
@@ -0,0 +1,28 @@
+!-----------------------------------------------------------------------
+&namzgr_sco    !   s-coordinate or hybrid z-s-coordinate
+!-----------------------------------------------------------------------
+   ln_s_sh94   = .false.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)|
+   ln_s_sf12   = .true.   !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied
+   ln_sigcrit  = .true.    !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
+                           !  stretching coefficients for all functions
+   rn_jpk      =  51       ! Number of S levels
+   ln_eq_taper = .false.   !  Tapering of S coords near equator
+   cn_coord_hgr = 'coordinates.nc'  ! File containing gphit (latitude) coordinate for use if ln_eq_taper=.true.
+   rn_sbot_min =   10.0    !  minimum depth of s-bottom surface (>0) (m)
+   rn_sbot_max = 7000.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m)
+   rn_hc       =   50.0    !  critical depth for transition to stretched coordinates
+                        !!!!!!!  Envelop bathymetry
+   rn_rmax     =    0.3    !  maximum cut-off r-value allowed (0<r_max<1)
+                        !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.)
+   rn_theta    =    6.0    !  surface control parameter (0<=theta<=20)
+   rn_bb       =    0.8    !  stretching with SH94 s-sigma   
+                        !!!!!!!  SF12 stretching coefficient  (ln_s_sf12 = .true.)
+   rn_alpha    =    4.4    !  stretching with SF12 s-sigma
+   rn_efold    =    0.0    !  efold length scale for transition to stretched coord
+   rn_zs       =    1.0    !  depth of surface grid box
+                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
+   rn_zb_a     =    0.024  !  bathymetry scaling factor for calculating Zb
+   rn_zb_b     =   -0.2    !  offset for calculating Zb
+                        !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above]
+   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1) 
+/
Index: /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/SCOORD_GEN/src/scoord_gen.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/SCOORD_GEN/src/scoord_gen.F90	(revision 6071)
+++ /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/SCOORD_GEN/src/scoord_gen.F90	(revision 6071)
@@ -0,0 +1,664 @@
+PROGRAM SCOORD_GEN
+      !!----------------------------------------------------------------------
+      !!                  ***  PROGRAM scoord_gen  ***
+      !!                     
+      !! ** Purpose :   define the s-coordinate system
+      !!
+      !! ** Method  :   s-coordinate
+      !!         The depth of model levels is defined as the product of an
+      !!      analytical function by the local bathymetry, while the vertical
+      !!      scale factors are defined as the product of the first derivative
+      !!      of the analytical function by the bathymetry.
+      !!      (this solution save memory as depth and scale factors are not
+      !!      3d fields)
+      !!          - Read bathymetry (in meters) at t-point and compute the
+      !!         bathymetry at u-, v-, and f-points.
+      !!            hbatu = mi( hbatt )
+      !!            hbatv = mj( hbatt )
+      !!            hbatf = mi( mj( hbatt ) )
+      !!          - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
+      !!         function and its derivative given as function.
+      !!            z_gsigt(k) = fssig (k    )
+      !!            z_gsigw(k) = fssig (k-0.5)
+      !!            z_esigt(k) = fsdsig(k    )
+      !!            z_esigw(k) = fsdsig(k-0.5)
+      !!      Three options for stretching are give, and they can be modified
+      !!      following the users requirements. Nevertheless, the output as
+      !!      well as the way to compute the model levels and scale factors
+      !!      must be respected in order to insure second order accuracy
+      !!      schemes.
+      !!
+      !!      The three methods for stretching available are:
+      !! 
+      !!           s_sh94 (Song and Haidvogel 1994)
+      !!                a sinh/tanh function that allows sigma and stretched sigma
+      !!
+      !!           s_sf12 (Siddorn and Furner 2012?)
+      !!                allows the maintenance of fixed surface and or
+      !!                bottom cell resolutions (cf. geopotential coordinates) 
+      !!                within an analytically derived stretched S-coordinate framework.
+      !! 
+      !!          s_tanh  (Madec et al 1996)
+      !!                a cosh/tanh function that gives stretched coordinates        
+      !!
+      !! ** History: 2015: Tim Graham - Code created based on online zdf_sco routine
+      !!
+      !!
+      !!----------------------------------------------------------------------
+      !
+      USE utils
+      IMPLICIT NONE
+
+
+     !!----------------------------------------------------------------------
+      !
+      !
+      OPEN( UNIT=numnam, FILE='namelist', FORM='FORMATTED', STATUS='OLD' )
+      READ( numnam, namzgr_sco )
+      CLOSE( numnam )
+      jpk = INT(rn_jpk)
+
+      WRITE(*,*)
+      WRITE(*,*) 'scoord_gen : s-coordinate or hybrid z-s-coordinate'
+      WRITE(*,*) '~~~~~~~~~~~'
+      WRITE(*,*) '   Namelist namzgr_sco'
+      WRITE(*,*) '     stretching coeffs '
+      WRITE(*,*) '        Number of levels                             rn_jpk        = ',rn_jpk
+      WRITE(*,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max
+      WRITE(*,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min
+      WRITE(*,*) '        Critical depth                               rn_hc         = ',rn_hc
+      WRITE(*,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax
+      WRITE(*,*) '        Tapering in vicinity of equator              ln_eq_taper   = ',ln_eq_taper
+      WRITE(*,*) '        Horizontal Coordinate File                   cn_coord_hgr  = ',cn_coord_hgr
+      WRITE(*,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94
+      WRITE(*,*) '        Song and Haidvogel 1994 stretching coefficients'
+      WRITE(*,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta
+      WRITE(*,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb
+      WRITE(*,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb
+      WRITE(*,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12
+      WRITE(*,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit
+      WRITE(*,*) '        Siddorn and Furner 2012 stretching coefficients'
+      WRITE(*,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha
+      WRITE(*,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold
+      WRITE(*,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs
+      WRITE(*,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a
+      WRITE(*,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b
+      WRITE(*,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'
+
+      
+      ! Read in bathy, jpi and jpj from bathy.nc
+      CALL read_bathy()
+
+      !Allocate all other allocatable variables
+      ios = dom_oce_alloc()
+      IF( ios .NE. 0) THEN
+         WRITE(0,*) 'Unable to allocate all arrays'
+         STOP 1
+      ENDIF
+
+      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
+      hifu(:,:) = rn_sbot_min
+      hifv(:,:) = rn_sbot_min
+      hiff(:,:) = rn_sbot_min
+
+      !                                        ! set maximum ocean depth
+      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+           IF( bathy(ji,jj) > 0. )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
+         END DO
+      END DO
+      !                                        ! =============================
+      !                                        ! Define the envelop bathymetry   (hbatt)
+      !                                        ! =============================
+      ! use r-value to create hybrid coordinates
+      zenv(:,:) = bathy(:,:)
+      !
+      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+           IF( bathy(ji,jj) == 0. ) THEN
+             iip1 = MIN( ji+1, jpi )
+             ijp1 = MIN( jj+1, jpj )
+             iim1 = MAX( ji-1, 1 )
+             ijm1 = MAX( jj-1, 1 )
+             IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) +              &
+        &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0. ) THEN
+               zenv(ji,jj) = rn_sbot_min
+             ENDIF
+           ENDIF
+         END DO
+      END DO
+      ! 
+      ! smooth the bathymetry (if required)
+      scosrf(:,:) = 0.             ! ocean surface depth (here zero: no under ice-shelf sea)
+      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
+      !
+      jl = 0
+      zrmax = 1.
+      !   
+      !     
+      ! set scaling factor used in reducing vertical gradients
+      zrfact = ( 1. - rn_rmax ) / ( 1. + rn_rmax )
+      !
+      ! initialise temporary evelope depth arrays
+      ztmpi1(:,:) = zenv(:,:)
+      ztmpi2(:,:) = zenv(:,:)
+      ztmpj1(:,:) = zenv(:,:)
+      ztmpj2(:,:) = zenv(:,:)
+      !
+      ! initialise temporary r-value arrays
+      zri(:,:) = 1.
+      zrj(:,:) = 1.
+
+      !                                                            ! ================ !
+      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8 ) !  Iterative loop  !
+         !                                                         ! ================ !
+         jl = jl + 1
+         zrmax = 0.
+         ! we set zrmax from previous r-values (zri and zrj) first
+         ! if set after current r-value calculation (as previously)
+         ! we could exit DO WHILE prematurely before checking r-value
+         ! of current zenv
+         DO jj = 1, jpj 
+            DO ji = 1, jpi
+               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
+            END DO
+         END DO
+         zri(:,:) = 0.
+         zrj(:,:) = 0.
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               iip1 = MIN(ji+1,jpi) 
+               ijp1 = MIN(jj+1,jpj)     
+               IF( (zenv(ji,jj) > 0.) .AND. (zenv(iip1,jj) > 0.)) THEN
+                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
+               END IF
+               IF( (zenv(ji,jj) > 0.) .AND. (zenv(ji,ijp1) > 0.)) THEN
+                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
+               END IF
+               IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
+               IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact
+               IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
+               IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
+            END DO
+         END DO
+         !
+         WRITE(*,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
+         !
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
+            END DO
+         END DO
+         !                                                  ! ================ !
+      END DO                                                !     End loop     !
+      !                                                     ! ================ !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
+         END DO
+      END DO
+      !
+      ! Envelope bathymetry saved in hbatt
+      hbatt(:,:) = zenv(:,:) 
+      IF( ln_eq_taper) THEN
+        CALL READ_GPHIT()
+        IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0. ) THEN
+          WRITE(*,*) 's-coordinates are tapered in vicinity of the Equator' 
+          DO jj = 1, jpj
+             DO ji = 1, jpi
+                ztaper = EXP( -(gphit(ji,jj)/8.)**2. )
+                hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1. - ztaper )
+             END DO
+          END DO
+        ENDIF
+      ENDIF
+      !
+      !                                        ! ==============================
+      !                                        !   hbatu, hbatv, hbatf fields
+      !                                        ! ==============================
+      WRITE(*,*)
+      WRITE(*,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
+      hbatu(:,:) = rn_sbot_min
+      hbatv(:,:) = rn_sbot_min
+      hbatf(:,:) = rn_sbot_min
+      DO jj = 1, jpj-1
+        DO ji = 1, jpi-1   ! NO vector opt.
+           hbatu(ji,jj) = 0.50 * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
+           hbatv(ji,jj) = 0.50 * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
+           hbatf(ji,jj) = 0.25 * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
+              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
+        END DO
+      END DO
+      ! 
+      zhbat(:,:) = hbatu(:,:)   
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF( hbatu(ji,jj) == 0. ) THEN
+               IF( zhbat(ji,jj) == 0. )   hbatu(ji,jj) = rn_sbot_min
+               IF( zhbat(ji,jj) /= 0. )   hbatu(ji,jj) = zhbat(ji,jj)
+            ENDIF
+         END DO
+      END DO
+      zhbat(:,:) = hbatv(:,:)  
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF( hbatv(ji,jj) == 0. ) THEN
+               IF( zhbat(ji,jj) == 0. )   hbatv(ji,jj) = rn_sbot_min
+               IF( zhbat(ji,jj) /= 0. )   hbatv(ji,jj) = zhbat(ji,jj)
+            ENDIF
+         END DO
+      END DO
+      zhbat(:,:) = hbatf(:,:)   
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF( hbatf(ji,jj) == 0. ) THEN
+               IF( zhbat(ji,jj) == 0. )   hbatf(ji,jj) = rn_sbot_min
+               IF( zhbat(ji,jj) /= 0. )   hbatf(ji,jj) = zhbat(ji,jj)
+            ENDIF
+         END DO
+      END DO
+
+!!bug:  key_helsinki a verifer
+      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
+      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
+      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
+      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
+
+         WRITE(*,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
+            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
+         WRITE(*,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
+            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
+         WRITE(*,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
+            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
+         WRITE(*,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
+            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
+      
+      ! Create the output file
+      CALL make_coord_file()
+
+      !                                            ! =======================
+      !                                            !   s-coordinate fields     (gdep., e3.)
+      !                                            ! =======================
+      !
+      ! non-dimensional "sigma" for model level depth at w- and t-levels
+
+
+!========================================================================
+! Song and Haidvogel  1994 (ln_s_sh94=T)
+! Siddorn and Furner 2012 (ln_sf12=T)
+! or  tanh function       (both false)                    
+! To reduce memory loop over jpk and write each level to file
+!========================================================================
+      IF      ( ln_s_sh94 ) THEN 
+                           CALL s_sh94()
+      ELSE IF ( ln_s_sf12 ) THEN
+                           CALL s_sf12()
+      ELSE                 
+                           CALL s_tanh()
+      ENDIF
+
+!Write all 2D variables to output file
+      CALL write_netcdf_2d_vars()
+      CALL check_nf90( nf90_close(ncout) )
+
+ 
+      !
+END PROGRAM SCOORD_GEN
+
+!!======================================================================
+   SUBROUTINE s_sh94()
+
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE s_sh94  ***
+      !!                     
+      !! ** Purpose :   stretch the s-coordinate system
+      !!
+      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994
+      !!                mixed S/sigma coordinate
+      !!
+      !! Reference : Song and Haidvogel 1994. 
+      !!----------------------------------------------------------------------
+      !
+      USE utils
+      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
+      !
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigt3, z_esigw3, z_esigtu3 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
+
+      ALLOCATE( z_gsigw3(jpi,jpj), z_gsigt3(jpi,jpj), z_gsi3w3(jpi,jpj) )
+      ALLOCATE( z_esigt3(jpi,jpj), z_esigw3(jpi,jpj), z_esigtu3(jpi,jpj)) 
+      ALLOCATE( z_esigtv3(jpi,jpj), z_esigtf3(jpi,jpj), z_esigwu3(jpi,jpj)) 
+      ALLOCATE( z_esigwv3(jpi,jpj), z_gsigw3p1(jpi,jpj), z_gsigt3m1(jpi,jpj) )
+      ALLOCATE( z_gsi3w3m1(jpi,jpj) )
+ 
+      z_gsigw3  = 0.   ;   z_gsigt3  = 0.   ;   z_gsi3w3  = 0.
+      z_esigt3  = 0.   ;   z_esigw3  = 0. 
+      z_esigtu3 = 0.   ;   z_esigtv3 = 0.   ;   z_esigtf3 = 0.
+      z_esigwu3 = 0.   ;   z_esigwv3 = 0.
+
+      DO jk = 1,jpk
+         DO ji = 1, jpi
+            DO jj = 1, jpj
+
+               IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
+                     z_gsigw3(ji,jj) = -fssig1( REAL(jk,wp)-0.5, rn_bb )
+                     z_gsigw3p1(ji,jj) = -fssig1( REAL(jk+1,wp)-0.5, rn_bb )
+                     z_gsigt3(ji,jj) = -fssig1( REAL(jk,wp)    , rn_bb )
+               ELSE ! shallow water, uniform sigma
+                     z_gsigw3(ji,jj) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
+                     z_gsigw3p1(ji,jj) =   REAL(jk,wp)            / REAL(jpk-1,wp)
+                     z_gsigt3(ji,jj) = ( REAL(jk-1,wp) + 0.5 ) / REAL(jpk-1,wp)
+               ENDIF
+               !
+             !gsi3w3m1 & gsigt3m1 only used if jk /= 1 and is set at the end of the loop over jk
+               IF( jk .EQ. 1) THEN
+                  z_esigw3(ji,jj ) = 2. * ( z_gsigt3(ji,jj ) - z_gsigw3(ji,jj ) )
+                  z_gsi3w3(ji,jj) = 0.5 * z_esigw3(ji,jj)
+               ELSE
+                  z_esigw3(ji,jj) = z_gsigt3(ji,jj) - z_gsigt3m1(ji,jj)
+                  z_gsi3w3(ji,jj) = z_gsi3w3m1(ji,jj) + z_esigw3(ji,jj)
+               ENDIF
+               IF( jk .EQ. jpk) THEN
+                  z_esigt3(ji,jj) = 2. * ( z_gsigt3(ji,jj) - z_gsigw3(ji,jj) )
+               ELSE
+                  z_esigt3(ji,jj ) = z_gsigw3p1(ji,jj) - z_gsigw3(ji,jj)
+               ENDIF
+               !
+
+               zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp)
+               zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp)
+               gdept_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj)+rn_hc*zcoeft )
+               gdepw_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj)+rn_hc*zcoefw )
+               gdep3w_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj)+rn_hc*zcoeft )
+              !
+            END DO   ! for all jj's
+         END DO    ! for all ji's
+
+         DO ji = 1, jpi-1
+            DO jj = 1, jpj-1
+               z_esigtu3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) )   &
+                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+               z_esigtv3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1) )   &
+                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+               z_esigtf3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj)     &
+                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1) )   &
+                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
+               z_esigwu3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj) )   &
+                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+               z_esigwv3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1) )   &
+                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+               !
+               e3t_0(ji,jj) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj) + rn_hc/REAL(jpk-1,wp) )
+               e3u_0(ji,jj) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
+               e3v_0(ji,jj) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
+               e3f_0(ji,jj) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
+               !
+               e3w_0 (ji,jj) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj) + rn_hc/REAL(jpk-1,wp) )
+               e3uw_0(ji,jj) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
+               e3vw_0(ji,jj) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
+            END DO
+         END DO
+
+         z_gsigt3m1 = z_gsigt3
+         z_gsi3w3m1 = z_gsi3w3
+ 
+         where (e3t_0   (:,:).eq.0.0)  e3t_0(:,:) = 1.0
+         where (e3u_0   (:,:).eq.0.0)  e3u_0(:,:) = 1.0
+         where (e3v_0   (:,:).eq.0.0)  e3v_0(:,:) = 1.0
+         where (e3f_0   (:,:).eq.0.0)  e3f_0(:,:) = 1.0
+         where (e3w_0   (:,:).eq.0.0)  e3w_0(:,:) = 1.0
+         where (e3uw_0  (:,:).eq.0.0)  e3uw_0(:,:) = 1.0
+         where (e3vw_0  (:,:).eq.0.0)  e3vw_0(:,:) = 1.0
+        
+         CALL write_netcdf_3d_vars(jk)
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+                  IF( scobot(ji,jj) >= gdept_0(ji,jj) )   mbathy(ji,jj) = MAX( 2, jk )
+                  IF( scobot(ji,jj) == 0.             )   mbathy(ji,jj) = 0
+            END DO
+         END DO
+      END DO !End of loop over jk
+
+      DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3, z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1)
+      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
+
+   END SUBROUTINE s_sh94
+
+   SUBROUTINE s_sf12
+
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE s_sf12 *** 
+      !!                     
+      !! ** Purpose :   stretch the s-coordinate system
+      !!
+      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
+      !!                mixed S/sigma/Z coordinate
+      !!
+      !!                This method allows the maintenance of fixed surface and or
+      !!                bottom cell resolutions (cf. geopotential coordinates) 
+      !!                within an analytically derived stretched S-coordinate framework.
+      !!
+      !!
+      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
+      !!----------------------------------------------------------------------
+      !
+      USE utils
+      REAL(wp) ::   zsmth               ! smoothing around critical depth
+      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
+      !
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigt3, z_esigw3, z_esigtu3 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
+
+      ALLOCATE( z_gsigw3(jpi,jpj), z_gsigt3(jpi,jpj), z_gsi3w3(jpi,jpj) )
+      ALLOCATE( z_esigt3(jpi,jpj), z_esigw3(jpi,jpj), z_esigtu3(jpi,jpj)) 
+      ALLOCATE( z_esigtv3(jpi,jpj), z_esigtf3(jpi,jpj), z_esigwu3(jpi,jpj)) 
+      ALLOCATE( z_esigwv3(jpi,jpj), z_gsigw3p1(jpi,jpj), z_gsigt3m1(jpi,jpj) )
+      ALLOCATE( z_gsi3w3m1(jpi,jpj) )
+
+      z_gsigw3  = 0.   ;   z_gsigt3  = 0.   ;   z_gsi3w3  = 0.
+      z_esigt3  = 0.   ;   z_esigw3  = 0. 
+      z_esigtu3 = 0.   ;   z_esigtv3 = 0.   ;   z_esigtf3 = 0.
+      z_esigwu3 = 0.   ;   z_esigwv3 = 0.
+      
+
+      DO jk = 1,jpk
+         DO ji = 1, jpi
+            DO jj = 1, jpj
+             IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
+              
+                zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
+                                                     ! could be changed by users but care must be taken to do so carefully
+                zzb = 1.0-(zzb/hbatt(ji,jj))
+            
+                zzs = rn_zs / hbatt(ji,jj) 
+              
+                IF (rn_efold /= 0.0) THEN
+                   zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
+                ELSE
+                   zsmth = 1.0 
+                ENDIF
+                  
+                z_gsigw3(ji,jj) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
+                z_gsigw3p1(ji,jj) =  REAL(jk,wp)        /REAL(jpk-1,wp)
+                z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
+                z_gsigw3(ji,jj) = fgamma( z_gsigw3(ji,jj), zzb, zzs, zsmth  )
+                z_gsigw3p1(ji,jj) = fgamma( z_gsigw3p1(ji,jj), zzb, zzs, zsmth  )
+                z_gsigt3(ji,jj) = fgamma( z_gsigt3(ji,jj), zzb, zzs, zsmth  )
+ 
+             ELSE IF(ln_sigcrit) THEN ! shallow water, uniform sigma
+
+                z_gsigw3(ji,jj) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
+                z_gsigw3p1(ji,jj) =  REAL(jk,wp)     /REAL(jpk-1,wp)
+                z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
+
+             ELSE  ! shallow water, z coordinates
+
+               z_gsigw3(ji,jj) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
+               z_gsigw3p1(ji,jj) =  REAL(jk,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
+               z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
+
+             ENDIF
+
+             !gsi3w3m1 & z_gsigt3m1 only used if jk /= 1 and is set at the end of the loop over jk
+             IF( jk .EQ. 1) THEN
+                z_esigw3(ji,jj  ) = 2.0 * (z_gsigt3(ji,jj ) - z_gsigw3(ji,jj ))
+                z_gsi3w3(ji,jj) = 0.5 * z_esigw3(ji,jj) 
+             ELSE
+                z_esigw3(ji,jj) = z_gsigt3(ji,jj) - z_gsigt3m1(ji,jj)
+                z_gsi3w3(ji,jj) = z_gsi3w3m1(ji,jj) + z_esigw3(ji,jj)
+             ENDIF
+             IF( jk .EQ. jpk) THEN
+                z_esigt3(ji,jj) = 2.0 * (z_gsigt3(ji,jj) - z_gsigw3(ji,jj))
+             ELSE
+                z_esigt3(ji,jj) = z_gsigw3p1(ji,jj) - z_gsigw3(ji,jj)
+             ENDIF
+
+             gdept_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj)
+             gdepw_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj)
+             gdep3w_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj)
+
+            ENDDO   ! for all jj's
+         ENDDO    ! for all ji's
+
+         DO ji=1,jpi-1
+            DO jj=1,jpj-1
+
+               z_esigtu3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) ) / &
+                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+               z_esigtv3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1) ) / &
+                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+               z_esigtf3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) +  &
+                                       hbatt(ji,jj+1)*z_esigt3(ji,jj+1)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1) ) / &
+                                     ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
+               z_esigwu3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj) ) / &
+                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+               z_esigwv3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1) ) / &
+                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+
+               e3t_0(ji,jj)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj)
+               e3u_0(ji,jj)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj)
+               e3v_0(ji,jj)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj)
+               e3f_0(ji,jj)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj)
+               !
+               e3w_0(ji,jj)=hbatt(ji,jj)*z_esigw3(ji,jj)
+               e3uw_0(ji,jj)=hbatu(ji,jj)*z_esigwu3(ji,jj)
+               e3vw_0(ji,jj)=hbatv(ji,jj)*z_esigwv3(ji,jj)
+            ENDDO
+         ENDDO
+         ! Keep some arrays for next level
+         z_gsigt3m1 = z_gsigt3
+         z_gsi3w3m1 = z_gsi3w3
+ 
+         where (e3t_0   (:,:).eq.0.0)  e3t_0(:,:) = 1.0
+         where (e3u_0   (:,:).eq.0.0)  e3u_0(:,:) = 1.0
+         where (e3v_0   (:,:).eq.0.0)  e3v_0(:,:) = 1.0
+         where (e3f_0   (:,:).eq.0.0)  e3f_0(:,:) = 1.0
+         where (e3w_0   (:,:).eq.0.0)  e3w_0(:,:) = 1.0
+         where (e3uw_0  (:,:).eq.0.0)  e3uw_0(:,:) = 1.0
+         where (e3vw_0  (:,:).eq.0.0)  e3vw_0(:,:) = 1.0
+         
+         CALL write_netcdf_3d_vars(jk)
+
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+                  IF( scobot(ji,jj) >= gdept_0(ji,jj) )   mbathy(ji,jj) = MAX( 2, jk )
+                  IF( scobot(ji,jj) == 0.             )   mbathy(ji,jj) = 0
+            END DO
+         END DO
+
+      ENDDO ! End of loop over jk
+
+      DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3, z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1)
+      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
+
+   END SUBROUTINE s_sf12
+
+   SUBROUTINE s_tanh()
+
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE s_tanh*** 
+      !!                     
+      !! ** Purpose :   stretch the s-coordinate system
+      !!
+      !! ** Method  :   s-coordinate stretch 
+      !!
+      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
+      !!----------------------------------------------------------------------
+
+      USE utils
+      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
+
+      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
+      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_esigt, z_esigw
+
+      ALLOCATE( z_gsigw(jpk), z_gsigt(jpk), z_gsi3w(jpk) )
+      ALLOCATE( z_esigt(jpk), z_esigw(jpk) )
+
+      z_gsigw  = 0.   ;   z_gsigt  = 0.   ;   z_gsi3w  = 0.
+      z_esigt  = 0.   ;   z_esigw  = 0. 
+
+      DO jk = 1, jpk
+        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5 )
+        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
+      END DO
+      WRITE(*,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk)
+      !
+      ! Coefficients for vertical scale factors at w-, t- levels
+!!gm bug :  define it from analytical function, not like juste bellow....
+!!gm        or betteroffer the 2 possibilities....
+      DO jk = 1, jpk-1
+         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
+         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
+      END DO
+      z_esigw( 1 ) = 2. * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
+      z_esigt(jpk) = 2. * ( z_gsigt(jpk) - z_gsigw(jpk) )
+      !
+      ! Coefficients for vertical depth as the sum of e3w scale factors
+      z_gsi3w(1) = 0.5 * z_esigw(1)
+      DO jk = 2, jpk
+         z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
+      END DO
+!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+              zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp)
+              zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp)
+              gdept_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsigt(jk) + hift(ji,jj)*zcoeft )
+              gdepw_0(:,:) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsigw(jk) + hift(ji,jj)*zcoefw )
+              gdep3w_0(:,:) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsi3w(jk) + hift(ji,jj)*zcoeft )
+              e3t_0(ji,jj) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpk-1,wp) )
+              e3u_0(ji,jj) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpk-1,wp) )
+              e3v_0(ji,jj) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpk-1,wp) )
+              e3f_0(ji,jj) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpk-1,wp) )
+              !
+              e3w_0(ji,jj) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpk-1,wp) )
+              e3uw_0(ji,jj) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpk-1,wp) )
+              e3vw_0(ji,jj) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpk-1,wp) )
+              IF( scobot(ji,jj) >= gdept_0(ji,jj) )   mbathy(ji,jj) = MAX( 2, jk )
+              IF( scobot(ji,jj) == 0.             )   mbathy(ji,jj) = 0
+            END DO
+         END DO
+         where (e3t_0   (:,:).eq.0.0)  e3t_0(:,:) = 1.0
+         where (e3u_0   (:,:).eq.0.0)  e3u_0(:,:) = 1.0
+         where (e3v_0   (:,:).eq.0.0)  e3v_0(:,:) = 1.0
+         where (e3f_0   (:,:).eq.0.0)  e3f_0(:,:) = 1.0
+         where (e3w_0   (:,:).eq.0.0)  e3w_0(:,:) = 1.0
+         where (e3uw_0  (:,:).eq.0.0)  e3uw_0(:,:) = 1.0
+         where (e3vw_0  (:,:).eq.0.0)  e3vw_0(:,:) = 1.0
+ 
+         CALL write_netcdf_3d_vars(jk)
+      ENDDO ! End of loop over jk
+
+
+      DEALLOCATE( z_gsigw, z_gsigt, z_gsi3w )
+      DEALLOCATE( z_esigt, z_esigw )
+
+   END SUBROUTINE s_tanh
Index: /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/SCOORD_GEN/src/utils.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/SCOORD_GEN/src/utils.F90	(revision 6071)
+++ /branches/2015/dev_merge_2015/NEMOGCM/TOOLS/SCOORD_GEN/src/utils.F90	(revision 6071)
@@ -0,0 +1,330 @@
+MODULE utils
+
+   USE netcdf
+
+   IMPLICIT NONE
+   PUBLIC             ! allows the acces to par_oce when dom_oce is used
+   !                  ! exception to coding rules... to be suppressed ???
+
+!   PUBLIC dom_oce_alloc
+
+   INTEGER, PARAMETER   :: dp=8 , sp=4, wp=dp
+
+   !! All coordinates
+   !! ---------------
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gdep3w_0           !: depth of t-points (sum of e3w) (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gdept_0, gdepw_0   !: analytical (time invariant) depth at t-w  points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gphit              !: latitude at t points
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3v_0  , e3f_0     !: analytical (time invariant) vertical scale factors at  v-f
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3t_0  , e3u_0     !:                                      t-u  points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3vw_0             !: analytical (time invariant) vertical scale factors at  vw
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3w_0  , e3uw_0    !:                                      w-uw points (m)
+
+   !! s-coordinate and hybrid z-s-coordinate
+   !! =----------------======---------------
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hbatv , hbatf      !: ocean depth at the vertical of  v--f
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hbatt , hbatu      !:                                 t--u points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   scosrf, scobot     !: ocean surface and bottom topographies 
+   !                                                                           !  (if deviating from coordinate surfaces in HYBRID)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hifv  , hiff       !: interface depth between stretching at v--f
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hift  , hifu       !: and quasi-uniform spacing             t--u points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rx1                !: Maximum grid stiffness ratio
+
+   !!----------------------------------------------------------------------
+   !! masks, bathymetry
+   !! ---------------------------------------------------------------------
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbathy             !: number of ocean level (=0, 1, ... , jpk-1)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy              !: ocean depth (meters - read from file)
+
+   !! Other variables needed by scoord_gen
+   INTEGER  ::   jpi, jpj, jpk            ! Size of the domain - read from bathy or namelist?
+   INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
+   INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
+   INTEGER  ::   ios                      ! Local integer output status for namelist read and allocation
+   INTEGER,PARAMETER  ::   numnam=8       ! File handle for namelist 
+   REAL(wp) ::   zrmax, ztaper   ! temporary scalars
+   REAL(wp) ::   zrfact
+   !
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zenv, ztmp, zmsk, zri, zrj, zhbat
+
+   !Namelist variables
+   REAL(wp) :: rn_jpk, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax, rn_theta
+   REAL(wp) :: rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
+   LOGICAL :: ln_s_sh94, ln_s_sf12, ln_sigcrit, ln_eq_taper
+   CHARACTER(len=50) :: cn_coord_hgr
+
+   NAMELIST/namzgr_sco/rn_jpk, ln_s_sh94, ln_s_sf12, ln_sigcrit, ln_eq_taper, & 
+                  &      cn_coord_hgr, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
+                  &      rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
+
+  ! IDs for output netcdf file
+  INTEGER :: id_x, id_y, id_z
+  INTEGER :: ncout
+  INTEGER, DIMENSION(20) :: var_ids  !Array to contain all variable IDs
+
+   CONTAINS
+
+   INTEGER FUNCTION dom_oce_alloc()
+      !!----------------------------------------------------------------------
+      INTEGER, DIMENSION(4) :: ierr
+      !!----------------------------------------------------------------------
+      ierr(:) = 0
+      !
+      ALLOCATE( zenv(jpi,jpj), ztmp(jpi,jpj), zmsk(jpi,jpj), zri(jpi,jpj), zrj(jpi,jpj), &
+         &      zhbat(jpi,jpj) , ztmpi1(jpi,jpj), ztmpi2(jpi,jpj), ztmpj1(jpi,jpj), ztmpj2(jpi,jpj), STAT=ierr(1) )
+         !
+      ALLOCATE( gdep3w_0(jpi,jpj) , e3v_0(jpi,jpj) , e3f_0(jpi,jpj) ,                         &
+         &      gdept_0(jpi,jpj) , e3t_0(jpi,jpj) , e3u_0 (jpi,jpj) ,                         &
+         &      gdepw_0(jpi,jpj) , e3w_0(jpi,jpj) , e3vw_0(jpi,jpj) ,                         &
+         &      gphit(jpi,jpj)   , e3uw_0(jpi,jpj) , STAT=ierr(2) )
+         !
+         !
+      ALLOCATE( hbatv (jpi,jpj) , hbatf (jpi,jpj) ,     &
+         &      hbatt (jpi,jpj) , hbatu (jpi,jpj) ,     &
+         &      scosrf(jpi,jpj) , scobot(jpi,jpj) ,     &
+         &      hifv  (jpi,jpj) , hiff  (jpi,jpj) ,     &
+         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(3) )
+
+      ALLOCATE( mbathy(jpi,jpj) , STAT=ierr(4) )
+     !
+      dom_oce_alloc = MAXVAL(ierr)
+      !
+   END FUNCTION dom_oce_alloc
+ 
+
+   SUBROUTINE read_bathy()
+     !! Read bathymetry from input netcdf file
+     INTEGER :: var_id, ncin
+
+     CALL check_nf90( nf90_open('bathy_meter.nc', NF90_NOWRITE, ncin), 'Error opening bathy_meter.nc file' )
+
+     ! Find the size of the input bathymetry
+     CALL dimlen(ncin, 'lon', jpi)    
+     CALL dimlen(ncin, 'lat', jpj)    
+     
+     ALLOCATE( bathy(jpi, jpj) )
+     
+     ! Read the bathymetry variable from file
+     CALL check_nf90( nf90_inq_varid( ncin, 'Bathymetry', var_id ), 'Cannot get variable ID for bathymetry')
+     CALL check_nf90( nf90_get_var( ncin, var_id, bathy, (/ 1,1 /), (/ jpi, jpj /) ) )
+
+     CALL check_nf90( nf90_close(ncin), 'Error closing bathy.nc file' )
+
+   END SUBROUTINE read_bathy
+
+   SUBROUTINE read_gphit()
+   !! Read gphit from horizontal coordinate file if required
+     INTEGER :: var_id, ncin
+
+     CALL check_nf90( nf90_open(cn_coord_hgr, NF90_NOWRITE, ncin), 'Error opening horizontal coordinate file' )
+
+     ! Read gphit variable from file
+     CALL check_nf90( nf90_inq_varid( ncin, 'gphit', var_id ), 'Cannot get variable ID for bathymetry')
+     CALL check_nf90( nf90_get_var( ncin, var_id, gphit, (/ 1,1 /), (/ jpi, jpj /) ) )
+
+     CALL check_nf90( nf90_close(ncin), 'Error closing horizontal coordinate file' )
+
+   END SUBROUTINE read_gphit
+
+   SUBROUTINE dimlen( ncid, dimname, len )
+     ! Determine the length of dimension dimname
+     INTEGER, INTENT(in)          :: ncid
+     CHARACTER(LEN=*), INTENT(in) :: dimname
+     INTEGER, INTENT(out)         :: len
+     ! Local variables
+     INTEGER :: id_var
+
+     id_var = 1
+     CALL check_nf90( nf90_inq_dimid(ncid, dimname, id_var), 'Dimension not found in file')
+     CALL check_nf90( nf90_inquire_dimension(ncid,id_var,len=len))
+
+   END SUBROUTINE dimlen
+  
+ 
+   SUBROUTINE make_coord_file()
+     ! Create new coordinates file and define dimensions and variables ready for
+     ! writing
+     
+
+     !Create the file
+     CALL check_nf90( nf90_create('coord_zgr.nc', NF90_NETCDF4, ncout), 'Could not create output file')
+     !
+     !Define dimensions
+     CALL check_nf90( nf90_def_dim(ncout, 'x', jpi, id_x) )
+     CALL check_nf90( nf90_def_dim(ncout, 'y', jpj, id_y) )
+     CALL check_nf90( nf90_def_dim(ncout, 'z', jpk, id_z) )
+     !
+     !Define variables - include all varibles that would be put into the mesh
+     !mask file
+     CALL check_nf90( nf90_def_var(ncout, 'gdept_0', nf90_double, (/id_x, id_y,id_z/), var_ids(1)) )
+     CALL check_nf90( nf90_def_var(ncout, 'gdepw_0', nf90_double, (/id_x, id_y,id_z/), var_ids(2)) )
+     CALL check_nf90( nf90_def_var(ncout, 'gdep3w_0', nf90_double, (/id_x, id_y,id_z/), var_ids(3)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3f_0', nf90_double, (/id_x, id_y,id_z/), var_ids(4)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3t_0', nf90_double, (/id_x, id_y,id_z/), var_ids(5)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3u_0', nf90_double, (/id_x, id_y,id_z/), var_ids(6)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3v_0', nf90_double, (/id_x, id_y,id_z/), var_ids(7)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3w_0', nf90_double, (/id_x, id_y,id_z/), var_ids(8)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3uw_0', nf90_double, (/id_x, id_y,id_z/), var_ids(9)) )
+     CALL check_nf90( nf90_def_var(ncout, 'e3vw_0', nf90_double, (/id_x, id_y,id_z/), var_ids(10)) )
+     ! 2D fields
+     CALL check_nf90( nf90_def_var(ncout, 'mbathy', nf90_double, (/id_x, id_y/), var_ids(11)) )
+     CALL check_nf90( nf90_def_var(ncout, 'hbatt', nf90_double, (/id_x, id_y/), var_ids(12)) )
+     CALL check_nf90( nf90_def_var(ncout, 'hbatu', nf90_double, (/id_x, id_y/), var_ids(13)) )
+     CALL check_nf90( nf90_def_var(ncout, 'hbatv', nf90_double, (/id_x, id_y/), var_ids(14)) )
+     CALL check_nf90( nf90_def_var(ncout, 'hbatf', nf90_double, (/id_x, id_y/), var_ids(15)) )
+     CALL check_nf90( nf90_def_var(ncout, 'rx1', nf90_double, (/id_x, id_y/), var_ids(16)) )
+
+     
+     ! End define mode
+     CALL check_nf90( nf90_enddef(ncout) )
+     
+     WRITE(*,*) 'Opened coord_zgr.nc file and defined variables'
+
+   END SUBROUTINE make_coord_file
+
+   SUBROUTINE write_netcdf_2d_vars()
+
+     CALL check_nf90( nf90_put_var(ncout, var_ids(11), mbathy, (/ 1,1 /), (/ jpi, jpj /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(12), hbatt, (/ 1,1 /), (/ jpi, jpj /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(13), hbatu, (/ 1,1 /), (/ jpi, jpj /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(14), hbatv, (/ 1,1 /), (/ jpi, jpj /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(15), hbatf, (/ 1,1 /), (/ jpi, jpj /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(16), rx1, (/ 1,1 /), (/ jpi, jpj /) ) )
+
+   END SUBROUTINE write_netcdf_2d_vars
+
+   SUBROUTINE write_netcdf_3d_vars(kk)
+   ! Write  variables to the netcdf file at level kk
+     INTEGER, INTENT(in) :: kk
+
+     CALL check_nf90( nf90_put_var(ncout, var_ids(1), gdept_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(2), gdepw_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(3), gdep3w_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(4), e3f_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(5), e3t_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(6), e3u_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(7), e3v_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(8), e3w_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(9), e3uw_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
+     CALL check_nf90( nf90_put_var(ncout, var_ids(10), e3vw_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
+
+   END SUBROUTINE write_netcdf_3d_vars
+
+   SUBROUTINE check_nf90( istat, message )
+      !Check for netcdf errors
+      INTEGER, INTENT(in) :: istat
+      CHARACTER(LEN=*), INTENT(in), OPTIONAL :: message
+
+      IF (istat /= nf90_noerr) THEN
+         WRITE(*,*) 'ERROR! : '//TRIM(nf90_strerror(istat))
+         IF ( PRESENT(message) ) THEN ; WRITE(*,*) message ; ENDIF
+         STOP
+      ENDIF
+
+   END SUBROUTINE check_nf90
+   FUNCTION fssig( pk ) RESULT( pf )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE fssig ***
+      !!       
+      !! ** Purpose :   provide the analytical function in s-coordinate
+      !!          
+      !! ** Method  :   the function provide the non-dimensional position of
+      !!                T and W (i.e. between 0 and 1)
+      !!                T-points at integer values (between 1 and jpk)
+      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
+      !!----------------------------------------------------------------------
+!     USE utils, ONLY : wp,rn_theta,rn_thetb,jpk
+      IMPLICIT NONE
+      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
+      REAL(wp)             ::   pf   ! sigma value
+      !!----------------------------------------------------------------------
+      !
+      pf =   (   TANH( rn_theta * ( -(pk-0.5) / REAL(jpk-1) + rn_thetb )  )   &
+         &     - TANH( rn_thetb * rn_theta                                )  )   &
+         & * (   COSH( rn_theta                           )                      &
+         &     + COSH( rn_theta * ( 2. * rn_thetb - 1. ) )  )              &
+         & / ( 2. * SINH( rn_theta ) )
+      !
+   END FUNCTION fssig
+
+   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE fssig1 ***
+      !!
+      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
+      !!
+      !! ** Method  :   the function provides the non-dimensional position of
+      !!                T and W (i.e. between 0 and 1)
+      !!                T-points at integer values (between 1 and jpk)
+      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
+      !!----------------------------------------------------------------------
+!     USE utils, ONLY : wp, jpk, rn_theta
+      IMPLICIT NONE
+      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
+      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
+      REAL(wp)             ::   pf1   ! sigma value
+      !!----------------------------------------------------------------------
+      !
+      IF ( rn_theta == 0 ) then      ! uniform sigma
+         pf1 = - ( pk1 - 0.5 ) / REAL( jpk-1 )
+      ELSE                        ! stretched sigma
+         pf1 =   ( 1. - pbb ) * ( SINH( rn_theta*(-(pk1-0.5)/REAL(jpk-1)) ) ) / SINH( rn_theta )              &
+            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5)/REAL(jpk-1)) + 0.5) ) - TANH( 0.5 * rn_theta )  )  &
+            &        / ( 2. * TANH( 0.5 * rn_theta ) )  )
+      ENDIF
+      !
+   END FUNCTION fssig1
+
+   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE fgamma  ***
+      !!
+      !! ** Purpose :   provide analytical function for the s-coordinate
+      !!
+      !! ** Method  :   the function provides the non-dimensional position of
+      !!                T and W (i.e. between 0 and 1)
+      !!                T-points at integer values (between 1 and jpk)
+      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
+      !!
+      !!                This method allows the maintenance of fixed surface and or
+      !!                bottom cell resolutions (cf. geopotential coordinates) 
+      !!                within an analytically derived stretched S-coordinate framework.
+      !!
+      !! Reference  :   Siddorn and Furner, in prep
+      !!----------------------------------------------------------------------
+!     USE utils, ONLY : jpk,wp,rn_alpha
+      IMPLICIT NONE
+      REAL(wp), INTENT(in   ) ::   pk1           ! continuous "k" coordinate
+      REAL(wp)                ::   p_gamma       ! stretched coordinate
+      REAL(wp), INTENT(in   ) ::   pzb           ! Bottom box depth
+      REAL(wp), INTENT(in   ) ::   pzs           ! surface box depth
+      REAL(wp), INTENT(in   ) ::   psmth         ! Smoothing parameter
+      REAL(wp)                ::   za1,za2,za3   ! local variables
+      REAL(wp)                ::   zn1,zn2       ! local variables
+      REAL(wp)                ::   za,zb,zx      ! local variables
+      !!----------------------------------------------------------------------
+      !
+
+      zn1  =  1./(jpk-1.)
+      zn2  =  1. -  zn1
+
+      za1 = (rn_alpha+2.0)*zn1**(rn_alpha+1.0)-(rn_alpha+1.0)*zn1**(rn_alpha+2.0) 
+      za2 = (rn_alpha+2.0)*zn2**(rn_alpha+1.0)-(rn_alpha+1.0)*zn2**(rn_alpha+2.0)
+      za3 = (zn2**3.0 - za2)/( zn1**3.0 - za1)
+     
+      za = pzb - za3*(pzs-za1)-za2
+      za = za/( zn2-0.5*(za2+zn2**2.0) - za3*(zn1-0.5*(za1+zn1**2.0) ) )
+      zb = (pzs - za1 - za*( zn1-0.5*(za1+zn1**2.0 ) ) ) / (zn1**3.0 - za1)
+      zx = 1.0-za/2.0-zb
+
+      p_gamma = za*(pk1*(1.0-pk1/2.0))+zb*pk1**3.0 +  &
+                  & zx*( (rn_alpha+2.0)*pk1**(rn_alpha+1.0)- &
+                  &      (rn_alpha+1.0)*pk1**(rn_alpha+2.0) )
+      p_gamma = p_gamma*psmth+pk1*(1.0-psmth)
+
+      !
+   END FUNCTION fgamma
+
+
+END MODULE utils
Index: /branches/2015/dev_merge_2015/NEMOGCM/fcm-make/inc/keys-orca2_lim.cfg
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/fcm-make/inc/keys-orca2_lim.cfg	(revision 6070)
+++ /branches/2015/dev_merge_2015/NEMOGCM/fcm-make/inc/keys-orca2_lim.cfg	(revision 6071)
@@ -1,2 +1,2 @@
 preprocess.prop{fpp.defs} = \
-  key_trabbl key_lim2 key_dynspg_flt key_diaeiv key_ldfslp key_traldf_c2d key_traldf_eiv key_dynldf_c3d key_zdftke key_zdfddm key_zdftmx key_iomput key_mpp_mpi key_diaobs key_asminc
+  key_trabbl key_lim2 key_diaeiv key_ldfslp key_traldf_c2d key_traldf_eiv key_dynldf_c3d key_zdftke key_zdfddm key_zdftmx key_iomput key_mpp_mpi key_asminc
Index: /branches/2015/dev_merge_2015/NEMOGCM/fcm-make/inc/keys-orca2_lim3.cfg
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/fcm-make/inc/keys-orca2_lim3.cfg	(revision 6070)
+++ /branches/2015/dev_merge_2015/NEMOGCM/fcm-make/inc/keys-orca2_lim3.cfg	(revision 6071)
@@ -1,2 +1,2 @@
 preprocess.prop{fpp.defs} = \
-  key_trabbl key_lim3 key_vvl key_dynspg_ts key_diaeiv key_ldfslp key_traldf_c2d key_traldf_eiv key_dynldf_c3d key_zdftke key_zdfddm key_zdftmx key_iomput key_mpp_mpi key_diaobs key_asminc
+  key_trabbl key_lim3 key_vvl key_diaeiv key_ldfslp key_traldf_c2d key_traldf_eiv key_dynldf_c3d key_zdftke key_zdfddm key_zdftmx key_iomput key_mpp_mpi key_asminc
