Index: /branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90
===================================================================
--- /branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 8641)
+++ /branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 8642)
@@ -62,5 +62,5 @@
 # endif
                                    zchd, zchn, zdin, zsil
-      USE dom_oce,           ONLY: e3t_0, e3t_n, gphit, tmask
+      USE dom_oce,           ONLY: e3t_0, e3t_n, gphit, tmask, mig, mjg
 # if defined key_iomput
       USE iom,               ONLY: lk_iomput
@@ -91,4 +91,6 @@
       USE trcoxy_medusa,     ONLY: trc_oxy_medusa
 # endif
+      USE lib_mpp,           ONLY: ctl_stop
+      USE trcstat 
 
    !!* Substitution
@@ -122,4 +124,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
@@ -134,8 +151,15 @@
          DO ji = 2,jpim1
             !! OPEN wet point IF..THEN loop
-            if (tmask(ji,jj,1) == 1) then
+            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)
+                  !!! 
+                  !!! Jpalm test on atm xCO2
+                  IF ( (f_xco2a(ji,jj) > 1000 ).OR.(f_xco2a(ji,jj) < 100 ) ) THEN
+                     IF(lwp) WRITE(numout,*) ' atm xCO2 = ',f_xco2a(ji,jj),   &
+                                             ' -- ji =', mig(ji),' jj = ', mjg(jj)
+                     CALL ctl_stop( 'MEDUSA - Air-Sea :', 'unrealistic atm xCO2 ' )
+                 ENDIF 
                ENDIF
                !!
@@ -162,19 +186,48 @@
                '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
       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...i
+               !!          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 : unrealistic surface Carb. Chemistry'
+                      CALL ctl_stop( 'MEDUSA - Air-Sea :', 'unrealistic surface Carb. Chemistry -- INPUTS' )
+                  ENDIF     
+               ENDIF     
                !!
                !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate
@@ -201,8 +254,46 @@
                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...i
+               !!          Check MOCSY inputs
+              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-2 ).OR.(f_co2flux(ji,jj) < -1.E-2 ) .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 : unrealistic surface Carb. Chemistry'
+                     CALL ctl_stop( 'MEDUSA - Air-Sea :', 'unrealistic surface Carb. Chemistry -- OUTPUTS' )
+                 ENDIF     
+              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   
 
Index: /branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90
===================================================================
--- /branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 8641)
+++ /branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 8642)
@@ -173,4 +173,23 @@
       fslowsinkc(:,:) = 0.0
 # endif      
+      !!
+      !! 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
+      zpho(:,:)  = 0.0
+      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
       !!
       !! allocate and initiate 2D diag
Index: /branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
===================================================================
--- /branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 8641)
+++ /branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 8642)
@@ -177,4 +177,8 @@
       !! temporary variables
       REAL(wp) ::    fq0,fq1,fq2,fq3,fq4
+      !!
+      !! T and S check temporary variable
+      REAL(wp) ::    sumtsn, tsnavg
+      INTEGER  ::    summask
       !!
       !!------------------------------------------------------------------
@@ -474,19 +478,106 @@
                   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, ',                   &
+                        ' trc_bio_medusa: T WARNING 3D, ',                   &
                         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
+                     !! temperature
                      ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem)
