Index: /branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/CONFIG/SHARED/field_def_bgc.xml
===================================================================
--- /branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/CONFIG/SHARED/field_def_bgc.xml	(revision 7330)
+++ /branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/CONFIG/SHARED/field_def_bgc.xml	(revision 7331)
@@ -327,7 +327,7 @@
        <field id= "MED_XZE"    long_name="Euphotic depth"                            unit="m"           />
        <field id= "WIND"       long_name="Surface scalar wind"                       unit="m/s"         />
-       <field id= "ATM_PCO2"   long_name="Atmospheric pCO2"                          unit="ppmv"        />
+       <field id= "ATM_PCO2"   long_name="Atmospheric pCO2"                          unit="uatm"        />
        <field id= "OCN_PH"     long_name="Surface ocean pH"                          unit="-"           />
-       <field id= "OCN_PCO2"   long_name="Surface ocean pCO2"                        unit="ppmv"        />
+       <field id= "OCN_PCO2"   long_name="Surface ocean pCO2"                        unit="uatm"        />
        <field id= "OCNH2CO3"   long_name="Surface ocean H2CO3 conc"                  unit="mmolC/m3"    />
        <field id= "OCN_HCO3"   long_name="Surface ocean HCO3 conc"                   unit="mmolC/m3"    />
@@ -340,5 +340,5 @@
        <field id= "TALK"       long_name="Total alkalinity"                          unit="ueq/kg"      />
        <field id= "KW660"      long_name="Piston velocity"                           unit="m/s"         />
-       <field id= "ATM_PP0"    long_name="Atmospheric pressure"                      unit="bar"         />
+       <field id= "ATM_PP0"    long_name="Atmospheric pressure"                      unit="atm"         />
        <field id= "O2FLUX"     long_name="Air-sea O2 flux"                           unit="mmolO2/m2/d" />
        <field id= "O2SAT"      long_name="Surface ocean O2 saturation"               unit="mmolO2/m3"   />
@@ -417,4 +417,14 @@
        <field id= "DMS_ARAN"   long_name="DMS Surface Concentration, AranamiT"       unit="nmol/L"      />
        <field id= "DMS_HALL"   long_name="DMS Surface Concentration, Halloran"       unit="nmol/L"      />
+       <field id= "ATM_XCO2"   long_name="Atmospheric xCO2"                          unit="ppm"         />
+       <field id= "OCN_FCO2"   long_name="Surface ocean fCO2"                        unit="uatm"        />
+       <field id= "ATM_FCO2"   long_name="Atmospheric fCO2"                          unit="uatm"        />
+       <field id= "OCN_RHOSW"  long_name="MOCSY rhoSW"                               unit="kg/m3"       />
+       <field id= "OCN_SCHCO2" long_name="MOCSY Schmidt CO2"                         unit="-"           />
+       <field id= "OCN_KWCO2"  long_name="MOCSY kw CO2"                              unit="m/s"         />
+       <field id= "OCN_K0"     long_name="MOCSY k0"                                  unit="(mol/kg)/atm" />
+       <field id= "CO2STARAIR" long_name="MOCSY CO2star air"                         unit="mmol/m3"     />
+       <field id= "OCN_DPCO2"  long_name="MOCSY delta pCO2"                          unit="uatm"        />
+       
 
        <!-- AXY (08/11/16): add new 2D CMIP6 diagnostics -->
@@ -431,7 +441,8 @@
        <field id= "FBDDTALK"   long_name="Int. change of Alkalinity due to Biology"  unit="meq/m2/d"    />
        <field id= "FBDDTDIC"   long_name="Int. change of DI Carbon due to Biology"   unit="mmolC/m2/d"  />
-       <field id= "FBDDTDIFE"  long_name="Int. change of DI Iron due to Biology"     unit="mmolFe/m2/d" />
+       <field id= "FBDDTIFE"   long_name="Int. change of DI Iron due to Biology"     unit="mmolFe/m2/d" />
        <field id= "FBDDTDIN"   long_name="Int. change of DI Nitrogen due to Biology" unit="mmolN/m2/d"  />
        <field id= "FBDDTDISI"  long_name="Int. change of DI Silicon due to Biology"  unit="mmolSi/m2/d" />
+       
 
        <!-- MEDUSA 3D diagnostics -->
