Index: /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90
===================================================================
--- /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 7958)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 7958)
@@ -0,0 +1,530 @@
+MODULE air_sea_mod
+   !!======================================================================
+   !!                         ***  MODULE air_sea_mod  ***
+   !! Calculate the carbon chemistry for the whole ocean
+   !!======================================================================
+   !! History :
+   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!----------------------------------------------------------------------
+#if defined key_medusa
+   !!----------------------------------------------------------------------
+   !!                                                   MEDUSA bio-model
+   !!----------------------------------------------------------------------
+
+   IMPLICIT NONE
+   PRIVATE
+      
+   PUBLIC   air_sea        ! 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 air_sea( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE air_sea  ***
+      !! This called from TRC_BIO_MEDUSA and 
+      !!  - calculate air-sea gas exchange
+      !!  - river inputs
+      !!----------------------------------------------------------------------
+      USE bio_medusa_mod,    ONLY: dms_andr, dms_andr2d, dms_aran,        &
+                                   dms_aran2d, dms_hall, dms_hall2d,      &
+                                   dms_simo, dms_simo2d, dms_surf,        &
+                                   dms_surf2d, iters, f_BetaD,            &
+                                   f_co2flux, f_co2flux2d, f_co2starair,  &
+                                   f_co2starair_2d, f_co3,                &
+                                   f_dcf, f_dpco2, f_fco2a_2d, f_fco2atm, &
+                                   f_fco2w, f_fco2w_2d, f_h2co3,          &
+                                   f_hco3, f_henry, f_insitut, f_K0,      &
+                                   f_kw660, f_kw6602d, f_kwco2, f_kwo2,   &
+                                   f_o2flux, f_o2flux2d,                  &
+                                   f_o2sat, f_o2sat2d, f_ocndpco2_2d,     &
+                                   f_ocnk0_2d, f_ocnkwco2_2d,             &
+                                   f_ocnrhosw_2d, f_ocnschco2_2d,         &
+                                   f_omarg, f_omcal, f_opres,             &
+                                   f_pco2a2d, f_pco2atm,                  &
+                                   f_pco2w, f_pco2w2d, f_ph, f_pp0,       &
+                                   f_pp02d, f_rhosw,                      &
+                                   f_riv_alk, f_riv_c, f_riv_n,           &
+                                   f_riv_si, f_runoff,                    & 
+                                   f_schmidtco2, f_TALK, f_TALK2d,        &
+                                   f_TDIC, f_TDIC2d, f_xco2a, f_xco2a_2d, &
+                                   fgco2,                                 &
+                                   zalk, zchd, zchn, zdic, zdin, zoxy,    &
+                                   zpho, zsal, zsil, ztmp
+      USE dom_oce,           ONLY: e3t_0, e3t_n, gphit, tmask
+      USE gastransfer,       ONLY: gas_transfer
+      USE iom,               ONLY: lk_iomput
+      USE in_out_manager,    ONLY: lwp, numout
+      USE mocsy_wrapper,     ONLY: mocsy_interface
+      USE oce,               ONLY: PCO2a_in_cpl
+      USE par_oce,           ONLY: jpim1, jpjm1, jpk
+      USE sbc_oce,           ONLY: fr_i, lk_oasis, qsr, wndm
+      USE sms_medusa,        ONLY: jdms, jdms_input, jdms_model,          &
+                                   jriver_alk, jriver_c,                  &
+                                   jriver_n, jriver_si,                   &
+                                   riv_alk, riv_c, riv_n, riv_si,         &
+                                   zn_dms_chd, zn_dms_chn, zn_dms_din,    &
+                                   zn_dms_mld, zn_dms_qsr 
+      USE trc,               ONLY: med_diag
+      USE trcco2_medusa,     ONLY: trc_co2_medusa
+      USE trcdms_medusa,     ONLY: trc_dms_medusa
+      USE trcoxy_medusa,     ONLY: trc_oxy_medusa
+      USE zdfmxl,            ONLY: hmld
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+
+      !! time (integer timestep)
+      INTEGER, INTENT( in ) ::    kt
+
+      !! jpalm 14-07-2016: convert CO2flux diag from mmol/m2/d to kg/m2/s
+      REAL, PARAMETER :: weight_CO2_mol = 44.0095  !! g / mol
+      REAL, PARAMETER :: secs_in_day    = 86400.0  !! s / d
+      REAL, PARAMETER :: CO2flux_conv   = (1.e-6 * weight_CO2_mol) / secs_in_day
+
+      INTEGER :: ji, jj
+
+# if defined key_roam
+      !!-----------------------------------------------------------
+      !! Air-sea gas exchange
+      !!-----------------------------------------------------------
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            !! OPEN wet point IF..THEN loop
+            if (tmask(ji,jj,1) == 1) then
+               IF (lk_oasis) THEN
+                  !! use 2D atm xCO2 from atm coupling
+                  f_xco2a(ji,jj) = PCO2a_in_cpl(ji,jj)
+               ENDIF
+               !!
+               !! AXY (23/06/15): as part of an effort to update the 
+               !!                 carbonate chemistry in MEDUSA, the gas 
+               !!                 transfer velocity used in the carbon
+               !!                 and oxygen cycles has been harmonised 
+               !!                 and is calculated by the same function 
+               !!                 here; this harmonisation includes
+               !!                 changes to the PML carbonate chemistry 
+               !!                 scheme so that it too makes use of the 
+               !!                 same gas transfer velocity; the
+               !!                 preferred parameterisation of this is 
+               !!                 Wanninkhof (2014), option 7
+               !!
+#   if defined key_debug_medusa
+               IF (lwp) write (numout,*) 'trc_bio_medusa: entering gas_transfer'
+               CALL flush(numout)
+#   endif
+               CALL gas_transfer( wndm(ji,jj), 1, 7,         &  ! inputs
+                                  f_kw660(ji,jj) )              ! outputs
+#   if defined key_debug_medusa
+               IF (lwp) write (numout,*) 'trc_bio_medusa: exiting gas_transfer'
+               CALL flush(numout)
+#   endif
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then
+               !! air pressure (atm); ultimately this will use air 
+               !! pressure at the base of the UKESM1 atmosphere 
+               !!                                     
+               f_pp0(ji,jj)   = 1.0
+               !!
+               !! IF(lwp) WRITE(numout,*) ' MEDUSA ztmp    =', ztmp(ji,jj)
+               !! IF(lwp) WRITE(numout,*) ' MEDUSA wndm    =', wndm(ji,jj)
+               !! IF(lwp) WRITE(numout,*) ' MEDUSA fr_i    =', fr_i(ji,jj)
+               !!
+#  if defined key_axy_carbchem
+#   if defined key_mocsy
+               !!
+               !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate
+               !!                 chemistry package; note that depth is set to
+               !!                 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),     &
+                                    f_pp0(ji,jj),0.0,                        &
+                                    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)
+            ENDIF
+         ENDDO
+      ENDDO
+
+#   else   
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then     
+               iters(ji,jj) = 0
+               !!
+               !! 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(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 (09/01/14): removed iteration and NaN checks; these have
+               !!                 been moved to trc_co2_medusa together with a
+               !!                 fudge that amends erroneous values (this is
+               !!                 intended to be a temporary fudge!); the
+               !!                 output warnings are retained here so that
+               !!                 failure position can be determined
+               if (iters(ji,jj) .eq. 25) then
+                  IF(lwp) WRITE(numout,*) ' trc_bio_medusa: ITERS WARNING, ', &
+                     iters(ji,jj), ' AT (', ji, ', ', jj, ', 1) AT ', kt
+               endif
+            ENDIF
+         ENDDO
+      ENDDO
+
+#   endif
+#  else
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then
+               !! AXY (18/04/13): switch off carbonate chemistry 
+               !!                 calculations; provide quasi-sensible 
+               !!                 alternatives
+               f_ph(ji,jj)           = 8.1
+               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)
+               f_co3(ji,jj)          = 0.130 * zdic(ji,jj)
+               f_omcal(ji,jj) = 4.
+               f_omarg(ji,jj) = 2.
+               f_co2flux(ji,jj)      = 0.
+               f_TDIC(ji,jj)         = zdic(ji,jj)
+               f_TALK(ji,jj)         = zalk(ji,jj)
+               f_dcf(ji,jj)          = 1.026
+               f_henry(ji,jj)        = 1.
+               !! AXY (23/06/15): add in some extra MOCSY diagnostics
+               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(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(ji,jj)
+               f_dpco2(ji,jj)        = 0.
+            ENDIF
+         ENDDO
+      ENDDO
+#  endif
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then
+               !! mmol/m2/s -> mmol/m3/d; correct for sea-ice; divide 
+               !! through by layer thickness
+               f_co2flux(ji,jj) = (1. - fr_i(ji,jj)) * f_co2flux(ji,jj) *    &
+                                  86400. / fse3t(ji,jj,1)
+               !!
+               !! oxygen (O2); OCMIP-2 code
+               !! AXY (23/06/15): amend input list for oxygen to account 
+               !!                 for common gas transfer velocity
+               !! Note that f_kwo2 is an about from the subroutine below,
+               !! which doesn't seem to be used - marc 10/4/17 
+               CALL trc_oxy_medusa(ztmp(ji,jj),zsal(ji,jj),f_kw660(ji,jj),   &
+                                   f_pp0(ji,jj),zoxy(ji,jj),                 &
+                                   f_kwo2(ji,jj),f_o2flux(ji,jj),            &
+                                   f_o2sat(ji,jj))
+               !!
+               !! mmol/m2/s -> mol/m3/d; correct for sea-ice; divide 
+               !! through by layer thickness
+               f_o2flux(ji,jj)  = (1. - fr_i(ji,jj)) * f_o2flux(ji,jj) *     &
+                                  86400. / fse3t(ji,jj,1)
+            ENDIF
+         ENDDO
+      ENDDO
+
+      !! Jpalm (08-2014)
+      !! DMS surface concentration calculation
+      !! initialy added for UKESM1 model.
+      !! using MET-OFFICE subroutine.
+      !! DMS module only needs Chl concentration and MLD
+      !! to get an aproximate value of DMS concentration.
+      !! air-sea fluxes are calculated by atmospheric chemitry model
+      !! from atm and oc-surface concentrations.
+      !!
+      !! AXY (13/03/15): this is amended to calculate all of the DMS
+      !!                 estimates examined during UKESM1 (see
+      !!                 comments in trcdms_medusa.F90)
+      !!
+      IF (jdms == 1) THEN
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               if (tmask(ji,jj,1) == 1) then
+                  !!
+                  !! feed in correct inputs
+                  if (jdms_input .eq. 0) then
+                     !! use instantaneous inputs
+                     CALL trc_dms_medusa(zchn(ji,jj),zchd(ji,jj),             &
+                                         hmld(ji,jj),qsr(ji,jj),              &
+                                         zdin(ji,jj),                         &
+                                         dms_andr(ji,jj),dms_simo(ji,jj),     &
+                                         dms_aran(ji,jj),dms_hall(ji,jj))
+                  else
+                     !! use diel-average inputs
+                     CALL trc_dms_medusa(zn_dms_chn(ji,jj),zn_dms_chd(ji,jj), &
+                                         zn_dms_mld(ji,jj),zn_dms_qsr(ji,jj), &
+                                         zn_dms_din(ji,jj),                   &
+                                         dms_andr(ji,jj),dms_simo(ji,jj),     &
+                                         dms_aran(ji,jj),dms_hall(ji,jj))
+                  endif
+                  !!
+                  !! assign correct output to variable passed to atmosphere
+                  if (jdms_model .eq. 1) then
+                     dms_surf(ji,jj) = dms_andr(ji,jj)
+                  elseif (jdms_model .eq. 2) then
+                     dms_surf(ji,jj) = dms_simo(ji,jj)
+                  elseif (jdms_model .eq. 3) then
+                     dms_surf(ji,jj) = dms_aran(ji,jj)
+                  elseif (jdms_model .eq. 4) then
+                     dms_surf(ji,jj) = dms_hall(ji,jj)
+                  endif
+                  !!
+                  !! 2D diag through iom_use
+                  IF( lk_iomput ) THEN
+                     IF( med_diag%DMS_SURF%dgsave ) THEN
+                        dms_surf2d(ji,jj) = dms_surf(ji,jj)
+                     ENDIF
+                     IF( med_diag%DMS_ANDR%dgsave ) THEN
+                        dms_andr2d(ji,jj) = dms_andr(ji,jj)
+                     ENDIF
+                     IF( med_diag%DMS_SIMO%dgsave ) THEN
+                        dms_simo2d(ji,jj) = dms_simo(ji,jj)
+                     ENDIF
+                     IF( med_diag%DMS_ARAN%dgsave ) THEN
+                        dms_aran2d(ji,jj) = dms_aran(ji,jj)
+                     ENDIF
+                     IF( med_diag%DMS_HALL%dgsave ) THEN
+                        dms_hall2d(ji,jj) = dms_hall(ji,jj)
+                     ENDIF
+#   if defined key_debug_medusa
+                     IF (lwp) write (numout,*)                                &
+                        'trc_bio_medusa: finish calculating dms'
+                     CALL flush(numout)
+#   endif 
+                  ENDIF      !! End iom
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF                  !! End IF (jdms == 1)
+
+
+      !!
+      !! store 2D outputs
+      !!
+      !! JPALM -- 17-11-16 -- put fgco2 out of diag request
+      !!       is needed for coupling; pass through restart
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,1) == 1) then
+               !! IF( med_diag%FGCO2%dgsave ) THEN
+               !! convert from  mol/m2/day to kg/m2/s
+               !! mmol-C/m3/d -> kg-CO2/m2/s
+               fgco2(ji,jj) = f_co2flux(ji,jj) * fse3t(ji,jj,1) *             &
+                              CO2flux_conv
+               !! ENDIF
+               IF ( lk_iomput ) THEN
+                  IF( med_diag%ATM_PCO2%dgsave ) THEN
+                     f_pco2a2d(ji,jj) = f_pco2atm(ji,jj)
+                  ENDIF
+                  IF( med_diag%OCN_PCO2%dgsave ) THEN
+                     f_pco2w2d(ji,jj) = f_pco2w(ji,jj)
+                  ENDIF
+                  IF( med_diag%CO2FLUX%dgsave ) THEN
+                     !! mmol/m3/d -> mmol/m2/d
+                     f_co2flux2d(ji,jj) = f_co2flux(ji,jj) * fse3t(ji,jj,1)
+                  ENDIF
+                  IF( med_diag%TCO2%dgsave ) THEN
+                     f_TDIC2d(ji,jj) = f_TDIC(ji,jj)
+                  ENDIF
+                  IF( med_diag%TALK%dgsave ) THEN
+                     f_TALK2d(ji,jj) = f_TALK(ji,jj)
+                  ENDIF
+                  IF( med_diag%KW660%dgsave ) THEN
+                      f_kw6602d(ji,jj) = f_kw660(ji,jj)
+                  ENDIF
+                  IF( med_diag%ATM_PP0%dgsave ) THEN
+                      f_pp02d(ji,jj) = f_pp0(ji,jj)
+                  ENDIF
+                  IF( med_diag%O2FLUX%dgsave ) THEN
+                     f_o2flux2d(ji,jj) = f_o2flux(ji,jj)
+                  ENDIF
+                  IF( med_diag%O2SAT%dgsave ) THEN
+                     f_o2sat2d(ji,jj) = f_o2sat(ji,jj)
+                  ENDIF
+                  !! AXY (24/11/16): add in extra MOCSY diagnostics
+                  IF( med_diag%ATM_XCO2%dgsave ) THEN
+                     f_xco2a_2d(ji,jj) = f_xco2a(ji,jj)
+                  ENDIF
+                  IF( med_diag%OCN_FCO2%dgsave ) THEN
+                     f_fco2w_2d(ji,jj) = f_fco2w(ji,jj)
+                  ENDIF
+                  IF( med_diag%ATM_FCO2%dgsave ) THEN
+                     f_fco2a_2d(ji,jj) = f_fco2atm(ji,jj)
+                  ENDIF
+                  IF( med_diag%OCN_RHOSW%dgsave ) THEN
+                     f_ocnrhosw_2d(ji,jj) = f_rhosw(ji,jj)
+                  ENDIF
+                  IF( med_diag%OCN_SCHCO2%dgsave ) THEN
+                     f_ocnschco2_2d(ji,jj) = f_schmidtco2(ji,jj)
+                  ENDIF
+                  IF( med_diag%OCN_KWCO2%dgsave ) THEN
+                     f_ocnkwco2_2d(ji,jj) = f_kwco2(ji,jj)
+                  ENDIF
+                  IF( med_diag%OCN_K0%dgsave ) THEN
+                     f_ocnk0_2d(ji,jj) = f_K0(ji,jj)
+                  ENDIF
+                  IF( med_diag%CO2STARAIR%dgsave ) THEN
+                     f_co2starair_2d(ji,jj) = f_co2starair(ji,jj)
+                  ENDIF
+                  IF( med_diag%OCN_DPCO2%dgsave ) THEN
+                     f_ocndpco2_2d(ji,jj) = f_dpco2(ji,jj)
+                  ENDIF
+               ENDIF
+            ENDIF
+         ENDDO
+      ENDDO
+
+# endif
+
+      !!-----------------------------------------------------------------
+      !! River inputs
+      !!-----------------------------------------------------------------
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            !! OPEN wet point IF..THEN loop
+            if (tmask(ji,jj,1) == 1) then
+               !!
+               !! runoff comes in as        kg / m2 / s
+               !! used and written out as   m3 / m2 / d (= m / d)
+               !! where                     1000 kg / m2 / d = 
+               !!                             1 m3 / m2 / d = 1 m / d
+               !!
+               !! AXY (17/07/14): the compiler doesn't like this line for 
+               !!                 some reason; as MEDUSA doesn't even use 
+               !!                 runoff for riverine inputs, a temporary 
+               !!                 solution is to switch off runoff entirely
+               !!                 here; again, this change is one of several 
+               !!                 that will need revisiting once MEDUSA has 
+               !!                 bedded down in UKESM1; particularly so if 
+               !!                 the land scheme provides information
+               !!                 concerning nutrient fluxes
+               !!
+               !! f_runoff(ji,jj) = sf_rnf(1)%fnow(ji,jj,1) / 1000. * 60. *   &
+               !!                   60. * 24.
+               f_runoff(ji,jj) = 0.0
+               !!
+               !! nutrients are added via rivers to the model in one of 
+               !! two ways:
+               !!   1. via river concentration; i.e. the average nutrient 
+               !!      concentration of a river water is described by a 
+               !!      spatial file, and this is multiplied by runoff to 
+               !!      give a nutrient flux
+               !!   2. via direct river flux; i.e. the average nutrient 
+               !!      flux due to rivers is described by a spatial file, 
+               !!      and this is simply applied as a direct nutrient 
+               !!      flux (i.e. it does not relate or respond to model 
+               !!      runoff) nutrient fields are derived from the 
+               !!      GlobalNEWS 2 database; carbon and alkalinity are 
+               !!      derived from continent-scale DIC estimates (Huang et 
+               !!      al., 2012) and some Arctic river alkalinity 
+               !!      estimates (Katya?)
+               !! 
+               !! as of 19/07/12, riverine nutrients can now be spread 
+               !! vertically across several grid cells rather than just 
+               !! poured into the surface box; this block of code is still 
+               !! executed, however, to set up the total amounts of 
+               !! nutrient entering via rivers
+               !!
+               !! nitrogen
+               if (jriver_n .eq. 1) then
+                  !! river concentration specified; use runoff to 
+                  !! calculate input
+                  f_riv_n(ji,jj) = f_runoff(ji,jj) * riv_n(ji,jj)
+               elseif (jriver_n .eq. 2) then
+                  !! river flux specified; independent of runoff
+                  f_riv_n(ji,jj) = riv_n(ji,jj)
+               endif
+               !!
+               !! silicon
+               if (jriver_si .eq. 1) then
+                  !! river concentration specified; use runoff to 
+                  !! calculate input
+                  f_riv_si(ji,jj) = f_runoff(ji,jj) * riv_si(ji,jj)
+               elseif (jriver_si .eq. 2) then
+                  !! river flux specified; independent of runoff
+                  f_riv_si(ji,jj) = riv_si(ji,jj)
+               endif
+               !!
+               !! carbon
+               if (jriver_c .eq. 1) then
+                  !! river concentration specified; use runoff to 
+                  !! calculate input
+                  f_riv_c(ji,jj) = f_runoff(ji,jj) * riv_c(ji,jj)
+               elseif (jriver_c .eq. 2) then
+                  !! river flux specified; independent of runoff
+                  f_riv_c(ji,jj) = riv_c(ji,jj)
+               endif
+               !!
+               !! alkalinity
+               if (jriver_alk .eq. 1) then
+                  !! river concentration specified; use runoff to 
+                  !! calculate input
+                  f_riv_alk(ji,jj) = f_runoff(ji,jj) * riv_alk(ji,jj)
+               elseif (jriver_alk .eq. 2) then
+                  !! river flux specified; independent of runoff
+                  f_riv_alk(ji,jj) = riv_alk(ji,jj)
+               endif
+            ENDIF
+         ENDDO
+      ENDDO
+
+   END SUBROUTINE air_sea
+
+#else
+   !!======================================================================
+   !!  Dummy module :                                   No MEDUSA bio-model
+   !!======================================================================
+CONTAINS
+   SUBROUTINE air_sea( )                    ! Empty routine
+      WRITE(*,*) 'air_sea: You should not have seen this print! error?'
+   END SUBROUTINE air_sea
+#endif 
+
+   !!======================================================================
+END MODULE air_sea_mod
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 7957)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90	(revision 7958)
@@ -94,4 +94,9 @@
    !! Oxygen production and consumption (and non-consumption)
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: foxy_prod,foxy_cons,foxy_anox
+
+   !! Add DMS in MEDUSA for UKESM1 model
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dms_surf
+   !! AXY (13/03/15): add in other DMS calculations
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dms_andr,dms_simo,dms_aran,dms_hall
 
    !! Benthic fluxes
