Index: branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90
===================================================================
--- branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90	(revision 7937)
+++ branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90	(revision 7938)
@@ -19,4 +19,14 @@
    PUBLIC   bio_medusa_alloc     ! called by trcini.F90
 
+   !! model state variables
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::    zchn,zchd,zphn,zphd,zpds,zzmi
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::    zzme,zdet,zdtc,zdin,zsil,zfer
+# if defined key_roam
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::    zdic, zalk, zoxy
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::    ztmp, zsal
+# endif
+# if defined key_mocsy
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::    zpho
+# endif
 
    !! AXY (01/03/10): add in mixed layer PP diagnostics
@@ -62,5 +72,19 @@
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fsil_prod,fsil_cons
 
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_xco2a
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_ph,f_pco2w,f_h2co3
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_hco3,f_co3,f_co2flux
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_TDIC,f_TALK,f_dcf,f_henry
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_pp0
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_kw660,f_o2flux,f_o2sat
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_omcal,f_omarg
+
+   !! AXY (23/06/15): additional diagnostics for MOCSY and oxygen
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_fco2w,f_BetaD,f_rhosw,f_opres
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_insitut,f_pco2atm,f_fco2atm
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_schmidtco2,f_kwco2,f_K0
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_co2starair,f_dpco2,f_kwo2
+
+   INTEGER, ALLOCATABLE, DIMENSION(:,:)  :: iters 
 
    !! Carbon, alkalinity production and consumption
@@ -172,5 +196,12 @@
       !!-------------------------------------------------------------------
       !