@@ -709,4 +720,13 @@
       <field field_ref= "INT_DET"    name="INT_DET"    />
       <field field_ref= "INT_DTC"    name="INT_DTC"    />
+      <field field_ref= "ATM_XCO2"   name="ATM_XCO2"   />
+      <field field_ref= "OCN_FCO2"   name="OCN_FCO2"   />
+      <field field_ref= "ATM_FCO2"   name="ATM_FCO2"   />
+      <field field_ref= "OCN_RHOSW"  name="OCN_RHOSW"  />
+      <field field_ref= "OCN_SCHCO2" name="OCN_SCHCO2" />
+      <field field_ref= "OCN_KWCO2"  name="OCN_KWCO2"  />
+      <field field_ref= "OCN_K0"     name="OCN_K0"     />
+      <field field_ref= "CO2STARAIR" name="CO2STARAIR" />
+      <field field_ref= "OCN_DPCO2"  name="OCN_DPCO2"  />
     </field_group>
 
Index: /branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
===================================================================
--- /branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 7330)
+++ /branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 7331)
@@ -175,4 +175,6 @@
       REAL(wp)                     ::    fthetan,faln,fchn1,fchn,fjln,fprn,frn
       REAL(wp)                     ::    fthetad,fald,fchd1,fchd,fjld,fprd,frd
+      !! AXY (23/11/16): add in light-only limitation term (normalised 0-1 range)
+      REAL(wp)                     ::    fjlim_pn, fjlim_pd
       !! AXY (03/02/11): add in Liebig terms
       REAL(wp) ::    fpnlim, fpdlim
@@ -296,5 +298,7 @@
       !!
       !! ROAM air-sea flux and diagnostic parameters
-      REAL(wp) ::    f_wind, f_pco2a
+      REAL(wp) ::    f_wind
+      !! AXY (24/11/16): add xCO2 variable for atmosphere (what we actually have)
+      REAL(wp) ::    f_xco2a
       REAL(wp) ::    f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_co2flux
       REAL(wp) ::    f_TDIC, f_TALK, f_dcf, f_henry
@@ -372,4 +376,8 @@
       REAL(wp), POINTER, DIMENSION(:,:  ) :: oben_fe2d, oben_c2d, oben_si2d, oben_ca2d, sfr_ocal2d
       REAL(wp), POINTER, DIMENSION(:,:  ) :: sfr_oarg2d, lyso_ca2d 
+      !! AXY (23/11/16): extra MOCSY diagnostics
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: f_xco2a_2d, f_fco2w_2d, f_fco2a_2d
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: f_ocnrhosw_2d, f_ocnschco2_2d, f_ocnkwco2_2d
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: f_ocnk0_2d, f_co2starair_2d, f_ocndpco2_2d
 # endif
       !!
@@ -986,4 +994,42 @@
             ENDIF
          ENDIF   
+         !!
+         !! AXY (24/11/16): extra MOCSY diagnostics, 2D
+         IF( med_diag%ATM_XCO2%dgsave ) THEN
+            CALL wrk_alloc( jpi, jpj, f_xco2a_2d      )
+            f_xco2a_2d(:,:)      = 0.0 !!
+         ENDIF
+         IF( med_diag%OCN_FCO2%dgsave ) THEN
+            CALL wrk_alloc( jpi, jpj, f_fco2w_2d      )
+            f_fco2w_2d(:,:)      = 0.0 !!
+         ENDIF
+         IF( med_diag%ATM_FCO2%dgsave ) THEN
+            CALL wrk_alloc( jpi, jpj, f_fco2a_2d      )
+            f_fco2a_2d(:,:)      = 0.0 !!
+         ENDIF
+         IF( med_diag%OCN_RHOSW%dgsave ) THEN
+            CALL wrk_alloc( jpi, jpj, f_ocnrhosw_2d   )
+            f_ocnrhosw_2d(:,:)      = 0.0 !!
+         ENDIF
+         IF( med_diag%OCN_SCHCO2%dgsave ) THEN
+            CALL wrk_alloc( jpi, jpj, f_ocnschco2_2d  )
+            f_ocnschco2_2d(:,:)      = 0.0 !!
+         ENDIF
+         IF( med_diag%OCN_KWCO2%dgsave ) THEN
+            CALL wrk_alloc( jpi, jpj, f_ocnkwco2_2d   )
+            f_ocnkwco2_2d(:,:)      = 0.0 !!
+         ENDIF
+         IF( med_diag%OCN_K0%dgsave ) THEN
+            CALL wrk_alloc( jpi, jpj, f_ocnk0_2d      )
+            f_ocnk0_2d(:,:)      = 0.0 !!
+         ENDIF
+         IF( med_diag%CO2STARAIR%dgsave ) THEN
+            CALL wrk_alloc( jpi, jpj, f_co2starair_2d )
+            f_co2starair_2d(:,:)      = 0.0 !!
+         ENDIF
+         IF( med_diag%OCN_DPCO2%dgsave ) THEN
+            CALL wrk_alloc( jpi, jpj, f_ocndpco2_2d   )
+            f_ocndpco2_2d(:,:)      = 0.0 !!
+         ENDIF
 # endif  
          IF( med_diag%TPP3%dgsave ) THEN