+                      IF(lwp) WRITE(numout,*)                                 &
+                     ' trc_bio_medusa: Abnormal T suroundings Temp'
+                      IF(lwp) WRITE(numout,*)                                 &
+                     tsn(ji-1,jj+1,jk,jp_tem), tsn(ji,jj+1,jk,jp_tem), tsn(ji+1,jj+1,jk,jp_tem)
+                      IF(lwp) WRITE(numout,*)                                 &
+                     tsn(ji-1,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), tsn(ji+1,jj,jk,jp_tem)
+                      IF(lwp) WRITE(numout,*)                                 &
+                     tsn(ji-1,jj-1,jk,jp_tem), tsn(ji,jj-1,jk,jp_tem), tsn(ji+1,jj-1,jk,jp_tem)
+                     sumtsn = ( tmask(ji-1,jj+1,jk) * tsn(ji-1,jj+1,jk,jp_tem) ) + & 
+                              ( tmask(ji  ,jj+1,jk) * tsn(ji  ,jj+1,jk,jp_tem) ) + & 
+                              ( tmask(ji+1,jj+1,jk) * tsn(ji+1,jj+1,jk,jp_tem) ) + & 
+                              ( tmask(ji-1,jj  ,jk) * tsn(ji-1,jj  ,jk,jp_tem) ) + & 
+                              ( tmask(ji+1,jj  ,jk) * tsn(ji+1,jj  ,jk,jp_tem) ) + & 
+                              ( tmask(ji-1,jj-1,jk) * tsn(ji-1,jj-1,jk,jp_tem) ) + & 
+                              ( tmask(ji  ,jj-1,jk) * tsn(ji  ,jj-1,jk,jp_tem) ) + & 
+                              ( tmask(ji+1,jj-1,jk) * tsn(ji+1,jj-1,jk,jp_tem) )
+                      IF(lwp) WRITE(numout,*)                                 &
+                     ' trc_bio_medusa: Abnormal T suroundings Sal'
+                      IF(lwp) WRITE(numout,*)                                 &
+                     tsn(ji-1,jj+1,jk,jp_sal), tsn(ji,jj+1,jk,jp_sal), tsn(ji+1,jj+1,jk,jp_sal)
+                      IF(lwp) WRITE(numout,*)                                 &
+                     tsn(ji-1,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), tsn(ji+1,jj,jk,jp_sal)
+                      IF(lwp) WRITE(numout,*)                                 &
+                     tsn(ji-1,jj-1,jk,jp_sal), tsn(ji,jj-1,jk,jp_sal), tsn(ji+1,jj-1,jk,jp_sal)
+                     !!
+                      IF(lwp) WRITE(numout,*)                                 &
+                     ' trc_bio_medusa: Abnormal T suroundings DIC'
+                      IF(lwp) WRITE(numout,*)                                 &
+                     trn(ji-1,jj+1,jk,jpdic), trn(ji,jj+1,jk,jpdic), trn(ji+1,jj+1,jk,jpdic)
+                      IF(lwp) WRITE(numout,*)                                 &
+                     trn(ji-1,jj,jk,jpdic), trn(ji,jj,jk,jpdic), trn(ji+1,jj,jk,jpdic)
+                      IF(lwp) WRITE(numout,*)                                 &
+                     trn(ji-1,jj-1,jk,jpdic), trn(ji,jj-1,jk,jpdic), trn(ji+1,jj-1,jk,jpdic)
+                      IF(lwp) WRITE(numout,*)                                 &
+                     ' trc_bio_medusa: Abnormal T suroundings Alk'
+                      IF(lwp) WRITE(numout,*)                                 &
+                     trn(ji-1,jj+1,jk,jpalk), trn(ji,jj+1,jk,jpalk), trn(ji+1,jj+1,jk,jpalk)
+                      IF(lwp) WRITE(numout,*)                                 &
+                     trn(ji-1,jj,jk,jpalk), trn(ji,jj,jk,jpalk), trn(ji+1,jj,jk,jpalk)
+                      IF(lwp) WRITE(numout,*)                                 &
+                     trn(ji-1,jj-1,jk,jpalk), trn(ji,jj-1,jk,jpalk), trn(ji+1,jj-1,jk,jpalk)
+                      IF(lwp) WRITE(numout,*)                                 &
+                     ' trc_bio_medusa: Abnormal T suroundings tmask'
+                      IF(lwp) WRITE(numout,*)                                 &
+                     tmask(ji-1,jj+1,jk), tmask(ji,jj+1,jk), tmask(ji+1,jj+1,jk)
+                      IF(lwp) WRITE(numout,*)                                 &
+                     tmask(ji-1,jj,jk), tmask(ji,jj,jk), tmask(ji+1,jj,jk)
+                      IF(lwp) WRITE(numout,*)                                 &
+                     tmask(ji-1,jj-1,jk), tmask(ji,jj-1,jk), tmask(ji+1,jj-1,jk)
+                     summask = tmask(ji-1,jj+1,jk) + tmask(ji  ,jj+1,jk)   +  &
+                               tmask(ji+1,jj+1,jk) + tmask(ji-1,jj  ,jk)   +  &
+                               tmask(ji+1,jj  ,jk) + tmask(ji-1,jj-1,jk)   +  & 
+                               tmask(ji  ,jj-1,jk) + tmask(ji+1,jj-1,jk) 
+                     tsnavg = ( sumtsn / summask )
+                     !! Correct out of range values
+                     IF ( ( summask .EQ. 0.0 ) .OR. (tsnavg .LT. -3.0 ) .OR.   &
+                          ( tsnavg .GT. 40.0 ) ) THEN    
+                        IF (ztmp(ji,jj) .LT. -3.0 ) THEN
+                           IF(lwp) WRITE(numout,*)                           &
+                           ' trc_bio_medusa: T SWITCHING 3D, ',              &
+                           tsn(ji,jj,jk,jp_tem), ' -> -3.0 '
+                           ztmp(ji,jj) = -3.0
+                        ENDIF
+                        IF (ztmp(ji,jj) .GT. 40.0 ) THEN
+                           IF(lwp) WRITE(numout,*)                           &
+                           ' trc_bio_medusa: T SWITCHING 3D, ',              &
+                           tsn(ji,jj,jk,jp_tem), ' -> 40.0 '
+                           ztmp(ji,jj) = 40.0
+                        ENDIF
+                     ELSE 
+                        IF(lwp) WRITE(numout,*)                              &
+                        ' trc_bio_medusa: T SWITCHING 3D, ',                 &
+                        tsn(ji,jj,jk,jp_tem), ' -> surounding avg : ', tsnavg
+                        ztmp(ji,jj) = tsnavg
+                     ENDIF
                   endif
