Index: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 5711)
+++ branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 5726)
@@ -28,4 +28,8 @@
    USE iom
    USE daymod
+   !! AXY (05/11/13): need these for MEDUSA to input/output benthic reservoirs
+   USE sms_medusa
+   USE trcsms_medusa
+   !!
    IMPLICIT NONE
    PRIVATE
@@ -102,4 +106,6 @@
       !!----------------------------------------------------------------------
       INTEGER  ::  jn     
+      !! AXY (05/11/13): temporary variables
+      REAL(wp) ::    fq0,fq1,fq2
 
       !!----------------------------------------------------------------------
@@ -117,4 +123,75 @@
          CALL iom_get( numrtr, jpdom_autoglo, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) )
       END DO
+
+      !! AXY (09/06/14): the ARCHER version of MEDUSA-2 does not include an equivalent
+      !!                 call to MEDUSA-2 at this point; this suggests that the FCM
+      !!                 version of NEMO date significantly earlier than the current
+      !!                 version
+
+#if defined key_medusa
+      !! AXY (13/01/12): check if the restart contains sediment fields;
+      !!                 this is only relevant for simulations that include
+      !!                 biogeochemistry and are restarted from earlier runs
+      !!                 in which there was no sediment component
+      !!
+      IF( iom_varid( numrtr, 'B_SED_N', ldstop = .FALSE. ) > 0 ) THEN
+         !! YES; in which case read them
+         !!
+         IF(lwp) WRITE(numout,*) ' MEDUSA sediment fields present - reading in ...'
+         CALL iom_get( numrtr, jpdom_autoglo, 'B_SED_N',  zb_sed_n(:,:)  )
+         CALL iom_get( numrtr, jpdom_autoglo, 'N_SED_N',  zn_sed_n(:,:)  )
+         CALL iom_get( numrtr, jpdom_autoglo, 'B_SED_FE', zb_sed_fe(:,:) )
+         CALL iom_get( numrtr, jpdom_autoglo, 'N_SED_FE', zn_sed_fe(:,:) )
+         CALL iom_get( numrtr, jpdom_autoglo, 'B_SED_SI', zb_sed_si(:,:) )
+         CALL iom_get( numrtr, jpdom_autoglo, 'N_SED_SI', zn_sed_si(:,:) )
+         CALL iom_get( numrtr, jpdom_autoglo, 'B_SED_C',  zb_sed_c(:,:)  )
+         CALL iom_get( numrtr, jpdom_autoglo, 'N_SED_C',  zn_sed_c(:,:)  )
+         CALL iom_get( numrtr, jpdom_autoglo, 'B_SED_CA', zb_sed_ca(:,:) )
+         CALL iom_get( numrtr, jpdom_autoglo, 'N_SED_CA', zn_sed_ca(:,:) )
+      ELSE
+         !! NO; in which case set them to zero
+         !!
+         IF(lwp) WRITE(numout,*) ' MEDUSA sediment fields absent - setting to zero ...'
+         zb_sed_n(:,:)  = 0.0   !! organic N
+         zn_sed_n(:,:)  = 0.0
+         zb_sed_fe(:,:) = 0.0   !! organic Fe
+         zn_sed_fe(:,:) = 0.0
+         zb_sed_si(:,:) = 0.0   !! inorganic Si
+         zn_sed_si(:,:) = 0.0
+         zb_sed_c(:,:)  = 0.0   !! organic C
+         zn_sed_c(:,:)  = 0.0
+         zb_sed_ca(:,:) = 0.0   !! inorganic C
+         zn_sed_ca(:,:) = 0.0
+      ENDIF
+      !!
+      !! calculate stats on these fields
+      IF(lwp) WRITE(numout,*) ' MEDUSA sediment field stats (min, max, sum) ...'
+      fq0 = MINVAL(zn_sed_n(:,:))
+      fq1 = MAXVAL(zn_sed_n(:,:))
+      fq2 = SUM(zn_sed_n(:,:))
+      if (lwp) write (numout,'(a,3f15.5)') 'Sediment  N ', &
+         &        fq0, fq1, fq2
+      fq0 = MINVAL(zn_sed_fe(:,:))
+      fq1 = MAXVAL(zn_sed_fe(:,:))
+      fq2 = SUM(zn_sed_fe(:,:))
+      if (lwp) write (numout,'(a,3f15.5)') 'Sediment Fe ', &
+         &        fq0, fq1, fq2
+      fq0 = MINVAL(zn_sed_si(:,:))
+      fq1 = MAXVAL(zn_sed_si(:,:))
+      fq2 = SUM(zn_sed_si(:,:))
+      if (lwp) write (numout,'(a,3f15.5)') 'Sediment Si ', &
+         &        fq0, fq1, fq2
+      fq0 = MINVAL(zn_sed_c(:,:))
+      fq1 = MAXVAL(zn_sed_c(:,:))
+      fq2 = SUM(zn_sed_c(:,:))
+      if (lwp) write (numout,'(a,3f15.5)') 'Sediment  C ', &
+         &        fq0, fq1, fq2
+      fq0 = MINVAL(zn_sed_ca(:,:))
+      fq1 = MAXVAL(zn_sed_ca(:,:))
+      fq2 = SUM(zn_sed_ca(:,:))
+      if (lwp) write (numout,'(a,3f15.5)') 'Sediment Ca ', &
+         &        fq0, fq1, fq2
+#endif
+ 
       !
    END SUBROUTINE trc_rst_read
