Index: branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
===================================================================
--- branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 6856)
+++ branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 6857)
@@ -33,4 +33,5 @@
    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
@@ -47,4 +48,5 @@
    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
@@ -94,4 +96,6 @@
    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
@@ -183,4 +187,6 @@
       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
@@ -213,4 +219,6 @@
          &            ln_fco2, ln_fco2fb,                             &
          &            fco2files, fco2fbfiles,                         &
+         &            ln_pco2, ln_pco2fb,                             &
+         &            pco2files, pco2fbfiles,                         &
          &            ln_profb_enatim, ln_ignmis, ln_cl4,             &
          &            ln_sstbias, sstbias_files
@@ -240,4 +248,6 @@
       INTEGER :: jnumfco2
       INTEGER :: jnumfco2fb
+      INTEGER :: jnumpco2
+      INTEGER :: jnumpco2fb
       INTEGER :: ji
       INTEGER :: jset
@@ -255,4 +265,6 @@
       ln_fco2     = .FALSE.
       ln_fco2fb   = .FALSE.
+      ln_pco2     = .FALSE.
+      ln_pco2fb   = .FALSE.
       
       !Initalise all values in namelist arrays
@@ -281,4 +293,6 @@
       fco2files(:) = ''
       fco2fbfiles(:) = ''
+      pco2files(:) = ''
+      pco2fbfiles(:) = ''
       sstbias_files(:) = ''
       endailyavtypes(:) = -1
@@ -409,4 +423,14 @@
          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
       
@@ -445,4 +469,6 @@
          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
@@ -570,5 +596,5 @@
          IF (ln_fco2) THEN
             DO ji = 1, jnumfco2
-               WRITE(numout,'(1X,2A)') '             fco2 input observation file name  fco2files = ', &
+               WRITE(numout,'(1X,2A)') '             fco2 input observation file name            fco2files = ', &
                   TRIM(fco2files(ji))
             END DO
@@ -576,6 +602,18 @@
          IF (ln_fco2fb) THEN
             DO ji = 1, jnumfco2fb
-               WRITE(numout,'(1X,2A)') '             Feedback fco2 input observation file name  fco2fbfiles = ', &
+               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
@@ -615,9 +653,9 @@
          & ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. &
          & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ).AND.( .NOT. ln_logchl ).AND. &
-         & ( .NOT. ln_spm ).AND.( .NOT. ln_fco2 ) ) THEN
+         & ( .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,', &
-            &                    ' ln_logchl, ln_spm, ln_fco2 are all set to .FALSE.'
+            &                    ' ln_logchl, ln_spm, ln_fco2, ln_pco2 are all set to .FALSE.'
          nwarn = nwarn + 1
       ENDIF
@@ -1282,4 +1320,59 @@
  
       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
@@ -1302,4 +1395,5 @@
       !!               - Sea surface spm
       !!               - Sea surface fco2
+      !!               - Sea surface pco2
       !!
       !! ** Action  : 
@@ -1340,15 +1434,17 @@
 #endif
 #if defined key_hadocc
-      USE trc, ONLY :  &                ! HadOCC chlorophyll and fCO2
+      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 and fCO2
+      USE trc, ONLY :  &                ! MEDUSA chlorophyll, fCO2 and pCO2
          & MEDUSA_CHL, &
          & MEDUSA_FCO2, &
+         & MEDUSA_PCO2, &
          & MEDUSA_FILL_FLT
 #elif defined key_fabm
-      !USE ???                           ! ERSEM chlorophyll and fCO2
+      !USE ???                           ! ERSEM chlorophyll, fCO2 and pCO2
 #endif
 #if defined key_spm
@@ -1370,4 +1466,5 @@
       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
@@ -1385,4 +1482,8 @@
       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
       CHARACTER(LEN=20) :: datestr=" ",timestr=" "
@@ -1582,4 +1683,40 @@
       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(:,:)  =  ???                 ! pCO2 from ERSEM
+         CALL ctl_stop( ' Trying to run pCO2 observation operator', &
+            &           ' but not properly implemented for FABM-ERSEM yet' )
+#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) 
@@ -1617,4 +1754,5 @@
       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
@@ -1939,4 +2077,34 @@
             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
Index: branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
===================================================================
--- branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 6856)
+++ branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 6857)
@@ -28,4 +28,6 @@
    !!                    observations
    !!   obs_fco2_opt :   Compute the model counterpart of fco2
+   !!                    observations
+   !!   obs_pco2_opt :   Compute the model counterpart of pco2
    !!                    observations
    !!----------------------------------------------------------------------
@@ -72,5 +74,6 @@
       &   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_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
@@ -2541,4 +2544,164 @@
    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_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_pco2.F90
===================================================================
--- branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_pco2.F90	(revision 6857)
+++ branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_pco2.F90	(revision 6857)
@@ -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_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_pco2_io.F90
===================================================================
--- branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_pco2_io.F90	(revision 6857)
+++ branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_pco2_io.F90	(revision 6857)
@@ -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_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
===================================================================
--- branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 6856)
+++ branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 6857)
@@ -15,4 +15,5 @@
    !!   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
@@ -47,4 +48,5 @@
       & 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  
 
@@ -1741,4 +1743,187 @@
    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
+                              ! Global counters for observations that
+      INTEGER :: iotdobsmpp   !  - outside time domain
+      INTEGER :: iosdsobsmpp  !  - outside space domain
+      INTEGER :: ilansobsmpp  !  - within a model land cell
+      INTEGER :: inlasobsmpp  !  - close to land
+      INTEGER :: igrdobsmpp   !  - fail the grid search
+      LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
+         & llvalid            ! data selection
+      INTEGER :: jobs         ! Obs. loop variable
+      INTEGER :: jstp         ! Time loop variable
+      INTEGER :: inrc         ! Time index variable
+
+      IF (lwp) WRITE(numout,*)'obs_pre_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
+
+      ! -----------------------------------------------------------------------
+      ! 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           ) 
+         
+      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
+      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
+      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
+
+      ! -----------------------------------------------------------------------
+      ! Copy useful data from the 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,*) ' 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_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_pco2.F90
===================================================================
--- branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_pco2.F90	(revision 6857)
+++ branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_pco2.F90	(revision 6857)
@@ -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_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90
===================================================================
--- branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90	(revision 6856)
+++ branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90	(revision 6857)
@@ -14,4 +14,5 @@
    !!   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
    !!----------------------------------------------------------------------
@@ -51,4 +52,5 @@
       &   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
    
@@ -1383,4 +1385,153 @@
    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_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obspco2_io.h90
===================================================================
--- branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obspco2_io.h90	(revision 6857)
+++ branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obspco2_io.h90	(revision 6857)
@@ -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
+
+
