Index: /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90	(revision 9380)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90	(revision 9381)
@@ -62,5 +62,7 @@
       & asm_pco2_bal
    USE par_medusa             ! MEDUSA parameters
-   USE mocsy_mainmod, ONLY: & ! MEDUSA carbonate system
+   USE mocsy_wrapper, ONLY: & ! MEDUSA carbonate system
+      & mocsy_interface
+   USE mocsy_mainmod, ONLY: & ! fCO2 to pCO2 conversion
       & f2pCO2
    USE sms_medusa, ONLY:    & ! MEDUSA diagnostic variables
@@ -177,4 +179,6 @@
 #endif
 
+#  include "domzgr_substitute.h90"
+
 CONTAINS
 
@@ -215,9 +219,13 @@
       ENDIF
 
-      IF ( ( ln_balwri ).AND.( .NOT. ln_phytobal ).AND. &
-         & ( .NOT. ln_spco2inc ).AND.( .NOT. ln_sfco2inc ).AND.( .NOT. ln_pphinc ) ) THEN
-         CALL ctl_warn( ' Balancing increments are only calculated for ocean colour', &
-            &           ' with ln_phytobal, or for pCO2, fCO2, or pH.',               &
-            &           ' Not the case, so ln_balwri will be set to .false.')
+      IF ( ( ln_balwri ).AND.                                         &
+         & ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. &
+         & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_schltotinc  ).AND. &
+         & ( .NOT. ln_slphytotinc ).AND.( .NOT. ln_slphydiainc ).AND. &
+         & ( .NOT. ln_slphynoninc ).AND.( .NOT. ln_plchltotinc ).AND. &
+         & ( .NOT. ln_pchltotinc  ).AND.( .NOT. ln_pphinc      ).AND. &
+         & ( .NOT. ln_spco2inc    ).AND.( .NOT. ln_sfco2inc    ) ) THEN
+         CALL ctl_warn( ' Balancing increments not being calculated', &
+            &           ' so ln_balwri will be set to .false.')
          ln_balwri = .FALSE.
       ENDIF
@@ -664,32 +672,46 @@
          CALL iom_rstput( kt, kt, inum, 'rdastp' , zdate   )
 
-         IF ( ln_slchltotinc ) THEN
+         IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
+            & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
+            & ln_slphynoninc ) THEN
 #if defined key_medusa
-            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chn', phyto2d_balinc(:,:,:,jpchn) )
-            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chd', phyto2d_balinc(:,:,:,jpchd) )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_chn', phyto2d_balinc(:,:,:,jpchn) )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_chd', phyto2d_balinc(:,:,:,jpchd) )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_phn', phyto2d_balinc(:,:,:,jpphn) )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_phd', phyto2d_balinc(:,:,:,jpphd) )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_pds', phyto2d_balinc(:,:,:,jppds) )
             IF ( ln_phytobal ) THEN
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phn', phyto2d_balinc(:,:,:,jpphn) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phd', phyto2d_balinc(:,:,:,jpphd) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_pds', phyto2d_balinc(:,:,:,jppds) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zmi', phyto2d_balinc(:,:,:,jpzmi) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zme', phyto2d_balinc(:,:,:,jpzme) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_din', phyto2d_balinc(:,:,:,jpdin) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_sil', phyto2d_balinc(:,:,:,jpsil) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_fer', phyto2d_balinc(:,:,:,jpfer) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_det', phyto2d_balinc(:,:,:,jpdet) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dtc', phyto2d_balinc(:,:,:,jpdtc) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dic', phyto2d_balinc(:,:,:,jpdic) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_alk', phyto2d_balinc(:,:,:,jpalk) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_oxy', phyto2d_balinc(:,:,:,jpoxy) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_zmi', phyto2d_balinc(:,:,:,jpzmi) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_zme', phyto2d_balinc(:,:,:,jpzme) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_din', phyto2d_balinc(:,:,:,jpdin) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_sil', phyto2d_balinc(:,:,:,jpsil) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_fer', phyto2d_balinc(:,:,:,jpfer) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_det', phyto2d_balinc(:,:,:,jpdet) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_dtc', phyto2d_balinc(:,:,:,jpdtc) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_dic', phyto2d_balinc(:,:,:,jpdic) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jpalk) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_oxy', phyto2d_balinc(:,:,:,jpoxy) )
             ENDIF
 #elif defined key_hadocc
