Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 10302)
@@ -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
@@ -43,4 +44,6 @@
    USE sbc_oce, ONLY: lk_oasis 
    USE oce,     ONLY: CO2Flux_out_cpl, DMS_out_cpl, chloro_out_cpl  !! Coupling variable
+   USE trcstat
+   USE obs_const, ONLY: obfillflt  ! Observation operator fill value
 
    IMPLICIT NONE
@@ -52,6 +55,7 @@
    PUBLIC   trc_rst_cal
    PUBLIC   trc_rst_stat
-   PUBLIC   trc_rst_dia_stat
-   PUBLIC   trc_rst_tra_stat
+#if defined key_medusa && defined key_roam
+   PUBLIC   trc_rst_conserve
+#endif
 
    !! * Substitutions
@@ -277,13 +281,17 @@
       !!                     as proxy of org matter from the ocean
       !!                  -- needed for the coupling with atm
+      !!       07-12-2017 -- To make things cleaner, we want to store an  
+      !!                     unscaled Chl field in the restart and only 
+      !!                     scale it when reading it in.
+
       IF( iom_varid( numrtr, 'N_CHL_srf', ldstop = .FALSE. ) > 0 ) THEN
-         IF(lwp) WRITE(numout,*) 'Chl surf concentration - reading in ...'
+         IF(lwp) WRITE(numout,*) 'Chl cpl concentration - reading in ... - scale by ', scl_chl
          CALL iom_get( numrtr, jpdom_autoglo, 'N_CHL_srf',  zn_chl_srf(:,:)  )
       ELSE
-         IF(lwp) WRITE(numout,*) 'Chl surf concentration - setting to zero ...'
-         zn_chl_srf(:,:)  = (trn(:,:,1,jpchn) + trn(:,:,1,jpchd)) * 1.E-6
+         IF(lwp) WRITE(numout,*) 'set Chl coupled concentration - scaled by ', scl_chl
+         zn_chl_srf(:,:)  = MAX( 0.0, (trn(:,:,1,jpchn) + trn(:,:,1,jpchd)) * 1.E-6 )
       ENDIF
       IF (lk_oasis) THEN
-         chloro_out_cpl(:,:) = zn_chl_srf(:,:)        !! Coupling variable
+         chloro_out_cpl(:,:) = zn_chl_srf(:,:) * scl_chl        !! Coupling variable
       END IF
       !!
@@ -297,5 +305,7 @@
       call trc_rst_dia_stat(zn_dms_srf(:,:), 'DMS surf')
       call trc_rst_dia_stat(zn_co2_flx(:,:), 'CO2 flux')
-      call trc_rst_dia_stat(zn_chl_srf(:,:), 'CHL surf')
+      IF (lk_oasis) THEN
+         call trc_rst_dia_stat(chloro_out_cpl(:,:), 'CHL  cpl')
+      END IF
       !!  
       !! JPALM 14-06-2016 -- add Carbonate chenistry variables through the restart
@@ -329,6 +339,34 @@
          IF(lwp) WRITE(numout,*) 'Or don t start from uncomplete restart...' 
       ENDIF
+      !
+      IF ( ln_foam_medusa ) THEN
+         !! 2D fields of pCO2 and fCO2 for observation operator on first timestep
+         IF( iom_varid( numrtr, 'PCO2W', ldstop = .FALSE. ) > 0 ) THEN
+            IF(lwp) WRITE(numout,*) ' MEDUSA pCO2 present - reading in ...'
+            CALL iom_get( numrtr, jpdom_autoglo, 'PCO2W',  f2_pco2w(:,:)  )
+            CALL iom_get( numrtr, jpdom_autoglo, 'FCO2W',  f2_fco2w(:,:)  )
+         ELSE
+            IF(lwp) WRITE(numout,*) ' MEDUSA pCO2 absent - setting to fill ...'
+            f2_pco2w(:,:) = obfillflt * tmask(:,:,1)
+            f2_fco2w(:,:) = obfillflt * tmask(:,:,1)
+         ENDIF
+      ENDIF
 # endif
