Index: /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
===================================================================
--- /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 7712)
+++ /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 7713)
@@ -21,4 +21,5 @@
    USE par_oce
    USE dom_oce                  ! Ocean space and time domain variables
+   USE obs_const, ONLY: obfillflt ! Fill value
    USE obs_fbm, ONLY: ln_cl4    ! Class 4 diagnostic switch
    USE obs_read_prof            ! Reading and allocation of observations (Coriolis)
@@ -29,4 +30,8 @@
    USE obs_read_seaice          ! Reading and allocation of Sea Ice observations  
    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_read_fco2            ! Reading and allocation of fco2 observations
+   USE obs_read_pco2            ! Reading and allocation of pco2 observations
    USE obs_prep                 ! Preparation of obs. (grid search etc).
    USE obs_oper                 ! Observation operators
@@ -40,4 +45,8 @@
    USE obs_sst                  ! SST data storage
    USE obs_seaice               ! Sea Ice data storage
+   USE obs_logchl               ! logchl data storage
+   USE obs_spm                  ! spm data storage
+   USE obs_fco2                 ! fco2 data storage
+   USE obs_pco2                 ! pco2 data storage
    USE obs_types                ! Definitions for observation types
    USE mpp_map                  ! MPP mapping
@@ -81,4 +90,12 @@
    LOGICAL, PUBLIC :: ln_velhradcp   !: Logical switch for raw high freq netCDF ADCP vel. data 
    LOGICAL, PUBLIC :: ln_velfb       !: Logical switch for velocities from feedback files
+   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_fco2        !: Logical switch for fco2
+   LOGICAL, PUBLIC :: ln_fco2fb      !: Logical switch for fco2 from feedback files
+   LOGICAL, PUBLIC :: ln_pco2        !: Logical switch for pco2
+   LOGICAL, PUBLIC :: ln_pco2fb      !: Logical switch for pco2 from feedback files
    LOGICAL, PUBLIC :: ln_ssh         !: Logical switch for sea surface height
    LOGICAL, PUBLIC :: ln_sss         !: Logical switch for sea surface salinity
@@ -164,4 +181,12 @@
       CHARACTER(len=128) :: velhradcpfiles(MaxNumFiles)
       CHARACTER(len=128) :: velfbfiles(MaxNumFiles)
+      CHARACTER(len=128) :: logchlfiles(MaxNumFiles)
+      CHARACTER(len=128) :: logchlfbfiles(MaxNumFiles)
+      CHARACTER(len=128) :: spmfiles(MaxNumFiles)
+      CHARACTER(len=128) :: spmfbfiles(MaxNumFiles)
+      CHARACTER(len=128) :: fco2files(MaxNumFiles)
+      CHARACTER(len=128) :: fco2fbfiles(MaxNumFiles)
+      CHARACTER(len=128) :: pco2files(MaxNumFiles)
+      CHARACTER(len=128) :: pco2fbfiles(MaxNumFiles)
       CHARACTER(LEN=128) :: reysstname
       CHARACTER(LEN=12)  :: reysstfmt
@@ -189,4 +214,12 @@
          &            ln_velhradcp, velhradcpfiles,                   &
          &            ln_velfb, velfbfiles, ln_velfb_av,              &
+         &            ln_logchl, ln_logchlfb,                         &
+         &            logchlfiles, logchlfbfiles,                     &
+         &            ln_spm, ln_spmfb,                               &
+         &            spmfiles, spmfbfiles,                           &
+         &            ln_fco2, ln_fco2fb,                             &
+         &            fco2files, fco2fbfiles,                         &
+         &            ln_pco2, ln_pco2fb,                             &
+         &            pco2files, pco2fbfiles,                         &
          &            ln_profb_enatim, ln_ignmis, ln_cl4,             &
          &            ln_sstbias, sstbias_files
@@ -210,4 +243,12 @@
       INTEGER :: jnumvelhradcp   
       INTEGER :: jnumvelfb
+      INTEGER :: jnumlogchl
+      INTEGER :: jnumlogchlfb
+      INTEGER :: jnumspm
+      INTEGER :: jnumspmfb
+      INTEGER :: jnumfco2
+      INTEGER :: jnumfco2fb
+      INTEGER :: jnumpco2
+      INTEGER :: jnumpco2fb
       INTEGER :: ji
       INTEGER :: jset
@@ -218,4 +259,13 @@
       ! Read namelist parameters
       !-----------------------------------------------------------------------
+
+      ln_logchl   = .FALSE.
+      ln_logchlfb = .FALSE.
+      ln_spm      = .FALSE.
+      ln_spmfb    = .FALSE.
+      ln_fco2     = .FALSE.
+      ln_fco2fb   = .FALSE.
+      ln_pco2     = .FALSE.
+      ln_pco2fb   = .FALSE.
       
       !Initalise all values in namelist arrays
@@ -238,4 +288,12 @@
       velcurfiles(:) = ''
       veladcpfiles(:) = ''
+      logchlfiles(:) = ''
+      logchlfbfiles(:) = ''
+      spmfiles(:) = ''
+      spmfbfiles(:) = ''
+      fco2files(:) = ''
+      fco2fbfiles(:) = ''
+      pco2files(:) = ''
+      pco2fbfiles(:) = ''
       sstbias_files(:) = ''
       endailyavtypes(:) = -1
@@ -337,4 +395,44 @@
          jnumvelfb = COUNT(lmask)
          lmask(:) = .FALSE.
+      ENDIF
+      IF (ln_logchl) THEN
+         lmask(:) = .FALSE.
+         WHERE (logchlfiles(:) /= '') lmask(:) = .TRUE.
+         jnumlogchl = COUNT(lmask)
+      ENDIF
+      IF (ln_logchlfb) THEN
+         lmask(:) = .FALSE.
+         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
+      IF (ln_fco2) THEN
+         lmask(:) = .FALSE.
+         WHERE (fco2files(:) /= '') lmask(:) = .TRUE.
+         jnumfco2 = COUNT(lmask)
+      ENDIF
+      IF (ln_fco2fb) THEN
+         lmask(:) = .FALSE.
+         WHERE (fco2fbfiles(:) /= '') lmask(:) = .TRUE.
+         jnumfco2fb = COUNT(lmask)
+      ENDIF
+      IF (ln_pco2) THEN
+         lmask(:) = .FALSE.
+         WHERE (pco2files(:) /= '') lmask(:) = .TRUE.
+         jnumpco2 = COUNT(lmask)
+      ENDIF
+      IF (ln_pco2fb) THEN
+         lmask(:) = .FALSE.
+         WHERE (pco2fbfiles(:) /= '') lmask(:) = .TRUE.
+         jnumpco2fb = COUNT(lmask)
       ENDIF
       
@@ -368,4 +466,12 @@
          WRITE(numout,*) '             Logical switch for velocity high freq. ADCP  ln_velhradcp = ', ln_velhradcp
          WRITE(numout,*) '             Logical switch for feedback velocity data        ln_velfb = ', ln_velfb
+         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,*) '             Logical switch for fco2 observations              ln_fco2 = ', ln_fco2
+         WRITE(numout,*) '             Logical switch for pco2 observations              ln_pco2 = ', ln_pco2
+         WRITE(numout,*) '             Logical switch for feedback pco2 data           ln_pco2fb = ', ln_pco2fb
+         WRITE(numout,*) '             Logical switch for feedback fco2 data           ln_fco2fb = ', ln_fco2fb
          WRITE(numout,*) '             Global distribtion of observations         ln_grid_global = ',ln_grid_global
          WRITE(numout,*) &
@@ -464,4 +570,52 @@
                      TRIM(velfbfiles(ji))
                ENDIF
+            END DO
+         ENDIF
+         IF (ln_logchl) THEN
+            DO ji = 1, jnumlogchl
+               WRITE(numout,'(1X,2A)') '             logchl input observation file name        logchlfiles = ', &
+                  TRIM(logchlfiles(ji))
+            END DO
+         ENDIF
+         IF (ln_logchlfb) THEN
+            DO ji = 1, jnumlogchlfb
+               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
+         IF (ln_fco2) THEN
+            DO ji = 1, jnumfco2
+               WRITE(numout,'(1X,2A)') '             fco2 input observation file name            fco2files = ', &
+                  TRIM(fco2files(ji))
+            END DO
+         ENDIF
+         IF (ln_fco2fb) THEN
+            DO ji = 1, jnumfco2fb
+               WRITE(numout,'(1X,2A)') '            Feedback fco2 input observation file name  fco2fbfiles = ', &
+                  TRIM(fco2fbfiles(ji))
+            END DO
+         ENDIF
+         IF (ln_pco2) THEN
+            DO ji = 1, jnumpco2
+               WRITE(numout,'(1X,2A)') '             pco2 input observation file name            pco2files = ', &
+                  TRIM(pco2files(ji))
+            END DO
+         ENDIF
+         IF (ln_pco2fb) THEN
+            DO ji = 1, jnumpco2fb
+               WRITE(numout,'(1X,2A)') '            Feedback pco2 input observation file name  pco2fbfiles = ', &
+                  TRIM(pco2fbfiles(ji))
             END DO
          ENDIF
@@ -501,8 +655,10 @@
          & ( .NOT. ln_vel3d ).AND.                                         &
          & ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. &
-         & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ) ) THEN
+         & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ).AND.( .NOT. ln_logchl ).AND. &
+         & ( .NOT. ln_spm ).AND.( .NOT. ln_fco2 ).AND.( .NOT. ln_pco2 ) ) 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 are all set to .FALSE.'
+            &                    ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d,', &
+            &                    ' ln_logchl, ln_spm, ln_fco2, ln_pco2 are all set to .FALSE.'
          nwarn = nwarn + 1
       ENDIF
@@ -1002,4 +1158,224 @@
 
       ENDIF