-      ALLOCATE(fprn_ml(jpi,jpj),fprd_ml(jpi,jpj),                     &
+      ALLOCATE(zchn(jpi,jpj),zchd(jpi,jpj),zphn(jpi,jpj),             &
+               zphd(jpi,jpj),zpds(jpi,jpj),zzmi(jpi,jpj),             &
+               zzme(jpi,jpj),zdet(jpi,jpj),zdtc(jpi,jpj),             &
+               zdin(jpi,jpj),zsil(jpi,jpj),zfer(jpi,jpj),             &
+               zdic(jpi,jpj),zalk(jpi,jpj),zoxy(jpi,jpj),             &
+               ztmp(jpi,jpj),zsal(jpi,jpj),                           &
+               zpho(jpi,jpj),                                         &
+               fprn_ml(jpi,jpj),fprd_ml(jpi,jpj),                     &
                xFree(jpi,jpj),                                        &
                fslownflux(jpi,jpj),fslowcflux(jpi,jpj),               &
@@ -201,5 +232,19 @@
                fnit_prod(jpi,jpj),fnit_cons(jpi,jpj),                 &
                fsil_prod(jpi,jpj),fsil_cons(jpi,jpj),                 &
+               f_xco2a(jpi,jpj),                                      &
+               f_ph(jpi,jpj),f_pco2w(jpi,jpj),f_h2co3(jpi,jpj),       &
+               f_hco3(jpi,jpj),f_co3(jpi,jpj),f_co2flux(jpi,jpj),     &
+               f_TDIC(jpi,jpj),f_TALK(jpi,jpj),f_dcf(jpi,jpj),        &
+               f_henry(jpi,jpj),                                      &
+               f_pp0(jpi,jpj),                                        &
+               f_kw660(jpi,jpj),f_o2flux(jpi,jpj),f_o2sat(jpi,jpj),   &
                f_omcal(jpi,jpj),f_omarg(jpi,jpj),                     &
+               f_fco2w(jpi,jpj),f_BetaD(jpi,jpj),f_rhosw(jpi,jpj),    &
+               f_opres(jpi,jpj),f_insitut(jpi,jpj),                   &
+               f_pco2atm(jpi,jpj),f_fco2atm(jpi,jpj),                 &
+               f_schmidtco2(jpi,jpj),f_kwco2(jpi,jpj),f_K0(jpi,jpj),  &
+               f_co2starair(jpi,jpj),f_dpco2(jpi,jpj),                &
+               f_kwo2(jpi,jpj),                                       &
+               iters(jpi,jpj),                                        &
                fcomm_resp(jpi,jpj),                                   &
                fcar_prod(jpi,jpj),fcar_cons(jpi,jpj),                 &
Index: branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90
===================================================================
--- branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90	(revision 7938)
+++ branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90	(revision 7938)
@@ -0,0 +1,249 @@
+MODULE carb_chem_mod
+   !!======================================================================
+   !!                         ***  MODULE carb_chem_mod  ***
+   !! Calculate the carbon chemistry for the whole ocean
+   !!======================================================================
+   !! History :
+   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!----------------------------------------------------------------------
+#if defined key_roam
+   !!----------------------------------------------------------------------
+   !!                                                   key_roam?
+   !!----------------------------------------------------------------------
+
+   IMPLICIT NONE
+   PRIVATE
+      
+   PUBLIC   carb_chem        ! Called in trcbio_medusa.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Id$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE carb_chem( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE carb_chem  ***
+      !! This called from TRC_BIO_MEDUSA and 
+      !!  - ...
+      !!----------------------------------------------------------------------
+      USE bio_medusa_mod,    ONLY: iters, f_BetaD, f_dcf, f_fco2atm,      &
+                                   f_co2flux, f_co2starair, f_co3,        &
+                                   f_fco2w, f_dpco2, f_h2co3, f_hco3,     &
+                                   f_insitut, f_K0, f_kw660, f_kwco2,     &
+                                   f_omarg, f_omcal, f_opres, f_pco2atm,  &
+                                   f_pco2w, f_ph, f_pp0, f_rhosw,         &
+                                   f_schmidtco2, f_TALK, f_TDIC, f_xco2a, &
+                                   zalk, zdic, zpho, zsal, zsil, ztmp 
+      USE dom_oce,           ONLY: gdept_0, gdept_n, gdepw_n, gphit,      &
+                                   mbathy, tmask
+      USE in_out_manager,    ONLY: lwp, numout
+      USE mocsy_wrapper,     ONLY: mocsy_interface
+      USE oce,               ONLY: PCO2a_in_cpl, tsb, tsn
+      USE par_kind,          ONLY: wp
+      USE par_medusa,        ONLY: jpalk, jpdic, jpdin, jpsil 
+      USE par_oce,           ONLY: jp_sal, jp_tem, jpi, jpim1, jpj,       &
+                                   jpjm1, jpk
+      USE sbc_oce,           ONLY: lk_oasis
+      USE sms_medusa,        ONLY: f2_ccd_arg, f2_ccd_cal, f3_co3,        &
+                                   f3_h2co3, f3_hco3, f3_omarg,           &
+                                   f3_omcal, f3_pH
+      USE trc,               ONLY: trn
+      USE trcco2_medusa,     ONLY: trc_co2_medusa
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+ 
+      !! time (integer timestep)
+      INTEGER, INTENT( in ) ::    kt
+
+      !! Flags to help with calculating the position of the CCD
+      INTEGER, DIMENSION(jpi,jpj) ::     i2_omcal,i2_omarg
+
+      !! temporary variables
+      REAL(wp) ::    fq0,fq1,fq2,fq3,fq4
+
+      INTEGER :: ji, jj, jk
+
+      !!----------------------------------------------------------------------
+      !! Calculate the carbonate chemistry for the whole ocean on the first
+      !! simulation timestep and every month subsequently; the resulting 3D
+      !! field of omega calcite is used to determine the depth of the CCD
+      !!----------------------------------------------------------------------
+      !!
+!      IF(lwp) WRITE(numout,*)                                               &
+!              ' MEDUSA calculating all carbonate chemistry at kt =', kt
+!         CALL flush(numout)
+      !! blank flags
+      i2_omcal(:,:) = 0
+      i2_omarg(:,:) = 0
+
+      !! loop over 3D space
+      DO jk = 1,jpk
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               !! OPEN wet point IF..THEN loop
+               if (tmask(ji,jj,jk).eq.1) then
+                  IF (lk_oasis) THEN
+                     !! use 2D atm xCO2 from atm coupling
+                     f_xco2a(ji,jj) = PCO2a_in_cpl(ji,jj)
+                  ENDIF
+                  !! Do carbonate chemistry
+                  !!
+                  !! Set up required state variables.
+                  !! dissolved inorganic carbon
+                  zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic))
+                  !! alkalinity
+                  zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk))
+                  !! temperature
+                  ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem)
+                  !! salinity
+                  zsal(ji,jj) = tsn(ji,jj,jk,jp_sal)
+#  if defined key_mocsy
+                  !! silicic acid
+                  zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil))
+                  !! phosphate via DIN and Redfield
+                  zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 
+#  endif
+		  !!
+		  !! AXY (28/02/14): check input fields
+		  if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then
+                     IF (lwp) WRITE(numout,*) ' carb_chem: T WARNING 3D, ',  &
+                        tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem),          &
+                        ' at (', ji, ',', jj, ',', jk, ') at time', kt
+	             IF(lwp) WRITE(numout,*)                                 &
+                        ' carb_chem: T SWITCHING 3D, ',                 &
+			tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem)
+                     ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem)     !! temperature
+                  endif
+		  if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then
+                     IF (lwp) WRITE(numout,*) ' carb_chem: S WARNING 3D, ',  &
+                        tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal),          &
+                        ' at (', ji, ',', jj, ',', jk, ') at time', kt
+                  endif
+! MAYBE BUT A BREAK IN HERE - marc
+                  !!
+                  !! Blank input variables not used at this stage (they 
+                  !! relate to air-sea flux)
+                  f_kw660(ji,jj) = 1.0
+                  f_pp0(ji,jj)   = 1.0
+                  !!
+                  !! calculate carbonate chemistry at grid cell midpoint
+#  if defined key_mocsy
+                  !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate
+                  !!                 chemistry package
+                  CALL mocsy_interface(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,f_ph(ji,jj),        &
+                                       f_pco2w(ji,jj),f_fco2w(ji,jj),       &
+                                       f_h2co3(ji,jj),f_hco3(ji,jj),        &
+                                       f_co3(ji,jj),f_omarg(ji,jj),         &
+                                       f_omcal(ji,jj),f_BetaD(ji,jj),       &
+                                       f_rhosw(ji,jj),f_opres(ji,jj),       &
+                                       f_insitut(ji,jj),f_pco2atm(ji,jj),   &
+                                       f_fco2atm(ji,jj),f_schmidtco2(ji,jj),&
+                                       f_kwco2(ji,jj),f_K0(ji,jj),          &
+                                       f_co2starair(ji,jj),f_co2flux(ji,jj),& 
+                                       f_dpco2(ji,jj))
+                  !!
+                  !! mmol / m3 -> umol / kg
+                  f_TDIC(ji,jj) = (zdic(ji,jj) / f_rhosw(ji,jj)) * 1000.
+                  !! meq / m3 -> ueq / kg
+                  f_TALK(ji,jj) = (zalk(ji,jj) / f_rhosw(ji,jj)) * 1000.
+                  f_dcf(ji,jj)  = f_rhosw(ji,jj)
+#  else
+                  !! AXY (22/06/15): use old PML carbonate chemistry 
+                  !! package (the MEDUSA-2 default)
+                  CALL trc_co2_medusa(ztmp(ji,jj),zsal(ji,jj),zdic(ji,jj),  &
+                                      zalk(ji,jj),fsdept(ji,jj,jk),         &
+                                      f_kw660(ji,jj),f_xco2a(ji,jj),        &
+                                      f_ph(ji,jj),                          &
+                                      f_pco2w(ji,jj),f_h2co3(ji,jj),        &
+                                      f_hco3(ji,jj),f_co3(ji,jj),           &
+                                      f_omcal(ji,jj),f_omarg(ji,jj),        &
+                                      f_co2flux(ji,jj),f_TDIC(ji,jj),       &
+                                      f_TALK(ji,jj), f_dcf(ji,jj),          &
+                                      f_henry(ji,jj), iters(ji,jj))
+                  !! 
+                  !! AXY (28/02/14): check output fields
+                  IF (iters(ji,jj) .eq. 25) THEN
+                     IF(lwp) WRITE(numout,*)                                &
+                        ' carb_chem: 3D ITERS WARNING, ',                   &
+                        iters(ji,jj), ' AT (', ji, ', ', jj, ', ',          &
+                        jk, ') AT ', kt
+                  ENDIF
+#  endif
+                  !!
+                  !! store 3D outputs
+                  f3_pH(ji,jj,jk)    = f_ph(ji,jj)
+                  f3_h2co3(ji,jj,jk) = f_h2co3(ji,jj)
+                  f3_hco3(ji,jj,jk)  = f_hco3(ji,jj)
+                  f3_co3(ji,jj,jk)   = f_co3(ji,jj)
+                  f3_omcal(ji,jj,jk) = f_omcal(ji,jj)
+                  f3_omarg(ji,jj,jk) = f_omarg(ji,jj)
+                  !!
+                  !! CCD calculation: calcite
+                  if (i2_omcal(ji,jj) == 0 .and. f_omcal(ji,jj) < 1.0) then
+                     if (jk .eq. 1) then
+                        f2_ccd_cal(ji,jj) = fsdept(ji,jj,jk)
+                     else
+                        fq0 = f3_omcal(ji,jj,jk-1) - f_omcal(ji,jj)
+                        fq1 = f3_omcal(ji,jj,jk-1) - 1.0
+                        fq2 = fq1 / (fq0 + tiny(fq0))
+                        fq3 = fsdept(ji,jj,jk) - fsdept(ji,jj,jk-1)
+                        fq4 = fq2 * fq3
+                        f2_ccd_cal(ji,jj) = fsdept(ji,jj,jk-1) + fq4
+                     endif
+                     i2_omcal(ji,jj)   = 1
+                  endif
+                  if ( i2_omcal(ji,jj) == 0 .and. jk == mbathy(ji,jj) ) then
+                     !! reached seafloor and still no dissolution; set 
+                     !! to seafloor (W-point)
+                     f2_ccd_cal(ji,jj) = fsdepw(ji,jj,jk+1)
+                     i2_omcal(ji,jj)   = 1
+                  endif
+                  !!
+                  !! CCD calculation: aragonite
+                  if (i2_omarg(ji,jj) == 0 .and. f_omarg(ji,jj) < 1.0) then
+                     if (jk .eq. 1) then
+                        f2_ccd_arg(ji,jj) = fsdept(ji,jj,jk)
+                     else
+                        fq0 = f3_omarg(ji,jj,jk-1) - f_omarg(ji,jj)
+                        fq1 = f3_omarg(ji,jj,jk-1) - 1.0
+                        fq2 = fq1 / (fq0 + tiny(fq0))
+                        fq3 = fsdept(ji,jj,jk) - fsdept(ji,jj,jk-1)
+                        fq4 = fq2 * fq3
+                        f2_ccd_arg(ji,jj) = fsdept(ji,jj,jk-1) + fq4
+                     endif
+                     i2_omarg(ji,jj)   = 1
+                  endif
+                  if ( i2_omarg(ji,jj) == 0 .and. jk == mbathy(ji,jj) ) then
+                     !! reached seafloor and still no dissolution; set 
+                     !! to seafloor (W-point)
+                     f2_ccd_arg(ji,jj) = fsdepw(ji,jj,jk+1)
+                     i2_omarg(ji,jj)   = 1
+                  endif
+               endif
+            ENDDO
+         ENDDO
+      ENDDO
+
+   END SUBROUTINE carb_chem
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   None key_roam?
+   !!======================================================================
+CONTAINS
+   SUBROUTINE carb_chem( )                    ! Empty routine
+      WRITE(*,*) 'carb_chem: You should not have seen this print! error?'
+   END SUBROUTINE carb_chem
+#endif 
+
+   !!======================================================================
+END MODULE carb_chem_mod
Index: branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 7937)
+++ branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 7938)
@@ -98,6 +98,7 @@
       USE bio_medusa_mod
       USE bio_medusa_init_mod,        ONLY: bio_medusa_init