-            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phy', phyto2d_balinc(:,:,:,jp_had_phy) )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_phy', phyto2d_balinc(:,:,:,jp_had_phy) )
             IF ( ln_phytobal ) THEN
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_nut', phyto2d_balinc(:,:,:,jp_had_nut) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zoo', phyto2d_balinc(:,:,:,jp_had_zoo) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_det', phyto2d_balinc(:,:,:,jp_had_pdn) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dic', phyto2d_balinc(:,:,:,jp_had_dic) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_alk', phyto2d_balinc(:,:,:,jp_had_alk) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_nut', phyto2d_balinc(:,:,:,jp_had_nut) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_zoo', phyto2d_balinc(:,:,:,jp_had_zoo) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_det', phyto2d_balinc(:,:,:,jp_had_pdn) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_dic', phyto2d_balinc(:,:,:,jp_had_dic) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jp_had_alk) )
             ENDIF
+#endif
+         ENDIF
+
+         IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
+#if defined key_medusa
+            CALL iom_rstput( kt, kt, inum, 'phy3d_chn', phyto3d_balinc(:,:,:,jpchn) )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_chd', phyto3d_balinc(:,:,:,jpchd) )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_phn', phyto3d_balinc(:,:,:,jpphn) )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_phd', phyto3d_balinc(:,:,:,jpphd) )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_pds', phyto3d_balinc(:,:,:,jppds) )
+#elif defined key_hadocc
+            CALL iom_rstput( kt, kt, inum, 'phy3d_phy', phyto3d_balinc(:,:,:,jp_had_phy) )
 #endif
          ENDIF
@@ -697,17 +719,27 @@
          IF ( ln_spco2inc ) 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) )
+            CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jpdic) )
+            CALL iom_rstput( kt, kt, inum, 'pco2_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) )
+            CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jp_had_dic) )
+            CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jp_had_alk) )
 #endif
          ELSE IF ( ln_sfco2inc ) 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) )
+            CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jpdic) )
+            CALL iom_rstput( kt, kt, inum, 'fco2_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) )
+            CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jp_had_dic) )
+            CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jp_had_alk) )
+#endif
+         ENDIF
+
+         IF ( ln_pphinc ) THEN
+#if defined key_medusa
+            CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jpdic) )
+            CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jpalk) )
+#elif defined key_hadocc
+            CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jp_had_dic) )
+            CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jp_had_alk) )
 #endif
          ENDIF
@@ -944,11 +976,14 @@
       REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights
       !
-      INTEGER                          :: ji, jj, jk ! Loop counters
-      INTEGER                          :: it         ! Index
-      REAL(wp)                         :: zincwgt    ! IAU weight for time step
-      REAL(wp)                         :: zfrac_chn  ! Fraction of jpchn
-      REAL(wp)                         :: zfrac_chd  ! Fraction of jpchd
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: chl_inc    ! Chlorophyll increments
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: bkg_chl    ! Chlorophyll background
+      INTEGER                          :: ji, jj, jk   ! Loop counters
+      INTEGER                          :: it           ! Index
+      REAL(wp)                         :: zincwgt      ! IAU weight for timestep
+      REAL(wp)                         :: zfrac_chn    ! Fraction of jpchn
+      REAL(wp)                         :: zfrac_chd    ! Fraction of jpchd
+      REAL(wp)                         :: zrat_phn_chn ! jpphn:jpchn ratio
+      REAL(wp)                         :: zrat_phd_chd ! jpphd:jpchd ratio
+      REAL(wp)                         :: zrat_pds_chd ! jppds:jpchd ratio
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: chl_inc      ! Chlorophyll increments
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: bkg_chl      ! Chlorophyll background
       !!------------------------------------------------------------------------
 
@@ -1006,4 +1041,10 @@
                      phyto3d_balinc(ji,jj,jk,jpchn) = chl_inc(ji,jj,jk) * zfrac_chn
                      phyto3d_balinc(ji,jj,jk,jpchd) = chl_inc(ji,jj,jk) * zfrac_chd
