Changeset 7886

Timestamp:

2017-04-10T17:01:50+02:00 (7 years ago)

Author:

marc

Message:

Changing many of the scalars to 2d arrays in prepartion for implementing Maff's changes

File:

• branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 (modified) (135 diffs)

branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90

@@ -155 +155 @@
       !!
       !! model state variables
-      REAL(wp) ::    zchn,zchd,zphn,zphd,zpds,zzmi
-      REAL(wp) ::    zzme,zdet,zdtc,zdin,zsil,zfer
+      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) ::    zdic, zalk, zoxy
-      REAL(wp) ::    ztmp, zsal
+      REAL(wp), DIMENSION(jpi,jpj) ::    zdic, zalk, zoxy
+      REAL(wp), DIMENSION(jpi,jpj) ::    ztmp, zsal
 # endif
 # if defined key_mocsy
-      REAL(wp) ::    zpho
+      REAL(wp), DIMENSION(jpi,jpj) ::    zpho
 # endif
       !!
@@ -170 +170 @@
       !! AXY (23/08/13): changed from individual variables for each flux to
       !!                 an array that holds all fluxes
-      REAL(wp), DIMENSION(jp_medusa) ::    btra
+      REAL(wp), DIMENSION(jpi,jpj,jp_medusa) ::    btra
       !!
       !! primary production and chl related quantities      
-      REAL(wp)                     ::    fthetan,faln,fchn1,fchn,fjln,fprn,frn
-      REAL(wp)                     ::    fthetad,fald,fchd1,fchd,fjld,fprd,frd
+      REAL(wp), DIMENSION(jpi,jpj) ::    fthetan,faln,fchn1,fchn,fjln,fprn,frn
+      REAL(wp), DIMENSION(jpi,jpj) ::    fthetad,fald,fchd1,fchd,fjld,fprd,frd
       !! AXY (23/11/16): add in light-only limitation term (normalised 0-1 range)
-      REAL(wp)                     ::    fjlim_pn, fjlim_pd
+      REAL(wp), DIMENSION(jpi,jpj) ::    fjlim_pn, fjlim_pd
       !! AXY (03/02/11): add in Liebig terms
-      REAL(wp) ::    fpnlim, fpdlim
+      REAL(wp), DIMENSION(jpi,jpj) ::    fpnlim, fpdlim
       !! AXY (16/07/09): add in Eppley curve functionality
-      REAL(wp) ::    loc_T,fun_T,xvpnT,xvpdT
+      REAL(wp), DIMENSION(jpi,jpj) ::    fun_T,xvpnT,xvpdT
       INTEGER  ::    ieppley
       !! AXY (16/05/11): per Katya's prompting, add in new T-dependence
       !!                 for phytoplankton growth only (i.e. no change
       !!                 for remineralisation)
-      REAL(wp) ::    fun_Q10
+      REAL(wp), DIMENSION(jpi,jpj) ::    fun_Q10
       !! AXY (01/03/10): add in mixed layer PP diagnostics
-      REAL(wp), DIMENSION(jpi,jpj) ::  fprn_ml,fprd_ml
+      REAL(wp), DIMENSION(jpi,jpj) ::    fprn_ml,fprd_ml
       !!
       !! nutrient limiting factors
-      REAL(wp) ::    fnln,ffln            !! N and Fe
-      REAL(wp) ::    fnld,ffld,fsld,fsld2 !! N, Fe and Si
+      REAL(wp), DIMENSION(jpi,jpj) ::    fnln,ffln            !! N and Fe
+      REAL(wp), DIMENSION(jpi,jpj) ::    fnld,ffld,fsld,fsld2 !! N, Fe and Si
       !!
       !! silicon cycle
-      REAL(wp) ::    fsin,fnsi,fsin1,fnsi1,fnsi2,fprds,fsdiss
+      REAL(wp), DIMENSION(jpi,jpj) ::    fsin,fnsi,fprds,fsdiss
+      REAL(wp)                     ::    fsin1,fnsi1,fnsi2
       !!
       !! iron cycle; includes parameters for Parekh et al. (2005) iron scheme
-      REAL(wp) ::    ffetop,ffebot,ffescav
+      REAL(wp), DIMENSION(jpi,jpj) ::    ffetop,ffebot,ffescav
       REAL(wp) ::    xLgF, xFeT, xFeF, xFeL         !! state variables for iron-ligand system
       REAL(wp), DIMENSION(jpi,jpj) ::  xFree        !! state variables for iron-ligand system
@@ -214 +215 @@
       !!
       !! microzooplankton grazing
-      REAL(wp) ::    fmi1,fmi,fgmipn,fgmid,fgmidc
-      REAL(wp) ::    finmi,ficmi,fstarmi,fmith,fmigrow,fmiexcr,fmiresp
+      REAL(wp), DIMENSION(jpi,jpj) ::    fmi1,fmi,fgmipn,fgmid,fgmidc
+      REAL(wp), DIMENSION(jpi,jpj) ::    finmi,ficmi,fstarmi,fmith,fmigrow,fmiexcr,fmiresp
       !!
       !! mesozooplankton grazing
-      REAL(wp) ::    fme1,fme,fgmepn,fgmepd,fgmepds,fgmezmi,fgmed,fgmedc
-      REAL(wp) ::    finme,ficme,fstarme,fmeth,fmegrow,fmeexcr,fmeresp
+      REAL(wp), DIMENSION(jpi,jpj) ::    fme1,fme,fgmepn,fgmepd,fgmepds,fgmezmi,fgmed,fgmedc
+      REAL(wp), DIMENSION(jpi,jpj) ::    finme,ficme,fstarme,fmeth,fmegrow,fmeexcr,fmeresp
       !!
       !! mortality/Remineralisation (defunct parameter "fz" removed)
-      REAL(wp) ::    fdpn,fdpd,fdpds,fdzmi,fdzme,fdd
+      REAL(wp, DIMENSION(jpi,jpj)) ::    fdpn,fdpd,fdpds,fdzmi,fdzme,fdd
 # if defined key_roam
-      REAL(wp) ::    fddc
-# endif
-      REAL(wp) ::    fdpn2,fdpd2,fdpds2,fdzmi2,fdzme2
-      REAL(wp) ::    fslown, fslowc
+      REAL(wp, DIMENSION(jpi,jpj)) ::    fddc
+# endif
+      REAL(wp), DIMENSION(jpi,jpj) ::    fdpn2,fdpd2,fdpds2,fdzmi2,fdzme2
+      REAL(wp), DIMENSION(jpi,jpj) ::    fslown, fslowc
       REAL(wp), DIMENSION(jpi,jpj) ::    fslownflux, fslowcflux
-      REAL(wp) ::    fregen,fregensi
+      REAL(wp), DIMENSION(jpi,jpj) ::    fregen,fregensi
       REAL(wp), DIMENSION(jpi,jpj) ::    fregenfast,fregenfastsi
 # if defined key_roam
-      REAL(wp) ::    fregenc
+! Doesn't look like this is used - marc 10/4/17
+!      REAL(wp), DIMENSION(jpi,jpj) ::    fregenc
       REAL(wp), DIMENSION(jpi,jpj) ::    fregenfastc
 # endif
       !!
       !! particle flux
-      REAL(WP) ::    fthk,fdep,fdep1,fdep2,flat,fcaco3
-      REAL(WP) ::    ftempn,ftempsi,ftempfe,ftempc,ftempca
-      REAL(wp) ::    freminn,freminsi,freminfe,freminc,freminca
+      REAL(WP), DIMENSION(jpi,jpj) ::    fdep1,fcaco3
+      REAL(WP), DIMENSION(jpi,jpj) ::    ftempn,ftempsi,ftempfe,ftempc,ftempca
+      REAL(wp), DIMENSION(jpi,jpj) ::    freminn,freminsi,freminfe,freminc,freminca
       REAL(wp), DIMENSION(jpi,jpj) ::    ffastn,ffastsi,ffastfe,ffastc,ffastca
-      REAL(wp) ::    fleftn,fleftsi,fleftfe,fleftc,fleftca
-      REAL(wp) ::    fheren,fheresi,fherefe,fherec,fhereca
-      REAL(wp) ::    fprotf
+      REAL(wp), DIMENSION(jpi,jpj) ::    fprotf
       REAL(wp), DIMENSION(jpi,jpj) ::    fsedn,fsedsi,fsedfe,fsedc,fsedca
       REAL(wp), DIMENSION(jpi,jpj) ::    fccd
-      REAL(wp) ::    fccd_dep
-      !! AXY (28/11/16): fix mbathy bug
-      INTEGER  ::    jmbathy
+      REAL(wp), DIMENSION(jpi,jpj) ::    fccd_dep
       !!
       !! AXY (06/07/11): alternative fast detritus schemes
@@ -254 +252 @@
       !!
       !! AXY (08/07/11): fate of fast detritus reaching the seafloor
-      REAL(wp) ::    ffast2slown,ffast2slowfe,ffast2slowc
+! I don't think ffast2slowfe is used - marc 10/4/17
+!      REAL(wp), DIMENSION(jpi,jpj) ::    ffast2slown,ffast2slowfe,ffast2slowc
+      REAL(wp), DIMENSION(jpi,jpj) ::    ffast2slown,ffast2slow
       !!
       !! conservation law
@@ -299 +299 @@
       INTEGER, DIMENSION(jpi,jpj) ::     i2_omcal,i2_omarg
       !!
-      !! ROAM air-sea flux and diagnostic parameters
-      REAL(wp) ::    f_wind
       !! AXY (24/11/16): add xCO2 variable for atmosphere (what we actually have)
-      REAL(wp) ::    f_xco2a
-      REAL(wp) ::    f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_co2flux
-      REAL(wp) ::    f_TDIC, f_TALK, f_dcf, f_henry
-      REAL(wp) ::    f_uwind, f_vwind, f_pp0
-      REAL(wp) ::    f_kw660, f_o2flux, f_o2sat, f_o2sat3
+      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, f_o2sat3
       REAL(wp), DIMENSION(jpi,jpj) ::    f_omcal, f_omarg
       !!
       !! AXY (23/06/15): additional diagnostics for MOCSY and oxygen
-      REAL(wp) ::    f_fco2w, f_BetaD, f_rhosw, f_opres, f_insitut, f_pco2atm, f_fco2atm
-      REAL(wp) ::    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
@@ -318 +316 @@
 
       !!
-      INTEGER  ::    iters
+      INTEGER, DIMENSION(jpi,jpj)  ::    iters
       REAL(wp) ::    f_year
       INTEGER  ::    i_year
@@ -329 +327 @@
       !!
       !! oxygen production and consumption (and non-consumption)
-      REAL(wp) ::    fo2_prod, fo2_cons, fo2_ncons, fo2_ccons
+      REAL(wp), DIMENSION(jpi,jpj) ::    fo2_prod, fo2_cons, fo2_ncons, fo2_ccons
       REAL(wp), DIMENSION(jpi,jpj) ::    foxy_prod, foxy_cons, foxy_anox
       !! Jpalm (11-08-2014)
       !! add DMS in MEDUSA for UKESM1 model
-      REAL(wp) ::    dms_surf
+      REAL(wp), DIMENSION(jpi,jpj) ::    dms_surf
       !! AXY (13/03/15): add in other DMS calculations
-      REAL(wp) ::    dms_andr, dms_simo, dms_aran, dms_hall
+      REAL(wp), DIMENSION(jpi,jpj) ::    dms_andr, dms_simo, dms_aran, dms_hall
 
 # endif
@@ -352 +350 @@
       REAL(wp), DIMENSION(jpi,jpj) :: f_runoff, f_riv_n, f_riv_si, f_riv_c, f_riv_alk
       !! AXY (19/07/12): variables for local riverine fluxes to handle inputs below surface
-      REAL(wp) ::    f_riv_loc_n, f_riv_loc_si, f_riv_loc_c, f_riv_loc_alk
+      REAL(wp), DIMENSION(jpi,jpj) :: f_riv_loc_n, f_riv_loc_si, f_riv_loc_c, f_riv_loc_alk
       !!
       !! Jpalm -- 11-10-2015 -- adapt diag to iom_use
@@ -1359 +1357 @@
                      !! do carbonate chemistry
                      !!
-                     fdep2 = fsdept(ji,jj,jk)           !! set up level midpoint
-                     !! AXY (28/11/16): local seafloor depth
-                     !!                 previously mbathy(ji,jj) - 1, now mbathy(ji,jj)
-                     jmbathy = mbathy(ji,jj)
-                     !!
                      !! set up required state variables
-                     zdic = max(0.,trn(ji,jj,jk,jpdic)) !! dissolved inorganic carbon
-                     zalk = max(0.,trn(ji,jj,jk,jpalk)) !! alkalinity
-                     ztmp = tsn(ji,jj,jk,jp_tem)        !! temperature
-                     zsal = tsn(ji,jj,jk,jp_sal)        !! salinity
+                     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 = max(0.,trn(ji,jj,jk,jpsil))        !! silicic acid
-                     zpho = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 !! phosphate via DIN and Redfield
+                     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 .lt. -3.0 .or. ztmp .gt. 40.0 ) then
+           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 (',    &
@@ -1381 +1374 @@
          IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 3D, ', &
          tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem)
-                        ztmp = tsb(ji,jj,jk,jp_tem)     !! temperature
+                        ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem)     !! temperature
                      endif
-           if (zsal .lt. 0.0 .or. zsal .gt. 45.0 ) then
+           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 (',    &
@@ -1390 +1383 @@
                      !!
                      !! blank input variables not used at this stage (they relate to air-sea flux)
-                     f_kw660 = 1.0
-                     f_pp0   = 1.0
+                     f_kw660(ji,jj) = 1.0
+                     f_pp0(ji,jj)   = 1.0
                      !!
                      !! calculate carbonate chemistry at grid cell midpoint
@@ -1397 +1390 @@
                      !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate
                      !!                 chemistry package
-                     CALL mocsy_interface( ztmp, zsal, zalk, zdic, zsil, zpho,         &    ! inputs
-                     f_pp0, fdep2, gphit(ji,jj), f_kw660, f_xco2a, 1,                  &    ! inputs
-                     f_ph, f_pco2w, f_fco2w, f_h2co3, f_hco3, f_co3, f_omarg(ji,jj),   &    ! outputs
-                     f_omcal(ji,jj), f_BetaD, f_rhosw, f_opres, f_insitut,             &    ! outputs
-                     f_pco2atm, f_fco2atm, f_schmidtco2, f_kwco2, f_K0,                &    ! outputs
-                     f_co2starair, f_co2flux, f_dpco2 )                                     ! outputs
+                     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 = (zdic / f_rhosw) * 1000. ! mmol / m3 -> umol / kg
-                     f_TALK = (zalk / f_rhosw) * 1000. !  meq / m3 ->  ueq / kg
-                     f_dcf  = f_rhosw
+                     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, zsal, zdic, zalk, fdep2, f_kw660,      &    ! inputs
-                     f_xco2a, f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_omcal(ji,jj),   &    ! outputs
-                     f_omarg(ji,jj), f_co2flux, f_TDIC, f_TALK, f_dcf, f_henry, iters)      ! outputs
+                     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 .eq. 25) then
+                     if (iters(ji,jj) .eq. 25) then
                         IF(lwp) WRITE(numout,*) ' trc_bio_medusa: 3D ITERS WARNING, ', &
-                        iters, ' AT (', ji, ', ', jj, ', ', jk, ') AT ', kt
+                        iters(ji,jj), ' AT (', ji, ', ', jj, ', ', jk, ') AT ', kt
                      endif
 #  endif
                      !!
                      !! store 3D outputs
-                     f3_pH(ji,jj,jk)    = f_ph
-                     f3_h2co3(ji,jj,jk) = f_h2co3
-                     f3_hco3(ji,jj,jk)  = f_hco3
-                     f3_co3(ji,jj,jk)   = f_co3
+                     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)
@@ -1432 +1425 @@
                      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) = fdep2
+                           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 = fdep2 - fsdept(ji,jj,jk-1)
+                           fq3 = fsdept(ji,jj,jk) - fsdept(ji,jj,jk-1)
                            fq4 = fq2 * fq3
                            f2_ccd_cal(ji,jj) = fsdept(ji,jj,jk-1) + fq4
@@ -1443 +1436 @@
                         i2_omcal(ji,jj)   = 1
                      endif
-                     if ( i2_omcal(ji,jj) .eq. 0 .and. jk .eq. jmbathy ) then
+                     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)
@@ -1452 +1445 @@
                      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) = fdep2
+                           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 = fdep2 - fsdept(ji,jj,jk-1)
+                           fq3 = fsdept(ji,jj,jk) - fsdept(ji,jj,jk-1)
                            fq4 = fq2 * fq3
                            f2_ccd_arg(ji,jj) = fsdept(ji,jj,jk-1) + fq4
@@ -1463 +1456 @@
                         i2_omarg(ji,jj)   = 1
                      endif
-                     if ( i2_omarg(ji,jj) .eq. 0 .and. jk .eq. jmbathy ) then
+                     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)
@@ -1514 +1507 @@
                !!---------------------------------------------------------------------
                !!
-               !! AXY (11/12/08): set up level thickness
-               fthk  = fse3t(ji,jj,jk)
-               !! AXY (25/02/10): set up level depth (top of level)
-               fdep  = fsdepw(ji,jj,jk)
                !! AXY (01/03/10): set up level depth (bottom of level)
-               fdep1 = fdep + fthk
+               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)
-               jmbathy = mbathy(ji,jj)
+               mbathy(ji,jj) = mbathy(ji,jj)
                !!
                !! set up model tracers
                !! negative values of state variables are not allowed to
                !! contribute to the calculated fluxes
-               zchn = max(0.,trn(ji,jj,jk,jpchn)) !! non-diatom chlorophyll
-               zchd = max(0.,trn(ji,jj,jk,jpchd)) !! diatom chlorophyll
-               zphn = max(0.,trn(ji,jj,jk,jpphn)) !! non-diatoms
-               zphd = max(0.,trn(ji,jj,jk,jpphd)) !! diatoms
-               zpds = max(0.,trn(ji,jj,jk,jppds)) !! diatom silicon
+               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
-               if (zchn.eq.0.) zphn = 0.
-               if (zchd.eq.0.) zphd = 0.
-               if (zphn.eq.0.) zchn = 0.
-               if (zphd.eq.0.) zchd = 0.
+               if (zchn(ji,jj).eq.0.) zphn(ji,jj) = 0.
+               if (zchd(ji,jj).eq.0.) zphd(ji,jj) = 0.
+               if (zphn(ji,jj).eq.0.) zchn(ji,jj) = 0.
+               if (zphd(ji,jj).eq.0.) zchd(ji,jj) = 0.
           !! AXY (23/01/14): duh - why did I forget diatom silicon?
-          if (zpds.eq.0.) zphd = 0.
-          if (zphd.eq.0.) zpds = 0.
-               zzmi = max(0.,trn(ji,jj,jk,jpzmi)) !! microzooplankton
-               zzme = max(0.,trn(ji,jj,jk,jpzme)) !! mesozooplankton
-               zdet = max(0.,trn(ji,jj,jk,jpdet)) !! detrital nitrogen
-               zdin = max(0.,trn(ji,jj,jk,jpdin)) !! dissolved inorganic nitrogen
-               zsil = max(0.,trn(ji,jj,jk,jpsil)) !! dissolved silicic acid
-               zfer = max(0.,trn(ji,jj,jk,jpfer)) !! dissolved "iron"
+          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"
 # if defined key_roam
-               zdtc = max(0.,trn(ji,jj,jk,jpdtc)) !! detrital carbon
-               zdic = max(0.,trn(ji,jj,jk,jpdic)) !! dissolved inorganic carbon
-               zalk = max(0.,trn(ji,jj,jk,jpalk)) !! alkalinity
-               zoxy = max(0.,trn(ji,jj,jk,jpoxy)) !! oxygen
+               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
 #  if defined key_axy_carbchem && defined key_mocsy
-               zpho = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 !! phosphate via DIN and Redfield
+               zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 !! phosphate via DIN and Redfield
 #  endif
                !!
                !! also need physical parameters for gas exchange calculations
-               ztmp = tsn(ji,jj,jk,jp_tem)
-               zsal = tsn(ji,jj,jk,jp_sal)
+               ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem)
+               zsal(ji,jj) = tsn(ji,jj,jk,jp_sal)
                !!
           !! AXY (28/02/14): check input fields
-               if (ztmp .lt. -3.0 .or. ztmp .gt. 40.0 ) then
+               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 (',    &
@@ -1566 +1555 @@
         IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 2D, ', &
                   tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem)
-                  ztmp = tsb(ji,jj,jk,jp_tem) !! temperature
+                  ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem) !! temperature
                endif
-               if (zsal .lt. 0.0 .or. zsal .gt. 45.0 ) then
+               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 (',    &
@@ -1574 +1563 @@
                endif
 # else
-               zdtc = zdet * xthetad              !! implicit detrital carbon
+               zdtc(ji,jj) = zdet(ji,jj) * xthetad              !! implicit detrital carbon
 # endif
 # if defined key_debug_medusa
@@ -1612 +1601 @@
                if (idf.eq.1.AND.idfval.eq.1) then
                   IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'fthk(',jk,') = ', fthk
-                  IF (lwp) write (numout,*) 'zphn(',jk,') = ', zphn
-                  IF (lwp) write (numout,*) 'zphd(',jk,') = ', zphd
-                  IF (lwp) write (numout,*) 'zpds(',jk,') = ', zpds
-                  IF (lwp) write (numout,*) 'zzmi(',jk,') = ', zzmi
-                  IF (lwp) write (numout,*) 'zzme(',jk,') = ', zzme
-                  IF (lwp) write (numout,*) 'zdet(',jk,') = ', zdet
-                  IF (lwp) write (numout,*) 'zdin(',jk,') = ', zdin
-                  IF (lwp) write (numout,*) 'zsil(',jk,') = ', zsil
-                  IF (lwp) write (numout,*) 'zfer(',jk,') = ', zfer
+                  IF (lwp) write (numout,*) 'fthk(',jk,') = ', fse3t(ji,jj,jk)
+                  IF (lwp) write (numout,*) 'zphn(',jk,') = ', zphn(ji,jj)
+                  IF (lwp) write (numout,*) 'zphd(',jk,') = ', zphd(ji,jj)
+                  IF (lwp) write (numout,*) 'zpds(',jk,') = ', zpds(ji,jj)
+                  IF (lwp) write (numout,*) 'zzmi(',jk,') = ', zzmi(ji,jj)
+                  IF (lwp) write (numout,*) 'zzme(',jk,') = ', zzme(ji,jj)
+                  IF (lwp) write (numout,*) 'zdet(',jk,') = ', zdet(ji,jj)
+                  IF (lwp) write (numout,*) 'zdin(',jk,') = ', zdin(ji,jj)
+                  IF (lwp) write (numout,*) 'zsil(',jk,') = ', zsil(ji,jj)
+                  IF (lwp) write (numout,*) 'zfer(',jk,') = ', zfer(ji,jj)
 #  if defined key_roam
-                  IF (lwp) write (numout,*) 'zdtc(',jk,') = ', zdtc
-                  IF (lwp) write (numout,*) 'zdic(',jk,') = ', zdic
-                  IF (lwp) write (numout,*) 'zalk(',jk,') = ', zalk
-                  IF (lwp) write (numout,*) 'zoxy(',jk,') = ', zoxy                  
+                  IF (lwp) write (numout,*) 'zdtc(',jk,') = ', zdtc(ji,jj)
+                  IF (lwp) write (numout,*) 'zdic(',jk,') = ', zdic(ji,jj)
+                  IF (lwp) write (numout,*) 'zalk(',jk,') = ', zalk(ji,jj)
+                  IF (lwp) write (numout,*) 'zoxy(',jk,') = ', zoxy(ji,jj)                  
 #  endif
                endif
@@ -1642 +1631 @@
                   IF ( med_diag%INVTN%dgsave )   THEN
                      ftot_n(ji,jj)  = ftot_n(ji,jj) + &
-                             (fthk * ( zphn + zphd + zzmi + zzme + zdet + zdin ) )
+                             (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) + & 
-                             (fthk * ( zpds + zsil ) )
+                             (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) + & 
-                             (fthk * ( xrfn * ( zphn + zphd + zzmi + zzme + zdet ) + zfer ) )
+                             (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) + & 
-                             (fthk * ( (xthetapn * zphn) + (xthetapd * zphd) + &
-                             (xthetazmi * zzmi) + (xthetazme * zzme) + zdtc +   &
-                             zdic ) )
+                             (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) + (fthk * ( zalk ) )
+                     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) + (fthk * ( zoxy ) )
+                     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) + (fthk * zdic)
+                     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) + (fthk * zdin)
+                     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) + (fthk * zsil)
+                     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) + (fthk * zalk)
+                     inttalk(ji,jj) = inttalk(ji,jj) + (fse3t(ji,jj,jk) * zalk(ji,jj))
                   ENDIF
                   IF ( med_diag%O2min%dgsave ) THEN
-                     if ( zoxy < o2min(ji,jj) ) then
-                        o2min(ji,jj)  = zoxy
+                     if ( zoxy(ji,jj) < o2min(ji,jj) ) then
+                        o2min(ji,jj)  = zoxy(ji,jj)
                         IF ( med_diag%ZO2min%dgsave ) THEN
-                           zo2min(ji,jj) = (fdep + fdep1) / 2. !! layer midpoint
+                           zo2min(ji,jj) = (fsdepw(ji,jj,jk) + fdep1(ji,jj)) / 2. !! layer midpoint
                         ENDIF
                      endif
@@ -1701 +1690 @@
                   !! Air-sea gas exchange
                   !!----------------------------------------------------------------------
-                  !!
-                  !! AXY (17/07/14): zwind_i and zwind_j do not exist in this
-                  !!                 version of NEMO because it does not include
-                  !!                 the SBC changes that our local version has
-                  !!                 for accessing the HadGEM2 forcing; they 
-                  !!                 could be added, but an alternative approach
-                  !!                 is to make use of wndm from oce_trc.F90
-                  !!                 which is wind speed at 10m (which is what
-                  !!                 is required here; this may need to be
-                  !!                 revisited when MEDUSA properly interacts
-                  !!                 with UKESM1 physics
-                  !!
-                  f_wind  = wndm(ji,jj)
                   IF (lk_oasis) THEN
                      f_xco2a = PCO2a_in_cpl(ji,jj)        !! use 2D atm xCO2 from atm coupling
@@ -1731 +1707 @@
                      CALL flush(numout)
 #   endif
-                  CALL gas_transfer( f_wind, 1, 7, &  ! inputs
-                                     f_kw660 )        ! outputs
+                  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'
@@ -1741 +1717 @@
                   !! of the UKESM1 atmosphere 
                   !!                                     
-                  f_pp0   = 1.0
-                  !!
-                  !! IF(lwp) WRITE(numout,*) ' MEDUSA ztmp    =', ztmp
+                  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 f_wind  =', f_wind
+                  !! IF(lwp) WRITE(numout,*) ' MEDUSA wndm    =', wndm(ji,jj)
                   !! IF(lwp) WRITE(numout,*) ' MEDUSA fr_i    =', fr_i(ji,jj)
                   !!
@@ -1755 +1731 @@
                   !!                 chemistry package; note that depth is set to
                   !!                 zero in this call
-                  CALL mocsy_interface( ztmp, zsal, zalk, zdic, zsil, zpho,        &  ! inputs
-                  f_pp0, 0.0, gphit(ji,jj), f_kw660, f_xco2a, 1,                   &  ! inputs
-                  f_ph, f_pco2w, f_fco2w, f_h2co3, f_hco3, f_co3, f_omarg(ji,jj),  &  ! outputs
-                  f_omcal(ji,jj), f_BetaD, f_rhosw, f_opres, f_insitut,            &  ! outputs
-                  f_pco2atm, f_fco2atm, f_schmidtco2, f_kwco2, f_K0,               &  ! outputs
-                  f_co2starair, f_co2flux, f_dpco2 )                                  ! outputs
-                  !!
-                  f_TDIC = (zdic / f_rhosw) * 1000. ! mmol / m3 -> umol / kg
-                  f_TALK = (zalk / f_rhosw) * 1000. !  meq / m3 ->  ueq / kg
-                  f_dcf  = f_rhosw
+                  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_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 = 0
+                  iters(ji,jj) = 0
                   !!
                   !! carbon dioxide (CO2); Jerry Blackford code (ostensibly OCMIP-2, but not)
-                  CALL trc_co2_medusa( ztmp, zsal, zdic, zalk, 0.0, f_kw660, f_xco2a,  &  ! inputs
-                  f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_omcal(ji,jj),               &  ! outputs
-                  f_omarg(ji,jj), f_co2flux, f_TDIC, f_TALK, f_dcf, f_henry, iters )      ! outputs
+                  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
+                  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
@@ -1779 +1755 @@
                   !!                 output warnings are retained here so that
                   !!                 failure position can be determined
-                  if (iters .eq. 25) then
+                  if (iters(ji,jj) .eq. 25) then
                      IF(lwp) WRITE(numout,*) ' trc_bio_medusa: ITERS WARNING, ', &
-                     iters, ' AT (', ji, ', ', jj, ', ', jk, ') AT ', kt
+                     iters(ji,jj), ' AT (', ji, ', ', jj, ', ', jk, ') AT ', kt
                   endif
 #   endif
@@ -1787 +1763 @@
                   !! AXY (18/04/13): switch off carbonate chemistry calculations; provide
                   !!                 quasi-sensible alternatives
-                  f_ph           = 8.1
-                  f_pco2w        = f_xco2a
-                  f_h2co3        = 0.005 * zdic
-                  f_hco3         = 0.865 * zdic
-                  f_co3          = 0.130 * zdic
+                  f_ph(ji,jj)           = 8.1
+                  f_pco2w(ji,jj)        = f_xco2a
+                  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      = 0.
-                  f_TDIC         = zdic
-                  f_TALK         = zalk
-                  f_dcf          = 1.026
-                  f_henry        = 1.
+                  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        = f_xco2a
-                  f_BetaD        = 1.
-                  f_rhosw        = 1.026
-                  f_opres        = 0.
-                  f_insitut      = ztmp
-                  f_pco2atm      = f_xco2a
-                  f_fco2atm      = f_xco2a
-                  f_schmidtco2   = 660.
-                  f_kwco2        = 0.
-                  f_K0           = 0.
-                  f_co2starair   = f_xco2a
-                  f_dpco2        = 0.
+                  f_fco2w(ji,jj)        = f_xco2a
+                  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_schmidtco2(ji,jj)   = 660.
+                  f_kwco2(ji,jj)        = 0.
+                  f_K0(ji,jj)           = 0.
+                  f_co2starair(ji,jj)   = f_xco2a
+                  f_dpco2(ji,jj)        = 0.
 #  endif
                   !!
                   !! mmol/m2/s -> mmol/m3/d; correct for sea-ice; divide through by layer thickness
-                  f_co2flux = (1. - fr_i(ji,jj)) * f_co2flux * 86400. / fthk
+                  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
-                  !! CALL trc_oxy_medusa( ztmp, zsal, f_uwind, f_vwind, f_pp0, zoxy / 1000., fthk,  &  ! inputs
-                  !! f_kw660, f_o2flux, f_o2sat )                                                      ! outputs
-                  CALL trc_oxy_medusa( ztmp, zsal, f_kw660, f_pp0, zoxy,  &  ! inputs
-                  f_kwo2, f_o2flux, f_o2sat )                                ! outputs
+                  !! 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  = (1. - fr_i(ji,jj)) * f_o2flux * 86400. / fthk
+                  f_o2flux(ji,jj)  = (1. - fr_i(ji,jj)) * f_o2flux(ji,jj) * 86400. / fse3t(ji,jj,jk)
                   !!
                   !! Jpalm (08-2014)
@@ -1846 +1822 @@
                      if (jdms_input .eq. 0) then
                         !! use instantaneous inputs
-                        CALL trc_dms_medusa( zchn, zchd, hmld(ji,jj), qsr(ji,jj), zdin, &  ! inputs
-                        dms_andr, dms_simo, dms_aran, dms_hall )                           ! outputs
+                        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, dms_simo, dms_aran, dms_hall )                      ! outputs
+                        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 = dms_andr
+                        dms_surf(ji,jj) = dms_andr(ji,jj)
                      elseif (jdms_model .eq. 2) then
