Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/C14b/trcsms_c14b.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/C14b/trcsms_c14b.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/C14b/trcsms_c14b.F90	(revision 10302)
@@ -330,4 +330,6 @@
 CONTAINS
    SUBROUTINE trc_sms_c14b( kt )       ! Empty routine
+      IMPLICIT NONE
+      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
       WRITE(*,*) 'trc_freons: You should not have seen this print! error?', kt
    END SUBROUTINE trc_sms_c14b
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/CFC/trcini_cfc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/CFC/trcini_cfc.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/CFC/trcini_cfc.F90	(revision 10302)
@@ -44,5 +44,5 @@
       !! ** Method  : - Read the namcfc namelist and check the parameter values
       !!----------------------------------------------------------------------
-      INTEGER  ::  ji, jj, jn, jl, jm, js, io, ierr
+      INTEGER  ::  ji, jj, jn, jl, jm, js, io, iostatus, ierr
       INTEGER  ::  iskip = 7   ! number of 1st descriptor lines
       REAL(wp) ::  zyy, zyd
@@ -57,13 +57,19 @@
       
       CALL ctl_opn( inum, clname, 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
-      REWIND(inum)
       
       ! compute the number of year in the file
       ! file starts in 1931 do jn represent the year in the century
+      iostatus = 0
       jn = 31 
-      DO 
-        READ(inum,'(1x)',END=100) 
+      DO WHILE ( iostatus == 0 )
+        READ(inum,'(1x)', IOSTAT=iostatus, END=100)
         jn = jn + 1
-      END DO
+      ENDDO
+      IF( iostatus .NE. 0 ) THEN
+        !! Error while reading CFC input file 
+        CALL ctl_stop('trc_ini_cfc:  &
+                      & Error on the 1st reading of cfc1112sf6.atm')
+        RETURN
+      ENDIF
  100  jpyear = jn - 1 - iskip
       IF ( lwp) WRITE(numout,*) '    ', jpyear ,' years read'
@@ -72,5 +78,6 @@
       ALLOCATE( p_cfc(jpyear,jphem,3), STAT=ierr )
       IF( ierr > 0 ) THEN
-         CALL ctl_stop( 'trc_ini_cfc: unable to allocate p_cfc array' )   ;   RETURN
+         CALL ctl_stop( 'trc_ini_cfc: unable to allocate p_cfc array' )   
+         RETURN
       ENDIF
       IF( trc_sms_cfc_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'trc_ini_cfc: unable to allocate CFC arrays' )
@@ -101,11 +108,15 @@
       ! file starts in 1931 do jn represent the year in the century.jhh
       ! Read file till the end
-      jn = 31
-      DO 
+      DO jn = 31, jpyear 
+        !!READ(inum, '(F6.1,6F7.2)', IOSTAT=io) zyy, p_cfc(jn,1,1), p_cfc(jn,1,2), &
         READ(inum,*, IOSTAT=io) zyy, p_cfc(jn,1,1), p_cfc(jn,1,2), &
              & p_cfc(jn,1,3), p_cfc(jn,2,1),  &
              & p_cfc(jn,2,2), p_cfc(jn,2,3)
-        IF( io < 0 ) exit
-        jn = jn + 1
+        IF( io .NE.0 )  THEN
+          !! Error while reading CFC input file 
+          CALL ctl_stop('trc_ini_cfc:   &
+                        & Error on the 2nd reading of cfc1112sf6.atm')
+          RETURN
+        ENDIF
       END DO
 
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/CFC/trcsms_cfc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/CFC/trcsms_cfc.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/CFC/trcsms_cfc.F90	(revision 10302)
@@ -24,4 +24,5 @@
    USE trdtrc
    USE iom           ! I/O library
+   USE wrk_nemo
 
    IMPLICIT NONE
@@ -54,4 +55,7 @@
    REAL(wp) ::   xconv3 = 1.0e+3       ! conversion from mol/l/atm to mol/m3/atm
    REAL(wp) ::   xconv4 = 1.0e-12      ! conversion from mol/m3/atm to mol/m3/pptv 
+
+   !! trend temporary array:
+   REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrcfc
 
    !! * Substitutions
@@ -265,7 +269,10 @@
       !
       IF( l_trdtrc ) THEN
+          CALL wrk_alloc( jpi, jpj, jpk, ztrcfc )
           DO jn = jp_cfc0, jp_cfc1
-            CALL trd_trc( tra(:,:,:,jn), jn, jptra_sms, kt )   ! save trends
+             ztrcfc(:,:,:) = tra(:,:,:,jn)
+            CALL trd_trc( ztrcfc, jn, jptra_sms, kt )   ! save trends
           END DO
+          CALL wrk_dealloc( jpi, jpj, jpk, ztrcfc )
       END IF
       !
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 10302)
@@ -6,4 +6,5 @@
    !! History :
    !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!   -   ! 2017-08 (A. Yool)            Add air-sea flux kill switch
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -61,19 +62,21 @@
 # endif
                                    zchd, zchn, zdin, zsil
-      USE dom_oce,           ONLY: e3t_0, e3t_n, gphit, tmask
-# if defined key_iomput
+      USE dom_oce,           ONLY: e3t_0, gphit, tmask, mig, mjg
+# if defined key_vvl
+      USE dom_oce,           ONLY: e3t_n
+# endif
       USE iom,               ONLY: lk_iomput
-# endif
       USE in_out_manager,    ONLY: lwp, numout
-      USE oce,               ONLY: PCO2a_in_cpl
       USE par_kind,          ONLY: wp
       USE par_oce,           ONLY: jpi, jpim1, jpj, jpjm1
-      USE sbc_oce,           ONLY: fr_i, lk_oasis, qsr, wndm
+      USE sbc_oce,           ONLY: fr_i, qsr, wndm
       USE sms_medusa,        ONLY: jdms, jdms_input, jdms_model,          &
                                    jriver_alk, jriver_c,                  &
                                    jriver_n, jriver_si,                   &
+                                   ln_foam_medusa,                        &
                                    riv_alk, riv_c, riv_n, riv_si,         &
                                    zn_dms_chd, zn_dms_chn, zn_dms_din,    &
                                    zn_dms_mld, zn_dms_qsr,                &
+                                   f2_pco2w, f2_fco2w,                    &
                                    xnln, xnld 
       USE trc,               ONLY: med_diag
@@ -86,8 +89,11 @@
 #  else
       USE trcco2_medusa,     ONLY: trc_co2_medusa
+      USE mocsy_mainmod,     ONLY: p2fCO2
 #  endif
       USE trcdms_medusa,     ONLY: trc_dms_medusa
       USE trcoxy_medusa,     ONLY: trc_oxy_medusa
 # endif
+      USE lib_mpp,           ONLY: ctl_stop
+      USE trcstat,           ONLY: trc_rst_dia_stat 
 
    !!* Substitution
@@ -121,4 +127,19 @@
 
 # if defined key_roam
+      !! init
+      f_fco2w(:,:)       = 0.0
+      f_fco2atm(:,:)     = 0.0
+      f_schmidtco2(:,:)  = 0.0
+      f_kwco2(:,:)       = 0.0
+      f_co2starair(:,:)  = 0.0
+      f_dpco2(:,:)       = 0.0
+      f_rhosw(:,:)       = 0.0
+      f_K0(:,:)          = 0.0
+      !! air pressure (atm); ultimately this will use air 
+      !! pressure at the base of the UKESM1 atmosphere 
+      !!                                     
+      f_pp0(:,:)   = 1.0
+
+
       !!-----------------------------------------------------------
       !! Air-sea gas exchange
@@ -133,9 +154,5 @@
          DO ji = 2,jpim1
             !! OPEN wet point IF..THEN loop
-            if (tmask(ji,jj,1) == 1) then
-               IF (lk_oasis) THEN
-                  !! use 2D atm xCO2 from atm coupling
-                  f_xco2a(ji,jj) = PCO2a_in_cpl(ji,jj)
-               ENDIF
+            IF (tmask(ji,jj,1) == 1) then
                !!
                !! AXY (23/06/15): as part of an effort to update the 
@@ -161,19 +178,49 @@
                'air-sea: carb-chem kt = ', kt
                CALL flush(numout)
+               !! JPALM add carb print:
+               call trc_rst_dia_stat(f_xco2a(:,:), 'f_xco2a')
+               call trc_rst_dia_stat(wndm(:,:), 'wndm')
+               call trc_rst_dia_stat(f_kw660(:,:), 'f_kw660')
+               call trc_rst_dia_stat(ztmp(:,:), 'ztmp')
+               call trc_rst_dia_stat(zsal(:,:), 'zsal')
+               call trc_rst_dia_stat(zalk(:,:), 'zalk')
+               call trc_rst_dia_stat(zdic(:,:), 'zdic')
+               call trc_rst_dia_stat(zsil(:,:), 'zsil')
+               call trc_rst_dia_stat(zpho(:,:), 'zpho')
 #   endif
+#  if defined key_axy_carbchem
+#   if defined key_mocsy
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
             if (tmask(ji,jj,1) == 1) then
-               !! air pressure (atm); ultimately this will use air 
-               !! pressure at the base of the UKESM1 atmosphere 
-               !!                                     
-               f_pp0(ji,jj)   = 1.0
-               !!
-               !! IF(lwp) WRITE(numout,*) ' MEDUSA ztmp    =', ztmp(ji,jj)
-               !! IF(lwp) WRITE(numout,*) ' MEDUSA wndm    =', wndm(ji,jj)
-               !! IF(lwp) WRITE(numout,*) ' MEDUSA fr_i    =', fr_i(ji,jj)
-               !!
-#  if defined key_axy_carbchem
-#   if defined key_mocsy
+               !!
+               !! Jpalm -- 12-09-2017 -- add extra check after reccurent
+               !!          carbonate failure in the coupled run.
+               !!          must be associated to air-sea flux or air xCO2...
+               !!          Check MOCSY inputs
+               IF ( (zsal(ji,jj) > 75.0 ).OR.(zsal(ji,jj) < 0.0 ) .OR.        &
+                    (ztmp(ji,jj) > 50.0 ).OR.(ztmp(ji,jj) < -20.0 ) .OR.      &
+                    (zalk(ji,jj) > 35.0E2 ).OR.(zalk(ji,jj) <= 0.0 ) .OR.     &
+                    (zdic(ji,jj) > 35.0E2 ).OR.(zdic(ji,jj) <= 0.0 ) .OR.     &
+                    (f_kw660(ji,jj) > 1.0E-2 ).OR.(f_kw660(ji,jj) < 0.0 ) ) THEN
+                  IF(lwp) THEN 
+                      WRITE(numout,*) ' surface T = ',ztmp(ji,jj)
+                      WRITE(numout,*) ' surface S = ',zsal(ji,jj)
+                      WRITE(numout,*) ' surface ALK = ',zalk(ji,jj)
+                      WRITE(numout,*) ' surface DIC = ',zdic(ji,jj)
+                      WRITE(numout,*) ' KW660 = ',f_kw660(ji,jj)
+                      WRITE(numout,*) ' atm xCO2 = ',f_xco2a(ji,jj)   
+                      WRITE(numout,*) ' surface pco2w  = ',f_pco2w(ji,jj)
+                      WRITE(numout,*) ' surface fco2w  = ',f_fco2w(ji,jj)
+                      WRITE(numout,*) ' surface fco2a  = ',f_fco2atm(ji,jj)
+                      WRITE(numout,*) ' surface co2flx = ',f_co2flux(ji,jj)
+                      WRITE(numout,*) ' surface dpco2  = ',f_dpco2(ji,jj)
+                      WRITE(numout,*) ' MOCSY input: ji =', mig(ji),' jj = ', mjg(jj),  &
+                                       ' kt = ', kt 
+                      WRITE(numout,*) 'MEDUSA - Air-Sea INPUT: unrealistic surface Carb. Chemistry'
+                  ENDIF     
+                  CALL ctl_stop( 'MEDUSA - Air-Sea INPUT: ',             &
+                                 'unrealistic surface Carb. Chemistry -- INPUTS' )
+               ENDIF     
                !!
                !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate
@@ -200,8 +247,57 @@
                f_TALK(ji,jj) = (zalk(ji,jj) / f_rhosw(ji,jj)) * 1000.
                f_dcf(ji,jj)  = f_rhosw(ji,jj)
+               !! Jpalm -- 12-09-2017 -- add extra check after reccurent
+               !!          carbonate failure in the coupled run.
+               !!          must be associated to air-sea flux or air xCO2...
+               !!          Check MOCSY outputs
+               !!===================
+               !! Jpalm -- 19-02-2018 -- remove the cap - only check MOCSY inputs
+               !!       because of specific area in arabic sea where strangely
+               !!       with core 2 forcing, ALK is lower than DIC and result in 
+               !!       Enormous dpco2 - even if all carb chem caract are OK.
+               !!       and this check stops the model.
+               !!       --Input checks are already more than enough to stop the
+               !!       model if carb chem goes crazy. 
+               !!       we remove the mocsy output checks
+               !!===================
+               !!IF ( (f_pco2w(ji,jj) > 1.E4 ).OR.(f_pco2w(ji,jj) < 0.0 ) .OR.     &
+               !!    (f_fco2w(ji,jj) > 1.E4 ).OR.(f_fco2w(ji,jj) < 0.0 ) .OR.     &   
+               !!    (f_fco2atm(ji,jj) > 1.E4 ).OR.(f_fco2atm(ji,jj) < 0.0 ) .OR.     &
+               !!    (f_co2flux(ji,jj) > 1.E-1 ).OR.(f_co2flux(ji,jj) < -1.E-1 ) .OR.     &
+               !!    (f_dpco2(ji,jj) > 1.E4 ).OR.(f_dpco2(ji,jj) < -1.E4 ) ) THEN
+               !!  IF(lwp) THEN 
+               !!      WRITE(numout,*) ' surface T = ',ztmp(ji,jj)
+               !!      WRITE(numout,*) ' surface S = ',zsal(ji,jj)
+               !!      WRITE(numout,*) ' surface ALK = ',zalk(ji,jj)
+               !!      WRITE(numout,*) ' surface DIC = ',zdic(ji,jj)
+               !!      WRITE(numout,*) ' KW660 = ',f_kw660(ji,jj)
+               !!      WRITE(numout,*) ' atm xCO2 = ',f_xco2a(ji,jj)   
+               !!      WRITE(numout,*) ' surface pco2w  = ',f_pco2w(ji,jj)
+               !!      WRITE(numout,*) ' surface fco2w  = ',f_fco2w(ji,jj)
+               !!      WRITE(numout,*) ' surface fco2a  = ',f_fco2atm(ji,jj)
+               !!      WRITE(numout,*) ' surface co2flx = ',f_co2flux(ji,jj)
+               !!      WRITE(numout,*) ' surface dpco2  = ',f_dpco2(ji,jj)
+               !!      WRITE(numout,*) ' MOCSY output: ji =', mig(ji),' jj = ', mjg(jj),  &
+               !!                        ' kt = ', kt     
+               !!      WRITE(numout,*) 'MEDUSA - Air-Sea OUTPUT: unrealistic surface Carb. Chemistry'
+               !!  ENDIF     
+               !!  CALL ctl_stop( 'MEDUSA - Air-Sea OUTPUT: ',            &
+               !!                 'unrealistic surface Carb. Chemistry -- OUTPUTS' )
+               !!ENDIF     
             ENDIF
          ENDDO
       ENDDO
 
+#   if defined key_debug_medusa
+               !! JPALM add carb print:
+               call trc_rst_dia_stat(f_pco2w(:,:), 'f_pco2w')
+               call trc_rst_dia_stat(f_fco2w(:,:), 'f_fco2w')
+               call trc_rst_dia_stat(f_fco2atm(:,:), 'f_fco2atm')
+               call trc_rst_dia_stat(f_schmidtco2(:,:), 'f_schmidtco2')
+               call trc_rst_dia_stat(f_kwco2(:,:), 'f_kwco2')
+               call trc_rst_dia_stat(f_co2starair(:,:), 'f_co2starair')
+               call trc_rst_dia_stat(f_co2flux(:,:), 'f_co2flux')
+               call trc_rst_dia_stat(f_dpco2(:,:), 'f_dpco2')
+#   endif
 #   else   
 
@@ -234,4 +330,8 @@
                      iters, ' AT (', ji, ', ', jj, ', 1) AT ', kt
                endif
+               IF ( ln_foam_medusa ) THEN
+                  !! DAF (Aug 2017): calculate fCO2 for observation operator
+                  CALL p2fCO2( f_pco2w, ztmp, f_pp0, 0.0, 1, f_fco2w )
+               ENDIF
             ENDIF
          ENDDO
@@ -277,4 +377,9 @@
          ENDDO
       ENDDO
+#  endif
+
+#  if defined key_axy_killco2flux
+      !! AXY (18/08/17): single kill switch on air-sea CO2 flux for budget checking
+      f_co2flux(:,:) = 0.
 #  endif
 
@@ -408,4 +513,9 @@
                               CO2flux_conv
                !! ENDIF
+               IF ( ln_foam_medusa ) THEN
+                  !! DAF (Aug 2017): Save pCO2 and fCO2 for observation operator
+                  f2_pco2w(ji,jj) = f_pco2w(ji,jj)
+                  f2_fco2w(ji,jj) = f_pco2w(ji,jj)
+               ENDIF
                IF ( lk_iomput ) THEN
                   IF( med_diag%ATM_PCO2%dgsave ) THEN
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_iomput.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_iomput.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_iomput.F90	(revision 10302)
@@ -31,5 +31,8 @@
       !!-------------------------------------------------------------------
       USE bio_medusa_mod
-      USE dom_oce,           ONLY: e3t_0, e3t_n, mbathy, tmask
+      USE dom_oce,           ONLY: e3t_0, mbathy, tmask
+# if defined key_vvl
+      USE dom_oce,           ONLY: e3t_n
+# endif
       USE in_out_manager,    ONLY: lwp, numout
       USE par_oce,           ONLY: jpim1, jpjm1
@@ -718,4 +721,5 @@
 CONTAINS
    SUBROUTINE bio_med_diag_iomput( )                    ! Empty routine
+      IMPLICIT NONE
       WRITE(*,*) 'bio_med_diag_iomput: You should not have seen this print! error?'
    END SUBROUTINE bio_med_diag_iomput
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag.F90	(revision 10302)
@@ -32,10 +32,10 @@
       USE bio_med_diag_iomput_mod,  ONLY: bio_med_diag_iomput
       USE bio_medusa_mod
-      USE dom_oce,                  ONLY: e3t_0, e3t_n,                  &
-                                          gdepw_0, gdepw_n, tmask
+      USE dom_oce,                  ONLY: e3t_0, gdepw_0, tmask
+# if defined key_vvl
+      USE dom_oce,                  ONLY: e3t_n, gdepw_n
+# endif
       USE in_out_manager,           ONLY: lwp, numout
-# if defined key_iomput
       USE iom,                      ONLY: lk_iomput
-# endif
       USE par_oce,                  ONLY: jpim1, jpjm1
       USE sms_medusa,               ONLY: xrfn, xthetapd, xthetapn,      &
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag_slice.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag_slice.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag_slice.F90	(revision 10302)
@@ -35,10 +35,7 @@
       USE dom_oce,           ONLY: tmask
       USE in_out_manager,    ONLY: lwp, numout
-# if defined key_iomput
       USE iom,               ONLY: iom_put
-# endif
       USE lbclnk,            ONLY: lbc_lnk
-      USE trc,               ONLY: trn
-      USE oce,               ONLY: CO2Flux_out_cpl, DMS_out_cpl, chloro_out_cpl
+      USE oce,               ONLY: CO2Flux_out_cpl, DMS_out_cpl
       USE par_oce,           ONLY: jpi, jpj
       USE sbc_oce,           ONLY: lk_oasis, qsr, wndm
@@ -48,5 +45,5 @@
                                    jdms, ocal_ccd, xpar, xze,              &
                                    zb_co2_flx, zb_dms_srf,                 &
-                                   zn_co2_flx, zn_dms_srf, zn_chl_srf
+                                   zn_co2_flx, zn_dms_srf
       USE trc,               ONLY: med_diag
 
@@ -66,11 +63,4 @@
       !!
       IF (jk.eq.1) THEN
-         !! JPALM -- 02-06-2017 --
-         !! add Chl surf coupling
-         !! no need to output, just pass to cpl var
-         IF (lk_oasis) THEN
-            zn_chl_srf(:,:) = (trn(:,:,1,jpchd) + trn(:,:,1,jpchn)) * 1.0E-6  !! surf Chl in Kg-chl/m3 as needed for cpl
-            chloro_out_cpl(:,:) = zn_chl_srf(:,:)        !! Coupling Chl
-         END IF
          IF( med_diag%MED_QSR%dgsave ) THEN
             CALL iom_put( "MED_QSR"  , qsr ) !
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 10302)
@@ -6,4 +6,5 @@
    !! History :
    !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!   -   ! 2017-08 (A. Yool)            Amend bethic reservoir updating
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -32,17 +33,19 @@
       !!----------------------------------------------------------------------
       USE bio_medusa_mod
-      USE dom_oce,           ONLY: atfp, atfp1, neuler, rdt, e3t_n, tmask
+      USE dom_oce,           ONLY: atfp, atfp1, neuler, rdt, tmask
       USE in_out_manager,    ONLY: lwp, numout
-# if defined key_iomput
-      USE iom,               ONLY: iom_put, lk_iomput
-# endif
+      USE iom,               ONLY: iom_put
       USE lbclnk,            ONLY: lbc_lnk
+      USE oce,               ONLY: chloro_out_cpl 
       USE par_medusa,        ONLY: jp_medusa_2d, jp_medusa_3d,          &
-                                   jp_medusa_trd
+                                   jp_medusa_trd, jpchd, jpchn
       USE par_oce,           ONLY: jpi, jpim1, jpj, jpjm1, jpk
       USE phycst,            ONLY: rsmall
+      USE sbc_oce,           ONLY: lk_oasis
       USE sms_medusa,        ONLY: jinorgben, jorgben,                  &
                                    f3_co3, f3_h2co3, f3_hco3,           &
                                    f3_omarg, f3_omcal, f3_pH,           &
+                                   ln_foam_medusa, mld_max, pgrow_avg,  &
+                                   ploss_avg, phyt_avg,                 &
                                    za_sed_c, za_sed_ca, za_sed_fe,      &
                                    za_sed_n, za_sed_si,                 &
@@ -50,7 +53,9 @@
                                    zb_sed_n, zb_sed_si,                 &
                                    zn_sed_c, zn_sed_ca, zn_sed_fe,      &
-                                   zn_sed_n, zn_sed_si
-      USE trc,               ONLY: med_diag, nittrc000 
+                                   zn_sed_n, zn_sed_si, zn_chl_srf,     &
+                                   scl_chl, chl_out
+      USE trc,               ONLY: med_diag, nittrc000, trn 
       USE trcnam_trp,        ONLY: ln_trcadv_cen2, ln_trcadv_tvd
+      USE zdfmxl,            ONLY: hmld
  
       !! time (integer timestep)
@@ -62,159 +67,73 @@
       REAL(wp) :: fq0,fq1,fq2,fq3
 
+# if defined key_roam                     
+      !!----------------------------------------------------------------------
+      !! AXY (09/08/17): fix benthic submodel
       !!----------------------------------------------------------------------
       !! Process benthic in/out fluxes
       !! These can be handled outside of the 3D calculations since the
-      !! benthic pools (and fluxes) are 2D in nature; this code is
-      !! (shamelessly) borrowed from corresponding code in the LOBSTER
-      !! model
+      !! benthic pools (and fluxes) are 2D in nature; this code was
+      !! developed with help from George Nurser (NOC); it cannot be run
+      !! in a configuration with variable time-stepping with depth
       !!----------------------------------------------------------------------
       !!
-      !! IF(lwp) WRITE(numout,*) 'AXY: rdt = ', rdt
+      !! time-step calculation
       if (jorgben.eq.1) then
-         za_sed_n(:,:)  = zn_sed_n(:,:)  +                                 & 
-                          ( f_sbenin_n(:,:)  + f_fbenin_n(:,:)  -          &
-                            f_benout_n(:,:)  ) * (rdt / 86400.)
+         za_sed_n(:,:)  = zb_sed_n(:,:)  + ((2. * (rdt / 86400.)) * &
+                          ( f_sbenin_n(:,:)  + f_fbenin_n(:,:)  - f_benout_n(:,:)  ))
+         za_sed_fe(:,:) = zb_sed_fe(:,:) + ((2. * (rdt / 86400.)) * &
+                          ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) - f_benout_fe(:,:) ))
+         za_sed_c(:,:)  = zb_sed_c(:,:)  + ((2. * (rdt / 86400.)) * &
+                          ( f_sbenin_c(:,:)  + f_fbenin_c(:,:)  - f_benout_c(:,:)  ))
+      endif
+      if (jinorgben.eq.1) then
+         za_sed_si(:,:) = zb_sed_si(:,:) + ((2. * (rdt / 86400.)) * &
+                          ( f_fbenin_si(:,:) - f_benout_si(:,:) ))
+         za_sed_ca(:,:) = zb_sed_ca(:,:) + ((2. * (rdt / 86400.)) * &
+                          ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ))
+      endif
+      !!
+      !! time-level calculation
+      if (jorgben.eq.1) then
+         zb_sed_n(:,:)  = zn_sed_n(:,:)  + (atfp * &
+                          ( za_sed_n(:,:)  - (2. * zn_sed_n(:,:))  + zb_sed_n(:,:)  ))
          zn_sed_n(:,:)  = za_sed_n(:,:)
