Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 9257)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 9385)
@@ -81,5 +81,7 @@
                                             gdept_0, gdept_n,               &
                                             gdepw_0, gdepw_n,               &
-                                            nday_year, nsec_day, nyear,     &
+                                            nday_year, nsec_day,            &
+                                            nyear, nyear_len, ndastp,       &
+                                            nsec_month,                     &
                                             rdt, tmask, mig, mjg, nimpp,    &
                                             njmpp 
@@ -92,5 +94,5 @@
                                             mpp_min, mpp_minloc,            &
                                             ctl_stop, ctl_warn, lk_mpp  
-      USE oce,                        ONLY: tsb, tsn
+      USE oce,                        ONLY: tsb, tsn, PCO2a_in_cpl
       USE par_kind,                   ONLY: wp
       USE par_medusa,                 ONLY: jpalk, jpchd, jpchn, jpdet,     &
@@ -102,5 +104,6 @@
       !! 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, &
+                                            lk_pi_co2
       USE trc,                        ONLY: ln_rsttr, nittrc000, trn
       USE bio_medusa_init_mod,        ONLY: bio_medusa_init
@@ -114,4 +117,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
@@ -181,5 +185,6 @@
       !! 
       !! 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
@@ -293,44 +298,85 @@
       !!------------------------------------------------------------------
       !!
-      !! 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 
+      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
 
@@ -358,6 +404,23 @@
       !!          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 .
+      IF ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        &
+           ( nsec_month .LE. INT(rdt) ) )  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
@@ -502,5 +565,5 @@
          !!              Exceptionnal value did exist
          !!
-         Call trc_bio_check(kt)
+         Call trc_bio_check(kt, jk)
 
          !!================================================================
@@ -810,5 +873,5 @@
    END SUBROUTINE trc_bio_exceptionnal_fix 
 
-   SUBROUTINE trc_bio_check(kt)
+   SUBROUTINE trc_bio_check(kt, jk)
       !!-----------------------------------
       !! JPALM -- 14-12-2017 -- Still dealing with this micro-boil/carb failure
@@ -824,5 +887,5 @@
       INTEGER               :: ii,ij         ! temporary scalars 
       INTEGER, DIMENSION(2) :: ilocs         ! 
-      INTEGER, INTENT( in ) :: kt
+      INTEGER, INTENT( in ) :: kt, jk
       !!
       !!==========================
@@ -852,6 +915,6 @@
          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: DIC concentration > 4000 '
+            WRITE(numout,9600) kt, zmax, ii, ij, jk
             WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****'
          ENDIF
@@ -869,6 +932,6 @@
          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: DIC concentration <= 0 '
+            WRITE(numout,9700) kt, zmin, ii, ij, jk
             WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****'
          ENDIF
@@ -901,6 +964,6 @@
          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: ALK concentration > 4000 '
+            WRITE(numout,9800) kt, zmax, ii, ij, jk
             WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****'
          ENDIF
@@ -918,6 +981,6 @@
          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:  ALK concentration <= 0 '
+            WRITE(numout,9900) kt, zmin, ii, ij, jk
             WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****'
          ENDIF
@@ -925,8 +988,8 @@
 
 
-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)
+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