-                        dms_surf = dms_simo
+                        dms_surf(ji,jj) = dms_simo(ji,jj)
                      elseif (jdms_model .eq. 3) then
-                        dms_surf = dms_aran
+                        dms_surf(ji,jj) = dms_aran(ji,jj)
                      elseif (jdms_model .eq. 4) then
-                        dms_surf = dms_hall
+                        dms_surf(ji,jj) = dms_hall(ji,jj)
                      endif
                      !!
@@ -1869 +1845 @@
                      IF( lk_iomput ) THEN
                        IF( med_diag%DMS_SURF%dgsave ) THEN
-                         dms_surf2d(ji,jj) = dms_surf
+                         dms_surf2d(ji,jj) = dms_surf(ji,jj)
                        ENDIF
                        IF( med_diag%DMS_ANDR%dgsave ) THEN
-                         dms_andr2d(ji,jj) = dms_andr
+                         dms_andr2d(ji,jj) = dms_andr(ji,jj)
                        ENDIF
                        IF( med_diag%DMS_SIMO%dgsave ) THEN
-                         dms_simo2d(ji,jj) = dms_simo
+                         dms_simo2d(ji,jj) = dms_simo(ji,jj)
                        ENDIF
                        IF( med_diag%DMS_ARAN%dgsave ) THEN
-                         dms_aran2d(ji,jj) = dms_aran
+                         dms_aran2d(ji,jj) = dms_aran(ji,jj)
                        ENDIF
                        IF( med_diag%DMS_HALL%dgsave ) THEN
-                         dms_hall2d(ji,jj) = dms_hall
+                         dms_hall2d(ji,jj) = dms_hall(ji,jj)
                        ENDIF
 #   if defined key_debug_medusa
@@ -1898 +1874 @@
                   !! IF( med_diag%FGCO2%dgsave ) THEN
                      !! convert from  mol/m2/day to kg/m2/s
-                     fgco2(ji,jj) = f_co2flux * fthk * CO2flux_conv  !! mmol-C/m3/d -> kg-CO2/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
+                         f_pco2a2d(ji,jj) = f_pco2atm(ji,jj)
                       ENDIF
                       IF( med_diag%OCN_PCO2%dgsave ) THEN
-                         f_pco2w2d(ji,jj) = f_pco2w
+                         f_pco2w2d(ji,jj) = f_pco2w(ji,jj)
                       ENDIF
                       IF( med_diag%CO2FLUX%dgsave ) THEN
-                         f_co2flux2d(ji,jj) = f_co2flux * fthk           !! mmol/m3/d -> mmol/m2/d
+                         f_co2flux(ji,jj)2d(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
+                         f_TDIC2d(ji,jj) = f_TDIC(ji,jj)
                       ENDIF
                       IF( med_diag%TALK%dgsave ) THEN
-                         f_TALK2d(ji,jj) = f_TALK
+                         f_TALK2d(ji,jj) = f_TALK(ji,jj)
                       ENDIF
                       IF( med_diag%KW660%dgsave ) THEN
-                         f_kw6602d(ji,jj) = f_kw660
+                         f_kw6602d(ji,jj) = f_kw660(ji,jj)
                       ENDIF
                       IF( med_diag%ATM_PP0%dgsave ) THEN
-                         f_pp02d(ji,jj) = f_pp0
+                         f_pp02d(ji,jj) = f_pp0(ji,jj)
                       ENDIF
                       IF( med_diag%O2FLUX%dgsave ) THEN
-                         f_o2flux2d(ji,jj) = f_o2flux
+                         f_o2flux2d(ji,jj) = f_o2flux(ji,jj)
                       ENDIF
                       IF( med_diag%O2SAT%dgsave ) THEN
-                         f_o2sat2d(ji,jj) = f_o2sat
+                         f_o2sat2d(ji,jj) = f_o2sat(ji,jj)
                       ENDIF
                       !! AXY (24/11/16): add in extra MOCSY diagnostics
@@ -1933 +1909 @@
                       ENDIF
                       IF( med_diag%OCN_FCO2%dgsave ) THEN
-                         f_fco2w_2d(ji,jj) = f_fco2w
+                         f_fco2w_2d(ji,jj) = f_fco2w(ji,jj)
                       ENDIF
                       IF( med_diag%ATM_FCO2%dgsave ) THEN
-                         f_fco2a_2d(ji,jj) = f_fco2atm
+                         f_fco2a_2d(ji,jj) = f_fco2atm(ji,jj)
                       ENDIF
                       IF( med_diag%OCN_RHOSW%dgsave ) THEN
-                         f_ocnrhosw_2d(ji,jj) = f_rhosw
+                         f_ocnrhosw_2d(ji,jj) = f_rhosw(ji,jj)
                       ENDIF
                       IF( med_diag%OCN_SCHCO2%dgsave ) THEN
-                         f_ocnschco2_2d(ji,jj) = f_schmidtco2
+                         f_ocnschco2_2d(ji,jj) = f_schmidtco2(ji,jj)
                       ENDIF
                       IF( med_diag%OCN_KWCO2%dgsave ) THEN
-                         f_ocnkwco2_2d(ji,jj) = f_kwco2
+                         f_ocnkwco2_2d(ji,jj) = f_kwco2(ji,jj)
                       ENDIF
                       IF( med_diag%OCN_K0%dgsave ) THEN
-                         f_ocnk0_2d(ji,jj) = f_K0
+                         f_ocnk0_2d(ji,jj) = f_K0(ji,jj)
                       ENDIF
                       IF( med_diag%CO2STARAIR%dgsave ) THEN
-                         f_co2starair_2d(ji,jj) = f_co2starair
+                         f_co2starair_2d(ji,jj) = f_co2starair(ji,jj)
                       ENDIF
                       IF( med_diag%OCN_DPCO2%dgsave ) THEN
-                         f_ocndpco2_2d(ji,jj) = f_dpco2
+                         f_ocndpco2_2d(ji,jj) = f_dpco2(ji,jj)
                       ENDIF
                   ENDIF
@@ -1963 +1939 @@
                !! AXY (11/11/16): CMIP6 oxygen saturation 3D diagnostic
                IF ( med_diag%O2SAT3%dgsave ) THEN
-                  call oxy_sato( ztmp, zsal, f_o2sat3 )
+                  call oxy_sato( ztmp(ji,jj), zsal(ji,jj), f_o2sat3 )
                   o2sat3(ji, jj, jk) = f_o2sat3
                ENDIF
@@ -2049 +2025 @@
                !!
                !! non-diatoms
-          if (zphn.GT.rsmall) then
-                  fthetan = max(tiny(zchn), (zchn * xxi) / (zphn + tiny(zphn)))
-                  faln    = xaln * fthetan
+          if (zphn(ji,jj).GT.rsmall) then
+                  fthetan(ji,jj) = max(tiny(zchn(ji,jj)), (zchn(ji,jj) * xxi) / (zphn(ji,jj) + tiny(zphn(ji,jj))))
+                  faln(ji,jj)    = xaln * fthetan(ji,jj)
                else
-                  fthetan = 0.
-                  faln    = 0.
+                  fthetan(ji,jj) = 0.
+                  faln(ji,jj)    = 0.
                endif
                !!
                !! diatoms
-          if (zphd.GT.rsmall) then
-                  fthetad = max(tiny(zchd), (zchd * xxi) / (zphd + tiny(zphd)))
-                  fald    = xald * fthetad
+          if (zphd(ji,jj).GT.rsmall) then
+                  fthetad(ji,jj) = max(tiny(zchd(ji,jj)), (zchd(ji,jj) * xxi) / (zphd(ji,jj) + tiny(zphd(ji,jj))))
+                  fald(ji,jj)    = xald * fthetad(ji,jj)
                else
-                  fthetad = 0.
-                  fald    = 0.
+                  fthetad(ji,jj) = 0.
+                  fald(ji,jj)    = 0.
                endif
 
@@ -2070 +2046 @@
                if (idf.eq.1.AND.idfval.eq.1) then
                   IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'faln(',jk,') = ', faln
-                  IF (lwp) write (numout,*) 'fald(',jk,') = ', fald
+                  IF (lwp) write (numout,*) 'faln(',jk,') = ', faln(ji,jj)
+                  IF (lwp) write (numout,*) 'fald(',jk,') = ', fald(ji,jj)
                endif
 # endif
@@ -2092 +2068 @@
                !! AXY (16/07/09)
                !! temperature for new Eppley style phytoplankton growth
-               loc_T   = tsn(ji,jj,jk,jp_tem)
-               fun_T   = 1.066**(1.0 * loc_T)
+               fun_T(ji,jj)   = 1.066**(1.0 * tsn(ji,jj,jk,jp_tem))
                !! AXY (16/05/11): add in new Q10 (1.5, not 2.0) for
                !phytoplankton
                !!                 growth; remin. unaffected
-               fun_Q10 = jq10**((loc_T - 0.0) / 10.0)
+               fun_Q10(ji,jj) = jq10**((tsn(ji,jj,jk,jp_tem) - 0.0) / 10.0)
                if (jphy.eq.1) then
-                  xvpnT = xvpn * fun_T
-                  xvpdT = xvpd * fun_T
+                  xvpnT(ji,jj) = xvpn * fun_T(ji,jj)
+                  xvpdT(ji,jj) = xvpd * fun_T(ji,jj)
                elseif (jphy.eq.2) then
-                  xvpnT = xvpn * fun_Q10
-                  xvpdT = xvpd * fun_Q10
+                  xvpnT(ji,jj) = xvpn * fun_Q10(ji,jj)
+                  xvpdT(ji,jj) = xvpd * fun_Q10(ji,jj)
                else
-                  xvpnT = xvpn
-                  xvpdT = xvpd
+                  xvpnT(ji,jj) = xvpn
+                  xvpdT(ji,jj) = xvpd
                endif
                !!
                !! non-diatoms
-               fchn1   = (xvpnT * xvpnT) + (faln * faln * xpar(ji,jj,jk) * xpar(ji,jj,jk))
-               if (fchn1.GT.rsmall) then
-                  fchn    = xvpnT / (sqrt(fchn1) + tiny(fchn1))
+               fchn1(ji,jj)   = (xvpnT(ji,jj) * xvpnT(ji,jj)) + (faln(ji,jj) * faln(ji,jj) * xpar(ji,jj,jk) * xpar(ji,jj,jk))
+               if (fchn1(ji,jj).GT.rsmall) then
+                  fchn(ji,jj)    = xvpnT(ji,jj) / (sqrt(fchn1(ji,jj)) + tiny(fchn1(ji,jj)))
                else
-                  fchn    = 0.
+                  fchn(ji,jj)    = 0.
                endif
-               fjln    = fchn * faln * xpar(ji,jj,jk) !! non-diatom J term
-               fjlim_pn = fjln / xvpnT
+               fjln(ji,jj)    = fchn(ji,jj) * faln(ji,jj) * xpar(ji,jj,jk) !! non-diatom J term
+               fjlim_pn(ji,jj) = fjln(ji,jj) / xvpnT(ji,jj)
                !!
                !! diatoms
-               fchd1   = (xvpdT * xvpdT) + (fald * fald * xpar(ji,jj,jk) * xpar(ji,jj,jk))
-               if (fchd1.GT.rsmall) then
-                  fchd    = xvpdT / (sqrt(fchd1) + tiny(fchd1))
+               fchd1(ji,jj)   = (xvpdT(ji,jj) * xvpdT(ji,jj)) + (fald(ji,jj) * fald(ji,jj) * xpar(ji,jj,jk) * xpar(ji,jj,jk))
+               if (fchd1(ji,jj).GT.rsmall) then
+                  fchd(ji,jj)    = xvpdT(ji,jj) / (sqrt(fchd1(ji,jj)) + tiny(fchd1(ji,jj)))
                else
-                  fchd    = 0.
+                  fchd(ji,jj)    = 0.
                endif
-               fjld    = fchd * fald * xpar(ji,jj,jk) !! diatom J term
-               fjlim_pd = fjld / xvpdT
+               fjld(ji,jj)    = fchd(ji,jj) * fald(ji,jj) * xpar(ji,jj,jk) !! diatom J term
+               fjlim_pd(ji,jj) = fjld(ji,jj) / xvpdT(ji,jj)
       
 # if defined key_debug_medusa
@@ -2133 +2108 @@
                if (idf.eq.1.AND.idfval.eq.1) then
                   IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'fchn(',jk,') = ', fchn
-                  IF (lwp) write (numout,*) 'fchd(',jk,') = ', fchd
-                  IF (lwp) write (numout,*) 'fjln(',jk,') = ', fjln
-                  IF (lwp) write (numout,*) 'fjld(',jk,') = ', fjld
+                  IF (lwp) write (numout,*) 'fchn(',jk,') = ', fchn(ji,jj)
+                  IF (lwp) write (numout,*) 'fchd(',jk,') = ', fchd(ji,jj)
+                  IF (lwp) write (numout,*) 'fjln(',jk,') = ', fjln(ji,jj)
+                  IF (lwp) write (numout,*) 'fjld(',jk,') = ', fjld(ji,jj)
                endif
 # endif
@@ -2145 +2120 @@
                !!
                !! non-diatoms (N, Fe)
-               fnln = zdin / (zdin + xnln) !! non-diatom Qn term
-               ffln = zfer / (zfer + xfln) !! non-diatom Qf term
+               fnln(ji,jj) = zdin(ji,jj) / (zdin(ji,jj) + xnln) !! non-diatom Qn term
+               ffln(ji,jj) = zfer(ji,jj) / (zfer(ji,jj) + xfln) !! non-diatom Qf term
                !!
                !! diatoms (N, Si, Fe)
-               fnld = zdin / (zdin + xnld) !! diatom Qn term
-               fsld = zsil / (zsil + xsld) !! diatom Qs term
-               ffld = zfer / (zfer + xfld) !! diatom Qf term
+               fnld(ji,jj) = zdin(ji,jj) / (zdin(ji,jj) + xnld) !! diatom Qn term
+               fsld(ji,jj) = zsil(ji,jj) / (zsil(ji,jj) + xsld) !! diatom Qs term
+               ffld(ji,jj) = zfer(ji,jj) / (zfer(ji,jj) + xfld) !! diatom Qf term
 
 # if defined key_debug_medusa
@@ -2157 +2132 @@
                if (idf.eq.1.AND.idfval.eq.1) then
                   IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'fnln(',jk,') = ', fnln
-                  IF (lwp) write (numout,*) 'fnld(',jk,') = ', fnld
-                  IF (lwp) write (numout,*) 'ffln(',jk,') = ', ffln
-                  IF (lwp) write (numout,*) 'ffld(',jk,') = ', ffld
-                  IF (lwp) write (numout,*) 'fsld(',jk,') = ', fsld
+                  IF (lwp) write (numout,*) 'fnln(',jk,') = ', fnln(ji,jj)
+                  IF (lwp) write (numout,*) 'fnld(',jk,') = ', fnld(ji,jj)
+                  IF (lwp) write (numout,*) 'ffln(',jk,') = ', ffln(ji,jj)
+                  IF (lwp) write (numout,*) 'ffld(',jk,') = ', ffld(ji,jj)
+                  IF (lwp) write (numout,*) 'fsld(',jk,') = ', fsld(ji,jj)
                endif
 # endif
@@ -2172 +2147 @@
                if (jliebig .eq. 0) then
                   !! multiplicative nutrient limitation
-                  fpnlim = fnln * ffln
+                  fpnlim(ji,jj) = fnln(ji,jj) * ffln(ji,jj)
                elseif (jliebig .eq. 1) then
                   !! Liebig Law (= most limiting) nutrient limitation
-                  fpnlim = min(fnln, ffln)
+                  fpnlim(ji,jj) = min(fnln(ji,jj), ffln(ji,jj))
                endif
-               fprn = fjln * fpnlim
+               fprn(ji,jj) = fjln(ji,jj) * fpnlim(ji,jj)
 
                !!----------------------------------------------------------------------
@@ -2191 +2166 @@
                if (jliebig .eq. 0) then
                   !! multiplicative nutrient limitation
-                  fpdlim = fnld * ffld
+                  fpdlim(ji,jj) = fnld(ji,jj) * ffld(ji,jj)
                elseif (jliebig .eq. 1) then
                   !! Liebig Law (= most limiting) nutrient limitation
-                  fpdlim = min(fnld, ffld)
+                  fpdlim(ji,jj) = min(fnld(ji,jj), ffld(ji,jj))
                endif
                !!
-          if (zphd.GT.rsmall .AND. zpds.GT.rsmall) then
+          if (zphd(ji,jj).GT.rsmall .AND. zpds(ji,jj).GT.rsmall) then
                   !! "intracellular" elemental ratios
-                  ! fsin  = zpds / (zphd + tiny(zphd))
-                  ! fnsi  = zphd / (zpds + tiny(zpds))
-                  fsin = 0.0
-                  IF( zphd .GT. rsmall) fsin  = zpds / zphd
-                  fnsi = 0.0
-                  IF( zpds .GT. rsmall) fnsi  = zphd / zpds
+                  ! fsin(ji,jj)  = zpds(ji,jj) / (zphd(ji,jj) + tiny(zphd(ji,jj)))
+                  ! fnsi(ji,jj)  = zphd(ji,jj) / (zpds(ji,jj) + tiny(zpds(ji,jj)))
+                  fsin(ji,jj) = 0.0
+                  IF( zphd(ji,jj) .GT. rsmall) fsin(ji,jj)  = zpds(ji,jj) / zphd(ji,jj)
+                  fnsi(ji,jj) = 0.0
+                  IF( zpds(ji,jj) .GT. rsmall) fnsi(ji,jj)  = zphd(ji,jj) / zpds(ji,jj)
                   !! AXY (23/02/10): these next variables derive from Mongin et al. (2003)
                   fsin1 = 3.0 * xsin0 !! = 0.6
@@ -2212 +2187 @@
                   !! conditionalities based on ratios
                   !! nitrogen (and iron and carbon)
-                  if (fsin.le.xsin0) then
-                     fprd  = 0.0
-                     fsld2 = 0.0
-                  elseif (fsin.lt.fsin1) then
-                     fprd  = xuif * ((fsin - xsin0) / (fsin + tiny(fsin))) * (fjld * fpdlim)
-                     fsld2 = xuif * ((fsin - xsin0) / (fsin + tiny(fsin)))
-                  elseif (fsin.ge.fsin1) then
-                     fprd  = (fjld * fpdlim)
-                     fsld2 = 1.0
+                  if (fsin(ji,jj).le.xsin0) then
+                     fprd(ji,jj)  = 0.0
+                     fsld2(ji,jj) = 0.0
+                  elseif (fsin(ji,jj).lt.fsin1) then
+                     fprd(ji,jj)  = xuif * ((fsin(ji,jj) - xsin0) / (fsin(ji,jj) + tiny(fsin(ji,jj)))) * (fjld(ji,jj) * fpdlim(ji,jj))
+                     fsld2(ji,jj) = xuif * ((fsin(ji,jj) - xsin0) / (fsin(ji,jj) + tiny(fsin(ji,jj))))
+                  elseif (fsin(ji,jj).ge.fsin1) then
+                     fprd(ji,jj)  = (fjld(ji,jj) * fpdlim(ji,jj))
+                     fsld2(ji,jj) = 1.0
                   endif
                   !!
                   !! silicon
-                  if (fsin.lt.fnsi1) then
-                     fprds = (fjld * fsld)
-                  elseif (fsin.lt.fnsi2) then
-                     fprds = xuif * ((fnsi - xnsi0) / (fnsi + tiny(fnsi))) * (fjld * fsld)
+                  if (fsin(ji,jj).lt.fnsi1) then
+                     fprds(ji,jj) = (fjld(ji,jj) * fsld(ji,jj))
+                  elseif (fsin(ji,jj).lt.fnsi2) then
+                     fprds(ji,jj) = xuif * ((fnsi(ji,jj) - xnsi0) / (fnsi(ji,jj) + tiny(fnsi(ji,jj)))) * (fjld(ji,jj) * fsld(ji,jj))
                   else
-                     fprds = 0.0
+                     fprds(ji,jj) = 0.0
                   endif     
                else
-                  fsin  = 0.0
-                  fnsi  = 0.0
-                  fprd  = 0.0
-                  fsld2 = 0.0
-                  fprds = 0.0
+                  fsin(ji,jj)  = 0.0
+                  fnsi(ji,jj)  = 0.0
+                  fprd(ji,jj)  = 0.0
+                  fsld2(ji,jj) = 0.0
+                  fprds(ji,jj) = 0.0
                endif
 
@@ -2243 +2218 @@
                if (idf.eq.1.AND.idfval.eq.1) then
                   IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'fsin(',jk,')   = ', fsin
-                  IF (lwp) write (numout,*) 'fnsi(',jk,')   = ', fnsi
-                  IF (lwp) write (numout,*) 'fsld2(',jk,')  = ', fsld2
-                  IF (lwp) write (numout,*) 'fprn(',jk,')   = ', fprn
-                  IF (lwp) write (numout,*) 'fprd(',jk,')   = ', fprd
-                  IF (lwp) write (numout,*) 'fprds(',jk,')  = ', fprds
+                  IF (lwp) write (numout,*) 'fsin(',jk,')   = ', fsin(ji,jj)
+                  IF (lwp) write (numout,*) 'fnsi(',jk,')   = ', fnsi(ji,jj)
+                  IF (lwp) write (numout,*) 'fsld2(',jk,')  = ', fsld2(ji,jj)
+                  IF (lwp) write (numout,*) 'fprn(',jk,')   = ', fprn(ji,jj)
+                  IF (lwp) write (numout,*) 'fprd(',jk,')   = ', fprd(ji,jj)
+                  IF (lwp) write (numout,*) 'fprds(',jk,')  = ', fprds(ji,jj)
                endif
 # endif
@@ -2261 +2236 @@
                !!----------------------------------------------------------------------
                !!
-               if (fdep1.le.hmld(ji,jj)) then
+               if (fdep1(ji,jj).le.hmld(ji,jj)) then
                   !! this level is entirely in the mixed layer
                   fq0 = 1.0
-               elseif (fdep.ge.hmld(ji,jj)) then
+               elseif (fsdepw(ji,jj,jk).ge.hmld(ji,jj)) then
                   !! this level is entirely below the mixed layer
                   fq0 = 0.0
                else
                   !! this level straddles the mixed layer
-                  fq0 = (hmld(ji,jj) - fdep) / fthk
+                  fq0 = (hmld(ji,jj) - fsdepw(ji,jj,jk)) / fse3t(ji,jj,jk)
                endif
                !!
-               fprn_ml(ji,jj) = fprn_ml(ji,jj) + (fprn * zphn * fthk * fq0)
-               fprd_ml(ji,jj) = fprd_ml(ji,jj) + (fprd * zphd * fthk * fq0)
+               fprn_ml(ji,jj) = fprn_ml(ji,jj) + (fprn(ji,jj) * zphn(ji,jj) * fse3t(ji,jj,jk) * fq0)
+               fprd_ml(ji,jj) = fprd_ml(ji,jj) + (fprd(ji,jj) * zphd(ji,jj) * fse3t(ji,jj,jk) * fq0)
                
                !!----------------------------------------------------------------------
                !! Vertical Integral --
                !!----------------------------------------------------------------------
-               ftot_pn(ji,jj)  = ftot_pn(ji,jj)  + (zphn * fthk)   !! vertical integral non-diatom phytoplankton
-               ftot_pd(ji,jj)  = ftot_pd(ji,jj)  + (zphd * fthk)   !! vertical integral diatom phytoplankton
-               ftot_zmi(ji,jj) = ftot_zmi(ji,jj) + (zzmi * fthk)   !! vertical integral microzooplankton
-               ftot_zme(ji,jj) = ftot_zme(ji,jj) + (zzme * fthk)   !! vertical integral mesozooplankton
-               ftot_det(ji,jj) = ftot_det(ji,jj) + (zdet * fthk)   !! vertical integral slow detritus, nitrogen
-               ftot_dtc(ji,jj) = ftot_dtc(ji,jj) + (zdtc * fthk)   !! vertical integral slow detritus, carbon
+               ftot_pn(ji,jj)  = ftot_pn(ji,jj)  + (zphn(ji,jj) * fse3t(ji,jj,jk))   !! vertical integral non-diatom phytoplankton
+               ftot_pd(ji,jj)  = ftot_pd(ji,jj)  + (zphd(ji,jj) * fse3t(ji,jj,jk))   !! vertical integral diatom phytoplankton
+               ftot_zmi(ji,jj) = ftot_zmi(ji,jj) + (zzmi(ji,jj) * fse3t(ji,jj,jk))   !! vertical integral microzooplankton
+               ftot_zme(ji,jj) = ftot_zme(ji,jj) + (zzme(ji,jj) * fse3t(ji,jj,jk))   !! vertical integral mesozooplankton
+               ftot_det(ji,jj) = ftot_det(ji,jj) + (zdet(ji,jj) * fse3t(ji,jj,jk))   !! vertical integral slow detritus, nitrogen
+               ftot_dtc(ji,jj) = ftot_dtc(ji,jj) + (zdtc(ji,jj) * fse3t(ji,jj,jk))   !! vertical integral slow detritus, carbon
                
                !!----------------------------------------------------------------------
@@ -2289 +2264 @@
                !!----------------------------------------------------------------------
                !!
-               !! frn = (xthetam / fthetan) * (fprn / (fthetan * xpar(ji,jj,jk)))
-               !! frd = (xthetam / fthetad) * (fprd / (fthetad * xpar(ji,jj,jk)))
-               frn = (xthetam * fchn * fnln * ffln       ) / (fthetan + tiny(fthetan))
+               !! frn(ji,jj) = (xthetam / fthetan(ji,jj)) * (fprn(ji,jj) / (fthetan(ji,jj) * xpar(ji,jj,jk)))
+               !! frd(ji,jj) = (xthetam / fthetad(ji,jj)) * (fprd(ji,jj) / (fthetad(ji,jj) * xpar(ji,jj,jk)))
+               frn(ji,jj) = (xthetam * fchn(ji,jj) * fnln(ji,jj) * ffln(ji,jj)       ) / (fthetan(ji,jj) + tiny(fthetan(ji,jj)))
                !! AXY (12/05/09): there's potentially a problem here; fsld, silicic acid 
                !!   limitation, is used in the following line to regulate chlorophyll 
@@ -2297 +2272 @@
                !!   of biomass growth; the Mongin term term used in growth is more complex
                !!   than the simple multiplicative function used below
-               !! frd = (xthetam * fchd * fnld * ffld * fsld) / (fthetad + tiny(fthetad))
+               !! frd(ji,jj) = (xthetam * fchd(ji,jj) * fnld(ji,jj) * ffld(ji,jj) * fsld(ji,jj)) / (fthetad(ji,jj) + tiny(fthetad(ji,jj)))
                !! AXY (12/05/09): this replacement line uses the new variable, fsld2, to
                !!   regulate chlorophyll growth
-               frd = (xthetamd * fchd * fnld * ffld * fsld2) / (fthetad + tiny(fthetad))
+               frd(ji,jj) = (xthetamd * fchd(ji,jj) * fnld(ji,jj) * ffld(ji,jj) * fsld2(ji,jj)) / (fthetad(ji,jj) + tiny(fthetad(ji,jj)))
 
 # if defined key_debug_medusa
@@ -2306 +2281 @@
                if (idf.eq.1.AND.idfval.eq.1) then
                   IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'fthetan(',jk,') = ', fthetan
-                  IF (lwp) write (numout,*) 'fthetad(',jk,') = ', fthetad
-                  IF (lwp) write (numout,*) 'frn(',jk,')     = ', frn
-                  IF (lwp) write (numout,*) 'frd(',jk,')     = ', frd
+                  IF (lwp) write (numout,*) 'fthetan(',jk,') = ', fthetan(ji,jj)
+                  IF (lwp) write (numout,*) 'fthetad(',jk,') = ', fthetad(ji,jj)
+                  IF (lwp) write (numout,*) 'frn(',jk,')     = ', frn(ji,jj)
+                  IF (lwp) write (numout,*) 'frd(',jk,')     = ', frd(ji,jj)
                endif
 # endif
@@ -2329 +2304 @@
                !!----------------------------------------------------------------------
                !!
-               fmi1    = (xkmi * xkmi) + (xpmipn * zphn * zphn) + (xpmid * zdet * zdet)
-               fmi     = xgmi * zzmi / fmi1
-               fgmipn  = fmi * xpmipn * zphn * zphn   !! grazing on non-diatoms
-               fgmid   = fmi * xpmid  * zdet * zdet   !! grazing on detrital nitrogen
+               fmi1(ji,jj)    = (xkmi * xkmi) + (xpmipn * zphn(ji,jj) * zphn(ji,jj)) + (xpmid * zdet(ji,jj) * zdet(ji,jj))
+               fmi(ji,jj)     = xgmi * zzmi(ji,jj) / fmi1(ji,jj)
+               fgmipn(ji,jj)  = fmi(ji,jj) * xpmipn * zphn(ji,jj) * zphn(ji,jj)   !! grazing on non-diatoms
+               fgmid(ji,jj)   = fmi(ji,jj) * xpmid  * zdet(ji,jj) * zdet(ji,jj)   !! grazing on detrital nitrogen
 # if defined key_roam
-               fgmidc  = rsmall !acc
-               IF ( zdet .GT. rsmall ) fgmidc  = (zdtc / (zdet + tiny(zdet))) * fgmid  !! grazing on detrital carbon
+               fgmidc(ji,jj)  = rsmall !acc
+               IF ( zdet(ji,jj) .GT. rsmall ) fgmidc(ji,jj)  = (zdtc(ji,jj) / (zdet(ji,jj) + tiny(zdet(ji,jj)))) * fgmid(ji,jj)  !! grazing on detrital carbon
 # else
                !! AXY (26/11/08): implicit detrital carbon change
-               fgmidc  = xthetad * fgmid              !! grazing on detrital carbon
+               fgmidc(ji,jj)  = xthetad * fgmid(ji,jj)              !! grazing on detrital carbon
 # endif
                !!
                !! which translates to these incoming N and C fluxes
-               finmi   = (1.0 - xphi) * (fgmipn + fgmid)
-               ficmi   = (1.0 - xphi) * ((xthetapn * fgmipn) + fgmidc)
+               finmi(ji,jj)   = (1.0 - xphi) * (fgmipn(ji,jj) + fgmid(ji,jj))
+               ficmi(ji,jj)   = (1.0 - xphi) * ((xthetapn * fgmipn(ji,jj)) + fgmidc(ji,jj))
                !!
                !! the ideal food C:N ratio for microzooplankton
                !! xbetan = 0.77; xthetaz = 5.625; xbetac = 0.64; xkc = 0.80
-               fstarmi = (xbetan * xthetazmi) / (xbetac * xkc)
+               fstarmi(ji,jj) = (xbetan * xthetazmi) / (xbetac * xkc)
                !!
                !! process these to determine proportioning of grazed N and C
                !! (since there is no explicit consideration of respiration,
                !! only growth and excretion are calculated here)
-               fmith   = (ficmi / (finmi + tiny(finmi)))
-               if (fmith.ge.fstarmi) then
-                  fmigrow = xbetan * finmi
-                  fmiexcr = 0.0
+               fmith(ji,jj)   = (ficmi(ji,jj) / (finmi(ji,jj) + tiny(finmi(ji,jj))))
+               if (fmith(ji,jj).ge.fstarmi(ji,jj)) then
+                  fmigrow(ji,jj) = xbetan * finmi(ji,jj)
+                  fmiexcr(ji,jj) = 0.0
                else
-                  fmigrow = (xbetac * xkc * ficmi) / xthetazmi
-                  fmiexcr = ficmi * ((xbetan / (fmith + tiny(fmith))) - ((xbetac * xkc) / xthetazmi))
+                  fmigrow(ji,jj) = (xbetac * xkc * ficmi(ji,jj)) / xthetazmi
+                  fmiexcr(ji,jj) = ficmi(ji,jj) * ((xbetan / (fmith(ji,jj) + tiny(fmith(ji,jj)))) - ((xbetac * xkc) / xthetazmi))
                endif
 # if defined key_roam
-               fmiresp = (xbetac * ficmi) - (xthetazmi * fmigrow)
+               fmiresp(ji,jj) = (xbetac * ficmi(ji,jj)) - (xthetazmi * fmigrow(ji,jj))
 # endif
 
@@ -2368 +2343 @@
                if (idf.eq.1.AND.idfval.eq.1) then
                   IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'fmi1(',jk,')    = ', fmi1
-                  IF (lwp) write (numout,*) 'fmi(',jk,')     = ', fmi
-                  IF (lwp) write (numout,*) 'fgmipn(',jk,')  = ', fgmipn
-                  IF (lwp) write (numout,*) 'fgmid(',jk,')   = ', fgmid
-                  IF (lwp) write (numout,*) 'fgmidc(',jk,')  = ', fgmidc
-                  IF (lwp) write (numout,*) 'finmi(',jk,')   = ', finmi
-                  IF (lwp) write (numout,*) 'ficmi(',jk,')   = ', ficmi
-                  IF (lwp) write (numout,*) 'fstarmi(',jk,') = ', fstarmi
-                  IF (lwp) write (numout,*) 'fmith(',jk,')   = ', fmith
-                  IF (lwp) write (numout,*) 'fmigrow(',jk,') = ', fmigrow
-                  IF (lwp) write (numout,*) 'fmiexcr(',jk,') = ', fmiexcr
+                  IF (lwp) write (numout,*) 'fmi1(',jk,')    = ', fmi1(ji,jj)
+                  IF (lwp) write (numout,*) 'fmi(',jk,')     = ', fmi(ji,jj)
+                  IF (lwp) write (numout,*) 'fgmipn(',jk,')  = ', fgmipn(ji,jj)
+                  IF (lwp) write (numout,*) 'fgmid(',jk,')   = ', fgmid(ji,jj)
+                  IF (lwp) write (numout,*) 'fgmidc(',jk,')  = ', fgmidc(ji,jj)
+                  IF (lwp) write (numout,*) 'finmi(',jk,')   = ', finmi(ji,jj)
+                  IF (lwp) write (numout,*) 'ficmi(',jk,')   = ', ficmi(ji,jj)
+                  IF (lwp) write (numout,*) 'fstarmi(',jk,') = ', fstarmi(ji,jj)
+                  IF (lwp) write (numout,*) 'fmith(',jk,')   = ', fmith(ji,jj)
+                  IF (lwp) write (numout,*) 'fmigrow(',jk,') = ', fmigrow(ji,jj)
+                  IF (lwp) write (numout,*) 'fmiexcr(',jk,') = ', fmiexcr(ji,jj)
 #  if defined key_roam
-                  IF (lwp) write (numout,*) 'fmiresp(',jk,') = ', fmiresp
+                  IF (lwp) write (numout,*) 'fmiresp(',jk,') = ', fmiresp(ji,jj)
 #  endif
                endif
@@ -2389 +2364 @@
                !!----------------------------------------------------------------------
                !!
-               fme1    = (xkme * xkme) + (xpmepn * zphn * zphn) + (xpmepd * zphd * zphd) + & 
-                         (xpmezmi * zzmi * zzmi) + (xpmed * zdet * zdet)
-               fme     = xgme * zzme / fme1
-               fgmepn  = fme * xpmepn  * zphn * zphn  !! grazing on non-diatoms
-               fgmepd  = fme * xpmepd  * zphd * zphd  !! grazing on diatoms
-               fgmepds = fsin * fgmepd                !! grazing on diatom silicon
-               fgmezmi = fme * xpmezmi * zzmi * zzmi  !! grazing on microzooplankton
-               fgmed   = fme * xpmed   * zdet * zdet  !! grazing on detrital nitrogen
+               fme1(ji,jj)    = (xkme * xkme) + (xpmepn * zphn(ji,jj) * zphn(ji,jj)) + (xpmepd * zphd(ji,jj) * zphd(ji,jj)) + & 
+                         (xpmezmi * zzmi(ji,jj) * zzmi(ji,jj)) + (xpmed * zdet(ji,jj) * zdet(ji,jj))
+               fme(ji,jj)     = xgme * zzme(ji,jj) / fme1(ji,jj)
+               fgmepn(ji,jj)  = fme(ji,jj) * xpmepn  * zphn(ji,jj) * zphn(ji,jj)  !! grazing on non-diatoms
+               fgmepd(ji,jj)  = fme(ji,jj) * xpmepd  * zphd(ji,jj) * zphd(ji,jj)  !! grazing on diatoms
+               fgmepds(ji,jj) = fsin(ji,jj) * fgmepd(ji,jj)                !! grazing on diatom silicon
+               fgmezmi(ji,jj) = fme(ji,jj) * xpmezmi * zzmi(ji,jj) * zzmi(ji,jj)  !! grazing on microzooplankton
+               fgmed(ji,jj)   = fme(ji,jj) * xpmed   * zdet(ji,jj) * zdet(ji,jj)  !! grazing on detrital nitrogen
 # if defined key_roam
-               fgmedc  = rsmall !acc
-               IF ( zdet .GT. rsmall ) fgmedc  = (zdtc / (zdet + tiny(zdet))) * fgmed  !! grazing on detrital carbon
+               fgmedc(ji,jj)  = rsmall !acc
+               IF ( zdet(ji,jj) .GT. rsmall ) fgmedc(ji,jj)  = (zdtc(ji,jj) / (zdet(ji,jj) + tiny(zdet(ji,jj)))) * fgmed(ji,jj)  !! grazing on detrital carbon
 # else
                !! AXY (26/11/08): implicit detrital carbon change
-               fgmedc  = xthetad * fgmed              !! grazing on detrital carbon
+               fgmedc(ji,jj)  = xthetad * fgmed(ji,jj)              !! grazing on detrital carbon
 # endif
                !!
                !! which translates to these incoming N and C fluxes
-               finme   = (1.0 - xphi) * (fgmepn + fgmepd + fgmezmi + fgmed)
-               ficme   = (1.0 - xphi) * ((xthetapn * fgmepn) + (xthetapd * fgmepd) + &
-                        (xthetazmi * fgmezmi) + fgmedc)
+               finme(ji,jj)   = (1.0 - xphi) * (fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj))
+               ficme(ji,jj)   = (1.0 - xphi) * ((xthetapn * fgmepn(ji,jj)) + (xthetapd * fgmepd(ji,jj)) + &
+                        (xthetazmi * fgmezmi(ji,jj)) + fgmedc(ji,jj))
                !!
                !! the ideal food C:N ratio for mesozooplankton
                !! xbetan = 0.77; xthetaz = 5.625; xbetac = 0.64; xkc = 0.80
