Index: branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
===================================================================
--- branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 6854)
+++ branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 6855)
@@ -31,4 +31,5 @@
    USE obs_read_vel             ! Reading and allocation of velocity component observations
    USE obs_read_logchl          ! Reading and allocation of logchl observations
+   USE obs_read_spm             ! Reading and allocation of spm observations
    USE obs_prep                 ! Preparation of obs. (grid search etc).
    USE obs_oper                 ! Observation operators
@@ -43,4 +44,5 @@
    USE obs_seaice               ! Sea Ice data storage
    USE obs_logchl               ! logchl data storage
+   USE obs_spm                  ! spm data storage
    USE obs_types                ! Definitions for observation types
    USE mpp_map                  ! MPP mapping
@@ -86,4 +88,6 @@
    LOGICAL, PUBLIC :: ln_logchl      !: Logical switch for log10(chlorophyll)
    LOGICAL, PUBLIC :: ln_logchlfb    !: Logical switch for logchl from feedback files
+   LOGICAL, PUBLIC :: ln_spm         !: Logical switch for spm
+   LOGICAL, PUBLIC :: ln_spmfb       !: Logical switch for spm from feedback files
    LOGICAL, PUBLIC :: ln_ssh         !: Logical switch for sea surface height
    LOGICAL, PUBLIC :: ln_sss         !: Logical switch for sea surface salinity
@@ -171,4 +175,6 @@
       CHARACTER(len=128) :: logchlfiles(MaxNumFiles)
       CHARACTER(len=128) :: logchlfbfiles(MaxNumFiles)
+      CHARACTER(len=128) :: spmfiles(MaxNumFiles)
+      CHARACTER(len=128) :: spmfbfiles(MaxNumFiles)
       CHARACTER(LEN=128) :: reysstname
       CHARACTER(LEN=12)  :: reysstfmt
@@ -197,4 +203,6 @@
          &            ln_logchl, ln_logchlfb,                         &
          &            logchlfiles, logchlfbfiles,                     &
+         &            ln_spm, ln_spmfb,                               &
+         &            spmfiles, spmfbfiles,                           &
          &            ln_profb_enatim, ln_ignmis, ln_cl4,             &
          &            ln_sstbias, sstbias_files
@@ -220,4 +228,6 @@
       INTEGER :: jnumlogchl
       INTEGER :: jnumlogchlfb
+      INTEGER :: jnumspm
+      INTEGER :: jnumspmfb
       INTEGER :: ji
       INTEGER :: jset
@@ -231,4 +241,6 @@
       ln_logchl   = .FALSE.
       ln_logchlfb = .FALSE.
+      ln_spm      = .FALSE.
+      ln_spmfb    = .FALSE.
       
       !Initalise all values in namelist arrays
@@ -253,4 +265,6 @@
       logchlfiles(:) = ''
       logchlfbfiles(:) = ''
+      spmfiles(:) = ''
+      spmfbfiles(:) = ''
       sstbias_files(:) = ''
       endailyavtypes(:) = -1
@@ -361,4 +375,14 @@
          WHERE (logchlfbfiles(:) /= '') lmask(:) = .TRUE.
          jnumlogchlfb = COUNT(lmask)
+      ENDIF
+      IF (ln_spm) THEN
+         lmask(:) = .FALSE.
+         WHERE (spmfiles(:) /= '') lmask(:) = .TRUE.
+         jnumspm = COUNT(lmask)
+      ENDIF
+      IF (ln_spmfb) THEN
+         lmask(:) = .FALSE.
+         WHERE (spmfbfiles(:) /= '') lmask(:) = .TRUE.
+         jnumspmfb = COUNT(lmask)
       ENDIF
       
@@ -394,4 +418,6 @@
          WRITE(numout,*) '             Logical switch for logchl observations          ln_logchl = ', ln_logchl
          WRITE(numout,*) '             Logical switch for feedback logchl data       ln_logchlfb = ', ln_logchlfb
+         WRITE(numout,*) '             Logical switch for spm observations                ln_spm = ', ln_spm
+         WRITE(numout,*) '             Logical switch for feedback spm data             ln_spmfb = ', ln_spmfb
          WRITE(numout,*) '             Global distribtion of observations         ln_grid_global = ',ln_grid_global
          WRITE(numout,*) &
@@ -502,4 +528,16 @@
                WRITE(numout,'(1X,2A)') '        Feedback logchl input observation file name  logchlfbfiles = ', &
                   TRIM(logchlfbfiles(ji))
+            END DO
+         ENDIF
+         IF (ln_spm) THEN
+            DO ji = 1, jnumspm
+               WRITE(numout,'(1X,2A)') '             spm input observation file name              spmfiles = ', &
+                  TRIM(spmfiles(ji))
+            END DO
+         ENDIF
+         IF (ln_spmfb) THEN
+            DO ji = 1, jnumspmfb
+               WRITE(numout,'(1X,2A)') '             Feedback spm input observation file name   spmfbfiles = ', &
+                  TRIM(spmfbfiles(ji))
             END DO
          ENDIF
@@ -538,9 +576,10 @@
          & ( .NOT. ln_vel3d ).AND.                                         &
          & ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. &
-         & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ).AND.( .NOT. ln_logchl ) ) THEN
+         & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ).AND.( .NOT. ln_logchl ).AND. &
+         & ( .NOT. ln_spm ) ) THEN
          IF(lwp) WRITE(numout,cform_war)
          IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', &
             &                    ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d,', &
-            &                    ' ln_logchl are all set to .FALSE.'
+            &                    ' ln_logchl, ln_spm are all set to .FALSE.'
          nwarn = nwarn + 1
       ENDIF
@@ -1095,4 +1134,59 @@
  
       ENDIF
+
+      !  - spm
+      
+      IF ( ln_spm ) THEN
+
+         ! Set the number of variables for spm to 1
+         nspmvars = 1
+
+         ! Set the number of extra variables for spm to 0
+         nspmextr = 0
+         
+         IF ( ln_spmfb ) THEN
+            nspmsets = jnumspmfb
+         ELSE
+            nspmsets = 1
+         ENDIF
+
+         ALLOCATE(spmdata(nspmsets))
+         ALLOCATE(spmdatqc(nspmsets))
+         spmdata(:)%nsurf=0
+         spmdatqc(:)%nsurf=0
+
+         nspmsets = 0
+
+         IF ( ln_spmfb ) THEN             ! Feedback file format
+
+            DO jset = 1, jnumspmfb
+            
+               nspmsets = nspmsets + 1
+
+               CALL obs_rea_spm( 0, spmdata(nspmsets), 1, &
+                  &                 spmfbfiles(jset:jset), &
+                  &                 nspmvars, nspmextr, nitend-nit000+2, &
+                  &                 dobsini, dobsend, ln_ignmis, .FALSE. )
+
+               CALL obs_pre_spm( spmdata(nspmsets), spmdatqc(nspmsets), &
+                  &                 ln_spm, ln_nea )
+            
+            ENDDO
+
+         ELSE                              ! Original file format
+
+            nspmsets = nspmsets + 1
+
+            CALL obs_rea_spm( 1, spmdata(nspmsets), jnumspm, &
+               &                 spmfiles(1:jnumspm), &
+               &                 nspmvars, nspmextr, nitend-nit000+2, &
+               &                 dobsini, dobsend, ln_ignmis, .FALSE. )
+
+            CALL obs_pre_spm( spmdata(nspmsets), spmdatqc(nspmsets), &
+               &                 ln_spm, ln_nea )
+
+         ENDIF
+ 
+      ENDIF
      
    END SUBROUTINE dia_obs_init
@@ -1113,4 +1207,5 @@
       !!               - Velocity component (U,V) profiles
       !!               - Sea surface log10(chlorophyll)
+      !!               - Sea surface spm
       !!
       !! ** Action  : 