@@ -1025,5 +1071,5 @@
           IF( med_diag%O2min%dgsave ) THEN
              CALL wrk_alloc( jpi, jpj,   o2min    )
-             o2min(:,:)  = 1.e6 !! set to high value as we're looking for min(o2)
+             o2min(:,:)  = 1.e3 !! set to high value as we're looking for min(o2)
           ENDIF          
           IF( med_diag%ZO2min%dgsave ) THEN
@@ -1227,8 +1273,8 @@
       if (iyr1 .le. 1) then
          !! before 1860
-         f_pco2a = hist_pco2(1)
+         f_xco2a = hist_pco2(1)
       elseif (iyr2 .ge. 242) then
          !! after 2099
-         f_pco2a = hist_pco2(242)
+         f_xco2a = hist_pco2(242)
       else
          !! just right
@@ -1243,9 +1289,8 @@
 #  endif
          fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3)
-         f_pco2a = fq4
+         f_xco2a = fq4
       endif
 #  if defined key_axy_pi_co2
-      !! f_pco2a = hist_pco2(1)
-      f_pco2a = 284.725          !! OCMIP pre-industrial pCO2
+      f_xco2a = 284.725          !! OCMIP pre-industrial pCO2
 #  endif
       !! IF(lwp) WRITE(numout,*) ' MEDUSA nyear     =', nyear
@@ -1257,5 +1302,5 @@
       !! IF(lwp) WRITE(numout,*) ' MEDUSA fq2       =', fq2
       !! IF(lwp) WRITE(numout,*) ' MEDUSA fq3       =', fq3
-      IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2  =', f_pco2a
+      IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2  =', f_xco2a
 # endif
 
@@ -1301,5 +1346,5 @@
                   if (tmask(ji,jj,jk).eq.1) then
                      IF (lk_oasis) THEN
-                        f_pco2a = PCO2a_in_cpl(ji,jj)        !! use 2D atm pCO2 from atm coupling
+                        f_xco2a = PCO2a_in_cpl(ji,jj)        !! use 2D atm xCO2 from atm coupling
                      ENDIF
                      !! do carbonate chemistry
@@ -1341,5 +1386,5 @@
                      !!                 chemistry package
                      CALL mocsy_interface( ztmp, zsal, zalk, zdic, zsil, zpho,         &    ! inputs
-                     f_pp0, fdep2, flatx, f_kw660, f_pco2a, 1,                         &    ! inputs
+                     f_pp0, fdep2, flatx, f_kw660, f_xco2a, 1,                         &    ! inputs
                      f_ph, f_pco2w, f_fco2w, f_h2co3, f_hco3, f_co3, f_omarg(ji,jj),   &    ! outputs
                      f_omcal(ji,jj), f_BetaD, f_rhosw, f_opres, f_insitut,             &    ! outputs
@@ -1354,5 +1399,5 @@
                      !!                 MEDUSA-2 default)
                      CALL trc_co2_medusa( ztmp, zsal, zdic, zalk, fdep2, f_kw660,      &    ! inputs
-                     f_pco2a, f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_omcal(ji,jj),   &    ! outputs
+                     f_xco2a, f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_omcal(ji,jj),   &    ! outputs
                      f_omarg(ji,jj), f_co2flux, f_TDIC, f_TALK, f_dcf, f_henry, iters)      ! outputs
                      !! 
@@ -1655,5 +1700,5 @@
                   f_wind  = wndm(ji,jj)
                   IF (lk_oasis) THEN