-               fstarme = (xbetan * xthetazme) / (xbetac * xkc)
+               fstarme(ji,jj) = (xbetan * xthetazme) / (xbetac * xkc)
                !!
                !! process these to determine proportioning of grazed N and C
                !! (since there is no explicit consideration of respiration,
                !! only growth and excretion are calculated here)
-               fmeth   = (ficme / (finme + tiny(finme)))
-               if (fmeth.ge.fstarme) then
-                  fmegrow = xbetan * finme
-                  fmeexcr = 0.0
+               fmeth(ji,jj)   = (ficme(ji,jj) / (finme(ji,jj) + tiny(finme(ji,jj))))
+               if (fmeth(ji,jj).ge.fstarme(ji,jj)) then
+                  fmegrow(ji,jj) = xbetan * finme(ji,jj)
+                  fmeexcr(ji,jj) = 0.0
                else
-                  fmegrow = (xbetac * xkc * ficme) / xthetazme
-                  fmeexcr = ficme * ((xbetan / (fmeth + tiny(fmeth))) - ((xbetac * xkc) / xthetazme))
+                  fmegrow(ji,jj) = (xbetac * xkc * ficme(ji,jj)) / xthetazme
+                  fmeexcr(ji,jj) = ficme(ji,jj) * ((xbetan / (fmeth(ji,jj) + tiny(fmeth(ji,jj)))) - ((xbetac * xkc) / xthetazme))
                endif
 # if defined key_roam
-               fmeresp = (xbetac * ficme) - (xthetazme * fmegrow)
+               fmeresp(ji,jj) = (xbetac * ficme(ji,jj)) - (xthetazme * fmegrow(ji,jj))
 # endif
 
@@ -2433 +2408 @@
                if (idf.eq.1.AND.idfval.eq.1) then
                   IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'fme1(',jk,')    = ', fme1
-                  IF (lwp) write (numout,*) 'fme(',jk,')     = ', fme
-                  IF (lwp) write (numout,*) 'fgmepn(',jk,')  = ', fgmepn
-                  IF (lwp) write (numout,*) 'fgmepd(',jk,')  = ', fgmepd
-                  IF (lwp) write (numout,*) 'fgmepds(',jk,') = ', fgmepds
-                  IF (lwp) write (numout,*) 'fgmezmi(',jk,') = ', fgmezmi
-                  IF (lwp) write (numout,*) 'fgmed(',jk,')   = ', fgmed
-                  IF (lwp) write (numout,*) 'fgmedc(',jk,')  = ', fgmedc
-                  IF (lwp) write (numout,*) 'finme(',jk,')   = ', finme
-                  IF (lwp) write (numout,*) 'ficme(',jk,')   = ', ficme
-                  IF (lwp) write (numout,*) 'fstarme(',jk,') = ', fstarme
-                  IF (lwp) write (numout,*) 'fmeth(',jk,')   = ', fmeth
-                  IF (lwp) write (numout,*) 'fmegrow(',jk,') = ', fmegrow
-                  IF (lwp) write (numout,*) 'fmeexcr(',jk,') = ', fmeexcr
+                  IF (lwp) write (numout,*) 'fme1(',jk,')    = ', fme1(ji,jj)
+                  IF (lwp) write (numout,*) 'fme(',jk,')     = ', fme(ji,jj)
+                  IF (lwp) write (numout,*) 'fgmepn(',jk,')  = ', fgmepn(ji,jj)
+                  IF (lwp) write (numout,*) 'fgmepd(',jk,')  = ', fgmepd(ji,jj)
+                  IF (lwp) write (numout,*) 'fgmepds(',jk,') = ', fgmepds(ji,jj)
+                  IF (lwp) write (numout,*) 'fgmezmi(',jk,') = ', fgmezmi(ji,jj)
+                  IF (lwp) write (numout,*) 'fgmed(',jk,')   = ', fgmed(ji,jj)
+                  IF (lwp) write (numout,*) 'fgmedc(',jk,')  = ', fgmedc(ji,jj)
+                  IF (lwp) write (numout,*) 'finme(',jk,')   = ', finme(ji,jj)
+                  IF (lwp) write (numout,*) 'ficme(',jk,')   = ', ficme(ji,jj)
+                  IF (lwp) write (numout,*) 'fstarme(',jk,') = ', fstarme(ji,jj)
+                  IF (lwp) write (numout,*) 'fmeth(',jk,')   = ', fmeth(ji,jj)
+                  IF (lwp) write (numout,*) 'fmegrow(',jk,') = ', fmegrow(ji,jj)
+                  IF (lwp) write (numout,*) 'fmeexcr(',jk,') = ', fmeexcr(ji,jj)
 #  if defined key_roam
-                  IF (lwp) write (numout,*) 'fmeresp(',jk,') = ', fmeresp
+                  IF (lwp) write (numout,*) 'fmeresp(',jk,') = ', fmeresp(ji,jj)
 #  endif
                endif
 # endif
 
-               fzmi_i(ji,jj)  = fzmi_i(ji,jj)  + fthk * (  &
-                  fgmipn + fgmid )
-               fzmi_o(ji,jj)  = fzmi_o(ji,jj)  + fthk * (  &
-                  fmigrow + (xphi * (fgmipn + fgmid)) + fmiexcr + ((1.0 - xbetan) * finmi) )
-               fzme_i(ji,jj)  = fzme_i(ji,jj)  + fthk * (  &
-                  fgmepn + fgmepd + fgmezmi + fgmed )
-               fzme_o(ji,jj)  = fzme_o(ji,jj)  + fthk * (  &
-                  fmegrow + (xphi * (fgmepn + fgmepd + fgmezmi + fgmed)) + fmeexcr + ((1.0 - xbetan) * finme) )
+               fzmi_i(ji,jj)  = fzmi_i(ji,jj)  + fse3t(ji,jj,jk) * (  &
+                  fgmipn(ji,jj) + fgmid(ji,jj) )
+               fzmi_o(ji,jj)  = fzmi_o(ji,jj)  + fse3t(ji,jj,jk) * (  &
+                  fmigrow(ji,jj) + (xphi * (fgmipn(ji,jj) + fgmid(ji,jj))) + fmiexcr(ji,jj) + ((1.0 - xbetan) * finmi(ji,jj)) )
+               fzme_i(ji,jj)  = fzme_i(ji,jj)  + fse3t(ji,jj,jk) * (  &
+                  fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj) )
+               fzme_o(ji,jj)  = fzme_o(ji,jj)  + fse3t(ji,jj,jk) * (  &
+                  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)) )
 
                !!----------------------------------------------------------------------
@@ -2467 +2442 @@
                !!----------------------------------------------------------------------
                !!
-               fdpn2  = xmetapn  * zphn
-               fdpd2  = xmetapd  * zphd
-               fdpds2 = xmetapd  * zpds
-               fdzmi2 = xmetazmi * zzmi
-               fdzme2 = xmetazme * zzme
+               fdpn2(ji,jj)  = xmetapn  * zphn(ji,jj)
+               fdpd2(ji,jj)  = xmetapd  * zphd(ji,jj)
+               fdpds2(ji,jj) = xmetapd  * zpds(ji,jj)
+               fdzmi2(ji,jj) = xmetazmi * zzmi(ji,jj)
+               fdzme2(ji,jj) = xmetazme * zzme(ji,jj)
 
                !!----------------------------------------------------------------------
@@ -2480 +2455 @@
                !!
                !! non-diatom phytoplankton
-               if (jmpn.eq.1) fdpn = xmpn * zphn               !! linear
-               if (jmpn.eq.2) fdpn = xmpn * zphn * zphn        !! quadratic
-               if (jmpn.eq.3) fdpn = xmpn * zphn * &           !! hyperbolic
-                  (zphn / (xkphn + zphn))
-               if (jmpn.eq.4) fdpn = xmpn * zphn * &           !! sigmoid
-                  ((zphn * zphn) / (xkphn + (zphn * zphn)))
+               if (jmpn.eq.1) fdpn(ji,jj) = xmpn * zphn(ji,jj)               !! linear
+               if (jmpn.eq.2) fdpn(ji,jj) = xmpn * zphn(ji,jj) * zphn(ji,jj)        !! quadratic
+               if (jmpn.eq.3) fdpn(ji,jj) = xmpn * zphn(ji,jj) * &           !! hyperbolic
+                  (zphn(ji,jj) / (xkphn + zphn(ji,jj)))
+               if (jmpn.eq.4) fdpn(ji,jj) = xmpn * zphn(ji,jj) * &           !! sigmoid
+                  ((zphn(ji,jj) * zphn(ji,jj)) / (xkphn + (zphn(ji,jj) * zphn(ji,jj))))
                !!
                !! diatom phytoplankton
-               if (jmpd.eq.1) fdpd = xmpd * zphd               !! linear
-               if (jmpd.eq.2) fdpd = xmpd * zphd * zphd        !! quadratic
-               if (jmpd.eq.3) fdpd = xmpd * zphd * &           !! hyperbolic
-                  (zphd / (xkphd + zphd))
-               if (jmpd.eq.4) fdpd = xmpd * zphd * &           !! sigmoid
-                  ((zphd * zphd) / (xkphd + (zphd * zphd)))
-               fdpds = fdpd * fsin
+               if (jmpd.eq.1) fdpd(ji,jj) = xmpd * zphd(ji,jj)               !! linear
+               if (jmpd.eq.2) fdpd(ji,jj) = xmpd * zphd(ji,jj) * zphd(ji,jj)        !! quadratic
+               if (jmpd.eq.3) fdpd(ji,jj) = xmpd * zphd(ji,jj) * &           !! hyperbolic
+                  (zphd(ji,jj) / (xkphd + zphd(ji,jj)))
+               if (jmpd.eq.4) fdpd(ji,jj) = xmpd * zphd(ji,jj) * &           !! sigmoid
+                  ((zphd(ji,jj) * zphd(ji,jj)) / (xkphd + (zphd(ji,jj) * zphd(ji,jj))))
+               fdpds(ji,jj) = fdpd(ji,jj) * fsin(ji,jj)
                !!
                !! microzooplankton
-               if (jmzmi.eq.1) fdzmi = xmzmi * zzmi            !! linear
-               if (jmzmi.eq.2) fdzmi = xmzmi * zzmi * zzmi     !! quadratic
-               if (jmzmi.eq.3) fdzmi = xmzmi * zzmi * &        !! hyperbolic
-                  (zzmi / (xkzmi + zzmi))
-               if (jmzmi.eq.4) fdzmi = xmzmi * zzmi * &        !! sigmoid
-                  ((zzmi * zzmi) / (xkzmi + (zzmi * zzmi)))
+               if (jmzmi.eq.1) fdzmi(ji,jj) = xmzmi * zzmi(ji,jj)            !! linear
+               if (jmzmi.eq.2) fdzmi(ji,jj) = xmzmi * zzmi(ji,jj) * zzmi(ji,jj)     !! quadratic
+               if (jmzmi.eq.3) fdzmi(ji,jj) = xmzmi * zzmi(ji,jj) * &        !! hyperbolic
+                  (zzmi(ji,jj) / (xkzmi + zzmi(ji,jj)))
+               if (jmzmi.eq.4) fdzmi(ji,jj) = xmzmi * zzmi(ji,jj) * &        !! sigmoid
+                  ((zzmi(ji,jj) * zzmi(ji,jj)) / (xkzmi + (zzmi(ji,jj) * zzmi(ji,jj))))
                !!
                !! mesozooplankton
-               if (jmzme.eq.1) fdzme = xmzme * zzme            !! linear
-               if (jmzme.eq.2) fdzme = xmzme * zzme * zzme     !! quadratic
-               if (jmzme.eq.3) fdzme = xmzme * zzme * &        !! hyperbolic
-                  (zzme / (xkzme + zzme))
-               if (jmzme.eq.4) fdzme = xmzme * zzme * &        !! sigmoid
-                  ((zzme * zzme) / (xkzme + (zzme * zzme)))
+               if (jmzme.eq.1) fdzme(ji,jj) = xmzme * zzme(ji,jj)            !! linear
+               if (jmzme.eq.2) fdzme(ji,jj) = xmzme * zzme(ji,jj) * zzme(ji,jj)     !! quadratic
+               if (jmzme.eq.3) fdzme(ji,jj) = xmzme * zzme(ji,jj) * &        !! hyperbolic
+                  (zzme(ji,jj) / (xkzme + zzme(ji,jj)))
+               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))))
 
                !!----------------------------------------------------------------------
@@ -2519 +2494 @@
                if (jmd.eq.1) then
                   !! temperature-dependent
-                  fdd  = xmd  * fun_T * zdet
+                  fdd(ji,jj)  = xmd  * fun_T(ji,jj) * zdet(ji,jj)
 # if defined key_roam
-                  fddc = xmdc * fun_T * zdtc
+                  fddc(ji,jj) = xmdc * fun_T(ji,jj) * zdtc(ji,jj)
 # endif
                elseif (jmd.eq.2) then
                   !! AXY (16/05/13): add in Q10-based parameterisation (def in nmlst)
                   !! temperature-dependent
-                  fdd  = xmd  * fun_Q10 * zdet
+                  fdd(ji,jj)  = xmd  * fun_Q10(ji,jj) * zdet(ji,jj)
 #if defined key_roam
-                  fddc = xmdc * fun_Q10 * zdtc
+                  fddc(ji,jj) = xmdc * fun_Q10(ji,jj) * zdtc(ji,jj)
 #endif
                else
                   !! temperature-independent
-                  fdd  = xmd  * zdet
+                  fdd(ji,jj)  = xmd  * zdet(ji,jj)
 # if defined key_roam
-                  fddc = xmdc * zdtc
+                  fddc(ji,jj) = xmdc * zdtc(ji,jj)
 # endif
                endif
                !!
                !! AXY (22/07/09): accelerate detrital remineralisation in the bottom box
-               if ((jk.eq.jmbathy) .and. jsfd.eq.1) then
-                  fdd  = 1.0  * zdet
+               if ((jk.eq.mbathy(ji,jj)) .and. jsfd.eq.1) then
+                  fdd(ji,jj)  = 1.0  * zdet(ji,jj)
 # if defined key_roam
-                  fddc = 1.0  * zdtc
+                  fddc(ji,jj) = 1.0  * zdtc(ji,jj)
 # endif
                endif
@@ -2550 +2525 @@
                if (idf.eq.1.AND.idfval.eq.1) then
                   IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'fdpn2(',jk,') = ', fdpn2
-                  IF (lwp) write (numout,*) 'fdpd2(',jk,') = ', fdpd2
-                  IF (lwp) write (numout,*) 'fdpds2(',jk,')= ', fdpds2
-                  IF (lwp) write (numout,*) 'fdzmi2(',jk,')= ', fdzmi2
-                  IF (lwp) write (numout,*) 'fdzme2(',jk,')= ', fdzme2
-                  IF (lwp) write (numout,*) 'fdpn(',jk,')  = ', fdpn
-                  IF (lwp) write (numout,*) 'fdpd(',jk,')  = ', fdpd
-                  IF (lwp) write (numout,*) 'fdpds(',jk,') = ', fdpds
-                  IF (lwp) write (numout,*) 'fdzmi(',jk,') = ', fdzmi
-                  IF (lwp) write (numout,*) 'fdzme(',jk,') = ', fdzme
-                  IF (lwp) write (numout,*) 'fdd(',jk,')   = ', fdd
+                  IF (lwp) write (numout,*) 'fdpn2(',jk,') = ', fdpn2(ji,jj)
+                  IF (lwp) write (numout,*) 'fdpd2(',jk,') = ', fdpd2(ji,jj)
+                  IF (lwp) write (numout,*) 'fdpds2(',jk,')= ', fdpds2(ji,jj)
+                  IF (lwp) write (numout,*) 'fdzmi2(',jk,')= ', fdzmi2(ji,jj)
+                  IF (lwp) write (numout,*) 'fdzme2(',jk,')= ', fdzme2(ji,jj)
+                  IF (lwp) write (numout,*) 'fdpn(',jk,')  = ', fdpn(ji,jj)
+                  IF (lwp) write (numout,*) 'fdpd(',jk,')  = ', fdpd(ji,jj)
+                  IF (lwp) write (numout,*) 'fdpds(',jk,') = ', fdpds(ji,jj)
+                  IF (lwp) write (numout,*) 'fdzmi(',jk,') = ', fdzmi(ji,jj)
+                  IF (lwp) write (numout,*) 'fdzme(',jk,') = ', fdzme(ji,jj)
+                  IF (lwp) write (numout,*) 'fdd(',jk,')   = ', fdd(ji,jj)
 #  if defined key_roam
-                  IF (lwp) write (numout,*) 'fddc(',jk,')  = ', fddc
+                  IF (lwp) write (numout,*) 'fddc(',jk,')  = ', fddc(ji,jj)
 #  endif
                endif
@@ -2573 +2548 @@
                !!----------------------------------------------------------------------
                !!
-               if ((jk.eq.jmbathy) .and. jorgben.eq.1) then
+               if ((jk.eq.mbathy(ji,jj)) .and. jorgben.eq.1) then
                   !! this is the BOTTOM OCEAN BOX -> into the benthic pool!
                   !!
-                  f_sbenin_n(ji,jj)  = (zdet * vsed * 86400.)
-                  f_sbenin_fe(ji,jj) = (zdet * vsed * 86400. * xrfn)
+                  f_sbenin_n(ji,jj)  = (zdet(ji,jj) * vsed * 86400.)
+                  f_sbenin_fe(ji,jj) = (zdet(ji,jj) * vsed * 86400. * xrfn)
 # if defined key_roam
-                  f_sbenin_c(ji,jj)  = (zdtc * vsed * 86400.)
+                  f_sbenin_c(ji,jj)  = (zdtc(ji,jj) * vsed * 86400.)
 # else
-                  f_sbenin_c(ji,jj)  = (zdet * vsed * 86400. * xthetad)
+                  f_sbenin_c(ji,jj)  = (zdet(ji,jj) * vsed * 86400. * xthetad)
 # endif
                endif
@@ -2594 +2569 @@
                !!
                !! total iron concentration (mmol Fe / m3 -> umol Fe / m3)
-               xFeT        = zfer * 1.e3
+               xFeT        = zfer(ji,jj) * 1.e3
                !!
                !! calculate fractionation (based on Diat-HadOCC; in turn based on Parekh et al., 2004)
@@ -2608 +2583 @@
                !! "free" iron concentration (and convert to mmol Fe / m3)
                xFeF        = (xFeT - xFeL) * 1.e-3
-               xFree(ji,jj)= xFeF / (zfer + tiny(zfer))
+               xFree(ji,jj)= xFeF / (zfer(ji,jj) + tiny(zfer(ji,jj)))
                !!
                !! scavenging of iron (multiple schemes); I'm only really happy with the 
@@ -2621 +2596 @@
                   !!----------------------------------------------------------------------
                   !!
-                  ffescav     = xk_sc_Fe * xFeF                     ! = mmol/m3/d
+                  ffescav(ji,jj)     = xk_sc_Fe * xFeF                     ! = mmol/m3/d
                   !!
                   !!----------------------------------------------------------------------
@@ -2652 +2627 @@
                   !! scale of 1 day; seems reasonable to me ... (famous last words)
                   !!