+      USE carb_chem_mod,              ONLY: carb_chem
+      USE bio_medusa_diag_slice_mod,  ONLY: bio_medusa_diag_slice
       USE bio_medusa_fin_mod,         ONLY: bio_medusa_fin
-      USE bio_medusa_diag_slice_mod,  ONLY: bio_medusa_diag_slice
 
       IMPLICIT NONE
@@ -160,14 +161,15 @@
       !!
       !! model state variables
-      REAL(wp), DIMENSION(jpi,jpj) ::    zchn,zchd,zphn,zphd,zpds,zzmi
-      REAL(wp), DIMENSION(jpi,jpj) ::    zzme,zdet,zdtc,zdin,zsil,zfer
-      REAL(wp) ::    zage
-# if defined key_roam
-      REAL(wp), DIMENSION(jpi,jpj) ::    zdic, zalk, zoxy
-      REAL(wp), DIMENSION(jpi,jpj) ::    ztmp, zsal
-# endif
-# if defined key_mocsy
-      REAL(wp), DIMENSION(jpi,jpj) ::    zpho
-# endif
+!      REAL(wp), DIMENSION(jpi,jpj) ::    zchn,zchd,zphn,zphd,zpds,zzmi
+!      REAL(wp), DIMENSION(jpi,jpj) ::    zzme,zdet,zdtc,zdin,zsil,zfer
+! zage doesn't seem to be used - marc 19/4/17
+!      REAL(wp) ::    zage
+!# if defined key_roam
+!      REAL(wp), DIMENSION(jpi,jpj) ::    zdic, zalk, zoxy
+!      REAL(wp), DIMENSION(jpi,jpj) ::    ztmp, zsal
+!# endif
+!# if defined key_mocsy
+!      REAL(wp), DIMENSION(jpi,jpj) ::    zpho
+!# endif
       !!
       !! integrated source and sink terms