-         !!
-         za_sed_fe(:,:) = zn_sed_fe(:,:) +                                 &
-                          ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) -          &
-                            f_benout_fe(:,:) ) * (rdt / 86400.)
+         zb_sed_fe(:,:) = zn_sed_fe(:,:) + (atfp * &
+                          ( za_sed_fe(:,:) - (2. * zn_sed_fe(:,:)) + zb_sed_fe(:,:) ))
          zn_sed_fe(:,:) = za_sed_fe(:,:)
-         !!
-         za_sed_c(:,:)  = zn_sed_c(:,:)  +                                 &
-                          ( f_sbenin_c(:,:)  + f_fbenin_c(:,:)  -          &
-                            f_benout_c(:,:)  ) * (rdt / 86400.)
+         zb_sed_c(:,:)  = zn_sed_c(:,:)  + (atfp * &
+                          ( za_sed_c(:,:)  - (2. * zn_sed_c(:,:))  + zb_sed_c(:,:)  ))
          zn_sed_c(:,:)  = za_sed_c(:,:)
       endif
       if (jinorgben.eq.1) then
-         za_sed_si(:,:) = zn_sed_si(:,:) +                                 & 
-                          ( f_fbenin_si(:,:) - f_benout_si(:,:) ) *        &
-                          (rdt / 86400.)
+         zb_sed_si(:,:) = zn_sed_si(:,:) + (atfp * &
+                          ( za_sed_si(:,:) - (2. * zn_sed_si(:,:)) + zb_sed_si(:,:) ))
          zn_sed_si(:,:) = za_sed_si(:,:)
-         !!
-         za_sed_ca(:,:) = zn_sed_ca(:,:) +                                 &
-                          ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) *        &
-                          (rdt / 86400.)
+         zb_sed_ca(:,:) = zn_sed_ca(:,:) + (atfp * &
+                          ( za_sed_ca(:,:) - (2. * zn_sed_ca(:,:)) + zb_sed_ca(:,:) ))
          zn_sed_ca(:,:) = za_sed_ca(:,:)
       endif
-      !!
-      if (ibenthic.eq.2) then
-         !! The code below (in this if ... then ... endif loop) is
-         !! effectively commented out because it does not work as 
-         !! anticipated; it can be deleted at a later date
-         if (jorgben.eq.1) then
-            za_sed_n(:,:)  = ( f_sbenin_n(:,:)  + f_fbenin_n(:,:)  -       &
-                               f_benout_n(:,:)  ) * rdt
-            za_sed_fe(:,:) = ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) -       &
-                               f_benout_fe(:,:) ) * rdt
-            za_sed_c(:,:)  = ( f_sbenin_c(:,:)  + f_fbenin_c(:,:)  -       &
-                               f_benout_c(:,:)  ) * rdt
-         endif
-         if (jinorgben.eq.1) then
-            za_sed_si(:,:) = ( f_fbenin_si(:,:) - f_benout_si(:,:) ) * rdt
-            za_sed_ca(:,:) = ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) * rdt
-         endif
+# endif      
+
+      IF ( ln_foam_medusa ) THEN
+         !!----------------------------------------------------------------------
+         !! Diagnostics required for ocean colour assimilation:
+         !! Mixed layer average phytoplankton growth, loss and concentration
+         !! Maximum mixed layer depth
+         !!----------------------------------------------------------------------
          !!
-         !! Leap-frog scheme - only in explicit case, otherwise the 
-         !! time stepping is already being done in trczdf
-         !! IF( l_trczdf_exp .AND. (ln_trcadv_cen2 .OR. ln_trcadv_tvd) ) THEN
-         !!    zfact = 2. * rdttra(jk) * FLOAT( ndttrc )
-         !!    IF( neuler == 0 .AND. kt == nittrc000 )  zfact = rdttra(jk) * 
-         !!                                             FLOAT(ndttrc)
-         !!    if (jorgben.eq.1) then
-         !!       za_sed_n(:,:)  = zb_sed_n(:,:)  + ( zfact * za_sed_n(:,:)  )
-         !!      za_sed_fe(:,:) = zb_sed_fe(:,:) + ( zfact * za_sed_fe(:,:) )
-         !!       za_sed_c(:,:)  = zb_sed_c(:,:)  + ( zfact * za_sed_c(:,:)  )
-         !!    endif
-         !!    if (jinorgben.eq.1) then
-         !!       za_sed_si(:,:) = zb_sed_si(:,:) + ( zfact * za_sed_si(:,:) )
-         !!       za_sed_ca(:,:) = zb_sed_ca(:,:) + ( zfact * za_sed_ca(:,:) )
-         !!    endif
-         !! ENDIF
-         !! 
-         !! Time filter and swap of arrays
-         IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd  ) THEN ! centred or tvd scheme
-            IF( neuler == 0 .AND. kt == nittrc000 ) THEN
-               if (jorgben.eq.1) then
-                  zb_sed_n(:,:)  = zn_sed_n(:,:)
-                  zn_sed_n(:,:)  = za_sed_n(:,:)
-                  za_sed_n(:,:)  = 0.0
-                  !!
-                  zb_sed_fe(:,:) = zn_sed_fe(:,:)
-                  zn_sed_fe(:,:) = za_sed_fe(:,:)
-                  za_sed_fe(:,:) = 0.0
-                  !!
-                  zb_sed_c(:,:)  = zn_sed_c(:,:)
-                  zn_sed_c(:,:)  = za_sed_c(:,:)
-                  za_sed_c(:,:)  = 0.0
-               endif
-               if (jinorgben.eq.1) then
-                  zb_sed_si(:,:) = zn_sed_si(:,:)
-                  zn_sed_si(:,:) = za_sed_si(:,:)
-                  za_sed_si(:,:) = 0.0
-                  !!
-                  zb_sed_ca(:,:) = zn_sed_ca(:,:)
-                  zn_sed_ca(:,:) = za_sed_ca(:,:)
-                  za_sed_ca(:,:) = 0.0
-               endif
-            ELSE
-               if (jorgben.eq.1) then
-                  zb_sed_n(:,:)  = (atfp  *                                 &
-                                    ( zb_sed_n(:,:)  + za_sed_n(:,:)  )) +  &
-                                      (atfp1 * zn_sed_n(:,:) )
-                  zn_sed_n(:,:)  = za_sed_n(:,:)
-                  za_sed_n(:,:)  = 0.0
-                  !!
-                  zb_sed_fe(:,:) = (atfp  *                                 &
-                                    ( zb_sed_fe(:,:) + za_sed_fe(:,:) )) +  &
-                                      (atfp1 * zn_sed_fe(:,:))
-                  zn_sed_fe(:,:) = za_sed_fe(:,:)
-                  za_sed_fe(:,:) = 0.0
-                  !!
-                  zb_sed_c(:,:)  = (atfp  *                                 &
-                                    ( zb_sed_c(:,:)  + za_sed_c(:,:)  )) +  &
-                                      (atfp1 * zn_sed_c(:,:) )
-                  zn_sed_c(:,:)  = za_sed_c(:,:)
-                  za_sed_c(:,:)  = 0.0
-               endif
-               if (jinorgben.eq.1) then
-                  zb_sed_si(:,:) = (atfp  *                                 &
-                                    ( zb_sed_si(:,:) + za_sed_si(:,:) )) +  &
-                                      (atfp1 * zn_sed_si(:,:))
-                  zn_sed_si(:,:) = za_sed_si(:,:)
-                  za_sed_si(:,:) = 0.0
-                  !!
-                  zb_sed_ca(:,:) = (atfp  *                                 &
-                                    ( zb_sed_ca(:,:) + za_sed_ca(:,:) )) +  &
-                                      (atfp1 * zn_sed_ca(:,:))
-                  zn_sed_ca(:,:) = za_sed_ca(:,:)
-                  za_sed_ca(:,:) = 0.0
-               endif
-            ENDIF
-         ELSE                   !  case of smolar scheme or muscl
-            if (jorgben.eq.1) then
-               zb_sed_n(:,:)  = za_sed_n(:,:)
-               zn_sed_n(:,:)  = za_sed_n(:,:)
-               za_sed_n(:,:)  = 0.0
-               !!
-               zb_sed_fe(:,:) = za_sed_fe(:,:)
-               zn_sed_fe(:,:) = za_sed_fe(:,:)
-               za_sed_fe(:,:) = 0.0
-               !!
-               zb_sed_c(:,:)  = za_sed_c(:,:)
-               zn_sed_c(:,:)  = za_sed_c(:,:)
-               za_sed_c(:,:)  = 0.0
-            endif
-            if (jinorgben.eq.1) then
-               zb_sed_si(:,:) = za_sed_si(:,:)
-               zn_sed_si(:,:) = za_sed_si(:,:)
-               za_sed_si(:,:) = 0.0
-               !!
-               zb_sed_ca(:,:) = za_sed_ca(:,:)
-               zn_sed_ca(:,:) = za_sed_ca(:,:)
-               za_sed_ca(:,:) = 0.0
-            endif
-         ENDIF
-      endif
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF ( hmld(ji,jj) .GT. 0.0 ) THEN
+                  pgrow_avg(ji,jj) = pgrow_avg(ji,jj) / hmld(ji,jj)
+                  ploss_avg(ji,jj) = ploss_avg(ji,jj) / hmld(ji,jj)
+                  phyt_avg(ji,jj)  = phyt_avg(ji,jj)  / hmld(ji,jj)
+                  IF ( hmld(ji,jj) .GT. mld_max(ji,jj) ) THEN
+                     mld_max(ji,jj) = hmld(ji,jj)
+                  ENDIF
+               ENDIF
+            END DO
+         END DO
+      ENDIF
 
 #  if defined key_debug_medusa
@@ -253,7 +172,8 @@
                   fq1 = f_sbenin_n(ji,jj) + f_fbenin_n(ji,jj)
                   fq2 = fq0 + fq1
-                  IF (lwp) write (numout,'(a,2i3,a,3f15.10)')               &
-                     'AXY N   cons: (i,j)=',ji,jj,', (flx,ben,err)=',      &
-                     fq0,fq1,fq2
+                  fq3 = f_benout_n(ji,jj)
+                  if (lwp) write (numout,'a,2i3,a,4f15,5)')                   &
+                     'AXY N   cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',      &
+                     fq0,fq1,fq2,fq3
                ENDIF
             ENDDO
@@ -266,7 +186,8 @@
                   fq1 = f_fbenin_si(ji,jj)
                   fq2 = fq0 + fq1
-                  IF (lwp) write (numout,'(a,2i3,a,3f15.10)')               &
-                     'AXY Si  cons: (i,j)=',ji,jj,', (flx,ben,err)=',     &
-                     fq0,fq1,fq2
+                  fq3 = f_benout_si(ji,jj)
+                  if (lwp) write (numout,'a,2i3,a,4f15,5)')                   &
+                     'AXY Si  cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',     &
+                     fq0,fq1,fq2,fq3
                ENDIF
             ENDDO
@@ -278,22 +199,24 @@
                   fq0 = fflx_c(ji,jj)
                   fq1 = f_sbenin_c(ji,jj) + f_fbenin_c(ji,jj) + f_fbenin_ca(ji,jj)
-                  fq2 = f_co2flux(ji,jj) * e3t_n(ji,jj,1)
+                  fq2 = f_co2flux(ji,jj) * fse3t(ji,jj,1)
                   fq3 = fq0 + fq1
