Index: /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_iomput.F90
===================================================================
--- /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_iomput.F90	(revision 8061)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_iomput.F90	(revision 8062)
@@ -367,49 +367,51 @@
             ENDDO
          ENDDO
-      ELSE IF (jk.eq.mbathy(ji,jj)) THEN
+      ELSE
          DO jj = 2,jpjm1
             DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
-                  IF( med_diag%IBEN_N%dgsave ) THEN
-                     iben_n2d(ji,jj) = f_sbenin_n(ji,jj) +                   &
-                                       f_fbenin_n(ji,jj)
-                  ENDIF
-                  IF( med_diag%IBEN_FE%dgsave ) THEN
-                     iben_fe2d(ji,jj) = f_sbenin_fe(ji,jj) +                 &
-                                        f_fbenin_fe(ji,jj)
-                  ENDIF
-                  IF( med_diag%IBEN_C%dgsave ) THEN
-                     iben_c2d(ji,jj) = f_sbenin_c(ji,jj) +                   &
-                                       f_fbenin_c(ji,jj)
-                  ENDIF
-                  IF( med_diag%IBEN_SI%dgsave ) THEN
-                     iben_si2d(ji,jj) = f_fbenin_si(ji,jj)
-                  ENDIF
-                  IF( med_diag%IBEN_CA%dgsave ) THEN
-                     iben_ca2d(ji,jj) = f_fbenin_ca(ji,jj)
-                  ENDIF
-                  IF( med_diag%OBEN_N%dgsave ) THEN
-                     oben_n2d(ji,jj) = f_benout_n(ji,jj)
-                  ENDIF
-                  IF( med_diag%OBEN_FE%dgsave ) THEN
-                     oben_fe2d(ji,jj) = f_benout_fe(ji,jj)
-                  ENDIF
-                  IF( med_diag%OBEN_C%dgsave ) THEN
-                     oben_c2d(ji,jj) = f_benout_c(ji,jj)
-                  ENDIF
-                  IF( med_diag%OBEN_SI%dgsave ) THEN
-                     oben_si2d(ji,jj) = f_benout_si(ji,jj)
-                  ENDIF
-                  IF( med_diag%OBEN_CA%dgsave ) THEN
-                     oben_ca2d(ji,jj) = f_benout_ca(ji,jj)
-                  ENDIF
-                  IF( med_diag%SFR_OCAL%dgsave ) THEN
-                     sfr_ocal2d(ji,jj) = f3_omcal(ji,jj,jk)
-                  ENDIF
-                  IF( med_diag%SFR_OARG%dgsave ) THEN
-                     sfr_oarg2d(ji,jj) =  f3_omarg(ji,jj,jk)
-                  ENDIF
-                  IF( med_diag%LYSO_CA%dgsave ) THEN
-                     lyso_ca2d(ji,jj) = f_benout_lyso_ca(ji,jj)
+               IF (jk.eq.mbathy(ji,jj)) THEN
+                  IF (tmask(ji,jj,jk) == 1) THEN
+                     IF( med_diag%IBEN_N%dgsave ) THEN
+                        iben_n2d(ji,jj) = f_sbenin_n(ji,jj) +                &
+                                          f_fbenin_n(ji,jj)
+                     ENDIF
+                     IF( med_diag%IBEN_FE%dgsave ) THEN
+                        iben_fe2d(ji,jj) = f_sbenin_fe(ji,jj) +              &
+                                           f_fbenin_fe(ji,jj)
+                     ENDIF
+                     IF( med_diag%IBEN_C%dgsave ) THEN
+                        iben_c2d(ji,jj) = f_sbenin_c(ji,jj) +                &
+                                          f_fbenin_c(ji,jj)
+                     ENDIF
+                     IF( med_diag%IBEN_SI%dgsave ) THEN
+                        iben_si2d(ji,jj) = f_fbenin_si(ji,jj)
+                     ENDIF
+                     IF( med_diag%IBEN_CA%dgsave ) THEN
+                        iben_ca2d(ji,jj) = f_fbenin_ca(ji,jj)
+                     ENDIF
+                     IF( med_diag%OBEN_N%dgsave ) THEN
+                        oben_n2d(ji,jj) = f_benout_n(ji,jj)
+                     ENDIF
+                     IF( med_diag%OBEN_FE%dgsave ) THEN
+                        oben_fe2d(ji,jj) = f_benout_fe(ji,jj)
+                     ENDIF
+                     IF( med_diag%OBEN_C%dgsave ) THEN
+                        oben_c2d(ji,jj) = f_benout_c(ji,jj)
+                     ENDIF
+                     IF( med_diag%OBEN_SI%dgsave ) THEN
+                        oben_si2d(ji,jj) = f_benout_si(ji,jj)
+                     ENDIF
+                     IF( med_diag%OBEN_CA%dgsave ) THEN
+                        oben_ca2d(ji,jj) = f_benout_ca(ji,jj)
+                     ENDIF
+                     IF( med_diag%SFR_OCAL%dgsave ) THEN
+                        sfr_ocal2d(ji,jj) = f3_omcal(ji,jj,jk)
+                     ENDIF
+                     IF( med_diag%SFR_OARG%dgsave ) THEN
+                        sfr_oarg2d(ji,jj) =  f3_omarg(ji,jj,jk)
+                     ENDIF
+                     IF( med_diag%LYSO_CA%dgsave ) THEN
+                        lyso_ca2d(ji,jj) = f_benout_lyso_ca(ji,jj)
+                     ENDIF
                   ENDIF
                ENDIF