+
+      !  - log10(chlorophyll)
+      
+      IF ( ln_logchl ) THEN
+
+         ! Set the number of variables for logchl to 1
+         nlogchlvars = 1
+
+         ! Set the number of extra variables for logchl to 0
+         nlogchlextr = 0
+         
+         IF ( ln_logchlfb ) THEN
+            nlogchlsets = jnumlogchlfb
+         ELSE
+            nlogchlsets = 1
+         ENDIF
+
+         ALLOCATE(logchldata(nlogchlsets))
+         ALLOCATE(logchldatqc(nlogchlsets))
+         logchldata(:)%nsurf=0
+         logchldatqc(:)%nsurf=0
+
+         nlogchlsets = 0
+
+         IF ( ln_logchlfb ) THEN             ! Feedback file format
+
+            DO jset = 1, jnumlogchlfb
+            
+               nlogchlsets = nlogchlsets + 1
+
+               CALL obs_rea_logchl( 0, logchldata(nlogchlsets), 1, &
+                  &                 logchlfbfiles(jset:jset), &
+                  &                 nlogchlvars, nlogchlextr, nitend-nit000+2, &
+                  &                 dobsini, dobsend, ln_ignmis, .FALSE. )
+
+               CALL obs_pre_logchl( logchldata(nlogchlsets), logchldatqc(nlogchlsets), &
+                  &                 ln_logchl, ln_nea )
+            
+            ENDDO
+
+         ELSE                              ! Original file format
+
+            nlogchlsets = nlogchlsets + 1
+
+            CALL obs_rea_logchl( 1, logchldata(nlogchlsets), jnumlogchl, &
+               &                 logchlfiles(1:jnumlogchl), &
+               &                 nlogchlvars, nlogchlextr, nitend-nit000+2, &
+               &                 dobsini, dobsend, ln_ignmis, .FALSE. )
+
+            CALL obs_pre_logchl( logchldata(nlogchlsets), logchldatqc(nlogchlsets), &
+               &                 ln_logchl, ln_nea )
+
+         ENDIF
+ 
+      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
+
+      !  - fco2
+      
+      IF ( ln_fco2 ) THEN
+
+         ! Set the number of variables for fco2 to 1
+         nfco2vars = 1
+
+         ! Set the number of extra variables for fco2 to 0
+         nfco2extr = 0
+         
+         IF ( ln_fco2fb ) THEN
+            nfco2sets = jnumfco2fb
+         ELSE
+            nfco2sets = 1
+         ENDIF
+
+         ALLOCATE(fco2data(nfco2sets))
+         ALLOCATE(fco2datqc(nfco2sets))
+         fco2data(:)%nsurf=0
+         fco2datqc(:)%nsurf=0
+
+         nfco2sets = 0
+
+         IF ( ln_fco2fb ) THEN             ! Feedback file format
+
+            DO jset = 1, jnumfco2fb
+            
+               nfco2sets = nfco2sets + 1
+
+               CALL obs_rea_fco2( 0, fco2data(nfco2sets), 1, &
+                  &                 fco2fbfiles(jset:jset), &
+                  &                 nfco2vars, nfco2extr, nitend-nit000+2, &
+                  &                 dobsini, dobsend, ln_ignmis, .FALSE. )
+
+               CALL obs_pre_fco2( fco2data(nfco2sets), fco2datqc(nfco2sets), &
+                  &                 ln_fco2, ln_nea )
+            
+            ENDDO
+
+         ELSE                              ! Original file format
+
+            nfco2sets = nfco2sets + 1
+
+            CALL obs_rea_fco2( 1, fco2data(nfco2sets), jnumfco2, &
+               &                 fco2files(1:jnumfco2), &
+               &                 nfco2vars, nfco2extr, nitend-nit000+2, &
+               &                 dobsini, dobsend, ln_ignmis, .FALSE. )
+
+            CALL obs_pre_fco2( fco2data(nfco2sets), fco2datqc(nfco2sets), &
+               &                 ln_fco2, ln_nea )
+
+         ENDIF
+ 
+      ENDIF
+
+      !  - pco2
+      
+      IF ( ln_pco2 ) THEN
+
+         ! Set the number of variables for pco2 to 1
+         npco2vars = 1
+
+         ! Set the number of extra variables for pco2 to 0
+         npco2extr = 0
+         
+         IF ( ln_pco2fb ) THEN
+            npco2sets = jnumpco2fb
+         ELSE
+            npco2sets = 1
+         ENDIF
+
+         ALLOCATE(pco2data(npco2sets))
+         ALLOCATE(pco2datqc(npco2sets))
+         pco2data(:)%nsurf=0
+         pco2datqc(:)%nsurf=0
+
+         npco2sets = 0
+
+         IF ( ln_pco2fb ) THEN             ! Feedback file format
+
+            DO jset = 1, jnumpco2fb
+            
+               npco2sets = npco2sets + 1
+
+               CALL obs_rea_pco2( 0, pco2data(npco2sets), 1, &
+                  &                 pco2fbfiles(jset:jset), &
+                  &                 npco2vars, npco2extr, nitend-nit000+2, &
+                  &                 dobsini, dobsend, ln_ignmis, .FALSE. )
+
+               CALL obs_pre_pco2( pco2data(npco2sets), pco2datqc(npco2sets), &
+                  &                 ln_pco2, ln_nea )
+            
+            ENDDO
+
+         ELSE                              ! Original file format
+
+            npco2sets = npco2sets + 1
+
+            CALL obs_rea_pco2( 1, pco2data(npco2sets), jnumpco2, &
+               &                 pco2files(1:jnumpco2), &
+               &                 npco2vars, npco2extr, nitend-nit000+2, &
+               &                 dobsini, dobsend, ln_ignmis, .FALSE. )
+
+            CALL obs_pre_pco2( pco2data(npco2sets), pco2datqc(npco2sets), &
+               &                 ln_pco2, ln_nea )
+
+         ENDIF
+ 
+      ENDIF
      
    END SUBROUTINE dia_obs_init
@@ -1019,4 +1395,8 @@
       !!               - Sea surface salinity
       !!               - Velocity component (U,V) profiles
+      !!               - Sea surface log10(chlorophyll)
+      !!               - Sea surface spm
+      !!               - Sea surface fco2
+      !!               - Sea surface pco2
       !!
       !! ** Action  : 
@@ -1043,5 +1423,6 @@
          & tmask, umask, vmask                            
       USE phycst, ONLY : &              ! Physical constants
-         & rday                         
+         & rday, &
+         & rt0
       USE oce, ONLY : &                 ! Ocean dynamics and tracers variables
          & tsn,  &             
@@ -1056,4 +1437,26 @@
          & frld
 #endif
+#if defined key_hadocc
+      USE trc, ONLY :  &                ! HadOCC chlorophyll, fCO2 and pCO2
+         & HADOCC_CHL, &
+         & HADOCC_FCO2, &
+         & HADOCC_PCO2, &
+         & HADOCC_FILL_FLT
+#elif defined key_medusa && defined key_foam_medusa
+      USE trc, ONLY :  &                ! MEDUSA chlorophyll, fCO2 and pCO2
+         & MEDUSA_CHL, &
+         & MEDUSA_FCO2, &
+         & MEDUSA_PCO2, &
+         & MEDUSA_FILL_FLT
+#elif defined key_fabm
+      USE fabm
+      USE par_fabm
+#endif
+#if defined key_spm
+      USE par_spm, ONLY: &              ! ERSEM/SPM sediments
+         & jp_spm
+      USE trc, ONLY :  &
+         & trn
+#endif
       IMPLICIT NONE
 
@@ -1067,7 +1470,31 @@
       INTEGER :: jseaiceset             ! sea ice data set loop variable
       INTEGER :: jveloset               ! velocity profile data loop variable
+      INTEGER :: jlogchlset             ! logchl data set loop variable
+      INTEGER :: jspmset                ! spm data set loop variable
+      INTEGER :: jfco2set               ! fco2 data set loop variable
+      INTEGER :: jpco2set               ! pco2 data set loop variable
       INTEGER :: jvar                   ! Variable number    
 #if ! defined key_lim2 && ! defined key_lim3
       REAL(wp), POINTER, DIMENSION(:,:) :: frld   
+#endif
+      REAL(wp) :: tiny                  ! small number
+      REAL(wp), DIMENSION(jpi,jpj) :: &
+         logchl                         ! array for log chlorophyll
+      REAL(wp), DIMENSION(jpi,jpj) :: &
+         maskchl                        ! array for special chlorophyll mask
+      REAL(wp), DIMENSION(jpi,jpj) :: &
+         spm                            ! array for spm
+      REAL(wp), DIMENSION(jpi,jpj) :: &
+         fco2                           ! array for fco2
+      REAL(wp), DIMENSION(jpi,jpj) :: &
+         maskfco2                       ! array for special fco2 mask
+      REAL(wp), DIMENSION(jpi,jpj) :: &
+         pco2                           ! array for pco2
+      REAL(wp), DIMENSION(jpi,jpj) :: &
+         maskpco2                       ! array for special pco2 mask
+      INTEGER :: jn                     ! loop index
+#if defined key_fabm
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: logchl_3d
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: pco2_3d
 #endif
       CHARACTER(LEN=20) :: datestr=" ",timestr=" "
@@ -1180,4 +1607,144 @@
       ENDIF
 
+      IF ( ln_logchl ) THEN
+
+#if defined key_hadocc
+         logchl(:,:) = HADOCC_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC
+#elif defined key_medusa && defined key_foam_medusa
+         logchl(:,:) = MEDUSA_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC
+#elif defined key_fabm
+         logchl_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabmdia_chltot)
+         logchl(:,:) = logchl_3d(:,:,1)
+#else
+         CALL ctl_stop( ' Trying to run logchl observation operator', &
+            &           ' but no biogeochemical model appears to have been defined' )
+#endif
+
+         maskchl(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things
+
+         ! Take the log10 where we can, otherwise exclude
+         tiny = 1.0e-20
+         WHERE(logchl(:,:) > tiny .AND. logchl(:,:) /= obfillflt )
+            logchl(:,:)  = LOG10(logchl(:,:))
+         ELSEWHERE
+            logchl(:,:)  = obfillflt
+            maskchl(:,:) = 0
+         END WHERE
+
+         DO jlogchlset = 1, nlogchlsets
+             CALL obs_logchl_opt( logchldatqc(jlogchlset),             &
+               &                  kstp, jpi, jpj, nit000, logchl(:,:), &
+               &                  maskchl(:,:), n2dint )
+         END DO         
+      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 ( ln_fco2 ) THEN
+         maskfco2(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things
+#if defined key_hadocc
+         fco2(:,:) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC
+         IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ).AND.( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN
+            fco2(:,:) = obfillflt
+            maskfco2(:,:) = 0
+            CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', &
+               &           ' on timestep ' // TRIM(STR(kstp)),                              &
+               &           ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' )
+         ENDIF
+#elif defined key_medusa && defined key_foam_medusa
+         fco2(:,:) = MEDUSA_FCO2(:,:)    ! fCO2 from MEDUSA
+         IF ( ( MINVAL( MEDUSA_FCO2 ) == MEDUSA_FILL_FLT ).AND.( MAXVAL( MEDUSA_FCO2 ) == MEDUSA_FILL_FLT ) ) THEN
+            fco2(:,:) = obfillflt
+            maskfco2(:,:) = 0
+            CALL ctl_warn( ' MEDUSA fCO2 values masked out for observation operator', &
+               &           ' on timestep ' // TRIM(STR(kstp)),                              &
+               &           ' as MEDUSA_FCO2(:,:) == MEDUSA_FILL_FLT' )
+         ENDIF
+#elif defined key_fabm
+         ! First, get pCO2 from FABM
+         pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc)
+         pco2(:,:) = pco2_3d(:,:,1)
+         ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of:
+         ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems
+         ! and data reduction routines, Deep-Sea Research II, 56: 512-522.
+         ! and
+         ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas,
+         ! Marine Chemistry, 2: 203-215.
+         ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so
+         ! not explicitly included - atmospheric pressure is not necessarily available so this is
+         ! the best assumption.
+         ! Further, the (1-xCO2)^2 term has been neglected. This is common practice
+         ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes)
+         ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal
+         ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway.
+         fco2(:,:) = pco2(:,:) * EXP((-1636.75                                                                               + &
+            &                         12.0408      * (tsn(:,:,1,jp_tem)+rt0)                                                 - &
+            &                         0.0327957    * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)                         + &
+            &                         0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + &
+            &                         2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0)))                                        / &
+            &                        (82.0578 * (tsn(:,:,1,jp_tem)+rt0)))
+#else
+         CALL ctl_stop( ' Trying to run fco2 observation operator', &
+            &           ' but no biogeochemical model appears to have been defined' )
+#endif
+
+         DO jfco2set = 1, nfco2sets
+             CALL obs_fco2_opt( fco2datqc(jfco2set),                      &
+               &                kstp, jpi, jpj, nit000, fco2(:,:), &
+               &                maskfco2(:,:), n2dint )
+         END DO
+      ENDIF
+
+      IF ( ln_pco2 ) THEN
+         maskpco2(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things
+#if defined key_hadocc
+         pco2(:,:) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC
+         IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ).AND.( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN
+            pco2(:,:) = obfillflt
+            maskpco2(:,:) = 0
+            CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', &
+               &           ' on timestep ' // TRIM(STR(kstp)),                              &
+               &           ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' )
+         ENDIF
+#elif defined key_medusa && defined key_foam_medusa
+         pco2(:,:) = MEDUSA_PCO2(:,:)    ! pCO2 from MEDUSA
+         IF ( ( MINVAL( MEDUSA_PCO2 ) == MEDUSA_FILL_FLT ).AND.( MAXVAL( MEDUSA_PCO2 ) == MEDUSA_FILL_FLT ) ) THEN
+            pco2(:,:) = obfillflt
+            maskpco2(:,:) = 0
+            CALL ctl_warn( ' MEDUSA pCO2 values masked out for observation operator', &
+               &           ' on timestep ' // TRIM(STR(kstp)),                              &
+               &           ' as MEDUSA_PCO2(:,:) == MEDUSA_FILL_FLT' )
+         ENDIF
+#elif defined key_fabm
+         pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc)
+         pco2(:,:) = pco2_3d(:,:,1)
+#else
+         CALL ctl_stop( ' Trying to run pCO2 observation operator', &
+            &           ' but no biogeochemical model appears to have been defined' )
+#endif
+
+         DO jpco2set = 1, npco2sets
+             CALL obs_pco2_opt( pco2datqc(jpco2set),                      &
+               &                kstp, jpi, jpj, nit000, pco2(:,:), &
+               &                maskpco2(:,:), n2dint )
+         END DO
+      ENDIF
+
 #if ! defined key_lim2 && ! defined key_lim3
       CALL wrk_dealloc(jpi,jpj,frld) 
