Index: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 10028)
+++ branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 10045)
@@ -7,4 +7,5 @@
    !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
    !!   -   ! 2017-08 (A. Yool)            Add air-sea flux kill switch
+   !!   -   ! 2018-08 (A. Yool)            add OMIP preindustrial DIC
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -57,4 +58,18 @@
                                    f_xco2a, f_xco2a_2d,                   &
                                    zalk, zdic, zoxy, zsal, ztmp,          &
+#  if defined key_omip_dic
+                                   pi_fgco2,                               &
+                                   f_pi_co2flux, f_pi_co2flux2d,           &
+                                   f_pi_co2starair_2d, f_pi_co3,           &
+                                   f_pi_fco2a_2d, f_pi_fco2w_2d,           &
+                                   f_pi_h2co3, f_pi_hco3,                  &
+                                   f_pi_ocndpco2_2d,                       &
+                                   f_pi_omarg, f_pi_omcal,                 &
+                                   f_pi_pco2a2d, f_pi_pco2atm, f_pi_pco2w, &
+                                   f_pi_pco2w2d, f_pi_ph,                  &
+                                   f_pi_TDIC, f_pi_TDIC2d,                 &
+                                   f_pi_xco2a, f_pi_xco2a_2d,              &
+                                   zomd,                                   &
+#  endif                                   
 # endif
 # if defined key_mocsy
@@ -119,6 +134,9 @@
       REAL(wp) :: f_insitut_dum
       REAL(wp) :: f_kwo2_dum
+#  if defined key_omip_dic
+      !! AXY (06/08/18): output variables that are not different for OMIP DIC
+      REAL(wp) :: f_rhosw_dum, f_schmidtco2_dum, f_kwco2_dum,f_K0_dum
+#  endif
 # endif
-
 
 # if defined key_roam
@@ -132,4 +150,10 @@
       f_rhosw(:,:)       = 0.0
       f_K0(:,:)          = 0.0
+#  if defined key_omip_dic
+      f_pi_fco2w(:,:)       = 0.0
+      f_pi_fco2atm(:,:)     = 0.0
+      f_pi_co2starair(:,:)  = 0.0
+      f_pi_dpco2(:,:)       = 0.0      
+#  endif
       !! air pressure (atm); ultimately this will use air 
       !! pressure at the base of the UKESM1 atmosphere 
@@ -137,14 +161,13 @@
       f_pp0(:,:)   = 1.0
 
-
       !!-----------------------------------------------------------
       !! Air-sea gas exchange
       !!-----------------------------------------------------------
 
-#   if defined key_debug_medusa
+#  if defined key_debug_medusa
                IF (lwp) write (numout,*)                     & 
                'air-sea: gas_transfer kt = ', kt
                CALL flush(numout)
-#   endif
+#  endif
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
@@ -170,25 +193,29 @@
       ENDDO
 
-#   if defined key_debug_medusa
-               IF (lwp) write (numout,*)                     &
-               '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
+#  if defined key_debug_medusa
+      IF (lwp) write (numout,*)                     &
+           '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
+
+#  if defined key_axy_carbchem
+#   if defined key_mocsy
+      !!-----------------------------------------------------------
+      !! MOCSY carbonate chemistry
+      !!-----------------------------------------------------------
+      !!
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
             if (tmask(ji,jj,1) == 1) then
-               !!
-#  if defined key_axy_carbchem
-#   if defined key_mocsy
                !! Jpalm -- 12-09-2017 -- add extra check after reccurent
                !!          carbonate failure in the coupled run.
@@ -201,22 +228,22 @@
                     (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     
+                     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     
+               ENDIF
                !!
                !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate
@@ -281,21 +308,56 @@
                !!                 'unrealistic surface Carb. Chemistry -- OUTPUTS' )
                !!ENDIF     
+#    if defined key_omip_dic
+               !! AXY (06/08/18): repeat carbonate chemistry calculations but
+               !!                 using preindustrial DIC for OMIP; note that
+               !!                 some outputs below are dummy since they are
+               !!                 the same (or should be!) as for regular DIC
+               CALL mocsy_interface(ztmp(ji,jj),zsal(ji,jj),zalk(ji,jj),        &
+                                    zomd(ji,jj),zsil(ji,jj),zpho(ji,jj),        &
+                                    f_pp0(ji,jj),0.0,                           &
+                                    gphit(ji,jj),f_kw660(ji,jj),                &
+                                    f_pi_xco2a(ji,jj),1,f_pi_ph(ji,jj),         &
+                                    f_pi_pco2w(ji,jj),f_pi_fco2w(ji,jj),        &
+                                    f_pi_h2co3(ji,jj),f_pi_hco3(ji,jj),         &
+                                    f_pi_co3(ji,jj),f_pi_omarg(ji,jj),          &
+                                    f_pi_omcal(ji,jj),f_BetaD_dum,              &
+                                    f_rhosw_dum,f_opres_dum,                    &
+                                    f_insitut_dum,f_pi_pco2atm(ji,jj),          &
+                                    f_pi_fco2atm(ji,jj),f_schmidtco2_dum,       &
+                                    f_kwco2_dum,f_K0_dum,                       &
+                                    f_pi_co2starair(ji,jj),f_pi_co2flux(ji,jj), &
+                                    f_pi_dpco2(ji,jj))
+               !! mmol / m3 -> umol / kg
+               f_pi_TDIC(ji,jj) = (zomd(ji,jj) / f_rhosw(ji,jj)) * 1000.
+#    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
+#    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')
+#     if defined key_omip_dic
+      call trc_rst_dia_stat(f_pi_pco2w(:,:), 'f_pi_pco2w')
+      call trc_rst_dia_stat(f_pi_fco2w(:,:), 'f_pi_fco2w')
+      call trc_rst_dia_stat(f_pi_fco2atm(:,:), 'f_pi_fco2atm')
+      call trc_rst_dia_stat(f_pi_co2starair(:,:), 'f_pi_co2starair')
+      call trc_rst_dia_stat(f_pi_co2flux(:,:), 'f_pi_co2flux')
+      call trc_rst_dia_stat(f_pi_dpco2(:,:), 'f_pi_dpco2')      
+#     endif
+#    endif
+               
 #   else   