-                     f_pco2a = PCO2a_in_cpl(ji,jj)        !! use 2D atm pCO2 from atm coupling
+                     f_xco2a = PCO2a_in_cpl(ji,jj)        !! use 2D atm xCO2 from atm coupling
                   ENDIF
                   !!
@@ -1696,5 +1741,5 @@
                   !!                 zero in this call
                   CALL mocsy_interface( ztmp, zsal, zalk, zdic, zsil, zpho,        &  ! inputs
-                  f_pp0, 0.0, flatx, f_kw660, f_pco2a, 1,                          &  ! inputs
+                  f_pp0, 0.0, flatx, f_kw660, f_xco2a, 1,                          &  ! inputs
                   f_ph, f_pco2w, f_fco2w, f_h2co3, f_hco3, f_co3, f_omarg(ji,jj),  &  ! outputs
                   f_omcal(ji,jj), f_BetaD, f_rhosw, f_opres, f_insitut,            &  ! outputs
@@ -1709,5 +1754,5 @@
                   !!
                   !! carbon dioxide (CO2); Jerry Blackford code (ostensibly OCMIP-2, but not)
-                  CALL trc_co2_medusa( ztmp, zsal, zdic, zalk, 0.0, f_kw660, f_pco2a,  &  ! inputs
+                  CALL trc_co2_medusa( ztmp, zsal, zdic, zalk, 0.0, f_kw660, f_xco2a,  &  ! inputs
                   f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_omcal(ji,jj),               &  ! outputs
                   f_omarg(ji,jj), f_co2flux, f_TDIC, f_TALK, f_dcf, f_henry, iters )      ! outputs
@@ -1728,5 +1773,5 @@
                   !!                 quasi-sensible alternatives
                   f_ph           = 8.1
-                  f_pco2w        = f_pco2a
+                  f_pco2w        = f_xco2a
                   f_h2co3        = 0.005 * zdic
                   f_hco3         = 0.865 * zdic
@@ -1740,15 +1785,15 @@
                   f_henry        = 1.
                   !! AXY (23/06/15): add in some extra MOCSY diagnostics
-                  f_fco2w        = fpco2a
+                  f_fco2w        = f_xco2a
                   f_BetaD        = 1.
                   f_rhosw        = 1.026
                   f_opres        = 0.
                   f_insitut      = ztmp
-                  f_pco2atm      = fpco2a
-                  f_fco2atm      = fpco2a
+                  f_pco2atm      = f_xco2a
+                  f_fco2atm      = f_xco2a
                   f_schmidtco2   = 660.
                   f_kwco2        = 0.
                   f_K0           = 0.
-                  f_co2starair   = fpco2a
+                  f_co2starair   = f_xco2a
                   f_dpco2        = 0.
 #  endif
@@ -1842,5 +1887,5 @@
                   IF ( lk_iomput ) THEN
                       IF( med_diag%ATM_PCO2%dgsave ) THEN
-                         f_pco2a2d(ji,jj) = f_pco2a
+                         f_pco2a2d(ji,jj) = f_pco2atm
                       ENDIF
                       IF( med_diag%OCN_PCO2%dgsave ) THEN
@@ -1868,4 +1913,32 @@
                          f_o2sat2d(ji,jj) = f_o2sat
                       ENDIF
+                      !! AXY (24/11/16): add in extra MOCSY diagnostics
+                      IF( med_diag%ATM_XCO2%dgsave ) THEN
+                         f_xco2a_2d(ji,jj) = f_xco2a
+                      ENDIF
+                      IF( med_diag%OCN_FCO2%dgsave ) THEN
+                         f_fco2w_2d(ji,jj) = f_fco2w
+                      ENDIF
+                      IF( med_diag%ATM_FCO2%dgsave ) THEN
+                         f_fco2a_2d(ji,jj) = f_fco2atm
+                      ENDIF
+                      IF( med_diag%OCN_RHOSW%dgsave ) THEN
+                         f_ocnrhosw_2d(ji,jj) = f_rhosw
+                      ENDIF
+                      IF( med_diag%OCN_SCHCO2%dgsave ) THEN
+                         f_ocnschco2_2d(ji,jj) = f_schmidtco2
+                      ENDIF
+                      IF( med_diag%OCN_KWCO2%dgsave ) THEN
+                         f_ocnkwco2_2d(ji,jj) = f_kwco2
+                      ENDIF
+                      IF( med_diag%OCN_K0%dgsave ) THEN
+                         f_ocnk0_2d(ji,jj) = f_K0
+                      ENDIF
+                      IF( med_diag%CO2STARAIR%dgsave ) THEN
+                         f_co2starair_2d(ji,jj) = f_co2starair
+                      ENDIF
+                      IF( med_diag%OCN_DPCO2%dgsave ) THEN
+                         f_ocndpco2_2d(ji,jj) = f_dpco2
+                      ENDIF
                   ENDIF
                   !! 
@@ -2029,4 +2102,5 @@
                endif
                fjln    = fchn * faln * xpar(ji,jj,jk) !! non-diatom J term
+               fjlim_pn = fjln / xvpnT
                !!
                !! diatoms
@@ -2038,4 +2112,5 @@
                endif
                fjld    = fchd * fald * xpar(ji,jj,jk) !! diatom J term
+               fjlim_pd = fjld / xvpdT
       
 # if defined key_debug_medusa
@@ -4014,5 +4089,6 @@
                   ENDIF
                   IF( med_diag%PN_JLIM%dgsave ) THEN
-                      fjln2d(ji,jj) = fjln2d(ji,jj) + (fjln  * zphn * fthk) 
+                      ! fjln2d(ji,jj) = fjln2d(ji,jj) + (fjln  * zphn * fthk) 
+                      fjln2d(ji,jj) = fjln2d(ji,jj) + (fjlim_pn * zphn * fthk) 
                   ENDIF
                   IF( med_diag%PN_NLIM%dgsave ) THEN
@@ -4023,5 +4099,6 @@
                   ENDIF
                   IF( med_diag%PD_JLIM%dgsave ) THEN
-                      fjld2d(ji,jj) = fjld2d(ji,jj) + (fjld  * zphd * fthk) 
+                      ! fjld2d(ji,jj) = fjld2d(ji,jj) + (fjld  * zphd * fthk) 
+                      fjld2d(ji,jj) = fjld2d(ji,jj) + (fjlim_pd * zphd * fthk) 
                   ENDIF
                   IF( med_diag%PD_NLIM%dgsave ) THEN
@@ -4269,5 +4346,5 @@
                   ENDIF
 		  IF ( med_diag%FEDISS3%dgsave  ) THEN
-                     fescav3(ji,jj,jk) = ffetop
+                     fediss3(ji,jj,jk) = ffetop
                   ENDIF
 		  IF ( med_diag%FESCAV3%dgsave  ) THEN
@@ -4299,5 +4376,6 @@
                   ENDIF
 		  IF ( med_diag%PNLIMJ3%dgsave  ) THEN
-                     pnlimj3(ji,jj,jk)  = fjln
+                     ! pnlimj3(ji,jj,jk)  = fjln
+                     pnlimj3(ji,jj,jk)  = fjlim_pn
                   ENDIF
 		  IF ( med_diag%PNLIMN3%dgsave  ) THEN
@@ -4308,5 +4386,6 @@
                   ENDIF
 		  IF ( med_diag%PDLIMJ3%dgsave  ) THEN
-                     pdlimj3(ji,jj,jk)  = fjld
+                     ! pdlimj3(ji,jj,jk)  = fjld
+                     pdlimj3(ji,jj,jk)  = fjlim_pd
                   ENDIF
 		  IF ( med_diag%PDLIMN3%dgsave  ) THEN
@@ -4357,8 +4436,8 @@
                   trc2d(ji,jj,23) = trc2d(ji,jj,23) + (ffebot  * fthk)         !! seafloor iron addition
                   trc2d(ji,jj,24) = trc2d(ji,jj,24) + (ffescav * fthk)         !! "free" iron scavenging
-                  trc2d(ji,jj,25) = trc2d(ji,jj,25) + (fjln  * zphn * fthk)    !! non-diatom J  limitation term 
+                  trc2d(ji,jj,25) = trc2d(ji,jj,25) + (fjlim_pn * zphn * fthk) !! non-diatom J  limitation term 
                   trc2d(ji,jj,26) = trc2d(ji,jj,26) + (fnln  * zphn * fthk)    !! non-diatom N  limitation term 
                   trc2d(ji,jj,27) = trc2d(ji,jj,27) + (ffln  * zphn * fthk)    !! non-diatom Fe limitation term 
-                  trc2d(ji,jj,28) = trc2d(ji,jj,28) + (fjld  * zphd * fthk)    !! diatom     J  limitation term 
+                  trc2d(ji,jj,28) = trc2d(ji,jj,28) + (fjlim_pd * zphd * fthk) !! diatom     J  limitation term 
                   trc2d(ji,jj,29) = trc2d(ji,jj,29) + (fnld  * zphd * fthk)    !! diatom     N  limitation term 
                   trc2d(ji,jj,30) = trc2d(ji,jj,30) + (ffld  * zphd * fthk)    !! diatom     Fe limitation term 
@@ -4432,5 +4511,5 @@
                   if (jk .eq. 1) then
                      trc2d(ji,jj,91)  = trc2d(ji,jj,91)  + f_wind              !! surface wind
-                     trc2d(ji,jj,92)  = trc2d(ji,jj,92)  + f_pco2a             !! atmospheric pCO2
+                     trc2d(ji,jj,92)  = trc2d(ji,jj,92)  + f_pco2atm           !! atmospheric pCO2
                      trc2d(ji,jj,93)  = trc2d(ji,jj,93)  + f_ph                !! ocean pH
                      trc2d(ji,jj,94)  = trc2d(ji,jj,94)  + f_pco2w             !! ocean pCO2
@@ -4821,4 +4900,41 @@
                        ENDIF
                      ENDIF
+                     !! AXY (24/11/16): extra MOCSY diagnostics
+                     IF( med_diag%ATM_XCO2%dgsave ) THEN
+                        CALL iom_put( "ATM_XCO2"  ,   f_xco2a_2d      )
+                        CALL wrk_dealloc( jpi, jpj,   f_xco2a_2d      )
+                     ENDIF
+                     IF( med_diag%OCN_FCO2%dgsave ) THEN
+                        CALL iom_put( "OCN_FCO2"  ,   f_fco2w_2d      )
+                        CALL wrk_dealloc( jpi, jpj,   f_fco2w_2d      )
+                     ENDIF
+                     IF( med_diag%ATM_FCO2%dgsave ) THEN
+                        CALL iom_put( "ATM_FCO2"  ,   f_fco2a_2d      )
+                        CALL wrk_dealloc( jpi, jpj,   f_fco2a_2d      )
+                     ENDIF
+                     IF( med_diag%OCN_RHOSW%dgsave ) THEN
+                        CALL iom_put( "OCN_RHOSW"  ,  f_ocnrhosw_2d   )
+                        CALL wrk_dealloc( jpi, jpj,   f_ocnrhosw_2d   )
+                     ENDIF
+                     IF( med_diag%OCN_SCHCO2%dgsave ) THEN
+                        CALL iom_put( "OCN_SCHCO2"  , f_ocnschco2_2d  )
+                        CALL wrk_dealloc( jpi, jpj,   f_ocnschco2_2d  )
+                     ENDIF
+                     IF( med_diag%OCN_KWCO2%dgsave ) THEN
+                        CALL iom_put( "OCN_KWCO2"  ,  f_ocnkwco2_2d   )
+                        CALL wrk_dealloc( jpi, jpj,   f_ocnkwco2_2d   )
+                     ENDIF
+                     IF( med_diag%OCN_K0%dgsave ) THEN
+                        CALL iom_put( "OCN_K0"  ,     f_ocnk0_2d      )
+                        CALL wrk_dealloc( jpi, jpj,   f_ocnk0_2d      )
+                     ENDIF
+                     IF( med_diag%CO2STARAIR%dgsave ) THEN
+                        CALL iom_put( "CO2STARAIR"  , f_co2starair_2d )
+                        CALL wrk_dealloc( jpi, jpj,   f_co2starair_2d )
+                     ENDIF
+                     IF( med_diag%OCN_DPCO2%dgsave ) THEN
+                        CALL iom_put( "OCN_DPCO2"  ,  f_ocndpco2_2d   )
+                        CALL wrk_dealloc( jpi, jpj,   f_ocndpco2_2d   )
+                     ENDIF
 # endif                     
                  ELSE IF (jk.eq.i0100) THEN 
@@ -5851,4 +5967,8 @@
           !!
           !! AXY (09/11/16): 3D CMIP6 diagnostics
+          IF( med_diag%TPPD3%dgsave ) THEN
+             CALL iom_put( "TPPD3"     , tppd3 )
+             CALL wrk_dealloc( jpi, jpj, jpk, tppd3      )
+          ENDIF          
           IF( med_diag%BDDTALK3%dgsave ) THEN
              CALL iom_put( "BDDTALK3"  , bddtalk3 )
Index: /branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90
===================================================================
--- /branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 7330)
+++ /branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 7331)
@@ -2098,4 +2098,50 @@
           med_diag%DMS_HALL%dgsave = .FALSE.
       ENDIF