-                  IF (lwp) write (numout,'(a,2i3,a,4f15.10)')               &
-                    'AXY C   cons: (i,j)=',ji,jj,', (flx,ben,asf,err)=',  &
-                    fq0,fq1,fq2,fq3
-               ENDIF
-            ENDDO
-         ENDDO   
-         !! alkalinity
-         DO jj = 2,jpjm1
-            DO ji = 2,jpim1
-               if (tmask(ji,jj,1) == 1) then
-                  fq0 = fflx_a(ji,jj)
-                  fq1 = 2.0 * f_fbenin_ca(ji,jj)
-                  fq2 = fq0 + fq1
-                  IF (lwp) write (numout,'(a,2i3,a,3f15.10)')               &
-                     'AXY alk cons: (i,j)=',ji,jj,', (flx,ben,err)=',     &
-                     fq0,fq1,fq2
+                  fq4 = f_benout_c(ji,jj) + f_benout_ca(ji,jj)
+                  if (lwp) write (numout,'a,2i3,a,5f15,5)')                   &
+                     'AXY C   cons: (i,j)=',ji,jj,', (flx,ben,asf,err,out)=', &
+                     fq0,fq1,fq2,fq3,fq4
+                ENDIF
+             ENDDO
+          ENDDO   
+          !! alkalinity
+          DO jj = 2,jpjm1
+             DO ji = 2,jpim1
+                if (tmask(ji,jj,1) == 1) then
+                   fq0 = fflx_a(ji,jj)
+                   fq1 = 2.0 * f_fbenin_ca(ji,jj)
+                   fq2 = fq0 + fq1
+                   fq3 = 2.0 * f_benout_ca(ji,jj)
+                   if (lwp) write (numout,'a,2i3,a,4f15,5)')                   &
+                      'AXY alk cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',     &
+                      fq0,fq1,fq2,fq3
                ENDIF
             ENDDO
@@ -333,4 +256,24 @@
             ENDDO
          ENDDO
+
+         !!!---------------------------------------------------------------
+         !! Calculates Chl diag for UM coupling 
+         !!!---------------------------------------------------------------
+         !! JPALM -- 02-06-2017 --
+         !! add Chl surf coupling
+         !! no need to output, just pass to cpl var
+         IF (lk_oasis) THEN
+            IF (chl_out.eq.1) THEN
+               !! export and scale surface chl
+               zn_chl_srf(:,:) = MAX( 0.0, (trn(:,:,1,jpchd) + trn(:,:,1,jpchn)) * 1.0E-6 )
+                                 !! surf Chl in Kg-chl/m3 as needed for cpl
+            ELSEIF (chl_out.eq.2) THEN
+               !! export and scale mld chl
+               zn_chl_srf(:,:) = MAX( 0.0, fchl_ml(:,:) * 1.0E-6 )
+                                 !! mld Chl in Kg-chl/m3 as needed for cpl
+            ENDIF
+            chloro_out_cpl(:,:) = zn_chl_srf(:,:) * scl_chl        !! Coupling Chl
+         END IF
+
          !!----------------------------------------------------------------
          !! Add in XML diagnostics stuff
@@ -360,4 +303,12 @@
             CALL iom_put( "OCAL_LVL"  , fccd )
          ENDIF
+         IF ( med_diag%CHL_MLD%dgsave ) THEN
+            CALL iom_put( "CHL_MLD"  , fchl_ml )
+         ENDIF
+         IF (lk_oasis) THEN
+            IF ( med_diag%CHL_CPL%dgsave ) THEN
+               CALL iom_put( "CHL_CPL"  , chloro_out_cpl )
+            ENDIF
+         ENDIF
          IF ( med_diag%PN_JLIM%dgsave ) THEN
             CALL iom_put( "PN_JLIM"  , fjln2d )
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 10302)
@@ -6,4 +6,5 @@
    !! History :
    !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!   -   ! 2017-08 (A. Yool)            Add slow-sinking detrius variables
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -34,12 +35,10 @@
       USE bio_medusa_mod
       USE par_oce,           ONLY: jpi, jpj, jpk
-      USE sms_medusa,        ONLY: jdms
+      USE sms_medusa,        ONLY: jdms, pgrow_avg, ploss_avg, phyt_avg, mld_max
       USE trc,               ONLY: ln_diatrc, med_diag, nittrc000 
-      USE in_out_manager,    ONLY: lwp
-
-# if defined key_iomput
-      USE iom,               ONLY: lk_iomput, numout
+      USE in_out_manager,    ONLY: lwp, numout
+
+      USE iom,               ONLY: lk_iomput
       USE trcnam_medusa,     ONLY: trc_nam_iom_medusa
-# endif
 
       !! time (integer timestep)
@@ -161,7 +160,46 @@
       fprn_ml(:,:) = 0.0        !! mixed layer PP diagnostics
       fprd_ml(:,:) = 0.0        !! mixed layer PP diagnostics
+      !! AXY (16/08/17)
+      fchl_ml(:,:) = 0.0	!! mixed layer chlorophyll diagnostics
       !! 
       fslownflux(:,:) = 0.0 
       fslowcflux(:,:) = 0.0 
+      !!
+      !! JPALM -- 21-09-2017 -- needed to debug air-sea carb
+      f_xco2a(:,:)  = 0.0
+      f_pco2w(:,:)  = 0.0
+      f_ph(:,:)     = 0.0
+      f_kw660(:,:)  = 0.0
+      ztmp(:,:)  = 0.0
+      zsal(:,:)  = 0.0
+      zalk(:,:)  = 0.0
+      zdic(:,:)  = 0.0
+      zsil(:,:)  = 0.0
+# if defined key_mocsy
+      ! zpho is only defined if key_mocsy
+      ! is active, so we must protect this
+      ! initialisation accordingly. 
+      zpho(:,:)  = 0.0
+# endif
+      f_co2flux(:,:)  = 0.0 
+      f_pco2atm(:,:)  = 0.0
+      f_h2co3(:,:)    = 0.0
+      f_hco3(:,:)     = 0.0
+      f_co3(:,:)      = 0.0
+      f_omarg(:,:)    = 0.0
+      f_omcal(:,:)    = 0.0
+      !!
+      !! AXY (08/08/17): zero slow detritus fluxes
+      fslowsink(:,:)  = 0.0
+# if defined key_roam
+      fslowsinkc(:,:) = 0.0
+# endif      
+      !!
+      pgrow_avg(:,:) = 0.0
+      ploss_avg(:,:) = 0.0
+      phyt_avg(:,:)  = 0.0
+      IF( kt == nittrc000 ) THEN
+         mld_max(:,:) = 0.0
+      ENDIF
       !!
       !! allocate and initiate 2D diag
@@ -836,4 +874,5 @@
 CONTAINS
    SUBROUTINE bio_medusa_init( )                   ! Empty routine
+      IMPLICIT NONE
       WRITE(*,*) 'bio_medusa_init: You should not have seen this print! error?'
    END SUBROUTINE bio_medusa_init
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90	(revision 10302)
@@ -7,4 +7,5 @@
    !! History :
    !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!   -   ! 2017-08 (A. Yool)            Slow detritus, ML-avg chl variables
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -54,4 +55,6 @@
    !! AXY (01/03/10): add in mixed layer PP diagnostics
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fprn_ml,fprd_ml
+   !! AXY (16/08/17): add in mixed layer chlorophyll diagnostic
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fchl_ml
    !!
    !! nutrient limiting factors
@@ -94,5 +97,11 @@
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fregenfastc
 # endif
-
+   !!
+   !! AXY (08/08/17): sinking of detritus moved here
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::    fslowsink, fslowgain, fslowloss
+# if defined key_roam
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::    fslowsinkc, fslowgainc, fslowlossc
+# endif
+   !!
    !! Particle flux
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fdep1
@@ -287,5 +296,5 @@
                fjlim_pn(jpi,jpj),fjlim_pd(jpi,jpj),                   &
                fun_T(jpi,jpj),fun_Q10(jpi,jpj),                       &
-               fprn_ml(jpi,jpj),fprd_ml(jpi,jpj),                     &
+               fprn_ml(jpi,jpj),fprd_ml(jpi,jpj),fchl_ml(jpi,jpj),    &
                fnln(jpi,jpj),ffln2(jpi,jpj),                          &
                fnld(jpi,jpj),ffld(jpi,jpj),fsld(jpi,jpj),             &
@@ -316,4 +325,10 @@
                fregenfastc(jpi,jpj),                                  &
 # endif
+	       fslowsink(jpi,jpj),fslowgain(jpi,jpj),                 &
+               fslowloss(jpi,jpj),                                    &
+# if defined key_roam
+	       fslowsinkc(jpi,jpj),fslowgainc(jpi,jpj),               &
+               fslowlossc(jpi,jpj),                                   &
+# endif
                fdep1(jpi,jpj),                                        &
                ftempn(jpi,jpj),ftempsi(jpi,jpj),ftempfe(jpi,jpj),     &
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90	(revision 10302)
@@ -6,4 +6,6 @@
    !! History :
    !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!   -   ! 2017-08 (A. Yool)            Amend slow-detritus bug
+   !!   -   ! 2017-08 (A. Yool)            Reformatting for clarity
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -60,5 +62,8 @@
                                    fsil_cons, fsil_prod, fsdiss,             &
                                    ftempca, fthetad, fthetan,                &
+                                   fslowsink, fslowgain, fslowloss,          & ! AXY (22/08/17)
+                                   f_sbenin_n, f_sbenin_c,                   &
 # if defined key_roam
+                                   fslowsinkc, fslowgainc, fslowlossc,       & ! AXY (22/08/17)
                                    fcar_cons, fcar_prod, fcomm_resp,         &
                                    fddc, fflx_a, fflx_c, fflx_o2, zoxy,      &
@@ -66,9 +71,12 @@
 # endif
                                    zpds, zphd, zphn
-      USE dom_oce,           ONLY: e3t_0, e3t_n, gphit, mbathy, tmask
+      USE dom_oce,           ONLY: e3t_0, gphit, mbathy, tmask
+# if defined key_vvl
+      USE dom_oce,           ONLY: e3t_n
+# endif
       USE in_out_manager,    ONLY: lwp, numout
       USE lib_mpp,           ONLY: ctl_stop
       USE par_kind,          ONLY: wp
-      USE par_medusa,        ONLY: jp_medusa,                                &
+      USE par_medusa,        ONLY: jp_medusa, jp_msa0, jp_msa1,              &
                                    jpalk, jpchd, jpchn, jpdet, jpdic,        &
                                    jpdin, jpdtc, jpfer, jpoxy, jppds,        &
@@ -78,5 +86,5 @@
                                    jpoxy_lc, jppds_lc, jpphd_lc, jpphn_lc,   &
                                    jpsil_lc, jpzme_lc, jpzmi_lc
-      USE par_oce,           ONLY: jpi, jpim1, jpj, jpjm1
+      USE par_oce,           ONLY: jpi, jpim1, jpj, jpjm1, jpk
       USE par_trc,           ONLY: jptra
       USE sms_medusa,        ONLY: friver_dep,                               &
@@ -181,25 +189,21 @@
       ENDDO
 
-      DO jj = 2,jpjm1
-         DO ji = 2,jpim1
-            if (tmask(ji,jj,jk) == 1) then
-               !!
-               !!----------------------------------------------------------
-               !! detritus
-               btra(ji,jj,jpdet_lc) = b0 *                                      &
-                                          ! mort. losses 
-                                   (fdpn(ji,jj) + ((1.0 - xfdfrac1) *        &
-                                                   fdpd(ji,jj)) +            &
-                                    fdzmi(ji,jj) +                           &
-                                    ((1.0 - xfdfrac2) * fdzme(ji,jj)) +      &
-                                          ! assim. inefficiency
-                                    ((1.0 - xbetan) * (finmi(ji,jj) +        &
-                                                       finme(ji,jj))) -      &
-                                          ! grazing and remin.
-                                    fgmid(ji,jj) - fgmed(ji,jj) -            &
-                                    fdd(ji,jj) +                             &
-                                          ! seafloor fast->slow
-                                    ffast2slown(ji,jj))
-               !!
+      !!----------------------------------------------------------
+      !! detritus
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,jk) == 1) then
+               !!
+               btra(ji,jj,jpdet_lc) = b0 * (                           &
+                   fdpn(ji,jj)                                         & ! mort. losses 
+                 + ((1.0 - xfdfrac1) * fdpd(ji,jj))                    & ! mort. losses 
+                 + fdzmi(ji,jj)                                        & ! mort. losses
+                 + ((1.0 - xfdfrac2) * fdzme(ji,jj))                   & ! mort. losses
+                 + ((1.0 - xbetan) * (finmi(ji,jj) + finme(ji,jj)))    & ! assim. inefficiency
+                 - fgmid(ji,jj) - fgmed(ji,jj)                         & ! grazing
+                 - fdd(ji,jj)                                          & ! remin.
+                 + fslowgain(ji,jj) - fslowloss(ji,jj)                 & ! slow-sinking
+                 - (f_sbenin_n(ji,jj) / fse3t(ji,jj,jk))               & ! slow-sinking loss to seafloor
+                 + ffast2slown(ji,jj) )                                  ! seafloor fast->slow 
             ENDIF
          ENDDO
@@ -266,5 +270,5 @@
                              ! mort. loss
                          ((1.0 - xfdfrac1) * fdpds(ji,jj)) +                 &
-                             &  ! egestion of grazed Si
+                             ! egestion of grazed Si
                          ((1.0 - xfdfrac3) * fgmepds(ji,jj)) +               &
                              ! fast diss. and metab. losses
@@ -305,24 +309,26 @@
                                           ffetop(ji,jj) + ffebot(ji,jj) -    &
                                           ffescav(ji,jj) )
+            ENDIF
+         ENDDO
+      ENDDO
+
 # if defined key_roam
-               !!
-               !!----------------------------------------------------------
-               !! AXY (26/11/08): implicit detrital carbon change
-               btra(ji,jj,jpdtc_lc) = b0 * (                                    &
-                                            ! mort. losses
-                                         (xthetapn * fdpn(ji,jj)) +          &
-                                         ((1.0 - xfdfrac1) *                 &
-                                          (xthetapd * fdpd(ji,jj))) +        &
-                                         (xthetazmi * fdzmi(ji,jj)) +        &
-                                         ((1.0 - xfdfrac2) *                 &
-                                          (xthetazme * fdzme(ji,jj))) +      &
-                                             ! assim. inefficiency
-                                         ((1.0 - xbetac) *                   &
-                                          (ficmi(ji,jj) + ficme(ji,jj))) -   &
-                                             ! grazing and remin.
-                                         fgmidc(ji,jj) - fgmedc(ji,jj) -     &
-                                         fddc(ji,jj) +                       &
-                                             ! seafloor fast->slow
-                                         ffast2slowc(ji,jj) )
+      !!----------------------------------------------------------
+      !! AXY (26/11/08): implicit detrital carbon change
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,jk) == 1) then 
+               !!
+               btra(ji,jj,jpdtc_lc) = b0 * (                           &
+                   (xthetapn * fdpn(ji,jj))                            & ! mort. losses 
+                 + ((1.0 - xfdfrac1) * (xthetapd * fdpd(ji,jj)))       & ! mort. losses 
+                 + (xthetazmi * fdzmi(ji,jj))                          & ! mort. losses 
+                 + ((1.0 - xfdfrac2) * (xthetazme * fdzme(ji,jj)))     & ! mort. losses 
+                 + ((1.0 - xbetac) * (ficmi(ji,jj) + ficme(ji,jj)))    & ! assim. inefficiency
+                 - fgmidc(ji,jj) - fgmedc(ji,jj)                       & ! grazing
+                 - fddc(ji,jj)                                         & ! remin.
+                 + fslowgainc(ji,jj) - fslowlossc(ji,jj)               & ! slow-sinking
+                 - (f_sbenin_c(ji,jj) / fse3t(ji,jj,jk))               & ! slow-sinking loss to seafloor
+                 + ffast2slowc(ji,jj) )                                  ! seafloor fast->slow
             ENDIF
          ENDDO
@@ -575,8 +581,8 @@
                                                            f_o2flux(ji,jj))
                endif
+            ENDIF
+         ENDDO
+      ENDDO
 # endif
-            ENDIF
-         ENDDO
-      ENDDO
 
 # if defined key_debug_medusa
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90	(revision 10302)
@@ -40,6 +40,8 @@
 # endif
                                    zalk, zdic, zsal, zsil, ztmp 
-      USE dom_oce,           ONLY: gdept_0, gdept_n, gdepw_0, gdepw_n,    &
-                                   gphit, mbathy, tmask
+      USE dom_oce,           ONLY: gdept_0, gdepw_0, gphit, mbathy, tmask
+# if defined key_vvl
+      USE dom_oce,           ONLY: gdept_n, gdepw_n
+# endif
       USE in_out_manager,    ONLY: lwp, numout
       USE oce,               ONLY: PCO2a_in_cpl, tsb, tsn
@@ -103,8 +105,4 @@
                !! OPEN wet point IF..THEN loop
                IF (tmask(ji,jj,jk).eq.1) THEN
-                  IF (lk_oasis) THEN
-                     !! use 2D atm xCO2 from atm coupling
-                     f_xco2a(ji,jj) = PCO2a_in_cpl(ji,jj)
-                  ENDIF
                   !! Do carbonate chemistry
                   !!
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus.F90	(revision 10302)
@@ -6,4 +6,5 @@
    !! History :
    !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!   -   ! 2017-08 (A. Yool)            Revise slow-sinking of detritus
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -35,14 +36,24 @@
                                         f_sbenin_n, fdd,                   &
                                         idf, idfval,                       &   
-# if defined key_roam
+                                        fslowsink,                         &
+                                        fslowgain, fslowloss,              &
+# if defined key_roam
+                                        fslowsinkc,                        &
+                                        fslowgainc, fslowlossc,            &
                                         fddc,                              &
 # endif
                                         fun_T, fun_Q10, zdet, zdtc
       USE detritus_fast_sink_mod, ONLY: detritus_fast_sink
-      USE dom_oce,                ONLY: mbathy, tmask
+      USE dom_oce,                ONLY: mbathy, e3t_0, gphit, tmask
+# if defined key_vvl
+      USE dom_oce,                ONLY: e3t_n
+# endif
       USE in_out_manager,         ONLY: lwp, numout
       USE par_oce,                ONLY: jpim1, jpjm1
       USE sms_medusa,             ONLY: jmd, jorgben, jsfd, vsed,          &
                                         xrfn, xmd, xmdc, xthetad
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
 
       !! Level
