Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 6043)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 6044)
@@ -652,5 +652,5 @@
             ELSE
                CALL ctl_stop('DIA_OBS: Generalised vertical interpolation not'// &
-                         'yet working for velocity date (turn off velocity observations')
+                         'yet working for velocity data (turn off velocity observations')
             ENDIF
 
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/obs_fbm.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/obs_fbm.F90	(revision 6043)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/obs_fbm.F90	(revision 6044)
@@ -45,8 +45,4 @@
    INTEGER, PARAMETER    :: fbimdi = -99999   !: Integers
    REAL(fbsp), PARAMETER :: fbrmdi =  99999   !: Reals
-
-   ! Output stream choice
-   LOGICAL               :: ln_cl4 = .FALSE.  !: Logical switch for
-                                              !: class 4 file outputs
  
    ! Main data structure for observation feedback data.
@@ -1030,31 +1026,4 @@
 
    SUBROUTINE write_obfbdata( cdfilename, fbdata )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE write_obfbdata  ***
-      !!
-      !! ** Purpose :   Write an obfbdata structure into a netCDF file.
-      !!
-      !! ** Method  :   Decides which output wrapper to use. 
-      !!
-      !! ** Action  : 
-      !!
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      CHARACTER(len=*) :: cdfilename ! Output filename
-      TYPE(obfbdata)   :: fbdata     ! obsfbdata structure
-#if defined key_offobsoper
-      IF (ln_cl4) THEN
-          ! Class 4 file output stream
-          CALL write_obfbdata_cl( cdfilename, fbdata )
-      ELSE
-#endif
-          ! Standard feedback file output stream
-          CALL write_obfbdata_fb( cdfilename, fbdata )
-#if defined key_offobsoper
-      ENDIF
-#endif
-   END SUBROUTINE write_obfbdata
-
-   SUBROUTINE write_obfbdata_fb( cdfilename, fbdata )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE write_obfbdata  ***
@@ -1555,237 +1524,5 @@
 
       
-   END SUBROUTINE write_obfbdata_fb
-
-#if defined key_offobsoper
-   SUBROUTINE write_obfbdata_cl(cdfilename, fbdata)
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE write_obfbdata_cl  ***
-      !!
-      !! ** Purpose : Write an obfbdata structure into a class 4 file.
-      !!
-      !! ** Method  : 1. Allocate memory needed by ooo_write
-      !!              2. Map obfbdata into allocated memory
-      !!              3. Pass mapped data to ooo_write
-      !!              4. Deallocate memory
-      !!----------------------------------------------------------------------
-      USE dom_oce, ONLY: narea
-      USE ooo_write
-      USE ooo_data
-      !! * Arguments
-      CHARACTER(len=*) :: cdfilename ! Feedback filename
-      TYPE(obfbdata)   :: fbdata     ! obsfbdata structure
-      !! * Local variables
-      CHARACTER(len=17), PARAMETER :: cpname = 'write_obfbdata_cl'
-      CHARACTER(len=64) :: &
-              & cdate, &   !: class 4 file validity date 
-              & cconf, &   !: model configuration
-              & csys, &    !: model system
-              & ccont, &   !: contact email
-              & cinst, &   !: institution
-              & cversion   !: model version
-      CHARACTER(len=8) :: &
-              & ckind      !: observation kind
-      CHARACTER(len=3) :: cfield
-      INTEGER :: kobs, &   !: number of observations
-              &  kvars, &  !: number of physical variables
-              &  kdeps, &  !: number of observed depths
-              &  kfcst, &  !: number of forecasts
-              &  kifcst, & !: current forecast number
-              &  kproc     !: processor number
-      INTEGER, DIMENSION(:, :, :), ALLOCATABLE :: &
-              &  kqc       !: quality control counterpart
-      INTEGER(KIND=2), DIMENSION(:, :, :), ALLOCATABLE :: &
-              &  k2qc       !: quality control counterpart
-      REAL(kind=fbdp) :: &
-              &  pmodjuld  !: model Julian day
-      REAL(kind=fbdp), DIMENSION(:), ALLOCATABLE :: &
-              &  plead, &  !: forecast lead time
-              &  plam, &   !: longitude of observation
-              &  pphi, &   !: latitude of observation
-              &  ptim      !: time of observation
-      REAL(kind=fbdp), DIMENSION(:, :), ALLOCATABLE :: &
-              &  pdep      !: depths of observations
-      REAL(kind=fbdp), DIMENSION(:, :, :), ALLOCATABLE :: &
-              &  pob, &    !: observation counterpart
-              &  pextra    !: extra field counterpart
-      REAL(kind=fbdp), DIMENSION(:, :, :), ALLOCATABLE :: &
-              &  pmod      !: model counterpart
-      CHARACTER(len=128) :: &
-              &  clfilename  !: class 4 file name
-      CHARACTER(len=128), DIMENSION(:), ALLOCATABLE :: &
-              &  ctype       !: Instrument type
-      CHARACTER(len=nf90_max_name) :: &
-              & cdtmp        !: NetCDF variable name
-      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: &
-              &  cwmo, &     !: Instrument WMO ID
-              &  cunit, &    !: Instrument WMO ID
-              &  cvarname    !: Instrument WMO ID
-      INTEGER :: &
-              &  idep, &     !: Loop variable
-              &  ivar, &     !: Loop variable
-              &  iobs, &     !: Loop variable
-              &  ii, &       !: Loop variable
-              &  ij, &       !: Loop variable
-              &  ik, &       !: Loop variable
-              &  il          !: Loop variable
-      cconf = TRIM(cl4_cfg)
-      csys = TRIM(cl4_sys)
-      cversion = TRIM(cl4_vn)
-      ccont = TRIM(cl4_contact)
-      cinst = TRIM(cl4_inst)
-      cdate = TRIM(cl4_date)
-      CALL locate_kind(cdfilename, ckind)
-      kproc = narea
-      kfcst = cl4_fcst_len
-      kobs = fbdata%nobs
-      kdeps = fbdata%nlev
-      kvars = fbdata%nvar
-      IF (kobs .GT. 0) THEN
-         ALLOCATE(plam(kobs), &
-               &  pphi(kobs), &
-               &  ptim(kobs), &
-               &  plead(kfcst), &
-               &  pdep(kdeps, kobs), &
-               &  kqc(kdeps, kvars, kobs), &
-               &  k2qc(kdeps, kvars, kobs), &
-               &  pob(kdeps, kvars, kobs), &
-               &  pmod(kdeps, kvars, kobs), &
-               &  pextra(kdeps, kvars, kobs), &
-               &  ctype(kobs), &
-               &  cwmo(kobs), &
-               &  cunit(kvars), &
-               &  cvarname(kvars))
-         plam(:) = fbdata%plam(:)
-         pphi(:) = fbdata%pphi(:)
-         ptim(:) = fbdata%ptim(:)
-         pdep(:, :) = fbdata%pdep(:, :)
-         kqc(:,:,:) = 1.
-         DO ii = 1, kvars
-            cvarname(ii)  = fbdata%cname(ii)
-            cunit(ii)     = fbdata%cobunit(ii)
-         END DO
-
-         ! Quality control algorithm
-         k2qc(:,:,:) = NF90_FILL_SHORT
-         DO idep = 1,kdeps
-            DO ivar = 1, kvars
-               DO iobs = 1, kobs
-                  ! 1 symbolises good for fbdata
-                  ! fbimdi symbolises that qc has not been set
-                  ! Essentially, if any fbdata flag is not an element of {1, fbimdi}
-                  ! then set the class 4 flag to bad.
-                  ! Note: fbdata%ioqc is marked good if zero.
-                  IF (((fbdata%ioqc(iobs) /= 0) .AND. &
-                            & (fbdata%ioqc(iobs) /= fbimdi)) .OR. &
-                    & ((fbdata%ipqc(iobs) /= 1) .AND. &
-                            & (fbdata%ipqc(iobs) /= fbimdi)) .OR. &
-                    & ((fbdata%idqc(idep,iobs) /= 1) .AND. &
-                            & (fbdata%idqc(idep,iobs) /= fbimdi)) .OR. &
-                    & ((fbdata%ivqc(iobs,ivar) /= 1) .AND. &
-                            & (fbdata%ivqc(iobs,ivar) /= fbimdi)) .OR. &
-                    & ((fbdata%ivlqc(idep,iobs,ivar) /= 1) .AND. &
-                            & (fbdata%ivlqc(idep,iobs,ivar) /= fbimdi)) .OR. &
-                    & ((fbdata%itqc(iobs) /= 1) .AND. &
-                            & (fbdata%itqc(iobs) /= fbimdi))) THEN
-                     ! 1 symbolises bad for class 4 file
-                     k2qc(idep, ivar, iobs) = 1
-                  ELSE
-                     ! 0 symbolises good for class 4 file
-                     k2qc(idep, ivar, iobs) = 0
-                  END IF 
-               END DO
-            END DO
-         END DO
-
-         ! Permute observation dimensions
-         pob(:,:,:) = RESHAPE(fbdata%pob, (/kdeps, kvars, kobs/), &
-                            & ORDER=(/1, 3, 2/))
-
-         ! Explicit model counterpart dimension permutation
-         ! 1,2,3,4 --> 1,4,2,3
-         pmod(:,:,:) = fbrmdi
-         ij = cl4_fcst_idx(jimatch)
-         DO ii = 1,kdeps
-            DO ik = 1, kvars
-               DO il = 1, kobs
-                  pmod(ii,ik,il) = fbdata%padd(ii,il,1,ik)
-               END DO
-            END DO
-         END DO
-
-         ! Extra fields set to missing for now
-         pextra(:,:,:) = fbrmdi
-
-         ! Lead time of class 4 file is a global parameter
-         plead = cl4_leadtime(1:cl4_fcst_len)
-
-         ! Model Julian day
-         pmodjuld = cl4_modjuld
-
-         ! Observation types
-         ctype(:) = 'X'
-         DO ii = 1,kobs
-            ctype(ii) = fbdata%cdtyp(ii)
-         END DO
-
-         ! World Meteorology Organisation codes
-         cwmo(:) = fbdata%cdwmo(:)
-
-         ! Initialise class 4 file
-         CALL ooo_wri_init(cconf, csys, ckind, cversion, ccont, cinst, cdate, &
-                         & kproc, kobs, kvars, kdeps, kfcst, &
-                         & clfilename)
-
-         ! Write standard variables
-         CALL ooo_wri_default(clfilename, kobs, kvars, kfcst, kdeps, &
-                            & ctype, cwmo, cunit, cvarname, &
-                            & plam, pphi, pdep, ptim, pob, plead, &
-                            & k2qc, pmodjuld)
-         !! Write to optional variables
-         cdtmp = cl4_vars(jimatch)
-         IF ( (TRIM(cdtmp) == "forecast") .OR. &
-              (TRIM(cdtmp) == "persistence") ) THEN
-            !! 4D variables
-            CALL ooo_wri_extra(clfilename, TRIM(cdtmp), kdeps, kfcst, &
-                            &  kvars, kobs, (/ 1,ij,1,1 /), (/ kdeps,1,kvars,kobs /), pmod)
-         ELSE
-            !! 3D variables
-            CALL ooo_wri_extra(clfilename, TRIM(cdtmp), kdeps, &
-                            &  kvars, kobs, (/ 1,1,1 /), (/ kdeps,kvars,kobs /), pmod)
-         ENDIF
-
-         DEALLOCATE(plam, pphi, ptim, pdep, plead, kqc, k2qc, &
-                  & pob, pmod, pextra, ctype, cwmo, &
-                  & cunit, cvarname)
-      END IF
-   END SUBROUTINE write_obfbdata_cl
-#endif
-
-#if defined key_offobsoper
-   SUBROUTINE locate_kind(cdfilename, ckind)
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE locate_kind  ***
-      !!
-      !! ** Purpose : Detect which kind of class 4 file is being produced.
-      !!
-      !! ** Method  : 1. Inspect cdfilename for observation kind.
-      !!----------------------------------------------------------------------
-      CHARACTER(len=*) :: cdfilename ! Feedback filename
-      CHARACTER(len=8) :: ckind
-      IF (cdfilename(1:3) == 'sst') THEN
-         ckind = 'SST'
-      ELSE IF (cdfilename(1:3) == 'sla') THEN
-         ckind = 'SLA'
-      ELSE IF (cdfilename(1:3) == 'pro') THEN
-         ckind = 'profile'
-      ELSE IF (cdfilename(1:3) == 'ena') THEN
-         ckind = 'profile'
-      ELSE IF (cdfilename(1:3) == 'sea') THEN
-         ckind = 'seaice'
-      ELSE
-         ckind = 'unknown'
-      END IF
-   END SUBROUTINE locate_kind
-#endif
+   END SUBROUTINE write_obfbdata
 
    SUBROUTINE putvaratt_obfbdata( idfile, idvar, cdlongname, cdunits, &
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAO_SRC/nemogcm.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAO_SRC/nemogcm.F90	(revision 6044)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAO_SRC/nemogcm.F90	(revision 6044)
@@ -0,0 +1,647 @@
+MODULE nemogcm
+   !!======================================================================
+   !!                       ***  MODULE nemogcm   ***
+   !! Ocean system   : NEMO GCM (ocean dynamics, on-line tracers, biochemistry and sea-ice)
+   !!======================================================================
+   !! History :  OPA  ! 1990-10  (C. Levy, G. Madec)  Original code
+   !!            7.0  ! 1991-11  (M. Imbard, C. Levy, G. Madec)
+   !!            7.1  ! 1993-03  (M. Imbard, C. Levy, G. Madec, O. Marti, M. Guyon, A. Lazar,
+   !!                             P. Delecluse, C. Perigaud, G. Caniaux, B. Colot, C. Maes) release 7.1
+   !!             -   ! 1992-06  (L.Terray)  coupling implementation
+   !!             -   ! 1993-11  (M.A. Filiberti) IGLOO sea-ice
+   !!            8.0  ! 1996-03  (M. Imbard, C. Levy, G. Madec, O. Marti, M. Guyon, A. Lazar,
+   !!                             P. Delecluse, L.Terray, M.A. Filiberti, J. Vialar, A.M. Treguier, M. Levy) release 8.0
+   !!            8.1  ! 1997-06  (M. Imbard, G. Madec)
+   !!            8.2  ! 1999-11  (M. Imbard, H. Goosse)  LIM sea-ice model
+   !!                 ! 1999-12  (V. Thierry, A-M. Treguier, M. Imbard, M-A. Foujols)  OPEN-MP
+   !!                 ! 2000-07  (J-M Molines, M. Imbard)  Open Boundary Conditions  (CLIPPER)
+   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90: Free form and modules
+   !!             -   ! 2004-06  (R. Redler, NEC CCRLE, Germany) add OASIS[3/4] coupled interfaces
+   !!             -   ! 2004-08  (C. Talandier) New trends organization
+   !!             -   ! 2005-06  (C. Ethe) Add the 1D configuration possibility
+   !!             -   ! 2005-11  (V. Garnier) Surface pressure gradient organization
+   !!             -   ! 2006-03  (L. Debreu, C. Mazauric)  Agrif implementation
+   !!             -   ! 2006-04  (G. Madec, R. Benshila)  Step reorganization
+   !!             -   ! 2007-07  (J. Chanut, A. Sellar) Unstructured open boundaries (BDY)
+   !!            3.2  ! 2009-08  (S. Masson)  open/write in the listing file in mpp
+   !!            3.3  ! 2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
+   !!             -   ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
+   !!            3.3.1! 2011-01  (A. R. Porter, STFC Daresbury) dynamical allocation
+   !!            3.4  ! 2011-11  (C. Harris) decomposition changes for running with CICE
+   !!                 ! 2012-05  (C. Calone, J. Simeon, G. Madec, C. Ethe) Add grid coarsening 
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   nemo_gcm       : solve ocean dynamics, tracer, biogeochemistry and/or sea-ice
+   !!   nemo_init      : initialization of the NEMO system
+   !!   nemo_ctl       : initialisation of the contol print
+   !!   nemo_closefile : close remaining open files
+   !!   nemo_alloc     : dynamical allocation
+   !!   nemo_partition : calculate MPP domain decomposition
+   !!   factorise      : calculate the factors of the no. of MPI processes
+   !!----------------------------------------------------------------------
+   USE step_oce        ! module used in the ocean time stepping module
+   USE domcfg          ! domain configuration               (dom_cfg routine)
+   USE mppini          ! shared/distributed memory setting (mpp_init routine)
+   USE domain          ! domain initialization             (dom_init routine)
+#if defined key_nemocice_decomp
+   USE ice_domain_size, only: nx_global, ny_global
+#endif
+   USE istate          ! initial state setting          (istate_init routine)
+   USE phycst          ! physical constant                  (par_cst routine)
+   USE diaobs          ! Observation diagnostics       (dia_obs_init routine)
+   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
+   USE step            ! NEMO time-stepping                 (stp     routine)
+   USE cpl_oasis3      ! OASIS3 coupling
+   USE lib_mpp         ! distributed memory computing
+#if defined key_iomput
+   USE xios
+#endif
+   USE lbcnfd, ONLY: isendto, nsndto, nfsloop, nfeloop ! Setup of north fold exchanges 
+
+   ! Stand Alone Observation operator modules
+   USE sao_data
+   USE sao_intp
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   nemo_gcm    ! called by model.F90
+   PUBLIC   nemo_init   ! needed by AGRIF
+   PUBLIC   nemo_alloc  ! needed by TAM
+
+   CHARACTER(lc) ::   cform_aaa="( /, 'AAAAAAAA', / ) "     ! flag for output listing
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE nemo_gcm
+         !!----------------------------------------------------------------------
+         !!                    ***  SUBROUTINE offline_obs_oper ***
+         !!
+         !! ** Purpose : To use NEMO components to interpolate model fields
+         !!              to observation space.
+         !!
+         !! ** Method : 1. Initialise NEMO
+         !!             2. Initialise offline obs_oper
+         !!             3. Cycle through match ups
+         !!             4. Write results to file
+         !!
+         !!----------------------------------------------------------------------
+         !! Initialise NEMO
+         CALL nemo_init
+         !! Initialise Stand Alone Observation operator data
+         CALL sao_data_init
+         !! Initialise obs_oper
+         CALL dia_obs_init
+         !! Interpolate to observation space
+         CALL sao_interp
+         !! Pipe to output files
+         CALL dia_obs_wri
+         !! Reset the obs_oper between
+         CALL dia_obs_dealloc
+         !! Safely stop MPI
+         IF(lk_mpp) CALL mppstop  ! end mpp communications
+   END SUBROUTINE nemo_gcm
+
+
+   SUBROUTINE nemo_init
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE nemo_init  ***
+      !!
+      !! ** Purpose :   initialization of the NEMO GCM
+      !!----------------------------------------------------------------------
+      INTEGER ::   ji            ! dummy loop indices
+      INTEGER ::   ilocal_comm   ! local integer
+      INTEGER ::   ios
+      CHARACTER(len=80), DIMENSION(16) ::   cltxt
+      !
+      NAMELIST/namctl/ ln_ctl, nn_print, nn_ictls, nn_ictle,   &
+         &             nn_isplt, nn_jsplt, nn_jctls, nn_jctle,   &
+         &             nn_bench, nn_timing
+      NAMELIST/namcfg/ cp_cfg, cp_cfz, jp_cfg, jpidta, jpjdta, jpkdta, jpiglo, jpjglo, &
+         &             jpizoom, jpjzoom, jperio, ln_use_jattr
+      !!----------------------------------------------------------------------
+      !
+      cltxt = ''
+      cxios_context = 'nemo'
+      !
+      !                             ! Open reference namelist and configuration namelist files
+      CALL ctl_opn( numnam_ref, 'namelist_ref', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. )
+      CALL ctl_opn( numnam_cfg, 'namelist_cfg', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. )
+      !
+      REWIND( numnam_ref )              ! Namelist namctl in reference namelist : Control prints & Benchmark
+      READ  ( numnam_ref, namctl, IOSTAT = ios, ERR = 901 )
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namctl in reference namelist', .TRUE. )
+
+      REWIND( numnam_cfg )              ! Namelist namctl in confguration namelist : Control prints & Benchmark
+      READ  ( numnam_cfg, namctl, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namctl in configuration namelist', .TRUE. )
+
+      !
+      REWIND( numnam_ref )              ! Namelist namcfg in reference namelist : Control prints & Benchmark
+      READ  ( numnam_ref, namcfg, IOSTAT = ios, ERR = 903 )
+903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcfg in reference namelist', .TRUE. )
+
+      REWIND( numnam_cfg )              ! Namelist namcfg in confguration namelist : Control prints & Benchmark
+      READ  ( numnam_cfg, namcfg, IOSTAT = ios, ERR = 904 )
+904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcfg in configuration namelist', .TRUE. )   
+
+! Force values for AGRIF zoom (cf. agrif_user.F90)
+#if defined key_agrif
+   IF( .NOT. Agrif_Root() ) THEN
+      jpiglo  = nbcellsx + 2 + 2*nbghostcells
+      jpjglo  = nbcellsy + 2 + 2*nbghostcells
+      jpi     = ( jpiglo-2*jpreci + (jpni-1+0) ) / jpni + 2*jpreci
+      jpj     = ( jpjglo-2*jprecj + (jpnj-1+0) ) / jpnj + 2*jprecj
+      jpidta  = jpiglo
+      jpjdta  = jpjglo
+      jpizoom = 1
+      jpjzoom = 1
+      nperio  = 0
+      jperio  = 0
+      ln_use_jattr = .false.
+   ENDIF
+#endif
+      !
+      !                             !--------------------------------------------!
+      !                             !  set communicator & select the local node  !
+      !                             !  NB: mynode also opens output.namelist.dyn !
+      !                             !      on unit number numond on first proc   !
+      !                             !--------------------------------------------!
+#if defined key_iomput
+      IF( Agrif_Root() ) THEN
+         IF( lk_oasis ) THEN
+            CALL cpl_init( "oceanx", ilocal_comm )                     ! nemo local communicator given by oasis
+            CALL xios_initialize( "not used",local_comm=ilocal_comm )    ! send nemo communicator to xios
+         ELSE
+            CALL  xios_initialize( "for_xios_mpi_id",return_comm=ilocal_comm )    ! nemo local communicator given by xios
+         ENDIF
+      ENDIF
+      ! Nodes selection (control print return in cltxt)
+      narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop, ilocal_comm )
+#else
+      IF( lk_oasis ) THEN
+         IF( Agrif_Root() ) THEN
+            CALL cpl_init( "oceanx", ilocal_comm )                      ! nemo local communicator given by oasis
+         ENDIF
+         ! Nodes selection (control print return in cltxt)
+         narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop, ilocal_comm )
+      ELSE
+         ilocal_comm = 0
+         ! Nodes selection (control print return in cltxt)
+         narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop )
+      ENDIF
+#endif
+      narea = narea + 1                                     ! mynode return the rank of proc (0 --> jpnij -1 )
+
+      lwm = (narea == 1)                                    ! control of output namelists
+      lwp = (narea == 1) .OR. ln_ctl                        ! control of all listing output print
+
+      IF(lwm) THEN
+         ! write merged namelists from earlier to output namelist now that the
+         ! file has been opened in call to mynode. nammpp has already been
+         ! written in mynode (if lk_mpp_mpi)
+         WRITE( numond, namctl )
+         WRITE( numond, namcfg )
+      ENDIF
+
+      ! If dimensions of processor grid weren't specified in the namelist file
+      ! then we calculate them here now that we have our communicator size
+      IF( (jpni < 1) .OR. (jpnj < 1) )THEN
+#if   defined key_mpp_mpi
+         IF( Agrif_Root() ) CALL nemo_partition(mppsize)
+#else
+         jpni  = 1
+         jpnj  = 1
+         jpnij = jpni*jpnj
+#endif
+      END IF
+
+      ! Calculate domain dimensions given calculated jpni and jpnj
+      ! This used to be done in par_oce.F90 when they were parameters rather
+      ! than variables
+      IF( Agrif_Root() ) THEN
+#if defined key_nemocice_decomp
+         jpi = ( nx_global+2-2*jpreci + (jpni-1) ) / jpni + 2*jpreci ! first  dim.
+         jpj = ( ny_global+2-2*jprecj + (jpnj-1) ) / jpnj + 2*jprecj ! second dim. 
+#else
+         jpi = ( jpiglo-2*jpreci + (jpni-1) ) / jpni + 2*jpreci   ! first  dim.
+         jpj = ( jpjglo-2*jprecj + (jpnj-1) ) / jpnj + 2*jprecj   ! second dim.
+#endif
+      ENDIF
+         jpk = jpkdta                                             ! third dim
+         jpim1 = jpi-1                                            ! inner domain indices
+         jpjm1 = jpj-1                                            !   "           "
+         jpkm1 = jpk-1                                            !   "           "
+         jpij  = jpi*jpj                                          !  jpi x j
+
+      IF(lwp) THEN                            ! open listing units
+         !
+         CALL ctl_opn( numout, 'ocean.output', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
+         !
+         WRITE(numout,*)
+         WRITE(numout,*) '   CNRS - NERC - Met OFFICE - MERCATOR-ocean - INGV - CMCC'
+         WRITE(numout,*) '                       NEMO team'
+         WRITE(numout,*) '            Stand Alone Observation operator'
+         WRITE(numout,*) '                  version 1.0  (2015) '
+         WRITE(numout,*)
+         WRITE(numout,*)
+         DO ji = 1, SIZE(cltxt)
+            IF( TRIM(cltxt(ji)) /= '' )   WRITE(numout,*) cltxt(ji)      ! control print of mynode
+         END DO
+         WRITE(numout,cform_aaa)                                         ! Flag AAAAAAA
+         !
+      ENDIF
+
+      ! Now we know the dimensions of the grid and numout has been set we can
+      ! allocate arrays
+      CALL nemo_alloc()
+
+      !                             !-------------------------------!
+      !                             !  NEMO general initialization  !
+      !                             !-------------------------------!
+
+      CALL nemo_ctl                          ! Control prints & Benchmark
+
+      !                                      ! Domain decomposition
+      IF( jpni*jpnj == jpnij ) THEN   ;   CALL mpp_init      ! standard cutting out
+      ELSE                            ;   CALL mpp_init2     ! eliminate land processors
+      ENDIF
+      !
+      IF( nn_timing == 1 )  CALL timing_init
+      !
+      !                                      ! General initialization
+                            CALL     phy_cst    ! Physical constants
+                            CALL     eos_init   ! Equation of state
+                            CALL     dom_cfg    ! Domain configuration
+                            CALL     dom_init   ! Domain
+
+      IF( ln_nnogather )    CALL nemo_northcomms   ! Initialise the northfold neighbour lists (must be done after the masks are defined)
+
+      IF( ln_ctl        )   CALL prt_ctl_init   ! Print control
+
+                            CALL  istate_init   ! ocean initial state (Dynamics and tracers)
+   END SUBROUTINE nemo_init
+
+
+   SUBROUTINE nemo_ctl
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE nemo_ctl  ***
+      !!
+      !! ** Purpose :   control print setting
+      !!
+      !! ** Method  : - print namctl information and check some consistencies
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) THEN                  ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'nemo_ctl: Control prints & Benchmark'
+         WRITE(numout,*) '~~~~~~~ '
+         WRITE(numout,*) '   Namelist namctl'
+         WRITE(numout,*) '      run control (for debugging)     ln_ctl     = ', ln_ctl
+         WRITE(numout,*) '      level of print                  nn_print   = ', nn_print
+         WRITE(numout,*) '      Start i indice for SUM control  nn_ictls   = ', nn_ictls
+         WRITE(numout,*) '      End i indice for SUM control    nn_ictle   = ', nn_ictle
+         WRITE(numout,*) '      Start j indice for SUM control  nn_jctls   = ', nn_jctls
+         WRITE(numout,*) '      End j indice for SUM control    nn_jctle   = ', nn_jctle
+         WRITE(numout,*) '      number of proc. following i     nn_isplt   = ', nn_isplt
+         WRITE(numout,*) '      number of proc. following j     nn_jsplt   = ', nn_jsplt
+         WRITE(numout,*) '      benchmark parameter (0/1)       nn_bench   = ', nn_bench
+         WRITE(numout,*) '      timing activated    (0/1)       nn_timing  = ', nn_timing
+      ENDIF
+      !
+      nprint    = nn_print          ! convert DOCTOR namelist names into OLD names
+      nictls    = nn_ictls
+      nictle    = nn_ictle
+      njctls    = nn_jctls
+      njctle    = nn_jctle
+      isplt     = nn_isplt
+      jsplt     = nn_jsplt
+      nbench    = nn_bench
+
+      IF(lwp) THEN                  ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'namcfg  : configuration initialization through namelist read'
+         WRITE(numout,*) '~~~~~~~ '
+         WRITE(numout,*) '   Namelist namcfg'
+         WRITE(numout,*) '      configuration name              cp_cfg      = ', TRIM(cp_cfg)
+         WRITE(numout,*) '      configuration zoom name         cp_cfz      = ', TRIM(cp_cfz)
+         WRITE(numout,*) '      configuration resolution        jp_cfg      = ', jp_cfg
+         WRITE(numout,*) '      1st lateral dimension ( >= jpi ) jpidta     = ', jpidta
+         WRITE(numout,*) '      2nd    "         "    ( >= jpj ) jpjdta     = ', jpjdta
+         WRITE(numout,*) '      3nd    "         "               jpkdta     = ', jpkdta
+         WRITE(numout,*) '      1st dimension of global domain in i jpiglo  = ', jpiglo
+         WRITE(numout,*) '      2nd    -                  -    in j jpjglo  = ', jpjglo
+         WRITE(numout,*) '      left bottom i index of the zoom (in data domain) jpizoom = ', jpizoom
+         WRITE(numout,*) '      left bottom j index of the zoom (in data domain) jpizoom = ', jpjzoom
+         WRITE(numout,*) '      lateral cond. type (between 0 and 6) jperio = ', jperio   
+         WRITE(numout,*) '      use file attribute if exists as i/p j-start ln_use_jattr = ', ln_use_jattr
+      ENDIF
+      !                             ! Parameter control
+      !
+      IF( ln_ctl ) THEN                 ! sub-domain area indices for the control prints
+         IF( lk_mpp .AND. jpnij > 1 ) THEN
+            isplt = jpni   ;   jsplt = jpnj   ;   ijsplt = jpni*jpnj   ! the domain is forced to the real split domain
+         ELSE
+            IF( isplt == 1 .AND. jsplt == 1  ) THEN
+               CALL ctl_warn( ' - isplt & jsplt are equal to 1',   &
+                  &           ' - the print control will be done over the whole domain' )
+            ENDIF
+            ijsplt = isplt * jsplt            ! total number of processors ijsplt
+         ENDIF
+         IF(lwp) WRITE(numout,*)'          - The total number of processors over which the'
+         IF(lwp) WRITE(numout,*)'            print control will be done is ijsplt : ', ijsplt
+         !
+         !                              ! indices used for the SUM control
+         IF( nictls+nictle+njctls+njctle == 0 )   THEN    ! print control done over the default area
+            lsp_area = .FALSE.
+         ELSE                                             ! print control done over a specific  area
+            lsp_area = .TRUE.
+            IF( nictls < 1 .OR. nictls > jpiglo )   THEN
+               CALL ctl_warn( '          - nictls must be 1<=nictls>=jpiglo, it is forced to 1' )
+               nictls = 1
+            ENDIF
+            IF( nictle < 1 .OR. nictle > jpiglo )   THEN
+               CALL ctl_warn( '          - nictle must be 1<=nictle>=jpiglo, it is forced to jpiglo' )
+               nictle = jpiglo
+            ENDIF
+            IF( njctls < 1 .OR. njctls > jpjglo )   THEN
+               CALL ctl_warn( '          - njctls must be 1<=njctls>=jpjglo, it is forced to 1' )
+               njctls = 1
+            ENDIF
+            IF( njctle < 1 .OR. njctle > jpjglo )   THEN
+               CALL ctl_warn( '          - njctle must be 1<=njctle>=jpjglo, it is forced to jpjglo' )
+               njctle = jpjglo
+            ENDIF
+         ENDIF
+      ENDIF
+      !
+      IF( nbench == 1 ) THEN              ! Benchmark
+         SELECT CASE ( cp_cfg )
+         CASE ( 'gyre' )   ;   CALL ctl_warn( ' The Benchmark is activated ' )
+         CASE DEFAULT      ;   CALL ctl_stop( ' The Benchmark is based on the GYRE configuration:',   &
+            &                                 ' cp_cfg = "gyre" in namelist &namcfg or set nbench = 0' )
+         END SELECT
+      ENDIF
+      !
+      IF( 1_wp /= SIGN(1._wp,-0._wp)  )   CALL ctl_stop( 'nemo_ctl: The intrinsec SIGN function follows ',  &
+         &                                               'f2003 standard. '                              ,  &
+         &                                               'Compile with key_nosignedzero enabled' )
+      !
+   END SUBROUTINE nemo_ctl
+
+
+   SUBROUTINE nemo_closefile
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE nemo_closefile  ***
+      !!
+      !! ** Purpose :   Close the files
+      !!----------------------------------------------------------------------
+      !
+      IF( lk_mpp )   CALL mppsync
+      !
+      CALL iom_close                                 ! close all input/output files managed by iom_*
+      !
+      IF( numstp          /= -1 )   CLOSE( numstp          )   ! time-step file
+      IF( numsol          /= -1 )   CLOSE( numsol          )   ! solver file
+      IF( numnam_ref      /= -1 )   CLOSE( numnam_ref      )   ! oce reference namelist
+      IF( numnam_cfg      /= -1 )   CLOSE( numnam_cfg      )   ! oce configuration namelist
+      IF( lwm.AND.numond  /= -1 )   CLOSE( numond          )   ! oce output namelist
+      IF( numnam_ice_ref  /= -1 )   CLOSE( numnam_ice_ref  )   ! ice reference namelist
+      IF( numnam_ice_cfg  /= -1 )   CLOSE( numnam_ice_cfg  )   ! ice configuration namelist
+      IF( lwm.AND.numoni  /= -1 )   CLOSE( numoni          )   ! ice output namelist
+      IF( numevo_ice      /= -1 )   CLOSE( numevo_ice      )   ! ice variables (temp. evolution)
+      IF( numout          /=  6 )   CLOSE( numout          )   ! standard model output file
+      IF( numdct_vol      /= -1 )   CLOSE( numdct_vol      )   ! volume transports
+      IF( numdct_heat     /= -1 )   CLOSE( numdct_heat     )   ! heat transports
+      IF( numdct_salt     /= -1 )   CLOSE( numdct_salt     )   ! salt transports
+
+      !
+      numout = 6                                     ! redefine numout in case it is used after this point...
+      !
+   END SUBROUTINE nemo_closefile
+
+
+   SUBROUTINE nemo_alloc
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE nemo_alloc  ***
+      !!
+      !! ** Purpose :   Allocate all the dynamic arrays of the OPA modules
+      !!
+      !! ** Method  :
+      !!----------------------------------------------------------------------
+      USE diawri    , ONLY: dia_wri_alloc
+      USE dom_oce   , ONLY: dom_oce_alloc
+      !
+      INTEGER :: ierr
+      !!----------------------------------------------------------------------
+      !
+      ierr =        oce_alloc       ()          ! ocean
+      ierr = ierr + dia_wri_alloc   ()
+      ierr = ierr + dom_oce_alloc   ()          ! ocean domain
+      !
+      IF( lk_mpp    )   CALL mpp_sum( ierr )
+      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'nemo_alloc : unable to allocate standard ocean arrays' )
+      !
+   END SUBROUTINE nemo_alloc
+
+
+   SUBROUTINE nemo_partition( num_pes )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE nemo_partition  ***
+      !!
+      !! ** Purpose :
+      !!
+      !! ** Method  :
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) :: num_pes ! The number of MPI processes we have
+      !
+      INTEGER, PARAMETER :: nfactmax = 20
+      INTEGER :: nfact ! The no. of factors returned
+      INTEGER :: ierr  ! Error flag
+      INTEGER :: ji
+      INTEGER :: idiff, mindiff, imin ! For choosing pair of factors that are closest in value
+      INTEGER, DIMENSION(nfactmax) :: ifact ! Array of factors
+      !!----------------------------------------------------------------------
+      !
+      ierr = 0
+      !
+      CALL factorise( ifact, nfactmax, nfact, num_pes, ierr )
+      !
+      IF( nfact <= 1 ) THEN
+         WRITE (numout, *) 'WARNING: factorisation of number of PEs failed'
+         WRITE (numout, *) '       : using grid of ',num_pes,' x 1'
+         jpnj = 1
+         jpni = num_pes
+      ELSE
+         ! Search through factors for the pair that are closest in value
+         mindiff = 1000000
+         imin    = 1
+         DO ji = 1, nfact-1, 2
+            idiff = ABS( ifact(ji) - ifact(ji+1) )
+            IF( idiff < mindiff ) THEN
+               mindiff = idiff
+               imin = ji
+            ENDIF
+         END DO
+         jpnj = ifact(imin)
+         jpni = ifact(imin + 1)
+      ENDIF
+      !
+      jpnij = jpni*jpnj
+      !
+   END SUBROUTINE nemo_partition
+
+
+   SUBROUTINE factorise( kfax, kmaxfax, knfax, kn, kerr )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE factorise  ***
+      !!
+      !! ** Purpose :   return the prime factors of n.
+      !!                knfax factors are returned in array kfax which is of
+      !!                maximum dimension kmaxfax.
+      !! ** Method  :
+      !!----------------------------------------------------------------------
+      INTEGER                    , INTENT(in   ) ::   kn, kmaxfax
+      INTEGER                    , INTENT(  out) ::   kerr, knfax
+      INTEGER, DIMENSION(kmaxfax), INTENT(  out) ::   kfax
+      !
+      INTEGER :: ifac, jl, inu
+      INTEGER, PARAMETER :: ntest = 14
+      INTEGER :: ilfax(ntest)
+      !
+      ! lfax contains the set of allowed factors.
+      data (ilfax(jl),jl=1,ntest) / 16384, 8192, 4096, 2048, 1024, 512, 256,  &
+         &                            128,   64,   32,   16,    8,   4,   2  /
+      !!----------------------------------------------------------------------
+
+      ! Clear the error flag and initialise output vars
+      kerr = 0
+      kfax = 1
+      knfax = 0
+
+      ! Find the factors of n.
+      IF( kn == 1 )   GOTO 20
+
+      ! nu holds the unfactorised part of the number.
+      ! knfax holds the number of factors found.
+      ! l points to the allowed factor list.
+      ! ifac holds the current factor.
+
+      inu   = kn
+      knfax = 0
+
+      DO jl = ntest, 1, -1
+         !
+         ifac = ilfax(jl)
+         IF( ifac > inu )   CYCLE
+
+         ! Test whether the factor will divide.
+
+         IF( MOD(inu,ifac) == 0 ) THEN
+            !
+            knfax = knfax + 1            ! Add the factor to the list
+            IF( knfax > kmaxfax ) THEN
+               kerr = 6
+               write (*,*) 'FACTOR: insufficient space in factor array ', knfax
+               return
+            ENDIF
+            kfax(knfax) = ifac
+            ! Store the other factor that goes with this one
+            knfax = knfax + 1
+            kfax(knfax) = inu / ifac
+            !WRITE (*,*) 'ARPDBG, factors ',knfax-1,' & ',knfax,' are ', kfax(knfax-1),' and ',kfax(knfax)
+         ENDIF
+         !
+      END DO
+
+   20 CONTINUE      ! Label 20 is the exit point from the factor search loop.
+      !
+   END SUBROUTINE factorise
+
+#if defined key_mpp_mpi
+
+   SUBROUTINE nemo_northcomms
+      !!======================================================================
+      !!                     ***  ROUTINE  nemo_northcomms  ***
+      !! nemo_northcomms    :  Setup for north fold exchanges with explicit 
+      !!                       point-to-point messaging
+      !!=====================================================================
+      !!----------------------------------------------------------------------
+      !!
+      !! ** Purpose :   Initialization of the northern neighbours lists.
+      !!----------------------------------------------------------------------
+      !!    1.0  ! 2011-10  (A. C. Coward, NOCS & J. Donners, PRACE)
+      !!    2.0  ! 2013-06 Setup avoiding MPI communication (I. Epicoco, S. Mocavero, CMCC) 
+      !!----------------------------------------------------------------------
+
+      INTEGER  ::   sxM, dxM, sxT, dxT, jn
+      INTEGER  ::   njmppmax
+
+      njmppmax = MAXVAL( njmppt )
+    
+      !initializes the north-fold communication variables
+      isendto(:) = 0
+      nsndto = 0
+
+      !if I am a process in the north
+      IF ( njmpp == njmppmax ) THEN
+          !sxM is the first point (in the global domain) needed to compute the
+          !north-fold for the current process
+          sxM = jpiglo - nimppt(narea) - nlcit(narea) + 1
+          !dxM is the last point (in the global domain) needed to compute the
+          !north-fold for the current process
+          dxM = jpiglo - nimppt(narea) + 2
+
+          !loop over the other north-fold processes to find the processes
+          !managing the points belonging to the sxT-dxT range
+  
+          DO jn = 1, jpni
+                !sxT is the first point (in the global domain) of the jn
+                !process
+                sxT = nfiimpp(jn, jpnj)
+                !dxT is the last point (in the global domain) of the jn
+                !process
+                dxT = nfiimpp(jn, jpnj) + nfilcit(jn, jpnj) - 1
+                IF ((sxM .gt. sxT) .AND. (sxM .lt. dxT)) THEN
+                   nsndto = nsndto + 1
+                     isendto(nsndto) = jn
+                ELSEIF ((sxM .le. sxT) .AND. (dxM .ge. dxT)) THEN
+                   nsndto = nsndto + 1
+                     isendto(nsndto) = jn
+                ELSEIF ((dxM .lt. dxT) .AND. (sxT .lt. dxM)) THEN
+                   nsndto = nsndto + 1
+                     isendto(nsndto) = jn
+                END IF
+          END DO
+          nfsloop = 1
+          nfeloop = nlci
+          DO jn = 2,jpni-1
+           IF(nfipproc(jn,jpnj) .eq. (narea - 1)) THEN
+              IF (nfipproc(jn - 1 ,jpnj) .eq. -1) THEN
+                 nfsloop = nldi
+              ENDIF
+              IF (nfipproc(jn + 1,jpnj) .eq. -1) THEN
+                 nfeloop = nlei
+              ENDIF
+           ENDIF
+        END DO
+
+      ENDIF
+      l_north_nogather = .TRUE.
+   END SUBROUTINE nemo_northcomms
+#else
+   SUBROUTINE nemo_northcomms      ! Dummy routine
+      WRITE(*,*) 'nemo_northcomms: You should not have seen this print! error?'
+   END SUBROUTINE nemo_northcomms
+#endif
+
+   !!======================================================================
+END MODULE nemogcm
+
+
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAO_SRC/obs_fbm.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAO_SRC/obs_fbm.F90	(revision 6044)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAO_SRC/obs_fbm.F90	(revision 6044)
@@ -0,0 +1,1 @@
+link ../OPA_SRC/OBS/obs_fbm.F90
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAO_SRC/sao.nml
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAO_SRC/sao.nml	(revision 6044)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAO_SRC/sao.nml	(revision 6044)
@@ -0,0 +1,10 @@
+!----------------------------------------------------------------------
+!       namsao Stand Alone Observation operator namelist
+!----------------------------------------------------------------------
+!   sao_files     specifies the files containing the model counterpart
+!   nn_sao_idx    specifies the index within the model file
+!   nn_sao_freq   specifies the number of timesteps between file reads
+&namsao
+   sao_files = 'fcst.0.nc'
+   nn_sao_idx = 1
+/
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAO_SRC/sao_data.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAO_SRC/sao_data.F90	(revision 6044)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAO_SRC/sao_data.F90	(revision 6044)
@@ -0,0 +1,82 @@
+MODULE sao_data
+   !! =================================================================
+   !!                    *** MODULE sao_data ***
+   !! =================================================================
+   USE par_kind, ONLY: lc
+   USE lib_mpp         ! distributed memory computing
+
+   IMPLICIT NONE
+
+   !! Public data
+
+   INTEGER, PARAMETER :: MaxNumFiles = 1000
+
+   !! Stand Alone Observation operator settings
+   CHARACTER(len=lc) :: &
+      & sao_files(MaxNumFiles)         !: model files
+   INTEGER            :: &
+      & n_files, &                     !: number of files
+      & nn_sao_idx(MaxNumFiles), &     !: time_counter indices
+      & nn_sao_freq                    !: read frequency in time steps
+CONTAINS
+   SUBROUTINE sao_data_init()
+      !!----------------------------------------------------------------------
+      !!                    ***  SUBROUTINE sao_data_init ***
+      !!
+      !! ** Purpose : To read namelists and initialise offline_oper run.
+      !!
+      !!----------------------------------------------------------------------
+      USE in_out_manager
+      INTEGER            :: &
+         & jf                           !: file dummy loop index
+      LOGICAL :: lmask(MaxNumFiles)     !: Logical mask used for counting
+      INTEGER :: ios
+
+      ! Standard offline obs_oper information
+      NAMELIST/namsao/sao_files, nn_sao_idx, nn_sao_freq
+
+      ! Standard offline obs_oper initialisation
+      n_files = 0                   !: number of files to cycle through
+      sao_files(:) = ''             !: list of files to read in
+      nn_sao_idx(:) = 0             !: list of indices inside each file
+      nn_sao_freq = -1              !: input frequency in time steps
+
+      ! Standard offline obs_oper settings
+      REWIND( numnam_ref )              ! Namelist namctl in reference namelist : Control prints & Benchmark
+      READ  ( numnam_ref, namsao, IOSTAT = ios, ERR = 901 )
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsao in reference namelist', .TRUE. )
+
+      REWIND( numnam_cfg )              ! Namelist namctl in confguration namelist : Control prints & Benchmark
+      READ  ( numnam_cfg, namsao, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsao in configuration namelist', .TRUE. )
+
+
+      ! count input files
+      lmask(:) = .FALSE.
+      WHERE (sao_files(:) /= '') lmask(:) = .TRUE.
+      n_files = COUNT(lmask)
+
+      !! Initialise sub obs window frequency
+      IF (nn_sao_freq == -1) THEN
+         !! Run length
+         nn_sao_freq = nitend - nit000 + 1
+      ENDIF
+
+      !! Print summary of settings
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'offline obs_oper : Initialization'
+         WRITE(numout,*) '~~~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namsao : set stand alone obs_oper parameters'
+         DO jf = 1, n_files
+            WRITE(numout,'(1X,2A)') '   Input forecast file name          forecastfile = ', &
+               TRIM(sao_files(jf))
+            WRITE(numout,*) '   Input forecast file index        forecastindex = ', &
+               nn_sao_idx(jf)
+         END DO
+      END IF
+
+   END SUBROUTINE sao_data_init
+
+END MODULE sao_data
+
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAO_SRC/sao_intp.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAO_SRC/sao_intp.F90	(revision 6044)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAO_SRC/sao_intp.F90	(revision 6044)
@@ -0,0 +1,48 @@
+MODULE sao_intp
+   !!======================================================================
+   !!                    ***  MODULE sao_intp ***
+   !! ** Purpose : Run NEMO observation operator in offline mode
+   !!======================================================================
+   !! NEMO modules
+   USE in_out_manager
+   USE diaobs
+   !! Stand Alone Observation operator modules
+   USE sao_read
+   USE sao_data
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC sao_interp
+
+   CONTAINS
+
+      SUBROUTINE sao_interp
+         !!----------------------------------------------------------------------
+         !!                    ***  SUBROUTINE sao_interp ***
+         !!
+         !! ** Purpose : To interpolate the model as if it were running online.
+         !!
+         !! ** Method : 1. Populate model counterparts
+         !!             2. Call dia_obs at appropriate time steps
+         !!----------------------------------------------------------------------
+         INTEGER :: &
+            & istp, & ! time step index
+            & ifile   ! file index
+         istp = nit000 - 1
+         nstop = 0
+         ifile = 1
+         CALL sao_rea_dri(ifile)
+         DO WHILE ( istp <= nitend .AND. nstop == 0 )
+            IF (ifile <= n_files + 1) THEN
+               IF ( MOD(istp, nn_sao_freq) == nit000 ) THEN
+                  CALL sao_rea_dri(ifile)
+                  ifile = ifile + 1
+               ENDIF
+               CALL dia_obs(istp)
+            ENDIF
+            istp = istp + 1
+         END DO
+      END SUBROUTINE sao_interp
+
+END MODULE sao_intp
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAO_SRC/sao_read.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAO_SRC/sao_read.F90	(revision 6044)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/SAO_SRC/sao_read.F90	(revision 6044)
@@ -0,0 +1,172 @@
+MODULE sao_read
+   !!==================================================================
+   !!                    *** MODULE sao_read ***
+   !! Read routines : I/O for Stand Alone Observation operator
+   !!==================================================================
+   USE mppini
+   USE lib_mpp
+   USE in_out_manager
+   USE par_kind, ONLY: lc
+   USE netcdf
+   USE oce,     ONLY: tsn, sshn
+   USE dom_oce, ONLY: nlci, nlcj, nimpp, njmpp, tmask
+   USE par_oce, ONLY: jpi, jpj, jpk
+   USE obs_fbm, ONLY: fbimdi, fbrmdi, fbsp, fbdp
+   USE sao_data
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC sao_rea_dri
+
+CONTAINS
+   SUBROUTINE sao_rea_dri(kfile)
+      !!------------------------------------------------------------------------
+      !!             *** sao_rea_dri ***
+      !!
+      !! Purpose : To choose appropriate read method
+      !! Method  :
+      !!
+      !! Author  : A. Ryan Oct 2013
+      !!
+      !!------------------------------------------------------------------------
+      INTEGER, INTENT(IN) :: &
+              & kfile         !: File number
+      CHARACTER(len=lc) :: &
+              & cdfilename    !: File name
+      INTEGER :: &
+              & kindex        !: File index to read
+
+      cdfilename = TRIM(sao_files(kfile))
+      kindex = nn_sao_idx(kfile)
+      CALL sao_read_file(TRIM(cdfilename), kindex)
+
+   END SUBROUTINE sao_rea_dri
+
+   SUBROUTINE sao_read_file(filename, ifcst)
+      !!------------------------------------------------------------------------
+      !!             *** sao_read_file ***
+      !!
+      !! Purpose : To fill tn and sn with dailymean field from netcdf files
+      !! Method  : Use subdomain indices to create start and count matrices
+      !!           for netcdf read.
+      !!
+      !! Author  : A. Ryan Oct 2010
+      !!------------------------------------------------------------------------
+
+      INTEGER,          INTENT(IN) :: ifcst
+      CHARACTER(len=*), INTENT(IN) :: filename
+      INTEGER                      :: ncid, &
+                                    & varid,&
+                                    & istat,&
+                                    & ntimes,&
+                                    & tdim, &
+                                    & xdim, &
+                                    & ydim, &
+                                    & zdim
+      INTEGER                      :: ii, ij, ik
+      INTEGER, DIMENSION(4)        :: start_n, &
+                                    & count_n
+      INTEGER, DIMENSION(3)        :: start_s, &
+                                    & count_s
+      REAL(fbdp), DIMENSION(:,:,:),ALLOCATABLE :: temp_tn, &
+                                              & temp_sn
+      REAL(fbdp), DIMENSION(:,:),  ALLOCATABLE :: temp_sshn
+      REAL(fbdp)                     :: fill_val
+
+      ! DEBUG
+      INTEGER :: istage
+
+      IF (TRIM(filename) == 'nofile') THEN
+         tsn(:,:,:,:) = fbrmdi
+         sshn(:,:) = fbrmdi
+      ELSE
+         WRITE(numout,*) "Opening :", TRIM(filename)
+         ! Open Netcdf file to find dimension id
+         istat = nf90_open(path=TRIM(filename), mode=nf90_nowrite, ncid=ncid)
+         IF ( istat /= nf90_noerr ) THEN
+             WRITE(numout,*) "WARNING: Could not open ", trim(filename)
+             WRITE(numout,*) "ERROR: ", nf90_strerror(istat)
+         ENDIF
+         istat = nf90_inq_dimid(ncid,'x',xdim)
+         istat = nf90_inq_dimid(ncid,'y',ydim)
+         istat = nf90_inq_dimid(ncid,'deptht',zdim)
+         istat = nf90_inq_dimid(ncid,'time_counter',tdim)
+         istat = nf90_inquire_dimension(ncid, tdim, len=ntimes)
+         IF (ifcst .LE. ntimes) THEN
+            ! Allocate temporary temperature array
+            ALLOCATE(temp_tn(nlci,nlcj,jpk))
+            ALLOCATE(temp_sn(nlci,nlcj,jpk))
+            ALLOCATE(temp_sshn(nlci,nlcj))
+
+            ! Set temp_tn, temp_sn to 0.
+            temp_tn(:,:,:) = fbrmdi
+            temp_sn(:,:,:) = fbrmdi
+            temp_sshn(:,:) = fbrmdi
+
+            ! Create start and count arrays
+            start_n = (/ nimpp, njmpp, 1,   ifcst /)
+            count_n = (/ nlci,  nlcj,  jpk, 1     /)
+            start_s = (/ nimpp, njmpp, ifcst /)
+            count_s = (/ nlci,  nlcj,  1     /)
+
+            ! Read information into temporary arrays
+            ! retrieve varid and read in temperature
+            istat = nf90_inq_varid(ncid,'votemper',varid)
+            istat = nf90_get_att(ncid, varid, '_FillValue', fill_val)
+            istat = nf90_get_var(ncid, varid, temp_tn, start_n, count_n)
+            WHERE(temp_tn(:,:,:) == fill_val) temp_tn(:,:,:) = fbrmdi
+
+            ! retrieve varid and read in salinity
+            istat = nf90_inq_varid(ncid,'vosaline',varid)
+            istat = nf90_get_att(ncid, varid, '_FillValue', fill_val)
+            istat = nf90_get_var(ncid, varid, temp_sn, start_n, count_n)
+            WHERE(temp_sn(:,:,:) == fill_val) temp_sn(:,:,:) = fbrmdi
+
+            ! retrieve varid and read in SSH
+            istat = nf90_inq_varid(ncid,'sossheig',varid)
+            IF (istat /= nf90_noerr) THEN
+               ! Altimeter bias
+               istat = nf90_inq_varid(ncid,'altbias',varid)
+            END IF
+
+            istat = nf90_get_att(ncid, varid, '_FillValue', fill_val)
+            istat = nf90_get_var(ncid, varid, temp_sshn, start_s, count_s)
+            WHERE(temp_sshn(:,:) == fill_val) temp_sshn(:,:) = fbrmdi
+
+            ! Initialise tsn, sshn to fbrmdi
+            tsn(:,:,:,:) = fbrmdi
+            sshn(:,:) = fbrmdi
+
+            ! Mask out missing data index
+            tsn(1:nlci,1:nlcj,1:jpk,1) = temp_tn(:,:,:) * tmask(1:nlci,1:nlcj,1:jpk)
+            tsn(1:nlci,1:nlcj,1:jpk,2) = temp_sn(:,:,:) * tmask(1:nlci,1:nlcj,1:jpk)
+            sshn(1:nlci,1:nlcj)        = temp_sshn(:,:) * tmask(1:nlci,1:nlcj,1)
+
+            ! Remove halo from tmask, tsn, sshn to prevent double obs counting
+            IF (jpi > nlci) THEN
+                tmask(nlci+1:,:,:) = 0
+                tsn(nlci+1:,:,:,1) = 0
+                tsn(nlci+1:,:,:,2) = 0
+                sshn(nlci+1:,:) = 0
+            END IF
+            IF (jpj > nlcj) THEN
+                tmask(:,nlcj+1:,:) = 0
+                tsn(:,nlcj+1:,:,1) = 0
+                tsn(:,nlcj+1:,:,2) = 0
+                sshn(:,nlcj+1:) = 0
+            END IF
+
+            ! Deallocate arrays
+            DEALLOCATE(temp_tn, temp_sn, temp_sshn)
+         ELSE
+            ! Mark all as missing data
+            tsn(:,:,:,:) = fbrmdi
+            sshn(:,:) = fbrmdi
+         ENDIF
+         ! Close netcdf file
+         WRITE(numout,*) "Closing :", TRIM(filename)
+         istat = nf90_close(ncid)
+      END IF
+   END SUBROUTINE sao_read_file
+END MODULE sao_read
