Index: /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/CONFIG/SHARED/field_def_bgc.xml
===================================================================
--- /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/CONFIG/SHARED/field_def_bgc.xml	(revision 11725)
+++ /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/CONFIG/SHARED/field_def_bgc.xml	(revision 11726)
@@ -376,6 +376,7 @@
        <field id= "SEAFLRC"    long_name="Seafloor flux of C"                        unit="mmol-C/m2/d"  />
        <field id= "SEAFLRCA"   long_name="Seafloor flux of CaCO3"                    unit="mmol-C/m2/d"  />
-       <field id= "MED_QSR"    long_name="Sea surface radiation"                     unit="W/m2"        />
-       <field id= "MED_XPAR"   long_name="Top level radiation"                       unit="W/m2"        />
+       <field id= "GEOENG_ALK" long_name="Surface Geo-ingeneering flux of CaCO3"     unit="mmol-eq/m2/d" />
+       <field id= "MAED_QSR"    long_name="Sea surface radiation"                    unit="W/m2"         />
+       <field id= "MED_XPAR"   long_name="Top level radiation"                       unit="W/m2"         />
        <field id= "INTFLX_N"   long_name="Integrated N  cycle fluxes"                unit="mmol-N/m2/d"  />
        <field id= "INTFLX_SI"  long_name="Integrated Si cycle fluxes"                unit="mmol-Si/m2/d" />
@@ -2268,4 +2269,5 @@
       <field field_ref= "SEAFLRC"    name="SEAFLRC"    />
       <field field_ref= "SEAFLRCA"   name="SEAFLRCA"   />
+      <field field_ref= "GEOENG_ALK" name="GEOENG_ALK" />
       <field field_ref= "MED_QSR"    name="MED_QSR"    />
       <field field_ref= "MED_XPAR"   name="MED_XPAR"   />
Index: /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_iomput.F90
===================================================================
--- /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_iomput.F90	(revision 11725)
+++ /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_iomput.F90	(revision 11726)
@@ -445,4 +445,10 @@
                                      fse3t(ji,jj,jk))
                ENDIF
+               !! AXY (16/10/19)
+               IF( med_diag%GEOENG_ALK%dgsave ) THEN
+                  geoengalk2d(ji,jj) = geoengalk2d(ji,jj) +                  &
+                                       (f_geoeng_alk(ji,jj) *                &
+                                       fse3t(ji,jj,jk))
+               ENDIF                 
                IF( med_diag%DETC%dgsave ) THEN
                   fslowc2d(ji,jj) = fslowc2d(ji,jj) +                        &
Index: /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90
===================================================================
--- /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 11725)
+++ /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 11726)
@@ -173,5 +173,5 @@
                   fq2 = fq0 + fq1
                   fq3 = f_benout_n(ji,jj)
-                  if (lwp) write (numout,'a,2i3,a,4f15,5)')                   &
+                  if (lwp) write (numout,'(a,2i3,a,4f15.5)')                   &
                      'AXY N   cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',      &
                      fq0,fq1,fq2,fq3
@@ -187,5 +187,5 @@
                   fq2 = fq0 + fq1
                   fq3 = f_benout_si(ji,jj)
-                  if (lwp) write (numout,'a,2i3,a,4f15,5)')                   &
+                  if (lwp) write (numout,'(a,2i3,a,4f15.5)')                   &
                      'AXY Si  cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',     &
                      fq0,fq1,fq2,fq3
@@ -202,5 +202,5 @@
                   fq3 = fq0 + fq1
                   fq4 = f_benout_c(ji,jj) + f_benout_ca(ji,jj)
-                  if (lwp) write (numout,'a,2i3,a,5f15,5)')                   &
+                  if (lwp) write (numout,'(a,2i3,a,5f15.5)')                   &
                      'AXY C   cons: (i,j)=',ji,jj,', (flx,ben,asf,err,out)=', &
                      fq0,fq1,fq2,fq3,fq4