@@ -123,4 +134,39 @@
          DO ji = 2,jpim1
             if (tmask(ji,jj,jk) == 1) then
+               !!----------------------------------------------------------------------
+               !! Detritus sinking (AXY, 08/08/18)
+	       !! Replaces slow-sinking done in trcsed_medusa.F90
+               !!
+               !! Uses the fslowsink variable to carry slow-sinking detritus from one
+               !! grid level to the next, variable fslowgain to "add" detritus sinking
+               !! from above and variable fslowloss to "subtract" detritus sinking out
+               !! to below; these variables appear in the differential equations of
+               !! detrital nitrogen and carbon below
+               !!----------------------------------------------------------------------
+               !!
+               fslowgain(ji,jj)  = fslowsink(ji,jj) / fse3t(ji,jj,jk)                  ! = mmol N / m3 / d
+               if (jk.lt.mbathy(ji,jj)) then
+                  fslowloss(ji,jj)  = (zdet(ji,jj) * vsed * 86400.) / fse3t(ji,jj,jk)  ! = mmol N / m3 / d
+               else
+                  fslowloss(ji,jj)  = 0.                                               ! = mmol N / m3 / d
+               endif
+               fslowsink(ji,jj) = fslowloss(ji,jj) * fse3t(ji,jj,jk)                   ! = mmol N / m2 / d
+               !!
+#  if defined key_roam
+               fslowgainc(ji,jj) = fslowsinkc(ji,jj) / fse3t(ji,jj,jk)                 ! = mmol C / m3 / d
+               if (jk.lt.mbathy(ji,jj)) then
+                  fslowlossc(ji,jj) = (zdtc(ji,jj) * vsed * 86400.) / fse3t(ji,jj,jk)  ! = mmol C / m3 / d
+               else
+                  fslowlossc(ji,jj) = 0.                                               ! = mmol C / m3 / d
+               endif
+               fslowsinkc(ji,jj) = fslowlossc(ji,jj) * fse3t(ji,jj,jk)                 ! = mmol C / m2 / d
+#  endif
+            ENDIF
+         ENDDO
+      ENDDO
+
+      DO jj = 2,jpjm1
+         DO ji = 2,jpim1
+            if (tmask(ji,jj,jk) == 1) then
                !!---------------------------------------------------------
                !! Detritus addition to benthos
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus_fast_sink.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus_fast_sink.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus_fast_sink.F90	(revision 10302)
@@ -66,6 +66,8 @@
                                    idf, idfval,                            &
                                    zdet, zdtc
-      USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n,         &
-                                   gphit, mbathy, tmask
+      USE dom_oce,           ONLY: e3t_0, gdepw_0, gphit, mbathy, tmask
+# if defined key_vvl
+      USE dom_oce,           ONLY: e3t_n, gdepw_n
+# endif
       USE in_out_manager,    ONLY: lwp, numout
       USE oce,               ONLY: tsn
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/iron_chem_scav.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/iron_chem_scav.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/iron_chem_scav.F90	(revision 10302)
@@ -35,6 +35,8 @@
                                    zdet, zfer, zphd, zphn, zzme, zzmi,    &
                                    idf, idfval                          
-      USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n,        &
-                                   mbathy, tmask
+      USE dom_oce,           ONLY: e3t_0, gdepw_0, mbathy, tmask
+#if defined key_vvl
+      USE dom_oce,           ONLY: e3t_n, gdepw_n
+#endif
       USE par_kind,          ONLY: wp
       USE in_out_manager,    ONLY: lwp, numout
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90	(revision 10302)
@@ -6,4 +6,5 @@
    !! History :
    !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!   -   ! 2017-08 (A. Yool)            Mean mixed layer chlorophyll
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -42,6 +43,9 @@
                                    zchd, zchn, zdet, zdin, zdtc,         &
                                    zfer, zpds, zphd, zphn, zsil,         &
-                                   zzme, zzmi
-      USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n, tmask
+                                   zzme, zzmi, fchl_ml
+      USE dom_oce,           ONLY: e3t_0, gdepw_0, tmask
+#if defined key_vvl
+      USE dom_oce,           ONLY: e3t_n, gdepw_n
+#endif
       USE in_out_manager,    ONLY: lwp, numout
       USE oce,               ONLY: tsn
@@ -55,4 +59,5 @@
                                    xvpd, xvpn, xxi
       USE zdfmxl,            ONLY: hmld
+      USE lbclnk,            ONLY: lbc_lnk
 
    !!* Substitution
@@ -373,7 +378,11 @@
                fprd_ml(ji,jj) = fprd_ml(ji,jj) + (fprd(ji,jj) * zphd(ji,jj) * &
                                                   fse3t(ji,jj,jk) * fq0)
-            ENDIF
-         ENDDO
-      ENDDO
+	       !! AXY (16/08/17)
+	       fchl_ml(ji,jj) = fchl_ml(ji,jj) + ((zchn(ji,jj) + zchd(ji,jj)) * &
+	                                          (fse3t(ji,jj,jk) * fq0) / hmld(ji,jj))
+            ENDIF
+         ENDDO
+      ENDDO
+      CALL lbc_lnk(fchl_ml(:,:),'T',1. )
 
       DO jj = 2,jpjm1
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90	(revision 10302)
@@ -36,13 +36,22 @@
                                    fdpn, fdpn2, fdzme, fdzme2,             &
                                    fdzmi, fdzmi2, fsdiss, fsin,            &
+                                   fdep1, fprn, fprd,                      &
+                                   fgmepd, fgmepn, fgmipn,                 &
                                    zphd, zphn, zpds, zzme, zzmi
-      USE dom_oce,           ONLY: tmask
+      USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n, tmask
+      USE par_kind,          ONLY: wp
       USE par_oce,           ONLY: jpim1, jpjm1
       USE phytoplankton_mod, ONLY: phytoplankton
       USE sms_medusa,        ONLY: jmpd, jmpn, jmzme, jmzmi,               &
+                                   ln_foam_medusa,                         &
+                                   pgrow_avg, ploss_avg, phyt_avg,         &
                                    xkphd, xkphn, xkzme, xkzmi,             &
                                    xmetapd, xmetapn, xmetazme, xmetazmi,   &
                                    xmpd, xmpn, xmzme, xmzmi, xsdiss
+      USE zdfmxl,            ONLY: hmld
       USE zooplankton_mod,   ONLY: zooplankton
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
 
       !! Level
@@ -50,4 +59,6 @@
 
       INTEGER :: ji, jj
+
+      REAL(wp) :: fq0
 
       !!-------------------------------------------------------------------
@@ -188,4 +199,39 @@
       ENDDO
 
+      IF ( ln_foam_medusa ) THEN
+         !! Mixed layer averages for ocean colour assimilation
+         !!
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  if (fdep1(ji,jj).le.hmld(ji,jj)) then
+                     !! this level is entirely in the mixed layer
+                     fq0 = 1.0
+                  elseif (fsdepw(ji,jj,jk).ge.hmld(ji,jj)) then
+                     !! this level is entirely below the mixed layer
+                     fq0 = 0.0
+                  else
+                     !! this level straddles the mixed layer
+                     fq0 = (hmld(ji,jj) - fsdepw(ji,jj,jk)) / fse3t(ji,jj,jk)
+                  endif
+                  !!
+                  pgrow_avg(ji,jj) = pgrow_avg(ji,jj) +                   &
+                                     (((fprn(ji,jj) * zphn(ji,jj)) +      &
+                                       (fprd(ji,jj) * zphd(ji,jj))  ) *   &
+                                      fse3t(ji,jj,jk) * fq0)
+                  ploss_avg(ji,jj) = ploss_avg(ji,jj) +                   &
+                                     ((fgmepd(ji,jj) + fdpd(ji,jj) +      &
+                                       fdpd2(ji,jj)                +      &
+                                       fgmepn(ji,jj) + fdpn(ji,jj) +      &
+                                       fdpn2(ji,jj)  + fgmipn(ji,jj) ) *  &
+                                      fse3t(ji,jj,jk) * fq0)
+                  phyt_avg(ji,jj)  = phyt_avg(ji,jj)  +                   &
+                                     ((zphn(ji,jj) + zphd(ji,jj)) *       &
+                                      fse3t(ji,jj,jk) * fq0)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
    END SUBROUTINE plankton
 
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90	(revision 10302)
@@ -177,5 +177,12 @@
    INTEGER  ::  jdms_model   !: choice of DMS model passed to atmosphere
 !!                              1 = ANDR, 2 = SIMO, 3 = ARAN, 4 = HALL
-!!
+!! JPALM --19-12-2017 -- UM people need to tune the Anderson DMS scheme    
+   REAL(wp) ::  dmsmin       !: Anderson DMS scheme - DMS minimum value
+   REAL(wp) ::  dmscut       !: Anderson DMS scheme - DMS cutoff value 
+   REAL(wp) ::  dmsslp       !: Anderson DMS scheme - DMS slope 
+!! FOR UKESM   
+   REAL(wp) ::  scl_chl      !: scaling factor for tuned Chl passed to the UM 
+   INTEGER  ::  chl_out      !: select Chl field exported and scaled for UM:
+                             !: 1- Surface Chl ; 2- MLD Chl
 !!
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   remdmp   !: depth dependent damping coefficient of passive tracers 
@@ -205,4 +212,8 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: f2_ccd_arg  !: 2D aragonite CCD depth
 !!
+!! 2D fields of pCO2 and fCO2 for observation operator
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: f2_pco2w    !: 2D pCO2
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: f2_fco2w    !: 2D fCO2
+!!
 !! 2D fields of organic and inorganic material sedimented on the seafloor
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: zb_sed_n    !: 2D organic nitrogen   (before)
@@ -276,112 +287,21 @@
 #if defined key_roam
 !!----------------------------------------------------------------------
-!! Atmospheric pCO2 data (1859 to 2100 inclusive)
-!!----------------------------------------------------------------------
-!!
-   REAL(wp), DIMENSION(242)         ::   hist_pco2 !: pCO2
-
-# if defined key_rcp26
-      !! UKMO, run AJKKH + KAAEC, RCP 2.6, pCO2 time evolution
-      DATA hist_pco2 / 286.0230, 286.1730, 286.3230, 286.4480, 286.5730, &
-      & 286.7230, 286.8480, 286.9480, 287.0480, 287.1730, &
-      & 287.3230, 287.4730, 287.6480, 287.8480, 288.0730, &
-      & 288.3480, 288.6480, 288.9730, 289.3470, 289.7470, &
-      & 290.1730, 290.6470, 291.1470, 291.6220, 292.0720, &
-      & 292.5220, 292.9220, 293.2470, 293.5220, 293.7470, &
-      & 293.9470, 294.1220, 294.2720, 294.4220, 294.5470, &
-      & 294.6470, 294.7470, 294.8470, 294.9710, 295.1710, &
-      & 295.4460, 295.7470, 296.0720, 296.4210, 296.7710, &
-      & 297.1460, 297.5710, 298.0210, 298.4460, 298.8460, &
-      & 299.2460, 299.6450, 300.0210, 300.3710, 300.7200, &
-      & 301.0450, 301.3460, 301.6710, 302.0200, 302.3450, &
-      & 302.6450, 302.9700, 303.3450, 303.7200, 304.0700, &
-      & 304.4700, 304.9200, 305.3440, 305.7700, 306.2450, &
-      & 306.7190, 307.1700, 307.6440, 308.1190, 308.5440, &
-      & 308.9440, 309.3440, 309.6940, 309.9440, 310.1190, &
-      & 310.2440, 310.3190, 310.3190, 310.2440, 310.1440, &
-      & 310.0690, 310.0440, 310.0690, 310.1440, 310.2690, &
-      & 310.4440, 310.6940, 311.0430, 311.4440, 311.8690, &
-      & 312.3680, 312.9430, 313.5430, 314.1680, 314.7900, &
-      & 315.4430, 316.2150, 317.0170, 317.7370, 318.3400, &
-      & 318.8680, 319.5900, 320.5890, 321.5470, 322.5770, &
-      & 323.8440, 324.9260, 325.7960, 327.0810, 328.6180, &
-      & 329.6830, 330.5250, 331.6880, 333.2120, 334.7870, &
-      & 336.4640, 338.2990, 339.6660, 340.7310, 342.1360, &
-      & 343.7200, 345.2200, 346.7350, 348.5820, 350.6740, &
-      & 352.4230, 353.7910, 354.9530, 355.8210, 356.7130, &
-      & 358.0630, 359.7720, 361.3970, 363.0900, 365.2560, &
-      & 367.2810, 368.7980, 370.4000, 372.4550, 374.6920, &
-      & 376.7440, 378.7440, 380.7580, 382.7080, 384.7300, &
-      & 386.9310, 389.2150, 391.4910, 393.7710, 396.0460, &
-      & 398.3240, 400.6080, 402.8950, 405.1780, 407.4550, &
-      & 409.7260, 411.9930, 414.2500, 416.4410, 418.5280, &
-      & 420.5250, 422.4390, 424.2720, 426.0200, 427.6750, &
-      & 429.2360, 430.7050, 432.0850, 433.3580, 434.5140, &
-      & 435.5740, 436.5490, 437.4420, 438.2550, 438.9810, &
-      & 439.6110, 440.1430, 440.5770, 440.9450, 441.2660, &
-      & 441.5410, 441.7840, 442.0050, 442.2040, 442.3780, &
-      & 442.5210, 442.6200, 442.6720, 442.6810, 442.6540, &
-      & 442.5830, 442.4670, 442.3270, 442.1680, 441.9960, &
-      & 441.8060, 441.5930, 441.3440, 441.0540, 440.7230, &
-      & 440.3510, 439.9300, 439.4650, 438.9730, 438.4630, &
-      & 437.9400, 437.4020, 436.8400, 436.2640, 435.6850, &
-      & 435.1030, 434.5160, 433.9170, 433.3060, 432.7010, &
-      & 432.1110, 431.5380, 430.9810, 430.4320, 429.8860, &
-      & 429.3370, 428.7810, 428.2220, 427.6490, 427.0660, &
-      & 426.4890, 425.9270, 425.3840, 424.8610, 424.3540, &
-      & 423.8540, 423.3540, 422.8530, 422.3510, 421.8410, &
-      & 421.3250, 420.8190 /
-# else
-      !! UKMO, run AJKKH + KAAEF, RCP 8.5, pCO2 time evolution
-      DATA hist_pco2 / 286.0230, 286.1730, 286.3230, 286.4480, 286.5730, &
-      & 286.7230, 286.8480, 286.9480, 287.0480, 287.1730, &
-      & 287.3230, 287.4730, 287.6480, 287.8480, 288.0730, &
-      & 288.3480, 288.6480, 288.9730, 289.3470, 289.7470, &
-      & 290.1730, 290.6470, 291.1470, 291.6220, 292.0720, &
-      & 292.5220, 292.9220, 293.2470, 293.5220, 293.7470, &
-      & 293.9470, 294.1220, 294.2720, 294.4220, 294.5470, &
-      & 294.6470, 294.7470, 294.8470, 294.9710, 295.1710, &
-      & 295.4460, 295.7470, 296.0720, 296.4210, 296.7710, &
-      & 297.1460, 297.5710, 298.0210, 298.4460, 298.8460, &
-      & 299.2460, 299.6450, 300.0210, 300.3710, 300.7200, &
-      & 301.0450, 301.3460, 301.6710, 302.0200, 302.3450, &
-      & 302.6450, 302.9700, 303.3450, 303.7200, 304.0700, &
-      & 304.4700, 304.9200, 305.3440, 305.7700, 306.2450, &
-      & 306.7190, 307.1700, 307.6440, 308.1190, 308.5440, &
-      & 308.9440, 309.3440, 309.6940, 309.9440, 310.1190, &
-      & 310.2440, 310.3190, 310.3190, 310.2440, 310.1440, &
-      & 310.0690, 310.0440, 310.0690, 310.1440, 310.2690, &
-      & 310.4440, 310.6940, 311.0430, 311.4440, 311.8690, &
-      & 312.3680, 312.9430, 313.5430, 314.1680, 314.7900, &
-      & 315.4430, 316.2150, 317.0170, 317.7370, 318.3400, &
-      & 318.8680, 319.5900, 320.5890, 321.5470, 322.5770, &
-      & 323.8440, 324.9260, 325.7960, 327.0810, 328.6180, &
-      & 329.6830, 330.5250, 331.6880, 333.2120, 334.7870, &
-      & 336.4640, 338.2990, 339.6660, 340.7310, 342.1360, &
-      & 343.7200, 345.2200, 346.7350, 348.5820, 350.6740, &
-      & 352.4230, 353.7910, 354.9530, 355.8210, 356.7130, &
-      & 358.0630, 359.7720, 361.3970, 363.0900, 365.2560, &
-      & 367.2810, 368.7980, 370.4000, 372.4550, 374.6920, &
-      & 376.7440, 378.7440, 380.7580, 382.7080, 384.7300, &
-      & 386.9420, 389.2540, 391.5670, 393.9370, 396.3920, &
-      & 398.9320, 401.5550, 404.2550, 407.0220, 409.8530, &
-      & 412.7470, 415.7050, 418.7210, 421.7880, 424.9180, &
-      & 428.1200, 431.3970, 434.7470, 438.1650, 441.6410, &
-      & 445.1700, 448.7530, 452.3920, 456.0950, 459.8810, &
-      & 463.7680, 467.7660, 471.8750, 476.0960, 480.4210, &
-      & 484.8390, 489.3470, 493.9430, 498.6400, 503.4380, &
-      & 508.3410, 513.3630, 518.5160, 523.8050, 529.2290, &
-      & 534.7780, 540.4450, 546.2230, 552.1120, 558.1110, &
-      & 564.2110, 570.4130, 576.7390, 583.1990, 589.7980, &
-      & 596.5390, 603.4110, 610.4060, 617.4940, 624.6500, &
-      & 631.8800, 639.1750, 646.5360, 653.9800, 661.5230, &
-      & 669.1840, 676.9570, 684.8290, 692.7790, 700.7690, &
-      & 708.8050, 716.8870, 725.0020, 733.1770, 741.3900, &
-      & 749.6700, 758.0480, 766.5050, 775.0350, 783.6110, &
-      & 792.2200, 800.8740, 809.5680, 818.2760, 827.0090, &
-      & 835.8020, 844.6550, 853.5730, 862.5690, 871.6190, &
-      & 880.7020, 889.8240, 898.9590, 908.1270, 917.3080, &
-      & 926.4960, 935.7040 /
-# endif
+!! JPALM -- change hist_pco2 init
+!!----------------------------------------------------------------------
+!!
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   hist_pco2 !: pCO2
+   INTEGER  :: co2_rec
+   REAL(wp) :: co2_yinit, co2_yend    !: First and Last year read in the xCO2.atm file
+   REAL(wp) :: xobs_xco2a             !: Observed atmospheric xCO2, from namelist
+#endif
+
+!!----------------------------------------------------------------------
+!! JPALM -- PI CO2 key
+!!----------------------------------------------------------------------
+!!
+#if defined key_axy_pi_co2
+   LOGICAL , PUBLIC ::   lk_pi_co2 = .TRUE.  !: PI xCO2 used
+#else
+   LOGICAL , PUBLIC ::   lk_pi_co2 = .FALSE. !: PI xCO2 unused
 #endif
 