@@ -305,15 +307,15 @@
       !!
       !! AXY (24/11/16): add xCO2 variable for atmosphere (what we actually have)
-      REAL(wp)                     ::    f_xco2a
-      REAL(wp), DIMENSION(jpi,jpj) ::    f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_co2flux
-      REAL(wp), DIMENSION(jpi,jpj) ::    f_TDIC, f_TALK, f_dcf, f_henry
-      REAL(wp), DIMENSION(jpi,jpj) ::    f_pp0
-      REAL(wp), DIMENSION(jpi,jpj) ::    f_kw660, f_o2flux, f_o2sat
+!      REAL(wp)                     ::    f_xco2a
+!      REAL(wp), DIMENSION(jpi,jpj) ::    f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_co2flux
+!      REAL(wp), DIMENSION(jpi,jpj) ::    f_TDIC, f_TALK, f_dcf, f_henry
+!      REAL(wp), DIMENSION(jpi,jpj) ::    f_pp0
+!      REAL(wp), DIMENSION(jpi,jpj) ::    f_kw660, f_o2flux, f_o2sat
       REAL(wp)                     ::    f_o2sat3
 !      REAL(wp), DIMENSION(jpi,jpj) ::    f_omcal, f_omarg
       !!
       !! AXY (23/06/15): additional diagnostics for MOCSY and oxygen