-
+      IF ( ln_foam_medusa ) THEN
+         !! Fields for ocean colour assimilation on first timestep
+         IF( iom_varid( numrtr, 'pgrow_avg', ldstop = .FALSE. ) > 0 ) THEN
+            IF(lwp) WRITE(numout,*) ' MEDUSA pgrow_avg present - reading in ...'
+            CALL iom_get( numrtr, jpdom_autoglo, 'pgrow_avg',  pgrow_avg(:,:)  )
+            CALL iom_get( numrtr, jpdom_autoglo, 'ploss_avg',  ploss_avg(:,:)  )
+            CALL iom_get( numrtr, jpdom_autoglo, 'phyt_avg',   phyt_avg(:,:)   )
+            CALL iom_get( numrtr, jpdom_autoglo, 'mld_max',    mld_max(:,:)    )
+         ELSE
+            IF(lwp) WRITE(numout,*) ' MEDUSA pgrow_avg absent - setting to zero ...'
+            pgrow_avg(:,:) = 0.0
+            ploss_avg(:,:) = 0.0
+            phyt_avg(:,:)  = 0.0
+            mld_max(:,:)   = 0.0
+         ENDIF
+      ENDIF
 
 #endif
@@ -369,5 +407,5 @@
             !!
             IF(lwp) WRITE(numout,*) ' CFC averaged properties absent - setting to zero ...'
-            qint_cfc(:,:,jn)  = 0.0   !! CHN
+            qint_cfc(:,:,jl)  = 0.0   !! CHN
          ENDIF
          !!
@@ -457,4 +495,7 @@
       CALL iom_rstput( kt, nitrst, numrtw, 'B_CO2_flx',  zb_co2_flx(:,:)  )
       CALL iom_rstput( kt, nitrst, numrtw, 'N_CO2_flx',  zn_co2_flx(:,:)  )
+      !! JPALM 07-12-2017 -- To make things cleaner, we want to store an  
+      !!                     unscaled Chl field in the restart and only 
+      !!                     scale it when reading it in.
       CALL iom_rstput( kt, nitrst, numrtw, 'N_CHL_srf',  zn_chl_srf(:,:)  )
       !!
@@ -468,5 +509,5 @@
       call trc_rst_dia_stat(zn_dms_srf(:,:), 'DMS surf')
       call trc_rst_dia_stat(zn_co2_flx(:,:), 'CO2 flux')
-      call trc_rst_dia_stat(zn_chl_srf(:,:), 'CHL surf')
+      call trc_rst_dia_stat(zn_chl_srf(:,:), 'unscaled CHL cpl')
       !!
       IF(lwp) WRITE(numout,*) ' MEDUSA averaged prop. for dust and iron dep.'
@@ -498,5 +539,19 @@
       call trc_rst_dia_stat( f2_ccd_arg(:,:),'CCD_ARG')
       !!
+      IF ( ln_foam_medusa ) THEN
+         !! Fields for observation operator on first timestep
+         IF(lwp) WRITE(numout,*) ' MEDUSA OBS fields - writing out ...'
+         CALL iom_rstput( kt, nitrst, numrtw, 'PCO2W', f2_pco2w(:,:)  )
+         CALL iom_rstput( kt, nitrst, numrtw, 'FCO2W', f2_fco2w(:,:)  )
+      ENDIF
 # endif
+      IF ( ln_foam_medusa ) THEN
+         !! Fields for assimilation on first timestep
+         IF(lwp) WRITE(numout,*) ' MEDUSA ASM fields - writing out ...'
+         CALL iom_rstput( kt, nitrst, numrtw, 'pgrow_avg', pgrow_avg(:,:) )
+         CALL iom_rstput( kt, nitrst, numrtw, 'ploss_avg', ploss_avg(:,:) )
+         CALL iom_rstput( kt, nitrst, numrtw, 'phyt_avg',  phyt_avg(:,:)  )
+         CALL iom_rstput( kt, nitrst, numrtw, 'mld_max',   mld_max(:,:)   )
+      ENDIF
 !!
 #endif
@@ -531,4 +586,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
@@ -697,88 +755,149 @@
 
 
-   SUBROUTINE trc_rst_tra_stat
-      !!----------------------------------------------------------------------
-      !!                    ***  trc_rst_tra_stat  ***
-      !!
-      !! ** purpose  :   Compute tracers statistics - check where crazy values appears
-      !!----------------------------------------------------------------------
-      INTEGER  :: jk, jn
-      REAL(wp) :: ztraf, zmin, zmax, zmean, zdrift, areasf
-      REAL(wp), DIMENSION(jpi,jpj) :: zvol
-      !!----------------------------------------------------------------------
-
+# 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
+      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 ) THEN 
          WRITE(numout,*)