-                  if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then
+                  !! end T chack
+                  if (zsal(ji,jj) .lt. 1.0 .or. zsal(ji,jj) .gt. 47.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
+                     !! Correct out of range values
+                     IF (zsal(ji,jj) .LT. 1.0 ) THEN
+                        IF(lwp) WRITE(numout,*)                              &
+                        ' trc_bio_medusa: S SWITCHING 3D, ',                 &
+                        tsn(ji,jj,jk,jp_sal), ' -> 1.0 '
+                        zsal(ji,jj) = 1.0
+                     ENDIF
+                     IF (zsal(ji,jj) .GT. 47.0 ) THEN
+                        IF(lwp) WRITE(numout,*)                              &
+                        ' trc_bio_medusa: T SWITCHING 3D, ',                 &
+                        tsn(ji,jj,jk,jp_sal), ' -> 47.0 '
+                        zsal(ji,jj) = 47.0
+                     ENDIF
                   endif
+                  !! end S check
                ENDIF
             ENDDO
Index: /branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90
===================================================================
--- /branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 8641)
+++ /branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 8642)
@@ -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/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- /branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 8641)
+++ /branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 8642)
@@ -43,4 +43,5 @@
    USE sbc_oce, ONLY: lk_oasis 
    USE oce,     ONLY: CO2Flux_out_cpl, DMS_out_cpl, chloro_out_cpl  !! Coupling variable
+   USE trcstat
 
    IMPLICIT NONE
@@ -52,6 +53,4 @@
    PUBLIC   trc_rst_cal
    PUBLIC   trc_rst_stat
-   PUBLIC   trc_rst_dia_stat
-   PUBLIC   trc_rst_tra_stat
 
    !! * Substitutions
@@ -697,90 +696,4 @@
 
 
-   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
    !!----------------------------------------------------------------------
Index: /branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/trcstat.F90
===================================================================
--- /branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/trcstat.F90	(revision 8642)
+++ /branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/trcstat.F90	(revision 8642)
@@ -0,0 +1,154 @@
+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, e3t_a, e3t_n, gdepw_0, gdepw_n,     &
+                          tmask, e1e2t
+   !* 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/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/trcstp.F90
===================================================================
--- /branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 8641)
+++ /branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 8642)
@@ -19,4 +19,5 @@
    USE trcwri
    USE trcrst
+   USE trcstat
    USE trdtrc_oce
    USE trdmxl_trc