-                  ffescav     = ffescav + fdeltaFe                  ! = mmol/m3/d
+                  ffescav(ji,jj)     = ffescav(ji,jj) + fdeltaFe                  ! = mmol/m3/d
                   !!
 # if defined key_deep_fe_fix
@@ -2670 +2645 @@
                   !! justifying this in a paper though!
                   !!
-                  if ((fdep.gt.1000.) .and. (xFeT.lt.0.5)) then
-                     ffescav = 0.
+                  if ((fsdepw(ji,jj,jk).gt.1000.) .and. (xFeT.lt.0.5)) then
+                     ffescav(ji,jj) = 0.
                   endif
 # endif
@@ -2684 +2659 @@
                   !!
                   !! total iron concentration (mmol Fe / m3 -> umol Fe / m3)
-                  xFeT        = zfer * 1.e3
+                  xFeT        = zfer(ji,jj) * 1.e3
                   !!
                   !! this has a base scavenging rate (12% / y) which is modified by local
@@ -2713 +2688 @@
                   !! by a magic factor, 0.002, to get it into nmol C / cm2 / s
                   !!
-                  fscal_part = ((xthetapn * zphn) + (xthetapd * zphd) + (xthetazmi * zzmi) + &
-                  (xthetazme * zzme) + (xthetad * zdet)) * 0.002
+                  fscal_part = ((xthetapn * zphn(ji,jj)) + (xthetapd * zphd(ji,jj)) + (xthetazmi * zzmi(ji,jj)) + &
+                  (xthetazme * zzme(ji,jj)) + (xthetad * zdet(ji,jj))) * 0.002
                   !!
                   !! calculate scaling factor for base scavenging rate
@@ -2748 +2723 @@
                   !! apply the calculated scavenging rate ...
                   !!
-                  ffescav = fscal_scav * zfer
+                  ffescav(ji,jj) = fscal_scav * zfer(ji,jj)
                   !!
                elseif (jiron.eq.3) then
@@ -2762 +2737 @@
                   !!
                   !! total iron concentration (mmol Fe / m3 -> umol Fe / m3)
-                  xFeT        = zfer * 1.e3
+                  xFeT        = zfer(ji,jj) * 1.e3
                   !!
                   !! this has a base scavenging rate which is modified by local
@@ -2822 +2797 @@
                   !! apply the calculated scavenging rate ...
                   !!
-                  ffescav = fscal_scav * zfer
+                  ffescav(ji,jj) = fscal_scav * zfer(ji,jj)
                   !!
                elseif (jiron.eq.4) then
@@ -2861 +2836 @@
                   !!
                   !! sum these two terms together
-                  ffescav     = xORGscav + xINORGscav
+                  ffescav(ji,jj)     = xORGscav + xINORGscav
                else
                   !!----------------------------------------------------------------------
@@ -2869 +2844 @@
                   !! you fool!
                   !!----------------------------------------------------------------------
-                  ffescav     = 0.
+                  ffescav(ji,jj)     = 0.
                endif
 
@@ -2879 +2854 @@
                if (jk.eq.1) then
                   !! zirondep   is in mmol-Fe / m2 / day
-                  !! ffetop     is in mmol-dissolved-Fe / m3 / day
-                  ffetop  = zirondep(ji,jj) * xfe_sol / fthk 
+                  !! ffetop(ji,jj)     is in mmol-dissolved-Fe / m3 / day
+                  ffetop(ji,jj)  = zirondep(ji,jj) * xfe_sol / fse3t(ji,jj,jk) 
                else
-                  ffetop  = 0.0
+                  ffetop(ji,jj)  = 0.0
                endif
                !!
@@ -2888 +2863 @@
                !! AXY (10/07/12): amended to only apply sedimentary flux up to ~500 m down
                !! if (jk.eq.(mbathy(ji,jj)-1).AND.jk.lt.i1100) then
-               if ((jk.eq.jmbathy).AND.jk.le.i0500) then
+               if ((jk.eq.mbathy(ji,jj)).AND.jk.le.i0500) then
                   !! Moore et al. (2004) cite a coastal California value of 5 umol/m2/d, but adopt a
                   !! global value of 2 umol/m2/d for all areas < 1100 m; here we use this latter value
                   !! but apply it everywhere
                   !! AXY (21/07/09): actually, let's just apply it below 1100 m (levels 1-37)
-                  ffebot  = (xfe_sed / fthk)
+                  ffebot(ji,jj)  = (xfe_sed / fse3t(ji,jj,jk))
                else
-                  ffebot  = 0.0
+                  ffebot(ji,jj)  = 0.0
                endif
 
@@ -2903 +2878 @@
                !! iron addition at the seafloor is still permitted with the idea that this extra
                !! iron will be removed by the restoration away from the source
-               !! ffescav = 0.0
-               !! ffetop  = 0.0
-               !! ffebot  = 0.0
+               !! ffescav(ji,jj) = 0.0
+               !! ffetop(ji,jj)  = 0.0
+               !! ffebot(ji,jj)  = 0.0
 
 # if defined key_debug_medusa
@@ -2913 +2888 @@
                   IF (lwp) write (numout,*) 'xfe_sol  = ', xfe_sol
                   IF (lwp) write (numout,*) 'xfe_mass = ', xfe_mass
-                  IF (lwp) write (numout,*) 'ffetop(',jk,')  = ', ffetop
-                  IF (lwp) write (numout,*) 'ffebot(',jk,')  = ', ffebot
+                  IF (lwp) write (numout,*) 'ffetop(',jk,')  = ', ffetop(ji,jj)
+                  IF (lwp) write (numout,*) 'ffebot(',jk,')  = ', ffebot(ji,jj)
                   IF (lwp) write (numout,*) 'xFree(',jk,')   = ', xFree(ji,jj)
-                  IF (lwp) write (numout,*) 'ffescav(',jk,') = ', ffescav
+                  IF (lwp) write (numout,*) 'ffescav(',jk,') = ', ffescav(ji,jj)
                endif
 # endif
@@ -2925 +2900 @@
                !!
                !! diatom frustule dissolution
-               fsdiss  = xsdiss * zpds
+               fsdiss(ji,jj)  = xsdiss * zpds(ji,jj)
 
 # if defined key_debug_medusa
@@ -2931 +2906 @@
                if (idf.eq.1.AND.idfval.eq.1) then
                   IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'fsdiss(',jk,')  = ', fsdiss
+                  IF (lwp) write (numout,*) 'fsdiss(',jk,')  = ', fsdiss(ji,jj)
                endif
 # endif
@@ -2941 +2916 @@
                !! to allow the split between fast and slow detritus to be 
                !! diagnosed
-               fslown  = fdpn + fdzmi + ((1.0 - xfdfrac1) * fdpd) + &
-               ((1.0 - xfdfrac2) * fdzme) + ((1.0 - xbetan) * (finmi + finme))
+               fslown(ji,jj)  = fdpn(ji,jj) + fdzmi(ji,jj) + ((1.0 - xfdfrac1) * fdpd(ji,jj)) + &
+               ((1.0 - xfdfrac2) * fdzme(ji,jj)) + ((1.0 - xbetan) * (finmi(ji,jj) + finme(ji,jj)))
                !!
                !! this variable records the slow detrital sinking flux at this
@@ -2948 +2923 @@
                !! standard depths in the diagnostic outputs; needs to be
                !! adjusted from per second to per day because of parameter vsed
-               fslownflux(ji,jj) = zdet * vsed * 86400.
+               fslownflux(ji,jj) = zdet(ji,jj) * vsed * 86400.
 # if defined key_roam
                !!
                !! and the same for detrital carbon
-               fslowc  = (xthetapn * fdpn) + (xthetazmi * fdzmi) + &
-               (xthetapd * (1.0 - xfdfrac1) * fdpd) + &
-               (xthetazme * (1.0 - xfdfrac2) * fdzme) + &
-               ((1.0 - xbetac) * (ficmi + ficme))
+               fslowc(ji,jj)  = (xthetapn * fdpn(ji,jj)) + (xthetazmi * fdzmi(ji,jj)) + &
+               (xthetapd * (1.0 - xfdfrac1) * fdpd(ji,jj)) + &
+               (xthetazme * (1.0 - xfdfrac2) * fdzme(ji,jj)) + &
+               ((1.0 - xbetac) * (ficmi(ji,jj) + ficme(ji,jj)))
                !!
                !! this variable records the slow detrital sinking flux at this
@@ -2961 +2936 @@
                !! standard depths in the diagnostic outputs; needs to be
                !! adjusted from per second to per day because of parameter vsed
-               fslowcflux(ji,jj) = zdtc * vsed * 86400.
+               fslowcflux(ji,jj) = zdtc(ji,jj) * vsed * 86400.
 # endif
 
@@ -2975 +2950 @@
                !!
                !! nitrogen
-               fregen   = (( (xphi * (fgmipn + fgmid)) +                &  ! messy feeding
-               (xphi * (fgmepn + fgmepd + fgmezmi + fgmed)) +           &  ! messy feeding
-               fmiexcr + fmeexcr + fdd +                                &  ! excretion + D remin.
-               fdpn2 + fdpd2 + fdzmi2 + fdzme2) * fthk)                    ! linear mortality
+               fregen(ji,jj)   = (( (xphi * (fgmipn(ji,jj) + fgmid(ji,jj))) +                &  ! messy feeding
+               (xphi * (fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj))) +           &  ! messy feeding
+               fmiexcr(ji,jj) + fmeexcr(ji,jj) + fdd(ji,jj) +                                &  ! excretion + D remin.
+               fdpn2(ji,jj) + fdpd2(ji,jj) + fdzmi2(ji,jj) + fdzme2(ji,jj)) * fse3t(ji,jj,jk))                    ! linear mortality
                !!
                !! silicon
-               fregensi = (( fsdiss + ((1.0 - xfdfrac1) * fdpds) +      &  ! dissolution + non-lin. mortality
-               ((1.0 - xfdfrac3) * fgmepds) +                           &  ! egestion by zooplankton
-               fdpds2) * fthk)                                             ! linear mortality
+               fregensi(ji,jj) = (( fsdiss(ji,jj) + ((1.0 - xfdfrac1) * fdpds(ji,jj)) +      &  ! dissolution + non-lin. mortality
+               ((1.0 - xfdfrac3) * fgmepds(ji,jj)) +                           &  ! egestion by zooplankton
+               fdpds2(ji,jj)) * fse3t(ji,jj,jk))                                             ! linear mortality
 # if defined key_roam
                !!
                !! carbon
-               fregenc  = (( (xphi * ((xthetapn * fgmipn) + fgmidc)) +  &  ! messy feeding
-               (xphi * ((xthetapn * fgmepn) + (xthetapd * fgmepd) +     &  ! messy feeding
-               (xthetazmi * fgmezmi) + fgmedc)) +                       &  ! messy feeding
-               fmiresp + fmeresp + fddc +                               &  ! respiration + D remin.
-               (xthetapn * fdpn2) + (xthetapd * fdpd2) +                &  ! linear mortality
-               (xthetazmi * fdzmi2) + (xthetazme * fdzme2)) * fthk)        ! linear mortality
+! Doesn't look this is used - marc 10/4/17
+!               fregenc(ji,jj)  = (( (xphi * ((xthetapn * fgmipn(ji,jj)) + fgmidc(ji,jj))) +  &  ! messy feeding
+!               (xphi * ((xthetapn * fgmepn(ji,jj)) + (xthetapd * fgmepd(ji,jj)) +     &  ! messy feeding
+!               (xthetazmi * fgmezmi(ji,jj)) + fgmedc(ji,jj))) +                       &  ! messy feeding
+!               fmiresp(ji,jj) + fmeresp(ji,jj) + fddc(ji,jj) +                               &  ! respiration + D remin.
+!               (xthetapn * fdpn2(ji,jj)) + (xthetapd * fdpd2(ji,jj)) +                &  ! linear mortality
+!               (xthetazmi * fdzmi2(ji,jj)) + (xthetazme * fdzme2(ji,jj))) * fse3t(ji,jj,jk))        ! linear mortality
 # endif
 
@@ -3005 +2981 @@
                !!
                !! nitrogen:   diatom and mesozooplankton mortality
-               ftempn         = b0 * ((xfdfrac1 * fdpd)  + (xfdfrac2 * fdzme))
+               ftempn(ji,jj)         = b0 * ((xfdfrac1 * fdpd(ji,jj))  + (xfdfrac2 * fdzme(ji,jj)))
                !!
                !! silicon:    diatom mortality and grazed diatoms
-               ftempsi        = b0 * ((xfdfrac1 * fdpds) + (xfdfrac3 * fgmepds))
+               ftempsi(ji,jj)        = b0 * ((xfdfrac1 * fdpds(ji,jj)) + (xfdfrac3 * fgmepds(ji,jj)))
                !!
                !! iron:       diatom and mesozooplankton mortality
-               ftempfe        = b0 * (((xfdfrac1 * fdpd) + (xfdfrac2 * fdzme)) * xrfn)
+               ftempfe(ji,jj)        = b0 * (((xfdfrac1 * fdpd(ji,jj)) + (xfdfrac2 * fdzme(ji,jj))) * xrfn)
                !!
                !! carbon:     diatom and mesozooplankton mortality
-               ftempc         = b0 * ((xfdfrac1 * xthetapd * fdpd) + & 
-                                (xfdfrac2 * xthetazme * fdzme))
+               ftempc(ji,jj)         = b0 * ((xfdfrac1 * xthetapd * fdpd(ji,jj)) + & 
+                                (xfdfrac2 * xthetazme * fdzme(ji,jj)))
                !!
 # if defined key_roam
                if (jrratio.eq.0) then
                   !! CaCO3:      latitudinally-based fraction of total primary production
-                  !!               absolute latitude of current grid cell
-                  flat           = abs(gphit(ji,jj))
                   !!               0.10 at equator; 0.02 at pole
-                  fcaco3         = xcaco3a + ((xcaco3b - xcaco3a) * ((90.0 - flat) / 90.0))
+                  fcaco3(ji,jj)         = xcaco3a + ((xcaco3b - xcaco3a) * ((90.0 - abs(gphit(ji,jj))) / 90.0))
                elseif (jrratio.eq.1) then
                   !! CaCO3:      Ridgwell et al. (2007) submodel, version 1
@@ -3032 +3006 @@
                      fq1 = 0.
                   endif
-                  fcaco3 = xridg_r0 * fq1
+                  fcaco3(ji,jj) = xridg_r0 * fq1
                elseif (jrratio.eq.2) then
                   !! CaCO3:      Ridgwell et al. (2007) submodel, version 2
@@ -3041 +3015 @@
                      fq1 = 0.
                   endif
-                  fcaco3 = xridg_r0 * fq1
+                  fcaco3(ji,jj) = xridg_r0 * fq1
                endif
 # else
                !! CaCO3:      latitudinally-based fraction of total primary production
-               !!               absolute latitude of current grid cell
-               flat           = abs(gphit(ji,jj))
                !!               0.10 at equator; 0.02 at pole
-               fcaco3         = xcaco3a + ((xcaco3b - xcaco3a) * ((90.0 - flat) / 90.0))
+               fcaco3(ji,jj)         = xcaco3a + ((xcaco3b - xcaco3a) * ((90.0 - abs(gphit(ji,jj))) / 90.0))
 # endif
                !! AXY (09/03/09): convert CaCO3 production from function of 
@@ -3056 +3028 @@
                !! chlorophyll to an export flux for which they apply conversion 
                !! factors to estimate the various elemental fractions (Si, Ca)
-               ftempca        = ftempc * fcaco3
+               ftempca(ji,jj)        = ftempc(ji,jj) * fcaco3(ji,jj)
 
 # if defined key_debug_medusa
                !! integrate total fast detritus production
                if (idf.eq.1) then
-                  fifd_n(ji,jj)  = fifd_n(ji,jj)  + (ftempn  * fthk)
-                  fifd_si(ji,jj) = fifd_si(ji,jj) + (ftempsi * fthk)
-                  fifd_fe(ji,jj) = fifd_fe(ji,jj) + (ftempfe * fthk)
+                  fifd_n(ji,jj)  = fifd_n(ji,jj)  + (ftempn(ji,jj)  * fse3t(ji,jj,jk))
+                  fifd_si(ji,jj) = fifd_si(ji,jj) + (ftempsi(ji,jj) * fse3t(ji,jj,jk))
+                  fifd_fe(ji,jj) = fifd_fe(ji,jj) + (ftempfe(ji,jj) * fse3t(ji,jj,jk))
 #  if defined key_roam
-                  fifd_c(ji,jj)  = fifd_c(ji,jj)  + (ftempc  * fthk)
+                  fifd_c(ji,jj)  = fifd_c(ji,jj)  + (ftempc(ji,jj)  * fse3t(ji,jj,jk))
 #  endif
                endif
@@ -3072 +3044 @@
                if (idf.eq.1.AND.idfval.eq.1) then
                   IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'fdpd(',jk,')    = ', fdpd
-                  IF (lwp) write (numout,*) 'fdzme(',jk,')   = ', fdzme
-                  IF (lwp) write (numout,*) 'ftempn(',jk,')  = ', ftempn
-                  IF (lwp) write (numout,*) 'ftempsi(',jk,') = ', ftempsi
-                  IF (lwp) write (numout,*) 'ftempfe(',jk,') = ', ftempfe
-                  IF (lwp) write (numout,*) 'ftempc(',jk,')  = ', ftempc
-                  IF (lwp) write (numout,*) 'ftempca(',jk,') = ', ftempca
-                  IF (lwp) write (numout,*) 'flat(',jk,')    = ', flat
-                  IF (lwp) write (numout,*) 'fcaco3(',jk,')  = ', fcaco3
+                  IF (lwp) write (numout,*) 'fdpd(',jk,')    = ', fdpd(ji,jj)
+                  IF (lwp) write (numout,*) 'fdzme(',jk,')   = ', fdzme(ji,jj)
+                  IF (lwp) write (numout,*) 'ftempn(',jk,')  = ', ftempn(ji,jj)
+                  IF (lwp) write (numout,*) 'ftempsi(',jk,') = ', ftempsi(ji,jj)
+                  IF (lwp) write (numout,*) 'ftempfe(',jk,') = ', ftempfe(ji,jj)
+                  IF (lwp) write (numout,*) 'ftempc(',jk,')  = ', ftempc(ji,jj)
+                  IF (lwp) write (numout,*) 'ftempca(',jk,') = ', ftempca(ji,jj)
+                  IF (lwp) write (numout,*) 'flat(',jk,')    = ', abs(gphit(ji,jj))
+                  IF (lwp) write (numout,*) 'fcaco3(',jk,')  = ', fcaco3(ji,jj)
                endif
 # endif
@@ -3151 +3123 @@
                      !! this is the SURFACE OCEAN BOX (no remineralisation)
                      !!
-                     freminc  = 0.0
-                     freminn  = 0.0
-                     freminfe = 0.0
-                     freminsi = 0.0
-                     freminca = 0.0
-                  elseif (jk.le.jmbathy) then
+                     freminc(ji,jj)  = 0.0
+                     freminn(ji,jj)  = 0.0
+                     freminfe(ji,jj) = 0.0
+                     freminsi(ji,jj) = 0.0
+                     freminca(ji,jj) = 0.0
+                  elseif (jk.le.mbathy(ji,jj)) then
                      !! this is an OCEAN BOX (remineralise some material)
                      !!
@@ -3162 +3134 @@
                      if (jocalccd.eq.0) then
                         !! use default CCD field
-                        fccd_dep = ocal_ccd(ji,jj)
+                        fccd_dep(ji,jj) = ocal_ccd(ji,jj)
                      elseif (jocalccd.eq.1) then
                         !! use calculated CCD field
-                        fccd_dep = f2_ccd_cal(ji,jj)
+                        fccd_dep(ji,jj) = f2_ccd_cal(ji,jj)
                      endif
                      !!
@@ -3178 +3150 @@
                         !! it needs to be protected in case ALL C is protected
                         if (fq4.lt.fq1) then
-                           fprotf   = (fq4 / (fq1 + tiny(fq1)))      !! protected fraction of total organic C     (non-dim)
+                           fprotf(ji,jj)   = (fq4 / (fq1 + tiny(fq1)))      !! protected fraction of total organic C     (non-dim)
                         else
-                           fprotf   = 1.0                            !! all organic C is protected                (non-dim)
+                           fprotf(ji,jj)   = 1.0                            !! all organic C is protected                (non-dim)
                         endif
-                        fq5      = (1.0 - fprotf)                    !! unprotected fraction of total organic C   (non-dim)
+                        fq5      = (1.0 - fprotf(ji,jj))                    !! unprotected fraction of total organic C   (non-dim)
                         fq6      = (fq0 * fq5)                       !! how much organic C is unprotected         (mol)
-                        fq7      = (fq6 * exp(-(fthk / xfastc)))     !! how much unprotected C leaves this box    (mol)
-                        fq8      = (fq7 + (fq0 * fprotf))            !! how much total C leaves this box          (mol)
-                        freminc  = (fq0 - fq8) / fthk                !! C remineralisation in this box            (mol)
+                        fq7      = (fq6 * exp(-(fse3t(ji,jj,jk) / xfastc)))     !! how much unprotected C leaves this box    (mol)
+                        fq8      = (fq7 + (fq0 * fprotf(ji,jj)))            !! how much total C leaves this box          (mol)
+                        freminc(ji,jj)  = (fq0 - fq8) / fse3t(ji,jj,jk)                !! C remineralisation in this box            (mol)
                         ffastc(ji,jj) = fq8                          
 # if defined key_debug_medusa
@@ -3194 +3166 @@
                               IF (lwp) write (numout,*) 'totalC(',jk,')  = ', fq1
                               IF (lwp) write (numout,*) 'prtctC(',jk,')  = ', fq4
-                              IF (lwp) write (numout,*) 'fprotf(',jk,')  = ', fprotf
+                              IF (lwp) write (numout,*) 'fprotf(',jk,')  = ', fprotf(ji,jj)
                               IF (lwp) write (numout,*) '------------------------------'
                               IF (lwp) write (numout,*) 'IN   C(',jk,')  = ', fq0
-                              IF (lwp) write (numout,*) 'LOST C(',jk,')  = ', freminc * fthk
+                              IF (lwp) write (numout,*) 'LOST C(',jk,')  = ', freminc(ji,jj) * fse3t(ji,jj,jk)
                               IF (lwp) write (numout,*) 'OUT  C(',jk,')  = ', fq8
-                              IF (lwp) write (numout,*) 'NEW  C(',jk,')  = ', ftempc * fthk
+                              IF (lwp) write (numout,*) 'NEW  C(',jk,')  = ', ftempc(ji,jj) * fse3t(ji,jj,jk)
                            endif
 # endif
                         else
-                        fq1      = fq0 * exp(-(fthk / xfastc))       !! how much organic C leaves this box        (mol)
-                        freminc  = (fq0 - fq1) / fthk                !! C remineralisation in this box            (mol)
+                        fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastc))       !! how much organic C leaves this box        (mol)
+                        freminc(ji,jj)  = (fq0 - fq1) / fse3t(ji,jj,jk)                !! C remineralisation in this box            (mol)
                         ffastc(ji,jj)  = fq1
                      endif
@@ -3211 +3183 @@
                      fq0      = ffastn(ji,jj)                        !! how much organic N enters this box        (mol)
                      if (iball.eq.1) then
-                        fq5      = (1.0 - fprotf)                    !! unprotected fraction of total organic N   (non-dim)
+                        fq5      = (1.0 - fprotf(ji,jj))                    !! unprotected fraction of total organic N   (non-dim)
                         fq6      = (fq0 * fq5)                       !! how much organic N is unprotected         (mol)
-                        fq7      = (fq6 * exp(-(fthk / xfastc)))     !! how much unprotected N leaves this box    (mol)
-                        fq8      = (fq7 + (fq0 * fprotf))            !! how much total N leaves this box          (mol)
-                        freminn  = (fq0 - fq8) / fthk                !! N remineralisation in this box            (mol)
+                        fq7      = (fq6 * exp(-(fse3t(ji,jj,jk) / xfastc)))     !! how much unprotected N leaves this box    (mol)
+                        fq8      = (fq7 + (fq0 * fprotf(ji,jj)))            !! how much total N leaves this box          (mol)
+                        freminn(ji,jj)  = (fq0 - fq8) / fse3t(ji,jj,jk)                !! N remineralisation in this box            (mol)
                         ffastn(ji,jj)  = fq8
 # if defined key_debug_medusa
@@ -3223 +3195 @@
                            IF (lwp) write (numout,*) 'totalN(',jk,')  = ', fq1
                            IF (lwp) write (numout,*) 'prtctN(',jk,')  = ', fq4
-                           IF (lwp) write (numout,*) 'fprotf(',jk,')  = ', fprotf
+                           IF (lwp) write (numout,*) 'fprotf(',jk,')  = ', fprotf(ji,jj)
                            IF (lwp) write (numout,*) '------------------------------'
-                           if (freminn < 0.0) then
+                           if (freminn(ji,jj) < 0.0) then
                               IF (lwp) write (numout,*) '** FREMIN ERROR **'
                            endif
                            IF (lwp) write (numout,*) 'IN   N(',jk,')  = ', fq0
-                           IF (lwp) write (numout,*) 'LOST N(',jk,')  = ', freminn * fthk
+                           IF (lwp) write (numout,*) 'LOST N(',jk,')  = ', freminn(ji,jj) * fse3t(ji,jj,jk)
                            IF (lwp) write (numout,*) 'OUT  N(',jk,')  = ', fq8
-                           IF (lwp) write (numout,*) 'NEW  N(',jk,')  = ', ftempn * fthk
+                           IF (lwp) write (numout,*) 'NEW  N(',jk,')  = ', ftempn(ji,jj) * fse3t(ji,jj,jk)
                         endif
 # endif
                      else
-                        fq1      = fq0 * exp(-(fthk / xfastc))       !! how much organic N leaves this box        (mol)
-                        freminn  = (fq0 - fq1) / fthk                !! N remineralisation in this box            (mol)
+                        fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastc))       !! how much organic N leaves this box        (mol)
+                        freminn(ji,jj)  = (fq0 - fq1) / fse3t(ji,jj,jk)                !! N remineralisation in this box            (mol)
                         ffastn(ji,jj)  = fq1
                      endif
@@ -3243 +3215 @@
                      fq0      = ffastfe(ji,jj)                       !! how much organic Fe enters this box       (mol)
                      if (iball.eq.1) then
-                        fq5      = (1.0 - fprotf)                    !! unprotected fraction of total organic Fe  (non-dim)
+                        fq5      = (1.0 - fprotf(ji,jj))                    !! unprotected fraction of total organic Fe  (non-dim)
                         fq6      = (fq0 * fq5)                       !! how much organic Fe is unprotected        (mol)
-                        fq7      = (fq6 * exp(-(fthk / xfastc)))     !! how much unprotected Fe leaves this box   (mol)
-                        fq8      = (fq7 + (fq0 * fprotf))            !! how much total Fe leaves this box         (mol)            
-                        freminfe = (fq0 - fq8) / fthk                !! Fe remineralisation in this box           (mol)
+                        fq7      = (fq6 * exp(-(fse3t(ji,jj,jk) / xfastc)))     !! how much unprotected Fe leaves this box   (mol)
+                        fq8      = (fq7 + (fq0 * fprotf(ji,jj)))            !! how much total Fe leaves this box         (mol)            
+                        freminfe(ji,jj) = (fq0 - fq8) / fse3t(ji,jj,jk)                !! Fe remineralisation in this box           (mol)
                         ffastfe(ji,jj) = fq8
                      else
-                        fq1      = fq0 * exp(-(fthk / xfastc))       !! how much total Fe leaves this box         (mol)      
-                        freminfe = (fq0 - fq1) / fthk                !! Fe remineralisation in this box           (mol)
+                        fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastc))       !! how much total Fe leaves this box         (mol)      
+                        freminfe(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)                !! Fe remineralisation in this box           (mol)
                         ffastfe(ji,jj) = fq1
                      endif
@@ -3257 +3229 @@
                      !! === biogenic silicon ===
                      fq0      = ffastsi(ji,jj)                       !! how much  opal centers this box           (mol) 
-                     fq1      = fq0 * exp(-(fthk / xfastsi))         !! how much  opal leaves this box            (mol)
-                     freminsi = (fq0 - fq1) / fthk                   !! Si remineralisation in this box           (mol)
+                     fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastsi))         !! how much  opal leaves this box            (mol)
+                     freminsi(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)                   !! Si remineralisation in this box           (mol)
                      ffastsi(ji,jj) = fq1
                      !!
                      !! === biogenic calcium carbonate ===
                      fq0      = ffastca(ji,jj)                       !! how much CaCO3 enters this box            (mol)
-                     if (fdep.le.fccd_dep) then
+                     if (fsdepw(ji,jj,jk).le.fccd_dep(ji,jj)) then
                         !! whole grid cell above CCD
                         fq1      = fq0                               !! above lysocline, no Ca dissolves          (mol)
-                        freminca = 0.0                               !! above lysocline, no Ca dissolves          (mol)
+                        freminca(ji,jj) = 0.0                               !! above lysocline, no Ca dissolves          (mol)
                         fccd(ji,jj) = real(jk)                       !! which is the last level above the CCD?    (#)
-                     elseif (fdep.ge.fccd_dep) then
+                     elseif (fsdepw(ji,jj,jk).ge.fccd_dep(ji,jj)) then
                         !! whole grid cell below CCD
-                        fq1      = fq0 * exp(-(fthk / xfastca))      !! how much CaCO3 leaves this box            (mol)
-                        freminca = (fq0 - fq1) / fthk                !! Ca remineralisation in this box           (mol)
+                        fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastca))      !! how much CaCO3 leaves this box            (mol)
+                        freminca(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)                !! Ca remineralisation in this box           (mol)
                      else
                         !! partial grid cell below CCD
-                        fq2      = fdep1 - fccd_dep                  !! amount of grid cell below CCD             (m)
+                        fq2      = fdep1(ji,jj) - fccd_dep(ji,jj)                  !! amount of grid cell below CCD             (m)
                         fq1      = fq0 * exp(-(fq2 / xfastca))       !! how much CaCO3 leaves this box            (mol)
