From 645c5537d48e617fd787439e327fca3dbfe3c35b Mon Sep 17 00:00:00 2001
From: Harikrishnan Ramesh <101573810+hra063@users.noreply.github.com>
Date: Fri, 14 Nov 2025 11:04:47 +0100
Subject: [PATCH] Added the Source modifications needed for the SST pacemaker
experiments
---
.../n1850_SST_pacemaker_Kuroshio/Readme.md | 29 +
.../SourceMods/src.blom/Readme.md | 34 +
.../SourceMods/src.blom/geoenv_file.F | 867 ++++++++++++++++++
.../SourceMods/src.blom/mod_grid.F90 | 177 ++++
.../SourceMods/src.drv/Readme.md | 60 ++
.../SourceMods/src.drv/ocn_comp_mct.F90 | 527 +++++++++++
.../SourceMods/src.drv/sumsbuff_mct.F | 177 ++++
7 files changed, 1871 insertions(+)
create mode 100644 cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/Readme.md
create mode 100644 cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.blom/Readme.md
create mode 100644 cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.blom/geoenv_file.F
create mode 100644 cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.blom/mod_grid.F90
create mode 100644 cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.drv/Readme.md
create mode 100644 cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.drv/ocn_comp_mct.F90
create mode 100644 cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.drv/sumsbuff_mct.F
diff --git a/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/Readme.md b/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/Readme.md
new file mode 100644
index 00000000..67e550d2
--- /dev/null
+++ b/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/Readme.md
@@ -0,0 +1,29 @@
+**Important usage Notes**
+
+External data format:
+
+The external SST file and pcmask must match the model domain and resolution.
+
+External SST file dimensions: (y, x, time)
+
+Time Handling:
+
+SST is applied at the current model timestep via the send buffer. Thus the time index (nstep_strm) should be aligned. (line 451 in ocn_comp_mct.F90)
+
+Branch runs:
+
+If starting from a specific model year (e.g., year 1320), the time index nstep_strm must account for prior simulation years. That is,
+
+nstep_strm = int(time) - int(time0) - 481434 ! (corresponding to year 1319 * 365 - 1)
+
+Hybrid runs:
+
+No special alignment is required; the time index can be set relative to the model start
+
+nstep_strm = int(time) - int(time0) - 1
+
+The current implementation hard-codes the time index for reading SST. We are trying to fix this with a dynamic approach that calculates the correct time index automatically, avoiding the need for hard-coded values.
+
+Mask Replacement:
+
+To run SST pacemaker experiments in a different region, provide a new mask file and update the pcmask variable.
diff --git a/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.blom/Readme.md b/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.blom/Readme.md
new file mode 100644
index 00000000..3dcdb7df
--- /dev/null
+++ b/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.blom/Readme.md
@@ -0,0 +1,34 @@
+Small modifications were made to the BLOM source code to enable reading and initialization of the pacemaker mask variable (pcmask) from an external NetCDF file.
+These updates allow BLOM to include region-specific SST forcing (overiding the BLOM SST) directly during grid initialization.
+
+Files Modified:
+
+1. mod_grid.F90
+
+ Added declaration of the new variable pcmask to the grid module.
+
+ Ensures the mask is initialized together with other grid variables in inivar_grid().
+
+2. geoenv_file.F
+
+ Updated to read the pacemaker_mask variable from the pacemaker_mask.nc file.
+
+ Integrated pcmask loading alongside other grid fields.
+
+ Ensures the mask is available globally after grid initialization.
+
+Description of pcmask:
+
+ Variable name: pcmask
+
+ Source file: pacemaker_mask.nc
+
+The mask has 1.0 in the core of the pacemaker region to 0.0 at the edges of the region (smooth tapering).
+
+Current configuration:
+
+The mask is set for the Kuroshio region in the North Pacific.
+
+Location of the file:
+
+ /cluster/work/users/hra063/So_t/pacemaker_mask.nc
diff --git a/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.blom/geoenv_file.F b/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.blom/geoenv_file.F
new file mode 100644
index 00000000..9512739a
--- /dev/null
+++ b/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.blom/geoenv_file.F
@@ -0,0 +1,867 @@
+! ------------------------------------------------------------------------------
+! Copyright (C) 2015-2020 Mats Bentsen, Ping-Gin Chiu
+!
+! This file is part of BLOM.
+!
+! BLOM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU Lesser General Public License as published by the Free
+! Software Foundation, either version 3 of the License, or (at your option)
+! any later version.
+!
+! BLOM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+! FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
+! more details.
+!
+! You should have received a copy of the GNU Lesser General Public License
+! along with BLOM. If not, see .
+! ------------------------------------------------------------------------------
+
+ subroutine geoenv_file
+c
+c --- ------------------------------------------------------------------
+c --- Get bathymetry and grid specification from file and compute
+c --- Coriolis parameter
+c --- ------------------------------------------------------------------
+c
+ use mod_config, only: inst_suffix
+ use mod_constants, only: rearth, pi, radian
+ use mod_xc
+ use mod_grid, only: grfile, qclon, qclat, pclon, pclat, uclon,
+ . uclat, vclon, vclat, scqx, scqy, scpx, scpy,
+ . scux, scuy, scvx, scvy, scq2, scp2, scu2,
+ . scv2, qlon, qlat, plon, plat, ulon, ulat,
+ . vlon, vlat, depths, corioq, coriop, betafp,
+ . angle, cosang, sinang, nwp, pcmask
+ use netcdf
+c
+ implicit none
+c
+ integer cwmlen
+ parameter (cwmlen=100)
+c
+ character (len = 256) :: nlfnm
+ character (len=80), dimension(cwmlen) :: cwmtag
+ character (len=1), dimension(cwmlen) :: cwmedg
+ real, dimension(itdm,jtdm) :: tmpg
+ real, dimension(cwmlen) :: cwmwth
+ integer, dimension(cwmlen) :: cwmi,cwmj
+ integer, dimension(3) :: start,count
+ integer i,j,k,status,ncid,dimid,varid,ios,ncwm,l
+ logical fexist
+c
+ namelist /cwmod/ cwmtag,cwmedg,cwmi,cwmj,cwmwth
+c
+c --- ------------------------------------------------------------------
+c --- read grid information from grid file
+c --- ------------------------------------------------------------------
+c
+ if (mnproc.eq.1) then
+ write (lp,'(2a)') ' reading grid information from ',trim(grfile)
+ call flush(lp)
+c
+c --- - open netcdf file
+ status=nf90_open(grfile,nf90_nowrite,ncid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(4a)') ' nf90_open: ',trim(grfile),': ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+c
+c --- - check dimensions
+ status=nf90_inq_dimid(ncid,'x',dimid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_dimid: x: ',nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_inquire_dimension(ncid,dimid,len=i)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inquire_dimension: x: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_inq_dimid(ncid,'y',dimid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_dimid: y: ',nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_inquire_dimension(ncid,dimid,len=j)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inquire_dimension: y: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ if (i.ne.itdm.or.j.ne.jtdm) then
+ write (lp,'(2a)') ' wrong dimensions in ',trim(grfile)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+c
+c --- - read bathymetry
+ status=nf90_inq_varid(ncid,'pdepth',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: pdepth: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: pdepth: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+c
+c --- - count number of wet points for subsequent xcsum testing
+ nwp=0
+ do j=1,jtdm
+ do i=1,itdm
+ if (tmpg(i,j).gt.0.) nwp=nwp+1
+ enddo
+ enddo
+ endif
+c
+ call xcaput(tmpg,depths,1)
+c
+c --- read grid coordinates
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'qlon',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: qlon: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: qlon: ',nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,qlon,1)
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'qlat',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: qlat: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: qlat: ',nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,qlat,1)
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'plon',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: plon: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: plon: ',nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,plon,1)
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'plat',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: plat: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: plat: ',nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,plat,1)
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'ulon',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: ulon: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: ulon: ',nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,ulon,1)
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'ulat',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: ulat: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: ulat: ',nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,ulat,1)
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'vlon',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: vlon: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: vlon: ',nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,vlon,1)
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'vlat',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: vlat: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: vlat: ',nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,vlat,1)
+c
+ start(1)=1
+ start(2)=1
+ count(1)=itdm
+ count(2)=jtdm
+ count(3)=1
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'qclon',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: qclon: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ do k=1,4
+ if (mnproc.eq.1) then
+ start(3)=k
+ status=nf90_get_var(ncid,varid,tmpg,start,count)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: qclon: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,qclon(1-nbdy,1-nbdy,k),1)
+ enddo
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'qclat',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: qclat: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ do k=1,4
+ if (mnproc.eq.1) then
+ start(3)=k
+ status=nf90_get_var(ncid,varid,tmpg,start,count)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: qclat: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,qclat(1-nbdy,1-nbdy,k),1)
+ enddo
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'pclon',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: pclon: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ do k=1,4
+ if (mnproc.eq.1) then
+ start(3)=k
+ status=nf90_get_var(ncid,varid,tmpg,start,count)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: pclon: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,pclon(1-nbdy,1-nbdy,k),1)
+ enddo
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'pclat',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: pclat: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ do k=1,4
+ if (mnproc.eq.1) then
+ start(3)=k
+ status=nf90_get_var(ncid,varid,tmpg,start,count)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: pclat: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,pclat(1-nbdy,1-nbdy,k),1)
+ enddo
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'uclon',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: uclon: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ do k=1,4
+ if (mnproc.eq.1) then
+ start(3)=k
+ status=nf90_get_var(ncid,varid,tmpg,start,count)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: uclon: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,uclon(1-nbdy,1-nbdy,k),1)
+ enddo
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'uclat',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: uclat: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ do k=1,4
+ if (mnproc.eq.1) then
+ start(3)=k
+ status=nf90_get_var(ncid,varid,tmpg,start,count)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: uclat: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,uclat(1-nbdy,1-nbdy,k),1)
+ enddo
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'vclon',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: vclon: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ do k=1,4
+ if (mnproc.eq.1) then
+ start(3)=k
+ status=nf90_get_var(ncid,varid,tmpg,start,count)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: vclon: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,vclon(1-nbdy,1-nbdy,k),1)
+ enddo
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'vclat',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: vclat: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ do k=1,4
+ if (mnproc.eq.1) then
+ start(3)=k
+ status=nf90_get_var(ncid,varid,tmpg,start,count)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: vclat: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,vclat(1-nbdy,1-nbdy,k),1)
+ enddo
+c
+c --- read scale factors
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'qdx',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: qdx: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: qdx: ',nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,scqx,1)
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'qdy',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: qdy: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: qdy: ',nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,scqy,1)
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'pdx',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: pdx: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: pdx: ',nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,scpx,1)
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'pdy',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: pdy: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: pdy: ',nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,scpy,1)
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'udx',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: udx: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: udx: ',nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,scux,1)
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'udy',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: udy: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: udy: ',nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,scuy,1)
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'vdx',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: vdx: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: vdx: ',nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,scvx,1)
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'vdy',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: vdy: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: vdy: ',nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,scvy,1)
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'qarea',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: qarea: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: qarea: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,scq2,1)
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'parea',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: parea: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: parea: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,scp2,1)
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'uarea',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: uarea: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: uarea: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,scu2,1)
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'varea',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: varea: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: varea: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,scv2,1)
+c
+ if (mnproc.eq.1) then
+ status=nf90_inq_varid(ncid,'angle',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: angle: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ if (mnproc.eq.1) then
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: angle: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,angle,1)
+c
+c --- close grid information file
+c
+ if (mnproc.eq.1) then
+ status=nf90_close(ncid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(4a)') ' nf90_close: ',trim(grfile),': ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+c
+ if (mnproc.eq.1) then
+ write (lp,'(2a)') ' reading pacemaker mask from ',
+ . trim('pacemaker_mask.nc')
+ call flush(lp)
+ status=nf90_open(
+ . '/cluster/work/users/hra063/So_t/pacemaker_mask.nc',
+ . nf90_nowrite,ncid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(4a)') ' nf90_open: pc_mask: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_inq_varid(ncid,'pc_mask',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_inq_varid: mask: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ status=nf90_get_var(ncid,varid,tmpg)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(2a)') ' nf90_get_var: mask: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcaput(tmpg,pcmask,1)
+c
+c --- close pacemaker file
+c
+ if (mnproc.eq.1) then
+ status=nf90_close(ncid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(4a)') ' nf90_close: pc_mask: ',
+ . nf90_strerror(status)
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+
+c --- ------------------------------------------------------------------
+c --- Apply channel width modifications if specified in namelist
+c --- ------------------------------------------------------------------
+c
+ if (mnproc.eq.1) then
+ nlfnm='ocn_in'//trim(inst_suffix)
+ inquire(file=nlfnm,exist=fexist)
+ if (fexist) then
+ open (unit=nfu,file=nlfnm,status='old',action='read')
+ else
+ nlfnm='limits'//trim(inst_suffix)
+ inquire(file=nlfnm,exist=fexist)
+ if (fexist) then
+ open (unit=nfu,file=nlfnm,status='old',action='read')
+ else
+ write (lp,*) 'geoenv_file: could not find namelist file!'
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ cwmtag=''
+ read (unit=nfu,nml=cwmod,iostat=ios)
+ close (unit=nfu)
+ endif
+ call xcbcst(ios)
+ if (ios.ne.0) then
+ if (mnproc.eq.1) then
+ write (lp,*) 'geoenv_file: no valid channel width '//
+ . 'modifications found in namelist.'
+ endif
+ else
+ if (mnproc.eq.1) then
+ ncwm=1
+ do while (cwmtag(ncwm).ne.'')
+ ncwm=ncwm+1
+ if (ncwm.gt.cwmlen) exit
+ enddo
+ ncwm=ncwm-1
+ endif
+ call xcbcst(ncwm)
+ do l=1,ncwm
+ if (mnproc.eq.1) then
+ write (lp,'(a,i3,a)') ' channel width modification number ',
+ . l,':'
+ write (lp,'(3a)') ' cwmtag = ''',trim(cwmtag(l)),''''
+ write (lp,'(3a)') ' cwmedg = ''',cwmedg(l),''''
+ write (lp,'(a,i5)') ' cwmi =',cwmi(l)
+ write (lp,'(a,i5)') ' cwmj =',cwmj(l)
+ write (lp,'(a,e10.3)') ' cwmwth = ',cwmwth(l)
+ if (cwmedg(l).ne.'u'.and.cwmedg(l).ne.'v') then
+ write (lp,*)
+ .'geoenv_file: channel width modification grid cell edge '//
+ .'specification must be either ''u'' or ''v''!'
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ if (cwmi(l).lt.1.or.cwmi(l).gt.itdm.or.
+ . cwmj(l).lt.1.or.cwmj(l).gt.jtdm) then
+ write (lp,*)
+ .'geoenv_file: channel width modification indices out of bounds!'
+ call xchalt('(geoenv_file)')
+ stop '(geoenv_file)'
+ endif
+ endif
+ call xcbcst(cwmtag(l))
+ call xcbcst(cwmedg(l))
+ call xcbcst(cwmi(l))
+ call xcbcst(cwmj(l))
+ call xcbcst(cwmwth(l))
+ if (cwmi(l).gt.i0.and.cwmi(l).le.i0+ii.and.
+ . cwmj(l).gt.j0.and.cwmj(l).le.j0+jj) then
+ i=cwmi(l)-i0
+ j=cwmj(l)-j0
+ if (cwmedg(l).eq.'u') then
+ write (lp,*) 'blom: geoenv_file: ',trim(cwmtag(l)),
+ . i+i0,j+j0,'scuy:',scuy(i,j),'->',cwmwth(l)
+ scu2(i,j) = scu2(i,j)*cwmwth(l)/scuy(i,j)
+ scuy(i,j) = cwmwth(l)
+ else
+ write (lp,*) 'blom: geoenv_file: ',trim(cwmtag(l)),
+ . i+i0,j+j0,'scvx:',scvx(i,j),'->',cwmwth(l)
+ scv2(i,j) = scv2(i,j)*cwmwth(l)/scvx(i,j)
+ scvx(i,j) = cwmwth(l)
+ endif
+ endif
+ enddo
+ endif
+c
+c --- ------------------------------------------------------------------
+c --- Get correct units of scale factors, precompute cosine and sine of
+c --- local angle of i-direction and with eastward direction, and
+c --- compute Coriolis and beta plane parameter
+c --- ------------------------------------------------------------------
+c
+c$OMP PARALLEL DO PRIVATE(i)
+ do j=1,jj
+ do i=1,ii
+c
+ scqx(i,j)=scqx(i,j)*1.e2
+ scqy(i,j)=scqy(i,j)*1.e2
+ scpx(i,j)=scpx(i,j)*1.e2
+ scpy(i,j)=scpy(i,j)*1.e2
+ scux(i,j)=scux(i,j)*1.e2
+ scuy(i,j)=scuy(i,j)*1.e2
+ scvx(i,j)=scvx(i,j)*1.e2
+ scvy(i,j)=scvy(i,j)*1.e2
+ scq2(i,j)=scq2(i,j)*1.e4
+ scp2(i,j)=scp2(i,j)*1.e4
+ scu2(i,j)=scu2(i,j)*1.e4
+ scv2(i,j)=scv2(i,j)*1.e4
+c
+ cosang(i,j)=cos(angle(i,j))
+ sinang(i,j)=sin(angle(i,j))
+c
+ corioq(i,j)=sin(qlat(i,j)/radian)*4.*pi/86164.
+ coriop(i,j)=sin(plat(i,j)/radian)*4.*pi/86164.
+ betafp(i,j)=cos(plat(i,j)/radian)*4.*pi/(86164.*rearth)
+c
+ enddo
+ enddo
+c$OMP END PARALLEL DO
+c
+ return
+ end
diff --git a/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.blom/mod_grid.F90 b/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.blom/mod_grid.F90
new file mode 100644
index 00000000..5ce4d7e7
--- /dev/null
+++ b/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.blom/mod_grid.F90
@@ -0,0 +1,177 @@
+! ------------------------------------------------------------------------------
+! Copyright (C) 2020 Mats Bentsen
+!
+! This file is part of BLOM.
+!
+! BLOM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU Lesser General Public License as published by the Free
+! Software Foundation, either version 3 of the License, or (at your option)
+! any later version.
+!
+! BLOM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+! FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
+! more details.
+!
+! You should have received a copy of the GNU Lesser General Public License
+! along with BLOM. If not, see .
+! ------------------------------------------------------------------------------
+
+module mod_grid
+! ------------------------------------------------------------------------------
+! This module contains declaration of grid variables.
+! ------------------------------------------------------------------------------
+
+ use mod_types, only: r8
+ use mod_constants, only: spval
+ use mod_xc, only: idm, jdm, kdm, nbdy, ii, jj, kk
+
+ implicit none
+
+ private
+
+ ! Variable to be set in namelist:
+ character(len = 256) :: &
+ grfile ! Name of file containing grid specification.
+
+ real(r8), dimension(1 - nbdy:idm + nbdy, 1 - nbdy:jdm + nbdy, kdm) :: &
+ sigmar ! Reference potential density [g cm-3].
+
+ real(r8), dimension(1 - nbdy:idm + nbdy, 1 - nbdy:jdm + nbdy, 4) :: &
+ qclon, & ! Longitude of q-cell corners [degrees].
+ qclat, & ! Latitude of q-cell corners [degrees].
+ pclon, & ! Longitude of p-cell corners [degrees].
+ pclat, & ! Latitude of p-cell corners [degrees].
+ uclon, & ! Longitude of u-cell corners [degrees].
+ uclat, & ! Latitude of u-cell corners [degrees].
+ vclon, & ! Longitude of v-cell corners [degrees].
+ vclat ! Latitude of v-cell corners [degrees].
+
+ real(r8), dimension(1 - nbdy:idm + nbdy, 1 - nbdy:jdm + nbdy) :: &
+ scqx, & ! Grid size in x-direction centered at q-point [cm].
+ scqy, & ! Grid size in y-direction centered at q-point [cm].
+ scpx, & ! Grid size in x-direction centered at p-point [cm].
+ scpy, & ! Grid size in y-direction centered at p-point [cm].
+ scux, & ! Grid size in x-direction centered at u-point [cm].
+ scuy, & ! Grid size in y-direction centered at u-point [cm].
+ scvx, & ! Grid size in x-direction centered at v-point [cm].
+ scvy, & ! Grid size in y-direction centered at v-point [cm].
+ scq2, & ! Area of grid cell centered at q-point [cm2].
+ scp2, & ! Area of grid cell centered at p-point [cm2].
+ scu2, & ! Area of grid cell centered at u-point [cm2].
+ scv2, & ! Area of grid cell centered at v-point [cm2].
+ scq2i, & ! Multiplicative inverse of scq2 [cm-2].
+ scp2i, & ! Multiplicative inverse of scp2 [cm-2].
+ scuxi, & ! Multiplicative inverse of scux [cm-1].
+ scuyi, & ! Multiplicative inverse of scuy [cm-1].
+ scvxi, & ! Multiplicative inverse of scvx [cm-1].
+ scvyi, & ! Multiplicative inverse of scvy [cm-1].
+ qlon, & ! Longitude of q-point [degrees].
+ qlat, & ! Latitude of q-point [degrees].
+ plon, & ! Longitude of p-point [degrees].
+ plat, & ! Latitude of p-point [degrees].
+ ulon, & ! Longitude of u-point [degrees].
+ ulat, & ! Latitude of u-point [degrees].
+ vlon, & ! Longitude of v-point [degrees].
+ vlat, & ! Latitude of v-point [degrees].
+ depths, & ! Water depth [m].
+ corioq, & ! Coriolis parameter at q-point [s-1].
+ coriop, & ! Coriolis parameter at p-point [s-1].
+ betafp, & ! Derivative of Coriolis parameter with respect to meridional
+ ! distance at p-point [cm-1 s-1].
+ angle, & ! Local angle between x-direction and eastward direction at
+ ! p-points [radians].
+ cosang, & ! Cosine of local angle between x-direction and eastward
+ ! direction at p-points [].
+ sinang, & ! Sine of local angle between x-direction and eastward
+ ! direction at p-points [].
+ pcmask, & ! pacemaker mask
+ So_t ! pacemaker SST
+
+ real(r8) :: &
+ area ! Total grid area [cm2].
+
+ integer :: &
+ nwp ! Number of wet grid cells.
+
+ public :: grfile, sigmar, &
+ qclon, qclat, pclon, pclat, uclon, uclat, vclon, vclat, &
+ scqx, scqy, scpx, scpy, scux, scuy, scvx, scvy, &
+ scq2, scp2, scu2, scv2, scq2i, scp2i, &
+ scuxi, scuyi, scvxi, scvyi, &
+ qlon, qlat, plon, plat, ulon, ulat, vlon, vlat, &
+ depths, corioq, coriop, betafp, angle, cosang, sinang, &
+ area, nwp, pcmask, So_t, &
+ inivar_grid
+
+contains
+
+ subroutine inivar_grid
+ ! ---------------------------------------------------------------------------
+ ! Initialize arrays.
+ ! ---------------------------------------------------------------------------
+
+ integer :: i, j, k
+
+ !$omp parallel do private(i, k)
+ do j = 1 - nbdy, jj + nbdy
+ do k = 1, kk
+ do i = 1 - nbdy, ii + nbdy
+ sigmar(i, j, k) = spval
+ enddo
+ enddo
+ do k = 1, 4
+ do i = 1 - nbdy, ii + nbdy
+ qclon(i, j, k) = spval
+ qclat(i, j, k) = spval
+ pclon(i, j, k) = spval
+ pclat(i, j, k) = spval
+ uclon(i, j, k) = spval
+ uclat(i, j, k) = spval
+ vclon(i, j, k) = spval
+ vclat(i, j, k) = spval
+ enddo
+ enddo
+ do i = 1 - nbdy, ii + nbdy
+ scqx(i, j) = spval
+ scqy(i, j) = spval
+ scpx(i, j) = spval
+ scpy(i, j) = spval
+ scux(i, j) = spval
+ scuy(i, j) = spval
+ scvx(i, j) = spval
+ scvy(i, j) = spval
+ scq2(i, j) = spval
+ scp2(i, j) = spval
+ scu2(i, j) = spval
+ scv2(i, j) = spval
+ scq2i(i, j) = spval
+ scp2i(i, j) = spval
+ scuxi(i, j) = spval
+ scuyi(i, j) = spval
+ scvxi(i, j) = spval
+ scvyi(i, j) = spval
+ qlon(i, j) = spval
+ qlat(i, j) = spval
+ plon(i, j) = spval
+ plat(i, j) = spval
+ ulon(i, j) = spval
+ ulat(i, j) = spval
+ vlon(i, j) = spval
+ vlat(i, j) = spval
+ depths(i, j) = spval
+ corioq(i, j) = spval
+ coriop(i, j) = spval
+ betafp(i, j) = spval
+ angle(i, j) = spval
+ cosang(i, j) = spval
+ sinang(i, j) = spval
+ pcmask(i, j) = spval
+ So_t(i, j) = spval
+ enddo
+ enddo
+ !$omp end parallel do
+
+ end subroutine inivar_grid
+
+end module mod_grid
diff --git a/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.drv/Readme.md b/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.drv/Readme.md
new file mode 100644
index 00000000..243bd2e7
--- /dev/null
+++ b/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.drv/Readme.md
@@ -0,0 +1,60 @@
+Modifications were made to the DRV source code to override the BLOM Sea Surface Temperature (SST) over a defined pacemaker region with SSTs from a separate NetCDF file. These changes do not affect other fields or fluxes outside the pacemaker region.
+
+
+Key Changes
+
+1. External SST (So_t) and Pacemaker Mask (pcmask)
+
+SST is read from an external NetCDF file using the subroutine "get_So_t" (lines 425 - 524 in ocn_comp_mct.F90)
+
+The pacemaker mask (pcmask) defines where the SST override is applied. It has values range from 1 in the core of the pacemaker region to 0 at the edges for smooth spatial transition. It is currently set for the Kuroshio region; other regions can be used by updating the mask.
+
+Initialization of these variables are done in mod_grid.f90 /geoenv.F files in src.blom.
+
+2. Send Buffer Accumulation (sumsbuff_mct.F)
+
+Modified to apply external SST wherever pcmask > 0 (see sumsbuff_mct subroutine, lines 112–115).
+
+Inside the masked region, model SST is fully overridden by the external SST:
+
+ sbuff(i,j,index_o2x_So_t) = sbuff(i,j,index_o2x_So_t) + temp(i,j,k1n)*baclin + pcmask(i,j)*(So_t(i,j)-temp(i,j,k1n))*baclin
+
+Outside the masked region, model SST remains unchanged.
+
+Other prognostic fields are unaffected.
+
+3. Coupling Integration
+
+The SST override is applied in the MCT send buffer, ensuring proper integration with the BLOM coupling cycle.
+
+
+Important usage Notes
+
+External data format:
+
+The external SST file and pcmask must match the model domain and resolution.
+
+External SST file dimensions: (y, x, time)
+
+Time Handling:
+
+SST is applied at the current model timestep via the send buffer. Thus the time index (nstep_strm) should be aligned. (line 451 in ocn_comp_mct.F90)
+
+Branch runs:
+
+If starting from a specific model year (e.g., year 1320), the time index nstep_strm must account for prior simulation years. That is,
+
+nstep_strm = int(time) - int(time0) - 481434 ! (corresponding to year 1319 * 365 - 1)
+
+
+Hybrid runs:
+
+No special alignment is required; the time index can be set relative to the model start
+
+nstep_strm = int(time) - int(time0) - 1
+
+The current implementation hard-codes the time index for reading SST. We are trying to fix this with a dynamic approach that calculates the correct time index automatically, avoiding the need for hard-coded values.
+
+Mask Replacement:
+
+To run SST pacemaker experiments in a different region, provide a new mask file and update the pcmask variable.
diff --git a/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.drv/ocn_comp_mct.F90 b/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.drv/ocn_comp_mct.F90
new file mode 100644
index 00000000..8bbf267b
--- /dev/null
+++ b/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.drv/ocn_comp_mct.F90
@@ -0,0 +1,527 @@
+! ------------------------------------------------------------------------------
+! Copyright (C) 2008-2020 Mats Bentsen, Alok Kumar Gupta, Ping-Gin Chiu
+!
+! This file is part of BLOM.
+!
+! BLOM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU Lesser General Public License as published by the Free
+! Software Foundation, either version 3 of the License, or (at your option)
+! any later version.
+!
+! BLOM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+! FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
+! more details.
+!
+! You should have received a copy of the GNU Lesser General Public License
+! along with BLOM. If not, see .
+! ------------------------------------------------------------------------------
+
+module ocn_comp_mct
+
+ ! -------------------------------------------------------------------
+ ! BLOM interface module for the cesm cpl7 mct system
+ ! -------------------------------------------------------------------
+
+ ! CESM modules
+ use mct_mod
+ use esmf, only: ESMF_Clock
+ use seq_cdata_mod, only: seq_cdata, seq_cdata_setptrs
+ use seq_infodata_mod, only: &
+ seq_infodata_type, seq_infodata_getdata, &
+ seq_infodata_putdata, seq_infodata_start_type_cont, &
+ seq_infodata_start_type_brnch, seq_infodata_start_type_start
+ use seq_flds_mod
+ use seq_timemgr_mod, only: &
+ seq_timemgr_EClockGetData, seq_timemgr_RestartAlarmIsOn, &
+ seq_timemgr_EClockDateInSync
+ use seq_comm_mct, only: seq_comm_suffix, seq_comm_inst, seq_comm_name
+ use shr_file_mod, only: &
+ shr_file_getUnit, shr_file_setIO, &
+ shr_file_getLogUnit, shr_file_getLogLevel, &
+ shr_file_setLogUnit, shr_file_setLogLevel, &
+ shr_file_freeUnit
+ use shr_cal_mod, only: shr_cal_date2ymd
+ use shr_sys_mod, only: shr_sys_abort, shr_sys_flush
+ use perf_mod, only: t_startf, t_stopf
+
+ use mod_types, only: r8
+ use mod_config, only: inst_index, inst_name, inst_suffix
+ use mod_time, only: blom_time, time, time0, &
+ date0, date, nday_of_year
+ use mod_cesm, only: runid_cesm, runtyp_cesm, ocn_cpl_dt_cesm
+ use mod_xc
+ use mod_grid, only: So_t
+ use blom_cpl_indices
+ use netcdf
+
+ implicit none
+
+ public :: ocn_init_mct, ocn_run_mct, ocn_final_mct
+ private :: ocn_SetGSMap_mct
+
+ private
+
+ integer, dimension(:), allocatable :: &
+ perm ! Permutation array to reorder points
+
+ real(r8), dimension(:,:,:), allocatable :: &
+ sbuff ! Accumulated sum of send buffer quantities for
+ ! averaging before being sent
+! real(r8), dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: So_t
+ integer :: &
+ lsize, & ! Size of attribute vector
+ jjcpl, & ! y-dimension of local ocean domain to send/receive
+ ! fields to coupler
+ nsend, & ! Number of fields to be sent to coupler
+ nrecv ! Number of fields to be received from coupler
+
+ real(r8) :: &
+ tlast_coupled, & ! Time since last coupling
+ precip_fact ! Correction factor for precipitation and runoff
+
+ logical :: &
+ lsend_precip_fact ! Flag for sending precipitation/runoff factor
+
+ contains
+
+
+ subroutine ocn_init_mct(EClock, cdata_o, x2o_o, o2x_o, NLFilename)
+
+ ! Input/output arguments
+
+ type(ESMF_Clock) , intent(inout) :: EClock
+ type(seq_cdata) , intent(inout) :: cdata_o
+ type(mct_aVect) , intent(inout) :: x2o_o, o2x_o
+ character (len=*), optional , intent(in) :: NLFilename ! Namelist filename
+ ! Local variables
+
+ type(mct_gsMap), pointer :: gsMap_ocn
+ type(mct_gGrid), pointer :: dom_ocn
+ type(seq_infodata_type), pointer :: infodata ! Input init object
+ integer :: OCNID, mpicom_ocn, shrlogunit, shrloglev, &
+ start_ymd, start_tod, start_year, start_day, start_month
+ logical :: exists
+ character(len=32) :: starttype
+
+ ! Set cdata pointers
+ call seq_cdata_setptrs(cdata_o, ID = OCNID, mpicom = mpicom_ocn, &
+ gsMap = gsMap_ocn, dom = dom_ocn, &
+ infodata = infodata)
+
+ ! Set communicator to be used by blom
+ mpicom_external = mpicom_ocn
+
+ ! Get file unit
+ nfu = shr_file_getUnit()
+
+ ! Get multiple instance data
+ inst_name = seq_comm_name(OCNID)
+ inst_index = seq_comm_inst(OCNID)
+ inst_suffix = seq_comm_suffix(OCNID)
+
+ ! ----------------------------------------------------------------
+ ! Initialize the model run
+ ! ----------------------------------------------------------------
+
+ call blom_cpl_indices_set()
+
+ call seq_infodata_GetData( infodata, case_name = runid_cesm )
+
+ call seq_infodata_GetData( infodata, start_type = starttype)
+
+ if (trim(starttype) == trim(seq_infodata_start_type_start)) then
+ runtyp_cesm = "initial"
+ elseif (trim(starttype) == trim(seq_infodata_start_type_cont) ) then
+ runtyp_cesm = "continue"
+ elseif (trim(starttype) == trim(seq_infodata_start_type_brnch)) then
+ runtyp_cesm = "branch"
+ else
+ write (lp,*) 'ocn_comp_mct ERROR: unknown starttype'
+ call shr_sys_flush(lp)
+ call shr_sys_abort()
+ endif
+
+ !-----------------------------------------------------------------
+ ! Get coupling time interval
+ !-----------------------------------------------------------------
+
+ call seq_timemgr_EClockGetData(EClock, dtime = ocn_cpl_dt_cesm)
+
+ ! ----------------------------------------------------------------
+ ! Initialize blom
+ ! ----------------------------------------------------------------
+
+ call t_startf('blom_init')
+ call blom_init
+ call t_stopf('blom_init')
+
+ ! ----------------------------------------------------------------
+ ! Reset shr logging to my log file
+ ! ----------------------------------------------------------------
+
+ call shr_file_getLogUnit (shrlogunit)
+ call shr_file_getLogLevel(shrloglev)
+ call shr_file_setLogUnit (lp)
+
+ call shr_sys_flush(lp)
+
+ ! ----------------------------------------------------------------
+ ! Check for consistency of BLOM calender information and EClock
+ ! ----------------------------------------------------------------
+
+ ! This must be completed!
+
+ if (runtyp_cesm == 'initial') then
+ call seq_timemgr_EClockGetData(EClock, &
+ start_ymd = start_ymd, &
+ start_tod = start_tod)
+ call shr_cal_date2ymd(start_ymd, start_year, start_month, start_day)
+ if (mnproc == 1) then
+ write (lp,'(a,i8,a2,i5,a2,i4.4,a1,i2.2,a1,i2.2)') &
+ ' cesm initial date: ',start_ymd,': ',start_tod,': ', &
+ start_year,'.',start_month,'.',start_day
+ call shr_sys_flush(lp)
+ endif
+ endif
+
+ ! ----------------------------------------------------------------
+ ! Initialize MCT attribute vectors and indices
+ ! ----------------------------------------------------------------
+
+ call t_startf ('blom_mct_init')
+
+ ! Initialize ocn gsMap
+
+ call ocn_SetGSMap_mct(mpicom_ocn, OCNID, gsMap_ocn)
+
+ ! Initialize mct ocn domain (needs ocn initialization info)
+
+ if (mnproc == 1) then
+ write (lp, *) 'blom: ocn_init_mct: lsize', lsize
+ endif
+
+ call domain_mct(gsMap_ocn, dom_ocn, lsize, perm, jjcpl)
+
+ ! Inialize mct attribute vectors
+
+ call mct_aVect_init(x2o_o, rList = seq_flds_x2o_fields, lsize = lsize)
+ call mct_aVect_zero(x2o_o)
+
+ call mct_aVect_init(o2x_o, rList = seq_flds_o2x_fields, lsize = lsize)
+ call mct_aVect_zero(o2x_o)
+
+ nsend = mct_avect_nRattr(o2x_o)
+ nrecv = mct_avect_nRattr(x2o_o)
+
+ !-----------------------------------------------------------------
+ ! Send intial state to driver
+ !-----------------------------------------------------------------
+
+ call getprecipfact_mct(lsend_precip_fact, precip_fact)
+ if ( lsend_precip_fact ) then
+ call seq_infodata_PutData( infodata, precip_fact=precip_fact)
+ endif
+ allocate(sbuff(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,nsend))
+ tlast_coupled = 0._r8
+ call get_So_t()
+ call sumsbuff_mct(nsend, sbuff, tlast_coupled)
+ call export_mct(o2x_o, lsize, perm, jjcpl, nsend, sbuff, tlast_coupled)
+ if (nreg == 2) then
+ call seq_infodata_PutData(infodata, ocn_prognostic = .true., &
+ ocnrof_prognostic = .true., &
+ ocn_nx = itdm , ocn_ny = jtdm-1)
+ else
+ call seq_infodata_PutData(infodata, ocn_prognostic = .true., &
+ ocnrof_prognostic = .true., &
+ ocn_nx = itdm , ocn_ny = jtdm)
+ endif
+
+ call t_stopf('blom_mct_init')
+
+
+ if (mnproc == 1) then
+ write (lp, *) 'blom: completed initialization!'
+ endif
+
+ !-----------------------------------------------------------------
+ ! Reset shr logging to original values
+ !-----------------------------------------------------------------
+
+ call shr_file_setLogUnit (shrlogunit)
+ call shr_file_setLogLevel(shrloglev)
+
+ call shr_sys_flush(lp)
+
+ end subroutine ocn_init_mct
+
+
+ subroutine ocn_run_mct(EClock, cdata_o, x2o_o, o2x_o)
+
+ ! Input/output arguments
+
+ type(ESMF_Clock), intent(inout) :: EClock
+ type(seq_cdata) , intent(inout) :: cdata_o
+ type(mct_aVect) , intent(inout) :: x2o_o
+ type(mct_aVect) , intent(inout) :: o2x_o
+
+ ! Local variables
+ type(seq_infodata_type), pointer :: infodata ! Input init object
+ integer :: shrlogunit, shrloglev, ymd, tod, ymd_sync, tod_sync
+
+ ! ----------------------------------------------------------------
+ ! Reset shr logging to my log file
+ ! ----------------------------------------------------------------
+
+ call shr_file_getLogUnit (shrlogunit)
+ call shr_file_getLogLevel(shrloglev)
+ call shr_file_setLogUnit (lp)
+
+ call seq_cdata_setptrs(cdata_o, infodata=infodata)
+
+ !-----------------------------------------------------------------
+ ! Advance the model in time over a coupling interval
+ !-----------------------------------------------------------------
+
+ blom_loop: do
+
+ if (nint(tlast_coupled) == 0) then
+ ! Obtain import state from driver
+ call import_mct(x2o_o, lsize, perm, jjcpl)
+ endif
+
+ ! Advance the model a time step
+ call blom_step
+ call get_So_t()
+ ! Add fields to send buffer sums
+ call sumsbuff_mct(nsend, sbuff, tlast_coupled)
+
+ if (nint(ocn_cpl_dt_cesm-tlast_coupled) == 0) then
+! call get_So_t(sbuff, ocn_cpl_dt_cesm)
+ ! Return export state to driver and exit integration loop
+ call export_mct(o2x_o, lsize, perm, jjcpl, nsend, sbuff, &
+ tlast_coupled)
+ exit blom_loop
+ endif
+
+ if (mnproc == 1) then
+ call shr_sys_flush(lp)
+ endif
+
+ enddo blom_loop
+
+ call getprecipfact_mct(lsend_precip_fact, precip_fact)
+ if ( lsend_precip_fact ) then
+ call seq_infodata_PutData( infodata, precip_fact=precip_fact)
+ endif
+
+ !-----------------------------------------------------------------
+ ! if requested, write restart file
+ !-----------------------------------------------------------------
+
+ if (seq_timemgr_RestartAlarmIsOn(EClock)) then
+ call restart_wt
+ endif
+
+ !-----------------------------------------------------------------
+ ! check that internal clock is in sync with master clock
+ !-----------------------------------------------------------------
+
+ if (mnproc == 1) then
+ call blom_time(ymd, tod)
+ if (.not. seq_timemgr_EClockDateInSync(EClock, ymd, tod )) then
+ call seq_timemgr_EClockGetData(EClock, curr_ymd=ymd_sync, &
+ curr_tod=tod_sync )
+ write(lp,*)' blom ymd=',ymd ,' blom tod= ',tod
+ write(lp,*)' sync ymd=',ymd_sync,' sync tod= ',tod_sync
+ call shr_sys_abort( 'ocn_run_mct'// &
+ ":: Internal blom clock not in sync with Sync Clock")
+ endif
+ endif
+
+ !-----------------------------------------------------------------
+ ! Reset shr logging to original values
+ !-----------------------------------------------------------------
+
+ call shr_file_setLogUnit (shrlogunit)
+ call shr_file_setLogLevel(shrloglev)
+
+ end subroutine ocn_run_mct
+
+
+ subroutine ocn_final_mct( EClock, cdata_o, x2o_o, o2x_o)
+
+ type(ESMF_Clock) , intent(inout) :: EClock
+ type(seq_cdata) , intent(inout) :: cdata_o
+ type(mct_aVect) , intent(inout) :: x2o_o
+ type(mct_aVect) , intent(inout) :: o2x_o
+
+ deallocate(perm)
+ deallocate(sbuff)
+
+ end subroutine ocn_final_mct
+
+
+ subroutine ocn_SetGSMap_mct(mpicom_ocn, OCNID, gsMap_ocn)
+
+ ! Input/output arguments
+
+ integer , intent(in) :: mpicom_ocn
+ integer , intent(in) :: OCNID
+ type(mct_gsMap), intent(inout) :: gsMap_ocn
+
+ ! Local variables
+
+ integer, allocatable :: gindex(:)
+ integer :: i, j, n, gsize
+
+ ! ----------------------------------------------------------------
+ ! Build the BLOM grid numbering for MCT
+ ! NOTE: Numbering scheme is: West to East and South to North
+ ! starting at south pole. Should be the same as what's used
+ ! in SCRIP
+ ! ----------------------------------------------------------------
+
+ if (nreg == 2 .and. nproc == jpr) then
+ jjcpl = jj - 1
+ else
+ jjcpl = jj
+ endif
+
+ lsize = ii*jjcpl
+
+ if (nreg == 2) then
+ gsize = itdm*(jtdm-1)
+ else
+ gsize = itdm*jtdm
+ endif
+
+ allocate(gindex(lsize))
+
+ n = 0
+ do j = 1, jjcpl
+ do i = 1, ii
+ n = n + 1
+ gindex(n) = (j0 + j - 1)*itdm + i0 + i
+ enddo
+ enddo
+
+ ! ----------------------------------------------------------------
+ ! reorder gindex to be in ascending order.
+ ! initialize a permutation array and sort gindex in-place
+ ! ----------------------------------------------------------------
+
+ allocate(perm(lsize))
+
+ call mct_indexset(perm)
+ call mct_indexsort(lsize, perm, gindex)
+ call mct_permute(gindex, perm, lsize)
+ call mct_gsMap_init(gsMap_ocn, gindex, mpicom_ocn, OCNID, lsize, gsize)
+
+ deallocate(gindex)
+
+ end subroutine ocn_SetGSMap_mct
+
+ subroutine get_So_t()
+ implicit none
+ ! Local variables
+ integer i, j, l, k, m, n, mm, nn, k1m, k1n
+ integer status,ncid,varid,start(3),count(3),stride(3)
+ integer :: nyear0,nyear1 ! model first year
+ integer :: ymd,tod,dimlen
+ integer :: nstrm,nmod,strm_year,nstep_strm
+ integer :: cplhist_year_align=1, &
+ strm_year_start=53, strm_year_end=56
+ real(r8), dimension(itdm,jtdm) :: tmp2d
+! real(r8), dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
+! intent(inout) :: So_t
+ character(len=100) :: filepath,filename,dimname
+! integer, intent(in) :: ocn_cpl_dt_cesm
+
+ if (mnproc == 1) then
+ nyear0=date0%year
+ nyear1=date%year
+ if (nyear0.ne.cplhist_year_align) then
+ write(*,*) 'So_t year is not aligned to model year'
+ endif
+ nstrm=strm_year_end-strm_year_start+1
+ nmod=mod(nyear1-nyear0,nstrm)
+ strm_year=strm_year_start+nmod
+ call blom_time(ymd,tod) ! tod = nint(mod(nstep, nstep_in_day)*baclin)
+ nstep_strm=int(time)-int(time0)-481434 !-int(nstep_strm/3650)*3650
+ start=(/nstep_strm,1,1/)
+ count=(/1,itdm,jtdm/)
+ stride=(/1,1,1/)
+ write(lp,*) start,count,stride
+ write(lp,*) 'time:',time,'time0',time0
+ !
+ filepath='/cluster/work/users/hra063/So_t/SST_new_mx_yr_36_50_orig.nc'
+ !
+ status=nf90_open(trim(filepath),nf90_nowrite,ncid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(4a)') 'nf90_open: sst_mx_orig.nc', &
+ nf90_strerror(status)
+ call xchalt('(get_So_t)')
+ stop '(get_So_t)'
+ endif
+
+ status=nf90_inq_varid(ncid,'sst',varid)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(4a)') 'nf90_inq_varid: So_t', &
+ nf90_strerror(status)
+ call xchalt('(get_So_t)')
+ stop '(get_So_t)'
+ endif
+ !
+ if (nstep_strm.eq.1) then
+ status=nf90_inquire_dimension(ncid, 1, dimname, dimlen)
+ write(lp,*) 'dim 1 name, size:',dimname,dimlen
+ status=nf90_inquire_dimension(ncid, 2, dimname, dimlen)
+ write(lp,*) 'dim 2 name, size:',dimname,dimlen
+ status=nf90_inquire_dimension(ncid, 3, dimname, dimlen)
+ write(lp,*) 'dim 3 name, size:',dimname,dimlen
+ endif
+ !
+ status=nf90_get_var(ncid,varid,tmp2d,start,count)
+ if (status.ne.nf90_noerr) then
+ write(lp,'(4a)') 'nf90_get_var: ','So_t',': ', &
+ nf90_strerror(status)
+ call xchalt('(get_So_t)')
+ stop '(get_So_t)'
+ endif
+ ! CLOSE
+! status=nf90_close(ncid)
+! if (status.ne.nf90_noerr) then
+! write(lp,'(4a)') 'nf90_close: ',trim(filename),': ', &
+! nf90_strerror(status)
+! call xchalt('(get_So_t)')
+! stop '(get_So_t)'
+! endif
+ endif
+ call xcaput(tmp2d,So_t,1)
+ if (mnproc == 1) then
+ write(lp,*) So_t(250,260)
+ endif
+ !if (mnproc == 1) then
+ !do j = 1, jtdm
+ !do i = 1, itdm
+ !write(lp,*) 'tmp2d(',i,j,'):',tmp2d(i,j)
+ !enddo
+ !enddo
+ !!write(lp,*) 'tmp2d()',shape(So_t2)
+ !endif
+ !(i0,j0)=(52,30)
+ !(i,j)=
+! do j = 1, jj
+! do l = 1, isp(j)
+! do i = max(1,ifp(j,l)), min(ii,ilp(j,l))
+! sbuff(i,j,index_o2x_So_t) = &
+! (1._r8-pcmask(i,j))*sbuff(i,j,index_o2x_So_t) + pcmask(i,j)*So_t(i,j)*ocn_cpl_dt_cesm
+! !sbuff(i,j,index_o2x_So_t) = So_t(i,j)*ocn_cpl_dt_cesm
+! enddo
+
+
+ end subroutine get_So_t
+
+
+end module ocn_comp_mct
diff --git a/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.drv/sumsbuff_mct.F b/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.drv/sumsbuff_mct.F
new file mode 100644
index 00000000..6e77d092
--- /dev/null
+++ b/cime_config/usermods_dirs/n1850_SST_pacemaker_Kuroshio/SourceMods/src.drv/sumsbuff_mct.F
@@ -0,0 +1,177 @@
+! ------------------------------------------------------------------------------
+! Copyright (C) 2008-2020 Mats Bentsen, Jerry Tjiputra
+!
+! This file is part of BLOM.
+!
+! BLOM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU Lesser General Public License as published by the Free
+! Software Foundation, either version 3 of the License, or (at your option)
+! any later version.
+!
+! BLOM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+! FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
+! more details.
+!
+! You should have received a copy of the GNU Lesser General Public License
+! along with BLOM. If not, see .
+! ------------------------------------------------------------------------------
+
+ subroutine sumsbuff_mct(nsend, sbuff, tlast_coupled)
+
+ ! Uses modules
+
+ use mod_types, only: r8
+ use mod_time, only: nstep, baclin, delt1, dlt, time0
+ use mod_xc
+ use mod_grid, only: scuy, scvx, scuxi, scvyi, pcmask, So_t
+ use mod_state, only: u, v, temp, saln, pbu, pbv, ubflxs, vbflxs,
+ . sealv
+ use mod_forcing, only: flxco2, flxdms
+ use mod_cesm, only: frzpot, ocn_cpl_dt_cesm
+ use blom_cpl_indices
+
+ implicit none
+
+ ! Input/output arguments
+
+ integer, intent(in) :: nsend
+ real(r8), dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,nsend),
+ . intent(inout) :: sbuff
+ real(r8), intent(inout) :: tlast_coupled
+
+ ! Local variables
+
+ integer i, j, l, k, m, n, mm, nn, k1m, k1n, nstep_strm
+ ! External functions
+
+ !-----------------------------------------------------------------
+ ! Set send buffer to zero if this is the first call after a
+ ! coupling interval
+ !-----------------------------------------------------------------
+
+ if (tlast_coupled == 0._r8) then
+ do k = 1, nsend
+ do j = 1-nbdy, jdm+nbdy
+ do i = 1-nbdy, idm+nbdy
+ sbuff(i,j,k) = 0._r8
+ enddo
+ enddo
+ enddo
+ endif
+
+ !-----------------------------------------------------------------
+ ! Accumulate fields in send buffer
+ !-----------------------------------------------------------------
+
+ m=mod(nstep+1,2)+1
+ n=mod(nstep ,2)+1
+ mm=(m-1)*kk
+ nn=(n-1)*kk
+ k1m=1+mm
+ k1n=1+nn
+
+ call xctilr(sealv, 1,1, 1,1, halo_ps)
+
+ do j = 1, jj
+ do l = 1, isu(j)
+ do i = max(1,ifu(j,l)), min(ii,ilu(j,l))
+ sbuff(i,j,index_o2x_So_u) =
+ . sbuff(i,j,index_o2x_So_u)
+ . + ( u(i,j,k1n)
+ . + (ubflxs(i,j,m) + ubflxs(i,j,n))*dlt
+ . /(pbu(i,j,n)*scuy(i,j)*delt1))*baclin
+ sbuff(i,j,index_o2x_So_dhdx) =
+ . sbuff(i,j,index_o2x_So_dhdx)
+ . + (sealv(i,j)-sealv(i-1,j))*scuxi(i,j)*baclin
+ enddo
+ enddo
+ enddo
+
+ do j = 1, jj
+ do l = 1, isv(j)
+ do i = max(1,ifv(j,l)), min(ii,ilv(j,l))
+ sbuff(i,j,index_o2x_So_v) =
+ . sbuff(i,j,index_o2x_So_v)
+ . + ( v(i,j,k1n)
+ . + (vbflxs(i,j,m) + vbflxs(i,j,n))*dlt
+ . /(pbv(i,j,n)*scvx(i,j)*delt1))*baclin
+ sbuff(i,j,index_o2x_So_dhdy) =
+ . sbuff(i,j,index_o2x_So_dhdy)
+ . + (sealv(i,j)-sealv(i,j-1))*scvyi(i,j)*baclin
+ enddo
+ enddo
+ enddo
+
+ do j = 1, jj
+ do l = 1, isp(j)
+ do i = max(1,ifp(j,l)), min(ii,ilp(j,l))
+! if ((plat(i,j).gt.35._r8.and.plat(i,j).lt.45._r8).and.
+! . (plon(i,j).gt.140._r8.and.plon(i,j).lt.160._r8)) then
+! if (pcmask(i,j).gt.0._r8) then
+ sbuff(i,j,index_o2x_So_t) =
+ . sbuff(i,j,index_o2x_So_t) +
+ . temp(i,j,k1n)*baclin +
+ . pcmask(i,j)*(So_t(i,j)-temp(i,j,k1n))*baclin
+ !write(lp,*) i+i0, j+j0, So_t(i,j)
+! else
+! sbuff(i,j,index_o2x_So_t) =
+! . sbuff(i,j,index_o2x_So_t) + temp(i,j,k1n)*baclin
+! endif
+ sbuff(i,j,index_o2x_So_s) =
+ . sbuff(i,j,index_o2x_So_s) + saln(i,j,k1n)*baclin
+ sbuff(i,j,index_o2x_Fioo_q) =
+ . sbuff(i,j,index_o2x_Fioo_q) + frzpot(i,j)
+ enddo
+ enddo
+ enddo
+
+! do j = 1, jj
+! do l = 1, isp(j)
+! do i = max(1,ifp(j,l)), min(ii,ilp(j,l))
+! sbuff(i,j,index_o2x_So_t) =
+! . sbuff(i,j,index_o2x_So_t) + temp(i,j,k1n)*baclin
+! sbuff(i,j,index_o2x_So_t) =
+! . sbuff(i,j,index_o2x_So_t) +
+! . max(1._r8-pcmask(i,j),0._r8)*temp(i,j,k1n)*baclin
+! . + max(pcmask(i,j),0._r8)*So_t(i,j)*baclin
+! sbuff(i,j,index_o2x_So_s) =
+! . sbuff(i,j,index_o2x_So_s) + saln(i,j,k1n)*baclin
+! sbuff(i,j,index_o2x_Fioo_q) =
+! . sbuff(i,j,index_o2x_Fioo_q) + frzpot(i,j)
+! enddo
+! enddo
+! enddo
+
+ if (index_o2x_Faoo_fco2_ocn > 0) then
+ do j = 1, jj
+ do l = 1, isp(j)
+ do i = max(1,ifp(j,l)), min(ii,ilp(j,l))
+ sbuff(i,j,index_o2x_Faoo_fco2_ocn) =
+ . sbuff(i,j,index_o2x_Faoo_fco2_ocn)
+ . + flxco2(i,j)*baclin
+ enddo
+ enddo
+ enddo
+ endif
+
+ if (index_o2x_Faoo_fdms_ocn > 0) then
+ do j = 1, jj
+ do l = 1, isp(j)
+ do i = max(1,ifp(j,l)), min(ii,ilp(j,l))
+ sbuff(i,j,index_o2x_Faoo_fdms_ocn) =
+ . sbuff(i,j,index_o2x_Faoo_fdms_ocn)
+ . + flxdms(i,j)*baclin
+ enddo
+ enddo
+ enddo
+ endif
+
+
+ !-----------------------------------------------------------------
+ ! Increment time since last coupling
+ !-----------------------------------------------------------------
+
+ tlast_coupled = tlast_coupled + baclin
+
+ end subroutine sumsbuff_mct