-      REAL(wp), DIMENSION(jpi,jpj) ::    f_fco2w, f_BetaD, f_rhosw, f_opres, f_insitut, f_pco2atm, f_fco2atm
-      REAL(wp), DIMENSION(jpi,jpj) ::    f_schmidtco2, f_kwco2, f_K0, f_co2starair, f_dpco2, f_kwo2
+!      REAL(wp), DIMENSION(jpi,jpj) ::    f_fco2w, f_BetaD, f_rhosw, f_opres, f_insitut, f_pco2atm, f_fco2atm
+!      REAL(wp), DIMENSION(jpi,jpj) ::    f_schmidtco2, f_kwco2, f_K0, f_co2starair, f_dpco2, f_kwo2
       !! jpalm 14-07-2016: convert CO2flux diag from mmol/m2/d to kg/m2/s
       REAL, PARAMETER :: weight_CO2_mol = 44.0095  !! g / mol
@@ -322,5 +324,5 @@
 
       !!
-      INTEGER, DIMENSION(jpi,jpj)  ::    iters
+!      INTEGER, DIMENSION(jpi,jpj)  ::    iters
       REAL(wp) ::    f_year
       INTEGER  ::    i_year
@@ -484,8 +486,8 @@
       if (iyr1 .le. 1) then
          !! before 1860
-         f_xco2a = hist_pco2(1)
+         f_xco2a(:,:) = hist_pco2(1)
       elseif (iyr2 .ge. 242) then
          !! after 2099
-         f_xco2a = hist_pco2(242)
+         f_xco2a(:,:) = hist_pco2(242)
       else
          !! just right
@@ -500,8 +502,8 @@
 #  endif
          fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3)
-         f_xco2a = fq4
+         f_xco2a(:,:) = fq4
       endif
 #  if defined key_axy_pi_co2
-      f_xco2a = 284.725          !! OCMIP pre-industrial pCO2
+      f_xco2a(:,:) = 284.725          !! OCMIP pre-industrial pCO2
 #  endif
       !! IF(lwp) WRITE(numout,*) ' MEDUSA nyear     =', nyear
@@ -513,5 +515,5 @@
       !! IF(lwp) WRITE(numout,*) ' MEDUSA fq2       =', fq2
       !! IF(lwp) WRITE(numout,*) ' MEDUSA fq3       =', fq3