@@ -251,4 +256,7 @@
                foxy_prod(jpi,jpj), foxy_cons(jpi,jpj),                &
                foxy_anox(jpi,jpj),                                    &
+               dms_surf(jpi,jpj),                                     &
+               dms_andr(jpi,jpj),dms_simo(jpi,jpj),                   &
+               dms_aran(jpi,jpj),dms_hall(jpi,jpj),                   &
                f_sbenin_n(jpi,jpj),f_sbenin_fe(jpi,jpj),              &
                f_sbenin_c(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 7957)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90	(revision 7958)
@@ -31,7 +31,8 @@
       !!  - ...
       !!----------------------------------------------------------------------
-      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,     &
+      USE bio_medusa_mod,    ONLY: iters, f_BetaD, f_co2flux,             &
+                                   f_co2starair, f_co3, f_dcf, f_dpco2,   &
+                                   f_fco2atm, f_fco2w, f_h2co3, f_hco3,   &
+                                   f_henry,                               &
                                    f_insitut, f_K0, f_kw660, f_kwco2,     &
                                    f_omarg, f_omcal, f_opres, f_pco2atm,  &
@@ -82,10 +83,11 @@
       i2_omarg(:,:) = 0
 
-      !! loop over 3D space
+      !! Loop over levels
       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 (tmask(ji,jj,jk).eq.1) THEN
                   IF (lk_oasis) THEN
                      !! use 2D atm xCO2 from atm coupling