-
+      !!-----------------------------------------------------------
+      !! Blackford et al. (2007) carbonate chemistry
+      !!-----------------------------------------------------------
+      !!
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
@@ -304,5 +366,5 @@
                !!
                !! carbon dioxide (CO2); Jerry Blackford code (ostensibly 
-               !! OCMIP-2, but not)
+               !! OCMIP-2, but not); this should be removed in the future
                CALL trc_co2_medusa(ztmp(ji,jj),zsal(ji,jj),zdic(ji,jj),      &
                                    zalk(ji,jj),0.0,                          &
@@ -324,5 +386,5 @@
                if (iters .eq. 25) then
                   IF(lwp) WRITE(numout,*) 'air-sea: ITERS WARNING, ',       &
-                     iters, ' AT (', ji, ', ', jj, ', 1) AT ', kt
+                       iters, ' AT (', ji, ', ', jj, ', 1) AT ', kt
                endif
             ENDIF
@@ -332,5 +394,8 @@
 #   endif
 #  else
-
+      !!-----------------------------------------------------------
+      !! No carbonate chemistry
+      !!-----------------------------------------------------------
+      !!
       DO jj = 2,jpjm1
          DO ji = 2,jpim1
@@ -374,4 +439,8 @@
       !! AXY (18/08/17): single kill switch on air-sea CO2 flux for budget checking
       f_co2flux(:,:) = 0.
+#   if defined key_omip_dic
+      !! AXY (06/08/18): applies to OMIP DIC too
+      f_pi_co2flux(:,:) = 0.
+#   endif      
 #  endif
 
@@ -562,4 +631,41 @@
                      f_ocndpco2_2d(ji,jj) = f_dpco2(ji,jj)
                   ENDIF
+#   if defined key_omip_dic
+                  !! AXY (06/08/18): diagnostics for OMIP PI DIC
+                  IF( med_diag%PI_ATM_PCO2%dgsave ) THEN
+                     f_pi_pco2a2d(ji,jj) = f_pi_pco2atm(ji,jj)
+                  ENDIF
+                  IF( med_diag%OCN_PCO2%dgsave ) THEN
+                     f_pi_pco2w2d(ji,jj) = f_pi_pco2w(ji,jj)
+                  ENDIF
+                  IF( med_diag%PI_CO2FLUX%dgsave ) THEN
+                     !! mmol/m3/d -> mmol/m2/d
+                     f_pi_co2flux2d(ji,jj) = f_pi_co2flux(ji,jj) * fse3t(ji,jj,1)
+                  ENDIF
+                  IF( med_diag%PI_TCO2%dgsave ) THEN
+                     f_PI_TDIC2d(ji,jj) = f_PI_TDIC(ji,jj)
+                  ENDIF
+                  IF( med_diag%PI_ATM_XCO2%dgsave ) THEN
+                     f_pi_xco2a_2d(ji,jj) = f_pi_xco2a(ji,jj)
+                  ENDIF
+                  IF( med_diag%PI_OCN_FCO2%dgsave ) THEN
+                     f_pi_fco2w_2d(ji,jj) = f_pi_fco2w(ji,jj)
+                  ENDIF
+                  IF( med_diag%PI_ATM_FCO2%dgsave ) THEN
+                     f_pi_fco2a_2d(ji,jj) = f_pi_fco2atm(ji,jj)
+                  ENDIF
+                  IF( med_diag%PI_CO2STARAIR%dgsave ) THEN
+                     f_pi_co2starair_2d(ji,jj) = f_pi_co2starair(ji,jj)
+                  ENDIF
+                  IF( med_diag%PI_OCN_DPCO2%dgsave ) THEN
+                     f_pi_ocndpco2_2d(ji,jj) = f_pi_dpco2(ji,jj)
+                  ENDIF                  
+                  IF( med_diag%PI_FGCO2%dgsave ) THEN
+                     !! convert from  mol/m2/day to kg/m2/s
+                     !! mmol-C/m3/d -> kg-CO2/m2/s
+                     pi_fgco2(ji,jj) = f_pi_co2flux(ji,jj) * fse3t(ji,jj,1) * &
+                          CO2flux_conv
+                  ENDIF
+#   endif                  
                ENDIF
             ENDIF
@@ -567,5 +673,5 @@
       ENDDO
 
-# endif
+#  endif
 
       !!-----------------------------------------------------------------
@@ -660,5 +766,6 @@
          ENDDO
       ENDDO
-
+# endif
+      
    END SUBROUTINE air_sea
 
Index: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag_slice.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag_slice.F90	(revision 10028)
+++ branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag_slice.F90	(revision 10045)
@@ -6,4 +6,5 @@
    !! History :
    !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!   -   ! 2018-08 (A. Yool)            add OMIP preindustrial DIC
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -236,4 +237,77 @@
             DEALLOCATE( f_ocndpco2_2d )
          ENDIF
+#  if defined key_omip_dic
+      !!----------------------------------------------------------------------
+      !! AXY (03/08/18): additional diagnostics for OMIP preindustrial DIC;
+      !!                 basically a subset of the normal diagnostics with
+      !!                 "PI_" stuck on the front; only includes those CO2
+      !!                 diagnostics that should differ under PI conditions
+      !!----------------------------------------------------------------------
+         IF( med_diag%PI_ATM_PCO2%dgsave ) THEN
+            CALL iom_put( "PI_ATM_PCO2"  , f_pi_pco2a2d )
+            DEALLOCATE( f_pi_pco2a2d )
+         ENDIF
+         IF( med_diag%PI_OCN_PH%dgsave ) THEN
+            CALL iom_put( "PI_OCN_PH"  , f_pi_ph2d )
+            DEALLOCATE( f_pi_ph2d )
+         ENDIF
+         IF( med_diag%PI_OCN_PCO2%dgsave ) THEN
+            CALL iom_put( "PI_OCN_PCO2"  , f_pi_pco2w2d )
+            DEALLOCATE( f_pi_pco2w2d )
+         ENDIF
+         IF( med_diag%PI_OCNH2CO3%dgsave ) THEN
+            CALL iom_put( "PI_OCNH2CO3"  , f_pi_h2co32d )
+            DEALLOCATE( f_pi_h2co32d )
+         ENDIF
+         IF( med_diag%PI_OCN_HCO3%dgsave ) THEN
+            CALL iom_put( "PI_OCN_HCO3"  , f_pi_hco32d )
+            DEALLOCATE( f_pi_hco32d )
+         ENDIF
+         IF( med_diag%PI_OCN_CO3%dgsave ) THEN
+            CALL iom_put( "PI_OCN_CO3"  , f_pi_co32d )
+            DEALLOCATE( f_pi_co32d )
+         ENDIF
+         IF( med_diag%PI_CO2FLUX%dgsave ) THEN
+            CALL iom_put( "PI_CO2FLUX"  , f_pi_co2flux2d )
+            DEALLOCATE( f_pi_co2flux2d )
+         ENDIF
+         IF( med_diag%PI_FGCO2%dgsave ) THEN
+            CALL iom_put( "PI_FGCO2"  , pi_fgco2 )
+            DEALLOCATE( pi_fgco2 )
+         ENDIF
+         IF( med_diag%PI_OM_CAL%dgsave ) THEN
+            CALL iom_put( "PI_OM_CAL"  , f_pi_omcal )
+            !! AXY: WHY IS THIS NOT DEALLOCATED?
+         ENDIF
+         IF( med_diag%PI_OM_ARG%dgsave ) THEN
+            CALL iom_put( "PI_OM_ARG"  , f_pi_omarg )
+            !! AXY: WHY IS THIS NOT DEALLOCATED?
+         ENDIF
+         IF( med_diag%PI_TCO2%dgsave ) THEN
+            CALL iom_put( "PI_TCO2"  , f_pi_TDIC2d )
+            DEALLOCATE( f_pi_TDIC2d )
+         ENDIF
+         IF( med_diag%PI_ATM_XCO2%dgsave ) THEN
+            CALL iom_put( "PI_ATM_XCO2"  ,   f_pi_xco2a_2d      )
+            DEALLOCATE( f_pi_xco2a_2d )
+         ENDIF
+         IF( med_diag%PI_OCN_FCO2%dgsave ) THEN
+            CALL iom_put( "PI_OCN_FCO2"  ,   f_pi_fco2w_2d      )
+            DEALLOCATE( f_pi_fco2w_2d )
+         ENDIF
+         IF( med_diag%PI_ATM_FCO2%dgsave ) THEN
+            CALL iom_put( "PI_ATM_FCO2"  ,   f_pi_fco2a_2d      )
+            DEALLOCATE( f_pi_fco2a_2d )
+         ENDIF
+         IF( med_diag%PI_CO2STARAIR%dgsave ) THEN
+            CALL iom_put( "PI_CO2STARAIR"  , f_pi_co2starair_2d )
+            DEALLOCATE( f_pi_co2starair_2d )
+         ENDIF
+         IF( med_diag%PI_OCN_DPCO2%dgsave ) THEN
+            CALL iom_put( "PI_OCN_DPCO2"  ,  f_pi_ocndpco2_2d   )
+            DEALLOCATE( f_pi_ocndpco2_2d )
+         ENDIF
+#  endif
+         
 # endif                     
       ELSE IF (jk.eq.i0100) THEN 
Index: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 10028)
+++ branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 10045)
@@ -7,4 +7,5 @@
    !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
    !!   -   ! 2017-08 (A. Yool)            Add slow-sinking detrius variables