@@ -130,4 +207,6 @@
       INTEGER  :: jn
       REAL(wp) :: zarak0
+      !! AXY (05/11/13): temporary variables
+      REAL(wp) ::    fq0,fq1,fq2
       !!----------------------------------------------------------------------
       !
@@ -142,4 +221,58 @@
          CALL iom_rstput( kt, nitrst, numrtw, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) )
       END DO
+
+      !! AXY (09/06/14): the ARCHER version of MEDUSA-2 does not include an equivalent
+      !!                 call to MEDUSA-2 at this point; this suggests that the FCM
+      !!                 version of NEMO date significantly earlier than the current
+      !!                 version
+
+#if defined key_medusa
+      !! AXY (13/01/12): write out "before" and "now" state of seafloor
+      !!                 sediment pools into restart; this happens
+      !!                 whether or not the pools are to be used by
+      !!                 MEDUSA (which is controlled by a switch in the
+      !!                 namelist_top file)
+      !!
+      IF(lwp) WRITE(numout,*) ' MEDUSA sediment fields - writing out ...'
+      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_N',  zb_sed_n(:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_N',  zn_sed_n(:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_FE', zb_sed_fe(:,:) )
+      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_FE', zn_sed_fe(:,:) )
+      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_SI', zb_sed_si(:,:) )
+      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_SI', zn_sed_si(:,:) )
+      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_C',  zb_sed_c(:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_C',  zn_sed_c(:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'B_SED_CA', zb_sed_ca(:,:) )
+      CALL iom_rstput( kt, nitrst, numrtw, 'N_SED_CA', zn_sed_ca(:,:) )
+      !!
+      !! calculate stats on these fields
+      IF(lwp) WRITE(numout,*) ' MEDUSA sediment field stats (min, max, sum) ...'
+      fq0 = MINVAL(zn_sed_n(:,:))
+      fq1 = MAXVAL(zn_sed_n(:,:))
+      fq2 = SUM(zn_sed_n(:,:))
+      if (lwp) write (numout,'(a,3f15.5)') 'Sediment  N ', &
+         &        fq0, fq1, fq2
+      fq0 = MINVAL(zn_sed_fe(:,:))
+      fq1 = MAXVAL(zn_sed_fe(:,:))
+      fq2 = SUM(zn_sed_fe(:,:))
+      if (lwp) write (numout,'(a,3f15.5)') 'Sediment Fe ', &
+         &        fq0, fq1, fq2
+      fq0 = MINVAL(zn_sed_si(:,:))
+      fq1 = MAXVAL(zn_sed_si(:,:))
+      fq2 = SUM(zn_sed_si(:,:))
+      if (lwp) write (numout,'(a,3f15.5)') 'Sediment Si ', &
+         &        fq0, fq1, fq2
+      fq0 = MINVAL(zn_sed_c(:,:))
+      fq1 = MAXVAL(zn_sed_c(:,:))
+      fq2 = SUM(zn_sed_c(:,:))
+      if (lwp) write (numout,'(a,3f15.5)') 'Sediment  C ', &
+         &        fq0, fq1, fq2
+      fq0 = MINVAL(zn_sed_ca(:,:))
+      fq1 = MAXVAL(zn_sed_ca(:,:))
+      fq2 = SUM(zn_sed_ca(:,:))
+      if (lwp) write (numout,'(a,3f15.5)') 'Sediment Ca ', &
+         &        fq0, fq1, fq2
+#endif
+
       !
       IF( kt == nitrst ) THEN