-      IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2  =', f_xco2a
+      IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2  =', f_xco2a(1,1)
 # endif
 
@@ -544,129 +546,6 @@
          !! field of omega calcite is used to determine the depth of the CCD
          !!----------------------------------------------------------------------
-         !!
-         IF(lwp) WRITE(numout,*) ' MEDUSA calculating all carbonate chemistry at kt =', kt
-         CALL flush(numout)
-         !! blank flags
-         i2_omcal(:,:) = 0
-         i2_omarg(:,:) = 0
-         !! loop over 3D space
-         DO jk = 1,jpk
-            DO jj = 2,jpjm1
-               DO ji = 2,jpim1
-                  !! OPEN wet point IF..THEN loop
-                  if (tmask(ji,jj,jk).eq.1) then
-                     IF (lk_oasis) THEN
-                        f_xco2a = PCO2a_in_cpl(ji,jj)        !! use 2D atm xCO2 from atm coupling
-                     ENDIF
-                     !! do carbonate chemistry
-                     !!
-                     !! set up required state variables
-                     zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic)) !! dissolved inorganic carbon
-                     zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk)) !! alkalinity
-                     ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem)        !! temperature
-                     zsal(ji,jj) = tsn(ji,jj,jk,jp_sal)        !! salinity
-#  if defined key_mocsy
-                     zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil))        !! silicic acid
-                     zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 !! phosphate via DIN and Redfield
-#  endif
-		     !!
-		     !! AXY (28/02/14): check input fields
-		     if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then
-                        IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T WARNING 3D, ', &
-                        tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), ' at (',    &
-                        ji, ',', jj, ',', jk, ') at time', kt
-			IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 3D, ', &
-			tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem)
-                        ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem)     !! temperature
-                     endif
-		     if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then
-                        IF(lwp) WRITE(numout,*) ' trc_bio_medusa: S WARNING 3D, ', &
-                        tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (',    &
-                        ji, ',', jj, ',', jk, ') at time', kt
-                     endif
-                     !!
-                     !! blank input variables not used at this stage (they relate to air-sea flux)
-                     f_kw660(ji,jj) = 1.0
-                     f_pp0(ji,jj)   = 1.0
-                     !!
-                     !! calculate carbonate chemistry at grid cell midpoint
-#  if defined key_mocsy
-                     !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate
-                     !!                 chemistry package
-                     CALL mocsy_interface( ztmp(ji,jj), zsal(ji,jj), zalk(ji,jj), zdic(ji,jj), zsil(ji,jj), zpho(ji,jj),         &    ! inputs
-                     f_pp0(ji,jj), fsdept(ji,jj,jk), gphit(ji,jj), f_kw660(ji,jj), f_xco2a, 1,                  &    ! inputs
-                     f_ph(ji,jj), f_pco2w(ji,jj), f_fco2w(ji,jj), f_h2co3(ji,jj), f_hco3(ji,jj), f_co3(ji,jj), f_omarg(ji,jj),   &    ! outputs
-                     f_omcal(ji,jj), f_BetaD(ji,jj), f_rhosw(ji,jj), f_opres(ji,jj), f_insitut(ji,jj),             &    ! outputs
-                     f_pco2atm(ji,jj), f_fco2atm(ji,jj), f_schmidtco2(ji,jj), f_kwco2(ji,jj), f_K0(ji,jj),                &    ! outputs
-                     f_co2starair(ji,jj), f_co2flux(ji,jj), f_dpco2(ji,jj) )                                     ! outputs
-                     !!
-                     f_TDIC(ji,jj) = (zdic(ji,jj) / f_rhosw(ji,jj)) * 1000. ! mmol / m3 -> umol / kg
-                     f_TALK(ji,jj) = (zalk(ji,jj) / f_rhosw(ji,jj)) * 1000. !  meq / m3 ->  ueq / kg
-                     f_dcf(ji,jj)  = f_rhosw(ji,jj)
-#  else
-                     !! AXY (22/06/15): use old PML carbonate chemistry package (the
-                     !!                 MEDUSA-2 default)
-                     CALL trc_co2_medusa( ztmp(ji,jj), zsal(ji,jj), zdic(ji,jj), zalk(ji,jj), fsdept(ji,jj,jk), f_kw660(ji,jj),      &    ! inputs
-                     f_xco2a, f_ph(ji,jj), f_pco2w(ji,jj), f_h2co3(ji,jj), f_hco3(ji,jj), f_co3(ji,jj), f_omcal(ji,jj),   &    ! outputs
-                     f_omarg(ji,jj), f_co2flux(ji,jj), f_TDIC(ji,jj), f_TALK(ji,jj), f_dcf(ji,jj), f_henry(ji,jj), iters(ji,jj))      ! outputs
-                     !! 
-                     !! AXY (28/02/14): check output fields
-                     if (iters(ji,jj) .eq. 25) then
-                        IF(lwp) WRITE(numout,*) ' trc_bio_medusa: 3D ITERS WARNING, ', &
-                        iters(ji,jj), ' AT (', ji, ', ', jj, ', ', jk, ') AT ', kt
-                     endif
-#  endif
-                     !!
-                     !! store 3D outputs
-                     f3_pH(ji,jj,jk)    = f_ph(ji,jj)
-                     f3_h2co3(ji,jj,jk) = f_h2co3(ji,jj)
-                     f3_hco3(ji,jj,jk)  = f_hco3(ji,jj)
-                     f3_co3(ji,jj,jk)   = f_co3(ji,jj)
-                     f3_omcal(ji,jj,jk) = f_omcal(ji,jj)
-                     f3_omarg(ji,jj,jk) = f_omarg(ji,jj)
-                     !!
-                     !! CCD calculation: calcite
-                     if (i2_omcal(ji,jj) .eq. 0 .and. f_omcal(ji,jj) .lt. 1.0) then
-                        if (jk .eq. 1) then
-                           f2_ccd_cal(ji,jj) = fsdept(ji,jj,jk)
-                        else
-                           fq0 = f3_omcal(ji,jj,jk-1) - f_omcal(ji,jj)
-                           fq1 = f3_omcal(ji,jj,jk-1) - 1.0
-                           fq2 = fq1 / (fq0 + tiny(fq0))
-                           fq3 = fsdept(ji,jj,jk) - fsdept(ji,jj,jk-1)
-                           fq4 = fq2 * fq3
-                           f2_ccd_cal(ji,jj) = fsdept(ji,jj,jk-1) + fq4
-                        endif
-                        i2_omcal(ji,jj)   = 1
-                     endif
-                     if ( i2_omcal(ji,jj) .eq. 0 .and. jk .eq. mbathy(ji,jj) ) then
-                        !! reached seafloor and still no dissolution; set to seafloor (W-point)
-                        f2_ccd_cal(ji,jj) = fsdepw(ji,jj,jk+1)
-                        i2_omcal(ji,jj)   = 1
-                     endif
-                     !!
-                     !! CCD calculation: aragonite
-                     if (i2_omarg(ji,jj) .eq. 0 .and. f_omarg(ji,jj) .lt. 1.0) then
-                        if (jk .eq. 1) then
-                           f2_ccd_arg(ji,jj) = fsdept(ji,jj,jk)
-                        else
-                           fq0 = f3_omarg(ji,jj,jk-1) - f_omarg(ji,jj)
-                           fq1 = f3_omarg(ji,jj,jk-1) - 1.0
-                           fq2 = fq1 / (fq0 + tiny(fq0))
-                           fq3 = fsdept(ji,jj,jk) - fsdept(ji,jj,jk-1)
-                           fq4 = fq2 * fq3
-                           f2_ccd_arg(ji,jj) = fsdept(ji,jj,jk-1) + fq4
-                        endif
-                        i2_omarg(ji,jj)   = 1
-                     endif
-                     if ( i2_omarg(ji,jj) .eq. 0 .and. jk .eq. mbathy(ji,jj) ) then
-                        !! reached seafloor and still no dissolution; set to seafloor (W-point)
-                        f2_ccd_arg(ji,jj) = fsdepw(ji,jj,jk+1)
-                        i2_omarg(ji,jj)   = 1
-                     endif
-                  endif
-               ENDDO
-            ENDDO
-         ENDDO
+         CALL carb_chem( kt )
+
       ENDIF
 # endif
@@ -895,5 +774,5 @@
                   !!----------------------------------------------------------------------
                   IF (lk_oasis) THEN
-                     f_xco2a = PCO2a_in_cpl(ji,jj)        !! use 2D atm xCO2 from atm coupling
+                     f_xco2a(ji,jj) = PCO2a_in_cpl(ji,jj)        !! use 2D atm xCO2 from atm coupling
                   ENDIF
                   !!
@@ -936,5 +815,5 @@
                   !!                 zero in this call
                   CALL mocsy_interface( ztmp(ji,jj), zsal(ji,jj), zalk(ji,jj), zdic(ji,jj), zsil(ji,jj), zpho(ji,jj),        &  ! inputs
-                  f_pp0(ji,jj), 0.0, gphit(ji,jj), f_kw660(ji,jj), f_xco2a, 1,                   &  ! inputs
+                  f_pp0(ji,jj), 0.0, gphit(ji,jj), f_kw660(ji,jj), f_xco2a(ji,jj), 1,                   &  ! inputs
                   f_ph(ji,jj), f_pco2w(ji,jj), f_fco2w(ji,jj), f_h2co3(ji,jj), f_hco3(ji,jj), f_co3(ji,jj), f_omarg(ji,jj),  &  ! outputs
                   f_omcal(ji,jj), f_BetaD(ji,jj), f_rhosw(ji,jj), f_opres(ji,jj), f_insitut(ji,jj),            &  ! outputs
@@ -949,5 +828,5 @@
                   !!
                   !! carbon dioxide (CO2); Jerry Blackford code (ostensibly OCMIP-2, but not)
-                  CALL trc_co2_medusa( ztmp(ji,jj), zsal(ji,jj), zdic(ji,jj), zalk(ji,jj), 0.0, f_kw660(ji,jj), f_xco2a,  &  ! inputs
+                  CALL trc_co2_medusa( ztmp(ji,jj), zsal(ji,jj), zdic(ji,jj), zalk(ji,jj), 0.0, f_kw660(ji,jj), f_xco2a(ji,jj),  &  ! inputs
                   f_ph(ji,jj), f_pco2w(ji,jj), f_h2co3(ji,jj), f_hco3(ji,jj), f_co3(ji,jj), f_omcal(ji,jj),               &  ! outputs
                   f_omarg(ji,jj), f_co2flux(ji,jj), f_TDIC(ji,jj), f_TALK(ji,jj), f_dcf(ji,jj), f_henry(ji,jj), iters(ji,jj) )      ! outputs
@@ -968,5 +847,5 @@
                   !!                 quasi-sensible alternatives
                   f_ph(ji,jj)           = 8.1
-                  f_pco2w(ji,jj)        = f_xco2a
+                  f_pco2w(ji,jj)        = f_xco2a(ji,jj)
                   f_h2co3(ji,jj)        = 0.005 * zdic(ji,jj)
                   f_hco3(ji,jj)         = 0.865 * zdic(ji,jj)
@@ -980,15 +859,15 @@
                   f_henry(ji,jj)        = 1.
                   !! AXY (23/06/15): add in some extra MOCSY diagnostics
-                  f_fco2w(ji,jj)        = f_xco2a
+                  f_fco2w(ji,jj)        = f_xco2a(ji,jj)
                   f_BetaD(ji,jj)        = 1.
                   f_rhosw(ji,jj)        = 1.026
                   f_opres(ji,jj)        = 0.
                   f_insitut(ji,jj)      = ztmp(ji,jj)
-                  f_pco2atm(ji,jj)      = f_xco2a
-                  f_fco2atm(ji,jj)      = f_xco2a
+                  f_pco2atm(ji,jj)      = f_xco2a(ji,jj)
+                  f_fco2atm(ji,jj)      = f_xco2a(ji,jj)
                   f_schmidtco2(ji,jj)   = 660.
                   f_kwco2(ji,jj)        = 0.
                   f_K0(ji,jj)           = 0.
-                  f_co2starair(ji,jj)   = f_xco2a
+                  f_co2starair(ji,jj)   = f_xco2a(ji,jj)
                   f_dpco2(ji,jj)        = 0.
 #  endif
@@ -1110,5 +989,5 @@
                       !! AXY (24/11/16): add in extra MOCSY diagnostics
                       IF( med_diag%ATM_XCO2%dgsave ) THEN
-                         f_xco2a_2d(ji,jj) = f_xco2a
+                         f_xco2a_2d(ji,jj) = f_xco2a(ji,jj)
                       ENDIF
                       IF( med_diag%OCN_FCO2%dgsave ) THEN
