Index: branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8455)
+++ branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8456)
@@ -1295,4 +1295,6 @@
     ln_logchlinc = .false. !  Logical switch for applying logchl increments
     ln_logchlbal = .false. !  Logical switch for applying logchl multivariate balancing
+    ln_fco2inc = .false.   !  Logical switch for applying fCO2 increments
+    ln_pco2inc = .false.   !  Logical switch for applying pCO2 increments
     ln_temnofreeze = .false. !  Logical to not add increments if temperature would fall below freezing
     nitbkg    = 0          !  Timestep of background in [0,nitend-nit000-1]
Index: branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbal.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbal.F90	(revision 8455)
+++ branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbal.F90	(revision 8456)
@@ -74,5 +74,5 @@
 
          IF ( ln_logchlinc ) THEN
-#if defined key_medusa && defined key_foam_medusa
+#if defined key_medusa
             CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chn', logchl_balinc(:,:,:,jpchn) )
             CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chd', logchl_balinc(:,:,:,jpchd) )
@@ -104,4 +104,22 @@
          ENDIF
 
+         IF ( ln_pco2inc ) THEN
+#if defined key_medusa
+            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_dic', pco2_balinc(:,:,:,jpdic) )
+            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_alk', pco2_balinc(:,:,:,jpalk) )
+#elif defined key_hadocc
+            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_dic', pco2_balinc(:,:,:,jp_had_dic) )
+            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_alk', pco2_balinc(:,:,:,jp_had_alk) )
+#endif
+         ELSE IF ( ln_fco2inc ) THEN
+#if defined key_medusa
+            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', pco2_balinc(:,:,:,jpdic) )
+            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', pco2_balinc(:,:,:,jpalk) )
+#elif defined key_hadocc
+            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', pco2_balinc(:,:,:,jp_had_dic) )
+            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', pco2_balinc(:,:,:,jp_had_alk) )
+#endif
+         ENDIF
+
          CALL iom_close( inum )
       ELSE
Index: branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 8455)
+++ branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 8456)
@@ -23,4 +23,5 @@
    !!   seaice_asm_inc : Apply the seaice increment
    !!   logchl_asm_inc : Apply the logchl increment
+   !!   pco2_asm_inc   : Apply the pCO2/fCO2 increment
    !!----------------------------------------------------------------------
    USE wrk_nemo         ! Memory Allocation
@@ -64,8 +65,13 @@
    USE asmlogchlbal_medusa, ONLY: &
       & asm_logchl_bal_medusa
+   USE asmpco2bal, ONLY: &
+      & asm_pco2_bal
    USE par_medusa
+   USE mocsy_mainmod
 #elif defined key_hadocc
    USE asmlogchlbal_hadocc, ONLY: &
       & asm_logchl_bal_hadocc
+   USE asmpco2bal, ONLY: &
+      & asm_pco2_bal
    USE par_hadocc
 #endif
@@ -81,4 +87,5 @@
    PUBLIC   seaice_asm_inc !: Apply the seaice increment
    PUBLIC   logchl_asm_inc !: Apply the logchl increment
+   PUBLIC   pco2_asm_inc   !: Apply the pCO2/fCO2 increment
 
 #if defined key_asminc
@@ -97,4 +104,6 @@
    LOGICAL, PUBLIC :: ln_logchlinc = .FALSE.   !: No log10(chlorophyll) increment
    LOGICAL, PUBLIC :: ln_logchlbal = .FALSE.   !: Don't apply multivariate log10(chlorophyll) balancing
+   LOGICAL, PUBLIC :: ln_pco2inc = .FALSE.     !: No pCO2 increment
+   LOGICAL, PUBLIC :: ln_fco2inc = .FALSE.     !: No fCO2 increment
    LOGICAL, PUBLIC :: ln_salfix = .FALSE.      !: Apply minimum salinity check
    LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
@@ -123,6 +132,8 @@
 
    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: logchl_bkginc !: Increment to background logchl
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: pco2_bkginc   !: Increment to background pCO2
 #if defined key_top
    REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: logchl_balinc  !: Increment to BGC variables from logchl assim
+   REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: pco2_balinc    !: Increment to BGC variables from pCO2 assim
 #endif
 #if defined key_hadocc || (defined key_medusa && defined key_foam_medusa)
@@ -192,4 +203,5 @@
          &                 ln_trainc, ln_dyninc, ln_sshinc,                &
          &                 ln_logchlinc, ln_logchlbal,                     &
+         &                 ln_pco2inc, ln_fco2inc,                         &
          &                 ln_asmdin, ln_asmiau,                           &
          &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
@@ -229,4 +241,6 @@
          WRITE(numout,*) '      Logical switch for applying logchl increments         ln_logchlinc = ', ln_logchlinc
          WRITE(numout,*) '      Logical switch for logchl multivariate balancing      ln_logchlbal = ', ln_logchlbal