@@ -1161,4 +1256,8 @@
       !USE ???                           ! ERSEM chlorophyll
 #endif
+#if defined key_spm
+      USE par_spm, ONLY: &              ! ERSEM/SPM sediments
+         & jp_spm
+#endif
       IMPLICIT NONE
 
@@ -1173,4 +1272,5 @@
       INTEGER :: jveloset               ! velocity profile data loop variable
       INTEGER :: jlogchlset             ! logchl data set loop variable
+      INTEGER :: jspmset                ! spm data set loop variable
       INTEGER :: jvar                   ! Variable number    
 #if ! defined key_lim2 && ! defined key_lim3
@@ -1182,4 +1282,7 @@
       REAL(wp), DIMENSION(jpi,jpj) :: &
          maskchl                        ! array for special chlorophyll mask
+      REAL(wp), DIMENSION(jpi,jpj) :: &
+         spm                            ! array for spm
+      INTEGER :: jn                     ! loop index
       CHARACTER(LEN=20) :: datestr=" ",timestr=" "
  
@@ -1324,4 +1427,22 @@
       ENDIF 
 
+      IF ( ln_spm ) THEN
+#if defined key_spm
+         spm(:,:) = 0.0
+         DO jn = 1, jp_spm
+            spm(:,:) = spm(:,:) + trn(:,:,1,jn)   ! sum SPM sizes
+         END DO
+#else
+         CALL ctl_stop( ' Trying to run spm observation operator', &
+            &           ' but no spm model appears to have been defined' )
+#endif
+
+         DO jspmset = 1, nspmsets
+             CALL obs_spm_opt( spmdatqc(jspmset),                &
+               &               kstp, jpi, jpj, nit000, spm(:,:), &
+               &               tmask(:,:,1), n2dint )
+         END DO         
+      ENDIF
+
 #if ! defined key_lim2 && ! defined key_lim3
       CALL wrk_dealloc(jpi,jpj,frld) 
@@ -1357,4 +1478,5 @@
       INTEGER :: jseaiceset               ! Sea Ice data set loop variable
       INTEGER :: jlogchlset               ! logchl data set loop variable
+      INTEGER :: jspmset                  ! spm data set loop variable
       INTEGER :: jset
       INTEGER :: jfbini
@@ -1628,4 +1750,25 @@
             WRITE(cdtmp,'(A,I2.2)')'logchlfb_',jlogchlset
             CALL obs_wri_logchl( cdtmp, logchldata(jlogchlset) )
+
+         END DO
+
+      ENDIF
+
+      !  - spm
+      IF ( ln_spm ) THEN
+
+         ! Copy data from spmdatqc to spmdata structures
+         DO jspmset = 1, nspmsets
+
+            CALL obs_surf_decompress( spmdatqc(jspmset), &
+                 &                    spmdata(jspmset), .TRUE., numout )
+
+         END DO
+
+         ! Write the spm data
+         DO jspmset = 1, nspmsets
+      
+            WRITE(cdtmp,'(A,I2.2)')'spmfb_',jspmset
+            CALL obs_wri_spm( cdtmp, spmdata(jspmset) )
 
          END DO
Index: branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
===================================================================
--- branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 6854)
+++ branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 6855)
@@ -24,4 +24,6 @@
    !!                    components of velocity from observations.
    !!   obs_logchl_opt : Compute the model counterpart of log10(chlorophyll)
+   !!                    observations
+   !!   obs_spm_opt :    Compute the model counterpart of spm
    !!                    observations
    !!----------------------------------------------------------------------
@@ -66,5 +68,6 @@
       &   obs_seaice_opt, &
       &   obs_vel_opt, &  ! Compute the model counterpart of velocity profile data
-      &   obs_logchl_opt  ! Compute the model counterpart of logchl data
+      &   obs_logchl_opt, & ! Compute the model counterpart of logchl data
+      &   obs_spm_opt     ! Compute the model counterpart of spm data
 
    INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types
@@ -2215,4 +2218,164 @@
    END SUBROUTINE obs_logchl_opt
 
+   SUBROUTINE obs_spm_opt( spmdatqc, kt, kpi, kpj, kit000, &
+      &                    pspmn, pspmmask, k2dint )
+
+      !!-----------------------------------------------------------------------
+      !!
+      !!                     ***  ROUTINE obs_spm_opt  ***
+      !!
+      !! ** Purpose : Compute the model counterpart of spm
+      !!              data by interpolating from the model grid to the 
+      !!              observation point.
+      !!
+      !! ** Method  : Linearly interpolate to each observation point using 
+      !!              the model values at the corners of the surrounding grid box.
+      !!
+      !!    The now model spm is first computed at the obs (lon, lat) point.
+      !!
+      !!    Several horizontal interpolation schemes are available:
+      !!        - distance-weighted (great circle) (k2dint = 0)
+      !!        - distance-weighted (small angle)  (k2dint = 1)
+      !!        - bilinear (geographical grid)     (k2dint = 2)
+      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
+      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
+      !!
+      !!
+      !! ** Action  :
+      !!
+      !! History :
+      !!      
+      !!-----------------------------------------------------------------------
+
+      !! * Modules used
+      USE obs_surf_def  ! Definition of storage space for surface observations
+
+      IMPLICIT NONE
+
+      !! * Arguments
+      TYPE(obs_surf), INTENT(INOUT) :: spmdatqc     ! Subset of surface data not failing screening
+      INTEGER, INTENT(IN) :: kt       ! Time step
+      INTEGER, INTENT(IN) :: kpi      ! Model grid parameters
+      INTEGER, INTENT(IN) :: kpj
+      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step 
+                                      !   (kit000-1 = restart time)
+      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header)
+      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
+         & pspmn,  &    ! Model spm field
+         & pspmmask     ! Land-sea mask
+         
+      !! * Local declarations
+      INTEGER :: ji
+      INTEGER :: jj
+      INTEGER :: jobs
+      INTEGER :: inrc
+      INTEGER :: ispm
+      INTEGER :: iobs
+       
+      REAL(KIND=wp) :: zlam
+      REAL(KIND=wp) :: zphi
+      REAL(KIND=wp) :: zext(1), zobsmask(1)
+      REAL(kind=wp), DIMENSION(2,2,1) :: &
+         & zweig
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
+         & zmask, &
+         & zspml, &
+         & zglam, &
+         & zgphi
+      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
+         & igrdi, &
+         & igrdj
+
+      !------------------------------------------------------------------------
+      ! Local initialization 
+      !------------------------------------------------------------------------
+      ! ... Record and data counters
+      inrc = kt - kit000 + 2
+      ispm = spmdatqc%nsstp(inrc)
+
+      ! Get the data for interpolation
+      
+      ALLOCATE( &
+         & igrdi(2,2,ispm), &
+         & igrdj(2,2,ispm), &
+         & zglam(2,2,ispm), &
+         & zgphi(2,2,ispm), &
+         & zmask(2,2,ispm), &
+         & zspml(2,2,ispm)  &
+         & )
+      
+      DO jobs = spmdatqc%nsurfup + 1, spmdatqc%nsurfup + ispm
+         iobs = jobs - spmdatqc%nsurfup
+         igrdi(1,1,iobs) = spmdatqc%mi(jobs)-1
+         igrdj(1,1,iobs) = spmdatqc%mj(jobs)-1
+         igrdi(1,2,iobs) = spmdatqc%mi(jobs)-1
+         igrdj(1,2,iobs) = spmdatqc%mj(jobs)
+         igrdi(2,1,iobs) = spmdatqc%mi(jobs)
+         igrdj(2,1,iobs) = spmdatqc%mj(jobs)-1
+         igrdi(2,2,iobs) = spmdatqc%mi(jobs)
+         igrdj(2,2,iobs) = spmdatqc%mj(jobs)
+      END DO
+      
+      CALL obs_int_comm_2d( 2, 2, ispm, &
+         &                  igrdi, igrdj, glamt, zglam )
+      CALL obs_int_comm_2d( 2, 2, ispm, &
+         &                  igrdi, igrdj, gphit, zgphi )
+      CALL obs_int_comm_2d( 2, 2, ispm, &
+         &                  igrdi, igrdj, pspmmask, zmask )
+      CALL obs_int_comm_2d( 2, 2, ispm, &
+         &                  igrdi, igrdj, pspmn, zspml )
+      
+      DO jobs = spmdatqc%nsurfup + 1, spmdatqc%nsurfup + ispm
+         
+         iobs = jobs - spmdatqc%nsurfup
+         
+         IF ( kt /= spmdatqc%mstp(jobs) ) THEN
+            
+            IF(lwp) THEN
+               WRITE(numout,*)
+               WRITE(numout,*) ' E R R O R : Observation',              &
+                  &            ' time step is not consistent with the', &
+                  &            ' model time step'
+               WRITE(numout,*) ' ========='
+               WRITE(numout,*)
+               WRITE(numout,*) ' Record  = ', jobs,                &
+                  &            ' kt      = ', kt,                  &
+                  &            ' mstp    = ', spmdatqc%mstp(jobs), &
+                  &            ' ntyp    = ', spmdatqc%ntyp(jobs)
+            ENDIF
+            CALL ctl_stop( 'obs_spm_opt', 'Inconsistent time' )
+            
+         ENDIF
+         
+         zlam = spmdatqc%rlam(jobs)
+         zphi = spmdatqc%rphi(jobs)
+         
+         ! Get weights to interpolate the model spm to the observation point
+         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
+            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
+            &                   zmask(:,:,iobs), zweig, zobsmask )
+         
+         ! ... Interpolate the model spm to the observation point
+         CALL obs_int_h2d( 1, 1,      &
+            &              zweig, zspml(:,:,iobs),  zext )
+         
+         spmdatqc%rmod(jobs,1) = zext(1)
+         
+      END DO
+      
+      ! Deallocate the data for interpolation
+      DEALLOCATE( &
+         & igrdi,    &
+         & igrdj,    &
+         & zglam,    &
+         & zgphi,    &
+         & zmask,    &
+         & zspml  &
+         & )
+      
+      spmdatqc%nsurfup = spmdatqc%nsurfup + ispm
+
+   END SUBROUTINE obs_spm_opt
+
 END MODULE obs_oper
 