+                     zrat_phn_chn = tracer_bkg(ji,jj,jk,jpphn) / tracer_bkg(ji,jj,jk,jpchn)
+                     zrat_phd_chd = tracer_bkg(ji,jj,jk,jpphd) / tracer_bkg(ji,jj,jk,jpchd)
+                     zrat_pds_chd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpchd)
+                     phyto3d_balinc(ji,jj,jk,jpphn) = phyto3d_balinc(ji,jj,jk,jpchn) * zrat_phn_chn
+                     phyto3d_balinc(ji,jj,jk,jpphd) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_phd_chd
+                     phyto3d_balinc(ji,jj,jk,jppds) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_pds_chd
                   ENDIF
                END DO
@@ -1192,6 +1233,6 @@
 
 #if defined key_medusa
-         ! Account for logchl balancing if required
-         IF ( ln_slchltotinc .AND. ln_phytobal ) THEN
+         ! Account for phytoplankton balancing if required
+         IF ( ln_phytobal ) THEN
             dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto2d_balinc(:,:,1,jpdic)
             alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto2d_balinc(:,:,1,jpalk)
@@ -1206,6 +1247,6 @@
 
 #elif defined key_hadocc
-         ! Account for slchltot balancing if required
-         IF ( ln_slchltotinc .AND. ln_phytobal ) THEN
+         ! Account for phytoplankton balancing if required
+         IF ( ln_phytobal ) THEN
             dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + phyto2d_balinc(:,:,1,jp_had_dic)
             alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + phyto2d_balinc(:,:,1,jp_had_alk)
@@ -1379,44 +1420,209 @@
       !! ** Purpose : Apply the pH assimilation increments.
       !!
-      !! ** Method  : Calculate increments to state variables using carbon
-      !!              balancing.
+      !! ** Method  : Calculate increments to DIC and alkalinity from pH
+      !!              Use a similar approach to the pCO2 scheme
       !!              Direct initialization or Incremental Analysis Updating.
       !!
       !! ** Action  : 
       !!------------------------------------------------------------------------