-         WRITE(numout,*) '           ----SURFACE TRA STAT----             '
+         WRITE(numout,*) ' check : cvol    : ', zsum3d
+         WRITE(numout,*) ' check : carea   : ', zsum2d
          WRITE(numout,*)
       ENDIF
       !
-      zvol(:,:) = e1e2t(:,:) * fse3t_a(:,:,1) * tmask(:,:,1)
-      areasf = glob_sum(zvol(:,:))
-      DO jn = 1, jptra
-         ztraf = glob_sum( tra(:,:,1,jn) * zvol(:,:) )
-         zmin  = MINVAL( tra(:,:,1,jn), mask= ((tmask(:,:,1).NE.0.)) )
-         zmax  = MAXVAL( tra(:,:,1,jn), mask= ((tmask(:,:,1).NE.0.)) )
-         IF( lk_mpp ) THEN
-            CALL mpp_min( zmin )      ! min over the global domain
-            CALL mpp_max( zmax )      ! max over the global domain
-         END IF
-         zmean  = ztraf / areasf
-         IF(lwp) WRITE(numout,9001) jn, TRIM( ctrcnm(jn) ), zmean, zmin, zmax
-      END DO
-      IF(lwp) WRITE(numout,*)
-9001  FORMAT(' tracer nb :',i2,'    name :',a10,'    mean :',e18.10,'    min :',e18.10, &
-      &      '    max :',e18.10)
-      !
-   END SUBROUTINE trc_rst_tra_stat
-
-
-
-   SUBROUTINE trc_rst_dia_stat( dgtr, names)
-      !!----------------------------------------------------------------------
-      !!                    ***  trc_rst_dia_stat  ***
-      !!
-      !! ** purpose  :   Compute tracers statistics
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) ::   dgtr      ! 2D diag var
-      CHARACTER(len=*)             , INTENT(in) ::   names     ! 2D diag name
-      !!---------------------------------------------------------------------
-      INTEGER  :: jk, jn
-      CHARACTER (LEN=18) :: text_zmean
-      REAL(wp) :: ztraf, zmin, zmax, zmean, areasf
-      REAL(wp), DIMENSION(jpi,jpj) :: zvol
-      !!----------------------------------------------------------------------
-
-      IF( lwp )  WRITE(numout,*) 'STAT- ', names
-      
-      ! fse3t_a will be undefined at the start of a run, but this routine
-      ! may be called at any stage! Hence we MUST make sure it is 
-      ! initialised to zero when allocated to enable us to test for 
-      ! zero content here and avoid potentially dangerous and non-portable 
-      ! operations (e.g. divide by zero, global sums of junk values etc.)   
-      zvol(:,:) = e1e2t(:,:) * fse3t_a(:,:,1) * tmask(:,:,1)
-      ztraf = glob_sum( dgtr(:,:) * zvol(:,:) )
-      !! areasf = glob_sum(e1e2t(:,:) * tmask(:,:,1) )
-      areasf = glob_sum(zvol(:,:))
-      zmin  = MINVAL( dgtr(:,:), mask= ((tmask(:,:,1).NE.0.)) )
-      zmax  = MAXVAL( dgtr(:,:), mask= ((tmask(:,:,1).NE.0.)) )
-      IF( lk_mpp ) THEN
-         CALL mpp_min( zmin )      ! min over the global domain
-         CALL mpp_max( zmax )      ! max over the global domain
-      END IF
-
-      text_zmean = "N/A"
-      ! Avoid divide by zero. areasf must be positive.
-      IF  (areasf > 0.0) THEN 
-         zmean = ztraf / areasf
-         WRITE(text_zmean,'(e18.10)') zmean
-      ENDIF
-
-      IF(lwp) WRITE(numout,9002) TRIM( names ), text_zmean, zmin, zmax
-
-  9002  FORMAT(' tracer name :',A,'    mean :',A,'    min :',e18.10, &
-      &      '    max :',e18.10 )
-      !
-   END SUBROUTINE trc_rst_dia_stat
+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
 
 