+   !!   -   ! 2018-08 (A. Yool)            add OMIP preindustrial DIC
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -186,4 +187,17 @@
       f_omarg(:,:)    = 0.0
       f_omcal(:,:)    = 0.0
+# if defined key_omip_dic
+      zomd(:,:)          = 0.0
+      f_pi_xco2a(:,:)    = 0.0
+      f_pi_pco2w(:,:)    = 0.0
+      f_pi_ph(:,:)       = 0.0
+      f_pi_co2flux(:,:)  = 0.0 
+      f_pi_pco2atm(:,:)  = 0.0
+      f_pi_h2co3(:,:)    = 0.0
+      f_pi_hco3(:,:)     = 0.0
+      f_pi_co3(:,:)      = 0.0
+      f_pi_omarg(:,:)    = 0.0
+      f_pi_omcal(:,:)    = 0.0
+# endif      
       !!
       !! AXY (08/08/17): zero slow detritus fluxes
@@ -382,4 +396,10 @@
             f_pco2a2d(:,:)      = 0.0 !!
          ENDIF
+#  if defined key_omip_dic         
+         IF( med_diag%PI_ATM_PCO2%dgsave ) THEN
+            ALLOCATE( f_pi_pco2a2d(1:jpi, 1:jpj) )
+            f_pi_pco2a2d(:,:)      = 0.0 !!
+         ENDIF
+#  endif
 !!
 !! skip OCN_PH diagnostic here
