Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 8521)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 9257)
@@ -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,           ONLY: trc_rst_dia_stat 
 
    !!* 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,17 @@
          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) > 1500.0 ).OR.(f_xco2a(ji,jj) < 100.0 ) ) THEN
+                    IF(lwp) THEN 
+                       WRITE(numout,*) ' atm xCO2 = ',f_xco2a(ji,jj),          &
+                                       ' -- ji =', mig(ji),' jj = ', mjg(jj)
+                    ENDIF
+                    CALL ctl_stop( 'MEDUSA - Air-Sea :', 'unrealistic atm xCO2 ' )
+                 ENDIF 
                ENDIF
                !!
@@ -162,19 +188,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
       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
@@ -201,8 +257,47 @@
                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
+              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 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   
 