@@ -116,5 +118,5 @@
                         ' at (', ji, ',', jj, ',', jk, ') at time', kt
 	             IF(lwp) WRITE(numout,*)                                 &
-                        ' carb_chem: T SWITCHING 3D, ',                 &
+                        ' 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
@@ -125,11 +127,16 @@
                         ' 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
-                  !!
+               ENDIF
+            ENDDO
+         ENDDO
+
+         !! Blank input variables not used at this stage (they 
+         !! relate to air-sea flux)
+         f_kw660(:,:) = 1.0
+         f_pp0(:,:)   = 1.0
+
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk).eq.1) THEN
                   !! calculate carbonate chemistry at grid cell midpoint
 #  if defined key_mocsy
@@ -168,6 +175,6 @@
                                       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))
+                                      f_TALK(ji,jj),f_dcf(ji,jj),           &
+                                      f_henry(ji,jj),iters(ji,jj))
                   !! 
                   !! AXY (28/02/14): check output fields
@@ -229,5 +236,5 @@
                      i2_omarg(ji,jj)   = 1
                   endif
-               endif
+               ENDIF
             ENDDO
          ENDDO
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 7957)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 7958)
@@ -99,4 +99,5 @@
       USE bio_medusa_init_mod,        ONLY: bio_medusa_init
       USE carb_chem_mod,              ONLY: carb_chem