@@ -434,4 +354,15 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: cmask       !: ???
 
+!!----------------------------------------------------------------------
+!! Parameters required for ocean colour assimilation
+!!----------------------------------------------------------------------
+!!
+   LOGICAL :: ln_foam_medusa                                 !: Flag to calculate and save diagnostics
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: pgrow_avg  !: Mixed layer average phytoplankton growth
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ploss_avg  !: Mixed layer average phytoplankton loss
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: phyt_avg   !: Mixed layer average phytoplankton
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: mld_max    !: Maximum mixed layer depth
+!!
+
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
@@ -446,5 +377,5 @@
       !!----------------------------------------------------------------------
       USE lib_mpp , ONLY: ctl_warn
-      INTEGER ::   ierr(8)        ! Local variables
+      INTEGER ::   ierr(9)        ! Local variables
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -456,4 +387,5 @@
       !* 2D and 3D fields of carbonate system parameters
       ALLOCATE( f2_ccd_cal(jpi,jpj)  , f2_ccd_arg(jpi,jpj)  ,       &
+                f2_pco2w(jpi,jpj)    , f2_fco2w(jpi,jpj)    ,       &
          &      f3_pH(jpi,jpj,jpk)   , f3_h2co3(jpi,jpj,jpk),       &
          &      f3_hco3(jpi,jpj,jpk) , f3_co3(jpi,jpj,jpk)  ,       &
@@ -504,4 +436,7 @@
          &      ffln(jpi,jpj,jpk)    , fflf(jpi,jpj,jpk)    ,       &
          &      ffls(jpi,jpj,jpk)    , cmask(jpi,jpj)       ,    STAT=ierr(8) ) 
+      !* Fields for ocean colour data assimilation
+      ALLOCATE( pgrow_avg(jpi,jpj)   , ploss_avg(jpi,jpj)   ,       &
+         &      phyt_avg(jpi,jpj)    , mld_max(jpi,jpj)     ,    STAT=ierr(9) )
 #endif
       !
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 10302)
@@ -23,4 +23,5 @@
    !!  -   !  2016-11  (A. Yool)              Updated diags for CMIP6
    !!  -   !  2017-05  (A. Yool)              Added extra DMS calculation
+   !!  -   !  2017-11  (J. Palm, A. Yool)     Diagnose tracer excursions
    !!----------------------------------------------------------------------
    !!
@@ -77,16 +78,21 @@
                                             zfer, zpds, zphd, zphn, zsil,   &
                                             zzme, zzmi
-      USE dom_oce,                    ONLY: e3t_0, e3t_n,                   &
-                                            gdept_0, gdept_n,               &
-                                            gdepw_0, gdepw_n,               &
-                                            nday_year, nsec_day, nyear,     &
-                                            rdt, tmask
+      USE dom_oce,                    ONLY: e3t_0, gdept_0,  gdepw_0,       &
+                                            nday_year, nsec_day,            & 
+                                            nyear, nyear_len, ndastp,       &
+                                            nsec_month,                     &
+                                            rdt, tmask, mig, mjg, nimpp,    &
+                                            njmpp 
+#if defined key_vvl
+      USE dom_oce,                    ONLY: e3t_n, gdept_n,  gdepw_n
+#endif
+
       USE in_out_manager,             ONLY: lwp, numout, nn_date0
-# if defined key_iomput
       USE iom,                        ONLY: lk_iomput
-# endif
       USE lbclnk,                     ONLY: lbc_lnk
-      USE lib_mpp,                    ONLY: ctl_stop
-      USE oce,                        ONLY: tsb, tsn
+      USE lib_mpp,                    ONLY: mpp_max, mpp_maxloc,            & 
+                                            mpp_min, mpp_minloc,            &
+                                            ctl_stop, ctl_warn, lk_mpp  
+      USE oce,                        ONLY: tsb, tsn, PCO2a_in_cpl
       USE par_kind,                   ONLY: wp
       USE par_medusa,                 ONLY: jpalk, jpchd, jpchn, jpdet,     &
@@ -98,5 +104,11 @@
       !! JPALM (27-06-2016): add lk_oasis for CO2 and DMS coupling with atm
       USE sbc_oce,                    ONLY: lk_oasis
-      USE sms_medusa,                 ONLY: hist_pco2
+      USE sms_medusa,                 ONLY: hist_pco2, co2_yinit, co2_yend, &
+# if defined key_roam
+                                            xobs_xco2a,                     &
+# endif
+                                            pgrow_avg,                      &
+                                            ploss_avg, phyt_avg, mld_max,   &
+                                            lk_pi_co2, ln_foam_medusa
       USE trc,                        ONLY: ln_rsttr, nittrc000, trn
       USE bio_medusa_init_mod,        ONLY: bio_medusa_init
@@ -110,4 +122,5 @@
       USE bio_medusa_diag_slice_mod,  ONLY: bio_medusa_diag_slice
       USE bio_medusa_fin_mod,         ONLY: bio_medusa_fin
+      USE trcstat,                    ONLY: trc_rst_dia_stat
 
       IMPLICIT NONE
@@ -115,4 +128,5 @@
       
       PUBLIC   trc_bio_medusa    ! called in trcsms_medusa.F90
+      PUBLIC   trc_bio_exceptionnal_fix ! here 
 
    !!* Substitution
@@ -176,5 +190,11 @@
       !! 
       !! temporary variables
-      REAL(wp) ::    fq0,fq1,fq2,fq3,fq4
+      REAL(wp) ::    fq3,fq4
+      REAL(wp) ::    this_y, this_d, this_s, fyear
+      !!
+      !! T and S check temporary variable
+      REAL(wp) ::    sumtsn, tsnavg
+      INTEGER  ::    summask
+      CHARACTER(40) :: charout, charout2, charout3, charout4, charout5
       !!
       !!------------------------------------------------------------------
@@ -283,44 +303,88 @@
       !!------------------------------------------------------------------
       !!
-      !! what's atmospheric pCO2 doing? (data start in 1859)
-      iyr1  = nyear - 1859 + 1
-      iyr2  = iyr1 + 1
-      if (iyr1 .le. 1) then
-         !! before 1860
-         f_xco2a(:,:) = hist_pco2(1)
-      elseif (iyr2 .ge. 242) then
-         !! after 2099
-         f_xco2a(:,:) = hist_pco2(242)
-      else
-         !! just right
-         fq0 = hist_pco2(iyr1)
-         fq1 = hist_pco2(iyr2)
-         fq2 = real(nsec_day) / (60.0 * 60.0 * 24.0)
-         !! AXY (14/06/12): tweaked to make more sense (and be correct)
-#  if defined key_bs_axy_yrlen
-         !! bugfix: for 360d year with HadGEM2-ES forcing
-         fq3 = (real(nday_year) - 1.0 + fq2) / 360.0  
-#  else
-         !! original use of 365 days (not accounting for leap year or 
-         !! 360d year)
-         fq3 = (real(nday_year) - 1.0 + fq2) / 365.0
-#  endif
-         fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3)
-         f_xco2a(:,:) = fq4
-      endif
-#  if defined key_axy_pi_co2
-      !! OCMIP pre-industrial pCO2
-      !! f_xco2a(:,:) = 284.725  !! CMIP5 pre-industrial pCO2
-      f_xco2a = 284.317          !! CMIP6 pre-industrial pCO2 
-#  endif
-      !! IF(lwp) WRITE(numout,*) ' MEDUSA nyear     =', nyear
-      !! IF(lwp) WRITE(numout,*) ' MEDUSA nsec_day  =', real(nsec_day)
-      !! IF(lwp) WRITE(numout,*) ' MEDUSA nday_year =', real(nday_year)
-      !! AXY (29/01/14): remove surplus diagnostics
-      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq0       =', fq0
-      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq1       =', fq1
-      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq2       =', fq2
-      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq3       =', fq3
-      IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2  =', f_xco2a(1,1)
+      IF (lk_oasis) THEN 
+         !! xCO2 from coupled
+         IF ( ( kt == nittrc000 ) .AND. lwp )     &
+             WRITE(numout,*) '** MEDUSA Atm xCO2 given by the UM **'
+         f_xco2a(:,:) = PCO2a_in_cpl(:,:)
+         !! Check the xCO2 from the UM is OK
+         !! piece of code moved from air-sea.F90
+         !!---
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               !! OPEN wet point IF..THEN loop
+               IF (tmask(ji,jj,1) == 1) then
+                  !!! Jpalm test on atm xCO2
+                  IF ( (f_xco2a(ji,jj) .GT. 10000.0 ).OR.     &
+                       (f_xco2a(ji,jj) .LE. 0.0 ) ) THEN
+                     IF(lwp) THEN
+                        WRITE(numout,*) ' atm xCO2 = ',f_xco2a(ji,jj),       &
+                                        ' -- ji =', mig(ji),' jj = ', mjg(jj)
+                     ENDIF
+                     CALL ctl_stop( 'MEDUSA - trc_bio :',     &
+                                    'unrealistic coupled atm xCO2 ' )
+                  ENDIF
+               ENDIF
+            ENDDO
+         ENDDO
+         !!--- 
+      ELSEIF (lk_pi_co2) THEN
+         !! OCMIP pre-industrial xCO2
+         IF ( ( kt == nittrc000 ) .AND. lwp )     &
+             WRITE(numout,*) '** MEDUSA Atm xCO2 fixed to pre-industrial value **'
+         !! f_xco2a(:,:) = 284.725       !! CMIP5 pre-industrial pCO2
+         f_xco2a(:,:) = 284.317          !! CMIP6 pre-industrial pCO2 
+      ELSEIF ( xobs_xco2a > 0.0 ) THEN
+         IF(lwp) WRITE(numout,*) ' using observed atm pCO2 = ', xobs_xco2a
+         f_xco2a(:,:) = xobs_xco2a
+      ELSE      
+         !! xCO2 from file
+         !! AXY - JPALM new interpolation scheme usinf nyear_len
+         this_y = real(nyear)
+         this_d = real(nday_year)
+         this_s = real(nsec_day)
+         !!
+         fyear = this_y + ((this_d - 1) + (this_s / (60. * 60. * 24.))) / real(nyear_len(1))
+         !!
+         IF ( ( kt == nittrc000 ) .AND. lwp ) THEN 
+            WRITE(numout,*) '** MEDUSA Atm xCO2 from file **'
+            WRITE(numout,*) ' MEDUSA year      =', this_y
+            WRITE(numout,*) ' Year length      =', real(nyear_len(1))
+            WRITE(numout,*) ' MEDUSA nday_year =', this_d
+            WRITE(numout,*) ' MEDUSA nsec_day  =', this_s
+         ENDIF
+         !! 
+         !! different case test 
+         IF (fyear .LE. co2_yinit) THEN
+            !! before first record -- pre-industrial value
+            f_xco2a(:,:) = hist_pco2(1)  
+         ELSEIF (fyear .GE. co2_yend) THEN
+            !! after last record - continue to use the last value
+            f_xco2a(:,:) = hist_pco2(int(co2_yend - co2_yinit + 1.) )
+         ELSE
+            !! just right
+            iyr1 = int(fyear - co2_yinit) + 1
+            iyr2 = iyr1 + 1 
+            fq3 = fyear - real(iyr1) - co2_yinit + 1.
+            fq4 = ((1 - fq3) * hist_pco2(iyr1)) + (fq3 * hist_pco2(iyr2))
+            f_xco2a(:,:) = fq4
+            !!
+            IF ( ( kt == nittrc000 ) .AND. lwp ) THEN 
+               WRITE(numout,*) ' MEDUSA year1     =', iyr1
+               WRITE(numout,*) ' MEDUSA year2     =', iyr2
+               WRITE(numout,*) ' xCO2 year1       =', hist_pco2(iyr1)
+               WRITE(numout,*) ' xCO2 year2       =', hist_pco2(iyr2)
+               WRITE(numout,*) ' Year2 weight     =', fq3
+            ENDIF
+         ENDIF
+      ENDIF
+ 
+      !! Writing xCO2 in output on start and on the 1st tsp of each  month
+      IF ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        &
+           ( nsec_month .LE. INT(rdt) ) )  THEN
+           IF ( lwp )  WRITE(numout,*) ' *** Atm xCO2 *** -- kt:', kt,       &
+                                       ';  current date:', ndastp
+          call trc_rst_dia_stat(f_xco2a(:,:), 'atm xCO2')
+      ENDIF
 # endif
 
@@ -348,6 +412,25 @@
       !!          x * 30d + 1*rdt(i.e: mod = rdt)   
       !!          ++ need to pass carb-chem output var through restarts
-      If ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        &
-           ((86400*mod(nn_date0,100) + mod(kt*rdt,2592000.)) == rdt) ) THEN
+      !!If ( (kt == nitt8rc000 .AND. .NOT.ln_rsttr) .OR.          &
+      !!     ( (mod(kt*rdt,2592000.)) == rdt) THEN
+      !!=============================
+      !! (Jpalm -- updated for restartability issues)
+      !! We want this to be start of month or if starting afresh from  
+      !! climatology - marc 20/6/17 
+      !!If ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                         & 
+      !!     ((86400*mod(nn_date0,100) + mod(kt*rdt,2592000.)) == rdt) ) THEN 
+      !!=============================
+      !! Jpalm -- 15-02-2018 -- need to change 3D carb-chem call freq again.
+      !! previous call did not work, probably the (86400*mod(nn_date0,100) part
+      !! should not be in...
+      !! now use the NEMO calendar tool : nsec_month to be sure to call 
+      !! at the beginning of a new month .
+      !! DAF: For FOAM we want to run daily
+      IF ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        &
+           ( nsec_month .LE. INT(rdt) )          .OR.                        &
+           ( nsec_day   .LE. INT(rdt) .AND. ln_foam_medusa ) )  THEN
+           IF ( lwp )  WRITE(numout,*)                                       &         
+                              ' *** 3D carb chem call *** -- kt:', kt,       &
+                              ';  current date:', ndastp 
          !!---------------------------------------------------------------
          !! Calculate the carbonate chemistry for the whole ocean on the first
@@ -450,4 +533,5 @@
 
 # if defined key_roam
+         !! extra MEDUSA-2 tracers
          DO jj = 2,jpjm1
             DO ji = 2,jpim1
@@ -456,7 +540,7 @@
                   zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc))
                   !! dissolved inorganic carbon
-                  zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic))
+                  zdic(ji,jj) = trn(ji,jj,jk,jpdic)
                   !! alkalinity
-                  zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk))
+                  zalk(ji,jj) = trn(ji,jj,jk,jpalk)
                   !! oxygen
                   zoxy(ji,jj) = max(0.,trn(ji,jj,jk,jpoxy))
@@ -470,30 +554,12 @@
                   ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem)
                   zsal(ji,jj) = tsn(ji,jj,jk,jp_sal)
-                  !!
-	          !! AXY (28/02/14): check input fields
-                  if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then
-                     IF(lwp) WRITE(numout,*)                                 &
-                        ' trc_bio_medusa: T WARNING 2D, ',                   &
-                        tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem),          &
-                        ' at (', ji, ',', jj, ',', jk, ') at time', kt
-		     IF(lwp) WRITE(numout,*)                                 &
-                        ' trc_bio_medusa: T SWITCHING 2D, ',                 &
-                        tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem)
-                     !! temperatur
-                     ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem)
-                  endif
-                  if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then
-                     IF(lwp) WRITE(numout,*)                                 &
-                        ' trc_bio_medusa: S WARNING 2D, ',                   &
-                        tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal),          &
-                        ' at (', ji, ',', jj, ',', jk, ') at time', kt
-                  endif
                ENDIF
             ENDDO
          ENDDO
 # else
+         !! diagnostic MEDUSA-1 detrital carbon tracer
          DO jj = 2,jpjm1
             DO ji = 2,jpim1
-               if (tmask(ji,jj,jk) == 1) then
+               IF (tmask(ji,jj,jk) == 1) THEN
                   !! implicit detrital carbon
                   zdtc(ji,jj) = zdet(ji,jj) * xthetad
@@ -502,4 +568,75 @@
          ENDDO
 # endif
