Index: branches/UKMO/dev_r5518_GO6M_dev/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6M_dev/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 8199)
+++ branches/UKMO/dev_r5518_GO6M_dev/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 8201)
@@ -27,5 +27,21 @@
    USE trcnam_trp
    USE iom
+   USE ioipsl, ONLY : ju2ymds    ! for calendar
    USE daymod
+   !! AXY (05/11/13): need these for MEDUSA to input/output benthic reservoirs
+   USE sms_medusa
+   USE trcsms_medusa
+   !!
+#if defined key_idtra
+   USE trcsms_idtra
+#endif
+   !!
+#if defined key_cfc
+   USE trcsms_cfc
+#endif
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE sbc_oce, ONLY: lk_oasis 
+   USE oce,     ONLY: CO2Flux_out_cpl, DMS_out_cpl, chloro_out_cpl  !! Coupling variable
+
    IMPLICIT NONE
    PRIVATE
@@ -35,4 +51,7 @@
    PUBLIC   trc_rst_wri       ! called by ???
    PUBLIC   trc_rst_cal
+   PUBLIC   trc_rst_stat
+   PUBLIC   trc_rst_dia_stat
+   PUBLIC   trc_rst_tra_stat
 
    !! * Substitutions
@@ -48,4 +67,7 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt       ! number of iteration
+      INTEGER             ::   iyear, imonth, iday
+      REAL (wp)           ::   zsec
+      REAL (wp)           ::   zfjulday
       !
       CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character
@@ -78,7 +100,26 @@
       ! except if we write tracer restart files every tracer time step or if a tracer restart file was writen at nitend - 2*nn_dttrc + 1
       IF( kt == nitrst - 2*nn_dttrc .OR. nstock == nn_dttrc .OR. ( kt == nitend - nn_dttrc .AND. .NOT. lrst_trc ) ) THEN