@@ -1212,8 +1779,12 @@
       INTEGER :: jsstset                  ! SST data set loop variable
       INTEGER :: jseaiceset               ! Sea Ice data set loop variable
+      INTEGER :: jlogchlset               ! logchl data set loop variable
+      INTEGER :: jspmset                  ! spm data set loop variable
+      INTEGER :: jfco2set                 ! fco2 data set loop variable
+      INTEGER :: jpco2set                 ! pco2 data set loop variable
       INTEGER :: jset
       INTEGER :: jfbini
       CHARACTER(LEN=20) :: datestr=" ",timestr=" "
-      CHARACTER(LEN=10) :: cdtmp
+      CHARACTER(LEN=20) :: cdtmp
       !-----------------------------------------------------------------------
       ! Depending on switches call various observation output routines
@@ -1456,4 +2027,115 @@
          ENDIF
          
+      ENDIF
+
+      !  - log10(chlorophyll)
+      IF ( ln_logchl ) THEN
+
+         ! Copy data from logchldatqc to logchldata structures
+         DO jlogchlset = 1, nlogchlsets
+
+            CALL obs_surf_decompress( logchldatqc(jlogchlset), &
+                 &                    logchldata(jlogchlset), .TRUE., numout )
+
+         END DO
+         
+         ! Mark as bad observations with no valid model counterpart due to activities in dia_obs
+         ! Seem to need to set to fill value rather than marking as bad to be effective, so do both
+         DO jlogchlset = 1, nlogchlsets
+            WHERE ( logchldata(jlogchlset)%rmod(:,1) == obfillflt )
+               logchldata(jlogchlset)%nqc(:)    = 1
+               logchldata(jlogchlset)%robs(:,1) = obfillflt
+            END WHERE
+         END DO
+
+         ! Write the logchl data
+         DO jlogchlset = 1, nlogchlsets
+      
+            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
+
+      ENDIF
+
+      !  - fco2
+      IF ( ln_fco2 ) THEN
+
+         ! Copy data from fco2datqc to fco2data structures
+         DO jfco2set = 1, nfco2sets
+
+            CALL obs_surf_decompress( fco2datqc(jfco2set), &
+                 &                    fco2data(jfco2set), .TRUE., numout )
+
+         END DO
+         
+         ! Mark as bad observations with no valid model counterpart due to fCO2 not being in the restart
+         ! Seem to need to set to fill value rather than marking as bad to be effective, so do both
+         DO jfco2set = 1, nfco2sets
+            WHERE ( fco2data(jfco2set)%rmod(:,1) == obfillflt )
+               fco2data(jfco2set)%nqc(:)    = 1
+               fco2data(jfco2set)%robs(:,1) = obfillflt
+            END WHERE
+         END DO
+
+         ! Write the fco2 data
+         DO jfco2set = 1, nfco2sets
+      
+            WRITE(cdtmp,'(A,I2.2)')'fco2fb_',jfco2set
+            CALL obs_wri_fco2( cdtmp, fco2data(jfco2set) )
+
+         END DO
+
+      ENDIF
+
+      !  - pco2
+      IF ( ln_pco2 ) THEN
+
+         ! Copy data from pco2datqc to pco2data structures
+         DO jpco2set = 1, npco2sets
+
+            CALL obs_surf_decompress( pco2datqc(jpco2set), &
+                 &                    pco2data(jpco2set), .TRUE., numout )
+
+         END DO
+         
+         ! Mark as bad observations with no valid model counterpart due to pco2 not being in the restart
+         ! Seem to need to set to fill value rather than marking as bad to be effective, so do both
+         DO jpco2set = 1, npco2sets
+            WHERE ( pco2data(jpco2set)%rmod(:,1) == obfillflt )
+               pco2data(jpco2set)%nqc(:)    = 1
+               pco2data(jpco2set)%robs(:,1) = obfillflt
+            END WHERE
+         END DO
+
+         ! Write the pco2 data
+         DO jpco2set = 1, npco2sets
+      
+            WRITE(cdtmp,'(A,I2.2)')'pco2fb_',jpco2set
+            CALL obs_wri_pco2( cdtmp, pco2data(jpco2set) )
+
+         END DO
+
       ENDIF
 
Index: /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_fco2.F90
===================================================================
--- /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_fco2.F90	(revision 7713)
+++ /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_fco2.F90	(revision 7713)
@@ -0,0 +1,35 @@
+MODULE obs_fco2
+   !!=====================================================================
+   !!                       ***  MODULE  obs_fco2  ***
+   !! Observation diagnostics: Storage space for fco2 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 nfco2vars, nfco2extr, nfco2sets, fco2data, fco2datqc
+
+   !! * Shared Module variables
+   INTEGER :: nfco2vars                               ! Number of fco2data variables
+   INTEGER :: nfco2extr                               ! Number of fco2data extra 
+                                                      ! variables
+   INTEGER :: nfco2sets                               ! Number of fco2data sets
+   TYPE(obs_surf), POINTER, DIMENSION(:) :: fco2data  ! Initial fco2 data
+   TYPE(obs_surf), POINTER, DIMENSION(:) :: fco2datqc ! Sea ice data after quality control
+
+END MODULE obs_fco2
+
Index: /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_fco2_io.F90
===================================================================
--- /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_fco2_io.F90	(revision 7713)
+++ /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_fco2_io.F90	(revision 7713)
@@ -0,0 +1,30 @@
+MODULE obs_fco2_io
+   !!======================================================================
+   !!                       ***  MODULE obs_fco2_io  ***
+   !! Observation operators : I/O for fco2 files
+   !!======================================================================
+   !! History : 
+   !!             !  09-01  (K. Mogensen) Initial version
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   read_fco2file    :  Read a obfbdata structure from a fco2 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 "obsfco2_io.h90"
+
+END MODULE obs_fco2_io
Index: /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_logchl.F90
===================================================================
--- /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_logchl.F90	(revision 7713)
+++ /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_logchl.F90	(revision 7713)
@@ -0,0 +1,35 @@
+MODULE obs_logchl
+   !!=====================================================================
+   !!                       ***  MODULE  obs_logchl  ***
+   !! Observation diagnostics: Storage space for logchl 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 nlogchlvars, nlogchlextr, nlogchlsets, logchldata, logchldatqc
+
+   !! * Shared Module variables
+   INTEGER :: nlogchlvars                               ! Number of logchldata variables
+   INTEGER :: nlogchlextr                               ! Number of logchldata extra 
+                                                        ! variables
+   INTEGER :: nlogchlsets                               ! Number of logchldata sets
+   TYPE(obs_surf), POINTER, DIMENSION(:) :: logchldata  ! Initial logchl data
+   TYPE(obs_surf), POINTER, DIMENSION(:) :: logchldatqc ! Sea ice data after quality control
+
+END MODULE obs_logchl
+
Index: /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_logchl_io.F90
===================================================================
--- /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_logchl_io.F90	(revision 7713)
+++ /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_logchl_io.F90	(revision 7713)
@@ -0,0 +1,30 @@
+MODULE obs_logchl_io
+   !!======================================================================
+   !!                       ***  MODULE obs_logchl_io  ***
+   !! Observation operators : I/O for logchl files
+   !!======================================================================
+   !! History : 
+   !!             !  09-01  (K. Mogensen) Initial version
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   read_logchlfile    :  Read a obfbdata structure from a logchl 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 "obslogchl_io.h90"
+
+END MODULE obs_logchl_io
Index: /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
===================================================================
--- /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 7712)
+++ /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 7713)
@@ -23,4 +23,12 @@
    !!   obs_vel_opt :    Compute the model counterpart of zonal and meridional
    !!                    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
+   !!   obs_fco2_opt :   Compute the model counterpart of fco2
+   !!                    observations
+   !!   obs_pco2_opt :   Compute the model counterpart of pco2
+   !!                    observations
    !!----------------------------------------------------------------------
 
@@ -63,5 +71,9 @@
       &   obs_sss_opt, &  ! Compute the model counterpart of SSS observations
       &   obs_seaice_opt, &
-      &   obs_vel_opt     ! Compute the model counterpart of velocity profile data
+      &   obs_vel_opt, &  ! Compute the model counterpart of velocity profile data
+      &   obs_logchl_opt, & ! Compute the model counterpart of logchl data
+      &   obs_spm_opt, &  ! Compute the model counterpart of spm data
+      &   obs_fco2_opt, & ! Compute the model counterpart of fco2 data
+      &   obs_pco2_opt    ! Compute the model counterpart of pco2 data
 
    INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types
@@ -2052,4 +2064,644 @@
    END SUBROUTINE obs_vel_opt
 