Index: /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_trc.F90
===================================================================
--- /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_trc.F90	(revision 8061)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_trc.F90	(revision 8062)
@@ -529,8 +529,8 @@
 
 
-      if (jk.eq.mbathy(ji,jj)) then
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               if (jk.eq.mbathy(ji,jj)) then
                   !! N  sedimentation flux
                   trc2d(ji,jj,69) = fsedn(ji,jj)
@@ -543,8 +543,8 @@
                   !! Ca sedimentation flux
                   trc2d(ji,jj,73) = fsedca(ji,jj)
-               ENDIF
-            ENDDO
-         ENDDO
-      endif
+               endif
+            ENDIF
+         ENDDO
+      ENDDO
 
       if (jk.eq.1) then
@@ -729,16 +729,16 @@
       endif
 
-      if (jk .eq. mbathy(ji,jj)) then
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               if (jk .eq. mbathy(ji,jj)) then
                   !! seafloor omega calcite
                   trc2d(ji,jj,109) = f3_omcal(ji,jj,jk)
                   !! seafloor omega aragonite
                   trc2d(ji,jj,110) = f3_omarg(ji,jj,jk)
-               ENDIF
-            ENDDO
-         ENDDO
-      endif
+               endif
+            ENDIF
+         ENDDO
+      ENDDO
 
       if (jk.eq.i0100) then
@@ -779,8 +779,8 @@
       endif
 
-      if (jk.eq.mbathy(ji,jj)) then
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            IF (tmask(ji,jj,jk) == 1) THEN
+               if (jk.eq.mbathy(ji,jj)) then
                   !! AXY (18/01/12): benthic flux diagnostics
                   trc2d(ji,jj,121) = f_sbenin_n(ji,jj)  + f_fbenin_n(ji,jj)
@@ -794,8 +794,8 @@
                   trc2d(ji,jj,129) = f_benout_si(ji,jj)
                   trc2d(ji,jj,130) = f_benout_ca(ji,jj)
-               ENDIF
-            ENDDO
-         ENDDO
-      endif
+               endif
+            ENDIF
+         ENDDO
+      ENDDO
 
       DO jj = 2,jpjm1
@@ -876,21 +876,9 @@
                !! oxygen     inventory
                trc2d(ji,jj,148)  = trc2d(ji,jj,148)  + ftot_o2(ji,jj)
-            ENDIF
-         ENDDO
-      ENDDO
-
-      if (jk.eq.mbathy(ji,jj)) then
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               IF (tmask(ji,jj,jk) == 1) THEN
+
+               if (jk.eq.mbathy(ji,jj)) then
                   trc2d(ji,jj,149) = f_benout_lyso_ca(ji,jj)
-               ENDIF
-            ENDDO
-         ENDDO
-      endif
-
-      DO jj = 2,jpjm1
-         DO ji = 2,jpim1
-            IF (tmask(ji,jj,jk) == 1) THEN
+               endif
+
                !! community respiration
                trc2d(ji,jj,150) = fcomm_resp(ji,jj) * fse3t(ji,jj,jk)