Index: branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
===================================================================
--- branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 6854)
+++ branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 6855)
@@ -13,4 +13,5 @@
    !!   obs_pre_vel  : First level check and screening of velocity obs.
    !!   obs_pre_logchl : First level check and screening of logchl obs.
+   !!   obs_pre_spm  : First level check and screening of spm obs.
    !!   obs_scr      : Basic screening of the observations
    !!   obs_coo_tim  : Compute number of time steps to the observation time
@@ -43,4 +44,5 @@
       & obs_pre_vel, &     ! First level check and screening of velocity profiles
       & obs_pre_logchl, & ! First level check and screening of logchl data
+      & obs_pre_spm, &    ! First level check and screening of spm data
       & calc_month_len     ! Calculate the number of days in the months of a year  
 
@@ -1371,4 +1373,187 @@
    END SUBROUTINE obs_pre_logchl
 
+   SUBROUTINE obs_pre_spm( spmdata, spmdatqc, ld_spm, ld_nea )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE obs_pre_spm  ***
+      !!
+      !! ** Purpose : First level check and screening of spm observations
+      !!
+      !! ** Method  : First level check and screening of spm observations
+      !!
+      !! ** Action  : 
+      !!
+      !! References :
+      !!   
+      !! History :
+      !!----------------------------------------------------------------------
+      !! * Modules used
+      USE domstp              ! Domain: set the time-step
+      USE par_oce             ! Ocean parameters
+      USE dom_oce, ONLY : &   ! Geographical information
+         & glamt,   &
+         & gphit,   &
+         & tmask
+      !! * Arguments
+      TYPE(obs_surf), INTENT(INOUT) :: spmdata     ! Full set of spm data
+      TYPE(obs_surf), INTENT(INOUT) :: spmdatqc    ! Subset of spm data not failing screening
+      LOGICAL :: ld_spm     ! Switch for spm data
+      LOGICAL :: ld_nea        ! Switch for rejecting observation near land
+      !! * Local declarations
+      INTEGER :: iyea0         ! Initial date
+      INTEGER :: imon0         !  - (year, month, day, hour, minute)
+      INTEGER :: iday0    
+      INTEGER :: ihou0    
+      INTEGER :: imin0
+      INTEGER :: icycle       ! Current assimilation cycle
+                              ! Counters for observations that
+      INTEGER :: iotdobs      !  - outside time domain
+      INTEGER :: iosdsobs     !  - outside space domain
+      INTEGER :: ilansobs     !  - within a model land cell
+      INTEGER :: inlasobs     !  - close to land
+      INTEGER :: igrdobs      !  - fail the grid search
+                              ! Global counters for observations that
+      INTEGER :: iotdobsmpp   !  - outside time domain
+      INTEGER :: iosdsobsmpp  !  - outside space domain
+      INTEGER :: ilansobsmpp  !  - within a model land cell
+      INTEGER :: inlasobsmpp  !  - close to land
+      INTEGER :: igrdobsmpp   !  - fail the grid search
+      LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
+         & llvalid            ! data selection
+      INTEGER :: jobs         ! Obs. loop variable
+      INTEGER :: jstp         ! Time loop variable
+      INTEGER :: inrc         ! Time index variable
+
+      IF (lwp) WRITE(numout,*)'obs_pre_spm : Preparing the spm observations...'
+
+      ! Initial date initialization (year, month, day, hour, minute)
+      iyea0 =   ndate0 / 10000
+      imon0 = ( ndate0 - iyea0 * 10000 ) / 100
+      iday0 =   ndate0 - iyea0 * 10000 - imon0 * 100
+      ihou0 = 0
+      imin0 = 0
+
+      icycle = no     ! Assimilation cycle
+
+      ! Diagnostics counters for various failures.
+
+      iotdobs  = 0
+      igrdobs  = 0
+      iosdsobs = 0
+      ilansobs = 0
+      inlasobs = 0
+
+      ! -----------------------------------------------------------------------
+      ! Find time coordinate for spm data
+      ! -----------------------------------------------------------------------
+
+      CALL obs_coo_tim( icycle, &
+         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
+         &              spmdata%nsurf,   spmdata%nyea, spmdata%nmon, &
+         &              spmdata%nday,    spmdata%nhou, spmdata%nmin, &
+         &              spmdata%nqc,     spmdata%mstp, iotdobs        )
+      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
+      ! -----------------------------------------------------------------------
+      ! Check for spm data failing the grid search
+      ! -----------------------------------------------------------------------
+
+      CALL obs_coo_grd( spmdata%nsurf,   spmdata%mi, spmdata%mj, &
+         &              spmdata%nqc,     igrdobs                         )
+      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
+
+      ! -----------------------------------------------------------------------
+      ! Check for land points. 
+      ! -----------------------------------------------------------------------
+
+      CALL obs_coo_spc_2d( spmdata%nsurf,                 &
+         &                 jpi,             jpj,             &
+         &                 spmdata%mi,   spmdata%mj,   & 
+         &                 spmdata%rlam, spmdata%rphi, &
+         &                 glamt,           gphit,           &
+         &                 tmask(:,:,1),    spmdata%nqc,  &
+         &                 iosdsobs,        ilansobs,        &
+         &                 inlasobs,        ld_nea           ) 
+         
+      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
+      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
+      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
+
+      ! -----------------------------------------------------------------------
+      ! Copy useful data from the spmdata data structure to
+      ! the spmdatqc data structure 
+      ! -----------------------------------------------------------------------
+
+      ! Allocate the selection arrays
+
+      ALLOCATE( llvalid(spmdata%nsurf) )
+      
+      ! We want all data which has qc flags <= 0
+
+      llvalid(:)  = ( spmdata%nqc(:)  <= 10 )
+
+      ! The actual copying
+
+      CALL obs_surf_compress( spmdata,     spmdatqc,       .TRUE.,  numout, &
+         &                    lvalid=llvalid )
+
+      ! Dellocate the selection arrays
+      DEALLOCATE( llvalid )
+
+      ! -----------------------------------------------------------------------
+      ! Print information about what observations are left after qc
+      ! -----------------------------------------------------------------------
+
+      ! Update the total observation counter array
+      
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'obs_pre_spm :'
+         WRITE(numout,*) '~~~~~~~~~~~'
+         WRITE(numout,*)
+         WRITE(numout,*) ' spm data outside time domain                  = ', &
+            &            iotdobsmpp
+         WRITE(numout,*) ' Remaining spm data that failed grid search    = ', &
+            &            igrdobsmpp
+         WRITE(numout,*) ' Remaining spm data outside space domain       = ', &
+            &            iosdsobsmpp
+         WRITE(numout,*) ' Remaining spm data at land points             = ', &
+            &            ilansobsmpp
+         IF (ld_nea) THEN
+            WRITE(numout,*) ' Remaining spm data near land points (removed) = ', &
+               &            inlasobsmpp
+         ELSE
+            WRITE(numout,*) ' Remaining spm data near land points (kept)    = ', &
+               &            inlasobsmpp
+         ENDIF
+         WRITE(numout,*) ' spm data accepted                             = ', &
+            &            spmdatqc%nsurfmpp
+
+         WRITE(numout,*)
+         WRITE(numout,*) ' Number of observations per time step :'
+         WRITE(numout,*)
+         WRITE(numout,1997)
+         WRITE(numout,1998)
+      ENDIF
+      
+      DO jobs = 1, spmdatqc%nsurf
+         inrc = spmdatqc%mstp(jobs) + 2 - nit000
+         spmdatqc%nsstp(inrc)  = spmdatqc%nsstp(inrc) + 1
+      END DO
+      
+      CALL obs_mpp_sum_integers( spmdatqc%nsstp, spmdatqc%nsstpmpp, &
+         &                       nitend - nit000 + 2 )
+
+      IF ( lwp ) THEN
+         DO jstp = nit000 - 1, nitend
+            inrc = jstp - nit000 + 2
+            WRITE(numout,1999) jstp, spmdatqc%nsstpmpp(inrc)
+         END DO
+      ENDIF
+
+1997  FORMAT(10X,'Time step',5X,'spm data')
+1998  FORMAT(10X,'---------',5X,'------------')
+1999  FORMAT(10X,I9,5X,I17)
+      
+   END SUBROUTINE obs_pre_spm
+
    SUBROUTINE obs_coo_tim( kcycle, &
       &                    kyea0,   kmon0,   kday0,   khou0,   kmin0,     &
Index: branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_spm.F90
===================================================================
--- branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_spm.F90	(revision 6855)
+++ branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_spm.F90	(revision 6855)
@@ -0,0 +1,465 @@
+MODULE obs_read_spm
+   !!======================================================================
+   !!                       ***  MODULE obs_read_spm  ***
+   !! Observation diagnostics: Read the along track spm data from
+   !!                          GHRSST or any spm data from feedback files
+   !!======================================================================
+
+   !!----------------------------------------------------------------------
+   !!   obs_rea_spm : Driver for reading spm data from the feedback
+   !!----------------------------------------------------------------------
+
+   !! * Modules used   
+   USE par_kind                 ! Precision variables
+   USE in_out_manager           ! I/O manager
+   USE dom_oce                  ! Ocean space and time domain variables
+   USE obs_mpp                  ! MPP support routines for observation diagnostics
+   USE julian                   ! Julian date routines
+   USE obs_utils                ! Observation operator utility functions
+   USE obs_grid                 ! Grid search
+   USE obs_sort                 ! Sorting observation arrays
+   USE obs_surf_def             ! Surface observation definitions
+   USE obs_types                ! Observation type definitions
+   USE obs_spm_io            ! I/O for spm files
+   USE iom                      ! I/O
+   USE netcdf                   ! NetCDF library
+
+   IMPLICIT NONE
+
+   !! * Routine accessibility
+   PRIVATE
+
+   PUBLIC obs_rea_spm      ! Read the spm observations from the point data
+   
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE obs_rea_spm( kformat, &
+      &                    spmdata, knumfiles, cfilenames, &
+      &                    kvars, kextr, kstp, ddobsini, ddobsend, &
+      &                    ldignmis, ldmod )
+      !!---------------------------------------------------------------------
+      !!
+      !!                   *** ROUTINE obs_rea_spm ***
+      !!
+      !! ** Purpose : Read from file the spm data
+      !!
+      !! ** Method  : Depending on kformat either old or new style
+      !!              feedback data files are read
+      !!
+      !! ** Action  : 
+      !!
+      !!
+      !! History :  
+      !!      ! :  2009-01 (K. Mogensen) Initial version based on old versions
+      !!----------------------------------------------------------------------
+      !! * Modules used
+
+      !! * Arguments
+      INTEGER :: kformat   ! Format of input data 
+      !                    ! 0: New-style feedback
+      !                    ! 1: Old-style feedback
+      TYPE(obs_surf), INTENT(INOUT) :: &
+         & spmdata     ! spm data to be read
+      INTEGER, INTENT(IN) :: knumfiles   ! Number of corio format files to read in
+      CHARACTER(LEN=128), INTENT(IN) :: cfilenames(knumfiles) ! File names to read in
+      INTEGER, INTENT(IN) :: kvars    ! Number of variables in spmdata
+      INTEGER, INTENT(IN) :: kextr    ! Number of extra fields for each var in spmdata
+      INTEGER, INTENT(IN) :: kstp     ! Ocean time-step index
+      LOGICAL, INTENT(IN) :: ldignmis   ! Ignore missing files
+      LOGICAL, INTENT(IN) :: ldmod      ! Initialize model from input data
+      REAL(KIND=dp), INTENT(IN) :: ddobsini    ! Obs. ini time in YYYYMMDD.HHMMSS
+      REAL(KIND=dp), INTENT(IN) :: ddobsend    ! Obs. end time in YYYYMMDD.HHMMSS
+         
+      !! * Local declarations
+      CHARACTER(LEN=14), PARAMETER :: cpname='obs_rea_spm'
+      INTEGER :: ji
+      INTEGER :: jj
+      INTEGER :: jk
+      INTEGER :: iflag
+      INTEGER :: inobf
+      INTEGER :: i_file_id
+      INTEGER :: inowin
+      INTEGER :: iyea
+      INTEGER :: imon
+      INTEGER :: iday
+      INTEGER :: ihou
+      INTEGER :: imin
+      INTEGER :: isec
+      INTEGER, DIMENSION(knumfiles) :: &
+         & irefdate
+      INTEGER :: iobsmpp
+      INTEGER, PARAMETER :: ispmmaxtype = 1024
+      INTEGER, DIMENSION(0:ispmmaxtype) :: &
+         & ityp, &
+         & itypmpp
+      INTEGER, DIMENSION(:), ALLOCATABLE :: &
+         & iobsi,    &
+         & iobsj,    &
+         & iproc,    &
+         & iindx,    &
+         & ifileidx, &
+         & ispmidx
+      INTEGER :: itype
+      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
+         & zphi, &
+         & zlam
+      real(wp), DIMENSION(:), ALLOCATABLE :: &
+         & zdat
+      LOGICAL :: llvalprof
+      TYPE(obfbdata), POINTER, DIMENSION(:) :: &
+         & inpfiles
+      real(wp), DIMENSION(knumfiles) :: &
+         & djulini, &
+         & djulend
+      INTEGER :: iobs
+      INTEGER :: iobstot
+      INTEGER :: ios
+      INTEGER :: ioserrcount
+      CHARACTER(len=8) :: cl_refdate
+   
+      ! Local initialization
+      iobs = 0
+ 
+      !-----------------------------------------------------------------------
+      ! Check data the model part is just with feedback data files
+      !-----------------------------------------------------------------------
+      IF ( ldmod .AND. ( kformat /= 0 ) ) THEN
+         CALL ctl_stop( 'Model can only be read from feedback data' )
+         RETURN
+      ENDIF
+
+      !-----------------------------------------------------------------------
+      ! Count the number of files needed and allocate the obfbdata type
+      !-----------------------------------------------------------------------
+      
+      inobf = knumfiles
+      
+      ALLOCATE( inpfiles(inobf) )
+
+      spm_files : DO jj = 1, inobf
+          
+         !---------------------------------------------------------------------
+         ! Prints
+         !---------------------------------------------------------------------
+         IF(lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) ' obs_rea_spm : Reading from file = ', &
+               & TRIM( TRIM( cfilenames(jj) ) )
+            WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+            WRITE(numout,*)
+         ENDIF
+
+         !---------------------------------------------------------------------
+         !  Initialization: Open file and get dimensions only
+         !---------------------------------------------------------------------
+         
+         iflag = nf90_open( TRIM( TRIM( cfilenames(jj) ) ), nf90_nowrite, &
+            &                      i_file_id )
+         
+         IF ( iflag /= nf90_noerr ) THEN
+
+            IF ( ldignmis ) THEN
+               inpfiles(jj)%nobs = 0
+               CALL ctl_warn( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // &
+                  &           ' not found' )
+            ELSE 
+               CALL ctl_stop( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // &
+                  &           ' not found' )
+            ENDIF
+
+         ELSE 
+            
+            !------------------------------------------------------------------
+            !  Close the file since it is opened elsewhere
+            !------------------------------------------------------------------
+            
+            iflag = nf90_close( i_file_id )
+
+            !------------------------------------------------------------------
+            !  Read the file into inpfiles
+            !------------------------------------------------------------------
+            IF ( kformat == 0 ) THEN
+               CALL init_obfbdata( inpfiles(jj) )
+               IF(lwp) THEN
+                  WRITE(numout,*)
+                  WRITE(numout,*)'Reading from feedback file :', &
+                     &           TRIM( cfilenames(jj) )
+               ENDIF
+               CALL read_obfbdata( TRIM( cfilenames(jj) ), inpfiles(jj), &
+                  &                ldgrid = .TRUE. )
+               IF ( ldmod .AND. ( ( inpfiles(jj)%nadd == 0 ) .OR.&
+                  &               ( inpfiles(jj)%next < 2 ) ) ) THEN
+                  CALL ctl_stop( 'Model not in input data' )
+                  RETURN
+               ENDIF
+            ELSEIF ( kformat == 1) THEN
+               CALL read_spm( TRIM( cfilenames(jj) ), inpfiles(jj), &
+               &                 numout, lwp, .TRUE. )
+            ELSE
+               CALL ctl_stop( 'File format unknown' )
+            ENDIF
+
+            !------------------------------------------------------------------
+            !  Change longitude (-180,180)
+            !------------------------------------------------------------------
+
+            DO ji = 1, inpfiles(jj)%nobs 
+
+               IF ( inpfiles(jj)%plam(ji) < -180. ) &
+                  &   inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) + 360.
+
+               IF ( inpfiles(jj)%plam(ji) >  180. ) &
+                  &   inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) - 360.
+
+            END DO
+
+            !------------------------------------------------------------------
+            !  Calculate the date  (change eventually)
+            !------------------------------------------------------------------
+            cl_refdate=inpfiles(jj)%cdjuldref(1:8)
+            READ(cl_refdate,'(I8)') irefdate(jj)
+            
+            CALL ddatetoymdhms( ddobsini, iyea, imon, iday, ihou, imin, isec )
+            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulini(jj), &
+               &           krefdate = irefdate(jj) )
+            CALL ddatetoymdhms( ddobsend, iyea, imon, iday, ihou, imin, isec )
+            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulend(jj), &
+               &           krefdate = irefdate(jj) )
+            IF ( inpfiles(jj)%nobs > 0 ) THEN
+               inpfiles(jj)%iproc = -1
+               inpfiles(jj)%iobsi = -1
+               inpfiles(jj)%iobsj = -1
+            ENDIF
+            inowin = 0
+            DO ji = 1, inpfiles(jj)%nobs
+               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
+                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
+                  inowin = inowin + 1
+               ENDIF
+            END DO
+            ALLOCATE( zlam(inowin)  )
+            ALLOCATE( zphi(inowin)  )
+            ALLOCATE( iobsi(inowin) )
+            ALLOCATE( iobsj(inowin) )
+            ALLOCATE( iproc(inowin) )
+            inowin = 0
+            DO ji = 1, inpfiles(jj)%nobs
+               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
+                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
+                  inowin = inowin + 1
+                  zlam(inowin) = inpfiles(jj)%plam(ji)
+                  zphi(inowin) = inpfiles(jj)%pphi(ji)
+               ENDIF
+            END DO
+
+            CALL obs_grid_search( inowin, zlam, zphi, iobsi, iobsj, iproc, 'T' )
+
+            inowin = 0
+            DO ji = 1, inpfiles(jj)%nobs
+               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
+                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
+                  inowin = inowin + 1
+                  inpfiles(jj)%iproc(ji,1) = iproc(inowin)
+                  inpfiles(jj)%iobsi(ji,1) = iobsi(inowin)
+                  inpfiles(jj)%iobsj(ji,1) = iobsj(inowin)
+               ENDIF
+            END DO
+            DEALLOCATE( zlam, zphi, iobsi, iobsj, iproc )
+
+            DO ji = 1, inpfiles(jj)%nobs
+               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
+                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
+                  IF ( nproc == 0 ) THEN
+                     IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
+                  ELSE
+                     IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
+                  ENDIF
+                  llvalprof = .FALSE.
+                  IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. &
+                     & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN
+                     iobs = iobs + 1
+                  ENDIF
+               ENDIF
+            END DO
+
+         ENDIF
+
+      END DO spm_files
+
+      !-----------------------------------------------------------------------
+      ! Get the time ordered indices of the input data
+      !-----------------------------------------------------------------------
+
+      !---------------------------------------------------------------------
+      !  Loop over input data files to count total number of obs
+      !---------------------------------------------------------------------
+      iobstot = 0
+      DO jj = 1, inobf
+         DO ji = 1, inpfiles(jj)%nobs
+            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
+               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
+               iobstot = iobstot + 1
+            ENDIF
+         END DO
+      END DO
+
+      ALLOCATE( iindx(iobstot), ifileidx(iobstot), &
+         &      ispmidx(iobstot), zdat(iobstot) )
+      jk = 0
+      DO jj = 1, inobf
+         DO ji = 1, inpfiles(jj)%nobs
+            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
+               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
+               jk = jk + 1
+               ifileidx(jk) = jj
+               ispmidx(jk) = ji
+               zdat(jk)     = inpfiles(jj)%ptim(ji)
+            ENDIF
+         END DO
+      END DO
+      CALL sort_dp_indx( iobstot, &
+         &               zdat,     &
+         &               iindx   )
+      
+      CALL obs_surf_alloc( spmdata, iobs, & 
+                           kvars, kextr, kstp, jpi, jpj )
+      
+      ! * Read obs/positions, QC, all variable and assign to spmdata
+ 
+      iobs = 0
+
+      ityp   (:) = 0
+      itypmpp(:) = 0
+
+      ioserrcount=0      
+
+      DO jk = 1, iobstot
+         
+         jj = ifileidx(iindx(jk))
+         ji = ispmidx(iindx(jk))
+         IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  &
+            & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
+            
+            IF ( nproc == 0 ) THEN
+               IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
+            ELSE
+               IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
+            ENDIF
+            
+            ! Set observation information
+            
+            IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. &
+               & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN
+
+               iobs = iobs + 1
+
+               CALL jul2greg( isec,                   &
+                  &           imin,                   &
+                  &           ihou,                   &
+                  &           iday,                   &
+                  &           imon,                   &
+                  &           iyea,                   &
+                  &           inpfiles(jj)%ptim(ji), &
+                  &           irefdate(jj) )
+
+
+               ! spm time coordinates
+               spmdata%nyea(iobs) = iyea
+               spmdata%nmon(iobs) = imon
+               spmdata%nday(iobs) = iday
+               spmdata%nhou(iobs) = ihou
+               spmdata%nmin(iobs) = imin
+               
+               ! spm space coordinates
+               spmdata%rlam(iobs) = inpfiles(jj)%plam(ji)
+               spmdata%rphi(iobs) = inpfiles(jj)%pphi(ji)
+
+               ! Coordinate search parameters
+               spmdata%mi  (iobs) = inpfiles(jj)%iobsi(ji,1)
+               spmdata%mj  (iobs) = inpfiles(jj)%iobsj(ji,1)
+               
+               ! Instrument type
+               READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 901 ) itype
+901            IF ( ios /= 0 ) THEN
+                  IF (ioserrcount == 0) CALL ctl_warn ( 'Problem converting an instrument type to integer. Setting type to zero' ) 
+                  ioserrcount = ioserrcount + 1
+                  itype = 0
+               ENDIF
+               spmdata%ntyp(iobs) = itype
+               IF ( itype < ispmmaxtype + 1 ) THEN
+                  ityp(itype+1) = ityp(itype+1) + 1
+               ELSE
+                  IF(lwp)WRITE(numout,*)'WARNING:Increase ispmmaxtype in ',&
+                     &                  cpname
+               ENDIF
+
+               ! Bookkeeping data to match observations
+               spmdata%nsidx(iobs) = iobs
+               spmdata%nsfil(iobs) = iindx(jk)
+
+               ! QC flags
+               spmdata%nqc(iobs) = inpfiles(jj)%ivqc(ji,1)
+
+               ! Observed value
+               spmdata%robs(iobs,1) = inpfiles(jj)%pob(1,ji,1)
+
+
+               ! Model and MDT is set to fbrmdi unless read from file
+               IF ( ldmod ) THEN
+                  spmdata%rmod(iobs,1) = inpfiles(jj)%padd(1,ji,1,1)
+               ELSE
+                  spmdata%rmod(iobs,1) = fbrmdi
+               ENDIF
+            ENDIF
+         ENDIF
+
+      END DO
+
+      !-----------------------------------------------------------------------
+      ! Sum up over processors
+      !-----------------------------------------------------------------------
+      
+      CALL obs_mpp_sum_integer( iobs, iobsmpp )
+      
+      !-----------------------------------------------------------------------
+      ! Output number of observations.
+      !-----------------------------------------------------------------------
+      IF (lwp) THEN
+
+         WRITE(numout,*)
+         WRITE(numout,'(1X,A)')'spm data types'
+         WRITE(numout,'(1X,A)')'-----------------'
+         DO jj = 1,8
+            IF ( itypmpp(jj) > 0 ) THEN
+               WRITE(numout,'(1X,A4,I4,A3,I10)')'Type ', jj,' = ',itypmpp(jj)
+            ENDIF
+         END DO
+         WRITE(numout,'(1X,A50)')'--------------------------------------------------'
+         WRITE(numout,'(1X,A40,I10)')'Total                                 = ',iobsmpp
+         WRITE(numout,*)
+
+      ENDIF
+
+      !-----------------------------------------------------------------------
+      ! Deallocate temporary data
+      !-----------------------------------------------------------------------
+      DEALLOCATE( ifileidx, ispmidx, zdat )
+
+      !-----------------------------------------------------------------------
+      ! Deallocate input data
+      !-----------------------------------------------------------------------
+      DO jj = 1, inobf
+         CALL dealloc_obfbdata( inpfiles(jj) )
+      END DO
+      DEALLOCATE( inpfiles )
+
+   END SUBROUTINE obs_rea_spm
+
+END MODULE obs_read_spm
+
Index: branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_spm.F90
===================================================================
--- branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_spm.F90	(revision 6855)
+++ branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_spm.F90	(revision 6855)
@@ -0,0 +1,35 @@
+MODULE obs_spm
+   !!=====================================================================
+   !!                       ***  MODULE  obs_spm  ***
+   !! Observation diagnostics: Storage space for spm observations
+   !!                          arrays and additional flags etc.
+   !!=====================================================================
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+   
+   !! * Modules used 
+   USE obs_surf_def ! Definition of surface data types and tools
+
+   IMPLICIT NONE
+   
+   SAVE
+
+   !! * Routine accessibility
+   PRIVATE
+
+   PUBLIC nspmvars, nspmextr, nspmsets, spmdata, spmdatqc
+
+   !! * Shared Module variables
+   INTEGER :: nspmvars                               ! Number of spmdata variables
+   INTEGER :: nspmextr                               ! Number of spmdata extra 
+                                                     ! variables
+   INTEGER :: nspmsets                               ! Number of spmdata sets
+   TYPE(obs_surf), POINTER, DIMENSION(:) :: spmdata  ! Initial spm data
+   TYPE(obs_surf), POINTER, DIMENSION(:) :: spmdatqc ! Sea ice data after quality control
+
+END MODULE obs_spm
+
Index: branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_spm_io.F90
===================================================================
--- branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_spm_io.F90	(revision 6855)
+++ branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_spm_io.F90	(revision 6855)
@@ -0,0 +1,30 @@
+MODULE obs_spm_io
+   !!======================================================================
+   !!                       ***  MODULE obs_spm_io  ***
+   !! Observation operators : I/O for spm files
+   !!======================================================================
+   !! History : 
+   !!             !  09-01  (K. Mogensen) Initial version
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   read_spmfile    :  Read a obfbdata structure from a spm file
+   !!----------------------------------------------------------------------
+   USE par_kind
+   USE obs_utils
+   USE obs_fbm
+   USE julian
+   USE netcdf
+
+   IMPLICIT NONE
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+#include "obsspm_io.h90"
+
+END MODULE obs_spm_io
Index: branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90
===================================================================
--- branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90	(revision 6854)
+++ branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90	(revision 6855)
@@ -12,4 +12,5 @@
    !!   obs_wri_vel   : Write velocity observation diagnostics in NetCDF format
    !!   obs_wri_logchl: Write logchl observation related diagnostics