@@ -389,4 +409,10 @@
             f_pco2w2d(:,:)      = 0.0 !!
          ENDIF
+#  if defined key_omip_dic         
+         IF( med_diag%PI_OCN_PCO2%dgsave ) THEN
+            ALLOCATE( f_pi_pco2w2d(1:jpi, 1:jpj) )
+            f_pi_pco2w2d(:,:)   = 0.0 !!
+         ENDIF
+#  endif
 !!
 !! skip OCNH2CO3, OCN_HCO3, OCN_CO3 diagnostics here
@@ -396,4 +422,10 @@
             f_co2flux2d(:,:)      = 0.0 !!
          ENDIF
+#  if defined key_omip_dic         
+         IF( med_diag%PI_CO2FLUX%dgsave ) THEN
+            ALLOCATE( f_pi_co2flux2d(1:jpi, 1:jpj) )
+            f_pi_co2flux2d(:,:)   = 0.0 !!
+         ENDIF
+#  endif
 !!
 !! skip OM_XXX diagnostics here
@@ -403,4 +435,10 @@
             f_TDIC2d(:,:)      = 0.0 !!
          ENDIF
+#  if defined key_omip_dic         
+         IF( med_diag%PI_TCO2%dgsave ) THEN
+            ALLOCATE( f_pi_TDIC2d(1:jpi, 1:jpj) )
+            f_pi_TDIC2d(:,:)      = 0.0 !!
+         ENDIF
+#  endif
          IF( med_diag%TALK%dgsave ) THEN
             ALLOCATE( f_TALK2d(1:jpi, 1:jpj) )
@@ -662,4 +700,30 @@
             f_ocndpco2_2d(:,:)      = 0.0 !!
          ENDIF
+#  if defined key_omip_dic         
+         IF( med_diag%PI_ATM_XCO2%dgsave ) THEN
+            ALLOCATE( f_pi_xco2a_2d(1:jpi, 1:jpj) )
+            f_pi_xco2a_2d(:,:)      = 0.0 !!
+         ENDIF
+         IF( med_diag%PI_OCN_FCO2%dgsave ) THEN
+            ALLOCATE( f_pi_fco2w_2d(1:jpi, 1:jpj) )
+            f_pi_fco2w_2d(:,:)      = 0.0 !!
+         ENDIF
+         IF( med_diag%PI_ATM_FCO2%dgsave ) THEN
+            ALLOCATE( f_pi_fco2a_2d(1:jpi, 1:jpj) )
+            f_pi_fco2a_2d(:,:)      = 0.0 !!
+         ENDIF
+         IF( med_diag%PI_CO2STARAIR%dgsave ) THEN
+            ALLOCATE( f_pi_co2starair_2d(1:jpi, 1:jpj) )
+            f_pi_co2starair_2d(:,:)      = 0.0 !!
+         ENDIF
+         IF( med_diag%OCN_PI_DPCO2%dgsave ) THEN
+            ALLOCATE( f_pi_ocndpco2_2d(1:jpi, 1:jpj) )
+            f_pi_ocndpco2_2d(:,:)      = 0.0 !!
+         ENDIF
+         IF( med_diag%PI_FGCO2%dgsave ) THEN
+            ALLOCATE( pi_fgco2(1:jpi, 1:jpj) )
+            pi_fgco2(:,:)      = 0.0 !!
+         ENDIF
+#  endif
 # endif  
          IF( med_diag%TPP3%dgsave ) THEN
Index: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90	(revision 10028)
+++ branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90	(revision 10045)
@@ -8,4 +8,5 @@
    !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
    !!   -   ! 2017-08 (A. Yool)            Slow detritus, ML-avg chl variables
+   !!   -   ! 2018-08 (A. Yool)            add OMIP preindustrial DIC
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -37,4 +38,7 @@
 # if defined key_mocsy
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zpho
+# endif
+# if defined key_omip_dic
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zomd
 # endif
 
@@ -162,4 +166,14 @@
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_pco2atm
 
+#  if defined key_omip_dic
+   !! AXY (06/08/18): OMIP PI DIC additions
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_pi_xco2a
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_pi_ph,f_pi_pco2w,f_pi_h2co3
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_pi_hco3,f_pi_co3,f_pi_co2flux
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_pi_TDIC
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_pi_omcal,f_pi_omarg
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_pi_pco2atm
+#  endif
+   
    !! Carbon, alkalinity production and consumption
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fcomm_resp
@@ -237,4 +251,10 @@
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_ocnk0_2d, f_co2starair_2d
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_ocndpco2_2d
+#  if defined key_omip_dic
+   !! AXY (06/08/18): OMIP PI DIC additions
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_pi_xco2a_2d, f_pi_fco2w_2d, f_pi_fco2a_2d
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_pi_co2starair_2d
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_pi_ocndpco2_2d
+#  endif
 # endif
    !!