Index: /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag.F90
===================================================================
--- /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag.F90	(revision 8061)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag.F90	(revision 8062)
@@ -78,5 +78,5 @@
          DO jj = 2,jpjm1
             DO ji = 2,jpim1
-               if (tmask(ji,jj,1) == 1) then
+               if (tmask(ji,jj,jk) == 1) then
                   !! sum tracers for inventory checks
                   IF ( med_diag%INVTN%dgsave )   THEN
@@ -105,5 +105,5 @@
          DO jj = 2,jpjm1
             DO ji = 2,jpim1
-               if (tmask(ji,jj,1) == 1) then
+               if (tmask(ji,jj,jk) == 1) then
                   IF ( med_diag%INVTC%dgsave )  THEN
                      ftot_c(ji,jj)  = ftot_c(ji,jj) +                        & 
@@ -128,5 +128,5 @@
          DO jj = 2,jpjm1
             DO ji = 2,jpim1
-               if (tmask(ji,jj,1) == 1) then
+               if (tmask(ji,jj,jk) == 1) then
                   IF ( med_diag%INVTC%dgsave )  THEN
                      !!
@@ -155,5 +155,5 @@
          DO jj = 2,jpjm1
             DO ji = 2,jpim1
-               if (tmask(ji,jj,1) == 1) then
+               if (tmask(ji,jj,jk) == 1) then
                   IF ( med_diag%O2min%dgsave ) THEN
                      if ( zoxy(ji,jj) < o2min(ji,jj) ) then
Index: /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus_fast_sink.F90
===================================================================
--- /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus_fast_sink.F90	(revision 8061)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus_fast_sink.F90	(revision 8062)
@@ -116,5 +116,5 @@
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
-            IF (tmask(ji,jj,1) == 1) THEN
+            IF (tmask(ji,jj,jk) == 1) THEN
                !!
                !! this variable integrates the creation of slow sinking 
@@ -167,5 +167,5 @@
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
-            IF (tmask(ji,jj,1) == 1) THEN
+            IF (tmask(ji,jj,jk) == 1) THEN
                !!
                !! nitrogen
Index: /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/iron_chem_scav.F90
===================================================================
--- /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/iron_chem_scav.F90	(revision 8061)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/iron_chem_scav.F90	(revision 8062)
@@ -87,5 +87,5 @@
          DO ji = 2,jpim1
             !! OPEN wet point IF..THEN loop
-            if (tmask(ji,jj,1) == 1) then
+            if (tmask(ji,jj,jk) == 1) then
                !!
                !! total iron concentration (mmol Fe / m3 -> umol Fe / m3)
Index: /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90
===================================================================
--- /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90	(revision 8061)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90	(revision 8062)
@@ -121,5 +121,5 @@
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
-            if (tmask(ji,jj,1) == 1) then
+            if (tmask(ji,jj,jk) == 1) then
                !!----------------------------------------------------------
                !! Phytoplankton light limitation
@@ -197,5 +197,5 @@
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
-            if (tmask(ji,jj,1) == 1) then
+            if (tmask(ji,jj,jk) == 1) then
                !!----------------------------------------------------------
                !! Phytoplankton nutrient limitation
@@ -233,5 +233,5 @@
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
-            if (tmask(ji,jj,1) == 1) then
+            if (tmask(ji,jj,jk) == 1) then
                !!----------------------------------------------------------
                !! Primary production (non-diatoms)
@@ -254,5 +254,5 @@
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
-            if (tmask(ji,jj,1) == 1) then
+            if (tmask(ji,jj,jk) == 1) then
                !!----------------------------------------------------------
                !! Primary production (diatoms)
@@ -348,5 +348,5 @@
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
-            if (tmask(ji,jj,1) == 1) then
+            if (tmask(ji,jj,jk) == 1) then
                !!----------------------------------------------------------
                !! Mixed layer primary production
@@ -379,5 +379,5 @@
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
-            if (tmask(ji,jj,1) == 1) then
+            if (tmask(ji,jj,jk) == 1) then
                !!----------------------------------------------------------
                !! Vertical Integral --