-         ! beware of the format used to write kt (default is i8.8, that should be large enough)
-         IF( nitrst > 1.0e9 ) THEN   ;   WRITE(clkt,*       ) nitrst
-         ELSE                        ;   WRITE(clkt,'(i8.8)') nitrst
+         IF ( ln_rstdate ) THEN
+            !! JPALM -- 22-12-2015 -- modif to get the good date on restart trc file name
+            !!                     -- the condition to open the rst file is not the same than for the dynamic rst.
+            !!                     -- here it - for an obscure reason - is open 2 time-step before the restart writing process
+            !!                     instead of 1.
+            !!                     -- i am not sure if someone forgot +1 in the if loop condition as
+            !!                     it is writen in all comments nitrst -2*nn_dttrc + 1 and the condition is 
+            !!                     nitrst - 2*nn_dttrc
+            !!                     -- nevertheless we didn't wanted to broke something already working 
+            !!                     and just adapted the part we added.
+            !!                     -- So instead of calling ju2ymds( fjulday + (rdttra(1)) 
+            !!                     we call ju2ymds( fjulday + (2*rdttra(1)) 
+            !!--------------------------------------------------------------------      
+            zfjulday = fjulday + (2*rdttra(1)) / rday
+            IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday )   zfjulday = REAL(NINT(zfjulday),wp)   ! avoid truncation error
+            CALL ju2ymds( zfjulday + (2*rdttra(1)) / rday, iyear, imonth, iday, zsec )
+            WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday
+         ELSE
+            ! beware of the format used to write kt (default is i8.8, that should be large enough)
+            IF( nitrst > 1.0e9 ) THEN   ;   WRITE(clkt,*       ) nitrst
+            ELSE                        ;   WRITE(clkt,'(i8.8)') nitrst
+            ENDIF
          ENDIF
          ! create the file
@@ -101,5 +142,7 @@
       !! ** purpose  :   read passive tracer fields in restart files
       !!----------------------------------------------------------------------
-      INTEGER  ::  jn     
+      INTEGER  ::  jn, jl     
+      !! AXY (05/11/13): temporary variables
+      REAL(wp) ::    fq0,fq1,fq2
 
       !!----------------------------------------------------------------------
@@ -112,9 +155,226 @@
       DO jn = 1, jptra
          CALL iom_get( numrtr, jpdom_autoglo, 'TRN'//ctrcnm(jn), trn(:,:,:,jn) )
+         trn(:,:,:,jn) = trn(:,:,:,jn) * tmask(:,:,:)
       END DO
 
       DO jn = 1, jptra
          CALL iom_get( numrtr, jpdom_autoglo, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) )
-      END DO
+         trb(:,:,:,jn) = trb(:,:,:,jn) * tmask(:,:,:)
+      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) ...'
+      call trc_rst_dia_stat(zn_sed_n(:,:), 'Sediment  N')
+      call trc_rst_dia_stat(zn_sed_fe(:,:), 'Sediment Fe')
+      call trc_rst_dia_stat(zn_sed_si(:,:), 'Sediment Si')
+      call trc_rst_dia_stat(zn_sed_c(:,:), 'Sediment C')
+      call trc_rst_dia_stat(zn_sed_ca(:,:), 'Sediment Ca')
+      !!
+      !! AXY (07/07/15): read in temporally averaged fields for DMS
+      !!                 calculations
+      !!
+      IF( iom_varid( numrtr, 'B_DMS_CHN', ldstop = .FALSE. ) > 0 ) THEN
+         !! YES; in which case read them
+         !!
+         IF(lwp) WRITE(numout,*) ' MEDUSA averaged properties for DMS present - reading in ...'
+         CALL iom_get( numrtr, jpdom_autoglo, 'B_DMS_CHN',  zb_dms_chn(:,:)  )
+         CALL iom_get( numrtr, jpdom_autoglo, 'N_DMS_CHN',  zn_dms_chn(:,:)  )
+         CALL iom_get( numrtr, jpdom_autoglo, 'B_DMS_CHD',  zb_dms_chd(:,:)  )
+         CALL iom_get( numrtr, jpdom_autoglo, 'N_DMS_CHD',  zn_dms_chd(:,:)  )
+         CALL iom_get( numrtr, jpdom_autoglo, 'B_DMS_MLD',  zb_dms_mld(:,:)  )
+         CALL iom_get( numrtr, jpdom_autoglo, 'N_DMS_MLD',  zn_dms_mld(:,:)  )
+         CALL iom_get( numrtr, jpdom_autoglo, 'B_DMS_QSR',  zb_dms_qsr(:,:)  )
+         CALL iom_get( numrtr, jpdom_autoglo, 'N_DMS_QSR',  zn_dms_qsr(:,:)  )
+         CALL iom_get( numrtr, jpdom_autoglo, 'B_DMS_DIN',  zb_dms_din(:,:)  )
+         CALL iom_get( numrtr, jpdom_autoglo, 'N_DMS_DIN',  zn_dms_din(:,:)  )
+      ELSE
+         !! NO; in which case set them to zero
+         !!
+         IF(lwp) WRITE(numout,*) ' MEDUSA averaged properties for DMS absent - setting to zero ...'
+         zb_dms_chn(:,:)  = 0.0   !! CHN
+         zn_dms_chn(:,:)  = 0.0
+         zb_dms_chd(:,:)  = 0.0   !! CHD
+         zn_dms_chd(:,:)  = 0.0
+         zb_dms_mld(:,:)  = 0.0   !! MLD
+         zn_dms_mld(:,:)  = 0.0
+         zb_dms_qsr(:,:)  = 0.0   !! QSR
+         zn_dms_qsr(:,:)  = 0.0
+         zb_dms_din(:,:)  = 0.0   !! DIN
+         zn_dms_din(:,:)  = 0.0
+      ENDIF
+      !!  
+      !! JPALM 14-06-2016 -- add CO2 flux and DMS surf through the restart
+      !!                  -- needed for the coupling with atm
+      IF( iom_varid( numrtr, 'N_DMS_srf', ldstop = .FALSE. ) > 0 ) THEN
+         IF(lwp) WRITE(numout,*) 'DMS surf concentration - reading in ...'
+         CALL iom_get( numrtr, jpdom_autoglo, 'B_DMS_srf',  zb_dms_srf(:,:)  )
+         CALL iom_get( numrtr, jpdom_autoglo, 'N_DMS_srf',  zn_dms_srf(:,:)  )
+      ELSE
+         IF(lwp) WRITE(numout,*) 'DMS surf concentration - setting to zero ...'
+         zb_dms_srf(:,:)  = 0.0   !! DMS
+         zn_dms_srf(:,:)  = 0.0
+      ENDIF
+      IF (lk_oasis) THEN
+         DMS_out_cpl(:,:) = zn_dms_srf(:,:)        !! Coupling variable
+      END IF
+      !!
+      IF( iom_varid( numrtr, 'B_CO2_flx', ldstop = .FALSE. ) > 0 ) THEN
+         IF(lwp) WRITE(numout,*) 'CO2 air-sea flux - reading in ...'
+         CALL iom_get( numrtr, jpdom_autoglo, 'B_CO2_flx',  zb_co2_flx(:,:)  )
+         CALL iom_get( numrtr, jpdom_autoglo, 'N_CO2_flx',  zn_co2_flx(:,:)  )
+      ELSE
+         IF(lwp) WRITE(numout,*) 'CO2 air-sea flux - setting to zero ...'
+         zb_co2_flx(:,:)  = 0.0   !! CO2 flx
+         zn_co2_flx(:,:)  = 0.0
+      ENDIF
+      IF (lk_oasis) THEN
+         CO2Flux_out_cpl(:,:) =  zn_co2_flx(:,:)   !! Coupling variable
+      END IF
+      !!
+      !! JPALM 02-06-2017 -- in complement to DMS surf 
+      !!                  -- the atm model needs surf Chl 
+      !!                     as proxy of org matter from the ocean
+      !!                  -- needed for the coupling with atm
+      IF( iom_varid( numrtr, 'N_CHL_srf', ldstop = .FALSE. ) > 0 ) THEN
+         IF(lwp) WRITE(numout,*) 'Chl surf concentration - reading in ...'
+         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
+      ENDIF
+      IF (lk_oasis) THEN
+         chloro_out_cpl(:,:) = zn_chl_srf(:,:)        !! Coupling variable
+      END IF
+      !!
+      !! calculate stats on these fields
+      IF(lwp) WRITE(numout,*) ' MEDUSA averaged properties for DMS stats (min, max, sum) ...'
+      call trc_rst_dia_stat(zn_dms_chn(:,:), 'DMS, CHN')
+      call trc_rst_dia_stat(zn_dms_chd(:,:), 'DMS, CHD')
+      call trc_rst_dia_stat(zn_dms_mld(:,:), 'DMS, MLD')
+      call trc_rst_dia_stat(zn_dms_qsr(:,:), 'DMS, QSR')
+      call trc_rst_dia_stat(zn_dms_din(:,:), 'DMS, DIN')
+      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')
+      !!  
+      !! JPALM 14-06-2016 -- add Carbonate chenistry variables through the restart
+      !!                  -- needed for monthly call of carb-chem routine and better reproducibility
+# if defined key_roam
+      IF( iom_varid( numrtr, 'pH_3D', ldstop = .FALSE. ) > 0 ) THEN
+         IF(lwp) WRITE(numout,*) 'Carbonate chem variable - reading in ...'
+         CALL iom_get( numrtr, jpdom_autoglo, 'pH_3D'   ,  f3_pH(:,:,:)     )
+         CALL iom_get( numrtr, jpdom_autoglo, 'h2CO3_3D',  f3_h2co3(:,:,:)  )
+         CALL iom_get( numrtr, jpdom_autoglo, 'hCO3_3D' ,  f3_hco3(:,:,:)   )
+         CALL iom_get( numrtr, jpdom_autoglo, 'CO3_3D'  ,  f3_co3(:,:,:)    )
+         CALL iom_get( numrtr, jpdom_autoglo, 'omcal_3D',  f3_omcal(:,:,:)  )
+         CALL iom_get( numrtr, jpdom_autoglo, 'omarg_3D',  f3_omarg(:,:,:)  )
+         CALL iom_get( numrtr, jpdom_autoglo, 'CCD_CAL' ,  f2_ccd_cal(:,:)  )
+         CALL iom_get( numrtr, jpdom_autoglo, 'CCD_ARG' ,  f2_ccd_arg(:,:)  )
+         !!
+         IF(lwp) WRITE(numout,*) ' MEDUSA averaged Carb-chem stats (min, max, sum) ...'
+      call trc_rst_dia_stat( f3_pH(:,:,1)   ,'pH 3D surf')
+      call trc_rst_dia_stat( f3_h2co3(:,:,1),'h2CO3 3D surf')
+      call trc_rst_dia_stat( f3_hco3(:,:,1) ,'hCO3 3D surf' )
+      call trc_rst_dia_stat( f3_co3(:,:,1)  ,'CO3 3D surf' )
+      call trc_rst_dia_stat( f3_omcal(:,:,1),'omcal 3D surf')
+      call trc_rst_dia_stat( f3_omarg(:,:,1),'omarg 3D surf')
+      call trc_rst_dia_stat( f2_ccd_cal(:,:),'CCD_CAL')
+      call trc_rst_dia_stat( f2_ccd_arg(:,:),'CCD_ARG')
+
+      ELSE
+         IF(lwp) WRITE(numout,*) 'WARNING : No Carbonate-chem variable in the restart.... '
+         IF(lwp) WRITE(numout,*) 'Is not a problem if start a month, but may be very problematic if not '
+         IF(lwp) WRITE(numout,*) 'Check if   mod(kt*rdt,2592000) == rdt' 
+         IF(lwp) WRITE(numout,*) 'Or don t start from uncomplete restart...' 
+      ENDIF
+# endif
+
+
+#endif
+      !
+#if defined key_idtra
+      !! JPALM -- 05-01-2016 -- mv idtra and CFC restart reading and 
+      !!                        writting here undre their key.
+      !!                        problems in CFC restart, maybe because of this...
+      !!                        and pb in idtra diag or diad-restart writing.
+      !!----------------------------------------------------------------------
+      IF( iom_varid( numrtr, 'qint_IDTRA', ldstop = .FALSE. ) > 0 ) THEN
+         !! YES; in which case read them
+         !!
+         IF(lwp) WRITE(numout,*) ' IDTRA averaged properties present - reading in ...'
+         CALL iom_get( numrtr, jpdom_autoglo, 'qint_IDTRA',  qint_idtra(:,:,1)  )
+      ELSE
+         !! NO; in which case set them to zero
+         !!
+         IF(lwp) WRITE(numout,*) ' IDTRA averaged properties absent - setting to zero ...'
+         qint_idtra(:,:,1)  = 0.0   !! CHN
+      ENDIF
+      !!
+      !! calculate stats on these fields
+      IF(lwp) WRITE(numout,*) ' IDTRA averaged properties stats (min, max, sum) ...'
+      call trc_rst_dia_stat(qint_idtra(:,:,1), 'qint_IDTRA')
+#endif
+      !
+#if defined key_cfc
+      DO jl = 1, jp_cfc
+         jn = jp_cfc0 + jl - 1
+         IF( iom_varid( numrtr, 'qint_'//ctrcnm(jn), ldstop = .FALSE. ) > 0 ) THEN
+            !! YES; in which case read them
+            !!
+            IF(lwp) WRITE(numout,*) ' CFC averaged properties present - reading in ...'
+            CALL iom_get( numrtr, jpdom_autoglo, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jl) )
+         ELSE
+            !! NO; in which case set them to zero
+            !!
+            IF(lwp) WRITE(numout,*) ' CFC averaged properties absent - setting to zero ...'
+            qint_cfc(:,:,jn)  = 0.0   !! CHN
+         ENDIF
+         !!
+         !! calculate stats on these fields
+         IF(lwp) WRITE(numout,*) ' CFC averaged properties stats (min, max, sum) ...'
+         call trc_rst_dia_stat(qint_cfc(:,:,jl), 'qint_'//ctrcnm(jn))
+      END DO
+#endif
       !
    END SUBROUTINE trc_rst_read
@@ -128,6 +388,8 @@
       INTEGER, INTENT( in ) ::   kt    ! ocean time-step index
       !!
-      INTEGER  :: jn
+      INTEGER  :: jn, jl
       REAL(wp) :: zarak0
+      !! AXY (05/11/13): temporary variables
+      REAL(wp) ::    fq0,fq1,fq2
       !!----------------------------------------------------------------------
       !
@@ -142,5 +404,129 @@
          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) ...'
+      call trc_rst_dia_stat(zn_sed_n(:,:), 'Sediment  N')
+      call trc_rst_dia_stat(zn_sed_fe(:,:), 'Sediment Fe')
+      call trc_rst_dia_stat(zn_sed_si(:,:), 'Sediment Si')
+      call trc_rst_dia_stat(zn_sed_c(:,:), 'Sediment C')
+      call trc_rst_dia_stat(zn_sed_ca(:,:), 'Sediment Ca')
+      !!
+      !! AXY (07/07/15): write out temporally averaged fields for DMS
+      !!                 calculations
+      !!
+      IF(lwp) WRITE(numout,*) ' MEDUSA averaged properties for DMS - writing out ...'
+      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_CHN',  zb_dms_chn(:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_CHN',  zn_dms_chn(:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_CHD',  zb_dms_chd(:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_CHD',  zn_dms_chd(:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_MLD',  zb_dms_mld(:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_MLD',  zn_dms_mld(:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_QSR',  zb_dms_qsr(:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_QSR',  zn_dms_qsr(:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_DIN',  zb_dms_din(:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_DIN',  zn_dms_din(:,:)  )
+         !! JPALM 14-06-2016 -- add CO2 flux and DMS surf through the restart
+         !!                  -- needed for the coupling with atm
+      CALL iom_rstput( kt, nitrst, numrtw, 'B_DMS_srf',  zb_dms_srf(:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'N_DMS_srf',  zn_dms_srf(:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'B_CO2_flx',  zb_co2_flx(:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'N_CO2_flx',  zn_co2_flx(:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'N_CHL_srf',  zn_chl_srf(:,:)  )
+      !!
+      !! calculate stats on these fields
+      IF(lwp) WRITE(numout,*) ' MEDUSA averaged properties for DMS stats (min, max, sum) ...'
+      call trc_rst_dia_stat(zn_dms_chn(:,:), 'DMS, CHN')
+      call trc_rst_dia_stat(zn_dms_chd(:,:), 'DMS, CHD')
+      call trc_rst_dia_stat(zn_dms_mld(:,:), 'DMS, MLD')
+      call trc_rst_dia_stat(zn_dms_qsr(:,:), 'DMS, QSR')
+      call trc_rst_dia_stat(zn_dms_din(:,:), 'DMS, DIN')
+      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(lwp) WRITE(numout,*) ' MEDUSA averaged prop. for dust and iron dep.'
+      call trc_rst_dia_stat(dust(:,:), 'Dust dep')
+      call trc_rst_dia_stat(zirondep(:,:), 'Iron dep')
+      !! 
+      !!  
+      !! JPALM 14-06-2016 -- add Carbonate chenistry variables through the restart
+      !!                  -- needed for monthly call of carb-chem routine and better reproducibility
+# if defined key_roam
+      IF(lwp) WRITE(numout,*) 'Carbonate chem variable - writing out ...'
+      CALL iom_rstput( kt, nitrst, numrtw, 'pH_3D'   ,  f3_pH(:,:,:)     )
+      CALL iom_rstput( kt, nitrst, numrtw, 'h2CO3_3D',  f3_h2co3(:,:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'hCO3_3D' ,  f3_hco3(:,:,:)   )
+      CALL iom_rstput( kt, nitrst, numrtw, 'CO3_3D'  ,  f3_co3(:,:,:)    )
+      CALL iom_rstput( kt, nitrst, numrtw, 'omcal_3D',  f3_omcal(:,:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'omarg_3D',  f3_omarg(:,:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'CCD_CAL' ,  f2_ccd_cal(:,:)  )
+      CALL iom_rstput( kt, nitrst, numrtw, 'CCD_ARG' ,  f2_ccd_arg(:,:)  )
+      !!
+      IF(lwp) WRITE(numout,*) ' MEDUSA averaged Carb-chem stats (min, max, sum) ...'
+      call trc_rst_dia_stat( f3_pH(:,:,1)   ,'pH 3D surf')
+      call trc_rst_dia_stat( f3_h2co3(:,:,1),'h2CO3 3D surf')
+      call trc_rst_dia_stat( f3_hco3(:,:,1) ,'hCO3 3D surf' )
+      call trc_rst_dia_stat( f3_co3(:,:,1)  ,'CO3 3D surf' )
+      call trc_rst_dia_stat( f3_omcal(:,:,1),'omcal 3D surf')
+      call trc_rst_dia_stat( f3_omarg(:,:,1),'omarg 3D surf')
+      call trc_rst_dia_stat( f2_ccd_cal(:,:),'CCD_CAL')
+      call trc_rst_dia_stat( f2_ccd_arg(:,:),'CCD_ARG')
+      !!
+# endif
+!!
+#endif
+      !
+#if defined key_idtra
+      !! JPALM -- 05-01-2016 -- mv idtra and CFC restart reading and 
+      !!                        writting here undre their key.
+      !!                        problems in CFC restart, maybe because of this...
+      !!                        and pb in idtra diag or diad-restart writing.
+      !!----------------------------------------------------------------------
+      IF(lwp) WRITE(numout,*) ' IDTRA averaged properties - writing out ...'
+      CALL iom_rstput( kt, nitrst, numrtw, 'qint_IDTRA',  qint_idtra(:,:,1) )
+      !!
+      !! calculate stats on these fields
+      IF(lwp) WRITE(numout,*) ' IDTRA averaged properties stats (min, max, sum) ...'
+      call trc_rst_dia_stat(qint_idtra(:,:,1), 'qint_IDTRA')
+#endif
+      !
+#if defined key_cfc
+      DO jl = 1, jp_cfc
+         jn = jp_cfc0 + jl - 1
+         IF(lwp) WRITE(numout,*) ' CFC averaged properties - writing out ...'
+         CALL iom_rstput( kt, nitrst, numrtw, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jl) )
+         !!
+         !! calculate stats on these fields
+         IF(lwp) WRITE(numout,*) ' CFC averaged properties stats (min, max, sum) ...'
+         call trc_rst_dia_stat(qint_cfc(:,:,jl), 'qint_'//ctrcnm(jn))
+      END DO
+#endif
+      !
+
       IF( kt == nitrst ) THEN
           CALL trc_rst_stat            ! statistics
@@ -304,9 +690,84 @@
          IF(lwp) WRITE(numout,9000) jn, TRIM( ctrcnm(jn) ), zmean, zmin, zmax, zdrift
       END DO
-      WRITE(numout,*) 
+      IF(lwp) WRITE(numout,*) 
 9000  FORMAT(' tracer nb :',i2,'    name :',a10,'    mean :',e18.10,'    min :',e18.10, &
       &      '    max :',e18.10,'    drift :',e18.10, ' %')
       !
    END SUBROUTINE trc_rst_stat
+
+
+   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( lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '           ----SURFACE TRA STAT----             '
+         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
+      REAL(wp) :: ztraf, zmin, zmax, zmean, areasf
+      REAL(wp), DIMENSION(jpi,jpj) :: zvol
+      !!----------------------------------------------------------------------
+
+      IF( lwp )  WRITE(numout,*) 'STAT- ', names
+      !
+      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
+      zmean  = ztraf / areasf
+      IF(lwp) WRITE(numout,9002) TRIM( names ), zmean, zmin, zmax
+      !
+      IF(lwp) WRITE(numout,*)
+9002  FORMAT(' tracer name :',a10,'    mean :',e18.10,'    min :',e18.10, &
+      &      '    max :',e18.10 )
+      !
+   END SUBROUTINE trc_rst_dia_stat
+
 
 #else