@@ -250,4 +270,7 @@
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fbddtalk,fbddtdic,fbddtdife
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fbddtdin,fbddtdisi
+#  if defined key_omip_dic
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: pi_fgco2
+#  endif
    !!
    !! 3D var for new CMIP6 diagnostics
@@ -291,4 +314,7 @@
 # if defined key_mocsy
                zpho(jpi,jpj),                                         &
+# endif
+# if defined key_omip_dic
+               zomd(jpi,jpj),                                         &
 # endif
                fthetan(jpi,jpj),fprn(jpi,jpj),frn(jpi,jpj),           &
@@ -370,5 +396,13 @@
                f_kw660(jpi,jpj),f_o2flux(jpi,jpj),f_o2sat(jpi,jpj),   &
                f_omcal(jpi,jpj),f_omarg(jpi,jpj),                     &
-               f_pco2atm(jpi,jpj),                 &
+               f_pco2atm(jpi,jpj),                                    &
+#  if defined key_omip_dic
+               f_pi_xco2a(jpi,jpj),                                        &
+               f_pi_ph(jpi,jpj),f_pi_pco2w(jpi,jpj),f_pi_h2co3(jpi,jpj),   &
+               f_pi_hco3(jpi,jpj),f_pi_co3(jpi,jpj),f_pi_co2flux(jpi,jpj), &
+               f_pi_TDIC(jpi,jpj),                                         &
+               f_pi_omcal(jpi,jpj),f_pi_omarg(jpi,jpj),                    &
+               f_pi_pco2atm(jpi,jpj),                                      &
+#  endif               
                fcomm_resp(jpi,jpj),                                   &
                fcar_prod(jpi,jpj),fcar_cons(jpi,jpj),                 &
Index: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90	(revision 10028)
+++ branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90	(revision 10045)
@@ -8,4 +8,5 @@
    !!   -   ! 2017-08 (A. Yool)            Amend slow-detritus bug
    !!   -   ! 2017-08 (A. Yool)            Reformatting for clarity
+   !!   -   ! 2018-08 (A. Yool)            add OMIP preindustrial DIC
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -38,4 +39,7 @@
                                    f_benout_si,                              &
                                    f_co2flux, f_o2flux,                      &
+# if defined key_omip_dic
+                                   f_pi_co2flux,                             & ! AXY (06/08/18)
+# endif                                   
                                    f_riv_loc_alk, f_riv_loc_c,               &
                                    f_riv_loc_n, f_riv_loc_si,                & 
@@ -395,4 +399,10 @@
                fcar_prod(ji,jj) = fcar_prod(ji,jj) + ( fse3t(ji,jj,jk) *     &
                                                        fc_prod )
+#  if defined key_omip_dic
+               !! AXY (06/08/18): OMIP PI DIC has the same BGC fluxes as
+               !!                 normal DIC, with the exception of its
+               !!                 air-sea exchange; see below
+               btra(ji,jj,jpomd_lc) = btra(ji,jj,jpdic_lc)
+#  endif               
             ENDIF
          ENDDO
@@ -573,4 +583,10 @@
                   btra(ji,jj,jpdic_lc) = btra(ji,jj,jpdic_lc) + (b0 *              &
                                                            f_co2flux(ji,jj))
+#  if defined key_omip_dic
+                  !! AXY (06/08/18): air-sea CO2 flux is the only difference
+                  !!                 between DIC and OMIP PI DIC tracers
+                  btra(ji,jj,jpomd_lc) = btra(ji,jj,jpomd_lc) + (b0 *              &
+                                                           f_pi_co2flux(ji,jj))                  
+#  endif                   
                   !!
                   !! O2 flux (mol/m3/s -> mmol/m3/d)
@@ -624,4 +640,8 @@
                   IF (lwp) write (numout,*) 'btra(ji,jj,jpoxy_lc)(',jk,')  = ', &
                                             btra(ji,jj,jpoxy_lc)
+#   if defined key_omip_dic
+                  IF (lwp) write (numout,*) 'btra(ji,jj,jpomd_lc)(',jk,')  = ', &
+                                            btra(ji,jj,jpomd_lc)
+#   endif
 #  endif
                endif
@@ -731,4 +751,8 @@
                   tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) +                &
                                         (btra(ji,jj,jpoxy_lc) / 86400.)
+#  if defined key_omip_dic
+                  tra(ji,jj,jk,jpomd) = tra(ji,jj,jk,jpomd) +                &
+                                        (btra(ji,jj,jpomd_lc) / 86400.)
+#  endif
 # endif
                endif
Index: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90	(revision 10028)
+++ branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90	(revision 10045)
@@ -6,4 +6,5 @@
    !! History :
    !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
+   !!   -   ! 2018-08 (A. Yool)            no OMIP preindustrial DIC changed needed
    !!----------------------------------------------------------------------
 #if defined key_roam
Index: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus.F90	(revision 10028)
+++ branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus.F90	(revision 10045)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                         ***  MODULE detritus_mod  ***
-   !! Calculates detritus processes and fast-sinking detritus
+   !! Calculates slow detritus remineralisation and sinking
    !!======================================================================
    !! History :