+   SUBROUTINE obs_logchl_opt( logchldatqc, kt, kpi, kpj, kit000, &
+      &                    plogchln, plogchlmask, k2dint )
+
+      !!-----------------------------------------------------------------------
+      !!
+      !!                     ***  ROUTINE obs_logchl_opt  ***
+      !!
+      !! ** Purpose : Compute the model counterpart of logchl
+      !!              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 logchl 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) :: logchldatqc     ! 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) :: &
+         & plogchln,  &    ! Model logchl field
+         & plogchlmask     ! Land-sea mask
+         
+      !! * Local declarations
+      INTEGER :: ji
+      INTEGER :: jj
+      INTEGER :: jobs
+      INTEGER :: inrc
+      INTEGER :: ilogchl
+      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, &
+         & zlogchll, &
+         & zglam, &
+         & zgphi
+      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
+         & igrdi, &
+         & igrdj
+
+      !------------------------------------------------------------------------
+      ! Local initialization 
+      !------------------------------------------------------------------------
+      ! ... Record and data counters
+      inrc = kt - kit000 + 2
+      ilogchl = logchldatqc%nsstp(inrc)
+
+      ! Get the data for interpolation
+      
+      ALLOCATE( &
+         & igrdi(2,2,ilogchl), &
+         & igrdj(2,2,ilogchl), &
+         & zglam(2,2,ilogchl), &
+         & zgphi(2,2,ilogchl), &
+         & zmask(2,2,ilogchl), &
+         & zlogchll(2,2,ilogchl)  &
+         & )
+      
+      DO jobs = logchldatqc%nsurfup + 1, logchldatqc%nsurfup + ilogchl
+         iobs = jobs - logchldatqc%nsurfup
+         igrdi(1,1,iobs) = logchldatqc%mi(jobs)-1
+         igrdj(1,1,iobs) = logchldatqc%mj(jobs)-1
+         igrdi(1,2,iobs) = logchldatqc%mi(jobs)-1
+         igrdj(1,2,iobs) = logchldatqc%mj(jobs)
+         igrdi(2,1,iobs) = logchldatqc%mi(jobs)
+         igrdj(2,1,iobs) = logchldatqc%mj(jobs)-1
+         igrdi(2,2,iobs) = logchldatqc%mi(jobs)
+         igrdj(2,2,iobs) = logchldatqc%mj(jobs)
+      END DO
+      
+      CALL obs_int_comm_2d( 2, 2, ilogchl, &
+         &                  igrdi, igrdj, glamt, zglam )
+      CALL obs_int_comm_2d( 2, 2, ilogchl, &
+         &                  igrdi, igrdj, gphit, zgphi )
+      CALL obs_int_comm_2d( 2, 2, ilogchl, &
+         &                  igrdi, igrdj, plogchlmask, zmask )
+      CALL obs_int_comm_2d( 2, 2, ilogchl, &
+         &                  igrdi, igrdj, plogchln, zlogchll )
+      
+      DO jobs = logchldatqc%nsurfup + 1, logchldatqc%nsurfup + ilogchl
+         
+         iobs = jobs - logchldatqc%nsurfup
+         
+         IF ( kt /= logchldatqc%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    = ', logchldatqc%mstp(jobs), &
+                  &            ' ntyp    = ', logchldatqc%ntyp(jobs)
+            ENDIF
+            CALL ctl_stop( 'obs_logchl_opt', 'Inconsistent time' )
+            
+         ENDIF
+         
+         zlam = logchldatqc%rlam(jobs)
+         zphi = logchldatqc%rphi(jobs)
+         
+         ! Get weights to interpolate the model logchl 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 logchl to the observation point
+         CALL obs_int_h2d( 1, 1,      &
+            &              zweig, zlogchll(:,:,iobs),  zext )
+         
+         logchldatqc%rmod(jobs,1) = zext(1)
+         
+      END DO
+      
+      ! Deallocate the data for interpolation
+      DEALLOCATE( &
+         & igrdi,    &
+         & igrdj,    &
+         & zglam,    &
+         & zgphi,    &
+         & zmask,    &
+         & zlogchll  &
+         & )
+      
+      logchldatqc%nsurfup = logchldatqc%nsurfup + ilogchl
+
+   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
+
+   SUBROUTINE obs_fco2_opt( fco2datqc, kt, kpi, kpj, kit000, &
+      &                    pfco2n, pfco2mask, k2dint )
+
+      !!-----------------------------------------------------------------------
+      !!
+      !!                     ***  ROUTINE obs_fco2_opt  ***
+      !!
+      !! ** Purpose : Compute the model counterpart of fco2
+      !!              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 fco2 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) :: fco2datqc     ! 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) :: &
+         & pfco2n,  &    ! Model fco2 field
+         & pfco2mask     ! Land-sea mask
+         
+      !! * Local declarations
+      INTEGER :: ji
+      INTEGER :: jj
+      INTEGER :: jobs
+      INTEGER :: inrc
+      INTEGER :: ifco2
+      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, &
+         & zfco2l, &
+         & zglam, &
+         & zgphi
+      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
+         & igrdi, &
+         & igrdj
+
+      !------------------------------------------------------------------------
+      ! Local initialization 
+      !------------------------------------------------------------------------
+      ! ... Record and data counters
+      inrc = kt - kit000 + 2
+      ifco2 = fco2datqc%nsstp(inrc)
+
+      ! Get the data for interpolation
+      
+      ALLOCATE( &
+         & igrdi(2,2,ifco2), &
+         & igrdj(2,2,ifco2), &
+         & zglam(2,2,ifco2), &
+         & zgphi(2,2,ifco2), &
+         & zmask(2,2,ifco2), &
+         & zfco2l(2,2,ifco2)  &
+         & )
+      
+      DO jobs = fco2datqc%nsurfup + 1, fco2datqc%nsurfup + ifco2
+         iobs = jobs - fco2datqc%nsurfup
+         igrdi(1,1,iobs) = fco2datqc%mi(jobs)-1
+         igrdj(1,1,iobs) = fco2datqc%mj(jobs)-1
+         igrdi(1,2,iobs) = fco2datqc%mi(jobs)-1
+         igrdj(1,2,iobs) = fco2datqc%mj(jobs)
+         igrdi(2,1,iobs) = fco2datqc%mi(jobs)
+         igrdj(2,1,iobs) = fco2datqc%mj(jobs)-1
+         igrdi(2,2,iobs) = fco2datqc%mi(jobs)
+         igrdj(2,2,iobs) = fco2datqc%mj(jobs)
+      END DO
+      
+      CALL obs_int_comm_2d( 2, 2, ifco2, &
+         &                  igrdi, igrdj, glamt, zglam )
+      CALL obs_int_comm_2d( 2, 2, ifco2, &
+         &                  igrdi, igrdj, gphit, zgphi )
+      CALL obs_int_comm_2d( 2, 2, ifco2, &
+         &                  igrdi, igrdj, pfco2mask, zmask )
+      CALL obs_int_comm_2d( 2, 2, ifco2, &
+         &                  igrdi, igrdj, pfco2n, zfco2l )
+      
+      DO jobs = fco2datqc%nsurfup + 1, fco2datqc%nsurfup + ifco2
+         
+         iobs = jobs - fco2datqc%nsurfup
+         
+         IF ( kt /= fco2datqc%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    = ', fco2datqc%mstp(jobs), &
+                  &            ' ntyp    = ', fco2datqc%ntyp(jobs)
+            ENDIF
+            CALL ctl_stop( 'obs_fco2_opt', 'Inconsistent time' )
+            
+         ENDIF
+         
+         zlam = fco2datqc%rlam(jobs)
+         zphi = fco2datqc%rphi(jobs)
+         
+         ! Get weights to interpolate the model fco2 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 fco2 to the observation point
+         CALL obs_int_h2d( 1, 1,      &
+            &              zweig, zfco2l(:,:,iobs),  zext )
+         
+         fco2datqc%rmod(jobs,1) = zext(1)
+         
+      END DO
+      
+      ! Deallocate the data for interpolation
+      DEALLOCATE( &
+         & igrdi,    &
+         & igrdj,    &
+         & zglam,    &
+         & zgphi,    &
+         & zmask,    &
+         & zfco2l  &
+         & )
+      
+      fco2datqc%nsurfup = fco2datqc%nsurfup + ifco2
+
+   END SUBROUTINE obs_fco2_opt
+
+   SUBROUTINE obs_pco2_opt( pco2datqc, kt, kpi, kpj, kit000, &
+      &                    ppco2n, ppco2mask, k2dint )
+
+      !!-----------------------------------------------------------------------
+      !!
+      !!                     ***  ROUTINE obs_pco2_opt  ***
+      !!
+      !! ** Purpose : Compute the model counterpart of pco2
+      !!              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 pco2 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) :: pco2datqc     ! 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) :: &
+         & ppco2n,  &    ! Model pco2 field
+         & ppco2mask     ! Land-sea mask
+         
+      !! * Local declarations
+      INTEGER :: ji
+      INTEGER :: jj
+      INTEGER :: jobs
+      INTEGER :: inrc
+      INTEGER :: ipco2
+      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, &
+         & zpco2l, &
+         & zglam, &
+         & zgphi
+      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
+         & igrdi, &
+         & igrdj
+
+      !------------------------------------------------------------------------
+      ! Local initialization 
+      !------------------------------------------------------------------------
+      ! ... Record and data counters
+      inrc = kt - kit000 + 2
+      ipco2 = pco2datqc%nsstp(inrc)
+
+      ! Get the data for interpolation
+      
+      ALLOCATE( &
+         & igrdi(2,2,ipco2), &
+         & igrdj(2,2,ipco2), &
+         & zglam(2,2,ipco2), &
+         & zgphi(2,2,ipco2), &
+         & zmask(2,2,ipco2), &
+         & zpco2l(2,2,ipco2)  &
+         & )
+      
+      DO jobs = pco2datqc%nsurfup + 1, pco2datqc%nsurfup + ipco2
+         iobs = jobs - pco2datqc%nsurfup
+         igrdi(1,1,iobs) = pco2datqc%mi(jobs)-1
+         igrdj(1,1,iobs) = pco2datqc%mj(jobs)-1
+         igrdi(1,2,iobs) = pco2datqc%mi(jobs)-1
+         igrdj(1,2,iobs) = pco2datqc%mj(jobs)
+         igrdi(2,1,iobs) = pco2datqc%mi(jobs)
+         igrdj(2,1,iobs) = pco2datqc%mj(jobs)-1
+         igrdi(2,2,iobs) = pco2datqc%mi(jobs)
+         igrdj(2,2,iobs) = pco2datqc%mj(jobs)
+      END DO
+      
+      CALL obs_int_comm_2d( 2, 2, ipco2, &
+         &                  igrdi, igrdj, glamt, zglam )
+      CALL obs_int_comm_2d( 2, 2, ipco2, &
+         &                  igrdi, igrdj, gphit, zgphi )
+      CALL obs_int_comm_2d( 2, 2, ipco2, &
+         &                  igrdi, igrdj, ppco2mask, zmask )
+      CALL obs_int_comm_2d( 2, 2, ipco2, &
+         &                  igrdi, igrdj, ppco2n, zpco2l )
+      
+      DO jobs = pco2datqc%nsurfup + 1, pco2datqc%nsurfup + ipco2
+         
+         iobs = jobs - pco2datqc%nsurfup
+         
+         IF ( kt /= pco2datqc%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    = ', pco2datqc%mstp(jobs), &
+                  &            ' ntyp    = ', pco2datqc%ntyp(jobs)
+            ENDIF
+            CALL ctl_stop( 'obs_pco2_opt', 'Inconsistent time' )
+            
+         ENDIF
+         
+         zlam = pco2datqc%rlam(jobs)
+         zphi = pco2datqc%rphi(jobs)
+         
+         ! Get weights to interpolate the model pco2 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 pco2 to the observation point
+         CALL obs_int_h2d( 1, 1,      &
+            &              zweig, zpco2l(:,:,iobs),  zext )
+         
+         pco2datqc%rmod(jobs,1) = zext(1)
+         
+      END DO
+      
+      ! Deallocate the data for interpolation
+      DEALLOCATE( &
+         & igrdi,    &
+         & igrdj,    &
+         & zglam,    &
+         & zgphi,    &
+         & zmask,    &
+         & zpco2l  &
+         & )
+      
+      pco2datqc%nsurfup = pco2datqc%nsurfup + ipco2
+
+   END SUBROUTINE obs_pco2_opt
+
 END MODULE obs_oper
 
Index: /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_pco2.F90
===================================================================
--- /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_pco2.F90	(revision 7713)
+++ /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_pco2.F90	(revision 7713)
@@ -0,0 +1,35 @@
+MODULE obs_pco2
+   !!=====================================================================
+   !!                       ***  MODULE  obs_pco2  ***
+   !! Observation diagnostics: Storage space for pco2 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 npco2vars, npco2extr, npco2sets, pco2data, pco2datqc
+
+   !! * Shared Module variables
+   INTEGER :: npco2vars                               ! Number of pco2data variables
+   INTEGER :: npco2extr                               ! Number of pco2data extra 
+                                                      ! variables
+   INTEGER :: npco2sets                               ! Number of pco2data sets
+   TYPE(obs_surf), POINTER, DIMENSION(:) :: pco2data  ! Initial pco2 data
+   TYPE(obs_surf), POINTER, DIMENSION(:) :: pco2datqc ! Sea ice data after quality control
+
+END MODULE obs_pco2
+
Index: /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_pco2_io.F90
===================================================================
--- /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_pco2_io.F90	(revision 7713)
+++ /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_pco2_io.F90	(revision 7713)
@@ -0,0 +1,30 @@
+MODULE obs_pco2_io
+   !!======================================================================
+   !!                       ***  MODULE obs_pco2_io  ***
+   !! Observation operators : I/O for pco2 files
+   !!======================================================================
+   !! History : 
+   !!             !  09-01  (K. Mogensen) Initial version
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   read_pco2file    :  Read a obfbdata structure from a pco2 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 "obspco2_io.h90"
+
+END MODULE obs_pco2_io
Index: /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
===================================================================
--- /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 7712)
+++ /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 7713)
@@ -12,4 +12,8 @@
    !!   obs_pre_seaice : First level check and screening of sea ice observations
    !!   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_pre_fco2 : First level check and screening of fco2 obs.
+   !!   obs_pre_pco2 : First level check and screening of pco2 obs.
    !!   obs_scr      : Basic screening of the observations
    !!   obs_coo_tim  : Compute number of time steps to the observation time
@@ -45,4 +49,8 @@
       & obs_pre_seaice, & ! First level check and screening of sea ice data
       & 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
+      & obs_pre_fco2, &   ! First level check and screening of fco2 data
+      & obs_pre_pco2, &   ! First level check and screening of pco2 data
       & calc_month_len     ! Calculate the number of days in the months of a year  
 