+      USE air_sea_mod,                ONLY: air_sea
       USE bio_medusa_diag_slice_mod,  ONLY: bio_medusa_diag_slice
       USE bio_medusa_fin_mod,         ONLY: bio_medusa_fin
@@ -304,5 +305,6 @@
       !!
       !! flags to help with calculating the position of the CCD
-      INTEGER, DIMENSION(jpi,jpj) ::     i2_omcal,i2_omarg
+! Moved into carb_chem.F90 - marc 20/4/17
+!      INTEGER, DIMENSION(jpi,jpj) ::     i2_omcal,i2_omarg
       !!
       !! AXY (24/11/16): add xCO2 variable for atmosphere (what we actually have)
@@ -339,7 +341,7 @@
       !! Jpalm (11-08-2014)
       !! add DMS in MEDUSA for UKESM1 model
-      REAL(wp), DIMENSION(jpi,jpj) ::    dms_surf
+!      REAL(wp), DIMENSION(jpi,jpj) ::    dms_surf
       !! AXY (13/03/15): add in other DMS calculations
-      REAL(wp), DIMENSION(jpi,jpj) ::    dms_andr, dms_simo, dms_aran, dms_hall
+!      REAL(wp), DIMENSION(jpi,jpj) ::    dms_andr, dms_simo, dms_aran, dms_hall
 
 # endif
@@ -578,20 +580,23 @@
          DO ji = 2,jpim1
             !! OPEN wet point IF..THEN loop
-            if (tmask(ji,jj,jk).eq.1) then               
-               !!======================================================================
+            if (tmask(ji,jj,jk) == 1) then               
+               !!===========================================================
                !! SETUP LOCAL GRID CELL
-               !!======================================================================
-               !!
-               !!---------------------------------------------------------------------
+               !!===========================================================
+               !!
+               !!-----------------------------------------------------------
                !! Some notes on grid vertical structure
-               !! - fsdepw(ji,jj,jk) is the depth of the upper surface of level jk
-               !! - fsde3w(ji,jj,jk) is *approximately* the midpoint of level jk
+               !! - fsdepw(ji,jj,jk) is the depth of the upper surface of 
+               !!   level jk
+               !! - fsde3w(ji,jj,jk) is *approximately* the midpoint of 
+               !!   level jk
                !! - fse3t(ji,jj,jk)  is the thickness of level jk
-               !!---------------------------------------------------------------------
+               !!-----------------------------------------------------------
                !!
                !! AXY (01/03/10): set up level depth (bottom of level)
                fdep1(ji,jj) = fsdepw(ji,jj,jk) + fse3t(ji,jj,jk)
                !! AXY (28/11/16): local seafloor depth
-               !!                 previously mbathy(ji,jj) - 1, now mbathy(ji,jj)
+               !!                 previously mbathy(ji,jj) - 1, now 
+               !!                 mbathy(ji,jj)
                mbathy(ji,jj) = mbathy(ji,jj)
                !!
@@ -599,10 +604,16 @@
                !! negative values of state variables are not allowed to
                !! contribute to the calculated fluxes
-               zchn(ji,jj) = max(0.,trn(ji,jj,jk,jpchn)) !! non-diatom chlorophyll
-               zchd(ji,jj) = max(0.,trn(ji,jj,jk,jpchd)) !! diatom chlorophyll
-               zphn(ji,jj) = max(0.,trn(ji,jj,jk,jpphn)) !! non-diatoms
-               zphd(ji,jj) = max(0.,trn(ji,jj,jk,jpphd)) !! diatoms
-               zpds(ji,jj) = max(0.,trn(ji,jj,jk,jppds)) !! diatom silicon
-               !! AXY (28/01/10): probably need to take account of chl/biomass connection
+               !! non-diatom chlorophyll
+               zchn(ji,jj) = max(0.,trn(ji,jj,jk,jpchn))
+               !! diatom chlorophyll
+               zchd(ji,jj) = max(0.,trn(ji,jj,jk,jpchd))
+               !! non-diatoms
+               zphn(ji,jj) = max(0.,trn(ji,jj,jk,jpphn))
+               !! diatoms
+               zphd(ji,jj) = max(0.,trn(ji,jj,jk,jpphd))
+               !! diatom silicon
+               zpds(ji,jj) = max(0.,trn(ji,jj,jk,jppds))
+               !! AXY (28/01/10): probably need to take account of 
+               !! chl/biomass connection
                if (zchn(ji,jj).eq.0.) zphn(ji,jj) = 0.
                if (zchd(ji,jj).eq.0.) zphd(ji,jj) = 0.
@@ -612,17 +623,28 @@
 	       if (zpds(ji,jj).eq.0.) zphd(ji,jj) = 0.
 	       if (zphd(ji,jj).eq.0.) zpds(ji,jj) = 0.
-               zzmi(ji,jj) = max(0.,trn(ji,jj,jk,jpzmi)) !! microzooplankton
-               zzme(ji,jj) = max(0.,trn(ji,jj,jk,jpzme)) !! mesozooplankton
-               zdet(ji,jj) = max(0.,trn(ji,jj,jk,jpdet)) !! detrital nitrogen
-               zdin(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) !! dissolved inorganic nitrogen
-               zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil)) !! dissolved silicic acid
-               zfer(ji,jj) = max(0.,trn(ji,jj,jk,jpfer)) !! dissolved "iron"
+               !! microzooplankton
+               zzmi(ji,jj) = max(0.,trn(ji,jj,jk,jpzmi))
+               !! mesozooplankton
+               zzme(ji,jj) = max(0.,trn(ji,jj,jk,jpzme))
+               !! detrital nitrogen
+               zdet(ji,jj) = max(0.,trn(ji,jj,jk,jpdet))
+               !! dissolved inorganic nitrogen
+               zdin(ji,jj) = max(0.,trn(ji,jj,jk,jpdin))
+               !! dissolved silicic acid
+               zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil))
+               !! dissolved "iron"
+               zfer(ji,jj) = max(0.,trn(ji,jj,jk,jpfer))
 # if defined key_roam