Index: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus_fast_sink.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus_fast_sink.F90	(revision 10028)
+++ branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus_fast_sink.F90	(revision 10045)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                         ***  MODULE detritus_fast_sink_mod  ***
-   !! Calculates fast-sinking detritus
+   !! Calculates fast-sinking detritus processes (plus other diagnostics)
    !!======================================================================
    !! History :
Index: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/par_medusa.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/par_medusa.F90	(revision 10028)
+++ branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/par_medusa.F90	(revision 10045)
@@ -10,4 +10,5 @@
    !!              -   !  2011-04  (A. Yool) updated for ROAM project
    !!		   -   !  2013-03  (A. Yool) updated for v3.5 NEMO
+   !!              -   !  2018-08  (A. Yool) add OMIP preindustrial DIC
    !!----------------------------------------------------------------------
    !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
@@ -33,5 +34,9 @@
    LOGICAL, PUBLIC, PARAMETER ::   lk_medusa     = .TRUE.   !: PTS flag 
 # if defined key_roam
+#  if defined key_omip_dic
+   INTEGER, PUBLIC, PARAMETER ::   jp_medusa     =  16      !: number of PTS tracers
+#  else
    INTEGER, PUBLIC, PARAMETER ::   jp_medusa     =  15      !: number of PTS tracers
+#  endif
    INTEGER, PUBLIC, PARAMETER ::   jp_medusa_2d  =  225     !: additional 2d output arrays (used if ln_diatrc=T)
    INTEGER, PUBLIC, PARAMETER ::   jp_medusa_3d  =  5       !: additional 3d output arrays (used if ln_diatrc=T)
@@ -61,4 +66,7 @@
    INTEGER, PUBLIC, PARAMETER ::   jpalk  =  jp_lm + 14     !: alkalinity
    INTEGER, PUBLIC, PARAMETER ::   jpoxy  =  jp_lm + 15     !: dissolved oxygen concentration
+#  if defined key_omip_dic
+   INTEGER, PUBLIC, PARAMETER ::   jpomd  =  jp_lm + 16     !: OMIP preindustrial DIC
+#  end   
 # endif
 
@@ -80,4 +88,7 @@
    INTEGER, PUBLIC, PARAMETER ::   jpalk_lc  =  14     !: alkalinity
    INTEGER, PUBLIC, PARAMETER ::   jpoxy_lc  =  15     !: dissolved oxygen concentration
+#  if defined key_omip_dic
+   INTEGER, PUBLIC, PARAMETER ::   jpomd_lc  =  16     !: OMIP preindustrial DIC
+#  end   
 # endif
 
Index: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90	(revision 10028)
+++ branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90	(revision 10045)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                         ***  MODULE phytoplankton_mod  ***
-   !! Calculates the phytoplankton growth
+   !! Calculates phytoplankton growth
    !!======================================================================
    !! History :
Index: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90	(revision 10028)
+++ branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90	(revision 10045)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                         ***  MODULE plankton_mod  ***
-   !! Calculate the carbon chemistry for the whole ocean
+   !! Calculates phytoplankton and zooplankton terms
    !!======================================================================
    !! History :
Index: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 10028)
+++ branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 10045)
@@ -23,5 +23,7 @@
    !!  -   !  2016-11  (A. Yool)              Updated diags for CMIP6
    !!  -   !  2017-05  (A. Yool)              Added extra DMS calculation
+   !!  -   !  2017-05  (M. Stringer)          Code reorganisation
    !!  -   !  2017-11  (J. Palm, A. Yool)     Diagnose tracer excursions
+   !!  -   !  2018-08  (A. Yool)              add OMIP preindustrial DIC
    !!----------------------------------------------------------------------
    !!
@@ -72,4 +74,7 @@
                                             zalk, zdic, zoxy, zsal, ztmp,   &
 # endif
+# if defined key_omip_dic
+                                            zomd,                           &
+# endif                                            
 # if defined key_mocsy
                                             zpho,                           &
@@ -96,8 +101,14 @@
       USE oce,                        ONLY: tsb, tsn, PCO2a_in_cpl
       USE par_kind,                   ONLY: wp
-      USE par_medusa,                 ONLY: jpalk, jpchd, jpchn, jpdet,     &
-                                            jpdic, jpdin, jpdtc, jpfer,     &
-                                            jpoxy, jppds, jpphd, jpphn,     &
-                                            jpsil, jpzme, jpzmi
+      !! AXY (06/08/28): fragmented to permit tracer options
+      USE par_medusa,                 ONLY: jpchd, jpchn, jpphn, jpphd,     &
+                                            jppds, jpzmi, jpzme, jpdet,     &
+# if defined key_roam
+                                            jpdtc, jpdic, jpalk, jpoxy,     &
+#  if defined key_omip_dic
+                                            jpomd,                          &
+#  endif
+# endif
+                                            jpdin, jpsil, jpfer
       USE par_oce,                    ONLY: jp_sal, jp_tem, jpi, jpim1,     &
                                             jpj, jpjm1, jpk
@@ -329,5 +340,5 @@
              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 
+         f_xco2a(:,:) = 284.317          !! CMIP6 pre-industrial pCO2
       ELSE      
          !! xCO2 from file
@@ -371,5 +382,9 @@
          ENDIF
       ENDIF