@@ -216,5 +216,5 @@
                    fq2 = fq0 + fq1
                    fq3 = 2.0 * f_benout_ca(ji,jj)
-                   if (lwp) write (numout,'a,2i3,a,4f15,5)')                   &
+                   if (lwp) write (numout,'(a,2i3,a,4f15.5)')                   &
                       'AXY alk cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',     &
                       fq0,fq1,fq2,fq3
@@ -593,4 +593,9 @@
             DEALLOCATE( rivalk2d )
          ENDIF
+         !! AXY (16/10/19)
+         IF( med_diag%GEOENG_ALK%dgsave ) THEN
+            CALL iom_put( "GEOENG_ALK"  , geoengalk2d )
+            DEALLOCATE( geoengalk2d )
+         ENDIF
          IF( med_diag%DETC%dgsave ) THEN
             CALL iom_put( "DETC"  , fslowc2d )
Index: /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90
===================================================================
--- /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 11725)
+++ /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 11726)
@@ -506,4 +506,9 @@
             ALLOCATE( rivalk2d(1:jpi, 1:jpj) )
             rivalk2d(:,:)      = 0.0 !!
+         ENDIF
+         !! AXY (16/10/19)
+         IF( med_diag%GEOENG_ALK%dgsave ) THEN
+            ALLOCATE( geoengalk2d(1:jpi, 1:jpj) )
+            geoengalk2d(:,:)   = 0.0 !!
          ENDIF
          IF( med_diag%DETC%dgsave ) THEN
Index: /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90
===================================================================
--- /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90	(revision 11725)
+++ /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90	(revision 11726)
@@ -194,4 +194,8 @@
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_riv_loc_n,f_riv_loc_si
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_riv_loc_c, f_riv_loc_alk
+
+   !! AXY (16/10/2019)
+   !! geoengineering flux
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: f_geoeng_alk
 
    !! Jpalm -- 11-10-2015 -- adapt diag to iom_use
@@ -390,4 +394,5 @@
                f_riv_loc_n(jpi,jpj),f_riv_loc_si(jpi,jpj),            &
                f_riv_loc_c(jpi,jpj),f_riv_loc_alk(jpi,jpj),           &
+               f_geoeng_alk(jpi,jpj),                                 &
                STAT = bio_medusa_alloc)
 
Index: /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90
===================================================================
--- /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90	(revision 11725)
+++ /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90	(revision 11726)
@@ -71,5 +71,6 @@
 # endif
                                    zpds, zphd, zphn
-      USE dom_oce,           ONLY: e3t_0, gphit, mbathy, tmask
+      USE dom_oce,           ONLY: e3t_0, gphit, mbathy, tmask, nyear          ! AXY (16/10/19)
+
 # if defined key_vvl
       USE dom_oce,           ONLY: e3t_n
@@ -232,4 +233,6 @@
                !! 
                !! riverine flux
+               !! (NOTE: this is "protected" for near-surface addition of riverine
+               !! by array friver_dep which is zero below a defined model level)
                if ( jriver_n .gt. 0 ) then
                   f_riv_loc_n(ji,jj) = f_riv_n(ji,jj) *                      &
@@ -276,4 +279,6 @@
                !! 
                !! riverine flux
+               !! (NOTE: this is "protected" for near-surface addition of riverine
+               !! by array friver_dep which is zero below a defined model level)
                if ( jriver_si .gt. 0 ) then
                   f_riv_loc_si(ji,jj) = f_riv_si(ji,jj) *                    &
@@ -362,4 +367,6 @@
                !! 
                !! riverine flux
+               !! (NOTE: this is "protected" for near-surface addition of riverine
+               !! by array friver_dep which is zero below a defined model level)
                if ( jriver_c .gt. 0 ) then
                   f_riv_loc_c(ji,jj) = f_riv_c(ji,jj) *                      &
