Index: /branches/NERC/dev_r5518_GO6_conserv_check_up/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- /branches/NERC/dev_r5518_GO6_conserv_check_up/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 9207)
+++ /branches/NERC/dev_r5518_GO6_conserv_check_up/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 9208)
@@ -851,4 +851,11 @@
       CALL wrk_alloc( jpi,jpj, rgt33 )
       !
+#if defined key_coare_bulk_patch
+     rgt33 = MIN(zw10, 22.)
+     !! cd_neutral_10m = 1.e-3 * ( 2.7/rgt33  + 0.142 + 0.06*rgt33  +
+     !!                  0.0025*rgt33 **2 - 1.25e-9*rgt33 **6)
+     !! A new version May 31st!
+      cd_neutral_10m = 1.e-3 * ( 2.4/rgt33  + 0.15 + 0.095*rgt33  + 0.0008*rgt33 **2 - 1.0e-9*rgt33 **6)
+#else
       !! When wind speed > 33 m/s => Cyclone conditions => special treatment
       rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) )   ! If zw10 < 33. => 0, else => 1  
@@ -856,4 +863,5 @@
          &       (1._wp - rgt33)*( 2.7_wp/zw10 + 0.142_wp + zw10/13.09_wp - 3.14807E-10*zw10**6) & ! zw10< 33.
          &      + rgt33         *      2.34   )                                                    ! zw10 >= 33.
+#endif
       !
       CALL wrk_dealloc( jpi,jpj, rgt33)
Index: /branches/NERC/dev_r5518_GO6_conserv_check_up/NEMOGCM/NEMO/TOP_SRC/trc.F90
===================================================================
--- /branches/NERC/dev_r5518_GO6_conserv_check_up/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 9207)
+++ /branches/NERC/dev_r5518_GO6_conserv_check_up/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 9208)
@@ -33,4 +33,9 @@
    !! --------------------------------------------------
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)               ::  trai           !: initial total tracer
+# if defined key_medusa && key_roam 
+   !! AXY (17/11/2017): elemental cycle initial totals
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)               ::  cycletot       !: initial elemental cycle total
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)               ::  cycletot2      !: initial elemental cycle total excl. halo in mpp_sum
+# endif
    REAL(wp), PUBLIC                                                ::  areatot        !: total volume 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:  )         ::  cvol           !: volume correction -degrad option- 
@@ -266,4 +271,7 @@
          &      cvol(jpi,jpj,jpk)     , rdttrc(jpk)           , trai(jptra)           ,       &
          &      ctrcnm(jptra)         , ctrcln(jptra)         , ctrcun(jptra)         ,       & 
+# if defined key_medusa && defined key_roam
+         &      cycletot(6), cycletot2(6)                                             ,       &
+# endif
          &      ln_trc_ini(jptra)     , ln_trc_wri(jptra)     , qsr_mean(jpi,jpj)     ,  STAT = trc_alloc  )  
 
Index: /branches/NERC/dev_r5518_GO6_conserv_check_up/NEMOGCM/NEMO/TOP_SRC/trcini.F90
===================================================================
--- /branches/NERC/dev_r5518_GO6_conserv_check_up/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 9207)
+++ /branches/NERC/dev_r5518_GO6_conserv_check_up/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 9208)
@@ -28,4 +28,5 @@
    USE trcini_idtra    ! idealize tracer initialisation
    USE trcini_medusa   ! MEDUSA   initialisation
+   USE par_medusa      ! MEDUSA   parameters (needed for elemental cycles)
    USE trcdta          ! initialisation from files
    USE daymod          ! calendar manager
@@ -35,5 +36,5 @@
    USE sbc_oce
    USE trcice          ! tracers in sea ice
- 
+   USE sms_medusa      ! MEDUSA   initialisation
    IMPLICIT NONE
    PRIVATE
@@ -62,5 +63,13 @@
       !!                or read data or analytical formulation
       !!---------------------------------------------------------------------