@@ -407,5 +407,5 @@
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
-            if (tmask(ji,jj,1) == 1) then
+            if (tmask(ji,jj,jk) == 1) then
                !!----------------------------------------------------------
                !! More chlorophyll calculations
Index: /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90
===================================================================
--- /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90	(revision 8061)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90	(revision 8062)
@@ -84,5 +84,5 @@
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
-            if (tmask(ji,jj,1) == 1) then
+            if (tmask(ji,jj,jk) == 1) then
                !!----------------------------------------------------------
                !! Plankton mortality losses
@@ -112,5 +112,5 @@
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
-            if (tmask(ji,jj,1) == 1) then
+            if (tmask(ji,jj,jk) == 1) then
                !!
                !! diatom phytoplankton
@@ -135,5 +135,5 @@
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
-            if (tmask(ji,jj,1) == 1) then
+            if (tmask(ji,jj,jk) == 1) then
                !!
                !! microzooplankton
@@ -157,5 +157,5 @@
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
-            if (tmask(ji,jj,1) == 1) then
+            if (tmask(ji,jj,jk) == 1) then
                !!
                !! mesozooplankton
Index: /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
===================================================================
--- /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 8061)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 8062)
@@ -432,5 +432,5 @@
          DO jj = 2,jpjm1
             DO ji = 2,jpim1
-               if (tmask(ji,jj,1) == 1) then
+               if (tmask(ji,jj,jk) == 1) then
                   !! microzooplankton
                   zzmi(ji,jj) = max(0.,trn(ji,jj,jk,jpzmi))
@@ -452,5 +452,5 @@
          DO jj = 2,jpjm1
             DO ji = 2,jpim1
-               if (tmask(ji,jj,1) == 1) then
+               if (tmask(ji,jj,jk) == 1) then
                   !! detrital carbon
                   zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc))
@@ -495,5 +495,5 @@
          DO jj = 2,jpjm1
             DO ji = 2,jpim1
-               if (tmask(ji,jj,1) == 1) then
+               if (tmask(ji,jj,jk) == 1) then
                   !! implicit detrital carbon
                   zdtc(ji,jj) = zdet(ji,jj) * xthetad
@@ -505,5 +505,5 @@
          DO jj = 2,jpjm1
             DO ji = 2,jpim1
-               if (tmask(ji,jj,1) == 1) then
+               if (tmask(ji,jj,jk) == 1) then
                   if (idf.eq.1) then
                      !! AXY (15/01/10)
@@ -560,5 +560,5 @@
 !            DO jj = 2,jpjm1
 !               DO ji = 2,jpim1
-!                  if (tmask(ji,jj,1) == 1) then
+!                  if (tmask(ji,jj,jk) == 1) then
 !                     !! report state variable values
 !                     IF (lwp) write (numout,*)                               &
@@ -592,5 +592,5 @@
 !            DO jj = 2,jpjm1
 !               DO ji = 2,jpim1
-!                  if (tmask(ji,jj,1) == 1) then
+!                  if (tmask(ji,jj,jk) == 1) then
 !                     IF (lwp) write (numout,*)                               &
 !                       '------------------------------'
Index: /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/zooplankton.F90
===================================================================
--- /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/zooplankton.F90	(revision 8061)
+++ /branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/zooplankton.F90	(revision 8062)
@@ -122,5 +122,5 @@
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
-            if (tmask(ji,jj,1) == 1) then
+            if (tmask(ji,jj,jk) == 1) then
                !!
                !! which translates to these incoming N and C fluxes
@@ -175,5 +175,5 @@
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
-            if (tmask(ji,jj,1) == 1) then
+            if (tmask(ji,jj,jk) == 1) then
                !!----------------------------------------------------------
                !! Mesozooplankton second
@@ -230,5 +230,5 @@
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
-            if (tmask(ji,jj,1) == 1) then
+            if (tmask(ji,jj,jk) == 1) then
                !!
                !! the ideal food C:N ratio for mesozooplankton
@@ -281,5 +281,5 @@
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
-            if (tmask(ji,jj,1) == 1) then
+            if (tmask(ji,jj,jk) == 1) then
                fzmi_i(ji,jj)  = fzmi_i(ji,jj)  + fse3t(ji,jj,jk) *          &
                                 ( fgmipn(ji,jj) + fgmid(ji,jj) )