@@ -387,7 +394,9 @@
                !! 
                !! riverine flux
-               if ( jk .eq. 1 .and. jriver_c .gt. 0 ) then
-                  fc_prod = fc_prod + f_riv_c(ji,jj)
-               endif
+               !! (NOTE: the lines immediately below are incorrect and riverine
+               !! flux is already handled above)
+               !! if ( jk .eq. 1 .and. jriver_c .gt. 0 ) then
+               !!    fc_prod = fc_prod + f_riv_c(ji,jj)
+               !! endif
                !!
                btra(ji,jj,jpdic_lc) = b0 * ( fc_prod + fc_cons )
@@ -414,4 +423,6 @@
                !! 
                !! riverine flux
+               !! (NOTE: this is "protected" for near-surface addition of riverine
+               !! by array friver_dep which is zero below a defined model level)
                if ( jriver_alk .gt. 0 ) then
                   f_riv_loc_alk(ji,jj) = f_riv_alk(ji,jj) *                  &
@@ -419,4 +430,16 @@
                                          fse3t(ji,jj,jk)
                   fa_prod = fa_prod + f_riv_loc_alk(ji,jj)
+               endif
+	       !!
+               !! AXY (16/10/19)
+               !! geoengineering flux
+               !! (NOTE: this is only added to the surface ocean grid cell, only
+               !! within a particular time window, and only within a particular
+               !! sub-polar geographical range; these parameters could be passed up
+               !! to ROSE if they need different values for different experiments)
+	       if ( jk .eq. 1 .and. nyear .gt. 2014 .and. nyear .lt. 2050 .and. &
+                   gphit(ji,jj) .ge .-60.0 .and. gphit(ji,jj) .le. 70.0 ) then
+                  f_geoeng_alk(ji,jj) = 1.17 / fse3t(ji,jj,jk)
+                  fa_prod = fa_prod + f_geoeng_alk(ji,jj)
                endif
                !!  
Index: /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90
===================================================================
--- /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 11725)
+++ /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 11726)
@@ -1858,4 +1858,10 @@
           med_diag%RIV_ALK%dgsave = .FALSE.
       ENDIF
+      !! AXY (16/10/19)
+      IF  (iom_use("GEOENG_ALK")) THEN 
+          med_diag%GEOENG_ALK%dgsave = .TRUE.
+      ELSE 
+          med_diag%GEOENG_ALK%dgsave = .FALSE.
+      ENDIF
       IF  (iom_use("DETC")) THEN 
           med_diag%DETC%dgsave = .TRUE.
Index: /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/trc.F90
===================================================================
--- /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 11725)
+++ /branches/NERC/dev_r5518_GO6_GeoAlk_MEDUSA/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 11726)
@@ -131,5 +131,7 @@
                   RR_1000, IBEN_N, IBEN_FE, IBEN_C, IBEN_SI, IBEN_CA, OBEN_N, OBEN_FE, OBEN_C,       &
                   OBEN_SI, OBEN_CA, BEN_N, BEN_FE, BEN_C, BEN_SI, BEN_CA, RUNOFF, RIV_N, RIV_SI,     &
-                  RIV_C, RIV_ALK, DETC, SDC__100, SDC__200, SDC__500, SDC_1000, INVTC, INVTALK,      &
+                  RIV_C, RIV_ALK,                                                                    &
+		  GEOENG_ALK,                                                                        & ! AXY (16/10/19)
+		  DETC, SDC__100, SDC __200, SDC__500, SDC_1000, INVTC, INVTALK,                     &
                   INVTO2, LYSO_CA, COM_RESP, PN_LLOSS, PD_LLOSS, ZI_LLOSS, ZE_LLOSS, ZI_MES_N,       &
                   ZI_MES_D, ZI_MES_C, ZI_MESDC, ZI_EXCR, ZI_RESP, ZI_GROW, ZE_MES_N, ZE_MES_D,       &