-      INTEGER ::   jk, jn, jl    ! dummy loop indices
+      INTEGER ::   ji, jj, jk, jn, jl    ! dummy loop indices
+# if defined key_medusa && defined key_roam
+      !! AXY (23/11/2017)
+      REAL(wp)                         :: zsum3d, zsum2d
+      REAL(wp)                         :: zq1, zq2, loc_vol, loc_area
+      REAL(wp), DIMENSION(6)           :: loc_cycletot3, loc_cycletot2
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztot3d
+      REAL(wp), DIMENSION(jpi,jpj)     :: ztot2d, carea
+# endif
       CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
@@ -98,4 +107,5 @@
       !                                                              ! total volume of the ocean 
       areatot = glob_sum( cvol(:,:,:) )
+      carea(:,:) = e1e2t(:,:) * tmask(:,:,1) 
 
       IF( lk_pisces  )       CALL trc_ini_pisces       ! PISCES  bio-model
@@ -192,4 +202,84 @@
       ENDIF
 
+# if defined key_medusa && defined key_roam
+      ! AXY (17/11/2017): calculate initial totals of elemental cycles
+      !
+      ! This is done in a very hard-wired way here; in future, this could be
+      ! replaced with loops and using a 2D array; one dimension would cover
+      ! the tracers, the other would be for the elements; each tracer would
+      ! have a factor for each element to say how much of that element was
+      ! in that tracer; for example, PHN would be 1.0 for N, xrfn for Fe and
+      ! xthetapn for C, with the other elements 0.0; the array entry for PHN
+      ! would then be (1. 0. xrfn xthetapn 0. 0.) for (N, Si, Fe, C, A, O2);
+      ! saving this for the next iteration
+      !
+      cycletot(:) = 0._wp
+      ! report elemental totals at initialisation as we go along
+      IF ( lwp ) WRITE(numout,*)
+      IF ( lwp ) WRITE(numout,*)    ' Elemental cycle totals: '
+      IF ( lwp ) CALL flush(numout)
+      ! nitrogen
+      ztot3d(:,:,:) = trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + &
+                      trn(:,:,:,jpzme) + trn(:,:,:,jpdet) + trn(:,:,:,jpdin)
+      ztot2d(:,:)   = zn_sed_n(:,:)
+      zsum3d        = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) )
+      zsum2d        = glob_sum( ztot2d(:,:) * carea(:,:) )
+      cycletot(1)   = zsum3d + zsum2d
+      IF ( lwp ) WRITE(numout,9010) 'nitrogen', zsum3d, zsum2d, cycletot(1)
+      IF ( lwp ) CALL flush(numout)
+      ! silicon
+      ztot3d(:,:,:) = trn(:,:,:,jppds) + trn(:,:,:,jpsil)
+      ztot2d(:,:)   = zn_sed_si(:,:)
+      zsum3d        = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) )
+      zsum2d        = glob_sum( ztot2d(:,:) * carea(:,:) )
+      cycletot(2)   = zsum3d + zsum2d
+      IF ( lwp ) WRITE(numout,9010) 'silicon', zsum3d, zsum2d, cycletot(2)
+      IF ( lwp ) CALL flush(numout)
+      ! iron
+      ztot3d(:,:,:) = ((trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + &
+                      trn(:,:,:,jpzme) + trn(:,:,:,jpdet)) * xrfn) + trn(:,:,:,jpfer)
+      ztot2d(:,:)   = zn_sed_fe(:,:)
+      zsum3d        = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) )
+      zsum2d        = glob_sum( ztot2d(:,:) * carea(:,:) )
+      cycletot(3)   = zsum3d + zsum2d
+      IF ( lwp ) WRITE(numout,9010) 'iron', zsum3d, zsum2d, cycletot(3)
+      IF ( lwp ) CALL flush(numout)
+      ! carbon (uses fixed C:N ratios on plankton tracers)
+      ztot3d(:,:,:) = (trn(:,:,:,jpphn) * xthetapn)  + (trn(:,:,:,jpphd) * xthetapd)  +  &
+                      (trn(:,:,:,jpzmi) * xthetazmi) + (trn(:,:,:,jpzme) * xthetazme) +  &
+                      trn(:,:,:,jpdtc) + trn(:,:,:,jpdic)
+      ztot2d(:,:)   = zn_sed_c(:,:) + zn_sed_ca(:,:)
+      zsum3d        = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) )
+      zsum2d        = glob_sum( ztot2d(:,:) * carea(:,:) )
+      cycletot(4)   = zsum3d + zsum2d
+      IF ( lwp ) WRITE(numout,9010) 'carbon', zsum3d, zsum2d, cycletot(4)
+      IF ( lwp ) CALL flush(numout)
+      ! alkalinity (note benthic correction)
+      ztot3d(:,:,:) = trn(:,:,:,jpalk)
+      ztot2d(:,:)   = zn_sed_ca(:,:) * 2._wp
+      zsum3d        = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) )
+      zsum2d        = glob_sum( ztot2d(:,:) * carea(:,:) )
+      cycletot(5)   = zsum3d + zsum2d
+      IF ( lwp ) WRITE(numout,9010) 'alkalinity', zsum3d, zsum2d, cycletot(5)
+      IF ( lwp ) CALL flush(numout)
+      ! oxygen (note no benthic)
+      ztot3d(:,:,:) = trn(:,:,:,jpoxy)
+      ztot2d(:,:)   = 0._wp
+      zsum3d        = glob_sum( ztot3d(:,:,:) * cvol(:,:,:) )
+      zsum2d        = glob_sum( ztot2d(:,:) * carea(:,:) )
+      cycletot(6)   = zsum3d + zsum2d
+      IF ( lwp ) WRITE(numout,9010) 'oxygen', zsum3d, zsum2d, cycletot(6)
+      IF ( lwp ) CALL flush(numout)
+      ! Check
+      zsum3d        = glob_sum( cvol(:,:,:) )
+      zsum2d        = glob_sum( carea(:,:) )      
+      IF ( lwp ) WRITE(numout,*)
+      IF ( lwp ) WRITE(numout,*) ' check : cvol    : ', zsum3d
+      IF ( lwp ) WRITE(numout,*) ' check : carea   : ', zsum2d
+      IF ( lwp ) WRITE(numout,*)
+      IF ( lwp ) CALL flush(numout)
+      !
+# endif
+
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) 'trc_init : passive tracer set up completed'
@@ -202,4 +292,7 @@
 
 9000  FORMAT(' tracer nb : ',i2,'      name :',a10,'      initial content :',e18.10)