+         WRITE(numout,*) '      Logical switch for applying pCO2 increments             ln_pco2inc = ', ln_pco2inc
+         WRITE(numout,*) '      Logical switch for applying fCO2 increments             ln_fco2inc = ', ln_fco2inc
          WRITE(numout,*) '      Logical switch for Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau
          WRITE(numout,*) '      Timestep of background in [0,nitend-nit000-1]            nitbkg    = ', nitbkg
@@ -287,7 +301,7 @@
       IF (      ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) &
          & .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ).OR. &
-         &        ( ln_logchlinc ) )) &
+         &        ( ln_logchlinc ).OR.( ln_pco2inc ).OR.( ln_fco2inc ) )) &
          & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', &
-         &                ' and ln_logchlinc is set to .true.', &
+         &                ' ln_logchlinc, ln_pco2inc, and ln_fco2inc is set to .true.', &
          &                ' but ln_asmdin and ln_asmiau are both set to .false. :', &
          &                ' Inconsistent options')
@@ -298,7 +312,7 @@
 
       IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) &
-         & .AND.( .NOT. ln_logchlinc ) )  &
+         & .AND.( .NOT. ln_logchlinc ).AND.( .NOT. ln_pco2inc ).AND.( .NOT. ln_fco2inc ) )  &
          & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', &
-         &                ' and ln_logchlinc are set to .false. :', &
+         &                ' ln_logchlinc, ln_pco2inc, and ln_fco2inc are set to .false. :', &
          &                ' The assimilation increments are not applied')
 
@@ -321,14 +335,18 @@
          &                ' the cycle interval')
 
-      IF ( ( ln_balwri ).AND.( .NOT. ( ln_logchlinc ) ) ) THEN
-         CALL ctl_warn( ' Balancing increments are only calculated for logchl', &
-            &           ' Not assimilating logchl, so ln_balwri will be set to .false.')
+      IF ( ( ln_balwri ).AND.( .NOT. ln_logchlinc ).AND.( .NOT. ln_pco2inc ).AND.( .NOT. ln_fco2inc ) ) THEN
+         CALL ctl_warn( ' Balancing increments are only calculated for logchl and pCO2/fCO2', &
+            &           ' Not assimilating logchl, pCO2 or fCO2, so ln_balwri will be set to .false.')
          ln_balwri = .FALSE.
       ENDIF
 
-      IF ( ( ln_logchlbal ).AND.( .NOT. ( ln_logchlinc ) ) ) THEN
+      IF ( ( ln_logchlbal ).AND.( .NOT. ln_logchlinc ) ) THEN
          CALL ctl_warn( ' Cannot calculate logchl balancing increments if not assimilating logchl', &
             &           ' Not assimilating logchl, so ln_logchlbal will be set to .false.')
          ln_logchlbal = .FALSE.
+      ENDIF
+
+      IF ( ( ln_pco2inc ).AND.( ln_fco2inc ) ) THEN
+         CALL ctl_stop( ' Can only assimilate pCO2 OR fCO2, not both' )
       ENDIF
 
@@ -436,6 +454,14 @@
 #endif
       ENDIF
+      IF ( ln_pco2inc .OR. ln_fco2inc ) THEN
+         ALLOCATE( pco2_bkginc(jpi,jpj) )
+         pco2_bkginc(:,:) = 0.0
+#if defined key_top
+         ALLOCATE( pco2_balinc(jpi,jpj,jpk,jptra) )
+         pco2_balinc(:,:,:,:) = 0.0
+#endif
+      ENDIF
       IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) &
-         &  .OR.( ln_logchlinc ) ) THEN
+         &  .OR.( ln_logchlinc ).OR.( ln_pco2inc ).OR.( ln_fco2inc ) ) THEN
 
          !--------------------------------------------------------------------
@@ -518,4 +544,17 @@
             ! to allow for differences in masks
             WHERE( ABS( logchl_bkginc(:,:) ) > 1.0e+10 ) logchl_bkginc(:,:) = 0.0
+         ENDIF
+
+         IF ( ln_pco2inc .OR. ln_fco2inc ) THEN
+            IF ( ln_pco2inc ) THEN
+               CALL iom_get( inum, jpdom_autoglo, 'bckinpco2', pco2_bkginc(:,:), 1 )
+            ELSE IF ( ln_fco2inc ) THEN
+               CALL iom_get( inum, jpdom_autoglo, 'bckinfco2', pco2_bkginc(:,:), 1 )
+            ENDIF
+            ! Apply the masks
+            pco2_bkginc(:,:) = pco2_bkginc(:,:) * tmask(:,:,1)
+            ! Set missing increments to 0.0 rather than 1e+20
+            ! to allow for differences in masks
+            WHERE( ABS( pco2_bkginc(:,:) ) > 1.0e+10 ) pco2_bkginc(:,:) = 0.0
          ENDIF
 
