Index: /NEMO/branches/NERC/dev_r4.0.2_NERC_Externals/src/TOP/trcrst.F90
===================================================================
--- /NEMO/branches/NERC/dev_r4.0.2_NERC_Externals/src/TOP/trcrst.F90	(revision 13149)
+++ /NEMO/branches/NERC/dev_r4.0.2_NERC_Externals/src/TOP/trcrst.F90	(revision 13150)
@@ -24,4 +24,5 @@
    USE daymod
    USE lib_mpp
+   USE sms_medusa
    
    IMPLICIT NONE
@@ -32,4 +33,7 @@
    PUBLIC   trc_rst_wri       ! called by ???
    PUBLIC   trc_rst_cal
+   PUBLIC   trc_rst_stat
+   PUBLIC   trc_rst_conserve  ! Conservation Checks 
+
 
    !!----------------------------------------------------------------------
@@ -149,4 +153,5 @@
       IF( kt == nitrst ) THEN
           CALL trc_rst_stat            ! statistics
+          CALL trc_rst_conserve        ! Conservation Checks 
           CALL iom_close( numrtw )     ! close the restart file (only at last time step)
 #if ! defined key_trdmxl_trc
@@ -338,4 +343,148 @@
    END SUBROUTINE trc_rst_stat
 
+   SUBROUTINE trc_rst_conserve
+      !!----------------------------------------------------------------------
+      !!                    ***  trc_rst_conserve  ***
+      !!
+      !! ** purpose  :   Compute tracers conservation statistics
+      !!
+      !! AXY (17/11/2017)
+      !! This routine calculates the "now" inventories of the elemental 
+      !! cycles of MEDUSA and compares them to those calculate when the
+      !! model was initialised / restarted; the cycles calculated are:
+      !!    nitrogen, silicon, iron, carbon, alkalinity and oxygen
+      !!----------------------------------------------------------------------
+      INTEGER  :: ji, jj, jk, jn
+      REAL(wp) :: zsum3d, zsum2d, zinvt, zdelta, zratio
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d, zvol
+      REAL(wp), DIMENSION(jpi,jpj)     :: z2d, zarea
+      REAL(wp), DIMENSION(6)           :: loc_cycletot3, loc_cycletot2
+      !!----------------------------------------------------------------------
+      !
+      IF( lwp ) THEN
+         WRITE(numout,*) 
+         WRITE(numout,*) '           ----TRACER CONSERVATION----             '
+         WRITE(numout,*) 
+      ENDIF
+      !
+      ! ocean volume
+      DO jk = 1, jpk
+         zvol(:,:,jk) = e1e2t(:,:) * e3t_a(:,:,jk) * tmask(:,:,jk)
+      END DO
+      !
+      ! ocean area (for sediments)
+      zarea(:,:)      = e1e2t(:,:) * tmask(:,:,1)
+      !
+      !----------------------------------------------------------------------
+      ! nitrogen
+      z3d(:,:,:) = trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + &
+                   trn(:,:,:,jpzme) + trn(:,:,:,jpdet) + trn(:,:,:,jpdin)
+      z2d(:,:)   = zn_sed_n(:,:)
+      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
+      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
+      ! total tracer, and delta
+      zinvt      = zsum3d + zsum2d
+      zdelta     = zinvt - cycletot(1)
+      zratio     = 1.0e2 * zdelta / cycletot(1)
+      !
+      IF ( lwp ) WRITE(numout,9010) 'nitrogen', zsum3d, zsum2d, zinvt,   &
+         cycletot(1), zdelta, zratio
+      IF ( lwp ) WRITE(numout,*) 
+      !
+      !----------------------------------------------------------------------
+      ! silicon
+      z3d(:,:,:) = trn(:,:,:,jppds) + trn(:,:,:,jpsil)
+      z2d(:,:)   = zn_sed_si(:,:)
+      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
+      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
+      ! total tracer, and delta
+      zinvt      = zsum3d + zsum2d
+      zdelta     = zinvt - cycletot(2)
+      zratio     = 1.0e2 * zdelta / cycletot(2)
+      !
+      IF ( lwp ) WRITE(numout,9010) 'silicon', zsum3d, zsum2d, zinvt,    &
+         cycletot(2), zdelta, zratio
+      IF ( lwp ) WRITE(numout,*) 
+      !
+      !----------------------------------------------------------------------
+      ! iron
+      z3d(:,:,:) = ((trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + &
+            trn(:,:,:,jpzme) + trn(:,:,:,jpdet)) * xrfn) + trn(:,:,:,jpfer)
+      z2d(:,:)   = zn_sed_fe(:,:)
+      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
+      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
+      ! total tracer, and delta
+      zinvt      = zsum3d + zsum2d
+      zdelta     = zinvt - cycletot(3)
+      zratio     = 1.0e2 * zdelta / cycletot(3)
+      !
+      IF ( lwp ) WRITE(numout,9010) 'iron', zsum3d, zsum2d, zinvt,       &
+         cycletot(3), zdelta, zratio
+      IF ( lwp ) WRITE(numout,*) 
+      !
+      !----------------------------------------------------------------------
+      ! carbon
+      z3d(:,:,:) = (trn(:,:,:,jpphn) * xthetapn)  + (trn(:,:,:,jpphd) * xthetapd)  + &
+                   (trn(:,:,:,jpzmi) * xthetazmi) + (trn(:,:,:,jpzme) * xthetazme) + &
+                   trn(:,:,:,jpdtc) + trn(:,:,:,jpdic)
+      z2d(:,:)   = zn_sed_c(:,:) + zn_sed_ca(:,:)
+      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
+      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
+      ! total tracer, and delta
+      zinvt      = zsum3d + zsum2d
+      zdelta     = zinvt - cycletot(4)
+      zratio     = 1.0e2 * zdelta / cycletot(4)
+      !
+      IF ( lwp ) WRITE(numout,9010) 'carbon', zsum3d, zsum2d, zinvt,     &
+         cycletot(4), zdelta, zratio
+      IF ( lwp ) WRITE(numout,*) 
+      !
+      !----------------------------------------------------------------------
+      ! alkalinity
+      z3d(:,:,:) = trn(:,:,:,jpalk)
+      z2d(:,:)   = zn_sed_ca(:,:) * 2.0
+      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
+      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
+      ! total tracer, and delta
+      zinvt      = zsum3d + zsum2d
+      zdelta     = zinvt - cycletot(5)
+      zratio     = 1.0e2 * zdelta / cycletot(5)
+      !
+      IF ( lwp ) WRITE(numout,9010) 'alkalinity', zsum3d, zsum2d, zinvt, &
+         cycletot(5), zdelta, zratio
+      IF ( lwp ) WRITE(numout,*) 
+      !
+      !----------------------------------------------------------------------
+      ! oxygen
+      z3d(:,:,:) = trn(:,:,:,jpoxy)
+      z2d(:,:)   = 0.0
+      zsum3d     = glob_sum( z3d(:,:,:) * zvol(:,:,:) )
+      zsum2d     = glob_sum( z2d(:,:) * zarea(:,:) )
+      ! total tracer, and delta
+      zinvt      = zsum3d + zsum2d
+      zdelta     = zinvt - cycletot(6)
+      zratio     = 1.0e2 * zdelta / cycletot(6)
+      !
+      IF ( lwp ) WRITE(numout,9010) 'oxygen', zsum3d, zsum2d, zinvt,     &
+         cycletot(6), zdelta, zratio
+      !
+      !----------------------------------------------------------------------
+      ! Check 
+      zsum3d        = glob_sum( zvol(:,:,:) )
+      zsum2d        = glob_sum( zarea(:,:) )
+      IF ( lwp ) THEN 
+         WRITE(numout,*)
+         WRITE(numout,*) ' check : cvol    : ', zsum3d
+         WRITE(numout,*) ' check : carea   : ', zsum2d
+         WRITE(numout,*)
+      ENDIF
+      !
+9010  FORMAT(' element:',a10,                     &
+             ' 3d sum:',e18.10,' 2d sum:',e18.10, &
+             ' total:',e18.10,' initial:',e18.10, &
+             ' delta:',e18.10,' %:',e18.10)
+      !
+   END SUBROUTINE trc_rst_conserve 
+
 #else
    !!----------------------------------------------------------------------