+
+# if defined key_roam
+         !! ---------------------------------------------
+         !! JPALM -- 14-12-2017 -- Here, before any exeptionnal crazy value is
+         !!              removed, we want to tell to the Master Processor that this 
+         !!              Exceptionnal value did exist
+         !!
+         Call trc_bio_check(kt, jk)
+
+         !!================================================================
+	 !! AXY (03/11/17): check input fields
+         !! tracer values that exceed thresholds can cause carbonate system
+         !! failures when passed to MOCSY; temporary temperature excursions
+         !! in recent UKESM0.8 runs trigger such failures but are too short
+         !! to have physical consequences; this section checks for such
+         !! values and amends them using neighbouring values
+         !! 
+         !! the check on temperature here is also carried out at the end of
+         !! each model time step and anomalies are reported in the master
+         !! ocean.output file; the error reporting below is strictly local
+         !! to the relevant ocean.output_XXXX file so will not be visible
+         !! unless all processors are reporting output
+         !!================================================================
+         !!
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               if (tmask(ji,jj,jk) == 1) then
+                  !! the thresholds for the four tracers are ...
+                  IF ( (ztmp(ji,jj) .LT. -3.0) .OR. (ztmp(ji,jj) .GT. 40.0  ) .OR.   &
+                       (zsal(ji,jj) .LE.  0.0) .OR. (zsal(ji,jj) .GT. 50.0  ) .OR.   &
+                       (zdic(ji,jj) .LE.  0.0) .OR. (zdic(ji,jj) .GT. 4.0E3 ) .OR.   &
+                       (zalk(ji,jj) .LE.  0.0) .OR. (zalk(ji,jj) .GT. 4.0E3 ) ) THEN
+                     !!
+                     !! all tracer values are reported in the event of any excursion
+                     IF (lwp) THEN
+                         WRITE(charout,*)  ' Tmp = ', ztmp(ji,jj)
+                         WRITE(charout2,*) ' Sal = ', zsal(ji,jj)
+                         WRITE(charout3,*) ' DIC = ', zdic(ji,jj)
+                         WRITE(charout4,*) ' Alk = ', zalk(ji,jj)
+                         WRITE(charout5,*) mig(ji), mjg(jj), jk, kt 
+                         CALL ctl_warn( 'trc_bio_medusa: carbonate chemistry WARNING:',  &
+                            TRIM(charout),TRIM(charout2),TRIM(charout3),TRIM(charout4),  & 
+                            'at i, j, k, kt:', TRIM(charout5) )
+                     ENDIF
+                     !!
+                     !! Detect, report and correct tracer excursions
+                     IF ( (ztmp(ji,jj) .LT. -3.0) .OR. (ztmp(ji,jj) .GT.  40.0) )             &
+                        CALL trc_bio_exceptionnal_fix(                                         &
+                        tsn(ji-1:ji+1,jj-1:jj+1,jk,jp_tem), tmask(ji-1:ji+1,jj-1:jj+1,jk),    &
+                        'Tmp', -3.0, 40.0, ztmp(ji,jj) )
+                     !!
+                     IF ( (zsal(ji,jj) .LE. 0.0) .OR. (zsal(ji,jj) .GT.  50.0)  )             &
+                        CALL trc_bio_exceptionnal_fix(                                         &
+                        tsn(ji-1:ji+1,jj-1:jj+1,jk,jp_sal), tmask(ji-1:ji+1,jj-1:jj+1,jk),    &
+                        'Sal', 1.0, 50.0, zsal(ji,jj) )
+                     !!
+                     IF ( (zdic(ji,jj) .LE. 0.0) .OR. (zdic(ji,jj) .GT. 4.0E3)  )             &
+                        CALL trc_bio_exceptionnal_fix(                                         &
+                        trn(ji-1:ji+1,jj-1:jj+1,jk,jpdic), tmask(ji-1:ji+1,jj-1:jj+1,jk),     &
+                        'DIC', 100.0, 4.0E3, zdic(ji,jj) )
+                     !!
+                     IF ( (zalk(ji,jj) .LE. 0.0) .OR. (zalk(ji,jj) .GT. 4.0E3)  )             &
+                        CALL trc_bio_exceptionnal_fix(                                         &
+                        trn(ji-1:ji+1,jj-1:jj+1,jk,jpalk), tmask(ji-1:ji+1,jj-1:jj+1,jk),     &
+                        'Alk', 100.0, 4.0E3, zalk(ji,jj) )
+                  ENDIF
+               ENDIF
+            ENDDO
+         ENDDO
+# endif
+
 # if defined key_debug_medusa
          DO jj = 2,jpjm1
@@ -657,4 +794,216 @@
    END SUBROUTINE trc_bio_medusa
 
+
+
+   SUBROUTINE trc_bio_exceptionnal_fix(tiny_var, tiny_mask, var_nm, mini, maxi, varout)
+      !! JPALM (27/10/17): This function is called only when abnormal values that 
+      !! could break the model's carbonate system are fed to MEDUSA
+      !! 
+      !! The function calculates an average tracer value based on the 3x3 cell
+      !! neighbourhood around the abnormal cell, and reports this back
+      !! 
+      !! Tracer array values are not modified, but MEDUSA uses "corrected" values
+      !! in its calculations
+      !!
+      !! temporary variables
+      REAL(wp), INTENT( in ), DIMENSION(3,3) :: tiny_var, tiny_mask
+      CHARACTER(3), INTENT( in )            :: var_nm
+      REAL(wp), INTENT( in )                 :: mini, maxi
+      REAL(wp), INTENT( out )                :: varout
+      REAL(wp)      :: sumtsn, tsnavg
+      INTEGER       :: summask
+      CHARACTER(25) :: charout1, charout2
+      CHARACTER(60) :: charout3, charout4
+      INTEGER       :: ii,ij
+    
+      !! point to the center of the 3*3 zoom-grid, to check around
+      ii = 2
+      ij = 2
+      !! Print surounding values to check if isolated Crazy value or 
+      !! General error
+      IF(lwp) THEN 
+          WRITE(numout,*)                                 &
+            '----------------------------------------------------------------------'
+          WRITE(numout,*)                                 &
+            'trc_bio_medusa: 3x3 neighbourhood surrounding abnormal ', TRIM(var_nm)
+          WRITE(numout,9100)                              &
+            3, tiny_var(ii-1,ij+1), tiny_var(ii  ,ij+1), tiny_var(ii+1,ij+1)
+          WRITE(numout,9100)                              &
+            2, tiny_var(ii-1,ij  ), tiny_var(ii  ,ij  ), tiny_var(ii+1,ij  )
+          WRITE(numout,9100)                              &
+            1, tiny_var(ii-1,ij-1), tiny_var(ii  ,ij-1), tiny_var(ii+1,ij-1)
+          WRITE(numout,*)                                 &
+            'trc_bio_medusa: 3x3 land-sea neighbourhood, tmask'
+          WRITE(numout,9100)                              &
+            3, tiny_mask(ii-1,ij+1), tiny_mask(ii  ,ij+1), tiny_mask(ii+1,ij+1)
+          WRITE(numout,9100)                              &
+            2, tiny_mask(ii-1,ij  ), tiny_mask(ii  ,ij  ), tiny_mask(ii+1,ij  )
+          WRITE(numout,9100)                              &
+            1, tiny_mask(ii-1,ij-1), tiny_mask(ii  ,ij-1), tiny_mask(ii+1,ij-1)
+      ENDIF
+      !! Correct out of range values
+      sumtsn = ( tiny_mask(ii-1,ij+1) * tiny_var(ii-1,ij+1) ) +  &
+               ( tiny_mask(ii  ,ij+1) * tiny_var(ii  ,ij+1) ) +  &
+               ( tiny_mask(ii+1,ij+1) * tiny_var(ii+1,ij+1) ) +  &
+               ( tiny_mask(ii-1,ij  ) * tiny_var(ii-1,ij  ) ) +  &
+               ( tiny_mask(ii+1,ij  ) * tiny_var(ii+1,ij  ) ) +  &
+               ( tiny_mask(ii-1,ij-1) * tiny_var(ii-1,ij-1) ) +  &
+               ( tiny_mask(ii  ,ij-1) * tiny_var(ii  ,ij-1) ) +  &
+               ( tiny_mask(ii+1,ij-1) * tiny_var(ii+1,ij-1) )
+      !!
+      summask = tiny_mask(ii-1,ij+1) + tiny_mask(ii  ,ij+1)   +  &
+                tiny_mask(ii+1,ij+1) + tiny_mask(ii-1,ij  )   +  &
+                tiny_mask(ii+1,ij  ) + tiny_mask(ii-1,ij-1)   +  &
+                tiny_mask(ii  ,ij-1) + tiny_mask(ii+1,ij-1)
+      !!
+      IF ( summask .GT. 0 ) THEN
+         tsnavg = ( sumtsn / summask )
+         varout = MAX( MIN( maxi, tsnavg), mini )
+      ELSE   
+         IF (ztmp(ii,ij) .LT. mini )  varout = mini
+         IF (ztmp(ii,ij) .GT. maxi )  varout = maxi
+      ENDIF
+      !!
+      IF (lwp) THEN 
+          WRITE(charout1,9200) tiny_var(ii,ij)
+          WRITE(charout2,9200) varout
+          WRITE(charout3,*) ' ', charout1, ' -> ', charout2
+          WRITE(charout4,*) ' Tracer: ', trim(var_nm)
+      !!
+          WRITE(numout,*) 'trc_bio_medusa: ** EXCEPTIONAL VALUE SWITCHING **'
+          WRITE(numout,*) charout4 
+          WRITE(numout,*) charout3
+          WRITE(numout,*) '----------------------------------------------------------------------'
+          WRITE(numout,*) ' '
+      ENDIF
+
+9100  FORMAT('Row:', i1, '  ', e16.6, ' ', e16.6, ' ', e16.6)
+9200  FORMAT(e16.6)
+
+   END SUBROUTINE trc_bio_exceptionnal_fix 
+
+   SUBROUTINE trc_bio_check(kt, jk)
+      !!-----------------------------------
+      !! JPALM -- 14-12-2017 -- Still dealing with this micro-boil/carb failure
+      !!                     problem. The model is now able to correct a local
+      !!                     crazy value. but does it silently.
+      !!                     We need to spread the word to the master processor. we
+      !!                     don't want the model  to correct values without telling us
+      !!                     This module will tell at least when crazy DIC or
+      !!                     ALK appears.
+      !!-------------------------------------
+      REAL(wp)              :: zmax, zmin    ! temporary scalars
+      INTEGER               :: ji,jj         ! dummy loop
+      INTEGER               :: ii,ij         ! temporary scalars 
+      INTEGER, DIMENSION(2) :: ilocs         ! 
+      INTEGER, INTENT( in ) :: kt, jk
+      !!
+      !!==========================
+      !! DIC Check
+      zmax =  -5.0  ! arbitrary  low maximum value
+      zmin =  4.0E4  ! arbitrary high minimum value
+      DO jj = 2, jpjm1
+         DO ji = 2,jpim1
+            IF( tmask(ji,jj,1) == 1) THEN
+               zmax = MAX(zmax,zdic(ji,jj))     ! find local maximum
+               zmin = MIN(zmin,zdic(ji,jj))     ! find local minimum
+            ENDIF
+         END DO
+      END DO
+      IF( lk_mpp )   CALL mpp_max( zmax )       ! max over the global domain
+      IF( lk_mpp )   CALL mpp_min( zmin )       ! min over the global domain
+      !
+      IF( zmax .GT. 4.0E3) THEN  ! we've got a problem
+         IF (lk_mpp) THEN
+            CALL mpp_maxloc ( zdic(:,:),tmask(:,:,1), zmax, ii,ij )
+         ELSE
+            ilocs = MAXLOC( zdic(:,:), mask = tmask(:,:,1) == 1. )
+            ii = ilocs(1) + nimpp - 1
+            ij = ilocs(2) + njmpp - 1
+         ENDIF
+         !
+         IF(lwp) THEN
+            WRITE(numout,*) 'trc_bio:tracer anomaly: *****    WARNING     *****'
+            WRITE(numout,*) 'trc_bio:tracer anomaly: DIC concentration > 4000 '
+            WRITE(numout,9600) kt, zmax, ii, ij, jk
+            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****'
+         ENDIF
+      ENDIF
+      !
+      IF( zmin .LE. 0.0) THEN  ! we've got a problem
+         IF (lk_mpp) THEN
+            CALL mpp_minloc ( zdic(:,:),tmask(:,:,1), zmin, ii,ij )
+         ELSE
+            ilocs = MINLOC( zdic(:,:), mask = tmask(:,:,1) == 1. )
+            ii = ilocs(1) + nimpp - 1
+            ij = ilocs(2) + njmpp - 1
+         ENDIF
+         !
+         IF(lwp) THEN
+            WRITE(numout,*) 'trc_bio:tracer anomaly: *****    WARNING     *****'
+            WRITE(numout,*) 'trc_bio:tracer anomaly: DIC concentration <= 0 '
+            WRITE(numout,9700) kt, zmin, ii, ij, jk
+            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****'
+         ENDIF
+      ENDIF
+      !!
+      !!==========================
+      !! ALKALINITY Check
+      zmax =  -5.0  ! arbitrary  low maximum value
+      zmin =  4.0E4  ! arbitrary high minimum value
+      DO jj = 2, jpjm1
+         DO ji = 2,jpim1
+            IF( tmask(ji,jj,1) == 1) THEN
+               zmax = MAX(zmax,zalk(ji,jj))     ! find local maximum
+               zmin = MIN(zmin,zalk(ji,jj))     ! find local minimum
+            ENDIF
+         END DO
+      END DO
+      IF( lk_mpp )   CALL mpp_max( zmax )       ! max over the global domain
+      IF( lk_mpp )   CALL mpp_min( zmin )       ! min over the global domain
+      !
+      IF( zmax .GT. 4.0E3) THEN  ! we've got a problem
+         IF (lk_mpp) THEN
+            CALL mpp_maxloc ( zalk(:,:),tmask(:,:,1), zmax, ii,ij )
+         ELSE
+            ilocs = MAXLOC( zalk(:,:), mask = tmask(:,:,1) == 1. )
+            ii = ilocs(1) + nimpp - 1
+            ij = ilocs(2) + njmpp - 1
+         ENDIF
+         !
+         IF(lwp) THEN
+            WRITE(numout,*) 'trc_bio:tracer anomaly: *****     WARNING    *****'
+            WRITE(numout,*) 'trc_bio:tracer anomaly: ALK concentration > 4000 '
+            WRITE(numout,9800) kt, zmax, ii, ij, jk
+            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****'
+         ENDIF
+      ENDIF
+      !
+      IF( zmin .LE. 0.0) THEN  ! we've got a problem
+         IF (lk_mpp) THEN
+            CALL mpp_minloc ( zalk(:,:),tmask(:,:,1), zmin, ii,ij )
+         ELSE
+            ilocs = MINLOC( zalk(:,:), mask = tmask(:,:,1) == 1. )
+            ii = ilocs(1) + nimpp - 1
+            ij = ilocs(2) + njmpp - 1
+         ENDIF
+         !
+         IF(lwp) THEN
+            WRITE(numout,*) 'trc_bio:tracer anomaly: *****    WARNING     *****'
+            WRITE(numout,*) 'trc_bio:tracer anomaly:  ALK concentration <= 0 '
+            WRITE(numout,9900) kt, zmin, ii, ij, jk
+            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****'
+         ENDIF
+      ENDIF
+
+
+9600  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' max DIC: ',f16.10,', i j k: ',3i5)
+9700  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' min DIC: ',f16.10,', i j k: ',3i5)
+9800  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' max ALK: ',f16.10,', i j k: ',3i5)
+9900  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' min ALK: ',f16.10,', i j k: ',3i5)
+
+   END SUBROUTINE trc_bio_check
+
+
 #else
    !!=====================================================================
@@ -663,4 +1012,5 @@
 CONTAINS
    SUBROUTINE trc_bio_medusa( kt )                   ! Empty routine
+      IMPLICIT NONE
       INTEGER, INTENT( in ) ::   kt
       WRITE(*,*) 'trc_bio_medusa: You should not have seen this print! error?', kt
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcdms_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcdms_medusa.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcdms_medusa.F90	(revision 10302)
@@ -113,4 +113,6 @@
 !
 ! AXY (13/03/15): Anderson et al. (2001)
+!! JPALM --19-12-2017-- Tunable through the namelist
+!!                      within dmsmin - dmscut - dmsslp
         Jterm = xqsr + 1.0e-6
         !! this next line makes a hard-coded assumption about the 
@@ -120,8 +122,8 @@
         Qterm = xdin / (xdin + 0.5)
         fq1 = log10(CHL * Jterm * Qterm)
-        if (fq1 > 1.72) then
-           dms_andr = (8.24 * (fq1 - 1.72)) + 2.29
-        else
-           dms_andr = 2.29
+        if (fq1 > dmscut) then
+           dms_andr = (dmsslp * (fq1 - dmscut)) + dmsmin
+        else
+           dms_andr = dmsmin
         endif
 !
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcini_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcini_medusa.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcini_medusa.F90	(revision 10302)
@@ -271,5 +271,5 @@
       zn_dms_srf(:,:)  = 0.0
       za_dms_srf(:,:)  = 0.0
-      zn_chl_srf(:,:)  = 2.0E-8 !! Chl srf
+      zn_chl_srf(:,:)  = 2.0E-8 !! Chl cpl - set first as surf
       !!
       IF(lwp) WRITE(numout,*) ' trc_ini_medusa: DMS and CO2 flux (UKESM) initialised to zero'
@@ -278,5 +278,5 @@
          CO2Flux_out_cpl(:,:) =  zn_co2_flx(:,:)   !! Coupling variable
          DMS_out_cpl(:,:)     =  zn_dms_srf(:,:)   !! Coupling variable
-         chloro_out_cpl(:,:)  =  zn_chl_srf(:,:)   !! Coupling variable
+         chloro_out_cpl(:,:)  =  zn_chl_srf(:,:) * scl_chl   !! Coupling variable
       END IF
       !!
@@ -324,4 +324,14 @@
       IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
       IF(lwp) CALL flush(numout)
+      !!
+      !!----------------------------------------------------------------------
+      !! JPALM (23-01-2017): new way to initialize CO2-atm for cmip6 
+      !!                     initially done in trcsms_medusa
+      !!----------------------------------------------------------------------
+      !! 
+      IF( ( .NOT.lk_oasis ) .AND. ( .NOT.lk_pi_co2 ) .AND. ( xobs_xco2a <= 0.0 ) ) THEN
+         IF(lwp) WRITE(numout,*) ' trc_ini_medusa: initialisating atm CO2 record'
+         CALL trc_ini_medusa_co2atm
+      ENDIF
 
    END SUBROUTINE trc_ini_medusa
@@ -480,4 +490,88 @@
    END SUBROUTINE trc_ini_medusa_river
    
