Index: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 9070)
+++ branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 9073)
@@ -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
    !!----------------------------------------------------------------------
    !!
@@ -81,5 +82,6 @@
                                             gdepw_0, gdepw_n,               &
                                             nday_year, nsec_day, nyear,     &
-                                            rdt, tmask
+                                            rdt, tmask, mig, mjg, nimpp,    &
+                                            njmpp 
       USE in_out_manager,             ONLY: lwp, numout, nn_date0
 # if defined key_iomput
@@ -87,5 +89,5 @@
 # endif
       USE lbclnk,                     ONLY: lbc_lnk
-      USE lib_mpp,                    ONLY: ctl_stop
+      USE lib_mpp  
       USE oce,                        ONLY: tsb, tsn
       USE par_kind,                   ONLY: wp
@@ -115,4 +117,5 @@
       
       PUBLIC   trc_bio_medusa    ! called in trcsms_medusa.F90
+      PUBLIC   trc_bio_exceptional_fix ! here 
 
    !!* Substitution
@@ -177,4 +180,9 @@
       !! temporary variables
       REAL(wp) ::    fq0,fq1,fq2,fq3,fq4
+      !!
+      !! T and S check temporary variable
+      REAL(wp) ::    sumtsn, tsnavg
+      INTEGER  ::    summask
+      CHARACTER(40) :: charout, charout2, charout3, charout4, charout5
       !!
       !!------------------------------------------------------------------
@@ -450,4 +458,5 @@
 
 # if defined key_roam
+         !! extra MEDUSA-2 tracers
          DO jj = 2,jpjm1
             DO ji = 2,jpim1
@@ -456,7 +465,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 +479,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 +493,73 @@
          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 Node that this 
+         !!              Exceptionnal value did exist
+         !!
+         Call trc_bio_check(kt)
+
+         !!================================================================
+	 !! 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 nodes 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
+                     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 
+                     IF (lwp) CALL ctl_warn( 'trc_bio_medusa: carbonate chemistry WARNING:',  &
+                        TRIM(charout),TRIM(charout2),TRIM(charout3),TRIM(charout4),           & 
+                        'at i, j, k, kt:', TRIM(charout5) )
+                     !!
+                     !! Detect, report and correct tracer excursions
+                     IF ( (ztmp(ji,jj) .LT. -3.0) .OR. (ztmp(ji,jj) .GT.  40.0) )             &
+                        CALL trc_bio_exceptional_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_exceptional_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_exceptional_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_exceptional_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 +717,213 @@
    END SUBROUTINE trc_bio_medusa
 
+
+
+   SUBROUTINE trc_bio_exceptional_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(25), 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) WRITE(numout,*)                                 &
+         '----------------------------------------------------------------------'
+      IF(lwp) WRITE(numout,*)                                 &
+         'trc_bio_medusa: 3x3 neighbourhood surrounding abnormal ', TRIM(var_nm)
+      IF(lwp) WRITE(numout,9100)                              &
+         3, tiny_var(ii-1,ij+1), tiny_var(ii  ,ij+1), tiny_var(ii+1,ij+1)
+      IF(lwp) WRITE(numout,9100)                              &
+         2, tiny_var(ii-1,ij  ), tiny_var(ii  ,ij  ), tiny_var(ii+1,ij  )
+      IF(lwp) WRITE(numout,9100)                              &
+         1, tiny_var(ii-1,ij-1), tiny_var(ii  ,ij-1), tiny_var(ii+1,ij-1)
+      IF(lwp) WRITE(numout,*)                                 &
+         'trc_bio_medusa: 3x3 land-sea neighbourhood, tmask'
+      IF(lwp) WRITE(numout,9100)                              &
+         3, tiny_mask(ii-1,ij+1), tiny_mask(ii  ,ij+1), tiny_mask(ii+1,ij+1)
+      IF(lwp) WRITE(numout,9100)                              &
+         2, tiny_mask(ii-1,ij  ), tiny_mask(ii  ,ij  ), tiny_mask(ii+1,ij  )
+      IF(lwp) WRITE(numout,9100)                              &
+         1, tiny_mask(ii-1,ij-1), tiny_mask(ii  ,ij-1), tiny_mask(ii+1,ij-1)
+
+      !! 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
+      !!
+      write(charout1,9200) tiny_var(ii,ij)
+      write(charout2,9200) varout
+      write(charout3,*) ' ', charout1, ' -> ', charout2
+      write(charout4,*) ' Tracer: ', trim(var_nm)
+      IF(lwp) WRITE(numout,*) 'trc_bio_medusa: ** EXCEPTIONAL VALUE SWITCHING **'
+      IF(lwp) WRITE(numout,*) charout4 
+      IF(lwp) WRITE(numout,*) charout3
+      IF(lwp) WRITE(numout,*)                                 &
+         '----------------------------------------------------------------------'
+      IF(lwp) WRITE(numout,*)                                 &
+         ' '
+9100  FORMAT('Row:', i1, '  ', e16.6, ' ', e16.6, ' ', e16.6)
+9200  FORMAT(e16.6)
+
+   END SUBROUTINE trc_bio_exceptional_fix 
+
+   SUBROUTINE trc_bio_check(kt)
+      !!-----------------------------------
+      !! 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 node. 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
+      !!
+      !!==========================
+      !! DIC Check
+      zmax =  -5.e0  ! arbitrary  low maximum value
+      zmin =  4.0E4  ! arbitrary high minimum value
+      DO jj = 2, jpjm1
+         DO ji = 1, jpi
+            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 > 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.e0 )
+            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: sea surface DIC > 4000 '
+            WRITE(numout,9600) kt, zmax, ii, ij
+            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****'
+         ENDIF
+      ENDIF
+      !
+      IF( zmin <= 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.e0 )
+            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: sea surface DIC <= 0 '
+            WRITE(numout,9700) kt, zmin, ii, ij
+            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****'
+         ENDIF
+      ENDIF
+      !!
+      !!==========================
+      !! ALKALINITY Check
+      zmax =  -5.e0  ! arbitrary  low maximum value
+      zmin =  4.0E4  ! arbitrary high minimum value
+      DO jj = 2, jpjm1
+         DO ji = 1, jpi
+            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 > 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.e0 )
+            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: sea surface Alkalinity > 4000 '
+            WRITE(numout,9800) kt, zmax, ii, ij
+            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****'
+         ENDIF
+      ENDIF
+      !
+      IF( zmin <= 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.e0 )
+            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:  sea surface Alkalinity <= 0 '
+            WRITE(numout,9900) kt, zmin, ii, ij
+            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: ',2i5)
+9700  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' min DIC: ',f16.10,', i j: ',2i5)
+9800  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' max ALK: ',f16.10,', i j: ',2i5)
+9900  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' min ALK: ',f16.10,', i j: ',2i5)
+
+   END SUBROUTINE trc_bio_check
+
+
 #else
    !!=====================================================================