@@ -703,5 +742,34 @@
             CALL iom_get( inum, jpdom_autoglo, 'medusa_oxy', tracer_bkg(:,:,:,jpoxy) )
 #endif
-         ENDIF
+         ELSE IF ( ln_pco2inc .OR. ln_fco2inc ) THEN
+#if defined key_hadocc
+            CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
+            CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
+#elif defined key_medusa
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
+#endif
+         ENDIF
+
+         CALL iom_close( inum )
+         
+         DO jt = 1, jptra
+            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
+         END DO
+      
+      ELSE IF ( ln_pco2inc .OR. ln_fco2inc ) THEN
+
+         ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
+         tracer_bkg(:,:,:,:) = 0.0
+
+         CALL iom_open( c_asmbkg, inum )
+         
+#if defined key_hadocc
+         CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
+         CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
+#elif defined key_medusa
+         CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
+         CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
+#endif
 
          CALL iom_close( inum )
@@ -1454,4 +1522,247 @@
       !
    END SUBROUTINE logchl_asm_inc
+
+
+   SUBROUTINE pco2_asm_inc( kt )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE pco2_asm_inc  ***
+      !!          
+      !! ** Purpose : Apply the pco2/fco2 assimilation increments.
+      !!
+      !! ** Method  : Calculate increments to state variables using carbon
+      !!              balancing.
+      !!              Direct initialization or Incremental Analysis Updating.
+      !!
+      !! ** Action  : 
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(IN)                   :: kt           ! Current time step
+      !
+      INTEGER                               :: ji, jj, jk   ! Loop counters
+      INTEGER                               :: it           ! Index
+      INTEGER                               :: jkmax        ! Index of mixed layer depth
+      REAL(wp)                              :: zincwgt      ! IAU weight for current time step
+      REAL(wp), DIMENSION(jpi,jpj)          :: zmld         ! Mixed layer depth for balancing
+      REAL(wp), DIMENSION(jpi,jpj)          :: pco2_bkginc_temp ! For fCO2 conversion
+      REAL(wp), DIMENSION(jpi,jpj)          :: dic_bkg_temp ! DIC background
+      REAL(wp), DIMENSION(jpi,jpj)          :: alk_bkg_temp ! Alkalinity background
+      REAL(wp), DIMENSION(jpi,jpj)          :: tem_bkg_temp ! Temperature background
+      REAL(wp), DIMENSION(jpi,jpj)          :: sal_bkg_temp ! Salinity background
+      REAL(wp), DIMENSION(1)                :: patm         ! Atmospheric pressure
+      REAL(wp), DIMENSION(1)                :: phyd         ! Hydrostatic pressure
+
+      ! Coefficients for fCO2 to pCO2 conversion
+      REAL(wp)                              :: zcoef_fco2_1 = -1636.75
+      REAL(wp)                              :: zcoef_fco2_2 = 12.0408
+      REAL(wp)                              :: zcoef_fco2_3 = 0.0327957
+      REAL(wp)                              :: zcoef_fco2_4 = 0.0000316528
+      REAL(wp)                              :: zcoef_fco2_5 = 2.0
+      REAL(wp)                              :: zcoef_fco2_6 = 57.7
+      REAL(wp)                              :: zcoef_fco2_7 = 0.118
+      REAL(wp)                              :: zcoef_fco2_8 = 82.0578
+      !!----------------------------------------------------------------------
+
+      IF ( kt <= nit000 ) THEN
+
+         !--------------------------------------------------------------------
+         ! DIC and alkalinity balancing
+         !--------------------------------------------------------------------
+
+         IF ( ln_fco2inc ) THEN
+#if defined key_medusa && defined key_roam
+            ! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine
+            patm(1) = 1.0
+            phyd(1) = 0.0
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  CALL f2pCO2( pco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) )
+               END DO
+            END DO
+#else
+            ! If assimilating fCO2, then convert to pCO2 using temperature
+            ! See flux_gas.F90 within HadOCC for details of calculation
+            pco2_bkginc_temp(:,:) = pco2_bkginc(:,:) /                                                             &
+               &                    EXP((zcoef_fco2_1                                                            + &
+               &                         zcoef_fco2_2 * (tsn(:,:,1,1)+rt0)                                       - &
+               &                         zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)                    + &
+               &                         zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + &
+               &                         zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0)))      / &
+               &                        (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0)))
+#endif
+         ELSE
+            pco2_bkginc_temp(:,:) = pco2_bkginc(:,:)
+         ENDIF
+
+         ! Account for physics assimilation if required
+         IF ( ln_trainc ) THEN
+            tem_bkg_temp(:,:) = tsn(:,:,1,1) + t_bkginc(:,:,1)
+            sal_bkg_temp(:,:) = tsn(:,:,1,2) + s_bkginc(:,:,1)
+         ELSE
+            tem_bkg_temp(:,:) = tsn(:,:,1,1)
+            sal_bkg_temp(:,:) = tsn(:,:,1,2)
+         ENDIF
+
+#if defined key_medusa
+         ! Account for logchl balancing if required
+         IF ( ln_logchlinc .AND. ln_logchlbal ) THEN
+            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + logchl_balinc(:,:,1,jpdic)
+            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + logchl_balinc(:,:,1,jpalk)
+         ELSE
+            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic)
+            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk)
+         ENDIF
+
+         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
+            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        &
+            &               pco2_balinc(:,:,1,jpdic), pco2_balinc(:,:,1,jpalk) )
+
+#elif defined key_hadocc
+         ! Account for logchl balancing if required
+         IF ( ln_logchlinc .AND. ln_logchlbal ) THEN
+            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + logchl_balinc(:,:,1,jp_had_dic)
+            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + logchl_balinc(:,:,1,jp_had_alk)
+         ELSE
+            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic)
+            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk)
+         ENDIF
+
+         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
+            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        &
+            &               pco2_balinc(:,:,1,jp_had_dic), pco2_balinc(:,:,1,jp_had_alk) )
+
+#else
+         CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', &
+            &           'but not defined a biogeochemical model' )
+#endif
+
+         ! Select mixed layer
+         SELECT CASE( mld_choice_bgc )
+         CASE ( 1 )                   ! Turbocline/mixing depth [W points]
+            zmld(:,:) = hmld(:,:)
+         CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
+            zmld(:,:) = hmlp(:,:)
+         CASE ( 3 )                   ! Kara MLD [Interpolated]
+#if defined key_karaml
+            IF ( ln_kara ) THEN
+               zmld(:,:) = hmld_kara(:,:)
+            ELSE
+               CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
+                  &           ' but ln_kara=.false.' )
+            ENDIF
+#else
+            CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
+               &           ' but is not defined' )
+#endif
+         CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points]
+            !zmld(:,:) = hmld_tref(:,:)
+            CALL ctl_stop( ' hmld_tref mixed layer requested for pCO2 assimilation,', &
+               &           ' but is not available in this version' )
+         CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
+            zmld(:,:) = hmlpt(:,:)
+         END SELECT
+
+         ! Propagate through mixed layer
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !
+               jkmax = jpk-1
+               DO jk = jpk-1, 1, -1
+                  IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
+                     & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
+                     zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
+                     jkmax = jk
+                  ENDIF
+               END DO
+               !
+               DO jk = 2, jkmax
+                  pco2_balinc(ji,jj,jk,:) = pco2_balinc(ji,jj,1,:)
+               END DO
+               !
+            END DO
+         END DO
+
+      ENDIF
+
+      IF ( ln_asmiau ) THEN
+
+         !--------------------------------------------------------------------
+         ! Incremental Analysis Updating
+         !--------------------------------------------------------------------
+
+         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
+
+            it = kt - nit000 + 1
+            zincwgt = wgtiau(it)      ! IAU weight for the current time step
+            ! note this is not a tendency so should not be divided by rdt
+
+            IF(lwp) THEN
+               WRITE(numout,*) 
+               IF ( ln_pco2inc ) THEN
+                  WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', &
+                     &  kt,' with IAU weight = ', wgtiau(it)
+               ELSE IF ( ln_fco2inc ) THEN
+                  WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', &
+                     &  kt,' with IAU weight = ', wgtiau(it)
+               ENDIF
+               WRITE(numout,*) '~~~~~~~~~~~~'
+            ENDIF
+
+            ! Update the biogeochemical variables
+            ! Add directly to trn and trb, rather than to tra, as not a tendency
+#if defined key_medusa && defined key_foam_medusa
+            DO jk = 1, jpkm1
+               trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + &
+                  &                          pco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
+               trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + &
+                  &                          pco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
+            END DO
+#elif defined key_hadocc
+            DO jk = 1, jpkm1
+               trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + &
+                  &                          pco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
+               trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + &
+                  &                          pco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
+            END DO
+#endif
+
+            ! Do not deallocate arrays - needed by asm_bal_wri
+            ! which is called at end of model run
+
+         ENDIF
+
+      ELSEIF ( ln_asmdin ) THEN 
+
+         !--------------------------------------------------------------------
+         ! Direct Initialization
+         !--------------------------------------------------------------------
+         
+         IF ( kt == nitdin_r ) THEN
+
+            neuler = 0                    ! Force Euler forward step
+
+#if defined key_medusa && defined key_foam_medusa
+            ! Initialize the now fields with the background + increment
+            ! Background currently is what the model is initialised with
+            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pCO2 assimilation', &
+               &           ' Background state is taken from model rather than background file' )
+            trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
+               &                         pco2_balinc(:,:,:,jp_msa0:jp_msa1)
+            trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
+#elif defined key_hadocc
+            ! Initialize the now fields with the background + increment
+            ! Background currently is what the model is initialised with
+            CALL ctl_warn( ' Doing direct initialisation of HadOCC with pCO2 assimilation', &
+               &           ' Background state is taken from model rather than background file' )
+            trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
+               &                         pco2_balinc(:,:,:,jp_had0:jp_had1)
+            trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
+#endif
+ 
+            ! Do not deallocate arrays - needed by asm_bal_wri
+            ! which is called at end of model run
+         ENDIF
+         !
+      ENDIF
+      !
+   END SUBROUTINE pco2_asm_inc
    
    !!======================================================================