-                        freminca = (fq0 - fq1) / fthk                !! Ca remineralisation in this box           (mol)
+                        freminca(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)                !! Ca remineralisation in this box           (mol)
                      endif
                      ffastca(ji,jj) = fq1 
                   else
                      !! this is BELOW THE LAST OCEAN BOX (do nothing)
-                     freminc  = 0.0
-                     freminn  = 0.0
-                     freminfe = 0.0
-                     freminsi = 0.0
-                     freminca = 0.0              
+                     freminc(ji,jj)  = 0.0
+                     freminn(ji,jj)  = 0.0
+                     freminfe(ji,jj) = 0.0
+                     freminsi(ji,jj) = 0.0
+                     freminca(ji,jj) = 0.0              
                   endif
 
@@ -3334 +3306 @@
                      !! this is the SURFACE OCEAN BOX (no remineralisation)
                      !!
-                     freminc  = 0.0
-                     freminn  = 0.0
-                     freminfe = 0.0
-                     freminsi = 0.0
-                     freminca = 0.0
-                  elseif (jk.le.jmbathy) then
+                     freminc(ji,jj)  = 0.0
+                     freminn(ji,jj)  = 0.0
+                     freminfe(ji,jj) = 0.0
+                     freminsi(ji,jj) = 0.0
+                     freminca(ji,jj) = 0.0
+                  elseif (jk.le.mbathy(ji,jj)) then
                      !! this is an OCEAN BOX (remineralise some material)
                      !!
                      !! === organic carbon ===
                      fq0      = ffastc(ji,jj)                        !! how much organic C enters this box        (mol)
-                     fq1      = fq0 * ((fdep1/fdep)**fb_val)         !! how much organic C leaves this box        (mol)
-                     freminc  = (fq0 - fq1) / fthk                   !! C remineralisation in this box            (mol)
+                     fq1      = fq0 * ((fdep1(ji,jj)/fsdepw(ji,jj,jk))**fb_val)         !! how much organic C leaves this box        (mol)
+                     freminc(ji,jj)  = (fq0 - fq1) / fse3t(ji,jj,jk)                   !! C remineralisation in this box            (mol)
                      ffastc(ji,jj)  = fq1
                      !!
                      !! === organic nitrogen ===
                      fq0      = ffastn(ji,jj)                        !! how much organic N enters this box        (mol)
-                     fq1      = fq0 * ((fdep1/fdep)**fb_val)         !! how much organic N leaves this box        (mol)
-                     freminn  = (fq0 - fq1) / fthk                   !! N remineralisation in this box            (mol)
+                     fq1      = fq0 * ((fdep1(ji,jj)/fsdepw(ji,jj,jk))**fb_val)         !! how much organic N leaves this box        (mol)
+                     freminn(ji,jj)  = (fq0 - fq1) / fse3t(ji,jj,jk)                   !! N remineralisation in this box            (mol)
                      ffastn(ji,jj)  = fq1
                      !!
                      !! === organic iron ===
                      fq0      = ffastfe(ji,jj)                       !! how much organic Fe enters this box       (mol)
-                     fq1      = fq0 * ((fdep1/fdep)**fb_val)         !! how much organic Fe leaves this box       (mol)
-                     freminfe = (fq0 - fq1) / fthk                   !! Fe remineralisation in this box           (mol)
+                     fq1      = fq0 * ((fdep1(ji,jj)/fsdepw(ji,jj,jk))**fb_val)         !! how much organic Fe leaves this box       (mol)
+                     freminfe(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)                   !! Fe remineralisation in this box           (mol)
                      ffastfe(ji,jj) = fq1
                      !!
                      !! === biogenic silicon ===
                      fq0      = ffastsi(ji,jj)                       !! how much  opal centers this box           (mol) 
-                     fq1      = fq0 * exp(-(fthk / xfastsi))         !! how much  opal leaves this box            (mol)
-                     freminsi = (fq0 - fq1) / fthk                   !! Si remineralisation in this box           (mol)
+                     fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastsi))         !! how much  opal leaves this box            (mol)
+                     freminsi(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)                   !! Si remineralisation in this box           (mol)
                      ffastsi(ji,jj) = fq1
                      !!
                      !! === biogenic calcium carbonate ===
                      fq0      = ffastca(ji,jj)                       !! how much CaCO3 enters this box            (mol)
-                     if (fdep.le.ocal_ccd(ji,jj)) then
+                     if (fsdepw(ji,jj,jk).le.ocal_ccd(ji,jj)) then
                         !! whole grid cell above CCD
                         fq1      = fq0                               !! above lysocline, no Ca dissolves          (mol)
-                        freminca = 0.0                               !! above lysocline, no Ca dissolves          (mol)
+                        freminca(ji,jj) = 0.0                               !! above lysocline, no Ca dissolves          (mol)
                         fccd(ji,jj) = real(jk)                       !! which is the last level above the CCD?    (#)
-                     elseif (fdep.ge.ocal_ccd(ji,jj)) then
+                     elseif (fsdepw(ji,jj,jk).ge.ocal_ccd(ji,jj)) then
                         !! whole grid cell below CCD
-                        fq1      = fq0 * exp(-(fthk / xfastca))      !! how much CaCO3 leaves this box            (mol)
-                        freminca = (fq0 - fq1) / fthk                !! Ca remineralisation in this box           (mol)
+                        fq1      = fq0 * exp(-(fse3t(ji,jj,jk) / xfastca))      !! how much CaCO3 leaves this box            (mol)
+                        freminca(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)                !! Ca remineralisation in this box           (mol)
                      else
                         !! partial grid cell below CCD
-                        fq2      = fdep1 - ocal_ccd(ji,jj)           !! amount of grid cell below CCD             (m)
+                        fq2      = fdep1(ji,jj) - ocal_ccd(ji,jj)           !! amount of grid cell below CCD             (m)
                         fq1      = fq0 * exp(-(fq2 / xfastca))       !! how much CaCO3 leaves this box            (mol)
-                        freminca = (fq0 - fq1) / fthk                !! Ca remineralisation in this box           (mol)
+                        freminca(ji,jj) = (fq0 - fq1) / fse3t(ji,jj,jk)                !! Ca remineralisation in this box           (mol)
                      endif
                      ffastca(ji,jj) = fq1 
                   else
                      !! this is BELOW THE LAST OCEAN BOX (do nothing)
-                     freminc  = 0.0
-                     freminn  = 0.0
-                     freminfe = 0.0
-                     freminsi = 0.0
-                     freminca = 0.0              
+                     freminc(ji,jj)  = 0.0
+                     freminn(ji,jj)  = 0.0
+                     freminfe(ji,jj) = 0.0
+                     freminsi(ji,jj) = 0.0
+                     freminca(ji,jj) = 0.0              
                   endif
 
@@ -3406 +3378 @@
                !!
                !! === organic carbon ===                          (diatom and mesozooplankton mortality)
-               ffastc(ji,jj)  = ffastc(ji,jj)  + (ftempc  * fthk)
+               ffastc(ji,jj)  = ffastc(ji,jj)  + (ftempc(ji,jj)  * fse3t(ji,jj,jk))
                !!
                !! === organic nitrogen ===                        (diatom and mesozooplankton mortality)
-               ffastn(ji,jj)  = ffastn(ji,jj)  + (ftempn  * fthk)
+               ffastn(ji,jj)  = ffastn(ji,jj)  + (ftempn(ji,jj)  * fse3t(ji,jj,jk))
                !!
                !! === organic iron ===                            (diatom and mesozooplankton mortality)
-               ffastfe(ji,jj) = ffastfe(ji,jj) + (ftempfe * fthk)
+               ffastfe(ji,jj) = ffastfe(ji,jj) + (ftempfe(ji,jj) * fse3t(ji,jj,jk))
                !!
                !! === biogenic silicon ===                        (diatom mortality and grazed diatoms)
-               ffastsi(ji,jj) = ffastsi(ji,jj) + (ftempsi * fthk)
+               ffastsi(ji,jj) = ffastsi(ji,jj) + (ftempsi(ji,jj) * fse3t(ji,jj,jk))
                !!
                !! === biogenic calcium carbonate ===              (latitudinally-based fraction of total primary production)
-               ffastca(ji,jj) = ffastca(ji,jj) + (ftempca * fthk)
+               ffastca(ji,jj) = ffastca(ji,jj) + (ftempca(ji,jj) * fse3t(ji,jj,jk))
 
                !!----------------------------------------------------------------------
@@ -3456 +3428 @@
                !!----------------------------------------------------------------------
                !! 
-               ffast2slowc  = 0.0
-               ffast2slown  = 0.0
-               ffast2slowfe = 0.0
-               !!
-               if (jk.eq.jmbathy) then
+               ffast2slowc(ji,jj)  = 0.0
+               ffast2slown(ji,jj)  = 0.0
+! I don't think this is used - marc 10/4/17
+!               ffast2slowfe(ji,jj) = 0.0
+               !!
+               if (jk.eq.mbathy(ji,jj)) then
                   !! this is the BOTTOM OCEAN BOX (remineralise everything)
                   !!
@@ -3467 +3440 @@
                   !! === organic carbon ===
                   if (jfdfate.eq.0 .and. jorgben.eq.0) then
-                     freminc  = freminc + (ffastc(ji,jj) / fthk)    !! C remineralisation in this box            (mol/m3)
+                     freminc(ji,jj)  = freminc(ji,jj) + (ffastc(ji,jj) / fse3t(ji,jj,jk))    !! C remineralisation in this box            (mol/m3)
                   elseif (jfdfate.eq.1 .and. jorgben.eq.0) then
-                     ffast2slowc = ffastc(ji,jj) / fthk             !! fast C -> slow C                          (mol/m3)
-                     fslowc      = fslowc + ffast2slowc
+                     ffast2slowc(ji,jj) = ffastc(ji,jj) / fse3t(ji,jj,jk)             !! fast C -> slow C                          (mol/m3)
+                     fslowc(ji,jj)      = fslowc(ji,jj) + ffast2slowc(ji,jj)
                   elseif (jfdfate.eq.0 .and. jorgben.eq.1) then
                      f_fbenin_c(ji,jj)  = ffastc(ji,jj)             !! fast C -> benthic C                       (mol/m2)
@@ -3479 +3452 @@
                   !! === organic nitrogen ===
                   if (jfdfate.eq.0 .and. jorgben.eq.0) then
-                     freminn  = freminn + (ffastn(ji,jj) / fthk)    !! N remineralisation in this box            (mol/m3)
+                     freminn(ji,jj)  = freminn(ji,jj) + (ffastn(ji,jj) / fse3t(ji,jj,jk))    !! N remineralisation in this box            (mol/m3)
                   elseif (jfdfate.eq.1 .and. jorgben.eq.0) then
-                     ffast2slown = ffastn(ji,jj) / fthk             !! fast N -> slow N                          (mol/m3)
-                     fslown      = fslown + ffast2slown
+                     ffast2slown(ji,jj) = ffastn(ji,jj) / fse3t(ji,jj,jk)             !! fast N -> slow N                          (mol/m3)
+                     fslown(ji,jj)      = fslown(ji,jj) + ffast2slown(ji,jj)
                   elseif (jfdfate.eq.0 .and. jorgben.eq.1) then
                      f_fbenin_n(ji,jj)  = ffastn(ji,jj)             !! fast N -> benthic N                       (mol/m2)
@@ -3491 +3464 @@
                   !! === organic iron ===
                   if (jfdfate.eq.0 .and. jorgben.eq.0) then
-                     freminfe = freminfe + (ffastfe(ji,jj) / fthk)  !! Fe remineralisation in this box           (mol/m3)
-                  elseif (jfdfate.eq.1 .and. jorgben.eq.0) then
-                     ffast2slowfe = ffastn(ji,jj) / fthk            !! fast Fe -> slow Fe                        (mol/m3)
+                     freminfe(ji,jj) = freminfe(ji,jj) + (ffastfe(ji,jj) / fse3t(ji,jj,jk))  !! Fe remineralisation in this box           (mol/m3)
+! I don't think ffast2slowfe is used - marc 10/4/17
+!                  elseif (jfdfate.eq.1 .and. jorgben.eq.0) then
+!                     ffast2slowfe(ji,jj) = ffastn(ji,jj) / fse3t(ji,jj,jk)            !! fast Fe -> slow Fe                        (mol/m3)
                   elseif (jfdfate.eq.0 .and. jorgben.eq.1) then
                      f_fbenin_fe(ji,jj) = ffastfe(ji,jj)            !! fast Fe -> benthic Fe                     (mol/m2)
@@ -3502 +3476 @@
                   !! === biogenic silicon ===
                   if (jinorgben.eq.0) then
-                     freminsi = freminsi + (ffastsi(ji,jj) / fthk)  !! Si remineralisation in this box           (mol/m3)
+                     freminsi(ji,jj) = freminsi(ji,jj) + (ffastsi(ji,jj) / fse3t(ji,jj,jk))  !! Si remineralisation in this box           (mol/m3)
                   elseif (jinorgben.eq.1) then
                      f_fbenin_si(ji,jj) = ffastsi(ji,jj)            !! fast Si -> benthic Si                     (mol/m2)
@@ -3511 +3485 @@
                   !! === biogenic calcium carbonate ===
                   if (jinorgben.eq.0) then
-                     freminca = freminca + (ffastca(ji,jj) / fthk)  !! Ca remineralisation in this box           (mol/m3) 
+                     freminca(ji,jj) = freminca(ji,jj) + (ffastca(ji,jj) / fse3t(ji,jj,jk))  !! Ca remineralisation in this box           (mol/m3) 
                   elseif (jinorgben.eq.1) then
                      f_fbenin_ca(ji,jj) = ffastca(ji,jj)            !! fast Ca -> benthic Ca                     (mol/m2)
@@ -3525 +3499 @@
                   !!----------------------------------------------------------------------
                   !!
-                  fofd_n(ji,jj)  = fofd_n(ji,jj)  + (freminn  * fthk)
-                  fofd_si(ji,jj) = fofd_si(ji,jj) + (freminsi * fthk)
-                  fofd_fe(ji,jj) = fofd_fe(ji,jj) + (freminfe * fthk)
+                  fofd_n(ji,jj)  = fofd_n(ji,jj)  + (freminn(ji,jj)  * fse3t(ji,jj,jk))
+                  fofd_si(ji,jj) = fofd_si(ji,jj) + (freminsi(ji,jj) * fse3t(ji,jj,jk))
+                  fofd_fe(ji,jj) = fofd_fe(ji,jj) + (freminfe(ji,jj) * fse3t(ji,jj,jk))
 #  if defined key_roam
-                  fofd_c(ji,jj)  = fofd_c(ji,jj)  + (freminc  * fthk)
+                  fofd_c(ji,jj)  = fofd_c(ji,jj)  + (freminc(ji,jj)  * fse3t(ji,jj,jk))
 #  endif
                endif
@@ -3539 +3513 @@
                !!
                !! update fast-sinking regeneration arrays
-               fregenfast(ji,jj)   = fregenfast(ji,jj)   + (freminn  * fthk)
-               fregenfastsi(ji,jj) = fregenfastsi(ji,jj) + (freminsi * fthk)
+               fregenfast(ji,jj)   = fregenfast(ji,jj)   + (freminn(ji,jj)  * fse3t(ji,jj,jk))
+               fregenfastsi(ji,jj) = fregenfastsi(ji,jj) + (freminsi(ji,jj) * fse3t(ji,jj,jk))
 # if defined key_roam
-               fregenfastc(ji,jj)  = fregenfastc(ji,jj)  + (freminc  * fthk)
+               fregenfastc(ji,jj)  = fregenfastc(ji,jj)  + (freminc(ji,jj)  * fse3t(ji,jj,jk))
 # endif
 
@@ -3549 +3523 @@
                !!----------------------------------------------------------------------
                !!
-               if (jk.eq.jmbathy) then
+               if (jk.eq.mbathy(ji,jj)) then
                   !!
                   !! organic components
@@ -3564 +3538 @@
                      !!
                      !! account for CaCO3 that dissolves when it shouldn't
-                     if ( fdep .le. fccd_dep ) then
+                     if ( fsdepw(ji,jj,jk) .le. fccd_dep(ji,jj) ) then
                         f_benout_lyso_ca(ji,jj) = xsedca * zn_sed_ca(ji,jj)
                      endif
@@ -3581 +3555 @@
                !!----------------------------------------------------------------------
                !! chlorophyll
-               btra(jpchn) = b0 * ( &
-                 + ((frn * fprn * zphn) - fgmipn - fgmepn - fdpn - fdpn2) * (fthetan / xxi) )
-               btra(jpchd) = b0 * ( &
-                 + ((frd * fprd * zphd) - fgmepd - fdpd - fdpd2) * (fthetad / xxi) )
+               btra(ji,jj,jpchn) = b0 * ( &
+                 + ((frn(ji,jj) * fprn(ji,jj) * zphn(ji,jj)) - fgmipn(ji,jj) - fgmepn(ji,jj) - fdpn(ji,jj) - fdpn2(ji,jj)) * (fthetan(ji,jj) / xxi) )
+               btra(ji,jj,jpchd) = b0 * ( &
+                 + ((frd(ji,jj) * fprd(ji,jj) * zphd(ji,jj)) - fgmepd(ji,jj) - fdpd(ji,jj) - fdpd2(ji,jj)) * (fthetad(ji,jj) / xxi) )
                !!
                !!----------------------------------------------------------------------
                !! phytoplankton
-               btra(jpphn) = b0 * ( &
-                 + (fprn * zphn) - fgmipn - fgmepn - fdpn - fdpn2 )
-               btra(jpphd) = b0 * ( &
-                 + (fprd * zphd) - fgmepd - fdpd - fdpd2 )
-               btra(jppds) = b0 * ( &
-                 + (fprds * zpds) - fgmepds - fdpds - fsdiss - fdpds2 )
+               btra(ji,jj,jpphn) = b0 * ( &
+                 + (fprn(ji,jj) * zphn(ji,jj)) - fgmipn(ji,jj) - fgmepn(ji,jj) - fdpn(ji,jj) - fdpn2(ji,jj) )
+               btra(ji,jj,jpphd) = b0 * ( &
+                 + (fprd(ji,jj) * zphd(ji,jj)) - fgmepd(ji,jj) - fdpd(ji,jj) - fdpd2(ji,jj) )
+               btra(ji,jj,jppds) = b0 * ( &
+                 + (fprds(ji,jj) * zpds(ji,jj)) - fgmepds(ji,jj) - fdpds(ji,jj) - fsdiss(ji,jj) - fdpds2(ji,jj) )
                !!
                !!----------------------------------------------------------------------
                !! zooplankton
-               btra(jpzmi) = b0 * ( &
-                 + fmigrow - fgmezmi - fdzmi - fdzmi2 )
-               btra(jpzme) = b0 * ( &
-                 + fmegrow - fdzme - fdzme2 )
+               btra(ji,jj,jpzmi) = b0 * ( &
+                 + fmigrow(ji,jj) - fgmezmi(ji,jj) - fdzmi(ji,jj) - fdzmi2(ji,jj) )
+               btra(ji,jj,jpzme) = b0 * ( &
+                 + fmegrow(ji,jj) - fdzme(ji,jj) - fdzme2(ji,jj) )
                !!
                !!----------------------------------------------------------------------
                !! detritus
-               btra(jpdet) = b0 * ( &
-                 + fdpn + ((1.0 - xfdfrac1) * fdpd)              &  ! mort. losses
-                 + fdzmi + ((1.0 - xfdfrac2) * fdzme)            &  ! mort. losses
-                 + ((1.0 - xbetan) * (finmi + finme))            &  ! assim. inefficiency
-                 - fgmid - fgmed - fdd                           &  ! grazing and remin.
-                 + ffast2slown )                                    ! seafloor fast->slow
+               btra(ji,jj,jpdet) = b0 * ( &
+                 + fdpn(ji,jj) + ((1.0 - xfdfrac1) * fdpd(ji,jj))              &  ! mort. losses
+                 + fdzmi(ji,jj) + ((1.0 - xfdfrac2) * fdzme(ji,jj))            &  ! mort. losses
+                 + ((1.0 - xbetan) * (finmi(ji,jj) + finme(ji,jj)))            &  ! assim. inefficiency
+                 - fgmid(ji,jj) - fgmed(ji,jj) - fdd(ji,jj)                           &  ! grazing and remin.
+                 + ffast2slown(ji,jj) )                                    ! seafloor fast->slow
                !!
                !!----------------------------------------------------------------------
                !! dissolved inorganic nitrogen nutrient
                fn_cons = 0.0  &
-                 - (fprn * zphn) - (fprd * zphd)                    ! primary production
+                 - (fprn(ji,jj) * zphn(ji,jj)) - (fprd(ji,jj) * zphd(ji,jj))                    ! primary production
                fn_prod = 0.0  &
-                 + (xphi * (fgmipn + fgmid))                     &  ! messy feeding remin.
-                 + (xphi * (fgmepn + fgmepd + fgmezmi + fgmed))  &  ! messy feeding remin.
-                 + fmiexcr + fmeexcr + fdd + freminn             &  ! excretion and remin.
-                 + fdpn2 + fdpd2 + fdzmi2 + fdzme2                  ! metab. losses
+                 + (xphi * (fgmipn(ji,jj) + fgmid(ji,jj)))                     &  ! messy feeding remin.
+                 + (xphi * (fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj)))  &  ! messy feeding remin.
+                 + fmiexcr(ji,jj) + fmeexcr(ji,jj) + fdd(ji,jj) + freminn(ji,jj)             &  ! excretion and remin.
+                 + fdpn2(ji,jj) + fdpd2(ji,jj) + fdzmi2(ji,jj) + fdzme2(ji,jj)                  ! metab. losses
                !! 
                !! riverine flux
                if ( jriver_n .gt. 0 ) then
-                  f_riv_loc_n = f_riv_n(ji,jj) * friver_dep(jk,jmbathy) / fthk
-                  fn_prod = fn_prod + f_riv_loc_n
+                  f_riv_loc_n(ji,jj) = f_riv_n(ji,jj) * friver_dep(jk,mbathy(ji,jj)) / fse3t(ji,jj,jk)
+                  fn_prod = fn_prod + f_riv_loc_n(ji,jj)
                endif
                !!  
                !! benthic remineralisation
-               if (jk.eq.jmbathy .and. jorgben.eq.1 .and. ibenthic.eq.1) then
-                  fn_prod = fn_prod + (f_benout_n(ji,jj) / fthk)
+               if (jk.eq.mbathy(ji,jj) .and. jorgben.eq.1 .and. ibenthic.eq.1) then
+                  fn_prod = fn_prod + (f_benout_n(ji,jj) / fse3t(ji,jj,jk))
                endif
                !!