@@ -1376,4 +1384,768 @@
    END SUBROUTINE obs_pre_vel
 
+   SUBROUTINE obs_pre_logchl( logchldata, logchldatqc, ld_logchl, ld_nea )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE obs_pre_logchl  ***
+      !!
+      !! ** Purpose : First level check and screening of logchl observations
+      !!
+      !! ** Method  : First level check and screening of logchl 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) :: logchldata     ! Full set of logchl data
+      TYPE(obs_surf), INTENT(INOUT) :: logchldatqc    ! Subset of logchl data not failing screening
+      LOGICAL :: ld_logchl     ! Switch for logchl 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
+      INTEGER :: ibdysobs     !  - close to open boundary
+                              ! 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
+      INTEGER :: ibdysobsmpp  !  - close to open boundary
+      LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
+         & llvalid            ! data selection
+      INTEGER :: jobs         ! Obs. loop variable
+      INTEGER :: jstp         ! Time loop variable
+      INTEGER :: inrc         ! Time index variable
+      INTEGER :: irec         ! Record index
+
+      IF (lwp) WRITE(numout,*)'obs_pre_logchl : Preparing the logchl 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
+      ibdysobs = 0
+
+      ! -----------------------------------------------------------------------
+      ! Find time coordinate for logchl data
+      ! -----------------------------------------------------------------------
+
+      CALL obs_coo_tim( icycle, &
+         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
+         &              logchldata%nsurf,   logchldata%nyea, logchldata%nmon, &
+         &              logchldata%nday,    logchldata%nhou, logchldata%nmin, &
+         &              logchldata%nqc,     logchldata%mstp, iotdobs        )
+      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
+      ! -----------------------------------------------------------------------
+      ! Check for logchl data failing the grid search
+      ! -----------------------------------------------------------------------
+
+      CALL obs_coo_grd( logchldata%nsurf,   logchldata%mi, logchldata%mj, &
+         &              logchldata%nqc,     igrdobs                         )
+      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
+
+      ! -----------------------------------------------------------------------
+      ! Check for land points. 
+      ! -----------------------------------------------------------------------
+
+      CALL obs_coo_spc_2d( logchldata%nsurf,                 &
+         &                 jpi,             jpj,             &
+         &                 logchldata%mi,   logchldata%mj,   & 
+         &                 logchldata%rlam, logchldata%rphi, &
+         &                 glamt,           gphit,           &
+         &                 tmask(:,:,1),    logchldata%nqc,  &
+         &                 iosdsobs,        ilansobs,        &
+         &                 inlasobs,        ld_nea,          &
+         &                 ibdysobs,        ln_bound_reject  ) 
+         
+      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
+      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
+      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
+      CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp )
+
+      ! -----------------------------------------------------------------------
+      ! Copy useful data from the logchldata data structure to
+      ! the logchldatqc data structure 
+      ! -----------------------------------------------------------------------
+
+      ! Allocate the selection arrays
+
+      ALLOCATE( llvalid(logchldata%nsurf) )
+      
+      ! We want all data which has qc flags <= 0
+
+      llvalid(:)  = ( logchldata%nqc(:)  <= 10 )
+
+      ! The actual copying
+
+      CALL obs_surf_compress( logchldata,     logchldatqc,       .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_logchl :'
+         WRITE(numout,*) '~~~~~~~~~~~'
+         WRITE(numout,*)
+         WRITE(numout,*) ' logchl data outside time domain                  = ', &
+            &            iotdobsmpp
+         WRITE(numout,*) ' Remaining logchl data that failed grid search    = ', &
+            &            igrdobsmpp
+         WRITE(numout,*) ' Remaining logchl data outside space domain       = ', &
+            &            iosdsobsmpp
+         WRITE(numout,*) ' Remaining logchl data at land points             = ', &
+            &            ilansobsmpp
+         IF (ld_nea) THEN
+            WRITE(numout,*) ' Remaining logchl data near land points (removed) = ', &
+               &            inlasobsmpp
+         ELSE
+            WRITE(numout,*) ' Remaining logchl data near land points (kept)    = ', &
+               &            inlasobsmpp
+         ENDIF
+         WRITE(numout,*) ' Remaining logchl data near open boundary (removed) = ', &
+           &            ibdysobsmpp
+         WRITE(numout,*) ' logchl data accepted                             = ', &
+            &            logchldatqc%nsurfmpp
+
+         WRITE(numout,*)
+         WRITE(numout,*) ' Number of observations per time step :'
+         WRITE(numout,*)
+         WRITE(numout,1997)
+         WRITE(numout,1998)
+      ENDIF
+      
+      DO jobs = 1, logchldatqc%nsurf
+         inrc = logchldatqc%mstp(jobs) + 2 - nit000
+         logchldatqc%nsstp(inrc)  = logchldatqc%nsstp(inrc) + 1
+      END DO
+      
+      CALL obs_mpp_sum_integers( logchldatqc%nsstp, logchldatqc%nsstpmpp, &
+         &                       nitend - nit000 + 2 )
+
+      IF ( lwp ) THEN
+         DO jstp = nit000 - 1, nitend
+            inrc = jstp - nit000 + 2
+            WRITE(numout,1999) jstp, logchldatqc%nsstpmpp(inrc)
+         END DO
+      ENDIF
+
+1997  FORMAT(10X,'Time step',5X,'logchl data')
+1998  FORMAT(10X,'---------',5X,'------------')
+1999  FORMAT(10X,I9,5X,I17)
+      
+   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
+      INTEGER :: ibdysobs     !  - close to open boundary
+                              ! 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
+      INTEGER :: ibdysobsmpp  !  - close to open boundary
+      LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
+         & llvalid            ! data selection
+      INTEGER :: jobs         ! Obs. loop variable
+      INTEGER :: jstp         ! Time loop variable
+      INTEGER :: inrc         ! Time index variable
+      INTEGER :: irec         ! Record index
+
+      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
+      ibdysobs = 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,          &
+         &                 ibdysobs,        ln_bound_reject  ) 
+         
+      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
+      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
+      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
+      CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp )
+
+      ! -----------------------------------------------------------------------
+      ! 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,*) ' Remaining spm data near open boundary (removed) = ', &
+            &            ibdysobsmpp
+         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_pre_fco2( fco2data, fco2datqc, ld_fco2, ld_nea )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE obs_pre_fco2  ***
+      !!
+      !! ** Purpose : First level check and screening of fco2 observations
+      !!
+      !! ** Method  : First level check and screening of fco2 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) :: fco2data     ! Full set of fco2 data
+      TYPE(obs_surf), INTENT(INOUT) :: fco2datqc    ! Subset of fco2 data not failing screening
+      LOGICAL :: ld_fco2     ! Switch for fco2 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
+      INTEGER :: ibdysobs     !  - close to open boundary
+                              ! 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
+      INTEGER :: ibdysobsmpp  !  - close to open boundary
+      LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
+         & llvalid            ! data selection
+      INTEGER :: jobs         ! Obs. loop variable
+      INTEGER :: jstp         ! Time loop variable
+      INTEGER :: inrc         ! Time index variable
+      INTEGER :: irec         ! Record index
+
+      IF (lwp) WRITE(numout,*)'obs_pre_fco2 : Preparing the fco2 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
+      ibdysobs = 0
+
+      ! -----------------------------------------------------------------------
+      ! Find time coordinate for fco2 data
+      ! -----------------------------------------------------------------------
+
+      CALL obs_coo_tim( icycle, &
+         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
+         &              fco2data%nsurf,   fco2data%nyea, fco2data%nmon, &
+         &              fco2data%nday,    fco2data%nhou, fco2data%nmin, &
+         &              fco2data%nqc,     fco2data%mstp, iotdobs        )
+      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
+      ! -----------------------------------------------------------------------
+      ! Check for fco2 data failing the grid search
+      ! -----------------------------------------------------------------------
+
+      CALL obs_coo_grd( fco2data%nsurf,   fco2data%mi, fco2data%mj, &
+         &              fco2data%nqc,     igrdobs                         )
+      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
+
+      ! -----------------------------------------------------------------------
+      ! Check for land points. 
+      ! -----------------------------------------------------------------------
+
+      CALL obs_coo_spc_2d( fco2data%nsurf,                 &
+         &                 jpi,             jpj,             &
+         &                 fco2data%mi,   fco2data%mj,   & 
+         &                 fco2data%rlam, fco2data%rphi, &
+         &                 glamt,           gphit,           &
+         &                 tmask(:,:,1),    fco2data%nqc,  &
+         &                 iosdsobs,        ilansobs,        &
+         &                 inlasobs,        ld_nea,          &
+         &                 ibdysobs,        ln_bound_reject  ) 
+         
+      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
+      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
+      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
+      CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp )
+
+      ! -----------------------------------------------------------------------
+      ! Copy useful data from the fco2data data structure to
+      ! the fco2datqc data structure 
+      ! -----------------------------------------------------------------------
+
+      ! Allocate the selection arrays
+
+      ALLOCATE( llvalid(fco2data%nsurf) )
+      
+      ! We want all data which has qc flags <= 0
+
+      llvalid(:)  = ( fco2data%nqc(:)  <= 10 )
+
+      ! The actual copying
+
+      CALL obs_surf_compress( fco2data,     fco2datqc,       .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_fco2 :'
+         WRITE(numout,*) '~~~~~~~~~~~'
+         WRITE(numout,*)
+         WRITE(numout,*) ' fco2 data outside time domain                  = ', &
+            &            iotdobsmpp
+         WRITE(numout,*) ' Remaining fco2 data that failed grid search    = ', &
+            &            igrdobsmpp
+         WRITE(numout,*) ' Remaining fco2 data outside space domain       = ', &
+            &            iosdsobsmpp
+         WRITE(numout,*) ' Remaining fco2 data at land points             = ', &
+            &            ilansobsmpp
+         IF (ld_nea) THEN
+            WRITE(numout,*) ' Remaining fco2 data near land points (removed) = ', &
+               &            inlasobsmpp
+         ELSE
+            WRITE(numout,*) ' Remaining fco2 data near land points (kept)    = ', &
+               &            inlasobsmpp
+         ENDIF
+         WRITE(numout,*) ' Remaining fco2 data near open boundary (removed) = ', &
+           &            ibdysobsmpp
+         WRITE(numout,*) ' fco2 data accepted                             = ', &
+            &            fco2datqc%nsurfmpp
+
+         WRITE(numout,*)
+         WRITE(numout,*) ' Number of observations per time step :'
+         WRITE(numout,*)
+         WRITE(numout,1997)
+         WRITE(numout,1998)
+      ENDIF
+      
+      DO jobs = 1, fco2datqc%nsurf
+         inrc = fco2datqc%mstp(jobs) + 2 - nit000
+         fco2datqc%nsstp(inrc)  = fco2datqc%nsstp(inrc) + 1
+      END DO
+      
+      CALL obs_mpp_sum_integers( fco2datqc%nsstp, fco2datqc%nsstpmpp, &
+         &                       nitend - nit000 + 2 )
+
+      IF ( lwp ) THEN
+         DO jstp = nit000 - 1, nitend
+            inrc = jstp - nit000 + 2
+            WRITE(numout,1999) jstp, fco2datqc%nsstpmpp(inrc)
+         END DO
+      ENDIF
+
+1997  FORMAT(10X,'Time step',5X,'fco2 data')
+1998  FORMAT(10X,'---------',5X,'------------')
+1999  FORMAT(10X,I9,5X,I17)
+      
+   END SUBROUTINE obs_pre_fco2
+
+   SUBROUTINE obs_pre_pco2( pco2data, pco2datqc, ld_pco2, ld_nea )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE obs_pre_pco2  ***
+      !!
+      !! ** Purpose : First level check and screening of pco2 observations
+      !!
+      !! ** Method  : First level check and screening of pco2 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) :: pco2data     ! Full set of pco2 data
+      TYPE(obs_surf), INTENT(INOUT) :: pco2datqc    ! Subset of pco2 data not failing screening
+      LOGICAL :: ld_pco2     ! Switch for pco2 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
+      INTEGER :: ibdysobs     !  - close to open boundary
+                              ! 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
+      INTEGER :: ibdysobsmpp  !  - close to open boundary
+      LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
+         & llvalid            ! data selection
+      INTEGER :: jobs         ! Obs. loop variable
+      INTEGER :: jstp         ! Time loop variable
+      INTEGER :: inrc         ! Time index variable
+      INTEGER :: irec         ! Record index
+
+      IF (lwp) WRITE(numout,*)'obs_pre_pco2 : Preparing the pco2 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
+      ibdysobs = 0
+
+      ! -----------------------------------------------------------------------
+      ! Find time coordinate for pco2 data
+      ! -----------------------------------------------------------------------
+
+      CALL obs_coo_tim( icycle, &
+         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
+         &              pco2data%nsurf,   pco2data%nyea, pco2data%nmon, &
+         &              pco2data%nday,    pco2data%nhou, pco2data%nmin, &
+         &              pco2data%nqc,     pco2data%mstp, iotdobs        )
+      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
+      ! -----------------------------------------------------------------------
+      ! Check for pco2 data failing the grid search
+      ! -----------------------------------------------------------------------
+
+      CALL obs_coo_grd( pco2data%nsurf,   pco2data%mi, pco2data%mj, &
+         &              pco2data%nqc,     igrdobs                         )
+      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
+
+      ! -----------------------------------------------------------------------
+      ! Check for land points. 
+      ! -----------------------------------------------------------------------
+
+      CALL obs_coo_spc_2d( pco2data%nsurf,                 &
+         &                 jpi,             jpj,             &
+         &                 pco2data%mi,   pco2data%mj,   & 
+         &                 pco2data%rlam, pco2data%rphi, &
+         &                 glamt,           gphit,           &
+         &                 tmask(:,:,1),    pco2data%nqc,  &
+         &                 iosdsobs,        ilansobs,        &
+         &                 inlasobs,        ld_nea,          &
+         &                 ibdysobs,        ln_bound_reject  ) 
+         
+      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
+      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
+      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
+      CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp )
+
+      ! -----------------------------------------------------------------------
+      ! Copy useful data from the pco2data data structure to
+      ! the pco2datqc data structure 
+      ! -----------------------------------------------------------------------
+
+      ! Allocate the selection arrays
+
+      ALLOCATE( llvalid(pco2data%nsurf) )
+      
+      ! We want all data which has qc flags <= 0
+
+      llvalid(:)  = ( pco2data%nqc(:)  <= 10 )
+
+      ! The actual copying
+
+      CALL obs_surf_compress( pco2data,     pco2datqc,       .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_pco2 :'
+         WRITE(numout,*) '~~~~~~~~~~~'
+         WRITE(numout,*)
+         WRITE(numout,*) ' pco2 data outside time domain                  = ', &
+            &            iotdobsmpp
+         WRITE(numout,*) ' Remaining pco2 data that failed grid search    = ', &
+            &            igrdobsmpp
+         WRITE(numout,*) ' Remaining pco2 data outside space domain       = ', &
+            &            iosdsobsmpp
+         WRITE(numout,*) ' Remaining pco2 data at land points             = ', &
+            &            ilansobsmpp
+         IF (ld_nea) THEN
+            WRITE(numout,*) ' Remaining pco2 data near land points (removed) = ', &
+               &            inlasobsmpp
+         ELSE
+            WRITE(numout,*) ' Remaining pco2 data near land points (kept)    = ', &
+               &            inlasobsmpp
+         ENDIF
+         WRITE(numout,*) ' Remaining pco2 data near open boundary (removed) = ', &
+           &            ibdysobsmpp
+         WRITE(numout,*) ' pco2 data accepted                             = ', &
+            &            pco2datqc%nsurfmpp
+
+         WRITE(numout,*)
+         WRITE(numout,*) ' Number of observations per time step :'
+         WRITE(numout,*)
+         WRITE(numout,1997)
+         WRITE(numout,1998)
+      ENDIF
+      
+      DO jobs = 1, pco2datqc%nsurf
+         inrc = pco2datqc%mstp(jobs) + 2 - nit000
+         pco2datqc%nsstp(inrc)  = pco2datqc%nsstp(inrc) + 1
+      END DO
+      
+      CALL obs_mpp_sum_integers( pco2datqc%nsstp, pco2datqc%nsstpmpp, &
+         &                       nitend - nit000 + 2 )
+
+      IF ( lwp ) THEN
+         DO jstp = nit000 - 1, nitend
+            inrc = jstp - nit000 + 2
+            WRITE(numout,1999) jstp, pco2datqc%nsstpmpp(inrc)
+         END DO
+      ENDIF
+
+1997  FORMAT(10X,'Time step',5X,'pco2 data')
+1998  FORMAT(10X,'---------',5X,'------------')
+1999  FORMAT(10X,I9,5X,I17)
+      
+   END SUBROUTINE obs_pre_pco2
+
    SUBROUTINE obs_coo_tim( kcycle, &
       &                    kyea0,   kmon0,   kday0,   khou0,   kmin0,     &
Index: /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_fco2.F90
===================================================================
--- /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_fco2.F90	(revision 7713)
+++ /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_fco2.F90	(revision 7713)
@@ -0,0 +1,465 @@
+MODULE obs_read_fco2
+   !!======================================================================
+   !!                       ***  MODULE obs_read_fco2  ***
+   !! Observation diagnostics: Read the along track fco2 data from
+   !!                          GHRSST or any fco2 data from feedback files
+   !!======================================================================
+
+   !!----------------------------------------------------------------------
+   !!   obs_rea_fco2 : Driver for reading fco2 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_fco2_io            ! I/O for fco2 files
+   USE iom                      ! I/O
+   USE netcdf                   ! NetCDF library
+
+   IMPLICIT NONE
+
+   !! * Routine accessibility
+   PRIVATE
+
+   PUBLIC obs_rea_fco2      ! Read the fco2 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_fco2( kformat, &
+      &                    fco2data, knumfiles, cfilenames, &
+      &                    kvars, kextr, kstp, ddobsini, ddobsend, &
+      &                    ldignmis, ldmod )
+      !!---------------------------------------------------------------------
+      !!
+      !!                   *** ROUTINE obs_rea_fco2 ***
+      !!
+      !! ** Purpose : Read from file the fco2 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) :: &
+         & fco2data     ! fco2 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 fco2data
+      INTEGER, INTENT(IN) :: kextr    ! Number of extra fields for each var in fco2data
+      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_fco2'
+      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 :: ifco2maxtype = 1024
+      INTEGER, DIMENSION(0:ifco2maxtype) :: &
+         & ityp, &
+         & itypmpp
+      INTEGER, DIMENSION(:), ALLOCATABLE :: &
+         & iobsi,    &
+         & iobsj,    &
+         & iproc,    &
+         & iindx,    &
+         & ifileidx, &
+         & ifco2idx
+      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) )
+
+      fco2_files : DO jj = 1, inobf
+          
+         !---------------------------------------------------------------------
+         ! Prints
+         !---------------------------------------------------------------------
+         IF(lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) ' obs_rea_fco2 : 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_fco2( 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 fco2_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), &
+         &      ifco2idx(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
+               ifco2idx(jk) = ji
+               zdat(jk)     = inpfiles(jj)%ptim(ji)
+            ENDIF
+         END DO
+      END DO
+      CALL sort_dp_indx( iobstot, &
+         &               zdat,     &
+         &               iindx   )
+      
+      CALL obs_surf_alloc( fco2data, iobs, & 
+                           kvars, kextr, kstp, jpi, jpj )
+      
+      ! * Read obs/positions, QC, all variable and assign to fco2data
+ 
+      iobs = 0
+
+      ityp   (:) = 0
+      itypmpp(:) = 0
+
+      ioserrcount=0      
+
+      DO jk = 1, iobstot
+         
+         jj = ifileidx(iindx(jk))
+         ji = ifco2idx(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) )
+
+
+               ! fco2 time coordinates
+               fco2data%nyea(iobs) = iyea
+               fco2data%nmon(iobs) = imon
+               fco2data%nday(iobs) = iday
+               fco2data%nhou(iobs) = ihou
+               fco2data%nmin(iobs) = imin
+               
+               ! fco2 space coordinates
+               fco2data%rlam(iobs) = inpfiles(jj)%plam(ji)
+               fco2data%rphi(iobs) = inpfiles(jj)%pphi(ji)
+
+               ! Coordinate search parameters
+               fco2data%mi  (iobs) = inpfiles(jj)%iobsi(ji,1)
+               fco2data%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
+               fco2data%ntyp(iobs) = itype
+               IF ( itype < ifco2maxtype + 1 ) THEN
+                  ityp(itype+1) = ityp(itype+1) + 1
+               ELSE
+                  IF(lwp)WRITE(numout,*)'WARNING:Increase ifco2maxtype in ',&
+                     &                  cpname
+               ENDIF
+
+               ! Bookkeeping data to match observations
+               fco2data%nsidx(iobs) = iobs
+               fco2data%nsfil(iobs) = iindx(jk)
+
+               ! QC flags
+               fco2data%nqc(iobs) = inpfiles(jj)%ivqc(ji,1)
+
+               ! Observed value
+               fco2data%robs(iobs,1) = inpfiles(jj)%pob(1,ji,1)
+
+
+               ! Model and MDT is set to fbrmdi unless read from file
+               IF ( ldmod ) THEN
+                  fco2data%rmod(iobs,1) = inpfiles(jj)%padd(1,ji,1,1)
+               ELSE
+                  fco2data%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)')'fco2 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, ifco2idx, zdat )
+
+      !-----------------------------------------------------------------------
+      ! Deallocate input data
+      !-----------------------------------------------------------------------
+      DO jj = 1, inobf
+         CALL dealloc_obfbdata( inpfiles(jj) )
+      END DO
+      DEALLOCATE( inpfiles )
+
+   END SUBROUTINE obs_rea_fco2
+
+END MODULE obs_read_fco2
+
Index: /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_logchl.F90
===================================================================
--- /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_logchl.F90	(revision 7713)
+++ /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_logchl.F90	(revision 7713)
@@ -0,0 +1,465 @@
+MODULE obs_read_logchl
+   !!======================================================================
+   !!                       ***  MODULE obs_read_logchl  ***
+   !! Observation diagnostics: Read the along track logchl data from
+   !!                          GHRSST or any logchl data from feedback files
+   !!======================================================================
+
+   !!----------------------------------------------------------------------
+   !!   obs_rea_logchl : Driver for reading logchl 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_logchl_io            ! I/O for logchl files
+   USE iom                      ! I/O
+   USE netcdf                   ! NetCDF library
+
+   IMPLICIT NONE
+
+   !! * Routine accessibility
+   PRIVATE
+
+   PUBLIC obs_rea_logchl      ! Read the logchl 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_logchl( kformat, &
+      &                    logchldata, knumfiles, cfilenames, &
+      &                    kvars, kextr, kstp, ddobsini, ddobsend, &
+      &                    ldignmis, ldmod )
+      !!---------------------------------------------------------------------
+      !!
+      !!                   *** ROUTINE obs_rea_logchl ***
+      !!
+      !! ** Purpose : Read from file the logchl 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) :: &
+         & logchldata     ! logchl 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 logchldata
+      INTEGER, INTENT(IN) :: kextr    ! Number of extra fields for each var in logchldata
+      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_logchl'
+      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 :: ilogchlmaxtype = 1024
+      INTEGER, DIMENSION(0:ilogchlmaxtype) :: &
+         & ityp, &
+         & itypmpp
+      INTEGER, DIMENSION(:), ALLOCATABLE :: &
+         & iobsi,    &
+         & iobsj,    &
+         & iproc,    &
+         & iindx,    &
+         & ifileidx, &
+         & ilogchlidx
+      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) )
+
+      logchl_files : DO jj = 1, inobf
+          
+         !---------------------------------------------------------------------
+         ! Prints
+         !---------------------------------------------------------------------
+         IF(lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) ' obs_rea_logchl : 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_logchl( 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 logchl_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), &
+         &      ilogchlidx(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
+               ilogchlidx(jk) = ji
+               zdat(jk)     = inpfiles(jj)%ptim(ji)
+            ENDIF
+         END DO
+      END DO
+      CALL sort_dp_indx( iobstot, &
+         &               zdat,     &
+         &               iindx   )
+      
+      CALL obs_surf_alloc( logchldata, iobs, & 
+                           kvars, kextr, kstp, jpi, jpj )
+      
+      ! * Read obs/positions, QC, all variable and assign to logchldata
+ 
+      iobs = 0
+
+      ityp   (:) = 0
+      itypmpp(:) = 0
+
+      ioserrcount=0      
+
+      DO jk = 1, iobstot
+         
+         jj = ifileidx(iindx(jk))
+         ji = ilogchlidx(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) )
+
+
+               ! logchl time coordinates
+               logchldata%nyea(iobs) = iyea
+               logchldata%nmon(iobs) = imon
+               logchldata%nday(iobs) = iday
+               logchldata%nhou(iobs) = ihou
+               logchldata%nmin(iobs) = imin
+               
+               ! logchl space coordinates
+               logchldata%rlam(iobs) = inpfiles(jj)%plam(ji)
+               logchldata%rphi(iobs) = inpfiles(jj)%pphi(ji)
+
+               ! Coordinate search parameters
+               logchldata%mi  (iobs) = inpfiles(jj)%iobsi(ji,1)
+               logchldata%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
+               logchldata%ntyp(iobs) = itype
+               IF ( itype < ilogchlmaxtype + 1 ) THEN
+                  ityp(itype+1) = ityp(itype+1) + 1
+               ELSE
+                  IF(lwp)WRITE(numout,*)'WARNING:Increase ilogchlmaxtype in ',&
+                     &                  cpname
+               ENDIF
+
+               ! Bookkeeping data to match observations
+               logchldata%nsidx(iobs) = iobs
+               logchldata%nsfil(iobs) = iindx(jk)
+
+               ! QC flags
+               logchldata%nqc(iobs) = inpfiles(jj)%ivqc(ji,1)
+
+               ! Observed value
+               logchldata%robs(iobs,1) = inpfiles(jj)%pob(1,ji,1)
+
+
+               ! Model and MDT is set to fbrmdi unless read from file
+               IF ( ldmod ) THEN
+                  logchldata%rmod(iobs,1) = inpfiles(jj)%padd(1,ji,1,1)
+               ELSE
+                  logchldata%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)')'logchl 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, ilogchlidx, zdat )
+
+      !-----------------------------------------------------------------------
+      ! Deallocate input data
+      !-----------------------------------------------------------------------
+      DO jj = 1, inobf
+         CALL dealloc_obfbdata( inpfiles(jj) )
+      END DO
+      DEALLOCATE( inpfiles )
+
+   END SUBROUTINE obs_rea_logchl
+
+END MODULE obs_read_logchl
+
Index: /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_pco2.F90
===================================================================
--- /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_pco2.F90	(revision 7713)
+++ /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_pco2.F90	(revision 7713)
@@ -0,0 +1,465 @@
+MODULE obs_read_pco2
+   !!======================================================================
+   !!                       ***  MODULE obs_read_pco2  ***
+   !! Observation diagnostics: Read the along track pco2 data from
+   !!                          GHRSST or any pco2 data from feedback files
+   !!======================================================================
+
+   !!----------------------------------------------------------------------
+   !!   obs_rea_pco2 : Driver for reading pco2 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_pco2_io            ! I/O for pco2 files
+   USE iom                      ! I/O
+   USE netcdf                   ! NetCDF library
+
+   IMPLICIT NONE
+
+   !! * Routine accessibility
+   PRIVATE
+
+   PUBLIC obs_rea_pco2      ! Read the pco2 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_pco2( kformat, &
+      &                    pco2data, knumfiles, cfilenames, &
+      &                    kvars, kextr, kstp, ddobsini, ddobsend, &
+      &                    ldignmis, ldmod )
+      !!---------------------------------------------------------------------
+      !!
+      !!                   *** ROUTINE obs_rea_pco2 ***
+      !!
+      !! ** Purpose : Read from file the pco2 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) :: &
+         & pco2data     ! pco2 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 pco2data
+      INTEGER, INTENT(IN) :: kextr    ! Number of extra fields for each var in pco2data
+      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_pco2'
+      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 :: ipco2maxtype = 1024
+      INTEGER, DIMENSION(0:ipco2maxtype) :: &
+         & ityp, &
+         & itypmpp
+      INTEGER, DIMENSION(:), ALLOCATABLE :: &
+         & iobsi,    &
+         & iobsj,    &
+         & iproc,    &
+         & iindx,    &
+         & ifileidx, &
+         & ipco2idx
+      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) )
+
+      pco2_files : DO jj = 1, inobf
+          
+         !---------------------------------------------------------------------
+         ! Prints
+         !---------------------------------------------------------------------
+         IF(lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) ' obs_rea_pco2 : 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_pco2( 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 pco2_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), &
+         &      ipco2idx(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
+               ipco2idx(jk) = ji
+               zdat(jk)     = inpfiles(jj)%ptim(ji)
+            ENDIF
+         END DO
+      END DO
+      CALL sort_dp_indx( iobstot, &
+         &               zdat,     &
+         &               iindx   )
+      
+      CALL obs_surf_alloc( pco2data, iobs, & 
+                           kvars, kextr, kstp, jpi, jpj )
+      
+      ! * Read obs/positions, QC, all variable and assign to pco2data
+ 
+      iobs = 0
+
+      ityp   (:) = 0
+      itypmpp(:) = 0
+
+      ioserrcount=0      
+
+      DO jk = 1, iobstot
+         
+         jj = ifileidx(iindx(jk))
+         ji = ipco2idx(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) )
+
+
+               ! pco2 time coordinates
+               pco2data%nyea(iobs) = iyea
+               pco2data%nmon(iobs) = imon
+               pco2data%nday(iobs) = iday
+               pco2data%nhou(iobs) = ihou
+               pco2data%nmin(iobs) = imin
+               
+               ! pco2 space coordinates
+               pco2data%rlam(iobs) = inpfiles(jj)%plam(ji)
+               pco2data%rphi(iobs) = inpfiles(jj)%pphi(ji)
+
+               ! Coordinate search parameters
+               pco2data%mi  (iobs) = inpfiles(jj)%iobsi(ji,1)
+               pco2data%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
+               pco2data%ntyp(iobs) = itype
+               IF ( itype < ipco2maxtype + 1 ) THEN
+                  ityp(itype+1) = ityp(itype+1) + 1
+               ELSE
+                  IF(lwp)WRITE(numout,*)'WARNING:Increase ipco2maxtype in ',&
+                     &                  cpname
+               ENDIF
+
+               ! Bookkeeping data to match observations
+               pco2data%nsidx(iobs) = iobs
+               pco2data%nsfil(iobs) = iindx(jk)
+
+               ! QC flags
+               pco2data%nqc(iobs) = inpfiles(jj)%ivqc(ji,1)
+
+               ! Observed value
+               pco2data%robs(iobs,1) = inpfiles(jj)%pob(1,ji,1)
+
+
+               ! Model and MDT is set to fbrmdi unless read from file
+               IF ( ldmod ) THEN
+                  pco2data%rmod(iobs,1) = inpfiles(jj)%padd(1,ji,1,1)
+               ELSE
+                  pco2data%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)')'pco2 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, ipco2idx, zdat )
+
+      !-----------------------------------------------------------------------
+      ! Deallocate input data
+      !-----------------------------------------------------------------------
+      DO jj = 1, inobf
+         CALL dealloc_obfbdata( inpfiles(jj) )
+      END DO
+      DEALLOCATE( inpfiles )
+
+   END SUBROUTINE obs_rea_pco2
+
+END MODULE obs_read_pco2
+
Index: /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_spm.F90
===================================================================
--- /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_spm.F90	(revision 7713)
+++ /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_spm.F90	(revision 7713)
@@ -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/NEMOGCM/NEMO/OPA_SRC/OBS/obs_spm.F90
===================================================================
--- /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_spm.F90	(revision 7713)
+++ /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_spm.F90	(revision 7713)
@@ -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/NEMOGCM/NEMO/OPA_SRC/OBS/obs_spm_io.F90
===================================================================
--- /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_spm_io.F90	(revision 7713)
+++ /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_spm_io.F90	(revision 7713)
@@ -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/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90
===================================================================
--- /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90	(revision 7712)
+++ /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90	(revision 7713)
@@ -11,4 +11,8 @@
    !!   obs_wri_seaice: Write seaice observation related diagnostics
    !!   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_fco2  : Write fco2 observation related diagnostics