-               zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc)) !! detrital carbon
-               zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic)) !! dissolved inorganic carbon
-               zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk)) !! alkalinity
-               zoxy(ji,jj) = max(0.,trn(ji,jj,jk,jpoxy)) !! oxygen
+               !! detrital carbon
+               zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc))
+               !! dissolved inorganic carbon
+               zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic))
+               !! alkalinity
+               zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk))
+               !! oxygen
+               zoxy(ji,jj) = max(0.,trn(ji,jj,jk,jpoxy))
 #  if defined key_axy_carbchem && defined key_mocsy
-               zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 !! phosphate via DIN and Redfield
+               !! phosphate via DIN and Redfield
+               zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0
 #  endif
                !!
@@ -633,47 +655,59 @@
 	       !! 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 2D, ', &
-                  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 2D, ', &
-                  tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem)
-                  ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem) !! temperature
+                  IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T WARNING 2D, ',  &
+                     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 2D, ',&
+                     tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem)
+                  !! temperatur
+                  ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem)
                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 2D, ', &
-                  tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (',    &
-                  ji, ',', jj, ',', jk, ') at time', kt
+                     tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (',    &
+                     ji, ',', jj, ',', jk, ') at time', kt
                endif
 # else
-               zdtc(ji,jj) = zdet(ji,jj) * xthetad              !! implicit detrital carbon
+               !! implicit detrital carbon
+               zdtc(ji,jj) = zdet(ji,jj) * xthetad
 # endif
 # if defined key_debug_medusa
                if (idf.eq.1) then
-               !! AXY (15/01/10)
+                  !! AXY (15/01/10)
                   if (trn(ji,jj,jk,jpdin).lt.0.) then
                      IF (lwp) write (numout,*) '------------------------------'
-                     IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR =', trn(ji,jj,jk,jpdin)
-                     IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR @', ji, jj, jk, kt
+                     IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR =',        &
+                        trn(ji,jj,jk,jpdin)
+                     IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR @',        &
+                        ji, jj, jk, kt
                   endif
                   if (trn(ji,jj,jk,jpsil).lt.0.) then
                      IF (lwp) write (numout,*) '------------------------------'
-                     IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR =', trn(ji,jj,jk,jpsil)
-                     IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR @', ji, jj, jk, kt
+                     IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR =',        &
+                        trn(ji,jj,jk,jpsil)
+                     IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR @',        &
+                        ji, jj, jk, kt
                   endif
 #  if defined key_roam
                   if (trn(ji,jj,jk,jpdic).lt.0.) then
                      IF (lwp) write (numout,*) '------------------------------'
-                     IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR =', trn(ji,jj,jk,jpdic)
-                     IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR @', ji, jj, jk, kt
+                     IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR =',        &
+                        trn(ji,jj,jk,jpdic)
+                     IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR @',        &
+                        ji, jj, jk, kt
                   endif
                   if (trn(ji,jj,jk,jpalk).lt.0.) then
                      IF (lwp) write (numout,*) '------------------------------'
-                     IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR =', trn(ji,jj,jk,jpalk)
-                     IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR @', ji, jj, jk, kt
+                     IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR =',        &
+                        trn(ji,jj,jk,jpalk)
+                     IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR @',        &
+                        ji, jj, jk, kt
                   endif
                   if (trn(ji,jj,jk,jpoxy).lt.0.) then
                      IF (lwp) write (numout,*) '------------------------------'
-                     IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR =', trn(ji,jj,jk,jpoxy)
-                     IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR @', ji, jj, jk, kt
+                     IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR =',        &
+                        trn(ji,jj,jk,jpoxy)
+                     IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR @',        &
+                        ji, jj, jk, kt
                   endif
 #  endif
@@ -713,41 +747,54 @@
                IF( lk_iomput ) THEN
                   IF ( med_diag%INVTN%dgsave )   THEN
-                     ftot_n(ji,jj)  = ftot_n(ji,jj) + &
-                             (fse3t(ji,jj,jk) * ( zphn(ji,jj) + zphd(ji,jj) + zzmi(ji,jj) + zzme(ji,jj) + zdet(ji,jj) + zdin(ji,jj) ) )
+                     ftot_n(ji,jj)  = ftot_n(ji,jj) +                        &
+                        (fse3t(ji,jj,jk) * ( zphn(ji,jj) + zphd(ji,jj) +     &
+                                             zzmi(ji,jj) + zzme(ji,jj) +     &
+                                             zdet(ji,jj) + zdin(ji,jj) ) )
                   ENDIF
                   IF ( med_diag%INVTSI%dgsave )  THEN
-                     ftot_si(ji,jj) = ftot_si(ji,jj) + & 
-                             (fse3t(ji,jj,jk) * ( zpds(ji,jj) + zsil(ji,jj) ) )
+                     ftot_si(ji,jj) = ftot_si(ji,jj) +                       & 
+                        (fse3t(ji,jj,jk) * ( zpds(ji,jj) + zsil(ji,jj) ) )
                   ENDIF
                   IF ( med_diag%INVTFE%dgsave )  THEN
-                     ftot_fe(ji,jj) = ftot_fe(ji,jj) + & 
-                             (fse3t(ji,jj,jk) * ( xrfn * ( zphn(ji,jj) + zphd(ji,jj) + zzmi(ji,jj) + zzme(ji,jj) + zdet(ji,jj) ) + zfer(ji,jj) ) )
+                     ftot_fe(ji,jj) = ftot_fe(ji,jj) +                       & 
+                        (fse3t(ji,jj,jk) * ( xrfn *                          &
+                                            ( zphn(ji,jj) + zphd(ji,jj) +    &
+                                              zzmi(ji,jj) + zzme(ji,jj) +    &
+                                              zdet(ji,jj) ) + zfer(ji,jj) ) )
                   ENDIF
 # if defined key_roam
                   IF ( med_diag%INVTC%dgsave )  THEN
                      ftot_c(ji,jj)  = ftot_c(ji,jj) + & 
-                             (fse3t(ji,jj,jk) * ( (xthetapn * zphn(ji,jj)) + (xthetapd * zphd(ji,jj)) + &
-                             (xthetazmi * zzmi(ji,jj)) + (xthetazme * zzme(ji,jj)) + zdtc(ji,jj) +   &
-                             zdic(ji,jj) ) )
+                        (fse3t(ji,jj,jk) * ( (xthetapn * zphn(ji,jj)) +      &
+                                             (xthetapd * zphd(ji,jj)) +      &
+                                             (xthetazmi * zzmi(ji,jj)) +     &
+                                             (xthetazme * zzme(ji,jj)) +     &
+                                             zdtc(ji,jj) + zdic(ji,jj) ) )
                   ENDIF
                   IF ( med_diag%INVTALK%dgsave ) THEN
-                     ftot_a(ji,jj)  = ftot_a(ji,jj) + (fse3t(ji,jj,jk) * ( zalk(ji,jj) ) )
+                     ftot_a(ji,jj)  = ftot_a(ji,jj) + (fse3t(ji,jj,jk) *     &
+                                                       zalk(ji,jj))
                   ENDIF
                   IF ( med_diag%INVTO2%dgsave )  THEN
-                     ftot_o2(ji,jj) = ftot_o2(ji,jj) + (fse3t(ji,jj,jk) * ( zoxy(ji,jj) ) )
+                     ftot_o2(ji,jj) = ftot_o2(ji,jj) + (fse3t(ji,jj,jk) *    &
+                                                        zoxy(ji,jj))
                   ENDIF
                   !!
                   !! AXY (10/11/16): CMIP6 diagnostics
                   IF ( med_diag%INTDISSIC%dgsave ) THEN
-                     intdissic(ji,jj) = intdissic(ji,jj) + (fse3t(ji,jj,jk) * zdic(ji,jj))
+                     intdissic(ji,jj) = intdissic(ji,jj) +                   &
+                                        (fse3t(ji,jj,jk) * zdic(ji,jj))
                   ENDIF
                   IF ( med_diag%INTDISSIN%dgsave ) THEN
-                     intdissin(ji,jj) = intdissin(ji,jj) + (fse3t(ji,jj,jk) * zdin(ji,jj))
+                     intdissin(ji,jj) = intdissin(ji,jj) +                   &
+                                        (fse3t(ji,jj,jk) * zdin(ji,jj))
                   ENDIF
                   IF ( med_diag%INTDISSISI%dgsave ) THEN
-                     intdissisi(ji,jj) = intdissisi(ji,jj) + (fse3t(ji,jj,jk) * zsil(ji,jj))
+                     intdissisi(ji,jj) = intdissisi(ji,jj) +                 &
+                                         (fse3t(ji,jj,jk) * zsil(ji,jj))
                   ENDIF
                   IF ( med_diag%INTTALK%dgsave ) THEN
-                     inttalk(ji,jj) = inttalk(ji,jj) + (fse3t(ji,jj,jk) * zalk(ji,jj))
+                     inttalk(ji,jj) = inttalk(ji,jj) +                       &
+                                      (fse3t(ji,jj,jk) * zalk(ji,jj))
                   ENDIF
                   IF ( med_diag%O2min%dgsave ) THEN
@@ -755,5 +802,7 @@
                         o2min(ji,jj)  = zoxy(ji,jj)
                         IF ( med_diag%ZO2min%dgsave ) THEN
-                           zo2min(ji,jj) = (fsdepw(ji,jj,jk) + fdep1(ji,jj)) / 2. !! layer midpoint
+                           !! layer midpoint
+                           zo2min(ji,jj) = (fsdepw(ji,jj,jk) +               &
+                                            fdep1(ji,jj)) / 2.0
                         ENDIF
                      endif
@@ -764,259 +813,24 @@
                CALL flush(numout)
 
-               !!======================================================================
-               !! LOCAL GRID CELL CALCULATIONS
-               !!======================================================================
-               !!
+            ENDIF
+         ENDDO
+         ENDDO
+
+         !!----------------------------------------------------------------
+         !! Calculate air-sea gas exchange and river inputs
+         !!----------------------------------------------------------------
+         IF ( jk == 1 ) THEN
+            call air_sea( kt )
+         END IF
+
 # if defined key_roam
-               if ( jk .eq. 1 ) then
-                  !!----------------------------------------------------------------------
-                  !! Air-sea gas exchange
-                  !!----------------------------------------------------------------------
-                  IF (lk_oasis) THEN
-                     f_xco2a(ji,jj) = PCO2a_in_cpl(ji,jj)        !! use 2D atm xCO2 from atm coupling
-                  ENDIF
-                  !!
-                  !! AXY (23/06/15): as part of an effort to update the carbonate chemistry
-                  !!                 in MEDUSA, the gas transfer velocity used in the carbon
-                  !!                 and oxygen cycles has been harmonised and is calculated
-                  !!                 by the same function here; this harmonisation includes
-                  !!                 changes to the PML carbonate chemistry scheme so that
-                  !!                 it too makes use of the same gas transfer velocity; the
-                  !!                 preferred parameterisation of this is Wanninkhof (2014),
-                  !!                 option 7
-                  !!
-#   if defined key_debug_medusa
-                     IF (lwp) write (numout,*) 'trc_bio_medusa: entering gas_transfer'
-                     CALL flush(numout)
-#   endif
-                  CALL gas_transfer( wndm(ji,jj), 1, 7, &  ! inputs
-                                     f_kw660(ji,jj) )        ! outputs
-#   if defined key_debug_medusa
-                     IF (lwp) write (numout,*) 'trc_bio_medusa: exiting gas_transfer'
-                     CALL flush(numout)
-#   endif
-                  !!
-                  !! air pressure (atm); ultimately this will use air pressure at the base
-                  !! of the UKESM1 atmosphere 
-                  !!                                     
-                  f_pp0(ji,jj)   = 1.0
-                  !!
-                  !! IF(lwp) WRITE(numout,*) ' MEDUSA ztmp    =', ztmp(ji,jj)
-                  !! IF(lwp) WRITE(numout,*) ' MEDUSA zwind_i =', zwind_i(ji,jj)
-                  !! IF(lwp) WRITE(numout,*) ' MEDUSA zwind_j =', zwind_j(ji,jj)
-                  !! IF(lwp) WRITE(numout,*) ' MEDUSA wndm    =', wndm(ji,jj)
-                  !! IF(lwp) WRITE(numout,*) ' MEDUSA fr_i    =', fr_i(ji,jj)
-                  !!
-#  if defined key_axy_carbchem
-#   if defined key_mocsy
-                  !!
-                  !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate
-                  !!                 chemistry package; note that depth is set to
-                  !!                 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(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
-                  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                  
-                  iters(ji,jj) = 0
-                  !!
-                  !! 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(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
-                  !!
-                  !! AXY (09/01/14): removed iteration and NaN checks; these have
-                  !!                 been moved to trc_co2_medusa together with a
-                  !!                 fudge that amends erroneous values (this is
-                  !!                 intended to be a temporary fudge!); the
-                  !!                 output warnings are retained here so that
-                  !!                 failure position can be determined
-                  if (iters(ji,jj) .eq. 25) then
-                     IF(lwp) WRITE(numout,*) ' trc_bio_medusa: ITERS WARNING, ', &
-                     iters(ji,jj), ' AT (', ji, ', ', jj, ', ', jk, ') AT ', kt
-                  endif
-#   endif
-#  else
-                  !! AXY (18/04/13): switch off carbonate chemistry calculations; provide
-                  !!                 quasi-sensible alternatives
-                  f_ph(ji,jj)           = 8.1
-                  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)
-                  f_co3(ji,jj)          = 0.130 * zdic(ji,jj)
-                  f_omcal(ji,jj) = 4.
-                  f_omarg(ji,jj) = 2.
-                  f_co2flux(ji,jj)      = 0.
-                  f_TDIC(ji,jj)         = zdic(ji,jj)
-                  f_TALK(ji,jj)         = zalk(ji,jj)
-                  f_dcf(ji,jj)          = 1.026
-                  f_henry(ji,jj)        = 1.
-                  !! AXY (23/06/15): add in some extra MOCSY diagnostics
-                  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(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(ji,jj)
-                  f_dpco2(ji,jj)        = 0.
-#  endif
-                  !!
-                  !! mmol/m2/s -> mmol/m3/d; correct for sea-ice; divide through by layer thickness
-                  f_co2flux(ji,jj) = (1. - fr_i(ji,jj)) * f_co2flux(ji,jj) * 86400. / fse3t(ji,jj,jk)
-                  !!
-                  !! oxygen (O2); OCMIP-2 code
-                  !! AXY (23/06/15): amend input list for oxygen to account for common gas
-                  !!                 transfer velocity
-                  !! Note that f_kwo2 is an about from the subroutine below,
-                  !! which doesn't seem to be used - marc 10/4/17 
-                  CALL trc_oxy_medusa( ztmp(ji,jj), zsal(ji,jj), f_kw660(ji,jj), f_pp0(ji,jj), zoxy(ji,jj),  &  ! inputs
-                  f_kwo2(ji,jj), f_o2flux(ji,jj), f_o2sat(ji,jj) )                                ! outputs
-                  !!
-                  !! mmol/m2/s -> mol/m3/d; correct for sea-ice; divide through by layer thickness
-                  f_o2flux(ji,jj)  = (1. - fr_i(ji,jj)) * f_o2flux(ji,jj) * 86400. / fse3t(ji,jj,jk)
-                  !!
-                  !! Jpalm (08-2014)
-                  !! DMS surface concentration calculation
-                  !! initialy added for UKESM1 model.
-                  !! using MET-OFFICE subroutine.
-                  !! DMS module only needs Chl concentration and MLD
-                  !! to get an aproximate value of DMS concentration.
-                  !! air-sea fluxes are calculated by atmospheric chemitry model
-                  !! from atm and oc-surface concentrations.
-                  !!
-                  !! AXY (13/03/15): this is amended to calculate all of the DMS
-                  !!                 estimates examined during UKESM1 (see comments
-                  !!                 in trcdms_medusa.F90)
-                  !!
-                  IF (jdms .eq. 1) THEN
-                     !!
-                     !! feed in correct inputs
-                     if (jdms_input .eq. 0) then
-                        !! use instantaneous inputs
-                        CALL trc_dms_medusa( zchn(ji,jj), zchd(ji,jj), hmld(ji,jj), qsr(ji,jj), zdin(ji,jj), &  ! inputs
-                        dms_andr(ji,jj), dms_simo(ji,jj), dms_aran(ji,jj), dms_hall(ji,jj) )                           ! outputs
-                     else
-                        !! use diel-average inputs
-                        CALL trc_dms_medusa( zn_dms_chn(ji,jj), zn_dms_chd(ji,jj), &  ! inputs
-                        zn_dms_mld(ji,jj), zn_dms_qsr(ji,jj), zn_dms_din(ji,jj),   &  ! inputs
-                        dms_andr(ji,jj), dms_simo(ji,jj), dms_aran(ji,jj), dms_hall(ji,jj) )                      ! outputs
-                     endif
-                     !!
-                     !! assign correct output to variable passed to atmosphere
-                     if     (jdms_model .eq. 1) then
-                        dms_surf(ji,jj) = dms_andr(ji,jj)
-                     elseif (jdms_model .eq. 2) then
-                        dms_surf(ji,jj) = dms_simo(ji,jj)
-                     elseif (jdms_model .eq. 3) then
-                        dms_surf(ji,jj) = dms_aran(ji,jj)
-                     elseif (jdms_model .eq. 4) then
-                        dms_surf(ji,jj) = dms_hall(ji,jj)
-                     endif
-                     !!
-                     !! 2D diag through iom_use
-                     IF( lk_iomput ) THEN
-                       IF( med_diag%DMS_SURF%dgsave ) THEN
-                         dms_surf2d(ji,jj) = dms_surf(ji,jj)
-                       ENDIF
-                       IF( med_diag%DMS_ANDR%dgsave ) THEN
-                         dms_andr2d(ji,jj) = dms_andr(ji,jj)
-                       ENDIF
-                       IF( med_diag%DMS_SIMO%dgsave ) THEN
-                         dms_simo2d(ji,jj) = dms_simo(ji,jj)
-                       ENDIF
-                       IF( med_diag%DMS_ARAN%dgsave ) THEN
-                         dms_aran2d(ji,jj) = dms_aran(ji,jj)
-                       ENDIF
-                       IF( med_diag%DMS_HALL%dgsave ) THEN
-                         dms_hall2d(ji,jj) = dms_hall(ji,jj)
-                       ENDIF
-#   if defined key_debug_medusa
-                       IF (lwp) write (numout,*) 'trc_bio_medusa: finish calculating dms'
-                     CALL flush(numout)
-#   endif 
-                     ENDIF
-                     !! End iom
-                  ENDIF
-                  !! End DMS Loop
-                  !!
-                  !! store 2D outputs
-                  !!
-                  !! JPALM -- 17-11-16 -- put fgco2 out of diag request
-                  !!                    is needed for coupling; pass through restart
-                  !! IF( med_diag%FGCO2%dgsave ) THEN
-                     !! convert from  mol/m2/day to kg/m2/s
-                     fgco2(ji,jj) = f_co2flux(ji,jj) * fse3t(ji,jj,jk) * CO2flux_conv  !! mmol-C/m3/d -> kg-CO2/m2/s
-                  !! ENDIF
-                  IF ( lk_iomput ) THEN
-                      IF( med_diag%ATM_PCO2%dgsave ) THEN
-                         f_pco2a2d(ji,jj) = f_pco2atm(ji,jj)
-                      ENDIF
-                      IF( med_diag%OCN_PCO2%dgsave ) THEN
-                         f_pco2w2d(ji,jj) = f_pco2w(ji,jj)
-                      ENDIF
-                      IF( med_diag%CO2FLUX%dgsave ) THEN
-                         f_co2flux2d(ji,jj) = f_co2flux(ji,jj) * fse3t(ji,jj,jk)           !! mmol/m3/d -> mmol/m2/d
-                      ENDIF
-                      IF( med_diag%TCO2%dgsave ) THEN
-                         f_TDIC2d(ji,jj) = f_TDIC(ji,jj)
-                      ENDIF
-                      IF( med_diag%TALK%dgsave ) THEN
-                         f_TALK2d(ji,jj) = f_TALK(ji,jj)
-                      ENDIF
-                      IF( med_diag%KW660%dgsave ) THEN
-                         f_kw6602d(ji,jj) = f_kw660(ji,jj)
-                      ENDIF
-                      IF( med_diag%ATM_PP0%dgsave ) THEN
-                         f_pp02d(ji,jj) = f_pp0(ji,jj)
-                      ENDIF
-                      IF( med_diag%O2FLUX%dgsave ) THEN
-                         f_o2flux2d(ji,jj) = f_o2flux(ji,jj)
-                      ENDIF
-                      IF( med_diag%O2SAT%dgsave ) THEN
-                         f_o2sat2d(ji,jj) = f_o2sat(ji,jj)
-                      ENDIF
-                      !! AXY (24/11/16): add in extra MOCSY diagnostics
-                      IF( med_diag%ATM_XCO2%dgsave ) THEN
-                         f_xco2a_2d(ji,jj) = f_xco2a(ji,jj)
-                      ENDIF
-                      IF( med_diag%OCN_FCO2%dgsave ) THEN
-                         f_fco2w_2d(ji,jj) = f_fco2w(ji,jj)
-                      ENDIF
-                      IF( med_diag%ATM_FCO2%dgsave ) THEN
-                         f_fco2a_2d(ji,jj) = f_fco2atm(ji,jj)
-                      ENDIF
-                      IF( med_diag%OCN_RHOSW%dgsave ) THEN
-                         f_ocnrhosw_2d(ji,jj) = f_rhosw(ji,jj)
-                      ENDIF
-                      IF( med_diag%OCN_SCHCO2%dgsave ) THEN
-                         f_ocnschco2_2d(ji,jj) = f_schmidtco2(ji,jj)
-                      ENDIF
-                      IF( med_diag%OCN_KWCO2%dgsave ) THEN
-                         f_ocnkwco2_2d(ji,jj) = f_kwco2(ji,jj)
-                      ENDIF
-                      IF( med_diag%OCN_K0%dgsave ) THEN
-                         f_ocnk0_2d(ji,jj) = f_K0(ji,jj)
-                      ENDIF
-                      IF( med_diag%CO2STARAIR%dgsave ) THEN
-                         f_co2starair_2d(ji,jj) = f_co2starair(ji,jj)
-                      ENDIF
-                      IF( med_diag%OCN_DPCO2%dgsave ) THEN
-                         f_ocndpco2_2d(ji,jj) = f_dpco2(ji,jj)
-                      ENDIF
-                  ENDIF
-                  !! 
-               endif
-               !! End jk = 1 loop within ROAM key 
+
+! Moved from above - marc 21/4/17
+! I think this should be moved into diagnostics at bottom - it
+! doesn't like it is used anyway else - marc 21/4/17
+         DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            !! OPEN wet point IF..THEN loop
+            if (tmask(ji,jj,jk) == 1) then
 
                !! AXY (11/11/16): CMIP6 oxygen saturation 3D diagnostic
@@ -1025,81 +839,13 @@
                   o2sat3(ji, jj, jk) = f_o2sat3
                ENDIF
-
-# endif
-
-               if ( jk .eq. 1 ) then
-                  !!----------------------------------------------------------------------
-                  !! River inputs
-                  !!----------------------------------------------------------------------
-                  !!
-                  !! runoff comes in as        kg / m2 / s
-                  !! used and written out as   m3 / m2 / d (= m / d)
-                  !! where                     1000 kg / m2 / d = 1 m3 / m2 / d = 1 m / d
-                  !!
-                  !! AXY (17/07/14): the compiler doesn't like this line for some reason;
-                  !!                 as MEDUSA doesn't even use runoff for riverine inputs,
-                  !!                 a temporary solution is to switch off runoff entirely
-                  !!                 here; again, this change is one of several that will 
-                  !!                 need revisiting once MEDUSA has bedded down in UKESM1;
-                  !!                 particularly so if the land scheme provides information
-                  !!                 concerning nutrient fluxes
-                  !!
-                  !! f_runoff(ji,jj) = sf_rnf(1)%fnow(ji,jj,1) / 1000. * 60. * 60. * 24.
-                  f_runoff(ji,jj) = 0.0
-                  !!
-                  !! nutrients are added via rivers to the model in one of two ways:
-                  !!   1. via river concentration; i.e. the average nutrient concentration
-                  !!      of a river water is described by a spatial file, and this is
-                  !!      multiplied by runoff to give a nutrient flux
-                  !!   2. via direct river flux; i.e. the average nutrient flux due to
-                  !!      rivers is described by a spatial file, and this is simply applied
-                  !!      as a direct nutrient flux (i.e. it does not relate or respond to
-                  !!      model runoff)
-                  !! nutrient fields are derived from the GlobalNEWS 2 database; carbon and
-                  !! alkalinity are derived from continent-scale DIC estimates (Huang et al., 
-                  !! 2012) and some Arctic river alkalinity estimates (Katya?)
-                  !! 
-                  !! as of 19/07/12, riverine nutrients can now be spread vertically across 
-                  !! several grid cells rather than just poured into the surface box; this
-                  !! block of code is still executed, however, to set up the total amounts
-                  !! of nutrient entering via rivers
-                  !!
-                  !! nitrogen
-                  if (jriver_n .eq. 1) then
-                     !! river concentration specified; use runoff to calculate input
-                     f_riv_n(ji,jj) = f_runoff(ji,jj) * riv_n(ji,jj)
-                  elseif (jriver_n .eq. 2) then
-                     !! river flux specified; independent of runoff
-                     f_riv_n(ji,jj) = riv_n(ji,jj)
-                  endif
-                  !!
-                  !! silicon
-                  if (jriver_si .eq. 1) then
-                     !! river concentration specified; use runoff to calculate input
-                     f_riv_si(ji,jj) = f_runoff(ji,jj) * riv_si(ji,jj)
-                  elseif (jriver_si .eq. 2) then
-                     !! river flux specified; independent of runoff
-                     f_riv_si(ji,jj) = riv_si(ji,jj)
-                  endif
-                  !!
-                  !! carbon
-                  if (jriver_c .eq. 1) then
-                     !! river concentration specified; use runoff to calculate input
-                     f_riv_c(ji,jj) = f_runoff(ji,jj) * riv_c(ji,jj)
-                  elseif (jriver_c .eq. 2) then
-                     !! river flux specified; independent of runoff
-                     f_riv_c(ji,jj) = riv_c(ji,jj)
-                  endif
-                  !!
-                  !! alkalinity
-                  if (jriver_alk .eq. 1) then
-                     !! river concentration specified; use runoff to calculate input
-                     f_riv_alk(ji,jj) = f_runoff(ji,jj) * riv_alk(ji,jj)
-                  elseif (jriver_alk .eq. 2) then
-                     !! river flux specified; independent of runoff
-                     f_riv_alk(ji,jj) = riv_alk(ji,jj)
-                  endif
-
-               endif
+            ENDIF
+         ENDDO
+         ENDDO
+# endif
+
+         DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            !! OPEN wet point IF..THEN loop
+            if (tmask(ji,jj,jk) == 1) then
 
                !!----------------------------------------------------------------------
@@ -1370,4 +1116,7 @@
                endif
 # endif
+
+! MAYBE BUT A BREAK IN HERE, ZOOPLANKTON GRAZING - marc 20/4/17 
+! (plankton growth is 281 lines)
 
                !!----------------------------------------------------------------------
@@ -1520,4 +1269,7 @@
                   fmegrow(ji,jj) + (xphi * (fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj))) + fmeexcr(ji,jj) + ((1.0 - xbetan) * finme(ji,jj)) )
 
+! MAYBE BUT A BREAK IN HERE, MISCELLANEOUS PLANKTON LOSSES - marc 20/4/17 
+! (zoo plankton grazing is 152 lines)
+
                !!----------------------------------------------------------------------
                !! Plankton metabolic losses
@@ -1569,4 +1321,7 @@
                if (jmzme.eq.4) fdzme(ji,jj) = xmzme * zzme(ji,jj) * &        !! sigmoid
                   ((zzme(ji,jj) * zzme(ji,jj)) / (xkzme + (zzme(ji,jj) * zzme(ji,jj))))
+
+! MAYBE BUT A BREAK IN HERE, DETRITUS PROCESS- marc 20/4/17 
+! (misc plankton loses is 53 lines - maybe we can link with another section)
 
                !!----------------------------------------------------------------------
@@ -1643,4 +1398,7 @@
                endif
 
+! MAYBE BUT A BREAK IN HERE, IRON CHEMISTRY AND SCAVENGING - marc 20/4/17 
+! (detritus processes is 74 lines)
+
                !!----------------------------------------------------------------------
                !! Iron chemistry and fractionation
@@ -1978,4 +1736,7 @@
 # endif
 
+! MAYBE BUT A BREAK IN HERE, MISCELLANEOUS PROCESSES - marc 20/4/17 
+! (iron chemistry and scavenging is 340 lines)
+
                !!----------------------------------------------------------------------
                !! Miscellaneous
@@ -2053,4 +1814,7 @@
 !               (xthetazmi * fdzmi2(ji,jj)) + (xthetazme * fdzme2(ji,jj))) * fse3t(ji,jj,jk))        ! linear mortality
 # endif
+
+! MAYBE BUT A BREAK IN HERE, FAST-SINKINIG DETRITUS - marc 20/4/17 
+! (miscellaneous processes is 79 lines)
 
                !!----------------------------------------------------------------------
@@ -2628,4 +2392,7 @@
                CALL flush(numout)
 
+! MAYBE BUT A BREAK IN HERE, BUSINESS AND UPDATING - marc 20/4/17 
+! (fast sinking detritus is 576 lines)
+
                !!======================================================================
                !! LOCAL GRID CELL TRENDS
@@ -3069,4 +2836,8 @@
                !!    endif
                !! endif     
+
+! MAYBE BUT A BREAK IN HERE - marc 20/4/17 
+! (this would make the previous section about 470 lines and the one below
+! about 700 lines)
 
 # if defined key_trc_diabio