- 
+
+#  if defined key_omip_dic
+         f_xco2a(:,:) = 284.317          !! CMIP6 pre-industrial pCO2
+#  endif
+      
       !! Writing xCO2 in output on start and on the 1st tsp of each  month
       IF ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        &
@@ -378,4 +393,7 @@
                                        ';  current date:', ndastp
           call trc_rst_dia_stat(f_xco2a(:,:), 'atm xCO2')
+#  if defined key_omip_dic
+          call trc_rst_dia_stat(f_pi_xco2a(:,:), 'PI atm xCO2')
+#  endif
       ENDIF
 # endif
@@ -456,4 +474,5 @@
       !! OPEN vertical loop
       DO jk = 1,jpk
+
          !! OPEN horizontal loops
          DO jj = 2,jpjm1
@@ -503,4 +522,5 @@
          ENDDO
 
+         !! OPEN horizontal loops
          DO jj = 2,jpjm1
             DO ji = 2,jpim1
@@ -524,4 +544,5 @@
 # if defined key_roam
          !! extra MEDUSA-2 tracers
+         !! OPEN horizontal loops
          DO jj = 2,jpjm1
             DO ji = 2,jpim1
@@ -559,4 +580,15 @@
 # endif
 
+# if defined key_roam && defined key_omip_dic
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  !! OMIP PI DIC
+                  zomd(ji,jj) = trn(ji,jj,jk,jpomd)
+               ENDIF
+            ENDDO
+         ENDDO
+# endif
+         
 # if defined key_roam
          !! ---------------------------------------------
@@ -676,4 +708,13 @@
                            ji, jj, jk, kt
                      endif
+#  if defined key_omip_dic
+                     if (trn(ji,jj,jk,jpomd).lt.0.) then
+                        IF (lwp) write (numout,*)                            &
+                           '------------------------------'
+                        IF (lwp) write (numout,*) 'NEGATIVE OMD ERROR =',    &
+                           trn(ji,jj,jk,jpomd)
+                        IF (lwp) write (numout,*) 'NEGATIVE OMD ERROR @',    &
+                           ji, jj, jk, kt
+                     endif
 #  endif
                   endif
Index: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcctl_medusa.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcctl_medusa.F90	(revision 10028)
+++ branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcctl_medusa.F90	(revision 10045)
@@ -9,6 +9,6 @@
    !!              -   !  2008-11  (A. Yool) continuing adaptation for MEDUSA
    !!              -   !  2010-03  (A. Yool) updated for branch inclusion
+   !!              -   !  2018-08  (A. Yool) add OMIP preindustrial DIC
    !!----------------------------------------------------------------------
-
 #if defined key_medusa
 
@@ -45,4 +45,16 @@
 ! -----------------------
 # if defined key_roam
+#  if defined key_omip_dic     
+      IF (jp_medusa /= 16) THEN 
+          IF (lwp) THEN 
+              WRITE (numout,*) ' ===>>>> : W A R N I N G '
+              WRITE (numout,*) ' =======   ============= '
+              WRITE (numout,*)                             &
+              &   ' STOP, change jp_medusa to 16 in '      &
+              &  ,' par_medusa.F90 '  
+          END IF 
+          STOP 'TRC_CTL'
+      END IF 
+#  else
       IF (jp_medusa /= 15) THEN 
           IF (lwp) THEN 
@@ -55,4 +67,5 @@
           STOP 'TRC_CTL'
       END IF 
+#  endif
 # else
       IF (jp_medusa /= 11) THEN 
Index: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcini_medusa.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcini_medusa.F90	(revision 10028)
+++ branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcini_medusa.F90	(revision 10045)
@@ -9,4 +9,5 @@
    !!              -   !  2010-03  (A. Yool) updated for branch inclusion
    !!              -   !  2011-04  (A. Yool) updated for ROAM project
+   !!              -   !  2018-08  (A. Yool) add OMIP preindustrial DIC
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -207,4 +208,9 @@
       !! dissolved oxygen                 (nominal average value; typically initialised from climatology)
       trn(:,:,:,jpoxy) = 175.
+#  if defined key_omip_dic
+      !!
+      !! OMIP preindustrial DIC           (nominal average value; typically initialised from climatology)
+      trn(:,:,:,jpomd) = 2330.
+#  endif
       !!
       IF(lwp) WRITE(numout,*) ' trc_ini_medusa: MEDUSA-2 fields initialised to defaults'
Index: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 10028)
+++ branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 10045)
@@ -12,4 +12,5 @@
    !!              -   !  2015-11  (J. Palmieri) added iom_use for diags
    !!              -   !  2016-11  (A. Yool) updated diags for CMIP6
+   !!              -   !  2018-08  (A. Yool) add OMIP preindustrial DIC
    !!----------------------------------------------------------------------
 #if defined key_medusa
@@ -490,4 +491,11 @@
          &   ' key_debug_medusa                                                       = INACTIVE'
 # endif
+# if defined key_omip_dic
+         WRITE(numout,*)     &
+         &   ' key_omip_dic                                                           = ACTIVE'
+#else
+         WRITE(numout,*)     &
+         &   ' key_omip_dic                                                           = INACTIVE'
+# endif
          WRITE(numout,*) ' '
 
@@ -2389,4 +2397,88 @@
       ENDIF
 