-               btra(jpdin) = b0 * ( &
+               btra(ji,jj,jpdin) = b0 * ( &
                  fn_prod + fn_cons )
                !!
-               fnit_cons(ji,jj) = fnit_cons(ji,jj) + ( fthk * (  &  ! consumption of dissolved nitrogen
+               fnit_cons(ji,jj) = fnit_cons(ji,jj) + ( fse3t(ji,jj,jk) * (  &  ! consumption of dissolved nitrogen
                  fn_cons ) )
-               fnit_prod(ji,jj) = fnit_prod(ji,jj) + ( fthk * (  &  ! production of dissolved nitrogen
+               fnit_prod(ji,jj) = fnit_prod(ji,jj) + ( fse3t(ji,jj,jk) * (  &  ! production of dissolved nitrogen
                  fn_prod ) )
                !!
@@ -3643 +3617 @@
                !! dissolved silicic acid nutrient
                fs_cons = 0.0  &
-                 - (fprds * zpds)                                   ! opal production
+                 - (fprds(ji,jj) * zpds(ji,jj))                                   ! opal production
                fs_prod = 0.0  &
-                 + fsdiss                                        &  ! opal dissolution
-                 + ((1.0 - xfdfrac1) * fdpds)                    &  ! mort. loss
-                 + ((1.0 - xfdfrac3) * fgmepds)                  &  ! egestion of grazed Si
-                 + freminsi + fdpds2                                ! fast diss. and metab. losses
+                 + fsdiss(ji,jj)                                        &  ! opal dissolution
+                 + ((1.0 - xfdfrac1) * fdpds(ji,jj))                    &  ! mort. loss
+                 + ((1.0 - xfdfrac3) * fgmepds(ji,jj))                  &  ! egestion of grazed Si
+                 + freminsi(ji,jj) + fdpds2(ji,jj)                                ! fast diss. and metab. losses
                !! 
                !! riverine flux
                if ( jriver_si .gt. 0 ) then
-                  f_riv_loc_si = f_riv_si(ji,jj) * friver_dep(jk,jmbathy) / fthk
-                  fs_prod = fs_prod + f_riv_loc_si
+                  f_riv_loc_si(ji,jj) = f_riv_si(ji,jj) * friver_dep(jk,mbathy(ji,jj)) / fse3t(ji,jj,jk)
+                  fs_prod = fs_prod + f_riv_loc_si(ji,jj)
                endif
                !!  
                !! benthic remineralisation
-               if (jk.eq.jmbathy .and. jinorgben.eq.1 .and. ibenthic.eq.1) then
-                  fs_prod = fs_prod + (f_benout_si(ji,jj) / fthk)
+               if (jk.eq.mbathy(ji,jj) .and. jinorgben.eq.1 .and. ibenthic.eq.1) then
+                  fs_prod = fs_prod + (f_benout_si(ji,jj) / fse3t(ji,jj,jk))
                endif
                !!
-               btra(jpsil) = b0 * ( &
+               btra(ji,jj,jpsil) = b0 * ( &
                  fs_prod + fs_cons )
                !!
-               fsil_cons(ji,jj) = fsil_cons(ji,jj) + ( fthk * (  &  ! consumption of dissolved silicon
+               fsil_cons(ji,jj) = fsil_cons(ji,jj) + ( fse3t(ji,jj,jk) * (  &  ! consumption of dissolved silicon
                  fs_cons ) )
-               fsil_prod(ji,jj) = fsil_prod(ji,jj) + ( fthk * (  &  ! production of dissolved silicon
+               fsil_prod(ji,jj) = fsil_prod(ji,jj) + ( fse3t(ji,jj,jk) * (  &  ! production of dissolved silicon
                  fs_prod ) )
                !!
                !!----------------------------------------------------------------------
                !! dissolved "iron" nutrient
-               btra(jpfer) = b0 * ( &
-               + (xrfn * btra(jpdin)) + ffetop + ffebot - ffescav )
+               btra(ji,jj,jpfer) = b0 * ( &
+               + (xrfn * btra(ji,jj,jpdin)) + ffetop(ji,jj) + ffebot(ji,jj) - ffescav(ji,jj) )
 
 # if defined key_roam
@@ -3678 +3652 @@
                !!----------------------------------------------------------------------
                !! AXY (26/11/08): implicit detrital carbon change
-               btra(jpdtc) = b0 * ( &
-                 + (xthetapn * fdpn) + ((1.0 - xfdfrac1) * (xthetapd * fdpd))      &  ! mort. losses
-                 + (xthetazmi * fdzmi) + ((1.0 - xfdfrac2) * (xthetazme * fdzme))  &  ! mort. losses
-                 + ((1.0 - xbetac) * (ficmi + ficme))                              &  ! assim. inefficiency
-                 - fgmidc - fgmedc - fddc                                          &  ! grazing and remin.
-                 + ffast2slowc )                                                      ! seafloor fast->slow
+               btra(ji,jj,jpdtc) = b0 * ( &
+                 + (xthetapn * fdpn(ji,jj)) + ((1.0 - xfdfrac1) * (xthetapd * fdpd(ji,jj)))      &  ! mort. losses
+                 + (xthetazmi * fdzmi(ji,jj)) + ((1.0 - xfdfrac2) * (xthetazme * fdzme(ji,jj)))  &  ! mort. losses
+                 + ((1.0 - xbetac) * (ficmi(ji,jj) + ficme(ji,jj)))                              &  ! assim. inefficiency
+                 - fgmidc(ji,jj) - fgmedc(ji,jj) - fddc(ji,jj)                                          &  ! grazing and remin.
+                 + ffast2slowc(ji,jj) )                                                      ! seafloor fast->slow
                !!
                !!----------------------------------------------------------------------
                !! dissolved inorganic carbon
                fc_cons = 0.0  &
-                 - (xthetapn * fprn * zphn) - (xthetapd * fprd * zphd)                ! primary production
+                 - (xthetapn * fprn(ji,jj) * zphn(ji,jj)) - (xthetapd * fprd(ji,jj) * zphd(ji,jj))                ! primary production
                fc_prod = 0.0  &
-                 + (xthetapn * xphi * fgmipn) + (xphi * fgmidc)                    &  ! messy feeding remin
-                 + (xthetapn * xphi * fgmepn) + (xthetapd * xphi * fgmepd)         &  ! messy feeding remin
-                 + (xthetazmi * xphi * fgmezmi) + (xphi * fgmedc)                  &  ! messy feeding remin
-                 + fmiresp + fmeresp + fddc + freminc + (xthetapn * fdpn2)         &  ! resp., remin., losses
-                 + (xthetapd * fdpd2) + (xthetazmi * fdzmi2)                       &  ! losses
-                 + (xthetazme * fdzme2)                                               ! losses
+                 + (xthetapn * xphi * fgmipn(ji,jj)) + (xphi * fgmidc(ji,jj))                    &  ! messy feeding remin
+                 + (xthetapn * xphi * fgmepn(ji,jj)) + (xthetapd * xphi * fgmepd(ji,jj))         &  ! messy feeding remin
+                 + (xthetazmi * xphi * fgmezmi(ji,jj)) + (xphi * fgmedc(ji,jj))                  &  ! messy feeding remin
+                 + fmiresp(ji,jj) + fmeresp(ji,jj) + fddc(ji,jj) + freminc(ji,jj) + (xthetapn * fdpn2(ji,jj))         &  ! resp., remin., losses
+                 + (xthetapd * fdpd2(ji,jj)) + (xthetazmi * fdzmi2(ji,jj))                       &  ! losses
+                 + (xthetazme * fdzme2(ji,jj))                                               ! losses
                !! 
                !! riverine flux
                if ( jriver_c .gt. 0 ) then
-                  f_riv_loc_c = f_riv_c(ji,jj) * friver_dep(jk,jmbathy) / fthk
-                  fc_prod = fc_prod + f_riv_loc_c
+                  f_riv_loc_c(ji,jj) = f_riv_c(ji,jj) * friver_dep(jk,mbathy(ji,jj)) / fse3t(ji,jj,jk)
+                  fc_prod = fc_prod + f_riv_loc_c(ji,jj)
                endif
                !!  
                !! benthic remineralisation
-               if (jk.eq.jmbathy .and. jorgben.eq.1 .and. ibenthic.eq.1) then
-                  fc_prod = fc_prod + (f_benout_c(ji,jj) / fthk)
+               if (jk.eq.mbathy(ji,jj) .and. jorgben.eq.1 .and. ibenthic.eq.1) then
+                  fc_prod = fc_prod + (f_benout_c(ji,jj) / fse3t(ji,jj,jk))
                endif
-               if (jk.eq.jmbathy .and. jinorgben.eq.1 .and. ibenthic.eq.1) then
-                  fc_prod = fc_prod + (f_benout_ca(ji,jj) / fthk)
+               if (jk.eq.mbathy(ji,jj) .and. jinorgben.eq.1 .and. ibenthic.eq.1) then
+                  fc_prod = fc_prod + (f_benout_ca(ji,jj) / fse3t(ji,jj,jk))
                endif
                !!
@@ -3715 +3689 @@
                !!
                !! CaCO3
-               fc_prod = fc_prod - ftempca + freminca
+               fc_prod = fc_prod - ftempca(ji,jj) + freminca(ji,jj)
                !! 
                !! riverine flux
@@ -3722 +3696 @@
                endif
                !!
-               btra(jpdic) = b0 * ( &
+               btra(ji,jj,jpdic) = b0 * ( &
                  fc_prod + fc_cons )
                !!
-               fcar_cons(ji,jj) = fcar_cons(ji,jj) + ( fthk * (  &  ! consumption of dissolved carbon
+               fcar_cons(ji,jj) = fcar_cons(ji,jj) + ( fse3t(ji,jj,jk) * (  &  ! consumption of dissolved carbon
                  fc_cons ) )
-               fcar_prod(ji,jj) = fcar_prod(ji,jj) + ( fthk * (  &  ! production of dissolved carbon
+               fcar_prod(ji,jj) = fcar_prod(ji,jj) + ( fse3t(ji,jj,jk) * (  &  ! production of dissolved carbon
                  fc_prod ) )
                !!
@@ -3733 +3707 @@
                !! alkalinity
                fa_prod = 0.0  &
-                 + (2.0 * freminca)                                                   ! CaCO3 dissolution
+                 + (2.0 * freminca(ji,jj))                                                   ! CaCO3 dissolution
                fa_cons = 0.0  &
-                 - (2.0 * ftempca)                                                    ! CaCO3 production
+                 - (2.0 * ftempca(ji,jj))                                                    ! CaCO3 production
                !! 
                !! riverine flux
                if ( jriver_alk .gt. 0 ) then
-                  f_riv_loc_alk = f_riv_alk(ji,jj) * friver_dep(jk,jmbathy) / fthk
-                  fa_prod = fa_prod + f_riv_loc_alk
+                  f_riv_loc_alk(ji,jj) = f_riv_alk(ji,jj) * friver_dep(jk,mbathy(ji,jj)) / fse3t(ji,jj,jk)
+                  fa_prod = fa_prod + f_riv_loc_alk(ji,jj)
                endif
                !!  
                !! benthic remineralisation
-               if (jk.eq.jmbathy .and. jinorgben.eq.1 .and. ibenthic.eq.1) then
-                  fa_prod = fa_prod + (2.0 * f_benout_ca(ji,jj) / fthk)
+               if (jk.eq.mbathy(ji,jj) .and. jinorgben.eq.1 .and. ibenthic.eq.1) then
+                  fa_prod = fa_prod + (2.0 * f_benout_ca(ji,jj) / fse3t(ji,jj,jk))
                endif
                !!
-               btra(jpalk) = b0 * ( &
+               btra(ji,jj,jpalk) = b0 * ( &
                  fa_prod + fa_cons )
                !!
                !!----------------------------------------------------------------------
                !! oxygen (has protection at low O2 concentrations; OCMIP-2 style)
-               fo2_prod = 0.0 &
-                 + (xthetanit * fprn * zphn)                                      & ! Pn primary production, N
-                 + (xthetanit * fprd * zphd)                                      & ! Pd primary production, N
-                 + (xthetarem * xthetapn * fprn * zphn)                           & ! Pn primary production, C
-                 + (xthetarem * xthetapd * fprd * zphd)                             ! Pd primary production, C
-               fo2_ncons = 0.0 &
-                 - (xthetanit * xphi * fgmipn)                                    & ! Pn messy feeding remin., N
-                 - (xthetanit * xphi * fgmid)                                     & ! D  messy feeding remin., N
-                 - (xthetanit * xphi * fgmepn)                                    & ! Pn messy feeding remin., N
-                 - (xthetanit * xphi * fgmepd)                                    & ! Pd messy feeding remin., N
-                 - (xthetanit * xphi * fgmezmi)                                   & ! Zi messy feeding remin., N
-                 - (xthetanit * xphi * fgmed)                                     & ! D  messy feeding remin., N
-                 - (xthetanit * fmiexcr)                                          & ! microzoo excretion, N
-                 - (xthetanit * fmeexcr)                                          & ! mesozoo  excretion, N
-                 - (xthetanit * fdd)                                              & ! slow detritus remin., N 
-                 - (xthetanit * freminn)                                          & ! fast detritus remin., N
-                 - (xthetanit * fdpn2)                                            & ! Pn  losses, N
-                 - (xthetanit * fdpd2)                                            & ! Pd  losses, N
-                 - (xthetanit * fdzmi2)                                           & ! Zmi losses, N
-                 - (xthetanit * fdzme2)                                             ! Zme losses, N
+               fo2_prod(ji,jj) = 0.0 &
+                 + (xthetanit * fprn(ji,jj) * zphn(ji,jj))                                      & ! Pn primary production, N
+                 + (xthetanit * fprd(ji,jj) * zphd(ji,jj))                                      & ! Pd primary production, N
+                 + (xthetarem * xthetapn * fprn(ji,jj) * zphn(ji,jj))                           & ! Pn primary production, C
+                 + (xthetarem * xthetapd * fprd(ji,jj) * zphd(ji,jj))                             ! Pd primary production, C
+               fo2_ncons(ji,jj) = 0.0 &
+                 - (xthetanit * xphi * fgmipn(ji,jj))                                    & ! Pn messy feeding remin., N
+                 - (xthetanit * xphi * fgmid(ji,jj))                                     & ! D  messy feeding remin., N
+                 - (xthetanit * xphi * fgmepn(ji,jj))                                    & ! Pn messy feeding remin., N
+                 - (xthetanit * xphi * fgmepd(ji,jj))                                    & ! Pd messy feeding remin., N
+                 - (xthetanit * xphi * fgmezmi(ji,jj))                                   & ! Zi messy feeding remin., N
+                 - (xthetanit * xphi * fgmed(ji,jj))                                     & ! D  messy feeding remin., N
+                 - (xthetanit * fmiexcr(ji,jj))                                          & ! microzoo excretion, N
+                 - (xthetanit * fmeexcr(ji,jj))                                          & ! mesozoo  excretion, N
+                 - (xthetanit * fdd(ji,jj))                                              & ! slow detritus remin., N 
+                 - (xthetanit * freminn(ji,jj))                                          & ! fast detritus remin., N
+                 - (xthetanit * fdpn2(ji,jj))                                            & ! Pn  losses, N
+                 - (xthetanit * fdpd2(ji,jj))                                            & ! Pd  losses, N
+                 - (xthetanit * fdzmi2(ji,jj))                                           & ! Zmi losses, N
+                 - (xthetanit * fdzme2(ji,jj))                                             ! Zme losses, N
                !!  
                !! benthic remineralisation
-               if (jk.eq.jmbathy .and. jorgben.eq.1 .and. ibenthic.eq.1) then
-                  fo2_ncons = fo2_ncons - (xthetanit * f_benout_n(ji,jj) / fthk)
+               if (jk.eq.mbathy(ji,jj) .and. jorgben.eq.1 .and. ibenthic.eq.1) then
+                  fo2_ncons(ji,jj) = fo2_ncons(ji,jj) - (xthetanit * f_benout_n(ji,jj) / fse3t(ji,jj,jk))
                endif
-               fo2_ccons = 0.0 &
-                 - (xthetarem * xthetapn * xphi * fgmipn)                         & ! Pn messy feeding remin., C
-                 - (xthetarem * xphi * fgmidc)                                    & ! D  messy feeding remin., C
-                 - (xthetarem * xthetapn * xphi * fgmepn)                         & ! Pn messy feeding remin., C
-                 - (xthetarem * xthetapd * xphi * fgmepd)                         & ! Pd messy feeding remin., C
-                 - (xthetarem * xthetazmi * xphi * fgmezmi)                       & ! Zi messy feeding remin., C
-                 - (xthetarem * xphi * fgmedc)                                    & ! D  messy feeding remin., C
-                 - (xthetarem * fmiresp)                                          & ! microzoo respiration, C
-                 - (xthetarem * fmeresp)                                          & ! mesozoo  respiration, C
-                 - (xthetarem * fddc)                                             & ! slow detritus remin., C
-                 - (xthetarem * freminc)                                          & ! fast detritus remin., C
-                 - (xthetarem * xthetapn * fdpn2)                                 & ! Pn  losses, C
-                 - (xthetarem * xthetapd * fdpd2)                                 & ! Pd  losses, C
-                 - (xthetarem * xthetazmi * fdzmi2)                               & ! Zmi losses, C
-                 - (xthetarem * xthetazme * fdzme2)                                 ! Zme losses, C
+               fo2_ccons(ji,jj) = 0.0 &
+                 - (xthetarem * xthetapn * xphi * fgmipn(ji,jj))                         & ! Pn messy feeding remin., C
+                 - (xthetarem * xphi * fgmidc(ji,jj))                                    & ! D  messy feeding remin., C
+                 - (xthetarem * xthetapn * xphi * fgmepn(ji,jj))                         & ! Pn messy feeding remin., C
+                 - (xthetarem * xthetapd * xphi * fgmepd(ji,jj))                         & ! Pd messy feeding remin., C
+                 - (xthetarem * xthetazmi * xphi * fgmezmi(ji,jj))                       & ! Zi messy feeding remin., C
+                 - (xthetarem * xphi * fgmedc(ji,jj))                                    & ! D  messy feeding remin., C
+                 - (xthetarem * fmiresp(ji,jj))                                          & ! microzoo respiration, C
+                 - (xthetarem * fmeresp(ji,jj))                                          & ! mesozoo  respiration, C
+                 - (xthetarem * fddc(ji,jj))                                             & ! slow detritus remin., C
+                 - (xthetarem * freminc(ji,jj))                                          & ! fast detritus remin., C
+                 - (xthetarem * xthetapn * fdpn2(ji,jj))                                 & ! Pn  losses, C
+                 - (xthetarem * xthetapd * fdpd2(ji,jj))                                 & ! Pd  losses, C
+                 - (xthetarem * xthetazmi * fdzmi2(ji,jj))                               & ! Zmi losses, C
+                 - (xthetarem * xthetazme * fdzme2(ji,jj))                                 ! Zme losses, C
                !!  
                !! benthic remineralisation
-               if (jk.eq.jmbathy .and. jorgben.eq.1 .and. ibenthic.eq.1) then
-                  fo2_ccons = fo2_ccons - (xthetarem * f_benout_c(ji,jj) / fthk)
+               if (jk.eq.mbathy(ji,jj) .and. jorgben.eq.1 .and. ibenthic.eq.1) then
+                  fo2_ccons(ji,jj) = fo2_ccons(ji,jj) - (xthetarem * f_benout_c(ji,jj) / fse3t(ji,jj,jk))
                endif
-               fo2_cons = fo2_ncons + fo2_ccons
+               fo2_cons(ji,jj) = fo2_ncons(ji,jj) + fo2_ccons(ji,jj)
                !!
                !! is this a suboxic zone?
-               if (zoxy.lt.xo2min) then  ! deficient O2; production fluxes only
-                  btra(jpoxy) = b0 * ( &
-                    fo2_prod )
-                  foxy_prod(ji,jj) = foxy_prod(ji,jj) + ( fthk * fo2_prod )
-                  foxy_anox(ji,jj) = foxy_anox(ji,jj) + ( fthk * fo2_cons )
+               if (zoxy(ji,jj).lt.xo2min) then  ! deficient O2; production fluxes only
+                  btra(ji,jj,jpoxy) = b0 * ( &
+                    fo2_prod(ji,jj) )
+                  foxy_prod(ji,jj) = foxy_prod(ji,jj) + ( fse3t(ji,jj,jk) * fo2_prod(ji,jj) )
+                  foxy_anox(ji,jj) = foxy_anox(ji,jj) + ( fse3t(ji,jj,jk) * fo2_cons(ji,jj) )
                else                      ! sufficient O2; production + consumption fluxes
-                  btra(jpoxy) = b0 * ( &
-                    fo2_prod + fo2_cons )
-                  foxy_prod(ji,jj) = foxy_prod(ji,jj) + ( fthk * fo2_prod )
-                  foxy_cons(ji,jj) = foxy_cons(ji,jj) + ( fthk * fo2_cons )
+                  btra(ji,jj,jpoxy) = b0 * ( &
+                    fo2_prod(ji,jj) + fo2_cons(ji,jj) )
+                  foxy_prod(ji,jj) = foxy_prod(ji,jj) + ( fse3t(ji,jj,jk) * fo2_prod(ji,jj) )
+                  foxy_cons(ji,jj) = foxy_cons(ji,jj) + ( fse3t(ji,jj,jk) * fo2_cons(ji,jj) )
                endif
                !!
@@ -3817 +3791 @@
                   !!
                   !! CO2 flux
-                  btra(jpdic) = btra(jpdic) + (b0 * f_co2flux)
+                  btra(ji,jj,jpdic) = btra(ji,jj,jpdic) + (b0 * f_co2flux(ji,jj))
                   !!
                   !! O2 flux (mol/m3/s -> mmol/m3/d)
-                  btra(jpoxy) = btra(jpoxy) + (b0 * f_o2flux)
+                  btra(ji,jj,jpoxy) = btra(ji,jj,jpoxy) + (b0 * f_o2flux(ji,jj))
                endif
 # endif
@@ -3828 +3802 @@
                if (idf.eq.1.AND.idfval.eq.1) then
                   IF (lwp) write (numout,*) '------------------------------'
-                  IF (lwp) write (numout,*) 'btra(jpchn)(',jk,')  = ', btra(jpchn)
-                  IF (lwp) write (numout,*) 'btra(jpchd)(',jk,')  = ', btra(jpchd)
-                  IF (lwp) write (numout,*) 'btra(jpphn)(',jk,')  = ', btra(jpphn)
-                  IF (lwp) write (numout,*) 'btra(jpphd)(',jk,')  = ', btra(jpphd)
-                  IF (lwp) write (numout,*) 'btra(jppds)(',jk,')  = ', btra(jppds)
-                  IF (lwp) write (numout,*) 'btra(jpzmi)(',jk,')  = ', btra(jpzmi)
-                  IF (lwp) write (numout,*) 'btra(jpzme)(',jk,')  = ', btra(jpzme)
-                  IF (lwp) write (numout,*) 'btra(jpdet)(',jk,')  = ', btra(jpdet)
-                  IF (lwp) write (numout,*) 'btra(jpdin)(',jk,')  = ', btra(jpdin)
-                  IF (lwp) write (numout,*) 'btra(jpsil)(',jk,')  = ', btra(jpsil)
-                  IF (lwp) write (numout,*) 'btra(jpfer)(',jk,')  = ', btra(jpfer)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpchn)(',jk,')  = ', btra(ji,jj,jpchn)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpchd)(',jk,')  = ', btra(ji,jj,jpchd)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpphn)(',jk,')  = ', btra(ji,jj,jpphn)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpphd)(',jk,')  = ', btra(ji,jj,jpphd)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jppds)(',jk,')  = ', btra(ji,jj,jppds)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpzmi)(',jk,')  = ', btra(ji,jj,jpzmi)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpzme)(',jk,')  = ', btra(ji,jj,jpzme)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdet)(',jk,')  = ', btra(ji,jj,jpdet)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdin)(',jk,')  = ', btra(ji,jj,jpdin)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpsil)(',jk,')  = ', btra(ji,jj,jpsil)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpfer)(',jk,')  = ', btra(ji,jj,jpfer)
 #  if defined key_roam
-                  IF (lwp) write (numout,*) 'btra(jpdtc)(',jk,')  = ', btra(jpdtc)
-                  IF (lwp) write (numout,*) 'btra(jpdic)(',jk,')  = ', btra(jpdic)
-                  IF (lwp) write (numout,*) 'btra(jpalk)(',jk,')  = ', btra(jpalk)
-                  IF (lwp) write (numout,*) 'btra(jpoxy)(',jk,')  = ', btra(jpoxy)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdtc)(',jk,')  = ', btra(ji,jj,jpdtc)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpdic)(',jk,')  = ', btra(ji,jj,jpdic)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpalk)(',jk,')  = ', btra(ji,jj,jpalk)
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpoxy)(',jk,')  = ', btra(ji,jj,jpoxy)
 #  endif
                endif
@@ -3854 +3828 @@
                !! === nitrogen ===
                fflx_n(ji,jj)  = fflx_n(ji,jj)  + &
-                  fthk * ( btra(jpphn) + btra(jpphd) + btra(jpzmi) + btra(jpzme) + btra(jpdet) + btra(jpdin) )
+                  fse3t(ji,jj,jk) * ( btra(ji,jj,jpphn) + btra(ji,jj,jpphd) + btra(ji,jj,jpzmi) + btra(ji,jj,jpzme) + btra(ji,jj,jpdet) + btra(ji,jj,jpdin) )
                !! === silicon ===
                fflx_si(ji,jj) = fflx_si(ji,jj) + &
-                  fthk * ( btra(jppds) + btra(jpsil) )
+                  fse3t(ji,jj,jk) * ( btra(ji,jj,jppds) + btra(ji,jj,jpsil) )
                !! === iron ===
                fflx_fe(ji,jj) = fflx_fe(ji,jj) + &
-                  fthk * ( ( xrfn * ( btra(jpphn) + btra(jpphd) + btra(jpzmi) + btra(jpzme) + btra(jpdet)) ) + btra(jpfer) )
+                  fse3t(ji,jj,jk) * ( ( xrfn * ( btra(ji,jj,jpphn) + btra(ji,jj,jpphd) + btra(ji,jj,jpzmi) + btra(ji,jj,jpzme) + btra(ji,jj,jpdet)) ) + btra(ji,jj,jpfer) )
 # if defined key_roam
                !! === carbon ===
                fflx_c(ji,jj)  = fflx_c(ji,jj)  + &
-                  fthk * ( (xthetapn * btra(jpphn)) + (xthetapd * btra(jpphd)) + &
-                  (xthetazmi * btra(jpzmi)) + (xthetazme * btra(jpzme)) + btra(jpdtc) + btra(jpdic) )
+                  fse3t(ji,jj,jk) * ( (xthetapn * btra(ji,jj,jpphn)) + (xthetapd * btra(ji,jj,jpphd)) + &
+                  (xthetazmi * btra(ji,jj,jpzmi)) + (xthetazme * btra(ji,jj,jpzme)) + btra(ji,jj,jpdtc) + btra(ji,jj,jpdic) )
                !! === alkalinity ===
                fflx_a(ji,jj)  = fflx_a(ji,jj)  + &
-                  fthk * ( btra(jpalk) )
+                  fse3t(ji,jj,jk) * ( btra(ji,jj,jpalk) )
                !! === oxygen ===
                fflx_o2(ji,jj) = fflx_o2(ji,jj) + &
-                  fthk * ( btra(jpoxy) )
+                  fse3t(ji,jj,jk) * ( btra(ji,jj,jpoxy) )
 # endif
 
@@ -3895 +3869 @@
 # endif               
                if (ibio_switch.eq.1) then
-                  tra(ji,jj,jk,jpchn) = tra(ji,jj,jk,jpchn) + (btra(jpchn) / 86400.)
-                  tra(ji,jj,jk,jpchd) = tra(ji,jj,jk,jpchd) + (btra(jpchd) / 86400.)
-                  tra(ji,jj,jk,jpphn) = tra(ji,jj,jk,jpphn) + (btra(jpphn) / 86400.)
-                  tra(ji,jj,jk,jpphd) = tra(ji,jj,jk,jpphd) + (btra(jpphd) / 86400.)
-                  tra(ji,jj,jk,jppds) = tra(ji,jj,jk,jppds) + (btra(jppds) / 86400.)
-                  tra(ji,jj,jk,jpzmi) = tra(ji,jj,jk,jpzmi) + (btra(jpzmi) / 86400.)
-                  tra(ji,jj,jk,jpzme) = tra(ji,jj,jk,jpzme) + (btra(jpzme) / 86400.)
-                  tra(ji,jj,jk,jpdet) = tra(ji,jj,jk,jpdet) + (btra(jpdet) / 86400.)
-                  tra(ji,jj,jk,jpdin) = tra(ji,jj,jk,jpdin) + (btra(jpdin) / 86400.)
-                  tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) + (btra(jpsil) / 86400.)
-                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + (btra(jpfer) / 86400.)
+                  tra(ji,jj,jk,jpchn) = tra(ji,jj,jk,jpchn) + (btra(ji,jj,jpchn) / 86400.)
+                  tra(ji,jj,jk,jpchd) = tra(ji,jj,jk,jpchd) + (btra(ji,jj,jpchd) / 86400.)
+                  tra(ji,jj,jk,jpphn) = tra(ji,jj,jk,jpphn) + (btra(ji,jj,jpphn) / 86400.)
+                  tra(ji,jj,jk,jpphd) = tra(ji,jj,jk,jpphd) + (btra(ji,jj,jpphd) / 86400.)
+                  tra(ji,jj,jk,jppds) = tra(ji,jj,jk,jppds) + (btra(ji,jj,jppds) / 86400.)
+                  tra(ji,jj,jk,jpzmi) = tra(ji,jj,jk,jpzmi) + (btra(ji,jj,jpzmi) / 86400.)
+                  tra(ji,jj,jk,jpzme) = tra(ji,jj,jk,jpzme) + (btra(ji,jj,jpzme) / 86400.)
+                  tra(ji,jj,jk,jpdet) = tra(ji,jj,jk,jpdet) + (btra(ji,jj,jpdet) / 86400.)
+                  tra(ji,jj,jk,jpdin) = tra(ji,jj,jk,jpdin) + (btra(ji,jj,jpdin) / 86400.)
+                  tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) + (btra(ji,jj,jpsil) / 86400.)
+                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + (btra(ji,jj,jpfer) / 86400.)
 # if defined key_roam
-                  tra(ji,jj,jk,jpdtc) = tra(ji,jj,jk,jpdtc) + (btra(jpdtc) / 86400.)
-                  tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + (btra(jpdic) / 86400.)
-                  tra(ji,jj,jk,jpalk) = tra(ji,jj,jk,jpalk) + (btra(jpalk) / 86400.)
-                  tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + (btra(jpoxy) / 86400.)
+                  tra(ji,jj,jk,jpdtc) = tra(ji,jj,jk,jpdtc) + (btra(ji,jj,jpdtc) / 86400.)
+                  tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + (btra(ji,jj,jpdic) / 86400.)
+                  tra(ji,jj,jk,jpalk) = tra(ji,jj,jk,jpalk) + (btra(ji,jj,jpalk) / 86400.)
+                  tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + (btra(ji,jj,jpoxy) / 86400.)
 # endif
                endif               
@@ -3916 +3890 @@
                !! AXY (18/11/16): CMIP6 diagnostics
                IF( med_diag%FBDDTALK%dgsave )  THEN
-                  fbddtalk(ji,jj)  =  fbddtalk(ji,jj)  + (btra(jpalk) * fthk)
+                  fbddtalk(ji,jj)  =  fbddtalk(ji,jj)  + (btra(ji,jj,jpalk) * fse3t(ji,jj,jk))
                ENDIF
                IF( med_diag%FBDDTDIC%dgsave )  THEN
-                  fbddtdic(ji,jj)  =  fbddtdic(ji,jj)  + (btra(jpdic) * fthk)
+                  fbddtdic(ji,jj)  =  fbddtdic(ji,jj)  + (btra(ji,jj,jpdic) * fse3t(ji,jj,jk))
                ENDIF
                IF( med_diag%FBDDTDIFE%dgsave ) THEN
-                  fbddtdife(ji,jj) =  fbddtdife(ji,jj) + (btra(jpfer) * fthk)
+                  fbddtdife(ji,jj) =  fbddtdife(ji,jj) + (btra(ji,jj,jpfer) * fse3t(ji,jj,jk))
                ENDIF
                IF( med_diag%FBDDTDIN%dgsave )  THEN
-                  fbddtdin(ji,jj)  =  fbddtdin(ji,jj)  + (btra(jpdin) * fthk)
+                  fbddtdin(ji,jj)  =  fbddtdin(ji,jj)  + (btra(ji,jj,jpdin) * fse3t(ji,jj,jk))
                ENDIF
                IF( med_diag%FBDDTDISI%dgsave ) THEN
-                  fbddtdisi(ji,jj) =  fbddtdisi(ji,jj) + (btra(jpsil) * fthk)
+                  fbddtdisi(ji,jj) =  fbddtdisi(ji,jj) + (btra(ji,jj,jpsil) * fse3t(ji,jj,jk))
                ENDIF
           !!
                IF( med_diag%BDDTALK3%dgsave )  THEN
-                  bddtalk3(ji,jj,jk)  =  btra(jpalk)
+                  bddtalk3(ji,jj,jk)  =  btra(ji,jj,jpalk)
                ENDIF
                IF( med_diag%BDDTDIC3%dgsave )  THEN
-                  bddtdic3(ji,jj,jk)  =  btra(jpdic)
+                  bddtdic3(ji,jj,jk)  =  btra(ji,jj,jpdic)
                ENDIF
                IF( med_diag%BDDTDIFE3%dgsave ) THEN
-                  bddtdife3(ji,jj,jk) =  btra(jpfer)
+                  bddtdife3(ji,jj,jk) =  btra(ji,jj,jpfer)
                ENDIF
                IF( med_diag%BDDTDIN3%dgsave )  THEN
-                  bddtdin3(ji,jj,jk)  =  btra(jpdin)
+                  bddtdin3(ji,jj,jk)  =  btra(ji,jj,jpdin)
                ENDIF
                IF( med_diag%BDDTDISI3%dgsave ) THEN
-                  bddtdisi3(ji,jj,jk) =  btra(jpsil)
+                  bddtdisi3(ji,jj,jk) =  btra(ji,jj,jpsil)
                ENDIF
 
@@ -3960 +3934 @@
                !!
                DO jn = 1,jptra
-                  fq0 = btra(jn)
+                  fq0 = btra(ji,jj,jn)
                   !! AXY (30/01/14): "isnan" problem on HECTOR
                   !! if (fq0 /= fq0 ) then
                   if ( ieee_is_nan( fq0 ) ) then
                      !! there's a NaN here
-                     if (lwp) write(numout,*) 'NAN detected in btra(', ji, ',', &
+                     if (lwp) write(numout,*) 'NAN detected in btra(ji,jj,', ji, ',', &
                      & jj, ',', jk, ',', jn, ') at time', kt
            CALL ctl_stop( 'trcbio_medusa, NAN in btra field' )
@@ -3993 +3967 @@
                !!----------------------------------------------------------------------
                !!
-               !! fnit0 = btra(jpphn) + btra(jpphd) + btra(jpzmi) + btra(jpzme) + btra(jpdet) + btra(jpdin)  ! + ftempn
-               !! fsil0 = btra(jppds) + btra(jpsil)                              ! + ftempsi
-               !! ffer0 = (xrfn * fnit0) + btra(jpfer)
+               !! fnit0 = btra(ji,jj,jpphn) + btra(ji,jj,jpphd) + btra(ji,jj,jpzmi) + btra(ji,jj,jpzme) + btra(ji,jj,jpdet) + btra(ji,jj,jpdin)  ! + ftempn(ji,jj)
+               !! fsil0 = btra(ji,jj,jppds) + btra(ji,jj,jpsil)                              ! + ftempsi(ji,jj)
+               !! ffer0 = (xrfn * fnit0) + btra(ji,jj,jpfer)
 # if defined key_roam
                !! fcar0 = 0.
@@ -4027 +4001 @@
                !!
                IF(lwp) WRITE(numout,*) ' MEDUSA does not support key_trc_diabio'
-               !!               trbio(ji,jj,jk, 1) = fprn
+               !!               trbio(ji,jj,jk, 1) = fprn(ji,jj)
 # endif
 
@@ -4041 +4015 @@
 #   endif
                   IF ( med_diag%PRN%dgsave ) THEN
-                      fprn2d(ji,jj) = fprn2d(ji,jj) + (fprn  * zphn * fthk) 
+                      fprn2d(ji,jj) = fprn2d(ji,jj) + (fprn(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF ( med_diag%MPN%dgsave ) THEN
-                      fdpn2d(ji,jj) = fdpn2d(ji,jj) + (fdpn         * fthk)
+                      fdpn2d(ji,jj) = fdpn2d(ji,jj) + (fdpn(ji,jj)         * fse3t(ji,jj,jk))
                   ENDIF
                   IF ( med_diag%PRD%dgsave ) THEN
-                      fprd2d(ji,jj) = fprd2d(ji,jj) + (fprd  * zphd * fthk)
+                      fprd2d(ji,jj) = fprd2d(ji,jj) + (fprd(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%MPD%dgsave ) THEN
-                      fdpd2d(ji,jj) = fdpd2d(ji,jj) + (fdpd         * fthk) 
+                      fdpd2d(ji,jj) = fdpd2d(ji,jj) + (fdpd(ji,jj)         * fse3t(ji,jj,jk)) 
                   ENDIF
                   !  IF( med_diag%DSED%dgsave ) THEN
@@ -4056 +4030 @@
                   !  ENDIF
                   IF( med_diag%OPAL%dgsave ) THEN
-                      fprds2d(ji,jj) = fprds2d(ji,jj) + (fprds * zpds * fthk) 
+                      fprds2d(ji,jj) = fprds2d(ji,jj) + (fprds(ji,jj) * zpds(ji,jj) * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF( med_diag%OPALDISS%dgsave ) THEN
-                      fsdiss2d(ji,jj) = fsdiss2d(ji,jj) + (fsdiss  * fthk)  
+                      fsdiss2d(ji,jj) = fsdiss2d(ji,jj) + (fsdiss(ji,jj)  * fse3t(ji,jj,jk))  
                   ENDIF
                   IF( med_diag%GMIPn%dgsave ) THEN
-                      fgmipn2d(ji,jj) = fgmipn2d(ji,jj) + (fgmipn  * fthk) 
+                      fgmipn2d(ji,jj) = fgmipn2d(ji,jj) + (fgmipn(ji,jj)  * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF( med_diag%GMID%dgsave ) THEN
-                      fgmid2d(ji,jj) = fgmid2d(ji,jj) + (fgmid   * fthk) 
+                      fgmid2d(ji,jj) = fgmid2d(ji,jj) + (fgmid(ji,jj)   * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF( med_diag%MZMI%dgsave ) THEN
-                      fdzmi2d(ji,jj) = fdzmi2d(ji,jj) + (fdzmi   * fthk) 
+                      fdzmi2d(ji,jj) = fdzmi2d(ji,jj) + (fdzmi(ji,jj)   * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF( med_diag%GMEPN%dgsave ) THEN
-                      fgmepn2d(ji,jj) = fgmepn2d(ji,jj) + (fgmepn  * fthk)
+                      fgmepn2d(ji,jj) = fgmepn2d(ji,jj) + (fgmepn(ji,jj)  * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%GMEPD%dgsave ) THEN
-                      fgmepd2d(ji,jj) = fgmepd2d(ji,jj) + (fgmepd  * fthk) 
+                      fgmepd2d(ji,jj) = fgmepd2d(ji,jj) + (fgmepd(ji,jj)  * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF( med_diag%GMEZMI%dgsave ) THEN
-                      fgmezmi2d(ji,jj) = fgmezmi2d(ji,jj) + (fgmezmi * fthk) 
+                      fgmezmi2d(ji,jj) = fgmezmi2d(ji,jj) + (fgmezmi(ji,jj) * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF( med_diag%GMED%dgsave ) THEN
-                      fgmed2d(ji,jj) = fgmed2d(ji,jj) + (fgmed   * fthk) 
+                      fgmed2d(ji,jj) = fgmed2d(ji,jj) + (fgmed(ji,jj)   * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF( med_diag%MZME%dgsave ) THEN
-                      fdzme2d(ji,jj) = fdzme2d(ji,jj) + (fdzme   * fthk) 
+                      fdzme2d(ji,jj) = fdzme2d(ji,jj) + (fdzme(ji,jj)   * fse3t(ji,jj,jk)) 
                   ENDIF
                   !  IF( med_diag%DEXP%dgsave ) THEN
@@ -4089 +4063 @@
                   !  ENDIF
                   IF( med_diag%DETN%dgsave ) THEN
-                      fslown2d(ji,jj) = fslown2d(ji,jj) + (fslown  * fthk)  
+                      fslown2d(ji,jj) = fslown2d(ji,jj) + (fslown(ji,jj)  * fse3t(ji,jj,jk))  
                   ENDIF
                   IF( med_diag%MDET%dgsave ) THEN
-                      fdd2d(ji,jj) = fdd2d(ji,jj) + (fdd     * fthk) 
+                      fdd2d(ji,jj) = fdd2d(ji,jj) + (fdd(ji,jj)     * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF( med_diag%AEOLIAN%dgsave ) THEN
-                      ffetop2d(ji,jj) = ffetop2d(ji,jj) + (ffetop  * fthk) 
+                      ffetop2d(ji,jj) = ffetop2d(ji,jj) + (ffetop(ji,jj)  * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF( med_diag%BENTHIC%dgsave ) THEN
-                      ffebot2d(ji,jj) = ffebot2d(ji,jj) + (ffebot  * fthk) 
+                      ffebot2d(ji,jj) = ffebot2d(ji,jj) + (ffebot(ji,jj)  * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF( med_diag%SCAVENGE%dgsave ) THEN
-                      ffescav2d(ji,jj) = ffescav2d(ji,jj) + (ffescav * fthk)  
+                      ffescav2d(ji,jj) = ffescav2d(ji,jj) + (ffescav(ji,jj) * fse3t(ji,jj,jk))  
                   ENDIF
                   IF( med_diag%PN_JLIM%dgsave ) THEN
-                      ! fjln2d(ji,jj) = fjln2d(ji,jj) + (fjln  * zphn * fthk) 
-                      fjln2d(ji,jj) = fjln2d(ji,jj) + (fjlim_pn * zphn * fthk) 
+                      ! fjln2d(ji,jj) = fjln2d(ji,jj) + (fjln(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk)) 
+                      fjln2d(ji,jj) = fjln2d(ji,jj) + (fjlim_pn(ji,jj) * zphn(ji,jj) * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF( med_diag%PN_NLIM%dgsave ) THEN
-                      fnln2d(ji,jj) = fnln2d(ji,jj) + (fnln  * zphn * fthk) 
+                      fnln2d(ji,jj) = fnln2d(ji,jj) + (fnln(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF( med_diag%PN_FELIM%dgsave ) THEN
-                      ffln2d(ji,jj) = ffln2d(ji,jj) + (ffln  * zphn * fthk) 
+                      ffln2d(ji,jj) = ffln2d(ji,jj) + (ffln(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF( med_diag%PD_JLIM%dgsave ) THEN
-                      ! fjld2d(ji,jj) = fjld2d(ji,jj) + (fjld  * zphd * fthk) 
-                      fjld2d(ji,jj) = fjld2d(ji,jj) + (fjlim_pd * zphd * fthk) 
+                      ! fjld2d(ji,jj) = fjld2d(ji,jj) + (fjld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk)) 
+                      fjld2d(ji,jj) = fjld2d(ji,jj) + (fjlim_pd(ji,jj) * zphd(ji,jj) * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF( med_diag%PD_NLIM%dgsave ) THEN
-                      fnld2d(ji,jj) = fnld2d(ji,jj) + (fnld  * zphd * fthk) 
+                      fnld2d(ji,jj) = fnld2d(ji,jj) + (fnld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF( med_diag%PD_FELIM%dgsave ) THEN
-                      ffld2d(ji,jj) = ffld2d(ji,jj) + (ffld  * zphd * fthk) 
+                      ffld2d(ji,jj) = ffld2d(ji,jj) + (ffld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF( med_diag%PD_SILIM%dgsave ) THEN
-                      fsld2d2(ji,jj) = fsld2d2(ji,jj) + (fsld2 * zphd * fthk) 
+                      fsld2d2(ji,jj) = fsld2d2(ji,jj) + (fsld2(ji,jj) * zphd(ji,jj) * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF( med_diag%PDSILIM2%dgsave ) THEN
-                      fsld2d(ji,jj) = fsld2d(ji,jj) + (fsld  * zphd * fthk)
+                      fsld2d(ji,jj) = fsld2d(ji,jj) + (fsld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   !! 
                   IF( med_diag%TOTREG_N%dgsave ) THEN
-                      fregen2d(ji,jj) = fregen2d(ji,jj) + fregen
+                      fregen2d(ji,jj) = fregen2d(ji,jj) + fregen(ji,jj)
                   ENDIF
                   IF( med_diag%TOTRG_SI%dgsave ) THEN
-                      fregensi2d(ji,jj) = fregensi2d(ji,jj) + fregensi
+                      fregensi2d(ji,jj) = fregensi2d(ji,jj) + fregensi(ji,jj)
                   ENDIF
                   !! 
                   IF( med_diag%FASTN%dgsave ) THEN
-                      ftempn2d(ji,jj) = ftempn2d(ji,jj) + (ftempn  * fthk)
+                      ftempn2d(ji,jj) = ftempn2d(ji,jj) + (ftempn(ji,jj)  * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%FASTSI%dgsave ) THEN
-                      ftempsi2d(ji,jj) = ftempsi2d(ji,jj) + (ftempsi * fthk)
+                      ftempsi2d(ji,jj) = ftempsi2d(ji,jj) + (ftempsi(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%FASTFE%dgsave ) THEN
-                      ftempfe2d(ji,jj) =ftempfe2d(ji,jj)  + (ftempfe * fthk)  
+                      ftempfe2d(ji,jj) =ftempfe2d(ji,jj)  + (ftempfe(ji,jj) * fse3t(ji,jj,jk))  
                   ENDIF
                   IF( med_diag%FASTC%dgsave ) THEN
-                      ftempc2d(ji,jj) = ftempc2d(ji,jj) + (ftempc  * fthk)
+                      ftempc2d(ji,jj) = ftempc2d(ji,jj) + (ftempc(ji,jj)  * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%FASTCA%dgsave ) THEN
-                      ftempca2d(ji,jj) = ftempca2d(ji,jj) + (ftempca * fthk)
+                      ftempca2d(ji,jj) = ftempca2d(ji,jj) + (ftempca(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   !! 
                   IF( med_diag%REMINN%dgsave ) THEN
-                      freminn2d(ji,jj) = freminn2d(ji,jj) + (freminn  * fthk)
+                      freminn2d(ji,jj) = freminn2d(ji,jj) + (freminn(ji,jj)  * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%REMINSI%dgsave ) THEN
-                      freminsi2d(ji,jj) = freminsi2d(ji,jj) + (freminsi * fthk)
+                      freminsi2d(ji,jj) = freminsi2d(ji,jj) + (freminsi(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%REMINFE%dgsave ) THEN
-                      freminfe2d(ji,jj)= freminfe2d(ji,jj) + (freminfe * fthk) 
+                      freminfe2d(ji,jj)= freminfe2d(ji,jj) + (freminfe(ji,jj) * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF( med_diag%REMINC%dgsave ) THEN
-                      freminc2d(ji,jj) = freminc2d(ji,jj) + (freminc  * fthk) 
+                      freminc2d(ji,jj) = freminc2d(ji,jj) + (freminc(ji,jj)  * fse3t(ji,jj,jk)) 
                   ENDIF
                   IF( med_diag%REMINCA%dgsave ) THEN
-                      freminca2d(ji,jj) = freminca2d(ji,jj) + (freminca * fthk) 
+                      freminca2d(ji,jj) = freminca2d(ji,jj) + (freminca(ji,jj) * fse3t(ji,jj,jk)) 
                   ENDIF
                   !!
@@ -4200 +4174 @@
                         ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)
                      ENDIF
-                  ELSE IF (jk.eq.jmbathy) THEN
+                  ELSE IF (jk.eq.mbathy(ji,jj)) THEN
                      IF( med_diag%IBEN_N%dgsave ) THEN
                         iben_n2d(ji,jj) = f_sbenin_n(ji,jj)  + f_fbenin_n(ji,jj)
@@ -4244 +4218 @@
                   !!
                   IF( med_diag%RIV_N%dgsave ) THEN
-                     rivn2d(ji,jj) = rivn2d(ji,jj) +  (f_riv_loc_n * fthk)
+                     rivn2d(ji,jj) = rivn2d(ji,jj) +  (f_riv_loc_n(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%RIV_SI%dgsave ) THEN
-                     rivsi2d(ji,jj) = rivsi2d(ji,jj) +  (f_riv_loc_si * fthk)
+                     rivsi2d(ji,jj) = rivsi2d(ji,jj) +  (f_riv_loc_si(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%RIV_C%dgsave ) THEN
-                     rivc2d(ji,jj) = rivc2d(ji,jj) +  (f_riv_loc_c * fthk)
+                     rivc2d(ji,jj) = rivc2d(ji,jj) +  (f_riv_loc_c(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%RIV_ALK%dgsave ) THEN
-                     rivalk2d(ji,jj) = rivalk2d(ji,jj) +  (f_riv_loc_alk * fthk)
+                     rivalk2d(ji,jj) = rivalk2d(ji,jj) +  (f_riv_loc_alk(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%DETC%dgsave ) THEN
-                     fslowc2d(ji,jj) = fslowc2d(ji,jj) + (fslowc  * fthk)   
+                     fslowc2d(ji,jj) = fslowc2d(ji,jj) + (fslowc(ji,jj)  * fse3t(ji,jj,jk))   
                   ENDIF
                   !! 
@@ -4262 +4236 @@
                   !!
                   IF( med_diag%PN_LLOSS%dgsave ) THEN
-                     fdpn22d(ji,jj) = fdpn22d(ji,jj) + (fdpn2  * fthk)
+                     fdpn22d(ji,jj) = fdpn22d(ji,jj) + (fdpn2(ji,jj)  * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%PD_LLOSS%dgsave ) THEN
-                     fdpd22d(ji,jj) = fdpd22d(ji,jj) + (fdpd2  * fthk)
+                     fdpd22d(ji,jj) = fdpd22d(ji,jj) + (fdpd2(ji,jj)  * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%ZI_LLOSS%dgsave ) THEN
-                     fdzmi22d(ji,jj) = fdzmi22d(ji,jj) + (fdzmi2 * fthk)
+                     fdzmi22d(ji,jj) = fdzmi22d(ji,jj) + (fdzmi2(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%ZE_LLOSS%dgsave ) THEN
-                     fdzme22d(ji,jj) = fdzme22d(ji,jj) + (fdzme2 * fthk)
+                     fdzme22d(ji,jj) = fdzme22d(ji,jj) + (fdzme2(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%ZI_MES_N%dgsave ) THEN
                      zimesn2d(ji,jj) = zimesn2d(ji,jj) +  &
-                     (xphi * (fgmipn + fgmid) * fthk)
+                     (xphi * (fgmipn(ji,jj) + fgmid(ji,jj)) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%ZI_MES_D%dgsave ) THEN
                      zimesd2d(ji,jj) = zimesd2d(ji,jj) + & 
-                     ((1. - xbetan) * finmi * fthk)
+                     ((1. - xbetan) * finmi(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%ZI_MES_C%dgsave ) THEN
                      zimesc2d(ji,jj) = zimesc2d(ji,jj) + &
-                     (xphi * ((xthetapn * fgmipn) + fgmidc) * fthk)
+                     (xphi * ((xthetapn * fgmipn(ji,jj)) + fgmidc(ji,jj)) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%ZI_MESDC%dgsave ) THEN
                      zimesdc2d(ji,jj) = zimesdc2d(ji,jj) + &
-                     ((1. - xbetac) * ficmi * fthk)
+                     ((1. - xbetac) * ficmi(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%ZI_EXCR%dgsave ) THEN
-                     ziexcr2d(ji,jj) = ziexcr2d(ji,jj) +  (fmiexcr * fthk)
+                     ziexcr2d(ji,jj) = ziexcr2d(ji,jj) +  (fmiexcr(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%ZI_RESP%dgsave ) THEN
-                     ziresp2d(ji,jj) = ziresp2d(ji,jj) +  (fmiresp * fthk)
+                     ziresp2d(ji,jj) = ziresp2d(ji,jj) +  (fmiresp(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%ZI_GROW%dgsave ) THEN
-                     zigrow2d(ji,jj) = zigrow2d(ji,jj) + (fmigrow * fthk)
+                     zigrow2d(ji,jj) = zigrow2d(ji,jj) + (fmigrow(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%ZE_MES_N%dgsave ) THEN
                      zemesn2d(ji,jj) = zemesn2d(ji,jj) + &
-                     (xphi * (fgmepn + fgmepd + fgmezmi + fgmed) * fthk)
+                     (xphi * (fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj)) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%ZE_MES_D%dgsave ) THEN
                      zemesd2d(ji,jj) = zemesd2d(ji,jj) + &
-                     ((1. - xbetan) * finme * fthk)
+                     ((1. - xbetan) * finme(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%ZE_MES_C%dgsave ) THEN
                      zemesc2d(ji,jj) = zemesc2d(ji,jj) +                         & 
-                     (xphi * ((xthetapn * fgmepn) + (xthetapd * fgmepd) +  &
-                     (xthetazmi * fgmezmi) + fgmedc) * fthk)
+                     (xphi * ((xthetapn * fgmepn(ji,jj)) + (xthetapd * fgmepd(ji,jj)) +  &
+                     (xthetazmi * fgmezmi(ji,jj)) + fgmedc(ji,jj)) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%ZE_MESDC%dgsave ) THEN
                      zemesdc2d(ji,jj) = zemesdc2d(ji,jj) +  &
-                     ((1. - xbetac) * ficme * fthk)
+                     ((1. - xbetac) * ficme(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%ZE_EXCR%dgsave ) THEN
-                     zeexcr2d(ji,jj) = zeexcr2d(ji,jj) + (fmeexcr * fthk)
+                     zeexcr2d(ji,jj) = zeexcr2d(ji,jj) + (fmeexcr(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%ZE_RESP%dgsave ) THEN
-                     zeresp2d(ji,jj) = zeresp2d(ji,jj) + (fmeresp * fthk)
+                     zeresp2d(ji,jj) = zeresp2d(ji,jj) + (fmeresp(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%ZE_GROW%dgsave ) THEN
-                     zegrow2d(ji,jj) = zegrow2d(ji,jj) + (fmegrow * fthk)
+                     zegrow2d(ji,jj) = zegrow2d(ji,jj) + (fmegrow(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%MDETC%dgsave ) THEN
-                     mdetc2d(ji,jj) = mdetc2d(ji,jj) + (fddc * fthk)
+                     mdetc2d(ji,jj) = mdetc2d(ji,jj) + (fddc(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%GMIDC%dgsave ) THEN
-                     gmidc2d(ji,jj) = gmidc2d(ji,jj) + (fgmidc * fthk)
+                     gmidc2d(ji,jj) = gmidc2d(ji,jj) + (fgmidc(ji,jj) * fse3t(ji,jj,jk))
                   ENDIF
                   IF( med_diag%GMEDC%dgsave ) THEN
-                     gmedc2d(ji,jj) = gmedc2d(ji,jj) + (fgmedc  * fthk)
+                     gmedc2d(ji,jj) = gmedc2d(ji,jj) + (fgmedc(ji,jj)  * fse3t(ji,jj,jk))
                   ENDIF
                   !!
@@ -4338 +4312 @@
                   !! ** 3D diagnostics
                   IF( med_diag%TPP3%dgsave ) THEN
-                     tpp3d(ji,jj,jk) =  (fprn * zphn) + (fprd * zphd)
+                     tpp3d(ji,jj,jk) =  (fprn(ji,jj) * zphn(ji,jj)) + (fprd(ji,jj) * zphd(ji,jj))
                      !CALL iom_put( "TPP3"  , tpp3d )
                   ENDIF
                   IF( med_diag%TPPD3%dgsave ) THEN
-                     tppd3(ji,jj,jk) =  (fprd * zphd)
+                     tppd3(ji,jj,jk) =  (fprd(ji,jj) * zphd(ji,jj))
                   ENDIF
                   
                   IF( med_diag%REMIN3N%dgsave ) THEN
-                     remin3dn(ji,jj,jk) = fregen + (freminn * fthk) !! remineralisation
+                     remin3dn(ji,jj,jk) = fregen(ji,jj) + (freminn(ji,jj) * fse3t(ji,jj,jk)) !! remineralisation
                      !CALL iom_put( "REMIN3N"  , remin3dn )
                   ENDIF
@@ -4358 +4332 @@
         !! AXY (09/11/16): CMIP6 diagnostics
         IF ( med_diag%DCALC3%dgsave   ) THEN
-                     dcalc3(ji,jj,jk) = freminca
+                     dcalc3(ji,jj,jk) = freminca(ji,jj)
                   ENDIF
         IF ( med_diag%FEDISS3%dgsave  ) THEN
-                     fediss3(ji,jj,jk) = ffetop
+                     fediss3(ji,jj,jk) = ffetop(ji,jj)
                   ENDIF
         IF ( med_diag%FESCAV3%dgsave  ) THEN
-                     fescav3(ji,jj,jk) = ffescav
+                     fescav3(ji,jj,jk) = ffescav(ji,jj)
                   ENDIF
         IF ( med_diag%MIGRAZP3%dgsave ) THEN
-                     migrazp3(ji,jj,jk) = fgmipn * xthetapn
+                     migrazp3(ji,jj,jk) = fgmipn(ji,jj) * xthetapn
                   ENDIF
         IF ( med_diag%MIGRAZD3%dgsave ) THEN
-                     migrazd3(ji,jj,jk) = fgmidc
+                     migrazd3(ji,jj,jk) = fgmidc(ji,jj)
                   ENDIF
         IF ( med_diag%MEGRAZP3%dgsave ) THEN
-                     megrazp3(ji,jj,jk) = (fgmepn * xthetapn) + (fgmepd * xthetapd)
+                     megrazp3(ji,jj,jk) = (fgmepn(ji,jj) * xthetapn) + (fgmepd(ji,jj) * xthetapd)
                   ENDIF
         IF ( med_diag%MEGRAZD3%dgsave ) THEN
-                     megrazd3(ji,jj,jk) = fgmedc
+                     megrazd3(ji,jj,jk) = fgmedc(ji,jj)
                   ENDIF
         IF ( med_diag%MEGRAZZ3%dgsave ) THEN
-                     megrazz3(ji,jj,jk) = (fgmezmi * xthetazmi)
+                     megrazz3(ji,jj,jk) = (fgmezmi(ji,jj) * xthetazmi)
                   ENDIF
         IF ( med_diag%PBSI3%dgsave    ) THEN
-                     pbsi3(ji,jj,jk)    = (fprds * zpds)
+                     pbsi3(ji,jj,jk)    = (fprds(ji,jj) * zpds(ji,jj))
                   ENDIF
         IF ( med_diag%PCAL3%dgsave    ) THEN
-                     pcal3(ji,jj,jk)    = ftempca
+                     pcal3(ji,jj,jk)    = ftempca(ji,jj)
                   ENDIF
         IF ( med_diag%REMOC3%dgsave   ) THEN
-                     remoc3(ji,jj,jk)   = freminc
+                     remoc3(ji,jj,jk)   = freminc(ji,jj)
                   ENDIF
         IF ( med_diag%PNLIMJ3%dgsave  ) THEN
-                     ! pnlimj3(ji,jj,jk)  = fjln
-                     pnlimj3(ji,jj,jk)  = fjlim_pn
+                     ! pnlimj3(ji,jj,jk)  = fjln(ji,jj)
+                     pnlimj3(ji,jj,jk)  = fjlim_pn(ji,jj)
                   ENDIF
         IF ( med_diag%PNLIMN3%dgsave  ) THEN
-                     pnlimn3(ji,jj,jk)  = fnln
+                     pnlimn3(ji,jj,jk)  = fnln(ji,jj)
                   ENDIF
         IF ( med_diag%PNLIMFE3%dgsave ) THEN
-                     pnlimfe3(ji,jj,jk) = ffln
+                     pnlimfe3(ji,jj,jk) = ffln(ji,jj)
                   ENDIF
         IF ( med_diag%PDLIMJ3%dgsave  ) THEN
-                     ! pdlimj3(ji,jj,jk)  = fjld
-                     pdlimj3(ji,jj,jk)  = fjlim_pd
+                     ! pdlimj3(ji,jj,jk)  = fjld(ji,jj)
+                     pdlimj3(ji,jj,jk)  = fjlim_pd(ji,jj)
                   ENDIF
         IF ( med_diag%PDLIMN3%dgsave  ) THEN
-                     pdlimn3(ji,jj,jk)  = fnld
+                     pdlimn3(ji,jj,jk)  = fnld(ji,jj)
                   ENDIF
         IF ( med_diag%PDLIMFE3%dgsave ) THEN
-                     pdlimfe3(ji,jj,jk) = ffld
+                     pdlimfe3(ji,jj,jk) = ffld(ji,jj)
                   ENDIF
         IF ( med_diag%PDLIMSI3%dgsave ) THEN
-                     pdlimsi3(ji,jj,jk) = fsld2
+                     pdlimsi3(ji,jj,jk) = fsld2(ji,jj)
                   ENDIF
                   !!
@@ -4430 +4404 @@
                   trc2d(ji,jj,2)  =  ftot_si(ji,jj)                            !! silicon  inventory
                   trc2d(ji,jj,3)  =  ftot_fe(ji,jj)                            !! iron     inventory
-                  trc2d(ji,jj,4)  = trc2d(ji,jj,4)  + (fprn  * zphn * fthk)    !! non-diatom production
-                  trc2d(ji,jj,5)  = trc2d(ji,jj,5)  + (fdpn         * fthk)    !! non-diatom non-grazing losses
-                  trc2d(ji,jj,6)  = trc2d(ji,jj,6)  + (fprd  * zphd * fthk)    !! diatom production
-                  trc2d(ji,jj,7)  = trc2d(ji,jj,7)  + (fdpd         * fthk)    !! diatom non-grazing losses
+                  trc2d(ji,jj,4)  = trc2d(ji,jj,4)  + (fprn(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk))    !! non-diatom production
+                  trc2d(ji,jj,5)  = trc2d(ji,jj,5)  + (fdpn(ji,jj)         * fse3t(ji,jj,jk))    !! non-diatom non-grazing losses
+                  trc2d(ji,jj,6)  = trc2d(ji,jj,6)  + (fprd(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))    !! diatom production
+                  trc2d(ji,jj,7)  = trc2d(ji,jj,7)  + (fdpd(ji,jj)         * fse3t(ji,jj,jk))    !! diatom non-grazing losses
                   !! diagnostic field  8 is (ostensibly) supplied by trcsed.F            
-                  trc2d(ji,jj,9)  = trc2d(ji,jj,9)  + (fprds * zpds * fthk)    !! diatom silicon production
-                  trc2d(ji,jj,10) = trc2d(ji,jj,10) + (fsdiss  * fthk)         !! diatom silicon dissolution
-                  trc2d(ji,jj,11) = trc2d(ji,jj,11) + (fgmipn  * fthk)         !! microzoo grazing on non-diatoms
-                  trc2d(ji,jj,12) = trc2d(ji,jj,12) + (fgmid   * fthk)         !! microzoo grazing on detrital nitrogen
-                  trc2d(ji,jj,13) = trc2d(ji,jj,13) + (fdzmi   * fthk)         !! microzoo non-grazing losses
-                  trc2d(ji,jj,14) = trc2d(ji,jj,14) + (fgmepn  * fthk)         !! mesozoo  grazing on non-diatoms
-                  trc2d(ji,jj,15) = trc2d(ji,jj,15) + (fgmepd  * fthk)         !! mesozoo  grazing on diatoms
-                  trc2d(ji,jj,16) = trc2d(ji,jj,16) + (fgmezmi * fthk)         !! mesozoo  grazing on microzoo
-                  trc2d(ji,jj,17) = trc2d(ji,jj,17) + (fgmed   * fthk)         !! mesozoo  grazing on detrital nitrogen
-                  trc2d(ji,jj,18) = trc2d(ji,jj,18) + (fdzme   * fthk)         !! mesozoo  non-grazing losses
+                  trc2d(ji,jj,9)  = trc2d(ji,jj,9)  + (fprds(ji,jj) * zpds(ji,jj) * fse3t(ji,jj,jk))    !! diatom silicon production
+                  trc2d(ji,jj,10) = trc2d(ji,jj,10) + (fsdiss(ji,jj)  * fse3t(ji,jj,jk))         !! diatom silicon dissolution
+                  trc2d(ji,jj,11) = trc2d(ji,jj,11) + (fgmipn(ji,jj)  * fse3t(ji,jj,jk))         !! microzoo grazing on non-diatoms
+                  trc2d(ji,jj,12) = trc2d(ji,jj,12) + (fgmid(ji,jj)   * fse3t(ji,jj,jk))         !! microzoo grazing on detrital nitrogen
+                  trc2d(ji,jj,13) = trc2d(ji,jj,13) + (fdzmi(ji,jj)   * fse3t(ji,jj,jk))         !! microzoo non-grazing losses
+                  trc2d(ji,jj,14) = trc2d(ji,jj,14) + (fgmepn(ji,jj)  * fse3t(ji,jj,jk))         !! mesozoo  grazing on non-diatoms
+                  trc2d(ji,jj,15) = trc2d(ji,jj,15) + (fgmepd(ji,jj)  * fse3t(ji,jj,jk))         !! mesozoo  grazing on diatoms
+                  trc2d(ji,jj,16) = trc2d(ji,jj,16) + (fgmezmi(ji,jj) * fse3t(ji,jj,jk))         !! mesozoo  grazing on microzoo
+                  trc2d(ji,jj,17) = trc2d(ji,jj,17) + (fgmed(ji,jj)   * fse3t(ji,jj,jk))         !! mesozoo  grazing on detrital nitrogen
+                  trc2d(ji,jj,18) = trc2d(ji,jj,18) + (fdzme(ji,jj)   * fse3t(ji,jj,jk))         !! mesozoo  non-grazing losses
                   !! diagnostic field 19 is (ostensibly) supplied by trcexp.F
-                  trc2d(ji,jj,20) = trc2d(ji,jj,20) + (fslown  * fthk)         !! slow sinking detritus N production
-                  trc2d(ji,jj,21) = trc2d(ji,jj,21) + (fdd     * fthk)         !! detrital remineralisation
-                  trc2d(ji,jj,22) = trc2d(ji,jj,22) + (ffetop  * fthk)         !! aeolian  iron addition
-                  trc2d(ji,jj,23) = trc2d(ji,jj,23) + (ffebot  * fthk)         !! seafloor iron addition
-                  trc2d(ji,jj,24) = trc2d(ji,jj,24) + (ffescav * fthk)         !! "free" iron scavenging
-                  trc2d(ji,jj,25) = trc2d(ji,jj,25) + (fjlim_pn * zphn * fthk) !! non-diatom J  limitation term 
-                  trc2d(ji,jj,26) = trc2d(ji,jj,26) + (fnln  * zphn * fthk)    !! non-diatom N  limitation term 
-                  trc2d(ji,jj,27) = trc2d(ji,jj,27) + (ffln  * zphn * fthk)    !! non-diatom Fe limitation term 
-                  trc2d(ji,jj,28) = trc2d(ji,jj,28) + (fjlim_pd * zphd * fthk) !! diatom     J  limitation term 
-                  trc2d(ji,jj,29) = trc2d(ji,jj,29) + (fnld  * zphd * fthk)    !! diatom     N  limitation term 
-                  trc2d(ji,jj,30) = trc2d(ji,jj,30) + (ffld  * zphd * fthk)    !! diatom     Fe limitation term 
-                  trc2d(ji,jj,31) = trc2d(ji,jj,31) + (fsld2 * zphd * fthk)    !! diatom     Si limitation term 
-                  trc2d(ji,jj,32) = trc2d(ji,jj,32) + (fsld  * zphd * fthk)    !! diatom     Si uptake limitation term
+                  trc2d(ji,jj,20) = trc2d(ji,jj,20) + (fslown(ji,jj)  * fse3t(ji,jj,jk))         !! slow sinking detritus N production
+                  trc2d(ji,jj,21) = trc2d(ji,jj,21) + (fdd(ji,jj)     * fse3t(ji,jj,jk))         !! detrital remineralisation
+                  trc2d(ji,jj,22) = trc2d(ji,jj,22) + (ffetop(ji,jj)  * fse3t(ji,jj,jk))         !! aeolian  iron addition
+                  trc2d(ji,jj,23) = trc2d(ji,jj,23) + (ffebot(ji,jj)  * fse3t(ji,jj,jk))         !! seafloor iron addition
+                  trc2d(ji,jj,24) = trc2d(ji,jj,24) + (ffescav(ji,jj) * fse3t(ji,jj,jk))         !! "free" iron scavenging
+                  trc2d(ji,jj,25) = trc2d(ji,jj,25) + (fjlim_pn(ji,jj) * zphn(ji,jj) * fse3t(ji,jj,jk)) !! non-diatom J  limitation term 
+                  trc2d(ji,jj,26) = trc2d(ji,jj,26) + (fnln(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk))    !! non-diatom N  limitation term 
+                  trc2d(ji,jj,27) = trc2d(ji,jj,27) + (ffln(ji,jj)  * zphn(ji,jj) * fse3t(ji,jj,jk))    !! non-diatom Fe limitation term 
+                  trc2d(ji,jj,28) = trc2d(ji,jj,28) + (fjlim_pd(ji,jj) * zphd(ji,jj) * fse3t(ji,jj,jk)) !! diatom     J  limitation term 
+                  trc2d(ji,jj,29) = trc2d(ji,jj,29) + (fnld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))    !! diatom     N  limitation term 
+                  trc2d(ji,jj,30) = trc2d(ji,jj,30) + (ffld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))    !! diatom     Fe limitation term 
+                  trc2d(ji,jj,31) = trc2d(ji,jj,31) + (fsld2(ji,jj) * zphd(ji,jj) * fse3t(ji,jj,jk))    !! diatom     Si limitation term 
+                  trc2d(ji,jj,32) = trc2d(ji,jj,32) + (fsld(ji,jj)  * zphd(ji,jj) * fse3t(ji,jj,jk))    !! diatom     Si uptake limitation term
                   if (jk.eq.i0100) trc2d(ji,jj,33) = fslownflux(ji,jj)         !! slow detritus flux at  100 m
                   if (jk.eq.i0200) trc2d(ji,jj,34) = fslownflux(ji,jj)         !! slow detritus flux at  200 m
                   if (jk.eq.i0500) trc2d(ji,jj,35) = fslownflux(ji,jj)         !! slow detritus flux at  500 m
                   if (jk.eq.i1000) trc2d(ji,jj,36) = fslownflux(ji,jj)         !! slow detritus flux at 1000 m
-                  trc2d(ji,jj,37) = trc2d(ji,jj,37) + fregen                   !! non-fast N  full column regeneration
-                  trc2d(ji,jj,38) = trc2d(ji,jj,38) + fregensi                 !! non-fast Si full column regeneration
+                  trc2d(ji,jj,37) = trc2d(ji,jj,37) + fregen(ji,jj)                   !! non-fast N  full column regeneration
+                  trc2d(ji,jj,38) = trc2d(ji,jj,38) + fregensi(ji,jj)                 !! non-fast Si full column regeneration
                   if (jk.eq.i0100) trc2d(ji,jj,39) = trc2d(ji,jj,37)           !! non-fast N  regeneration to  100 m
                   if (jk.eq.i0200) trc2d(ji,jj,40) = trc2d(ji,jj,37)           !! non-fast N  regeneration to  200 m
                   if (jk.eq.i0500) trc2d(ji,jj,41) = trc2d(ji,jj,37)           !! non-fast N  regeneration to  500 m
                   if (jk.eq.i1000) trc2d(ji,jj,42) = trc2d(ji,jj,37)           !! non-fast N  regeneration to 1000 m
-                  trc2d(ji,jj,43) = trc2d(ji,jj,43) + (ftempn  * fthk)         !! fast sinking detritus N production
-                  trc2d(ji,jj,44) = trc2d(ji,jj,44) + (ftempsi * fthk)         !! fast sinking detritus Si production
-                  trc2d(ji,jj,45) = trc2d(ji,jj,45) + (ftempfe * fthk)         !! fast sinking detritus Fe production
-                  trc2d(ji,jj,46) = trc2d(ji,jj,46) + (ftempc  * fthk)         !! fast sinking detritus C production
-                  trc2d(ji,jj,47) = trc2d(ji,jj,47) + (ftempca * fthk)         !! fast sinking detritus CaCO3 production
+                  trc2d(ji,jj,43) = trc2d(ji,jj,43) + (ftempn(ji,jj)  * fse3t(ji,jj,jk))         !! fast sinking detritus N production
+                  trc2d(ji,jj,44) = trc2d(ji,jj,44) + (ftempsi(ji,jj) * fse3t(ji,jj,jk))         !! fast sinking detritus Si production
+                  trc2d(ji,jj,45) = trc2d(ji,jj,45) + (ftempfe(ji,jj) * fse3t(ji,jj,jk))         !! fast sinking detritus Fe production
+                  trc2d(ji,jj,46) = trc2d(ji,jj,46) + (ftempc(ji,jj)  * fse3t(ji,jj,jk))         !! fast sinking detritus C production
+                  trc2d(ji,jj,47) = trc2d(ji,jj,47) + (ftempca(ji,jj) * fse3t(ji,jj,jk))         !! fast sinking detritus CaCO3 production
                   if (jk.eq.i0100) trc2d(ji,jj,48) = ffastn(ji,jj)             !! fast detritus N  flux at  100 m
                   if (jk.eq.i0200) trc2d(ji,jj,49) = ffastn(ji,jj)             !! fast detritus N  flux at  200 m
@@ -4490 +4464 @@
                   if (jk.eq.i0500) trc2d(ji,jj,62) = fregenfastsi(ji,jj)       !! Si regeneration to  500 m
                   if (jk.eq.i1000) trc2d(ji,jj,63) = fregenfastsi(ji,jj)       !! Si regeneration to 1000 m
-                  trc2d(ji,jj,64) = trc2d(ji,jj,64) + (freminn  * fthk)        !! sum of fast-sinking N  fluxes
-                  trc2d(ji,jj,65) = trc2d(ji,jj,65) + (freminsi * fthk)        !! sum of fast-sinking Si fluxes
-                  trc2d(ji,jj,66) = trc2d(ji,jj,66) + (freminfe * fthk)        !! sum of fast-sinking Fe fluxes
-                  trc2d(ji,jj,67) = trc2d(ji,jj,67) + (freminc  * fthk)        !! sum of fast-sinking C  fluxes
-                  trc2d(ji,jj,68) = trc2d(ji,jj,68) + (freminca * fthk)        !! sum of fast-sinking Ca fluxes
-                  if (jk.eq.jmbathy) then
+                  trc2d(ji,jj,64) = trc2d(ji,jj,64) + (freminn(ji,jj)  * fse3t(ji,jj,jk))        !! sum of fast-sinking N  fluxes
+                  trc2d(ji,jj,65) = trc2d(ji,jj,65) + (freminsi(ji,jj) * fse3t(ji,jj,jk))        !! sum of fast-sinking Si fluxes
+                  trc2d(ji,jj,66) = trc2d(ji,jj,66) + (freminfe(ji,jj) * fse3t(ji,jj,jk))        !! sum of fast-sinking Fe fluxes
+                  trc2d(ji,jj,67) = trc2d(ji,jj,67) + (freminc(ji,jj)  * fse3t(ji,jj,jk))        !! sum of fast-sinking C  fluxes
+                  trc2d(ji,jj,68) = trc2d(ji,jj,68) + (freminca(ji,jj) * fse3t(ji,jj,jk))        !! sum of fast-sinking Ca fluxes
+                  if (jk.eq.mbathy(ji,jj)) then
                      trc2d(ji,jj,69) = fsedn(ji,jj)                                   !! N  sedimentation flux                                  
                      trc2d(ji,jj,70) = fsedsi(ji,jj)                                  !! Si sedimentation flux
@@ -4504 +4478 @@
                   if (jk.eq.1)  trc2d(ji,jj,74) = qsr(ji,jj)
                   if (jk.eq.1)  trc2d(ji,jj,75) = xpar(ji,jj,jk)
-                  !! if (jk.eq.1)  trc2d(ji,jj,75) = real(iters)
+                  !! if (jk.eq.1)  trc2d(ji,jj,75) = real(iters(ji,jj))
                   !! diagnostic fields 76 to 80 calculated below
                   trc2d(ji,jj,81) = trc2d(ji,jj,81) + fprn_ml(ji,jj)           !! mixed layer non-diatom production
@@ -4525 +4499 @@
                   !! ROAM provisionally has access to a further 20 2D diagnostics
                   if (jk .eq. 1) then
-                     trc2d(ji,jj,91)  = trc2d(ji,jj,91)  + f_wind              !! surface wind
-                     trc2d(ji,jj,92)  = trc2d(ji,jj,92)  + f_pco2atm           !! atmospheric pCO2
-                     trc2d(ji,jj,93)  = trc2d(ji,jj,93)  + f_ph                !! ocean pH
-                     trc2d(ji,jj,94)  = trc2d(ji,jj,94)  + f_pco2w             !! ocean pCO2
-                     trc2d(ji,jj,95)  = trc2d(ji,jj,95)  + f_h2co3             !! ocean H2CO3 conc.
-                     trc2d(ji,jj,96)  = trc2d(ji,jj,96)  + f_hco3              !! ocean HCO3 conc.
-                     trc2d(ji,jj,97)  = trc2d(ji,jj,97)  + f_co3               !! ocean CO3 conc.
-                     trc2d(ji,jj,98)  = trc2d(ji,jj,98)  + f_co2flux           !! air-sea CO2 flux
+                     trc2d(ji,jj,91)  = trc2d(ji,jj,91)  + wndm(ji,jj)              !! surface wind
+                     trc2d(ji,jj,92)  = trc2d(ji,jj,92)  + f_pco2atm(ji,jj)           !! atmospheric pCO2
+                     trc2d(ji,jj,93)  = trc2d(ji,jj,93)  + f_ph(ji,jj)                !! ocean pH
+                     trc2d(ji,jj,94)  = trc2d(ji,jj,94)  + f_pco2w(ji,jj)             !! ocean pCO2
+                     trc2d(ji,jj,95)  = trc2d(ji,jj,95)  + f_h2co3(ji,jj)             !! ocean H2CO3 conc.
+                     trc2d(ji,jj,96)  = trc2d(ji,jj,96)  + f_hco3(ji,jj)              !! ocean HCO3 conc.
+                     trc2d(ji,jj,97)  = trc2d(ji,jj,97)  + f_co3(ji,jj)               !! ocean CO3 conc.
+                     trc2d(ji,jj,98)  = trc2d(ji,jj,98)  + f_co2flux(ji,jj)           !! air-sea CO2 flux
                      trc2d(ji,jj,99)  = trc2d(ji,jj,99)  + f_omcal(ji,jj)      !! ocean omega calcite 
                      trc2d(ji,jj,100) = trc2d(ji,jj,100) + f_omarg(ji,jj)      !! ocean omega aragonite
-                     trc2d(ji,jj,101) = trc2d(ji,jj,101) + f_TDIC              !! ocean TDIC
-                     trc2d(ji,jj,102) = trc2d(ji,jj,102) + f_TALK              !! ocean TALK
-                     trc2d(ji,jj,103) = trc2d(ji,jj,103) + f_kw660             !! surface kw660
-                     trc2d(ji,jj,104) = trc2d(ji,jj,104) + f_pp0               !! surface pressure
-                     trc2d(ji,jj,105) = trc2d(ji,jj,105) + f_o2flux            !! air-sea O2 flux
-                     trc2d(ji,jj,106) = trc2d(ji,jj,106) + f_o2sat             !! ocean O2 saturation
+                     trc2d(ji,jj,101) = trc2d(ji,jj,101) + f_TDIC(ji,jj)              !! ocean TDIC
+                     trc2d(ji,jj,102) = trc2d(ji,jj,102) + f_TALK(ji,jj)              !! ocean TALK
+                     trc2d(ji,jj,103) = trc2d(ji,jj,103) + f_kw660(ji,jj)             !! surface kw660
+                     trc2d(ji,jj,104) = trc2d(ji,jj,104) + f_pp0(ji,jj)               !! surface pressure
+                     trc2d(ji,jj,105) = trc2d(ji,jj,105) + f_o2flux(ji,jj)            !! air-sea O2 flux
+                     trc2d(ji,jj,106) = trc2d(ji,jj,106) + f_o2sat(ji,jj)             !! ocean O2 saturation
                      trc2d(ji,jj,107) = f2_ccd_cal(ji,jj)                      !! depth calcite CCD
                      trc2d(ji,jj,108) = f2_ccd_arg(ji,jj)                      !! depth aragonite CCD
                   endif
-                  if (jk .eq. jmbathy) then
+                  if (jk .eq. mbathy(ji,jj)) then
                      trc2d(ji,jj,109) = f3_omcal(ji,jj,jk)                     !! seafloor omega calcite
                      trc2d(ji,jj,110) = f3_omarg(ji,jj,jk)                     !! seafloor omega aragonite
@@ -4553 +4527 @@
                   if (jk.eq.i1000) trc2d(ji,jj,120) = ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall)  !! rain ratio at 1000 m
                   !! AXY (18/01/12): benthic flux diagnostics
-                  if (jk.eq.jmbathy) then
+                  if (jk.eq.mbathy(ji,jj)) then
                      trc2d(ji,jj,121) = f_sbenin_n(ji,jj)  + f_fbenin_n(ji,jj)
                      trc2d(ji,jj,122) = f_sbenin_fe(ji,jj) + f_fbenin_fe(ji,jj)
@@ -4568 +4542 @@
                   trc2d(ji,jj,136) = f_runoff(ji,jj)
                   !! AXY (19/07/12): amended to allow for riverine nutrient addition below surface
-                  trc2d(ji,jj,137) = trc2d(ji,jj,137) + (f_riv_loc_n * fthk)
-                  trc2d(ji,jj,138) = trc2d(ji,jj,138) + (f_riv_loc_si * fthk)
-                  trc2d(ji,jj,139) = trc2d(ji,jj,139) + (f_riv_loc_c * fthk)
-                  trc2d(ji,jj,140) = trc2d(ji,jj,140) + (f_riv_loc_alk * fthk)
-                  trc2d(ji,jj,141) = trc2d(ji,jj,141) + (fslowc  * fthk)       !! slow sinking detritus C production
+                  trc2d(ji,jj,137) = trc2d(ji,jj,137) + (f_riv_loc_n(ji,jj) * fse3t(ji,jj,jk))
+                  trc2d(ji,jj,138) = trc2d(ji,jj,138) + (f_riv_loc_si(ji,jj) * fse3t(ji,jj,jk))
+                  trc2d(ji,jj,139) = trc2d(ji,jj,139) + (f_riv_loc_c(ji,jj) * fse3t(ji,jj,jk))
+                  trc2d(ji,jj,140) = trc2d(ji,jj,140) + (f_riv_loc_alk(ji,jj) * fse3t(ji,jj,jk))
+                  trc2d(ji,jj,141) = trc2d(ji,jj,141) + (fslowc(ji,jj)  * fse3t(ji,jj,jk))       !! slow sinking detritus C production
                   if (jk.eq.i0100) trc2d(ji,jj,142) = fslowcflux(ji,jj)        !! slow detritus flux at  100 m
                   if (jk.eq.i0200) trc2d(ji,jj,143) = fslowcflux(ji,jj)        !! slow detritus flux at  200 m
@@ -4580 +4554 @@
                   trc2d(ji,jj,147)  = trc2d(ji,jj,147)  + ftot_a(ji,jj)        !! alkalinity inventory
                   trc2d(ji,jj,148)  = trc2d(ji,jj,148)  + ftot_o2(ji,jj)       !! oxygen     inventory
-                  if (jk.eq.jmbathy) then
+                  if (jk.eq.mbathy(ji,jj)) then
                      trc2d(ji,jj,149) = f_benout_lyso_ca(ji,jj)
                   endif
-                  trc2d(ji,jj,150) = fcomm_resp(ji,jj) * fthk                  !! community respiration
+                  trc2d(ji,jj,150) = fcomm_resp(ji,jj) * fse3t(ji,jj,jk)                  !! community respiration
         !!
         !! AXY (14/02/14): a Valentines Day gift to BASIN - a shedload of new
@@ -4593 +4567 @@
                   !! linear losses
                   !! non-diatom
-                  trc2d(ji,jj,151) = trc2d(ji,jj,151) + (fdpn2  * fthk)
+                  trc2d(ji,jj,151) = trc2d(ji,jj,151) + (fdpn2(ji,jj)  * fse3t(ji,jj,jk))
                   !! diatom
-                  trc2d(ji,jj,152) = trc2d(ji,jj,152) + (fdpd2  * fthk)
+                  trc2d(ji,jj,152) = trc2d(ji,jj,152) + (fdpd2(ji,jj)  * fse3t(ji,jj,jk))
                   !! microzooplankton
-                  trc2d(ji,jj,153) = trc2d(ji,jj,153) + (fdzmi2 * fthk)
+                  trc2d(ji,jj,153) = trc2d(ji,jj,153) + (fdzmi2(ji,jj) * fse3t(ji,jj,jk))
                   !! mesozooplankton
-                  trc2d(ji,jj,154) = trc2d(ji,jj,154) + (fdzme2 * fthk)
+                  trc2d(ji,jj,154) = trc2d(ji,jj,154) + (fdzme2(ji,jj) * fse3t(ji,jj,jk))
                   !! ----------------------------------------------------------------------
                   !! microzooplankton grazing
                   !! microzooplankton messy -> N
-                  trc2d(ji,jj,155) = trc2d(ji,jj,155) + (xphi * (fgmipn + fgmid) * fthk)
+                  trc2d(ji,jj,155) = trc2d(ji,jj,155) + (xphi * (fgmipn(ji,jj) + fgmid(ji,jj)) * fse3t(ji,jj,jk))
                   !! microzooplankton messy -> D
-                  trc2d(ji,jj,156) = trc2d(ji,jj,156) + ((1. - xbetan) * finmi * fthk)
+                  trc2d(ji,jj,156) = trc2d(ji,jj,156) + ((1. - xbetan) * finmi(ji,jj) * fse3t(ji,jj,jk))
                   !! microzooplankton messy -> DIC
-                  trc2d(ji,jj,157) = trc2d(ji,jj,157) + (xphi * ((xthetapn * fgmipn) + fgmidc) * fthk)
+                  trc2d(ji,jj,157) = trc2d(ji,jj,157) + (xphi * ((xthetapn * fgmipn(ji,jj)) + fgmidc(ji,jj)) * fse3t(ji,jj,jk))
                   !! microzooplankton messy -> Dc
-                  trc2d(ji,jj,158) = trc2d(ji,jj,158) + ((1. - xbetac) * ficmi * fthk)
+                  trc2d(ji,jj,158) = trc2d(ji,jj,158) + ((1. - xbetac) * ficmi(ji,jj) * fse3t(ji,jj,jk))
                   !! microzooplankton excretion
-                  trc2d(ji,jj,159) = trc2d(ji,jj,159) + (fmiexcr * fthk)
+                  trc2d(ji,jj,159) = trc2d(ji,jj,159) + (fmiexcr(ji,jj) * fse3t(ji,jj,jk))
                   !! microzooplankton respiration
-                  trc2d(ji,jj,160) = trc2d(ji,jj,160) + (fmiresp * fthk)
+                  trc2d(ji,jj,160) = trc2d(ji,jj,160) + (fmiresp(ji,jj) * fse3t(ji,jj,jk))
                   !! microzooplankton growth
-                  trc2d(ji,jj,161) = trc2d(ji,jj,161) + (fmigrow * fthk)
+                  trc2d(ji,jj,161) = trc2d(ji,jj,161) + (fmigrow(ji,jj) * fse3t(ji,jj,jk))
                   !! ----------------------------------------------------------------------
                   !! mesozooplankton grazing
                   !! mesozooplankton messy -> N
-                  trc2d(ji,jj,162) = trc2d(ji,jj,162) + (xphi * (fgmepn + fgmepd + fgmezmi + fgmed) * fthk)
+                  trc2d(ji,jj,162) = trc2d(ji,jj,162) + (xphi * (fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj)) * fse3t(ji,jj,jk))
                   !! mesozooplankton messy -> D
-                  trc2d(ji,jj,163) = trc2d(ji,jj,163) + ((1. - xbetan) * finme * fthk)
+                  trc2d(ji,jj,163) = trc2d(ji,jj,163) + ((1. - xbetan) * finme(ji,jj) * fse3t(ji,jj,jk))
                   !! mesozooplankton messy -> DIC
-                  trc2d(ji,jj,164) = trc2d(ji,jj,164) + (xphi * ((xthetapn * fgmepn) + (xthetapd * fgmepd) + &
-                  &                  (xthetazmi * fgmezmi) + fgmedc) * fthk)
+                  trc2d(ji,jj,164) = trc2d(ji,jj,164) + (xphi * ((xthetapn * fgmepn(ji,jj)) + (xthetapd * fgmepd(ji,jj)) + &
+                  &                  (xthetazmi * fgmezmi(ji,jj)) + fgmedc(ji,jj)) * fse3t(ji,jj,jk))
                   !! mesozooplankton messy -> Dc
-                  trc2d(ji,jj,165) = trc2d(ji,jj,165) + ((1. - xbetac) * ficme * fthk)
+                  trc2d(ji,jj,165) = trc2d(ji,jj,165) + ((1. - xbetac) * ficme(ji,jj) * fse3t(ji,jj,jk))
                   !! mesozooplankton excretion
-                  trc2d(ji,jj,166) = trc2d(ji,jj,166) + (fmeexcr * fthk)
+                  trc2d(ji,jj,166) = trc2d(ji,jj,166) + (fmeexcr(ji,jj) * fse3t(ji,jj,jk))
                   !! mesozooplankton respiration
-                  trc2d(ji,jj,167) = trc2d(ji,jj,167) + (fmeresp * fthk)
+                  trc2d(ji,jj,167) = trc2d(ji,jj,167) + (fmeresp(ji,jj) * fse3t(ji,jj,jk))
                   !! mesozooplankton growth
-                  trc2d(ji,jj,168) = trc2d(ji,jj,168) + (fmegrow * fthk)
+                  trc2d(ji,jj,168) = trc2d(ji,jj,168) + (fmegrow(ji,jj) * fse3t(ji,jj,jk))
                   !! ----------------------------------------------------------------------
                   !! miscellaneous
-                  trc2d(ji,jj,169) = trc2d(ji,jj,169) + (fddc    * fthk) !! detrital C remineralisation
-                  trc2d(ji,jj,170) = trc2d(ji,jj,170) + (fgmidc  * fthk) !! microzoo grazing on detrital carbon
-                  trc2d(ji,jj,171) = trc2d(ji,jj,171) + (fgmedc  * fthk) !! mesozoo  grazing on detrital carbon
+                  trc2d(ji,jj,169) = trc2d(ji,jj,169) + (fddc(ji,jj)    * fse3t(ji,jj,jk)) !! detrital C remineralisation
+                  trc2d(ji,jj,170) = trc2d(ji,jj,170) + (fgmidc(ji,jj)  * fse3t(ji,jj,jk)) !! microzoo grazing on detrital carbon
+                  trc2d(ji,jj,171) = trc2d(ji,jj,171) + (fgmedc(ji,jj)  * fse3t(ji,jj,jk)) !! mesozoo  grazing on detrital carbon
                   !!
                   !! ----------------------------------------------------------------------
@@ -4646 +4620 @@
         !!                 below this section)
         !!
-        !! extract fields at surface
-       !! if (jk .eq. 1) then
-                 !!    trc2d(ji,jj,172) = zchn              !! Pn chlorophyll
-                 !!    trc2d(ji,jj,173) = zphn              !! Pn biomass
-                 !!    trc2d(ji,jj,174) = fjln              !! Pn J-term
-                 !!    trc2d(ji,jj,175) = (fprn * zphn)     !! Pn PP
-                 !!    trc2d(ji,jj,176) = zchd              !! Pd chlorophyll
-                 !!    trc2d(ji,jj,177) = zphd              !! Pd biomass
-                 !!    trc2d(ji,jj,178) = fjld              !! Pd J-term
-                 !!    trc2d(ji,jj,179) = xpar(ji,jj,jk)    !! Pd PP
-                 !!    trc2d(ji,jj,180) = loc_T             !! local temperature
-                 !! endif
-       !! !!
-       !! !! extract fields at 50m (actually 44-50m)
-       !! if (jk .eq. 18) then
-                 !!    trc2d(ji,jj,181) = zchn              !! Pn chlorophyll
-                 !!    trc2d(ji,jj,182) = zphn              !! Pn biomass
-                 !!    trc2d(ji,jj,183) = fjln              !! Pn J-term
-                 !!    trc2d(ji,jj,184) = (fprn * zphn)     !! Pn PP
-                 !!    trc2d(ji,jj,185) = zchd              !! Pd chlorophyll
-                 !!    trc2d(ji,jj,186) = zphd              !! Pd biomass
-                 !!    trc2d(ji,jj,187) = fjld              !! Pd J-term
-                 !!    trc2d(ji,jj,188) = xpar(ji,jj,jk)    !! Pd PP
-                 !!    trc2d(ji,jj,189) = loc_T             !! local temperature
-                 !! endif
-       !! !!
-       !! !! extract fields at 100m
-       !! if (jk .eq. i0100) then
-                 !!    trc2d(ji,jj,190) = zchn              !! Pn chlorophyll
-                 !!    trc2d(ji,jj,191) = zphn              !! Pn biomass
-                 !!    trc2d(ji,jj,192) = fjln              !! Pn J-term
-                 !!    trc2d(ji,jj,193) = (fprn * zphn)     !! Pn PP
-                 !!    trc2d(ji,jj,194) = zchd              !! Pd chlorophyll
-                 !!    trc2d(ji,jj,195) = zphd              !! Pd biomass
-                 !!    trc2d(ji,jj,196) = fjld              !! Pd J-term
-                 !!    trc2d(ji,jj,197) = xpar(ji,jj,jk)    !! Pd PP
-                 !!    trc2d(ji,jj,198) = loc_T             !! local temperature
-                 !! endif
-                 !!
                   !! extract relevant BASIN fields at 150m
                   if (jk .eq. i0150) then
@@ -4738 +4673 @@
                   !!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
                   if ((jk .eq. 1) .and.( jdms.eq.1)) then
-                     trc2d(ji,jj,221) = dms_surf          !! DMS surface concentration 
+                     trc2d(ji,jj,221) = dms_surf(ji,jj)          !! DMS surface concentration 
                      !! AXY (13/03/15): add in other DMS estimates
-                     trc2d(ji,jj,222) = dms_andr          !! DMS surface concentration 
-                     trc2d(ji,jj,223) = dms_simo          !! DMS surface concentration 
-                     trc2d(ji,jj,224) = dms_aran          !! DMS surface concentration 
-                     trc2d(ji,jj,225) = dms_hall          !! DMS surface concentration 
+                     trc2d(ji,jj,222) = dms_andr(ji,jj)          !! DMS surface concentration 
+                     trc2d(ji,jj,223) = dms_simo(ji,jj)          !! DMS surface concentration 
+                     trc2d(ji,jj,224) = dms_aran(ji,jj)          !! DMS surface concentration 
+                     trc2d(ji,jj,225) = dms_hall(ji,jj)          !! DMS surface concentration 
                   endif
 # endif
@@ -4756 +4691 @@
                   !!----------------------------------------------------------------------
                   !!
-                  trc3d(ji,jj,jk,1)  = ((fprn + fprd) * zphn)     !! primary production  
+                  trc3d(ji,jj,jk,1)  = ((fprn(ji,jj) + fprd(ji,jj)) * zphn(ji,jj))     !! primary production  
                   trc3d(ji,jj,jk,2)  = fslownflux(ji,jj) + ffastn(ji,jj) !! detrital flux
-                  trc3d(ji,jj,jk,3)  = fregen + (freminn * fthk)  !! remineralisation
+                  trc3d(ji,jj,jk,3)  = fregen(ji,jj) + (freminn(ji,jj) * fse3t(ji,jj,jk))  !! remineralisation
 # if defined key_roam
                   trc3d(ji,jj,jk,4)  = f3_pH(ji,jj,jk)            !! pH