+   !!   obs_wri_spm   : Write spm observation related diagnostics
    !!   obs_wri_stats : Print basic statistics on the data being written out
    !!----------------------------------------------------------------------
@@ -47,4 +48,5 @@
       &   obs_wri_vel, &    ! Write velocity observation related diagnostics
       &   obs_wri_logchl, & ! Write logchl observation related diagnostics
+      &   obs_wri_spm, &    ! Write spm observation related diagnostics
       &   obswriinfo
    
@@ -1081,4 +1083,153 @@
    END SUBROUTINE obs_wri_logchl
 
+   SUBROUTINE obs_wri_spm( cprefix, spmdata, padd, pext )
+      !!-----------------------------------------------------------------------
+      !!
+      !!                     *** ROUTINE obs_wri_spm  ***
+      !!
+      !! ** Purpose : Write spm observation diagnostics
+      !!              related 
+      !!
+      !! ** Method  : NetCDF
+      !! 
+      !! ** Action  :
+      !!
+      !!-----------------------------------------------------------------------
+
+      !! * Modules used
+      IMPLICIT NONE
+
+      !! * Arguments
+      CHARACTER(LEN=*), INTENT(IN) :: cprefix       ! Prefix for output files
+      TYPE(obs_surf), INTENT(INOUT) :: spmdata   ! Full set of spm
+      TYPE(obswriinfo), OPTIONAL :: padd            ! Additional info for each variable
+      TYPE(obswriinfo), OPTIONAL :: pext            ! Extra info
+
+      !! * Local declarations 
+      TYPE(obfbdata) :: fbdata
+      CHARACTER(LEN=40) :: cfname             ! netCDF filename
+      CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_spm'
+      INTEGER :: jo
+      INTEGER :: ja
+      INTEGER :: je
+      INTEGER :: nadd
+      INTEGER :: next
+
+      IF ( PRESENT( padd ) ) THEN
+         nadd = padd%inum
+      ELSE
+         nadd = 0
+      ENDIF
+
+      IF ( PRESENT( pext ) ) THEN
+         next = pext%inum
+      ELSE
+         next = 0
+      ENDIF
+
+      CALL init_obfbdata( fbdata )
+
+      CALL alloc_obfbdata( fbdata, 1, spmdata%nsurf, 1, &
+         &                 1 + nadd, next, .TRUE. )
+
+      fbdata%cname(1)      = 'spm'
+      fbdata%coblong(1)    = 'spm'
+      fbdata%cobunit(1)    = 'g/m3'
+      DO je = 1, next
+         fbdata%cextname(je) = pext%cdname(je)
+         fbdata%cextlong(je) = pext%cdlong(je,1)
+         fbdata%cextunit(je) = pext%cdunit(je,1)
+      END DO
+      fbdata%caddname(1)   = 'Hx'
+      fbdata%caddlong(1,1) = 'Model interpolated spm'
+      fbdata%caddunit(1,1) = 'g/m3'
+      fbdata%cgrid(1)      = 'T'
+      DO ja = 1, nadd
+         fbdata%caddname(1+ja) = padd%cdname(ja)
+         fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
+         fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
+      END DO
+
+      WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*)'obs_wri_spm :'
+         WRITE(numout,*)'~~~~~~~~~~~~~~~~'
+         WRITE(numout,*)'Writing spm feedback file : ',TRIM(cfname)
+      ENDIF
+
+      ! Transform obs_prof data structure into obfbdata structure
+      fbdata%cdjuldref = '19500101000000'
+      DO jo = 1, spmdata%nsurf
+         fbdata%plam(jo)      = spmdata%rlam(jo)
+         fbdata%pphi(jo)      = spmdata%rphi(jo)
+         WRITE(fbdata%cdtyp(jo),'(I4)') spmdata%ntyp(jo)
+         fbdata%ivqc(jo,:)    = 0
+         fbdata%ivqcf(:,jo,:) = 0
+         IF ( spmdata%nqc(jo) > 10 ) THEN
+            fbdata%ioqc(jo)    = 4
+            fbdata%ioqcf(1,jo) = 0
+            fbdata%ioqcf(2,jo) = spmdata%nqc(jo) - 10
+         ELSE
+            fbdata%ioqc(jo)    = MAX(spmdata%nqc(jo),1)
+            fbdata%ioqcf(:,jo) = 0
+         ENDIF
+         fbdata%ipqc(jo)      = 0
+         fbdata%ipqcf(:,jo)   = 0
+         fbdata%itqc(jo)      = 0
+         fbdata%itqcf(:,jo)   = 0
+         fbdata%cdwmo(jo)     = ''
+         fbdata%kindex(jo)    = spmdata%nsfil(jo)
+         IF (ln_grid_global) THEN
+            fbdata%iobsi(jo,1) = spmdata%mi(jo)
+            fbdata%iobsj(jo,1) = spmdata%mj(jo)
+         ELSE
+            fbdata%iobsi(jo,1) = mig(spmdata%mi(jo))
+            fbdata%iobsj(jo,1) = mjg(spmdata%mj(jo))
+         ENDIF
+         CALL greg2jul( 0, &
+            &           spmdata%nmin(jo), &
+            &           spmdata%nhou(jo), &
+            &           spmdata%nday(jo), &
+            &           spmdata%nmon(jo), &
+            &           spmdata%nyea(jo), &
+            &           fbdata%ptim(jo),   &
+            &           krefdate = 19500101 )
+         fbdata%padd(1,jo,1,1) = spmdata%rmod(jo,1)
+         fbdata%pob(1,jo,1)    = spmdata%robs(jo,1)
+         fbdata%pdep(1,jo)     = 0.0
+         fbdata%idqc(1,jo)     = 0
+         fbdata%idqcf(:,1,jo)  = 0
+         IF ( spmdata%nqc(jo) > 10 ) THEN
+            fbdata%ivlqc(1,jo,1) = 4
+            fbdata%ivlqcf(1,1,jo,1) = 0
+            fbdata%ivlqcf(2,1,jo,1) = spmdata%nqc(jo) - 10
+         ELSE
+            fbdata%ivlqc(1,jo,1) = MAX(spmdata%nqc(jo),1)
+            fbdata%ivlqcf(:,1,jo,1) = 0
+         ENDIF
+         fbdata%iobsk(1,jo,1)  = 0
+         DO ja = 1, nadd
+            fbdata%padd(1,jo,1+ja,1) = &
+               & spmdata%rext(jo,padd%ipoint(ja))
+         END DO
+         DO je = 1, next
+            fbdata%pext(1,jo,je) = &
+               & spmdata%rext(jo,pext%ipoint(je))
+         END DO
+
+      END DO
+
+      ! Write the obfbdata structure
+      CALL write_obfbdata( cfname, fbdata )
+      
+      ! Output some basic statistics
+      CALL obs_wri_stats( fbdata )
+
+      CALL dealloc_obfbdata( fbdata )
+
+   END SUBROUTINE obs_wri_spm
+
    SUBROUTINE obs_wri_stats( fbdata )
       !!-----------------------------------------------------------------------