+      IF  (iom_use("ATM_XCO2")) THEN 
+          med_diag%ATM_XCO2%dgsave = .TRUE.
+      ELSE 
+          med_diag%ATM_XCO2%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("OCN_FCO2")) THEN 
+          med_diag%OCN_FCO2%dgsave = .TRUE.
+      ELSE 
+          med_diag%OCN_FCO2%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("ATM_FCO2")) THEN 
+          med_diag%ATM_FCO2%dgsave = .TRUE.
+      ELSE 
+          med_diag%ATM_FCO2%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("OCN_RHOSW")) THEN 
+          med_diag%OCN_RHOSW%dgsave = .TRUE.
+      ELSE 
+          med_diag%OCN_RHOSW%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("OCN_SCHCO2")) THEN 
+          med_diag%OCN_SCHCO2%dgsave = .TRUE.
+      ELSE 
+          med_diag%OCN_SCHCO2%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("OCN_KWCO2")) THEN 
+          med_diag%OCN_KWCO2%dgsave = .TRUE.
+      ELSE 
+          med_diag%OCN_KWCO2%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("OCN_K0")) THEN 
+          med_diag%OCN_K0%dgsave = .TRUE.
+      ELSE 
+          med_diag%OCN_K0%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("CO2STARAIR")) THEN 
+          med_diag%CO2STARAIR%dgsave = .TRUE.
+      ELSE 
+          med_diag%CO2STARAIR%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("OCN_DPCO2")) THEN 
+          med_diag%OCN_DPCO2%dgsave = .TRUE.
+      ELSE 
+          med_diag%OCN_DPCO2%dgsave = .FALSE.
+      ENDIF
+      !!
       IF  (iom_use("TPP3")) THEN 
           med_diag%TPP3%dgsave = .TRUE.