-      INTEGER, INTENT(IN)                   :: kt           ! Current time step
-      LOGICAL, INTENT(IN)                   :: ll_asmdin    ! Flag for direct initialisation
-      LOGICAL, INTENT(IN)                   :: ll_asmiau    ! Flag for incremental analysis update
-      INTEGER, INTENT(IN)                   :: kcycper      ! Dimension of pwgtiau
-      REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights
-      LOGICAL, INTENT(IN)                   :: ll_trainc    ! Flag for T/S increments
+      INTEGER,                          INTENT(IN) :: kt        ! Current time step
+      LOGICAL,                          INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
+      LOGICAL,                          INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
+      INTEGER,                          INTENT(IN) :: kcycper   ! Dimension of pwgtiau
+      REAL(wp), DIMENSION(kcycper),     INTENT(IN) :: pwgtiau   ! IAU weights
+      LOGICAL,                          INTENT(IN) :: ll_trainc ! Flag for T/S increments
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
-      !!------------------------------------------------------------------------
-      
-      CALL ctl_stop( ' pH balancing not yet implemented' )
-      
-      ! See solve_at_general line 281 of mocsy_phsolvers.F90
-      !
-      ! Or, call to mocsy_interface line 158 of carb_chem.F90
-      ! Input variables (rest are output) are:
-      ! ztmp(ji,jj),zsal(ji,jj),zalk(ji,jj), &
-      ! zdic(ji,jj),zsil(ji,jj),zpho(ji,jj), &
-      ! f_pp0(ji,jj),fsdept(ji,jj,jk),       &
-      ! gphit(ji,jj),f_kw660(ji,jj),         &
-      ! f_xco2a(ji,jj),1
-      !
-      ! ztmp    = tsn(:,:,:,jp_tem)
-      ! zsal    = tsn(:,:,:,jp_sal)
-      ! zalk    = trn(:,:,:,jpalk)
-      ! zdic    = trn(:,:,:,jpdic)
-      ! zsil    = trn(:,:,:,jpsil)
-      ! zpho    = trn(:,:,:,jpdin) / 16.0
-      ! f_pp0   = 1.0
-      ! fsdept  = fsdept(:,:,:)
-      ! gphit   = gphit(:,:)
-      ! f_kw660 = 1.0
-      ! f_xco2a = f_xco2a(:,:)
-      
+      
+      REAL(wp)                         :: zsearch = 10.0 ! Increment to DIC/alk in pH calculation
+      REAL(wp)                         :: DIC_IN, ALK_IN ! DIC/alk in pH calculation
+      REAL(wp)                         :: PH_OUT         ! pH from pH calculation
+      REAL(wp)                         :: ph_bkg         ! pH from pH calculation
+      REAL(wp)                         :: ph_dic         ! pH from pH calculation
+      REAL(wp)                         :: ph_alk         ! pH from pH calculation
+      REAL(wp)                         :: dph_ddic       ! Change in pH wrt DIC
+      REAL(wp)                         :: dph_dalk       ! Change in pH wrt alk
+      REAL(wp)                         :: weight         ! Increment weighting
+      REAL(wp)                         :: zincwgt        ! IAU weight for current time step
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: dic_bkg_temp   ! DIC background
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: alk_bkg_temp   ! Alkalinity background
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: din_bkg_temp   ! DIN background
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: sil_bkg_temp   ! Silicate background
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: tem_bkg_temp   ! Temperature background
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: sal_bkg_temp   ! Salinity background
+      REAL(wp), DIMENSION(19)          :: dummy_out      ! Dummy variables
+      INTEGER                          :: ji, jj, jk, jx ! Loop counters
+      INTEGER                          :: it             ! Index
+      !!------------------------------------------------------------------------
+
+#if ! (defined key_medusa && defined key_foam_medusa)
+      CALL ctl_stop( ' pH balancing only implemented for MEDUSA' )
+#else
+
+      IF ( kt <= nit000 ) THEN
+
+         !----------------------------------------------------------------------
+         ! DIC and alkalinity balancing
+         !----------------------------------------------------------------------
+
+         ! Account for physics assimilation if required
+         IF ( ll_trainc ) THEN
+            tem_bkg_temp(:,:,:) = tsn(:,:,:,1) + pt_bkginc(:,:,:)
+            sal_bkg_temp(:,:,:) = tsn(:,:,:,2) + ps_bkginc(:,:,:)
+         ELSE
+            tem_bkg_temp(:,:,:) = tsn(:,:,:,1)
+            sal_bkg_temp(:,:,:) = tsn(:,:,:,2)
+         ENDIF
+
+         ! Account for phytoplankton balancing if required
+         IF ( ln_phytobal ) THEN
+            dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) + phyto2d_balinc(:,:,:,jpdic)
+            alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) + phyto2d_balinc(:,:,:,jpalk)
+            din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) + phyto2d_balinc(:,:,:,jpdin)
+            sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) + phyto2d_balinc(:,:,:,jpsil)
+         ELSE
+            dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic)
+            alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk)
+            din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin)
+            sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil)
+         ENDIF
+         
+         ! Account for pCO2 balancing if required
+         IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
+            dic_bkg_temp(:,:,:) = dic_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpdic)
+            alk_bkg_temp(:,:,:) = alk_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpalk)
+         ENDIF
+         
+         ! Loop over grid points and calculate dpH/dDIC and dpH/dAlk
+         ! This requires three calls to the MOCSY carbonate package
+         ! One to find pH at background DIC and alkalinity
+         ! One to find pH when DIC is increased by 10
+         ! One to find pH when alkalinity is increased by 10
+         ! Then calculate the gradients
+
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+
+                  IF ( pph_bkginc(ji,jj,jk) /= 0.0 ) THEN
+
+                     DO jx = 1, 3
+                        IF ( jx == 1 ) THEN
+                           DIC_IN = dic_bkg_temp(ji,jj,jk)
+                           ALK_IN = alk_bkg_temp(ji,jj,jk)
+                        ELSE IF ( jx == 2 ) THEN
+                           DIC_IN = dic_bkg_temp(ji,jj,jk) + zsearch
+                           ALK_IN = alk_bkg_temp(ji,jj,jk)
+                        ELSE IF ( jx == 3 ) THEN
+                           DIC_IN = dic_bkg_temp(ji,jj,jk)
+                           ALK_IN = alk_bkg_temp(ji,jj,jk) + zsearch
+                        ENDIF
+
+                        CALL mocsy_interface(tsn(ji,jj,jk,jp_tem),         & ! IN:  temperature
+                           &                 tsn(ji,jj,jk,jp_sal),         & !      salinity
+                           &                 ALK_IN,                       & !      alkalinity
+                           &                 DIC_IN,                       & !      DIC
+                           &                 sil_bkg_temp(ji,jj,jk),       & !      silicate
+                           &                 din_bkg_temp(ji,jj,jk)/16.0,  & !      phosphate (Redfield from DIN)
+                           &                 1.0,                          & !      atmospheric pressure (dummy)
+                           &                 fsdept(ji,jj,jk),             & !      depth
+                           &                 gphit(ji,jj),                 & !      latitude
+                           &                 1.0,                          & !      gas transfer (dummy)
+                           &                 1.0,                          & !      atmospheric xCO2 (dummy)
+                           &                 1,                            & !      number of points
+                           &                 PH_OUT,                       & ! OUT: pH
+                           &                 dummy_out(1),  dummy_out(2),  & !      stuff we don't care about
+                           &                 dummy_out(3),  dummy_out(4),  &
+                           &                 dummy_out(5),  dummy_out(6),  &
+                           &                 dummy_out(7),  dummy_out(8),  &
+                           &                 dummy_out(9),  dummy_out(10), &
+                           &                 dummy_out(11), dummy_out(12), &
+                           &                 dummy_out(13), dummy_out(14), &
+                           &                 dummy_out(15), dummy_out(16), &
+                           &                 dummy_out(17), dummy_out(18), &
+                           &                 dummy_out(19) )
+
+                        IF ( jx == 1 ) THEN
+                           ph_bkg = PH_OUT
+                        ELSE IF ( jx == 2 ) THEN
+                           ph_dic = PH_OUT
+                        ELSE IF ( jx == 3 ) THEN
+                           ph_alk = PH_OUT
+                        ENDIF
+                     END DO
+
+                     dph_ddic = (ph_dic - ph_bkg) / zsearch
+                     dph_dalk = (ph_alk - ph_bkg) / zsearch
+                     weight   = pph_bkginc(ji,jj,jk) / ( (dph_ddic*dph_ddic) + (dph_dalk*dph_dalk) )
+
+                     ph_balinc(ji,jj,jk,jpdic) = weight * dph_ddic
+                     ph_balinc(ji,jj,jk,jpalk) = weight * dph_dalk
+                     
+                  ENDIF
+                  
+               END DO
+            END DO
+         END DO
+
+      ENDIF
+      
+      IF ( ll_asmiau ) THEN
+
+         !--------------------------------------------------------------------
+         ! Incremental Analysis Updating
+         !--------------------------------------------------------------------
+
+         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
+
+            it = kt - nit000 + 1
+            zincwgt = pwgtiau(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,*) 
+               WRITE(numout,*) 'ph_asm_inc : pH IAU at time step = ', &
+                  &  kt,' with IAU weight = ', pwgtiau(it)
+               WRITE(numout,*) '~~~~~~~~~~~~'
+            ENDIF
+
+            ! Update the biogeochemical variables
+            ! Add directly to trn and trb, rather than to tra, because tra gets
+            ! reset to zero at the start of trc_stp, called after this routine
+            DO jk = 1, jpkm1
+               trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + &
+                  &                          ph_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
+               trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + &
+                  &                          ph_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
+            END DO
+
+            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
+            ! which is called at end of model run
+
+         ENDIF
+
+      ELSEIF ( ll_asmdin ) THEN 
+
+         !--------------------------------------------------------------------
+         ! Direct Initialization
+         !--------------------------------------------------------------------
+         
+         IF ( kt == nitdin_r ) THEN
+
+            neuler = 0                    ! Force Euler forward step
+
+            ! 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 pH assimilation', &
+               &           ' Background state is taken from model rather than background file' )
+            trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
+               &                         ph_balinc(:,:,:,jp_msa0:jp_msa1)
+            trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
+ 
+            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
+            ! which is called at end of model run
+         ENDIF
+         !
+      ENDIF
+#endif      
       !
    END SUBROUTINE ph_asm_inc
Index: /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 9380)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 9381)
@@ -1239,5 +1239,5 @@
       
       ! Surface pCO2/fCO2 next
-      IF( ln_sfco2inc .OR. ln_spco2inc ) THEN
+      IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
          CALL pco2_asm_inc( kt, ln_asmdin, ln_asmiau, icycper, wgtiau, &
             &               ln_trainc, t_bkginc, s_bkginc )
@@ -1245,5 +1245,5 @@
       
       ! Profile pH next
-      IF( ln_pphinc ) THEN
+      IF ( ln_pphinc ) THEN
          CALL ph_asm_inc( kt, ln_asmdin, ln_asmiau, icycper, wgtiau, &
             &             ln_trainc, t_bkginc, s_bkginc )