Index: branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obsspm_io.h90
===================================================================
--- branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obsspm_io.h90	(revision 6855)
+++ branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obsspm_io.h90	(revision 6855)
@@ -0,0 +1,270 @@
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+   SUBROUTINE read_spm( cdfilename, inpfile, kunit, ldwp, ldgrid )
+      !!---------------------------------------------------------------------
+      !!
+      !!                     ** ROUTINE read_spm **
+      !!
+      !! ** Purpose : Read from file the spm observations.
+      !!
+      !! ** Method  : The data file is a NetCDF file. 
+      !!
+      !! ** Action  :
+      !!
+      !! References : 
+      !!
+      !! 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=12),PARAMETER :: cpname = 'read_spm'
+      INTEGER :: i_file_id    ! netcdf IDS
+      INTEGER :: i_time_id
+      INTEGER :: i_ni_id
+      INTEGER :: i_data_id
+      INTEGER :: i_var_id
+      INTEGER :: i_data       ! Number of data per parameter in current file
+      INTEGER :: i_time       ! Number of reference times in file
+      INTEGER, DIMENSION(:), POINTER :: &
+         & i_reftime        ! Reference time in file in seconds since 1/1/1981.
+      INTEGER, DIMENSION(:,:), POINTER :: &
+         & i_dtime, &       ! Offset in seconds since reference time
+         & i_qc,    &       ! Quality control flag.
+         & i_type           ! Type of spm measurement.            
+      REAL(wp), DIMENSION(:), POINTER :: &
+         & z_phi,   &       ! Latitudes
+         & z_lam            ! Longitudes
+      REAL(wp), DIMENSION(:,:), POINTER :: &
+         & z_spm         ! spm data     
+      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(KIND=wp) :: zsca      ! Scale factor
+      REAL(KIND=wp) :: zoff      ! Offset for data in netcdf file
+      REAL(KIND=wp) :: z_offset  ! Offset for time conversion
+      REAL(KIND=wp) :: zfill     ! Fill value in netcdf file
+      CHARACTER (len=33) ::creftime     ! Reference time of file
+      INTEGER :: i_refyear       ! Integer version of reference time
+      INTEGER :: i_refmonth
+      INTEGER :: i_refday
+      INTEGER :: i_refhour
+      INTEGER :: i_refmin
+      INTEGER :: i_refsec
+      INTEGER :: ichunk
+      INTEGER :: jtim
+      INTEGER :: jobs
+      INTEGER :: iobs
+
+      CALL chkerr( nf90_open( TRIM( cdfilename ), nf90_nowrite, &
+         &                i_file_id, chunksize=ichunk), cpname, __LINE__ )
+      
+      ! Get the netCDF dimensions
+      
+      CALL chkerr( nf90_inq_dimid( i_file_id, 'time', i_time_id ),  &
+         &       cpname, __LINE__ )
+      CALL chkerr( nf90_inquire_dimension( i_file_id, i_time_id, &
+         &                              len = i_time ),  &
+         &      cpname, __LINE__ )
+      
+      CALL chkerr( nf90_inq_dimid( i_file_id, 'ni', i_ni_id ),  &
+         &       cpname, __LINE__ )
+      CALL chkerr( nf90_inquire_dimension( i_file_id, i_ni_id, &
+         &                              len = i_data ),  &
+         &      cpname, __LINE__ )
+      
+      
+      ! Allocate NetCDF variables
+      
+      ALLOCATE( &
+         & i_reftime    ( i_time                 ), &
+         & i_dtime      ( i_data,i_time  ), &      
+         & i_qc         ( i_data,i_time  ), & 
+         & i_type       ( i_data,i_time  ), &     
+         & z_phi        ( i_data                 ), &   
+         & z_lam        ( i_data                 ), &  
+         & z_spm        ( i_data,i_time  )  &
+         & )
+      
+      ! Get reference time of file which is in seconds since 1981/1/1 00:00. 
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, 'time', i_var_id ), & 
+         &       cpname, __LINE__ )
+      idims(1) = i_time
+      CALL chkdim( i_file_id, i_var_id, 1, idims, cpname, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, i_reftime),&
+         &       cpname, __LINE__ )
+      IF (nf90_inquire_attribute( i_file_id, i_var_id, "units") &
+         &                      == nf90_noerr) THEN
+         CALL chkerr( nf90_get_att( i_file_id, i_var_id, &
+            &                     "units",creftime), cpname,  __LINE__ )
+      ELSE
+         creftime = "seconds since 1981-01-01 00:00:00"
+      ENDIF
+      READ(creftime(15:18),*)i_refyear
+      READ(creftime(20:21),*)i_refmonth
+      READ(creftime(23:24),*)i_refday
+      READ(creftime(26:27),*)i_refhour
+      READ(creftime(29:30),*)i_refmin
+      READ(creftime(32:33),*)i_refsec
+      !Work out offset in days between reference time and 1/1/1950.
+      CALL greg2jul( i_refsec, i_refmin, i_refhour, i_refday, &
+         &           i_refmonth, i_refyear, z_offset)
+      
+      ! Get list of times for each ob in seconds relative to reference time
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, 'spm_dtime', i_var_id ), & 
+         &       cpname, __LINE__ )
+      idims(1) = i_data
+      idims(2) = i_time
+      CALL chkdim( i_file_id, i_var_id, 2, idims, cpname, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, i_dtime),&
+         &       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
+      zoff = 0.0
+      IF (nf90_inquire_attribute( i_file_id, i_var_id, "add_offset") &
+         &                      == nf90_noerr) THEN
+         CALL chkerr( nf90_get_att( i_file_id, i_var_id, &
+            &                     "add_offset",zoff), cpname,  __LINE__ )
+      ENDIF
+      i_dtime(:,:) = NINT((zsca * FLOAT(i_dtime(:,:))) &
+         &                 + zoff)
+      
+      ! Get longitudes
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, 'lon', i_var_id ), & 
+         &         cpname, __LINE__ )
+      idims(1) = i_data
+      CALL chkdim( i_file_id, i_var_id, 1, idims, cpname, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, z_lam),   &
+         &         cpname, __LINE__ )
+      
+      ! Get latitudes
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, 'lat', i_var_id ), & 
+         &         cpname, __LINE__ )
+      idims(1) = i_data
+      CALL chkdim( i_file_id, i_var_id, 1, idims, cpname, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, z_phi),   &
+         &         cpname, __LINE__ )
+      
+      ! Get spm data
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, 'spm', &
+         &                         i_var_id ), & 
+         &         cpname, __LINE__ )
+      idims(1) = i_data
+      idims(2) = i_time
+      CALL chkdim( i_file_id, i_var_id, 2, idims, cpname, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, z_spm), &
+         &       cpname, __LINE__ )
+      zoff = 0.
+      IF (nf90_inquire_attribute( i_file_id, i_var_id, "add_offset") &
+         &                      == nf90_noerr) THEN
+         CALL chkerr( nf90_get_att( i_file_id, i_var_id, &
+            &                     "add_offset", zoff), cpname, __LINE__ )
+      ENDIF
+      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(z_spm(:,:) /=  zfill)
+         z_spm(:,:) = (zsca * z_spm(:,:)) + zoff
+      ELSEWHERE
+         z_spm(:,:) = fbrmdi
+      END WHERE
+      
+      ! Get QC flag
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, 'confidence_flag', i_var_id ), & 
+         &       cpname, __LINE__ )
+      idims(1) = i_data
+      idims(2) = i_time
+      CALL chkdim( i_file_id, i_var_id, 2, idims, cpname, __LINE__ )
+      CALL chkerr( nf90_get_var  ( i_file_id, i_var_id, i_qc),   &
+            &       cpname, __LINE__ )
+      
+      ! Get spm obs type
+      
+      i_type(:,:)=1
+      
+      ! Close the file
+      
+      CALL chkerr( nf90_close( i_file_id ), cpname, __LINE__ )
+
+      ! Fill the obfbdata structure
+
+      ! Allocate obfbdata
+      
+      iobs = i_data * i_time
+      CALL init_obfbdata( inpfile )
+      CALL alloc_obfbdata( inpfile, 1, iobs, 1, 0, 0, ldgrid )
+      inpfile%cname(1) = 'spm'
+
+      ! Fill the obfbdata structure from input data
+
+      inpfile%cdjuldref = "19500101000000"
+      iobs = 0
+      DO jtim = 1, i_time
+         DO jobs = 1, i_data
+            iobs = iobs + 1
+            ! Characters
+            WRITE(inpfile%cdwmo(iobs),'(A6,A2)') 'spm','  '
+            WRITE(inpfile%cdtyp(iobs),'(I4)') i_type(jobs,jtim)
+            ! Real values
+            inpfile%plam(iobs)         = z_lam(jobs)
+            inpfile%pphi(iobs)         = z_phi(jobs)
+            inpfile%pob(1,iobs,1)      = z_spm(jobs,jtim)
+            inpfile%ptim(iobs)         = &
+               & REAL(i_reftime(jtim))/(60.*60.*24.) + &
+               & z_offset + REAL(i_dtime(jobs,jtim))/(60.*60.*24.)
+            inpfile%pdep(1,iobs)       = 0.0
+            ! Integers
+            inpfile%kindex(iobs)       = iobs
+            IF ( z_spm(jobs,jtim) == fbrmdi ) THEN
+               inpfile%ioqc(iobs)      = 4
+               inpfile%ivqc(iobs,1)    = 4 
+               inpfile%ivlqc(1,iobs,1) = 4
+            ELSE 
+               inpfile%ioqc(iobs)      = i_qc(jobs,jtim)
+               inpfile%ivqc(iobs,1)    = i_qc(jobs,jtim)
+               inpfile%ivlqc(1,iobs,1) = 1
+            ENDIF
+            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
+
+   END SUBROUTINE read_spm
+
+