+   SUBROUTINE trc_ini_medusa_co2atm
+      !!----------------------------------------------------------------------
+      !!                     ***  trc_ini_medusa_co2atm  ***  
+      !!
+      !! ** Purpose :   initialization atmospheric co2 record
+      !!
+      !! ** Method  : - Read the xco2 file
+      !!----------------------------------------------------------------------
+      INTEGER                       ::  jn, jm, io, ierr, inum, iostatus
+      INTEGER, PARAMETER            ::  iskip = 4   ! number of 1st descriptor lines
+      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   zyy !: xCO2 record years
+      CHARACTER (len=10), PARAMETER ::  clname = 'xco2.atm'  !! atm CO2 record file
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' trc_ini_medusa_co2atm: initialisation of atm CO2 historical record'
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
+
+
+      IF(lwp) WRITE(numout,*) 'read of formatted file xco2.atm'
+
+      CALL ctl_opn( inum, clname, 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
+      !!!
+      ! -Compute the number of year in the file
+      ! -File starts in co2_yinit, jn represents the record number in the file.
+      ! -Remove the file head (iskip lines) to jn
+      ! -The year is jn + yinit - 1 
+      !! Determine the number of lines in xCO2 input file
+      iostatus = 0
+      jn = 1
+      DO WHILE ( iostatus == 0 )
+        READ(inum,'(1x)', IOSTAT=iostatus, END=100)
+        jn = jn + 1
+      ENDDO
+      IF( iostatus .NE. 0 ) THEN
+        !! Error while reading xCO2 input file 
+        CALL ctl_stop('trc_ini_medusa_co2atm: &
+                      & Error on the 1st reading of xco2.atm')
+        RETURN
+      ENDIF
+ 100  co2_rec = jn - 1 - iskip
+      IF ( lwp) WRITE(numout,*) '    ', co2_rec ,' years read in the file'
+      !                                ! Allocate CO2 hist arrays
+      ierr = 0 
+      ALLOCATE( hist_pco2(co2_rec),zyy(co2_rec), STAT=ierr )
+      IF( ierr > 0 ) THEN
+         CALL ctl_stop( 'trc_ini_medusa_co2atm: unable to allocate  array' )  
+         RETURN
+      ENDIF
+
+      REWIND(inum)
+
+      DO jm = 1, iskip        ! Skip over 1st six descriptor lines
+         READ(inum,'(1x)')
+      END DO
+      ! file starts in 1931 do jn represent the year in the century.jhh
+      ! Read file till the end
+      ! allocate start and end year of the file
+      DO jn = 1, co2_rec
+        READ(inum,'(F6.1,F12.7)', IOSTAT=io) zyy(jn), hist_pco2(jn)
+        IF( io .NE. 0 ) THEN
+          !! Error while reading xCO2 input file 
+          CALL ctl_stop('trc_ini_medusa_co2atm: &
+                        & Error on the 2nd reading of xco2.atm')
+          RETURN
+        ENDIF
+
+        IF(jn==1) co2_yinit = zyy(jn)
+      END DO
+      co2_yend = co2_yinit + real(co2_rec) - 1.
+
+      IF(lwp) THEN        ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'CO2 hist start year: ', co2_yinit
+         WRITE(numout,*) 'CO2 hist end   year: ', co2_yend
+         WRITE(numout,*) ' Year   xCO2 atm '
+         DO jn = 1, co2_rec
+            WRITE(numout, '(F6.1,F12.7)') zyy(jn), hist_pco2(jn)
+         END DO
+      ENDIF
+
+   END SUBROUTINE trc_ini_medusa_co2atm
+
+
 #else
    !!----------------------------------------------------------------------
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 10302)
@@ -24,4 +24,5 @@
    USE sms_medusa      ! sms trends
    USE iom             ! I/O manager
+   USE sbc_oce, ONLY: lk_oasis
    !!USE trc_nam_dia     ! JPALM 13-11-2015 -- if iom_use for diag
 
@@ -85,7 +86,10 @@
       &  jriver_n,jriver_si,jriver_c,jriver_alk,jriver_dep,   &
       &  xsdiss,                                              &
-      &  sedlam,sedlostpoc,jpkb,jdms,jdms_input,jdms_model
+      &  sedlam,sedlostpoc,jpkb,jdms,jdms_input,jdms_model,   &
+      &  scl_chl, chl_out, dmsmin, dmscut, dmsslp,            &
+      &  ln_foam_medusa
 #if defined key_roam
       NAMELIST/natroam/ xthetaphy,xthetazoo,xthetanit,        &
+      &    xobs_xco2a,                                        &
       &    xthetarem,xo2min 
 #endif
@@ -125,4 +129,6 @@
       ! 1.4 namelist natbio : biological parameters
       ! -------------------------------------------
+      !! Note: the default values below will all be overwritten by the
+      !!       input in the namelist natbio.
       
       xxi         = 0.
@@ -246,4 +252,11 @@
       jdms_input  = 0
       jdms_model  = 0
+      scl_chl     = 1.
+      chl_out     = 1
+      dmsmin      = 2.29 !! Anderson DMS default
+      dmscut      = 1.72 !! Anderson DMS default
+      dmsslp      = 8.24 !! Anderson DMS default
+!!
+      ln_foam_medusa = .FALSE.
             
       !REWIND(numnatm)
@@ -399,5 +412,15 @@
 !!       jdms_model  :  choice of DMS model passed to atmosphere
 !!                      1 = ANDR, 2 = SIMO, 3 = ARAN, 4 = HALL, 5 = ANDM
-!!
+!!       dmsmin      : DMS minimum value for DMS Anderson (ANDR) sheme ONLY
+!!       dmscut      : DMS cutoff value for DMS Anderson (ANDR) sheme ONLY
+!!       dmsslp      : DMS slope value for DMS Anderson (ANDR) sheme ONLY
+!! UKESM1 - exported Chl to UM
+!!       scl_chl     : scaling factor to tune the chl field sent to the UM
+!!       chl_out     : select the chl field to send at the UM:
+!!                     1- Surf Chl ; 2- MLD Chl 
+!!
+!! FOAM - observation operator and data assimilation
+!!       ln_foam_medusa : calculate required diagnostics
+
       IF(lwp) THEN
 !!
@@ -909,30 +932,58 @@
 !!
 !! UKESM1 - new diagnostics  !! Jpalm; AXY (08/07/15)
-         WRITE(numout,*) '=== UKESM1-related parameters'
-         WRITE(numout,*)     &
-         &   ' include DMS diagnostic?,                                   jdms        = ', jdms
-         if (jdms_input .eq. 0) then
-            WRITE(numout,*)     &
-            &   ' use instant (0) or diel-avg (1) inputs,                    jdms_input  = instantaneous'
-         else
-            WRITE(numout,*)     &
-            &   ' use instant (0) or diel-avg (1) inputs,                    jdms_input  = diel-average'
-         endif
-	 if (jdms_model .eq. 1) then
-            WRITE(numout,*)     &
-            &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Anderson et al. (2001)'
-	 elseif (jdms_model .eq. 2) then
-            WRITE(numout,*)     &
-            &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Simo & Dachs (2002)'
-	 elseif (jdms_model .eq. 3) then
-            WRITE(numout,*)     &
-            &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Aranami & Tsunogai (2004)'
-	 elseif (jdms_model .eq. 4) then
-            WRITE(numout,*)     &
-            &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Halloran et al. (2010)'
-	 elseif (jdms_model .eq. 5) then
-            WRITE(numout,*)     &
-            &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Anderson et al. (2001; modified)'
-         endif
+         WRITE(numout,*) '=== UKESM1-related parameters ==='
+         WRITE(numout,*) ' ---- --- ---'
+
+         IF (lk_oasis) THEN
+            WRITE(numout,*) '=== UKESM1 --  coupled DMS to the atmosphere'
+            WRITE(numout,*)     &
+            &   ' include DMS diagnostic?,                                   jdms        = ', jdms
+            if (jdms_input .eq. 0) then
+               WRITE(numout,*)     &
+               &   ' use instant (0) or diel-avg (1) inputs,                    jdms_input  = instantaneous'
+            else
+               WRITE(numout,*)     &
+               &   ' use instant (0) or diel-avg (1) inputs,                    jdms_input  = diel-average'
+            endif
+   	    if (jdms_model .eq. 1) then
+               WRITE(numout,*)     &
+               &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Anderson et al. (2001)'
+	    elseif (jdms_model .eq. 2) then
+               WRITE(numout,*)     &
+               &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Simo & Dachs (2002)'
+	    elseif (jdms_model .eq. 3) then
+               WRITE(numout,*)     &
+               &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Aranami & Tsunogai (2004)'
+	    elseif (jdms_model .eq. 4) then
+               WRITE(numout,*)     &
+               &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Halloran et al. (2010)'
+	    elseif (jdms_model .eq. 5) then
+               WRITE(numout,*)     &
+               &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Anderson et al. (2001; modified)'
+            endif
+            if (jdms_model .eq. 1) then
+               WRITE(numout,*)     &
+               &   ' Anderson DMS model tuned parameters:                       DMS minimum = ',dmsmin,'. -- Default = 2.29 '
+               WRITE(numout,*)     &
+               &   ' Anderson DMS model tuned parameters:                       DMS cutoff  = ',dmscut,'. -- Default = 1.72 '
+               WRITE(numout,*)     &
+               &   ' Anderson DMS model tuned parameters:                       DMS slope   = ',dmsslp,'. -- Default = 8.24 '
+            endif
+
+            WRITE(numout,*) '=== UKESM1 --  coupled Chl to the atmosphere'
+            WRITE(numout,*)        &
+               &   ' Scaling factor to export tuned Chl to the atmosphere       scl_chl  = ', scl_chl
+            IF (chl_out .eq. 1) THEN
+               WRITE(numout,*)        &
+               &   ' Chl field to be scaled and sent to the atmosphere:         chl_out  = Surface Chl field '
+            ELSEIF (chl_out .eq. 2) THEN
+               WRITE(numout,*)        &
+               &   ' Chl field to be scaled and sent to the atmosphere:         chl_out  = MLD Chl field '
+            ENDIF
+         ENDIF ! IF lk_oasis=true
+!! FOAM
+         WRITE(numout,*) '=== FOAM-related parameters'
+         WRITE(numout,*)     &
+         &   ' calculate diagnostics for data assimilation,            ln_foam_medusa = ', ln_foam_medusa
 !!
       ENDIF
@@ -990,4 +1041,5 @@
       xthetarem = 0.
       xo2min    = 0.
+      xobs_xco2a = 0.
 
       !READ(numnatm,natroam)
@@ -1009,4 +1061,5 @@
 !!       xthetarem :  oxygen consumption by carbon remineralisation
 !!       xo2min    :  oxygen minimum concentration
+!!       xobs_xco2a : observed atmospheric xCO2 (not used if <= 0)
 
       IF(lwp) THEN
@@ -1026,4 +1079,6 @@
           WRITE(numout,*)     &
           &   ' oxygen minimum concentration                               xo2min      = ', xo2min
+          WRITE(numout,*)     &
+          &   ' observed atmospheric xCO2 (not used if <= 0)               xobs_xco2a  = ', xobs_xco2a
        ENDIF
 
@@ -2053,5 +2108,16 @@
           med_diag%OCN_DPCO2%dgsave = .FALSE.
       ENDIF
-      !!
+      !! UKESM additional
+      IF  (iom_use("CHL_MLD")) THEN 
+          med_diag%CHL_MLD%dgsave = .TRUE.
+      ELSE 
+          med_diag%CHL_MLD%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("CHL_CPL")) THEN 
+          med_diag%CHL_CPL%dgsave = .TRUE.
+      ELSE 
+          med_diag%CHL_CPL%dgsave = .FALSE.
+      ENDIF
+      !! 3D
       IF  (iom_use("TPP3")) THEN 
           med_diag%TPP3%dgsave = .TRUE.
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsms_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsms_medusa.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsms_medusa.F90	(revision 10302)
@@ -8,4 +8,5 @@
    !!              -   !  2008-11  (A. Yool) continuing adaptation for MEDUSA
    !!              -   !  2010-03  (A. Yool) updated for branch inclusion
+   !!              -   !  2017-08  (A. Yool) amend for slow detritus bug
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -22,4 +23,9 @@
    USE trcsed_medusa
    USE trcavg_medusa
+   !! for SMS trends
+   USE par_medusa,    ONLY: jp_msa0, jp_msa1, jp_medusa
+   USE par_oce,       ONLY: jpi, jpj, jpk
+   USE trd_oce,       ONLY: jptra_sms, l_trdtrc
+   USE trdtrc
 
 
@@ -46,4 +52,9 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) :: kt   ! ocean time-step index
+      !! Loop variables
+      INTEGER :: jn
+      !! trend temporary array:
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrmed
+
 
 # if defined key_debug_medusa
@@ -57,4 +68,13 @@
        IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
       ENDIF
+
+      !! MEDUSA SMS trends:
+      IF( l_trdtrc ) THEN
+          CALL wrk_alloc( jpi, jpj, jpk, jp_medusa, ztrmed )
+          ztrmed(:,:,:,:)=0.0  
+          DO jn = 1, jp_medusa
+            ztrmed(:,:,:,jn) = tra(:,:,:,jp_msa0 + jn - 1)
+          END DO
+      END IF
 
       CALL trc_avg_medusa( kt ) ! rolling average module
@@ -88,10 +108,21 @@
 #  endif
       
-      CALL trc_sed_medusa( kt ) ! sedimentation model
-#  if defined key_debug_medusa
-      IF(lwp) WRITE(numout,*) ' MEDUSA done trc_sed_medusa'
-      CALL flush(numout)
-#  endif
+!! AXY (08/08/2017): remove call to buggy subroutine (now handled by detritus.F90)
+!!       CALL trc_sed_medusa( kt ) ! sedimentation model
+!! #  if defined key_debug_medusa
+!!       IF(lwp) WRITE(numout,*) ' MEDUSA done trc_sed_medusa'
+!!       CALL flush(numout)
+!! #  endif
 # endif
+
+      !! MEDUSA SMS trends:
+      IF( l_trdtrc ) THEN
+          DO jn = 1, jp_medusa
+            ztrmed(:,:,:,jn) = tra(:,:,:,jp_msa0 + jn - 1)-ztrmed(:,:,:,jn)
+            CALL trd_trc( ztrmed(:,:,:,jn), jn, jptra_sms, kt )   ! save trends
+          END DO
+          CALL wrk_dealloc( jpi, jpj, jpk, jp_medusa, ztrmed )
+      END IF
+
 
    END SUBROUTINE trc_sms_medusa
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/zooplankton.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/zooplankton.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/zooplankton.F90	(revision 10302)
@@ -41,5 +41,8 @@
                                    idf, idfval,                          &
                                    zdet, zdtc, zphd, zphn, zzme, zzmi
-      USE dom_oce,           ONLY: e3t_0, e3t_n, tmask
+      USE dom_oce,           ONLY: e3t_0, tmask
+#if defined key_vvl
+      USE dom_oce,           ONLY: e3t_n
+#endif
       USE par_kind,          ONLY: wp
       USE in_out_manager,    ONLY: lwp, numout
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/TRP/trcbbl.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/TRP/trcbbl.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/TRP/trcbbl.F90	(revision 10302)
@@ -51,7 +51,9 @@
       !!
       !!----------------------------------------------------------------------  
+      IMPLICIT NONE
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step 
       CHARACTER (len=22) :: charout
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrtrd
+      INTEGER :: jn                   ! Local loop index
       !!----------------------------------------------------------------------
       !
@@ -108,4 +110,6 @@
 CONTAINS
    SUBROUTINE trc_bbl( kt )              ! Empty routine
+      IMPLICIT NONE
+      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
       WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
    END SUBROUTINE trc_bbl
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/TRP/trcnxt.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/TRP/trcnxt.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/TRP/trcnxt.F90	(revision 10302)
@@ -27,7 +27,9 @@
    !!----------------------------------------------------------------------
    USE oce_trc         ! ocean dynamics and tracers variables
+   USE domvvl          ! variable volume  
    USE trc             ! ocean passive tracers variables
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE prtctl_trc      ! Print control for debbuging
+   USE trcnam_trp      ! passive tracers transport namelist variables
    USE trd_oce
    USE trdtra
@@ -45,4 +47,6 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   r2dt
 
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
@@ -124,5 +128,37 @@
       IF( l_trdtrc )  THEN
          CALL wrk_alloc( jpi, jpj, jpk, jptra, ztrdt )  !* store now fields before applying the Asselin filter
-         ztrdt(:,:,:,:)  = trn(:,:,:,:)
+         ztrdt(:,:,jpk,:) = 0._wp
+         IF( ln_trcldf_iso ) THEN                       ! diagnose the "pure" Kz diffusive trend 
+            DO jn = 1, jptra
+               CALL trd_tra( kt, 'TRC', jn, jptra_zdfp, ztrdt(:,:,:,jn) )
+            ENDDO
+         ENDIF
+         ! total trend for the non-time-filtered variables.
+         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn
+         ! cancel from tsn terms
+         IF( lk_vvl ) THEN
+            DO jn = 1, jptra
+               DO jk = 1, jpkm1
+                  zfact = 1.0 / rdttrc(jk)
+                  ztrdt(:,:,jk,jn) = ( tra(:,:,jk,jn)*fse3t_a(:,:,jk) / fse3t_n(:,:,jk) - &
+                                       trn(:,:,jk,jn) ) * zfact
+               END DO
+            END DO
+         ELSE
+            DO jn = 1, jptra
+               DO jk = 1, jpkm1
+                  zfact = 1.0 / rdttrc(jk)
+                  ztrdt(:,:,jk,jn) = ( tra(:,:,jk,jn) - trn(:,:,jk,jn) ) * zfact
+               END DO
+            END DO
+         END IF
+         DO jn = 1, jptra
+            CALL trd_tra( kt, 'TRC', jn, jptra_tot, ztrdt(:,:,:,jn) )
+         ENDDO
+         IF( .NOT.lk_vvl )  THEN
+            ! Store now fields before applying the Asselin filter 
+            ! in order to calculate Asselin filter trend later.
+            ztrdt(:,:,:,:)  = trn(:,:,:,:)
+         ENDIF
       ENDIF
       ! Leap-Frog + Asselin filter time stepping
@@ -134,4 +170,15 @@
             END DO
          END DO