+   !!   obs_wri_pco2  : Write fco2 observation related diagnostics
    !!   obs_wri_stats : Print basic statistics on the data being written out
    !!----------------------------------------------------------------------
@@ -45,4 +49,8 @@
       &   obs_wri_seaice, & ! Write seaice observation related diagnostics
       &   obs_wri_vel, &    ! Write velocity observation related diagnostics
+      &   obs_wri_logchl, & ! Write logchl observation related diagnostics
+      &   obs_wri_spm, &    ! Write spm observation related diagnostics
+      &   obs_wri_fco2, &   ! Write fco2 observation related diagnostics
+      &   obs_wri_pco2, &   ! Write pco2 observation related diagnostics
       &   obswriinfo
    
@@ -930,4 +938,600 @@
    END SUBROUTINE obs_wri_vel
 
+   SUBROUTINE obs_wri_logchl( cprefix, logchldata, padd, pext )
+      !!-----------------------------------------------------------------------
+      !!
+      !!                     *** ROUTINE obs_wri_logchl  ***
+      !!
+      !! ** Purpose : Write logchl 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) :: logchldata   ! Full set of logchl
+      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_logchl'
+      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, logchldata%nsurf, 1, &
+         &                 1 + nadd, next, .TRUE. )
+
+      fbdata%cname(1)      = 'LOGCHL'
+      fbdata%coblong(1)    = 'logchl concentration'
+      fbdata%cobunit(1)    = 'mg/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 LOGCHL'
+      fbdata%caddunit(1,1) = 'mg/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_logchl :'
+         WRITE(numout,*)'~~~~~~~~~~~~~~~~'
+         WRITE(numout,*)'Writing logchl feedback file : ',TRIM(cfname)
+      ENDIF
+
+      ! Transform obs_prof data structure into obfbdata structure
+      fbdata%cdjuldref = '19500101000000'
+      DO jo = 1, logchldata%nsurf
+         fbdata%plam(jo)      = logchldata%rlam(jo)
+         fbdata%pphi(jo)      = logchldata%rphi(jo)
+         WRITE(fbdata%cdtyp(jo),'(I4)') logchldata%ntyp(jo)
+         fbdata%ivqc(jo,:)    = 0
+         fbdata%ivqcf(:,jo,:) = 0
+         IF ( logchldata%nqc(jo) > 10 ) THEN
+            fbdata%ioqc(jo)    = 4
+            fbdata%ioqcf(1,jo) = 0
+            fbdata%ioqcf(2,jo) = logchldata%nqc(jo) - 10
+         ELSE
+            fbdata%ioqc(jo)    = MAX(logchldata%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)    = logchldata%nsfil(jo)
+         IF (ln_grid_global) THEN
+            fbdata%iobsi(jo,1) = logchldata%mi(jo)
+            fbdata%iobsj(jo,1) = logchldata%mj(jo)
+         ELSE
+            fbdata%iobsi(jo,1) = mig(logchldata%mi(jo))
+            fbdata%iobsj(jo,1) = mjg(logchldata%mj(jo))
+         ENDIF
+         CALL greg2jul( 0, &
+            &           logchldata%nmin(jo), &
+            &           logchldata%nhou(jo), &
+            &           logchldata%nday(jo), &
+            &           logchldata%nmon(jo), &
+            &           logchldata%nyea(jo), &
+            &           fbdata%ptim(jo),   &
+            &           krefdate = 19500101 )
+         fbdata%padd(1,jo,1,1) = logchldata%rmod(jo,1)
+         fbdata%pob(1,jo,1)    = logchldata%robs(jo,1)
+         fbdata%pdep(1,jo)     = 0.0
+         fbdata%idqc(1,jo)     = 0
+         fbdata%idqcf(:,1,jo)  = 0
+         IF ( logchldata%nqc(jo) > 10 ) THEN
+            fbdata%ivlqc(1,jo,1) = 4
+            fbdata%ivlqcf(1,1,jo,1) = 0
+            fbdata%ivlqcf(2,1,jo,1) = logchldata%nqc(jo) - 10
+         ELSE
+            fbdata%ivlqc(1,jo,1) = MAX(logchldata%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) = &
+               & logchldata%rext(jo,padd%ipoint(ja))
+         END DO
+         DO je = 1, next
+            fbdata%pext(1,jo,je) = &
+               & logchldata%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_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_fco2( cprefix, fco2data, padd, pext )
+      !!-----------------------------------------------------------------------
+      !!
+      !!                     *** ROUTINE obs_wri_fco2  ***
+      !!
+      !! ** Purpose : Write fco2 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) :: fco2data   ! Full set of fco2
+      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_fco2'
+      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, fco2data%nsurf, 1, &
+         &                 1 + nadd, next, .TRUE. )
+
+      fbdata%cname(1)      = 'fco2'
+      fbdata%coblong(1)    = 'fco2'
+      fbdata%cobunit(1)    = 'uatm'
+      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 fco2'
+      fbdata%caddunit(1,1) = 'uatm'
+      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_fco2 :'
+         WRITE(numout,*)'~~~~~~~~~~~~~~~~'
+         WRITE(numout,*)'Writing fco2 feedback file : ',TRIM(cfname)
+      ENDIF
+
+      ! Transform obs_prof data structure into obfbdata structure
+      fbdata%cdjuldref = '19500101000000'
+      DO jo = 1, fco2data%nsurf
+         fbdata%plam(jo)      = fco2data%rlam(jo)
+         fbdata%pphi(jo)      = fco2data%rphi(jo)
+         WRITE(fbdata%cdtyp(jo),'(I4)') fco2data%ntyp(jo)
+         fbdata%ivqc(jo,:)    = 0
+         fbdata%ivqcf(:,jo,:) = 0
+         IF ( fco2data%nqc(jo) > 10 ) THEN
+            fbdata%ioqc(jo)    = 4
+            fbdata%ioqcf(1,jo) = 0
+            fbdata%ioqcf(2,jo) = fco2data%nqc(jo) - 10
+         ELSE
+            fbdata%ioqc(jo)    = MAX(fco2data%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)    = fco2data%nsfil(jo)
+         IF (ln_grid_global) THEN
+            fbdata%iobsi(jo,1) = fco2data%mi(jo)
+            fbdata%iobsj(jo,1) = fco2data%mj(jo)
+         ELSE
+            fbdata%iobsi(jo,1) = mig(fco2data%mi(jo))
+            fbdata%iobsj(jo,1) = mjg(fco2data%mj(jo))
+         ENDIF
+         CALL greg2jul( 0, &
+            &           fco2data%nmin(jo), &
+            &           fco2data%nhou(jo), &
+            &           fco2data%nday(jo), &
+            &           fco2data%nmon(jo), &
+            &           fco2data%nyea(jo), &
+            &           fbdata%ptim(jo),   &
+            &           krefdate = 19500101 )
+         fbdata%padd(1,jo,1,1) = fco2data%rmod(jo,1)
+         fbdata%pob(1,jo,1)    = fco2data%robs(jo,1)
+         fbdata%pdep(1,jo)     = 0.0
+         fbdata%idqc(1,jo)     = 0
+         fbdata%idqcf(:,1,jo)  = 0
+         IF ( fco2data%nqc(jo) > 10 ) THEN
+            fbdata%ivlqc(1,jo,1) = 4
+            fbdata%ivlqcf(1,1,jo,1) = 0
+            fbdata%ivlqcf(2,1,jo,1) = fco2data%nqc(jo) - 10
+         ELSE
+            fbdata%ivlqc(1,jo,1) = MAX(fco2data%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) = &
+               & fco2data%rext(jo,padd%ipoint(ja))
+         END DO
+         DO je = 1, next
+            fbdata%pext(1,jo,je) = &
+               & fco2data%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_fco2
+
+   SUBROUTINE obs_wri_pco2( cprefix, pco2data, padd, pext )
+      !!-----------------------------------------------------------------------
+      !!
+      !!                     *** ROUTINE obs_wri_pco2  ***
+      !!
+      !! ** Purpose : Write pco2 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) :: pco2data   ! Full set of pco2
+      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_pco2'
+      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, pco2data%nsurf, 1, &
+         &                 1 + nadd, next, .TRUE. )
+
+      fbdata%cname(1)      = 'pco2'
+      fbdata%coblong(1)    = 'pco2'
+      fbdata%cobunit(1)    = 'uatm'
+      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 pco2'
+      fbdata%caddunit(1,1) = 'uatm'
+      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_pco2 :'
+         WRITE(numout,*)'~~~~~~~~~~~~~~~~'
+         WRITE(numout,*)'Writing pco2 feedback file : ',TRIM(cfname)
+      ENDIF
+
+      ! Transform obs_prof data structure into obfbdata structure
+      fbdata%cdjuldref = '19500101000000'
+      DO jo = 1, pco2data%nsurf
+         fbdata%plam(jo)      = pco2data%rlam(jo)
+         fbdata%pphi(jo)      = pco2data%rphi(jo)
+         WRITE(fbdata%cdtyp(jo),'(I4)') pco2data%ntyp(jo)
+         fbdata%ivqc(jo,:)    = 0
+         fbdata%ivqcf(:,jo,:) = 0
+         IF ( pco2data%nqc(jo) > 10 ) THEN
+            fbdata%ioqc(jo)    = 4
+            fbdata%ioqcf(1,jo) = 0
+            fbdata%ioqcf(2,jo) = pco2data%nqc(jo) - 10
+         ELSE
+            fbdata%ioqc(jo)    = MAX(pco2data%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)    = pco2data%nsfil(jo)
+         IF (ln_grid_global) THEN
+            fbdata%iobsi(jo,1) = pco2data%mi(jo)
+            fbdata%iobsj(jo,1) = pco2data%mj(jo)
+         ELSE
+            fbdata%iobsi(jo,1) = mig(pco2data%mi(jo))
+            fbdata%iobsj(jo,1) = mjg(pco2data%mj(jo))
+         ENDIF
+         CALL greg2jul( 0, &
+            &           pco2data%nmin(jo), &
+            &           pco2data%nhou(jo), &
+            &           pco2data%nday(jo), &
+            &           pco2data%nmon(jo), &
+            &           pco2data%nyea(jo), &
+            &           fbdata%ptim(jo),   &
+            &           krefdate = 19500101 )
+         fbdata%padd(1,jo,1,1) = pco2data%rmod(jo,1)
+         fbdata%pob(1,jo,1)    = pco2data%robs(jo,1)
+         fbdata%pdep(1,jo)     = 0.0
+         fbdata%idqc(1,jo)     = 0
+         fbdata%idqcf(:,1,jo)  = 0
+         IF ( pco2data%nqc(jo) > 10 ) THEN
+            fbdata%ivlqc(1,jo,1) = 4
+            fbdata%ivlqcf(1,1,jo,1) = 0
+            fbdata%ivlqcf(2,1,jo,1) = pco2data%nqc(jo) - 10
+         ELSE
+            fbdata%ivlqc(1,jo,1) = MAX(pco2data%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) = &
+               & pco2data%rext(jo,padd%ipoint(ja))
+         END DO
+         DO je = 1, next
+            fbdata%pext(1,jo,je) = &
+               & pco2data%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_pco2
+
    SUBROUTINE obs_wri_stats( fbdata )
       !!-----------------------------------------------------------------------