Index: branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmpco2bal.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmpco2bal.F90	(revision 8456)
+++ branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmpco2bal.F90	(revision 8456)
@@ -0,0 +1,366 @@
+MODULE asmpco2bal
+   !Contains subroutine to convert a global pco2 increment to independent DIC and alkalinity increments
+
+   USE par_kind, ONLY: wp
+   USE par_oce, ONLY: jpi, jpj, jpk
+   USE dom_oce, ONLY: tmask
+#if defined key_hadocc
+   USE had_bgc_const, ONLY: RHO_WATER_SI
+#else
+   USE phycst, ONLY: RHO_WATER_SI => rau0
+#endif
+   
+   IMPLICIT NONE
+   PRIVATE
+   
+   PUBLIC asm_pco2_bal
+   CONTAINS 
+
+   SUBROUTINE  asm_pco2_bal(pco2inc, tco2_tracer_bkg, alk_tracer_bkg, &
+      &                     tem_tracer_bkg, sal_tracer_bkg,           &
+      &                     tco2_tracer_inc, alk_tracer_inc)
+      !Converts a global pco2 increment to independent DIC and alkalinity increments
+      !Calculates the solution to the equation [DDIC, DAlk]=(DpCO2 / |Grad (pCO2)|)* Grad(pCO2); where Grad(pco2)=[dpCO2/dDIC, dpCO2/dalk]
+      !Gradient code pco2_grad is a modification of assimcd1.mf77 by Ian Totterdell.
+
+      !James While 06/09
+
+      IMPLICIT NONE
+
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) ::  pco2inc            !pCO2 increment
+     
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: tco2_tracer_bkg,      &    !In: tCO2 (DIC) background
+                                                  alk_tracer_bkg      !In alkalinity background
+     
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: tem_tracer_bkg,      &    !In: temperature background
+                                                  sal_tracer_bkg      !In salinity background
+     
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(OUT) :: tco2_tracer_inc,      &    !Out: tCO2 (DIC) increment
+                                                   alk_tracer_inc      !Out alkalinity
+						   
+      !local variables
+      REAL(wp), DIMENSION(jpi,jpj) ::  Dpco2_tco2, Dpco2_alk
+      REAL(wp) :: wgt_fct
+      INTEGER :: ji,jj
+      
+      
+      !calculate the gradients with respect to DIC and Alkalinity
+      CALL pco2_grad(tco2_tracer_bkg, alk_tracer_bkg, tem_tracer_bkg, sal_tracer_bkg, Dpco2_tco2, Dpco2_alk)
+
+      
+      !loop over the ocean grid calculating the increments to tCO2 (DIC) and alkalinity
+      DO jj=1, jpj
+        DO ji=1, jpi
+
+          wgt_fct=pco2inc(ji,jj)/( (Dpco2_tco2(ji,jj)*Dpco2_tco2(ji,jj)) + (Dpco2_alk(ji,jj)*Dpco2_alk(ji,jj)) )
+          tco2_tracer_inc(ji,jj)=wgt_fct * Dpco2_tco2(ji,jj) * tmask(ji,jj,1)
+          alk_tracer_inc(ji,jj)=wgt_fct * Dpco2_alk(ji,jj) * tmask(ji,jj,1)
+ 
+        ENDDO
+      ENDDO
+       						  
+                                    
+
+   END SUBROUTINE  asm_pco2_bal
+   
+  
+   SUBROUTINE pco2_grad( tco2_tracer_bkg, alk_tracer_bkg, tem_tracer_bkg, sal_tracer_bkg, DPDC, DPDA )
+!
+!-------M---------M---------M---------M---------M---------M---------M--
+!
+!  Ian J Totterdell
+!  2007 August 08
+!
+!  Code for calculating derivatives of pCO2 wrt DIC and Alk
+!  
+!  J While 
+!  June 2009
+!  Modiefied for use in NEMO
+
+!-------M---------M---------M---------M---------M---------M---------M--
+      IMPLICIT NONE
+!
+!
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: tco2_tracer_bkg,      &    !In: tCO2 (DIC) background
+                                                  alk_tracer_bkg      !In alkalinity background
+     
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: tem_tracer_bkg,      &    !In: temperature background
+                                                  sal_tracer_bkg      !In salinity background
+!
+      REAL(wp), INTENT(OUT), DIMENSION(jpi,jpj) :: &
+        DPDC,   &
+        DPDA
+!
+!
+!  Local variables.
+!
+      INTEGER :: &
+        mxiter,  &
+        ji,jj,   &
+        iter
+!
+      REAL(wp) :: &
+        a1_co2, a2_co2, a3_co2,               &
+        b1_co2, b2_co2, b3_co2,               &
+        a1_k1, a2_k1, a3_k1,                  &
+        b1_k1, b2_k1, b3_k1, b4_k1, b5_k1,    &
+        a1_k2, a2_k2, a3_k2,                  &
+        b1_k2, b2_k2, b3_k2, b4_k2, b5_k2,    &
+        a1_kb, a2_kb, a3_kb, a4_kb, a5_kb,    &
+        b1_kb, b2_kb, b3_kb,                  &
+        c1_kb, c2_kb, c3_kb, d_kb,            &
+        a1_kw, a2_kw, a3_kw,                  &
+        b1_kw, b2_kw, b3_kw, c_kw
+!
+      REAL(wp) :: &
+        tk, logtk, ttk,             &
+        sal, sqsal,                 &
+        kone, ktwo, kboric, kwater, &
+        alkprime, qa,qb,qc,         &
+        fxa, xtrma, xtrmb,          &
+        dxda, dxdc,                 &
+        dgda, dgdc
+!
+      REAL(wp), DIMENSION(jpi,jpj) :: &
+      co2_sol,     &
+      rek,         &
+      rg1,         &
+      rg2,         &
+      rg3,         &
+      boron,       &
+      tco2,        &
+      alk_t,       &
+      xs,          &
+      x1,          &
+      x2,          &
+      x2p,         &
+      x3,          &
+      xx,          &
+      xf
+!
+!  Set parameter and constant values.
+!
+      PARAMETER ( mxiter = 10 )
+!
+!  CO2 solubility, constants for the calculation of.
+!   From Weiss (1974)
+!
+      DATA &
+      a1_co2 / -60.2409 /,  a2_co2 / 93.4517 /,   a3_co2 / 23.3585 /, &
+      b1_co2 / 0.023517 /,  b2_co2 / -0.023656 /, b3_co2 / 0.0047036 /
+!
+!  kone, constants for the calculation of.
+!   From Roy et al (1993)
+!
+      DATA a1_k1 /-2307.1266/,  a2_k1 / 2.83655 /,   a3_k1 /-1.5529413/ &
+        , b1_k1 / -4.0484 /,   b2_k1 /-0.20760841/, b3_k1 /0.08468345/ &
+        , b4_k1 /-0.00654208/, b5_k1 / -0.001005 /
+!
+!  ktwo, constants for the calculation of.
+!   From Roy et al (1993)
+!
+      DATA a1_k2 /-3351.6106/,  a2_k2 / -9.226508 /, a3_k2 /-0.2005743/ &
+        , b1_k2 / -23.9722 /,  b2_k2 /-0.106901773/,b3_k2 /0.1130822/ &
+        , b4_k2 /-0.00846934/, b5_k2 / -0.001005 /
+!
+!  kboric, constants for the calculation of.
+!   Using data from Dickson (1990);
+!
+      DATA a1_kb / -8966.90 /,  a2_kb / -2890.53 /,  a3_kb / -77.942 / &
+        , a4_kb / 1.728 /,     a5_kb /-0.0996/ &
+        , b1_kb / 148.0248 /,  b2_kb / 137.1942 /,  b3_kb / 1.62142 / &
+        , c1_kb / -24.4344 /,  c2_kb / -25.085 /,   c3_kb / -0.2474 / &
+        , d_kb / 0.053105 /
+!
+!  kwater, constants for the calculation of.
+!   From Millero (1995, p670), using composite data.
+!
+      DATA a1_kw / -13847.26 /, a2_kw / 148.9652 /,  a3_kw / -23.6521 / &
+        , b1_kw / 118.67 /,    b2_kw / -5.977 /,    b3_kw / 1.0495 / &
+        , c_kw / -0.01615 / 
+!
+!  Executable code begins here.
+!
+!  Calculate the carbon system constants.
+!
+
+      DO ji=1,jpi
+        DO jj=1,jpj
+!
+          tk = 273.15 + MAX( -2.0, MIN( 40.0, tem_tracer_bkg(ji,jj) ) )
+          logtk = LOG(tk)
+          ttk = MAX( 2.71, 0.01 * ( 273.15 + tem_tracer_bkg(ji,jj) ) )
+          sal = MAX( 0.0, MIN( 45.0, sal_tracer_bkg(ji,jj) ) )
+          sqsal = SQRT(sal)
+!
+!  Calculate the CO2 solubility from Weiss (1974)
+!   Units are mol/(kg.soln)/atm .
+!   Check value (Zeebe and Wolf-G; tc=25,sal=35): ln(co2_sol) = -3.5617
+!
+          co2_sol(ji,jj) = EXP( &
+           a1_co2 + a2_co2 / ttk + a3_co2 * LOG(ttk) + &
+           sal * ( b1_co2 + ttk * ( b2_co2 + b3_co2 * ttk ) ) )
+!
+!  Calculate the 1st constant of carbonic acid from Roy et al. (1993)
+!   Units are mol/(kg.soln) on the total hydrogen ion scale.
+!   K1 = (HCO3-).(H+)/(CO2(g))
+!   Check value (Dickson and Goyet; tc=25,sal=35): ln(K1) = -13.4847
+!
+          kone &
+           = 1.0e-7 *( 1.0 + b5_k1*sal ) * EXP( a1_k1/tk + a2_k1 + a3_k1*logtk &
+           + ( b1_k1/tk + b2_k1 )*sqsal + b3_k1*sal + b4_k1*sal*sqsal )
+!
+!  Calculate the 2nd constant of carbonic acid from Roy et al. (1993)
+!   Units are mol/(kg.soln) on the total hydrogen ion scale.
+!   K2 = (CO3--).(H+)/(HCO3-)
+!   Check value (Dickson and Goyet; tc=25,sal=35): ln(K2) = -20.5504
+!
+          ktwo &
+           = 1.0e-7 *( 1.0 + b5_k2*sal ) * EXP( a1_k2/tk + a2_k2 + a3_k2*logtk &
+           + ( b1_k2/tk + b2_k2 )*sqsal + b3_k2*sal + b4_k2*sal*sqsal )
+!
+!  Calculate the equilibrium constant of boric acid from Dickson (1990).
+!   Units are mol/(kg.soln) on the total hydrogen ion scale.
+!   Kb = (B(OH)4-).(H+)/(B(OH)3) (with (Bt) = (B(OH)3) + (B(OH)4-) ).
+!   Check value (Dickson and Goyet; tc=25,sal=35): ln(Kb) = -19.7964
+!
+          kboric &
+           = 1.0e-7 * EXP( ( a1_kb + a2_kb*sqsal + a3_kb*sal + a4_kb*sal*sqsal &
+           + a5_kb*sal*sal )/tk + ( b1_kb + b2_kb*sqsal + b3_kb*sal ) &
+           + ( c1_kb + c2_kb*sqsal + c3_kb*sal )*logtk + d_kb*sqsal*tk )
+!
+!  Calculate the equilibrium constant of water from Millero (1994).
+!   Units are ( mol/(kg.soln) )^2 on the total hydrogen ion scale.
+!   Kw = (OH-).(H+)
+!   Check value (Dickson and Goyet; tc=25,sal=35): ln(Kw) = -30.434
+!
+          kwater &
+           = 1.0e-14  *EXP( ( a1_kw/tk + a2_kw + a3_kw*logtk ) &
+           + ( b1_kw/tk + b2_kw + b3_kw*logtk )*sqsal + c_kw*sal )
+!
+!  Calculate the combinations (as in PPCO2).
+!
+          rek(ji,jj) = SQRT( kone / ktwo )
+          rg3(ji,jj) = SQRT( kone * ktwo )
+          rg1(ji,jj) = rg3(ji,jj) / kboric
+          rg2(ji,jj) = kwater / rg3(ji,jj)
+!
+!  Set Boron concentration according to the surface salinity
+!   (Boron accounts for about 10% of the total Alkalinity).
+!   Use relationship used by Peng et al, Tellus 39b 1987,
+!   originally from Culkin 1965:
+!
+!         boron(ji,jj) = 4.16e-4 * sal / 35.0   ! Correct value but not in UM
+          boron(ji,jj) = 4.106e-4 * sal / 35.0
+!
+!  tco2 and alk_t are got from the modelled tracers at the surface
+!   (converted to units of mol/kg).
+!
+          tco2(ji,jj)  = MAX( 0.0, tco2_tracer_bkg(ji,jj) * 0.001 / RHO_WATER_SI )
+          alk_t(ji,jj) = MAX( 0.0, alk_tracer_bkg(ji,jj)  * 0.001 / RHO_WATER_SI )
+!
+!  Set initial values for iterative method.
+!
+          xs(ji,jj) = 1.0
+          x2p(ji,jj) = 0.0
+          x3(ji,jj) = 0.0
+          xx(ji,jj) = 0.0
+          xf(ji,jj) = xs(ji,jj)
+!
+        ENDDO
+      ENDDO
+
+!
+!  Use the secant method of similar triangles to calculate the H+ ion
+!   concentration, [H] (though strictly xf refers to  rg3/[H] ).
+!   To aid vectorisation this has been written to iterate a set
+!   number of times (mxiter), large enough to bring every example to
+!   within the 1e-5 fractional change criterion.
+!
+!  The method follows that of Bacastow (1990) in separating the
+!   alkalinity into two parts: the major part is the carbonate alk
+!   which uses the [H] that we are trying to find, while the minor
+!   part (xw) is sum of the borate and water dissociation terms which
+!   use the previous approximation to [H].
+!   This leaves an easy-to-solve quadratic.
+!
+      DO iter=1,mxiter
+!
+
+        DO ji=1,jpi
+          DO jj=1,jpj
+!
+            IF((abs(xf(ji,jj)-xx(ji,jj))/xf(ji,jj))<1.e-5) CYCLE   !prevent iteration collapsing to 0.0/0.0
+
+
+            x1(ji,jj) = x2p(ji,jj)
+            x2(ji,jj) = x3(ji,jj)
+            x2p(ji,jj) = xx(ji,jj)
+            x3(ji,jj) = xf(ji,jj)
+
+            IF (iter<=2) THEN
+              xx(ji,jj) = xf(ji,jj)
+            ELSE
+              xx(ji,jj) = x2p(ji,jj) + ( x3(ji,jj) - x2p(ji,jj) )*( x1(ji,jj) - x2p(ji,jj) ) &
+              /( x1(ji,jj) - x2(ji,jj) - x2p(ji,jj) + x3(ji,jj) )
+            ENDIF
+
+
+!
+!  All the above has been setting up the similar triangles for the
+!   secant method and using it to find an improved interpolated
+!   estimate; it has been completely independent of the equations
+!   we are trying to solve.
+!
+!  The next bit, however, IS dependent on the equations.
+!   alkprime is the carbonate alkalinity, which is calculated
+!   from the total alkalinity by subtracting the other, smaller, terms,
+!   which use the interpolated previous approximation to xf (xx).
+!
+
+            alkprime = alk_t(ji,jj) - rg2(ji,jj)*xx(ji,jj) + rg3(ji,jj)/xx(ji,jj) &
+            - boron(ji,jj)/( 1.0 + rg1(ji,jj)/xx(ji,jj) )	    
+	    
+!
+!  Set up the quadratic terms.
+!
+
+            qa = 2.0*tco2(ji,jj) - alkprime
+            qb = -rek(ji,jj)*( alkprime - tco2(ji,jj) )
+            qc = -alkprime
+
+            xf(ji,jj) = ( -qb + SQRT( qb*qb - 4.0*qa*qc ) )*0.5/qa
+
+
+          ENDDO
+        ENDDO
+!     
+      ENDDO
+
+!
+      DO ji=1,jpi
+        DO jj=1,jpj
+!
+          xtrma = 1.0/( 1.0 + rek(ji,jj)*xf(ji,jj) + xf(ji,jj)*xf(ji,jj) )
+          xtrmb = rg2(ji,jj) + rg3(ji,jj)/(xf(ji,jj)*xf(ji,jj)) &
+          + rg1(ji,jj)*boron(ji,jj)/(( rg1(ji,jj) + xf(ji,jj) )*( rg1(ji,jj) + xf(ji,jj) ))
+!
+          dxda = 1.0/( tco2(ji,jj)*( rek(ji,jj) + 4.0*xf(ji,jj) &
+          + rek(ji,jj)*xf(ji,jj)*xf(ji,jj) )*xtrma*xtrma + xtrmb )
+          dxdc = -dxda*xf(ji,jj)*( rek(ji,jj) + 2.0*xf(ji,jj) )*xtrma
+!
+          dgdc = xtrma - dxdc*tco2(ji,jj)*( rek(ji,jj) + 2.0*xf(ji,jj) )*xtrma*xtrma
+          dgda = ( 1.0 - dxda*( tco2(ji,jj)*( rek(ji,jj) + 4.0*xf(ji,jj) )*xtrma &
+          + xtrmb ) )/(xf(ji,jj)*( rek(ji,jj) + 2.0*xf(ji,jj) ))
+!
+          fxa = 1.0e-3*RHO_WATER_SI*co2_sol(ji,jj)
+          DPDC(ji,jj) = dgdc/fxa
+          DPDA(ji,jj) = dgda/fxa
+!
+        ENDDO
+      ENDDO
+         
+      RETURN
+      END SUBROUTINE pco2_grad
+
+END MODULE asmpco2bal
Index: branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 8455)
+++ branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 8456)
@@ -162,4 +162,5 @@
                 IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1 )    ! SSH
                 IF( ln_logchlinc ) CALL logchl_asm_inc( nit000 - 1 )
+                IF( ln_fco2inc .OR. ln_pco2inc ) CALL pco2_asm_inc( nit000 - 1 )
              ENDIF
           ENDIF
Index: branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8455)
+++ branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8456)
@@ -255,4 +255,6 @@
       IF( lk_asminc .AND. ln_asmiau .AND. ln_logchlinc ) &
          &               CALL logchl_asm_inc( kstp )  ! logchl assimilation
+      IF( lk_asminc .AND. ln_asmiau .AND. ( ln_fco2inc .OR. ln_pco2inc ) ) &
+         &               CALL pco2_asm_inc( kstp )    ! pCO2/fCO2 assimilation
                          CALL trc_stp( kstp )         ! time-stepping
 #endif