+9010  FORMAT(' element:',a10,                     &
+             ' 3d sum:',e18.10,' 2d sum:',e18.10, &
+             ' total:',e18.10)
       !
       IF( nn_timing == 1 )   CALL timing_stop('trc_init')
Index: /branches/NERC/dev_r5518_GO6_conserv_check_up/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- /branches/NERC/dev_r5518_GO6_conserv_check_up/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 9207)
+++ /branches/NERC/dev_r5518_GO6_conserv_check_up/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 9208)
@@ -30,4 +30,5 @@
    USE daymod
    !! AXY (05/11/13): need these for MEDUSA to input/output benthic reservoirs
+   USE par_medusa
    USE sms_medusa
    USE trcsms_medusa
@@ -53,4 +54,7 @@
    PUBLIC   trc_rst_cal
    PUBLIC   trc_rst_stat
+#if defined key_medusa && defined key_roam
+   PUBLIC   trc_rst_conserve
+#endif
 
    !! * Substitutions
@@ -539,4 +543,7 @@
       IF( kt == nitrst ) THEN
           CALL trc_rst_stat            ! statistics
+#if defined key_medusa && defined key_roam
+          CALL trc_rst_conserve        ! conservation check
+#endif
           CALL iom_close( numrtw )     ! close the restart file (only at last time step)
 #if ! defined key_trdmxl_trc
@@ -705,4 +712,151 @@
 
 
+# if defined key_medusa && defined key_roam
+   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, loc_vol, loc_are
+      REAL(wp) :: zq1, zq2, loc_vol, loc_area
+      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(:,:) * fse3t_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 ) WRITE(numout,*)
+      IF ( lwp ) WRITE(numout,*) ' check : cvol    : ', zsum3d
+      IF ( lwp ) WRITE(numout,*) ' check : carea   : ', zsum2d
+      IF ( lwp ) WRITE(numout,*)
+      IF ( lwp ) CALL flush(numout)
+      !
+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 
+# endif
+
+
 #else
    !!----------------------------------------------------------------------