@@ -2199,4 +2245,9 @@
           med_diag%FBDDTDIFE%dgsave = .FALSE.
       ENDIF
+      IF  (iom_use("FBDDTDIN")) THEN 
+          med_diag%FBDDTDIN%dgsave = .TRUE.
+      ELSE 
+          med_diag%FBDDTDIN%dgsave = .FALSE.
+      ENDIF
       IF  (iom_use("FBDDTDISI")) THEN 
           med_diag%FBDDTDISI%dgsave = .TRUE.
@@ -2225,4 +2276,9 @@
       ELSE 
           med_diag%BDDTDIFE3%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("BDDTDIN3")) THEN 
+          med_diag%BDDTDIN3%dgsave = .TRUE.
+      ELSE 
+          med_diag%BDDTDIN3%dgsave = .FALSE.
       ENDIF
       IF  (iom_use("BDDTDISI3")) THEN 
Index: /branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/trc.F90
===================================================================
--- /branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 7330)
+++ /branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 7331)
@@ -131,5 +131,6 @@
                   ZE_MES_C, ZE_MESDC, ZE_EXCR, ZE_RESP, ZE_GROW, MDETC, GMIDC, GMEDC,                &
                   INT_ZMI, INT_ZME, INT_DET, INT_DTC, DMS_SURF, DMS_ANDR, DMS_SIMO, DMS_ARAN,        &
-                  DMS_HALL,                                                                          & ! end of regular 2D
+                  DMS_HALL, ATM_XCO2, OCN_FCO2, ATM_FCO2, OCN_RHOSW, OCN_SCHCO2, OCN_KWCO2,          &
+                  OCN_K0, CO2STARAIR, OCN_DPCO2,                                                     & ! end of regular 2D
                   TPP3, DETFLUX3, REMIN3N, PH3, OM_CAL3,                                             & ! end of regular 3D
 ! AXY (11/11/16): additional CMIP6 2D diagnostics