+         IF (l_trdtrc.AND.lk_vvl) THEN      ! Zero Asselin filter contribution
+                                            ! must be explicitly written out since for vvl
+                                            ! Asselin filter is output by
+                                            ! tra_nxt_vvl that is not called on
+                                            ! this time step
+            ztrdt(:,:,:,:) = 0._wp
+            DO jn = 1, jptra
+               CALL trd_tra( kt, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
+            ENDDO
+         END IF
+
          !                                              
       ELSE
@@ -144,14 +191,15 @@
 
       ! trends computation
-      IF( l_trdtrc ) THEN                                      ! trends
+      IF( l_trdtrc.AND..NOT.lk_vvl) THEN                                      ! trends
          DO jn = 1, jptra
             DO jk = 1, jpkm1
                zfact = 1.e0 / r2dt(jk)  
                ztrdt(:,:,jk,jn) = ( trb(:,:,jk,jn) - ztrdt(:,:,jk,jn) ) * zfact 
-               CALL trd_tra( kt, 'TRC', jn, jptra_atf, ztrdt )
-            END DO
+            END DO
+            CALL trd_tra( kt, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
          END DO
-         CALL wrk_dealloc( jpi, jpj, jpk, jptra, ztrdt ) 
       END IF
+      !
+      IF( l_trdtrc)  CALL wrk_dealloc( jpi, jpj, jpk, jptra, ztrdt ) 
       !
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90	(revision 10302)
@@ -140,6 +140,9 @@
       DO jn = 1, jptra
          !
-         IF( l_trdtrc ) ztrtrd(:,:,:) = tra(:,:,:,jn)  ! save trends
-         !                                             ! add the trend to the general tracer trend
+         IF( l_trdtrc ) THEN
+            ztrtrd(:,:,:) = 0.0
+            ztrtrd(:,:,1) = tra(:,:,1,jn)  ! save surface trends
+         !                                 ! add the trend to the general tracer trend
+         ENDIF
 
          IF ( nn_ice_tr == -1 ) THEN  ! No tracers in sea ice (null concentration in sea ice)
@@ -184,5 +187,5 @@
          !
          IF( l_trdtrc ) THEN
-            ztrtrd(:,:,:) = tra(:,:,:,jn) - ztrtrd(:,:,:)
+            ztrtrd(:,:,1) = tra(:,:,1,jn) - ztrtrd(:,:,1)
             CALL trd_tra( kt, 'TRC', jn, jptra_nsr, ztrtrd )
          END IF
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 10302)
@@ -28,5 +28,5 @@
    USE zpshde          ! partial step: hor. derivative       (zps_hde routine)
 # if defined key_debug_medusa
-   USE trcrst
+   USE trcstat
 # endif
 
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90	(revision 10302)
@@ -15,4 +15,5 @@
    !!----------------------------------------------------------------------
    USE oce_trc         ! ocean dynamics and active tracers
+   USE domvvl          ! variable volume       
    USE trc             ! ocean passive tracers variables
    USE trcnam_trp      ! passive tracers transport namelist variables
@@ -98,8 +99,21 @@
 
       IF( l_trdtrc )   THEN                      ! save the vertical diffusive trends for further diagnostics
+         !! JPALM -- 18-08-2017 -- vvl case, do as done by G Nurser in trazdf 
+         IF( lk_vvl ) THEN
+            DO jn = 1, jptra
+               DO jk = 1, jpkm1
+                  ztrtrd(:,:,jk,jn) = ( ( tra(:,:,jk,jn)*fse3t_a(:,:,jk) - &
+                                          trb(:,:,jk,jn)*fse3t_b(:,:,jk) ) &
+                                       / (fse3t_n(:,:,jk)*r2dt(jk)) ) - ztrtrd(:,:,jk,jn)
+               END DO
+            END DO
+         ELSE
+            DO jn = 1, jptra
+               DO jk = 1, jpkm1
+                  ztrtrd(:,:,jk,jn) = ( ( tra(:,:,jk,jn) - trb(:,:,jk,jn) ) / r2dt(jk) ) - ztrtrd(:,:,jk,jn)
+               END DO
+            END DO
+         ENDIF 
          DO jn = 1, jptra
-            DO jk = 1, jpkm1
-               ztrtrd(:,:,jk,jn) = ( ( tra(:,:,jk,jn) - trb(:,:,jk,jn) ) / r2dt(jk) ) - ztrtrd(:,:,jk,jn)
-            END DO
             CALL trd_tra( kt, 'TRC', jn, jptra_zdf, ztrtrd(:,:,:,jn) )
          END DO
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/TRP/trdtrc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/TRP/trdtrc.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/TRP/trdtrc.F90	(revision 10302)
@@ -19,4 +19,7 @@
    USE trdmxl_trc        ! Mixed layer trends diag.
    USE iom               ! I/O library
+# if defined key_debug_medusa
+   USE trcstat,          ONLY: trc_rst_dia_stat     
+# endif
 
    IMPLICIT NONE
@@ -86,21 +89,70 @@
 
       IF( lk_trdtrc .AND. ln_trdtrc( kjn ) ) THEN
-         !
+      !! JPALM -- 17-08-2017 -- modif following trd_tra_iom as suggested by Georges
+      !!                     -- add jptra_tot; jptra_totad; jptra_zdfp
+      !!                     -- shange to output trends every 2 time-step, except tot.
+      !!                     -- move cltra and iomput inside the select case
+      !!                     So if an non-wanted case arrives here it will not go
+      !!                     through cltra (without value) and break iomput.
+      !!                     -- Add iom_use in prevision of not using All trends
+      !!                     for All passive tracers (will create a HUGE 3D file otherwise --
+      !!                     might be interested in very few of them : SMS and TOT probably)
+         !
+         SELECT CASE( ktrd )
+         !! tot - output every time-step:
+         CASE( jptra_tot  )       ;    WRITE (cltra,'("TOT_",4a)')
+                           cltra = TRIM(cltra)//TRIM(ctrcnm(kjn))
+                           CALL trd_trc_iomput( cltra, ptrtrd, kjn, kt )
+         END SELECT
+         !
+       IF( MOD( kt, 2 ) == 0 ) THEN
          SELECT CASE( ktrd )
          CASE( jptra_xad  )       ;    WRITE (cltra,'("XAD_",4a)')
+                           cltra = TRIM(cltra)//TRIM(ctrcnm(kjn))
+                           CALL trd_trc_iomput( cltra, ptrtrd, kjn, kt )
          CASE( jptra_yad  )       ;    WRITE (cltra,'("YAD_",4a)')
-         CASE( jptra_zad  )       ;    WRITE (cltra,'("ZAD_",4a)')
+                           cltra = TRIM(cltra)//TRIM(ctrcnm(kjn))
+                           CALL trd_trc_iomput( cltra, ptrtrd, kjn, kt )
+         CASE( jptra_zad  )       ;    WRITE (cltra,'("ZAD_",4a)')      !! care vvl case
+                           cltra = TRIM(cltra)//TRIM(ctrcnm(kjn))
+                           CALL trd_trc_iomput( cltra, ptrtrd, kjn, kt )
+         CASE( jptra_totad  )     ;    WRITE (cltra,'("TAD_",4a)')      !! total adv
+                           cltra = TRIM(cltra)//TRIM(ctrcnm(kjn))
+                           CALL trd_trc_iomput( cltra, ptrtrd, kjn, kt )
          CASE( jptra_ldf  )       ;    WRITE (cltra,'("LDF_",4a)')
+                           cltra = TRIM(cltra)//TRIM(ctrcnm(kjn))
+                           CALL trd_trc_iomput( cltra, ptrtrd, kjn, kt )
          CASE( jptra_bbl  )       ;    WRITE (cltra,'("BBL_",4a)')
+                           cltra = TRIM(cltra)//TRIM(ctrcnm(kjn))
+                           CALL trd_trc_iomput( cltra, ptrtrd, kjn, kt )
          CASE( jptra_nsr  )       ;    WRITE (cltra,'("FOR_",4a)')
+                           cltra = TRIM(cltra)//TRIM(ctrcnm(kjn))
+                           CALL trd_trc_iomput( cltra, ptrtrd, kjn, kt )
          CASE( jptra_zdf  )       ;    WRITE (cltra,'("ZDF_",4a)')
+                           cltra = TRIM(cltra)//TRIM(ctrcnm(kjn))
+                           CALL trd_trc_iomput( cltra, ptrtrd, kjn, kt )
+         CASE( jptra_zdfp )       ;    WRITE (cltra,'("ZDP_",4a)')
+                           cltra = TRIM(cltra)//TRIM(ctrcnm(kjn))
+                           CALL trd_trc_iomput( cltra, ptrtrd, kjn, kt )
          CASE( jptra_dmp  )       ;    WRITE (cltra,'("DMP_",4a)')
+                           cltra = TRIM(cltra)//TRIM(ctrcnm(kjn))
+                           CALL trd_trc_iomput( cltra, ptrtrd, kjn, kt )
          CASE( jptra_sms  )       ;    WRITE (cltra,'("SMS_",4a)')
+                           cltra = TRIM(cltra)//TRIM(ctrcnm(kjn))
+                           CALL trd_trc_iomput( cltra, ptrtrd, kjn, kt )
+         CASE( jptra_radb )       ;    WRITE (cltra,'("RDB_",4a)')
+                           cltra = TRIM(cltra)//TRIM(ctrcnm(kjn))
+                           CALL trd_trc_iomput( cltra, ptrtrd, kjn, kt )
+         CASE( jptra_radn )       ;    WRITE (cltra,'("RDN_",4a)')
+                           cltra = TRIM(cltra)//TRIM(ctrcnm(kjn))
+                           CALL trd_trc_iomput( cltra, ptrtrd, kjn, kt )
+         END SELECT
+       ELSE IF( MOD( kt, 2 ) == 1 ) THEN
+         SELECT CASE( ktrd )
          CASE( jptra_atf  )       ;    WRITE (cltra,'("ATF_",4a)')
-         CASE( jptra_radb )       ;    WRITE (cltra,'("RDB_",4a)')
-         CASE( jptra_radn )       ;    WRITE (cltra,'("RDN_",4a)')
-         END SELECT
-                                          cltra = TRIM(cltra)//TRIM(ctrcnm(kjn))
-                                          CALL iom_put( cltra,  ptrtrd(:,:,:) )
+                           cltra = TRIM(cltra)//TRIM(ctrcnm(kjn))
+                           CALL trd_trc_iomput( cltra, ptrtrd, kjn, kt )
+         END SELECT
+       END IF
          !
       END IF
@@ -123,4 +175,33 @@
 
    END SUBROUTINE trd_trc_bio
+
+   SUBROUTINE trd_trc_iomput( cltra, ptrtrd, kjn, kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_trc_iomput  ***
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in )  ::   kt                                  ! timestep
+      INTEGER, INTENT( in )  ::   kjn                                 ! biotrend index
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout )  ::   ptrtrd  ! var trend
+      CHARACTER (len=*),INTENT( in ) :: cltra                         ! trend name
+      !!----------------------------------------------------------------------
+
+
+      IF  (iom_use(cltra)) THEN
+# if defined key_debug_medusa
+         IF(lwp) WRITE(numout,*) ' TREND stats (min, max,sum) kt = ',kt ,' jn = ',kjn
+         CALL trc_rst_dia_stat( ptrtrd(:,:,1), cltra)
+# endif
+         CALL iom_put( cltra,  ptrtrd(:,:,:) )
+# if defined key_debug_medusa
+      ELSE
+         IF(lwp) WRITE(numout,*) &
+                      ' TREND -- No output asked for ',cltra,' kt = ',kt,' jn = ',kjn
+         CALL trc_rst_dia_stat( ptrtrd(:,:,1), cltra)
+# endif
+      ENDIF
+
+   END SUBROUTINE trd_trc_iomput
+
+
 #else
    !!----------------------------------------------------------------------
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/trc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 10302)
@@ -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- 
@@ -105,5 +110,5 @@
    END TYPE DIAG
 
-#if defined key_medusa && defined key_iomput
+#if defined key_medusa
    TYPE, PUBLIC :: BDIAG
       LOGICAL              :: dgsave
@@ -134,4 +139,6 @@
                   OCN_KWCO2, OCN_K0, CO2STARAIR, OCN_DPCO2,                                          & ! end of regular 2D
                   TPP3, DETFLUX3, REMIN3N, PH3, OM_CAL3,                                             & ! end of regular 3D
+! JPALM (01/09/17): additional UKESM 2D diag
+                  CHL_MLD, CHL_CPL,                                                                  &
 ! AXY (11/11/16): additional CMIP6 2D diagnostics
                   epC100, epCALC100, epN100, epSI100,                                                &
@@ -264,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/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/trcini.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 10302)
@@ -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,7 @@
    USE sbc_oce
    USE trcice          ! tracers in sea ice
- 
+# if defined key_medusa
+   USE sms_medusa      ! MEDUSA   initialisation
+# endif
    IMPLICIT NONE
    PRIVATE
@@ -62,5 +65,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 +109,7 @@
       !                                                              ! total volume of the ocean 
       areatot = glob_sum( cvol(:,:,:) )
+# if defined key_medusa && defined key_roam
+      carea(:,:) = e1e2t(:,:) * tmask(:,:,1) 
+# endif
 
       IF( lk_pisces  )       CALL trc_ini_pisces       ! PISCES  bio-model
@@ -192,8 +206,83 @@
       ENDIF
 
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'trc_init : passive tracer set up completed'
-      IF(lwp) WRITE(numout,*) '~~~~~~~'
-      IF(lwp) CALL flush(numout)
+# 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: '
+      ! 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)
+      ! 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)
+      ! 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)
+      ! 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)
+      ! 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)
+      ! 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)
+      ! Check
+      zsum3d        = glob_sum( cvol(:,:,:) )
+      zsum2d        = glob_sum( carea(:,:) )      
+      IF ( lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) ' check : cvol    : ', zsum3d
+         WRITE(numout,*) ' check : carea   : ', zsum2d
+         WRITE(numout,*)
+      ENDIF
+      !
+# endif
+
+      IF(lwp) THEN 
+          WRITE(numout,*)
+          WRITE(numout,*) 'trc_init : passive tracer set up completed'
+          WRITE(numout,*) '~~~~~~~'
+      ENDIF 
 # if defined key_debug_medusa
          CALL trc_rst_stat
@@ -202,4 +291,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/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/trcnam.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/trcnam.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/trcnam.F90	(revision 10302)
@@ -60,5 +60,16 @@
       !!                ( (PISCES, CFC, MY_TRC, MEDUSA, IDTRA, Age )
       !!---------------------------------------------------------------------
-      INTEGER  ::   jn, jk                     ! dummy loop indice
+      INTEGER ::  ierr
+#if defined key_trdmxl_trc  || defined key_trdtrc
+      NAMELIST/namtrc_trd/ nn_trd_trc, nn_ctls_trc, rn_ucf_trc, &
+         &                ln_trdmxl_trc_restart, ln_trdmxl_trc_instant, &
+         &                cn_trdrst_trc_in, cn_trdrst_trc_out, ln_trdtrc
+#endif
+
+      INTEGER  ::   jn, jk              ! dummy loop indice
+      INTEGER  ::   ios                 ! Local integer output status for namelist read
+      !!---------------------------------------------------------------------
+
+
       !                                        !   Parameters of the run 
       IF( .NOT. lk_offline ) CALL trc_nam_run
@@ -68,5 +79,5 @@
       
       !                                        !   Parameters of additional diagnostics
-      CALL trc_nam_dia
+      IF( .NOT. lk_iomput )  CALL trc_nam_dia
 
       !                                        !   namelist of transport
@@ -171,6 +182,6 @@
       ENDIF
 
-      IF( lk_c14b     ) THEN   ;   CALL trc_nam_c14b         ! C14 bomb     tracers
-      ELSE                    ;   IF(lwp) WRITE(numout,*) '          C14 not used'
+      IF( lk_c14b    ) THEN  ;   CALL trc_nam_c14b         ! C14 bomb     tracers
+      ELSE                   ;   IF(lwp) WRITE(numout,*) '          C14 not used'
       ENDIF
 
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
 
 
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/trcstat.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/trcstat.F90	(revision 10302)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/trcstat.F90	(revision 10302)
@@ -0,0 +1,156 @@
+MODULE trcstat
+   !!======================================================================
+   !!                         ***  MODULE trcrst  ***
+   !! TOP :   Manage the passive tracer restart
+   !!======================================================================
+   !! History :    -   !  1991-03  ()  original code
+   !!             1.0  !  2005-03 (O. Aumont, A. El Moussaoui) F90
+   !!              -   !  2005-10 (C. Ethe) print control
+   !!             2.0  !  2005-10 (C. Ethe, G. Madec) revised architecture
+   !!----------------------------------------------------------------------
+#if defined key_top
+   !!----------------------------------------------------------------------
+   !!   'key_top'                                                TOP models
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   trc_rst :   Restart for passive tracer
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   'key_top'                                                TOP models
+   !!----------------------------------------------------------------------
+   !!   trc_rst_opn    : open  restart file
+   !!   trc_rst_read   : read  restart file
+   !!   trc_rst_wri    : write restart file
+   !!----------------------------------------------------------------------
+   USE trc,               ONLY: tra, ctrcnm
+   USE par_kind,          ONLY: wp
+   USE in_out_manager,    ONLY: lwp, numout
+   USE par_oce,           ONLY: jpi, jpj
+   USE par_trc,           ONLY: jptra
+   USE dom_oce,           ONLY: e3t_0, gdepw_0, tmask, e1e2t
+#if defined key_vvl
+   USE dom_oce,           ONLY: e3t_a, e3t_n, gdepw_n
+#endif
+   !* MPP library                         
+   USE lib_mpp
+   !* Fortran utilities                         
+   USE lib_fortran
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   trc_rst_dia_stat
+   PUBLIC   trc_rst_tra_stat
+
+   !! * Substitutions
+#  include "top_substitute.h90"
+   
+CONTAINS
+   
+   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
+      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
+
+
+#else
+   !!----------------------------------------------------------------------
+   !!  Dummy module :                                     No passive tracer
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_rst_dia_stat                      ! Empty routines
+      WRITE(*,*) 'trc_rst_wri: You should not have seen this print! error?'
+   END SUBROUTINE trc_rst_dia_stat
+   SUBROUTINE trc_rst_dia_stat( dgtr, names)
+      REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) ::   dgtr      ! 2D diag var
+      CHARACTER(len=*)             , INTENT(in) ::   names     ! 2D diag name
+      WRITE(*,*) 'trc_rst_wri: You should not have seen this print! error?'
+   END SUBROUTINE trc_rst_dia_stat  
+#endif
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!======================================================================
+END MODULE trcstat
Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/trcstp.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 9287)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 10302)
@@ -19,4 +19,5 @@
    USE trcwri
    USE trcrst
+   USE trcstat
    USE trdtrc_oce
    USE trdmxl_trc
@@ -177,7 +178,7 @@
          !
          !                                            !* Restart: read in restart file
-         IF( ln_rsttr .AND. iom_varid( numrtr, 'qsr_mean' , ldstop = .FALSE. ) > 0 .AND. &
-                            iom_varid( numrtr, 'qsr_arr_1', ldstop = .FALSE. ) > 0 .AND. &
-                            iom_varid( numrtr, 'zsecfst'  , ldstop = .FALSE. ) > 0 ) THEN 
+         IF( ln_rsttr .AND. nn_rsttr /= 0 .AND. iom_varid( numrtr, 'qsr_mean' , ldstop = .FALSE. ) > 0 &
+                                          .AND. iom_varid( numrtr, 'qsr_arr_1', ldstop = .FALSE. ) > 0 &
+                                          .AND. iom_varid( numrtr, 'zsecfst'  , ldstop = .FALSE. ) > 0 ) THEN
             IF(lwp) WRITE(numout,*) 'trc_qsr_mean:   qsr_mean read in the restart file'
             CALL iom_get( numrtr, jpdom_autoglo, 'qsr_mean', qsr_mean )   !  A mean of qsr