+# if defined key_omip_dic
+      !!----------------------------------------------------------------------
+      !! AXY (03/08/18): additional diagnostics for OMIP preindustrial DIC; 
+      !!                 basically a subset of the normal diagnostics with
+      !!                 "PI_" stuck on the front; only includes those CO2
+      !!                 diagnostics that should differ under PI conditions
+      !!----------------------------------------------------------------------
+      IF  (iom_use("PI_ATM_PCO2")) THEN 
+          med_diag%PI_ATM_PCO2%dgsave = .TRUE.
+      ELSE 
+          med_diag%PI_ATM_PCO2%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("PI_OCN_PH")) THEN 
+          med_diag%PI_OCN_PH%dgsave = .TRUE.
+      ELSE 
+          med_diag%PI_OCN_PH%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("PI_OCN_PCO2")) THEN 
+          med_diag%PI_OCN_PCO2%dgsave = .TRUE.
+      ELSE 
+          med_diag%PI_OCN_PCO2%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("PI_OCNH2CO3")) THEN 
+          med_diag%PI_OCNH2CO3%dgsave = .TRUE.
+      ELSE 
+          med_diag%PI_OCNH2CO3%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("PI_OCN_HCO3")) THEN 
+          med_diag%PI_OCN_HCO3%dgsave = .TRUE.
+      ELSE 
+          med_diag%PI_OCN_HCO3%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("PI_OCN_CO3")) THEN 
+          med_diag%PI_OCN_CO3%dgsave = .TRUE.
+      ELSE 
+          med_diag%PI_OCN_CO3%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("PI_CO2FLUX")) THEN 
+          med_diag%PI_CO2FLUX%dgsave = .TRUE.
+      ELSE 
+          med_diag%PI_CO2FLUX%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("PI_OM_CAL")) THEN 
+          med_diag%PI_OM_CAL%dgsave = .TRUE.
+      ELSE 
+          med_diag%PI_OM_CAL%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("PI_OM_ARG")) THEN 
+          med_diag%PI_OM_ARG%dgsave = .TRUE.
+      ELSE 
+          med_diag%PI_OM_ARG%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("PI_TCO2")) THEN 
+          med_diag%PI_TCO2%dgsave = .TRUE.
+      ELSE 
+          med_diag%PI_TCO2%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("PI_ATM_XCO2")) THEN 
+          med_diag%PI_ATM_XCO2%dgsave = .TRUE.
+      ELSE 
+          med_diag%PI_ATM_XCO2%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("PI_OCN_FCO2")) THEN 
+          med_diag%PI_OCN_FCO2%dgsave = .TRUE.
+      ELSE 
+          med_diag%PI_OCN_FCO2%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("PI_ATM_FCO2")) THEN 
+          med_diag%PI_ATM_FCO2%dgsave = .TRUE.
+      ELSE 
+          med_diag%PI_ATM_FCO2%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("PI_CO2STARAIR")) THEN 
+          med_diag%PI_CO2STARAIR%dgsave = .TRUE.
+      ELSE 
+          med_diag%PI_CO2STARAIR%dgsave = .FALSE.
+      ENDIF
+      IF  (iom_use("PI_OCN_DPCO2")) THEN 
+          med_diag%PI_OCN_DPCO2%dgsave = .TRUE.
+      ELSE 
+          med_diag%PI_OCN_DPCO2%dgsave = .FALSE.
+      ENDIF
+# endif      
+      
    END SUBROUTINE   trc_nam_iom_medusa
    
Index: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcopt_medusa.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcopt_medusa.F90	(revision 10028)
+++ branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcopt_medusa.F90	(revision 10045)
@@ -12,4 +12,6 @@
    !!              -   !  2008-11  (A. Yool) continuing adaptation for MEDUSA
    !!              -   !  2010-03  (A. Yool) updated for branch inclusion
+   !!              -   !  2014-08  (A. Yool) euphotic depth correction
+   !!              -   !  2015-07  (A. Yool) include rolling average irradiance
    !!----------------------------------------------------------------------
 #if defined key_medusa
Index: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsed_medusa.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsed_medusa.F90	(revision 10028)
+++ branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsed_medusa.F90	(revision 10045)
@@ -11,4 +11,9 @@
    !!              -   !  2010-03  (A. Yool) updated for branch inclusion
    !!              -   !  2011-04  (A. Yool) updated for ROAM project
+   !!              -   !  2014-XX  (A. Yool) various minor fixes
+   !!              -   !  2015-11  (J. Palmieri) IOM
+   !!              -   !  2016-03  (J. Palmieri) namelist dust changes
+   !!              -   !  2016-06  (J. Palmieri) atm coupling 
+   !!              -   !  2016-07  (J. Palmieri) dust unit conversions
    !!----------------------------------------------------------------------
 #if defined key_medusa
Index: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsms_medusa.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsms_medusa.F90	(revision 10028)
+++ branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsms_medusa.F90	(revision 10045)
@@ -8,4 +8,5 @@
    !!              -   !  2008-11  (A. Yool) continuing adaptation for MEDUSA
    !!              -   !  2010-03  (A. Yool) updated for branch inclusion
+   !!              -   !  2017-02  (A. Yool, J. Palmieri) move dust call
    !!              -   !  2017-08  (A. Yool) amend for slow detritus bug
    !!----------------------------------------------------------------------
Index: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/zooplankton.F90
===================================================================
--- branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/zooplankton.F90	(revision 10028)
+++ branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/zooplankton.F90	(revision 10045)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                         ***  MODULE zooplankton_mod  ***
-   !! Calculates the zooplankton grazing
+   !! Calculates zooplankton grazing
    !!======================================================================
    !! History :