Index: /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obsfco2_io.h90
===================================================================
--- /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obsfco2_io.h90	(revision 7713)
+++ /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obsfco2_io.h90	(revision 7713)
@@ -0,0 +1,270 @@
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+   SUBROUTINE read_fco2( cdfilename, inpfile, kunit, ldwp, ldgrid )
+      !!---------------------------------------------------------------------
+      !!
+      !!                     ** ROUTINE read_fco2 **
+      !!
+      !! ** Purpose : Read from file the fco2 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_fco2'
+      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 fco2 measurement.            
+      REAL(wp), DIMENSION(:), POINTER :: &
+         & z_phi,   &       ! Latitudes
+         & z_lam            ! Longitudes
+      REAL(wp), DIMENSION(:,:), POINTER :: &
+         & z_fco2         ! fco2 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_fco2        ( 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, 'fco2_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 fco2 data
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, 'fco2', &
+         &                         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_fco2), &
+         &       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_fco2(:,:) /=  zfill)
+         z_fco2(:,:) = (zsca * z_fco2(:,:)) + zoff
+      ELSEWHERE
+         z_fco2(:,:) = 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 fco2 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) = 'fco2'
+
+      ! 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)') 'fco2','  '
+            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_fco2(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_fco2(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_fco2
+
+
Index: /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obslogchl_io.h90
===================================================================
--- /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obslogchl_io.h90	(revision 7713)
+++ /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obslogchl_io.h90	(revision 7713)
@@ -0,0 +1,270 @@
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+   SUBROUTINE read_logchl( cdfilename, inpfile, kunit, ldwp, ldgrid )
+      !!---------------------------------------------------------------------
+      !!
+      !!                     ** ROUTINE read_logchl **
+      !!
+      !! ** Purpose : Read from file the logchl 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_logchl'
+      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 logchl measurement.            
+      REAL(wp), DIMENSION(:), POINTER :: &
+         & z_phi,   &       ! Latitudes
+         & z_lam            ! Longitudes
+      REAL(wp), DIMENSION(:,:), POINTER :: &
+         & z_logchl         ! logchl 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_logchl        ( 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, 'LogChl_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 logchl data
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, 'LogChl', &
+         &                         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_logchl), &
+         &       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_logchl(:,:) /=  zfill)
+         z_logchl(:,:) = (zsca * z_logchl(:,:)) + zoff
+      ELSEWHERE
+         z_logchl(:,:) = 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 logchl 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) = 'LOGCHL'
+
+      ! 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)') 'logchl','  '
+            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_logchl(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_logchl(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_logchl
+
+
Index: /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obspco2_io.h90
===================================================================
--- /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obspco2_io.h90	(revision 7713)
+++ /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obspco2_io.h90	(revision 7713)
@@ -0,0 +1,270 @@
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+   SUBROUTINE read_pco2( cdfilename, inpfile, kunit, ldwp, ldgrid )
+      !!---------------------------------------------------------------------
+      !!
+      !!                     ** ROUTINE read_pco2 **
+      !!
+      !! ** Purpose : Read from file the pco2 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_pco2'
+      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 pco2 measurement.            
+      REAL(wp), DIMENSION(:), POINTER :: &
+         & z_phi,   &       ! Latitudes
+         & z_lam            ! Longitudes
+      REAL(wp), DIMENSION(:,:), POINTER :: &
+         & z_pco2         ! pco2 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_pco2        ( 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, 'pco2_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 pco2 data
+      
+      CALL chkerr( nf90_inq_varid( i_file_id, 'pco2', &
+         &                         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_pco2), &
+         &       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_pco2(:,:) /=  zfill)
+         z_pco2(:,:) = (zsca * z_pco2(:,:)) + zoff
+      ELSEWHERE
+         z_pco2(:,:) = 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 pco2 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) = 'pco2'
+
+      ! 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)') 'pco2','  '
+            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_pco2(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_pco2(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_pco2
+
+
Index: /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obsspm_io.h90
===================================================================
--- /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obsspm_io.h90	(revision 7713)
+++ /branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obsspm_io.h90	(revision 7713)
@@ -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
+
+
