Index: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/SAO_SRC/nemogcm.F90
===================================================================
--- branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/SAO_SRC/nemogcm.F90	(revision 4847)
+++ branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/SAO_SRC/nemogcm.F90	(revision 4848)
@@ -66,5 +66,4 @@
    ! Stand Alone Observation operator modules
    USE ooo_data
-   USE ooo_read
    USE sao_intp
 
Index: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/SAO_SRC/ooo_read.F90
===================================================================
--- branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/SAO_SRC/ooo_read.F90	(revision 4847)
+++ 	(revision )
@@ -1,347 +1,0 @@
-
-MODULE ooo_read
-   !!==================================================================
-   !!                    *** MODULE ooo_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 ooo_data
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC ooo_rea_dri
-
-CONTAINS
-   SUBROUTINE ooo_rea_dri(kfile)
-      !!------------------------------------------------------------------------
-      !!             *** ooo_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
-              & cmatchname    !: Match name
-      INTEGER :: &
-              & kindex       !: File index to read
- 
-      !! Filename, index and match-up kind
-      cdfilename = TRIM(sao_files(kfile))
-      cmatchname = TRIM(cl4_vars(kfile))
-      kindex = nn_sao_idx(kfile)
-
-      !! Update model fields
-      !! Class 4 variables: forecast, persistence,
-      !!                    nrt_analysis, best_estimate
-      !! Feedback variables: empty string
-      IF ( (TRIM(cmatchname) == 'forecast') .OR. &
-         & (TRIM(cmatchname) == 'persistence') .OR. &
-         & (TRIM(cmatchname) == 'nrt_analysis') .OR. &
-         & (TRIM(cmatchname) == 'best_estimate').OR. &
-         & (TRIM(cmatchname) == '') ) THEN
-         CALL ooo_read_file(TRIM(cdfilename), kindex)
-         CALL ooo_read_juld(TRIM(cdfilename), kindex, cl4_modjuld)
-      ELSE IF (TRIM(cmatchname) == 'climatology') THEN
-         WRITE(numout,*) 'Interpolating climatologies'
-      ELSE IF (TRIM(cmatchname) == 'altimeter') THEN
-         CALL ooo_read_altbias(TRIM(cdfilename))
-         CALL ooo_read_juld(TRIM(cdfilename), kindex, cl4_modjuld)
-      END IF
-
-   END SUBROUTINE ooo_rea_dri
-
-   SUBROUTINE ooo_read_altbias(filename)
-      !!------------------------------------------------------------------------
-      !!                      *** ooo_read_altbias ***
-      !!
-      !! Purpose : To read altimeter bias and set tn,sn to missing values
-      !! Method  : Use subdomain indices to create start and count matrices
-      !!           for netcdf read.
-      !!
-      !! Author  : A. Ryan Sept 2012
-      !!------------------------------------------------------------------------
-      CHARACTER(len=*), INTENT(IN) :: filename
-      INTEGER                      :: ncid, &
-                                    & varid,&
-                                    & istat,&
-                                    & ntimes,&
-                                    & tdim, &
-                                    & xdim, &
-                                    & ydim, &
-                                    & zdim
-      INTEGER                      :: ii, ij, ik
-      INTEGER, DIMENSION(3)        :: start_s, &
-                                    & count_s
-      REAL(fbdp), DIMENSION(:,:),  ALLOCATABLE :: temp_sshn
-      REAL(fbdp)                     :: fill_val
-
-      IF (TRIM(filename) == 'nofile') THEN
-         tsn(:,:,:,:) = fbrmdi
-         sshn(:,:) = fbrmdi 
-      ELSE
-         ! Open Netcdf file to find dimension id
-         istat = nf90_open(trim(filename),nf90_nowrite,ncid)
-         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',tdim)
-         istat = nf90_inquire_dimension(ncid, tdim, len=ntimes)
-         ! Allocate temporary temperature array
-         ALLOCATE(temp_sshn(nlci,nlcj))
-         ! Create start and count arrays
-         start_s = (/ nimpp, njmpp, 1 /)
-         count_s = (/ nlci,  nlcj,  1 /)
-          
-         ! Altimeter bias
-         istat = nf90_inq_varid(ncid,'altbias',varid)
-         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 
-
-         ! Fill sshn with altimeter bias 
-         sshn(1:nlci,1:nlcj) = temp_sshn(:,:) * tmask(1:nlci,1:nlcj,1)
-
-         ! Remove halo from tmask, sshn to prevent double obs counting
-         IF (jpi > nlci) THEN
-             tmask(nlci+1:,:,:) = 0
-             sshn(nlci+1:,:) = 0
-         END IF
-         IF (jpj > nlcj) THEN
-             tmask(:,nlcj+1:,:) = 0
-             sshn(:,nlcj+1:) = 0
-         END IF
-
-         ! Deallocate arrays
-         DEALLOCATE(temp_sshn)
-
-         ! Close netcdf file
-         istat = nf90_close(ncid)
-      END IF
-   
-   END SUBROUTINE ooo_read_altbias
-
-   SUBROUTINE ooo_read_file(filename, ifcst)
-      !!------------------------------------------------------------------------
-      !!             *** ooo_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 ooo_read_file
-
-   SUBROUTINE ooo_read_juld(filename, ifcst, julian)
-      USE calendar
-      !!--------------------------------------------------------------------
-      !!                 *** ooo_read_juld ***
-      !!
-      !!   Purpose : To read model julian day information from file
-      !!   Author  : A. Ryan Nov 2010
-      !!--------------------------------------------------------------------
-
-      !! Routine arguments
-      CHARACTER(len=*), INTENT(IN)  :: filename
-      INTEGER,          INTENT(IN)  :: ifcst
-      REAL,             INTENT(OUT) :: julian    !: Julian day
-
-      !! Local variables
-      INTEGER :: year,  &   !: Date information
-               & month, &
-               & day,   &
-               & hour,  &
-               & minute,&
-               & second
-      INTEGER :: istat, &   !: Netcdf variables
-               & ncid,  &
-               & dimid, &
-               & varid, &
-               & ntime      
-      REAL,DIMENSION(:),ALLOCATABLE :: r_sec     !: Remainder seconds
-      CHARACTER(len=120) :: time_str  !: time string
-
-      time_str=''
-      !! Read in time_counter and remainder seconds
-      istat = nf90_open(trim(filename),nf90_nowrite,ncid)
-      istat = nf90_inq_dimid(ncid,'time_counter',dimid)
-      IF (istat /= nf90_noerr) THEN
-          istat = nf90_inq_dimid(ncid,'time',dimid)
-      ENDIF
-      istat = nf90_inquire_dimension(ncid,dimid,len=ntime)
-      istat = nf90_inq_varid(ncid,'time_counter',varid)
-      IF (istat /= nf90_noerr) THEN
-          istat = nf90_inq_dimid(ncid,'time',dimid)
-      ENDIF
-      istat = nf90_get_att(ncid,varid,'units',time_str)
-      ALLOCATE(r_sec(ntime))
-      istat = nf90_get_var(ncid,varid, r_sec)
-      istat = nf90_close(ncid)
-
-      !! Fill yyyy-mm-dd hh:mm:ss
-      !! format(('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2))
-      100 format((14x, I4.4,1x,I2.2,1x,I2.2,1x,I2.2,1x,I2.2,1x,I2.2))
-      READ( time_str, 100 ) year, month, day, hour, minute, second
-
-      CALL ymds2ju(year, month, day, r_sec(ifcst), julian)
-
-      !! To take a comment from the ymds2ju subroutine
-
-   !- In the case of the Gregorian calendar we have chosen to use
-   !- the Lilian day numbers. This is the day counter which starts
-   !- on the 15th October 1582.
-   !- This is the day at which Pope Gregory XIII introduced the
-   !- Gregorian calendar.
-   !- Compared to the true Julian calendar, which starts some
-   !- 7980 years ago, the Lilian days are smaler and are dealt with
-   !- easily on 32 bit machines. With the true Julian days you can only
-   !- the fraction of the day in the real part to a precision of
-   !- a 1/4 of a day with 32 bits.
-      
-      !! The obs operator routines prefer to calculate Julian days since 
-      !! 01/01/1950 00:00:00
-      !! In order to convert to the 1950 version we must adjust by the number
-      !! of days between 15th October 1582 and 1st Jan 1950
-
-      julian = julian - 134123.
-      
-      DEALLOCATE(r_sec)
-      
-   END SUBROUTINE ooo_read_juld
-
-END MODULE ooo_read 
-
Index: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/SAO_SRC/sao_intp.F90
===================================================================
--- branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/SAO_SRC/sao_intp.F90	(revision 4847)
+++ branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/SAO_SRC/sao_intp.F90	(revision 4848)
@@ -8,5 +8,5 @@
    USE diaobs
    !! Stand Alone Observation operator modules
-   USE ooo_read
+   USE sao_read
    USE ooo_data
 
@@ -35,5 +35,5 @@
                IF ( MOD(istp, nn_sao_freq) == nit000 ) THEN
                   !! Read next model counterpart
-                  CALL ooo_rea_dri(jifile)
+                  CALL sao_rea_dri(jifile)
                   jifile = jifile + 1
                ENDIF
Index: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/SAO_SRC/sao_read.F90
===================================================================
--- branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/SAO_SRC/sao_read.F90	(revision 4848)
+++ branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/SAO_SRC/sao_read.F90	(revision 4848)
@@ -0,0 +1,347 @@
+
+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 ooo_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
+              & cmatchname    !: Match name
+      INTEGER :: &
+              & kindex       !: File index to read
+ 
+      !! Filename, index and match-up kind
+      cdfilename = TRIM(sao_files(kfile))
+      cmatchname = TRIM(cl4_vars(kfile))
+      kindex = nn_sao_idx(kfile)
+
+      !! Update model fields
+      !! Class 4 variables: forecast, persistence,
+      !!                    nrt_analysis, best_estimate
+      !! Feedback variables: empty string
+      IF ( (TRIM(cmatchname) == 'forecast') .OR. &
+         & (TRIM(cmatchname) == 'persistence') .OR. &
+         & (TRIM(cmatchname) == 'nrt_analysis') .OR. &
+         & (TRIM(cmatchname) == 'best_estimate').OR. &
+         & (TRIM(cmatchname) == '') ) THEN
+         CALL sao_read_file(TRIM(cdfilename), kindex)
+         CALL sao_read_juld(TRIM(cdfilename), kindex, cl4_modjuld)
+      ELSE IF (TRIM(cmatchname) == 'climatology') THEN
+         WRITE(numout,*) 'Interpolating climatologies'
+      ELSE IF (TRIM(cmatchname) == 'altimeter') THEN
+         CALL sao_read_altbias(TRIM(cdfilename))
+         CALL sao_read_juld(TRIM(cdfilename), kindex, cl4_modjuld)
+      END IF
+
+   END SUBROUTINE sao_rea_dri
+
+   SUBROUTINE sao_read_altbias(filename)
+      !!------------------------------------------------------------------------
+      !!                      *** sao_read_altbias ***
+      !!
+      !! Purpose : To read altimeter bias and set tn,sn to missing values
+      !! Method  : Use subdomain indices to create start and count matrices
+      !!           for netcdf read.
+      !!
+      !! Author  : A. Ryan Sept 2012
+      !!------------------------------------------------------------------------
+      CHARACTER(len=*), INTENT(IN) :: filename
+      INTEGER                      :: ncid, &
+                                    & varid,&
+                                    & istat,&
+                                    & ntimes,&
+                                    & tdim, &
+                                    & xdim, &
+                                    & ydim, &
+                                    & zdim
+      INTEGER                      :: ii, ij, ik
+      INTEGER, DIMENSION(3)        :: start_s, &
+                                    & count_s
+      REAL(fbdp), DIMENSION(:,:),  ALLOCATABLE :: temp_sshn
+      REAL(fbdp)                     :: fill_val
+
+      IF (TRIM(filename) == 'nofile') THEN
+         tsn(:,:,:,:) = fbrmdi
+         sshn(:,:) = fbrmdi 
+      ELSE
+         ! Open Netcdf file to find dimension id
+         istat = nf90_open(trim(filename),nf90_nowrite,ncid)
+         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',tdim)
+         istat = nf90_inquire_dimension(ncid, tdim, len=ntimes)
+         ! Allocate temporary temperature array
+         ALLOCATE(temp_sshn(nlci,nlcj))
+         ! Create start and count arrays
+         start_s = (/ nimpp, njmpp, 1 /)
+         count_s = (/ nlci,  nlcj,  1 /)
+          
+         ! Altimeter bias
+         istat = nf90_inq_varid(ncid,'altbias',varid)
+         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 
+
+         ! Fill sshn with altimeter bias 
+         sshn(1:nlci,1:nlcj) = temp_sshn(:,:) * tmask(1:nlci,1:nlcj,1)
+
+         ! Remove halo from tmask, sshn to prevent double obs counting
+         IF (jpi > nlci) THEN
+             tmask(nlci+1:,:,:) = 0
+             sshn(nlci+1:,:) = 0
+         END IF
+         IF (jpj > nlcj) THEN
+             tmask(:,nlcj+1:,:) = 0
+             sshn(:,nlcj+1:) = 0
+         END IF
+
+         ! Deallocate arrays
+         DEALLOCATE(temp_sshn)
+
+         ! Close netcdf file
+         istat = nf90_close(ncid)
+      END IF
+   
+   END SUBROUTINE sao_read_altbias
+
+   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
+
+   SUBROUTINE sao_read_juld(filename, ifcst, julian)
+      USE calendar
+      !!--------------------------------------------------------------------
+      !!                 *** sao_read_juld ***
+      !!
+      !!   Purpose : To read model julian day information from file
+      !!   Author  : A. Ryan Nov 2010
+      !!--------------------------------------------------------------------
+
+      !! Routine arguments
+      CHARACTER(len=*), INTENT(IN)  :: filename
+      INTEGER,          INTENT(IN)  :: ifcst
+      REAL,             INTENT(OUT) :: julian    !: Julian day
+
+      !! Local variables
+      INTEGER :: year,  &   !: Date information
+               & month, &
+               & day,   &
+               & hour,  &
+               & minute,&
+               & second
+      INTEGER :: istat, &   !: Netcdf variables
+               & ncid,  &
+               & dimid, &
+               & varid, &
+               & ntime      
+      REAL,DIMENSION(:),ALLOCATABLE :: r_sec     !: Remainder seconds
+      CHARACTER(len=120) :: time_str  !: time string
+
+      time_str=''
+      !! Read in time_counter and remainder seconds
+      istat = nf90_open(trim(filename),nf90_nowrite,ncid)
+      istat = nf90_inq_dimid(ncid,'time_counter',dimid)
+      IF (istat /= nf90_noerr) THEN
+          istat = nf90_inq_dimid(ncid,'time',dimid)
+      ENDIF
+      istat = nf90_inquire_dimension(ncid,dimid,len=ntime)
+      istat = nf90_inq_varid(ncid,'time_counter',varid)
+      IF (istat /= nf90_noerr) THEN
+          istat = nf90_inq_dimid(ncid,'time',dimid)
+      ENDIF
+      istat = nf90_get_att(ncid,varid,'units',time_str)
+      ALLOCATE(r_sec(ntime))
+      istat = nf90_get_var(ncid,varid, r_sec)
+      istat = nf90_close(ncid)
+
+      !! Fill yyyy-mm-dd hh:mm:ss
+      !! format(('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2))
+      100 format((14x, I4.4,1x,I2.2,1x,I2.2,1x,I2.2,1x,I2.2,1x,I2.2))
+      READ( time_str, 100 ) year, month, day, hour, minute, second
+
+      CALL ymds2ju(year, month, day, r_sec(ifcst), julian)
+
+      !! To take a comment from the ymds2ju subroutine
+
+   !- In the case of the Gregorian calendar we have chosen to use
+   !- the Lilian day numbers. This is the day counter which starts
+   !- on the 15th October 1582.
+   !- This is the day at which Pope Gregory XIII introduced the
+   !- Gregorian calendar.
+   !- Compared to the true Julian calendar, which starts some
+   !- 7980 years ago, the Lilian days are smaler and are dealt with
+   !- easily on 32 bit machines. With the true Julian days you can only
+   !- the fraction of the day in the real part to a precision of
+   !- a 1/4 of a day with 32 bits.
+      
+      !! The obs operator routines prefer to calculate Julian days since 
+      !! 01/01/1950 00:00:00
+      !! In order to convert to the 1950 version we must adjust by the number
+      !! of days between 15th October 1582 and 1st Jan 1950
+
+      julian = julian - 134123.
+      
+      DEALLOCATE(r_sec)
+      
+   END SUBROUTINE sao_read_juld
+
+END MODULE sao_read 
+
