diff --git a/dev/compute_interpolation_weights/README.txt b/dev/compute_interpolation_weights/README.txt
new file mode 100644
index 0000000..0004b0e
--- /dev/null
+++ b/dev/compute_interpolation_weights/README.txt
@@ -0,0 +1,23 @@
+Scripts to create interpolation weights and distance to boundary files for interpolating forcing to RWPS mesh.
+To run:
+
+$ cd RWPS/dev/compute_interpolation_weights/ush
+$ sh compute_interpolation_weights.sh oc_1500m_30km
+
+to generate interpolation files for mesh rwps.oc_1500m_30km.msh.  Files will be writen to directory:
+RWPS/interpolation_weights
+
+Interpolation weights are created for:
+
+nbm oc domain
+rrfs hi domain (used for wind)
+rrfs pr domain (used for wind)
+rrfs ak domain (used for wind)
+rrfs na domain (used for wind)
+rrfs conus domain (used for wind)
+nbm ak domain (used for ice concentration)
+rtofs glo domain (used for current without extrapolation)
+rtofs glo domain (used for ice with extrapolation)
+stofs domain (used for waterlevel and current)
+
+Currently only setup to work on wcoss2.
diff --git a/dev/compute_interpolation_weights/ecf/compute_interpolation_weights.ecf b/dev/compute_interpolation_weights/ecf/compute_interpolation_weights.ecf
new file mode 100644
index 0000000..a257891
--- /dev/null
+++ b/dev/compute_interpolation_weights/ecf/compute_interpolation_weights.ecf
@@ -0,0 +1,26 @@
+#!/bin/bash
+#PBS -N compute_interpolation_weights
+#PBS -j oe
+#PBS -S /bin/bash
+#PBS -q dev
+#PBS -A NWPS-DEV
+#PBS -l walltime=04:00:00
+#PBS -l select=1:ncpus=32:mem=500GB
+#PBS -l place=excl
+#PBS -l debug=true
+
+module reset
+module load PrgEnv-intel/8.5.0
+module load intel/19.1.3.304
+module load craype/2.7.17
+module load cray-mpich/8.1.19
+module load hdf5-C/1.14.0
+module load netcdf-C/4.9.2
+module load esmf-C/8.6.0
+module load ve/hafs/2.1
+
+pip list -v
+
+$HOMErwps/dev/compute_interpolation_weights/jobs/COMPUTE_INTERPOLATION_WEIGHTS
+
+
diff --git a/dev/compute_interpolation_weights/jobs/COMPUTE_INTERPOLATION_WEIGHTS b/dev/compute_interpolation_weights/jobs/COMPUTE_INTERPOLATION_WEIGHTS
new file mode 100755
index 0000000..b9a65e1
--- /dev/null
+++ b/dev/compute_interpolation_weights/jobs/COMPUTE_INTERPOLATION_WEIGHTS
@@ -0,0 +1,135 @@
+#!/bin/bash
+
+######################################################################
+# This job computes interpolation weights and distance to boundary for 
+# a RWPS mesh. Currently configured for:
+# nbm oc domain (wind)
+# nbm ak domain (ice)
+# rrfs hi domain (wind)
+# rrfs pr domain (wind)
+# rrfs ak domain (wind)
+# rrfs conus domain (wind)
+# rrfs na domain (wind)
+# rtofs glo domain (surface currents, ice)
+# stofs v3 (waterlevel, currents)
+#
+# Different sets of interpolation weights are created for rtofs ice 
+# and surface currents so that ice concentration can be extrapolated
+# by nearest node while currents are not extrapolated into coastal 
+# regions outside of rtofs marine domain.
+#
+# Files are created in local RWPS/fix directory.
+######################################################################
+
+cd $HOMErwps/dev/compute_interpolation_weights/ush
+
+# Sample input directory
+SID="/lfs/h2/emc/couple/noscrub/keston.smith/SampleInput"
+date="20260730"
+cycl="00"
+Nprocs=32
+
+######################################################################
+# Compute rrfs interpolation weights for interpolating wind.
+# no extrapolation for rrfs forecasts.
+######################################################################
+rrfs_hi=$SID/wind.$date.$cycl/rrfs.$date.$cycl.wind10m.hi.nc
+rrfs_pr=$SID/wind.$date.$cycl/rrfs.$date.$cycl.wind10m.pr.nc
+rrfs_ak=$SID/wind.$date.$cycl/rrfs.$date.$cycl.wind10m.ak.nc
+rrfs_na=$SID/wind.$date.$cycl/rrfs.$date.$cycl.wind10m.na.nc
+rrfs_conus=$SID/wind.$date.$cycl/rrfs.$date.$cycl.wind10m.conus.nc
+
+
+rrfsIntrp="$out/InterpolationWeights.$meshname.rrfs"
+rrfsDst="$out/DistToBndy.$meshname.rrfs"
+
+rm -rf $rrfsIntrp.hi.nc $rrfsDst.hi.nc __pycache__
+python ./compute_gridded_to_rwps_interp_weights.py $rrfs_hi $mesh $rrfsIntrp.hi.nc $rrfsDst.hi.nc -1
+
+rm -rf $rrfsIntrp.pr.nc $rrfsDst.pr.nc __pycache__
+python ./compute_gridded_to_rwps_interp_weights.py $rrfs_pr $mesh $rrfsIntrp.pr.nc $rrfsDst.pr.nc -1
+
+rm -rf $rrfsIntrp.ak.nc $rrfsDst.ak.nc __pycache__
+python ./compute_gridded_to_rwps_interp_weights.py $rrfs_ak $mesh $rrfsIntrp.ak.nc $rrfsDst.ak.nc -1
+
+rm -rf $rrfsIntrp.na.nc $rrfsDst.na.nc __pycache__
+python ./compute_gridded_to_rwps_interp_weights.py $rrfs_na $mesh $rrfsIntrp.na.nc $rrfsDst.na.nc -1
+
+rm -rf $rrfsIntrp.conus.nc $rrfsDst.conus.nc __pycache__
+python ./compute_gridded_to_rwps_interp_weights.py $rrfs_conus $mesh $rrfsIntrp.conus.nc $rrfsDst.conus.nc -1
+
+######################################################################
+# Compute nbm interpolation weights on oc domain for interpolating wind.
+# Extrapolate to cover minor edge differences on south and east boundaries.
+######################################################################
+nbm_oc=$SID/wind.$date.$cycl/nbm.$date.$cycl.wind10m.oc.nc
+nbmIntrp="$out/InterpolationWeights.$meshname.nbm"
+nbmDst="$out/DistToBndy.$meshname.nbm"
+
+#extrapolate nbm oc wind to all points
+rm $nbmIntrp.oc.nc $nbmDst.oc.nc
+python ./compute_gridded_to_rwps_interp_weights.py $nbm_oc $mesh $nbmIntrp.oc.nc $nbmDst.oc.nc 1
+
+
+######################################################################
+# Compute nbm interpolation weights on ak domain for interpolating ice.
+# No extrapolation outside the ak region
+######################################################################
+nbm_ak=$SID/ice.$date.$cycl/nbm.$date.$cycl.ice.ak.nc
+#do not extrapolate for NBM .ak. domain
+rm $nbmIntrp.ak.nc $nbmDst.ak.nc
+python ./compute_gridded_to_rwps_interp_weights.py $nbm_ak $mesh $nbmIntrp.ak.nc $nbmDst.ak.nc -1
+
+rtofs_ice=$SID/ice.$date.$cycl/rtofs.ice.$date.nc
+rtofsIntrp="$out/InterpolationWeights.$meshname.rtofs"
+rtofsDst="$out/DistToBndy.$meshname.rtofs"
+
+######################################################################
+# Compute rtofs interpolation weights for interpolating ice.
+# Do extrapolate rtofs .glo. for full ice coverage in near shore.
+######################################################################
+rm $rtofsIntrp.ice.nc $rtofsDst.ice.nc 
+python ./compute_gridded_to_rwps_interp_weights.py $rtofs_ice $mesh $rtofsIntrp.ice.nc $rtofsDst.ice.nc 1 
+
+
+######################################################################
+# Compute rtofs interpolation weights for interpolating currents.
+# Do not extrapolate into coastal regions
+######################################################################
+rtofs_cur=$SID/rtofs.$date/currents.$date.rtofs.nc
+#do not extrapolate rtofs .glo. currents 
+rm $rtofsIntrp.current.nc $rtofsDst.current.nc
+python ./compute_gridded_to_rwps_interp_weights.py $rtofs_cur $mesh $rtofsIntrp.current.nc $rtofsDst.current.nc -1 
+
+
+######################################################################
+# Compute stofs interpolation weights for stofs v3 output
+# All dates after ~ August 18 2026
+######################################################################
+stofs_glo="$SID/stofs.v3/stofs.20260819.00/stofs_2d_glo.t00z.fields.cwl.nc"
+
+mkdir -p  STOFSInterpWeights.$meshname
+
+# First construct parallel jobscript to compute weights using Nprocs division of destination domain(meshfile)
+python compute_unstr_to_rwps_interp_weights.py $stofs_glo $mesh $Nprocs
+
+# Now run the parallel job script
+qsub -W block=true jobcardcompute_unstr_to_rwps_interp_weightsPBS
+wait
+
+# combine text file output of parallel jobs
+cat STOFSInterpWeights.$meshname/Part.IntrpWghts.*.txt > InterpWeights.$meshname.stofs.txt
+
+# sort order of nodes and convert output weights to netcdf file 
+python convert_weights_to_netcdf.py $stofs_glo $mesh InterpWeights.$meshname.stofs.txt 1
+mv InterpWeights.$meshname.stofs.nc $out/InterpolationWeights.$meshname.stofs.nc
+mv DistToBndy.$meshname.stofs.nc $out/DistToBndy.$meshname.stofs.nc
+
+# clean up temporary directory and files
+rm -rf STOFSInterpWeights.$meshname
+rm InterpWeights.$meshname.stofs.txt 
+rm *.out
+rm ESMPy*
+rm jobcardcompute_unstr_to_rwps*
+rm -rf __pycache__
+
diff --git a/dev/compute_interpolation_weights/ush/compute_gridded_to_rwps_interp_weights.py b/dev/compute_interpolation_weights/ush/compute_gridded_to_rwps_interp_weights.py
new file mode 100644
index 0000000..44b6157
--- /dev/null
+++ b/dev/compute_interpolation_weights/ush/compute_gridded_to_rwps_interp_weights.py
@@ -0,0 +1,213 @@
+# This is a a set of routines that interpolate wind forecasts to an unstructured 
+# mesh from regular or curvilinear grids (using ESMPY). In addition to interpolation 
+# field the out file contains a spatially variable error variance estimate based on the 
+# distance to the forecast boundary. This error variance is used later to update 
+# dispirate forecasts in a bayesian manner to yield a spatially smooth estimate 
+# incorporating all forecasts. The understanding is that the more local forecast 
+# products are of higher accuracy than the coarser broader scale forecasts
+
+import numpy as np
+import os
+import netCDF4 as nc
+import sys
+import xarray as xr
+import scipy.sparse as sp
+import compute_interpolation_weights_utilities as iutil
+
+# Main program
+AddExtrapolationSupport=True
+
+nargin = len(sys.argv) - 1
+
+flin=sys.argv[1]
+mshfl=sys.argv[2]
+weights_file=sys.argv[3]
+dist2bnd_file=sys.argv[4]
+# Don't use nearest neighbor interpolation unless 6th positive integer argument present.
+# This may be needed on boundary of RWPS mesh if node alignment is outside NBM OC domain
+Extrapolate=False
+if nargin > 4:
+    if int(sys.argv[5])>0:
+        print("using nearest neighbor to extrapolate wind field beyond geometric coverage")
+        Extrapolate=True
+
+xi, yi, ei, zi = iutil.loadWW3Mesh(mshfl)
+nn=len(xi)
+#shift coords to[129,370]
+data = nc.Dataset(flin,"r")
+
+#read spaital dimensions and determine if input mesh is curvilinear or regular
+x1=np.asarray(data["longitude"][:])
+y1=np.asarray(data["latitude"][:])
+
+#   rwps        :[-231, 11]
+#   nbm         :[129, 370]
+#   rrfs,pr     :[ -75.5000  -62.5087]
+#   rrfs,hi     :[ -161.5250 -153.8690]
+#   rrfs,ak     :[150.2012  266.2886]
+#   rrfs,na     :[67.5721  427.0000]
+#   rrfs,conus  :[ 225.9045  299.0828]
+#   rtofs,glo   :[74.1552,434.0146]
+#   stofs, glo  :[-180,180]
+
+#Shift longitude coordinates to RWPS specs for various files
+if "hi" in flin:
+    x1=x1+360.
+if "pr" in flin:
+    x1=x1+360.
+
+if "rtofs" in flin: #remove bad geometry edges
+    x1=x1[1:-1,1:-1]
+    y1=y1[1:-1,1:-1]
+    x1=x1-360.
+
+if (len(x1.shape)==2 and len(y1.shape)==2):
+    IsCrvLn=True
+elif (len(x1.shape)==1 and len(y1.shape)==1):
+    #represent regular grid as curvilinear grid
+    IsCrvLn=False
+    nx=len(x1)
+    ny=len(y1)
+    x1 = np.tile(x1,(ny,1))
+    y1 = np.tile(y1,(nx,1)).T
+    IsCrvLn=True
+else:
+    print("input file spatial dimension is not recognized. ending program")
+    sys.exit()
+    
+#######################################
+# === Create weights  ===#
+#######################################
+meshslash=mshfl.rfind('/')+1
+dom=flin.split(".")
+dom=dom[len(dom)-2]
+
+print("interpolation weights will be written to file = "+ weights_file)
+
+nx=x1.shape[0]
+ny=x1.shape[1]
+n1=nx*ny
+
+print("Computing weights and saving to file: "+ weights_file)
+
+#Use esmpy to construct bilinear interpolation weights
+iutil.CurvilinearGridCreateInterpWeights(xi, yi, x1, y1, weights_file)
+
+#######################################
+# === read weights and  ===#
+#######################################
+with xr.open_dataset(weights_file) as ds_s:
+   # Standard sparse storage uses 'row', 'col', and 'S'==weights variables
+   row = ds_s['row'].values
+   col = ds_s['col'].values
+   weights = ds_s['S'].values
+   Nrows=ds_s.attrs.get('Nrows')
+   Ncols=ds_s.attrs.get('Ncols')
+print("nn = "+str(nn)+": Nrows = "+str(Nrows))
+print("n1 = "+str(n1)+": Ncols = "+str(Ncols))
+if not ((nn==Nrows) and (n1==Ncols)):
+    print("Wrong matrix weights: number of rows from "+ mshfl +" = "+str(nn)+
+    " but number of rows in "+ weights_file +" = "+str(Nrows)+ 
+    ", number of spatial points in "+ flin +" = "+str(n1)+ 
+    " but number of columns in "+ weights_file +" = "+str(Ncols)  )
+    print("  You probably need to remove file "+ weights_file +" and rerun to generate appropriate weights")
+matrix = sp.coo_matrix((weights, (row-1, col-1)), shape=(nn,n1)).tocsr()
+print("sparse interpolation matrix")
+print(matrix)
+
+##################################################################################
+# START: Extrapolate for nodes not covered by interpolator
+##################################################################################
+x1v=np.transpose(x1).reshape(n1) # vectorize src nodes, consistant with data to interpolate
+y1v=np.transpose(y1).reshape(n1)
+
+if Extrapolate:
+    from scipy.interpolate import NearestNDInterpolator
+    srcp = np.array((x1v,y1v)).T
+    srcv = 1.+x1v**2 + y1v**2 #dummy input field
+    #dstv = matrix @ srcv.T
+    row_sum = matrix.sum(axis=1)
+    j0=np.where( row_sum==0 ) # destination nodes with no coverage from interpolation matrix
+    j0=np.array(j0[0]).tolist()
+    dstp = np.array((xi[j0],yi[j0])).T
+    interpolator = NearestNDInterpolator(srcp, srcv)
+    distances, j0src = interpolator.tree.query(dstp)
+    weightsExtrp=weights.tolist().append([1.0] * len(j0) )
+    j0=np.array(j0)
+    j0src=np.array(j0src)
+    rowExtrp=np.concatenate( (row, j0+1    ) )
+    colExtrp=np.concatenate( (col, j0src+1 ) )
+    weightsExtrp=np.concatenate( (weights, np.array([1.0] * len(j0))) )
+    os.replace(weights_file, weights_file[0:-3]+".NoExtrap.nc")
+    iutil.WriteInterpolationWeightsToNetCDF(weights_file,rowExtrp,colExtrp,weightsExtrp,len(xi),len(x1v))
+##################################################################################
+# FINISHED: Extrapolate for nodes not covered by interpolator
+##################################################################################
+
+##################################################################################
+# START: Extrapolation support for NaN occurances in source field
+##################################################################################
+if AddExtrapolationSupport:
+    with nc.Dataset(weights_file, 'r+', format='NETCDF4') as ncadd:
+        ncadd.createDimension('nn_src' , len(x1v))
+        ncadd.createDimension('nn_dst' , len(xi))
+        
+        xsrc_var=ncadd.createVariable('x_src', 'f8', ('nn_src',))
+        xsrc_var.long_name     = 'interpolation source node longitude'
+        xsrc_var[:]=x1v[:]
+        
+        ysrc_var=ncadd.createVariable('y_src', 'f8', ('nn_src',))
+        ysrc_var.long_name     = 'interpolation source node latitude'
+        ysrc_var[:]=y1v[:]
+        
+        xdst_var=ncadd.createVariable('x_dst', 'f8', ('nn_dst',))
+        xdst_var.long_name     = 'interpolation destination node longitude'
+        xdst_var[:]=xi[:]
+    
+        ydst_var=ncadd.createVariable('y_dst', 'f8', ('nn_dst',))
+        ydst_var.long_name     = 'interpolation destination node latitude'
+        ydst_var[:]=yi[:]
+        
+##################################################################################
+# FINISHED: Extrapolation support for NaN occurances in source field
+##################################################################################
+with nc.Dataset(weights_file, 'r+', format='NETCDF4') as ncadd:
+    ncadd.setncattr("SrcFieldType", "gridded")
+    ncadd.setncattr("InputFile", flin)
+    ncadd.setncattr("MeshFile", mshfl)
+
+
+##################################################################################
+# START: Compute distance to boundary for each node in mesh:
+##################################################################################
+row_sum = matrix.sum(axis=1)
+j0=np.where( row_sum==0 ) # destination nodes with no coverage from interpolation matrix
+j0=np.array(j0[0]).tolist()
+u0=np.ones(xi.shape)
+nan=float("nan")
+u0[j0]=nan
+
+if not j0: #all points are interpolation envelope
+    dist2bnd=0.*xi + np.inf
+else:
+    dist2bnd=iutil.CalculateDistanceToInterpEnvelope(xi,yi,u0, 1.)
+
+if Extrapolate:
+    dist2bnd=0.*dist2bnd + np.inf #all points are inside boundary- No boundary with this type of extrapolation
+    print(np.mean(dist2bnd))
+
+
+with nc.Dataset(dist2bnd_file, 'w', format='NETCDF4') as ncout:
+    ncout.createDimension('node' , nn)
+    d_var=ncout.createVariable('dist2bnd', 'f4', ('node',))
+    d_var.long_name     = 'distance to boundary'
+    d_var.units         = 'km'
+    d_var.standard_name = 'distance to boundary'
+    d_var[:]=dist2bnd[:]
+    
+    z_var=ncout.createVariable('depth', 'f4', ('node',))
+    z_var.long_name     = 'mesh depth'
+    z_var.units         = 'm'
+    z_var.standard_name = 'depth'
+    z_var[:]=zi[:]
+
diff --git a/dev/compute_interpolation_weights/ush/compute_interpolation_weights.sh b/dev/compute_interpolation_weights/ush/compute_interpolation_weights.sh
new file mode 100755
index 0000000..37c8517
--- /dev/null
+++ b/dev/compute_interpolation_weights/ush/compute_interpolation_weights.sh
@@ -0,0 +1,38 @@
+#!/bin/bash
+
+# --------------------------------------------------------------------------- #
+# Launch scripts to compute interpolation weights for an RWPS mesh. To run:   #
+#                                                                             #
+# $sh ./compute_interpolation_weights.sh oc_500m_10km                         #
+#                                                                             #
+# to generate interpolation weights for unstructured mesh                     #
+# rwps.oc_500m_10km.msh.                                                      #
+#                                                                             #
+# Last Changed : 09-09-2026                                                   #
+# --------------------------------------------------------------------------- #
+
+echo 'setting paths...'
+
+export meshID=$1
+
+export HOMErwps=$(cd "$(dirname "$(readlink -f -n "${BASH_SOURCE[0]}")")" && git rev-parse --show-toplevel)
+
+export fix="$HOMErwps/fix"
+
+export out=$HOMErwps/interpolation_weights
+
+mkdir -p $out
+
+#machine dependend path to RWPS fix files
+export RWPSfix=/lfs/h2/emc/couple/noscrub/keston.smith/RWPS
+
+# link mesh corresponding to meshID to local fix directory
+$HOMErwps/sorc/link_workflow.sh
+export mesh="$HOMErwps/fix/$meshID/rwps.$meshID.msh"
+
+meshname="${mesh##*/}"
+export meshname="${meshname: 0: -4}"
+
+#Retrieve current and process for forecast cycle
+qsub -V $HOMErwps/dev/compute_interpolation_weights/ecf/compute_interpolation_weights.ecf
+
diff --git a/dev/compute_interpolation_weights/ush/compute_interpolation_weights_utilities.py b/dev/compute_interpolation_weights/ush/compute_interpolation_weights_utilities.py
new file mode 100644
index 0000000..c66034a
--- /dev/null
+++ b/dev/compute_interpolation_weights/ush/compute_interpolation_weights_utilities.py
@@ -0,0 +1,458 @@
+import numpy as np
+import netCDF4 as nc
+
+def loadWW3Mesh(fl):
+    print("mesh file="+fl)
+    f=open(fl, 'r')
+    header = f.readline() 
+    header = f.readline() 
+    header = f.readline() 
+    header = f.readline() 
+    header = f.readline() # number of nodes
+    nn=int(header)
+    print("nn = "+str(nn))
+    xi=np.zeros(nn)
+    yi=np.zeros(nn)
+    zi=np.zeros(nn)
+    k=0
+    for i in range(nn):
+        A = f.readline()
+        B=A.lstrip()
+        values = B.split(" ")
+#        print(values)
+        if len(values)>5:
+            xi[k]=values[2]
+            yi[k]=values[4]
+            zi[k]=values[6]
+        else:
+            xi[k]=values[1]
+            yi[k]=values[2]
+            zi[k]=values[3]
+        k=k+1
+    print("number of nodes read: "+str(k))
+    header = f.readline() 
+    header = f.readline() 
+    header = f.readline() # number of elements
+    ne=int(header)#includes boundary nodes and actual elements
+    print("ne="+str(ne)+" -includes boundary nodes")
+    nbnd=0
+    bnd=[]
+    eix=np.zeros((ne,3), dtype=int)
+    k=0
+    for i in range(ne):
+        A = f.readline()
+        values = A.split(" ")
+        if len(values) == 6:
+            if int(values[2])==2:
+                bnd.append(int(values[5]))
+                nbnd=nbnd+1
+        if len(values)>15:
+            eix[k,0]=int(values[12])
+            eix[k,1]=int(values[14])
+            eix[k,2]=int(values[16])
+            k=k+1
+        elif len(values)>7:
+            eix[k,0]=int(values[6])
+            eix[k,1]=int(values[7])
+            eix[k,2]=int(values[8])
+            k=k+1
+    ei=eix[range(k),:]
+    print("number of open boundary nodes read: "+str(nbnd))
+    print("number of elements read: "+str(k))
+    return xi, yi, ei, zi
+
+############################################################################################
+# BEGIN WIND TO RWPS INTERP ROUTINES 
+import esmpy
+
+def CurvilinearGridCreateInterpWeights(xi,yi,x1,y1, weights_file):
+# Compute interpolation weights to interpolate from curvilinear grid (x1,y1) to points (xi,yi)
+# and store in netcdf file using ESMPY
+    debuging_output=False
+    
+    nx,ny=x1.shape
+    nn=len(xi)
+    n1=nx*ny
+    src_lon = x1
+    src_lat = y1
+    dst_lon=np.zeros((1,nn))
+    dst_lat=np.zeros((1,nn))
+    print("nn="+str(nn))
+    dst_lon[0,:] = xi[:]
+    dst_lat[0,:] = yi[:]
+
+    esmpy.Manager()
+
+    src_grid = esmpy.Grid(
+      max_index=np.array([src_lon.shape[0], src_lon.shape[1]]),
+      staggerloc=esmpy.StaggerLoc.CENTER,
+      coord_sys=esmpy.CoordSys.SPH_DEG
+    )
+
+    dst_grid = esmpy.Grid(
+      max_index=np.array([nn, 1]),
+      staggerloc=esmpy.StaggerLoc.CENTER,
+      coord_sys=esmpy.CoordSys.SPH_DEG
+    )
+
+    # 4. Populate the coordinate data
+    src_lon_ptr = src_grid.get_coords(0)
+    src_lat_ptr = src_grid.get_coords(1)
+    src_lon_ptr[...] = src_lon
+    src_lat_ptr[...] = src_lat
+
+    dst_lon_ptr = dst_grid.get_coords(0)
+    dst_lat_ptr = dst_grid.get_coords(1)
+    dst_lon_ptr[...] = dst_lon.T
+    dst_lat_ptr[...] = dst_lat.T
+
+    # 5. Create esmpy Fields
+    src_field = esmpy.Field(src_grid, name="src_field")
+    dst_field = esmpy.Field(dst_grid, name="dst_field")
+    src_field.data[...]=np.sqrt(np.abs(x1/180))/(90+y1) # arbitrary function of x,y
+    
+    if debuging_output: 
+        np.savetxt('F.txt', src_field.data[...])
+        np.savetxt('X.txt', x1)
+        np.savetxt('Y.txt', y1)
+    
+    print(f"Creating weights: {weights_file}")
+    regrid = esmpy.Regrid(
+      src_field,
+      dst_field,
+      filename=weights_file,
+      regrid_method=esmpy.RegridMethod.BILINEAR,
+      ignore_degenerate=True, # <--- Add this parameter
+      unmapped_action=esmpy.UnmappedAction.IGNORE # Optional: Ignores missing/masked points
+    )
+#Add number of rows and columns to weights file for clarity when constructing sparse matrix for interpolation
+    with nc.Dataset(weights_file, mode="a") as ds:
+        ds.Nrows = nn
+        ds.Ncols = n1
+    
+    if debuging_output: 
+        np.savetxt('Fi.txt', dst_field.data[...])
+        np.savetxt('xi.txt', xi)
+        np.savetxt('yi.txt', yi)
+    
+    return
+
+
+def CalculateDistanceToInterpEnvelope(xi,yi,fi,SearchWidth):
+# Alternative distance to boundary calculation for use when interpolation envelope is 
+# distinctly interior to curvilinear grid boundary as happens for RRFS NA grid 
+#
+# Inputs:
+#   xi (nn): longitude of unstructured mesh nodes 
+#   yi (nn): latitude of unstructured mesh nodes 
+#   fi (nn): interpolated field on mesh nodes with 'nan' values outside of interpolation envelope
+#   SearchWidth: Computational speed up to remove extra search points in distance to boundary
+#
+# Outputs:
+#   dist2bnd (nn) : distance to edge of interpolation envelope.  dist2bnd[k]=0 if (xi[k],yi[k]) is outside
+#                   of the interpolation envelope
+    nn=len(xi)
+    #dist2bnd=np.zeros(nn)
+    dist2bnd=np.full(nn,np.nan)
+    jin = np.where(~np.isnan(fi))[0].tolist()#points inside interpolation envelope
+    jout = np.where(np.isnan(fi))[0].tolist() #points outside interpolation envelope
+    xin=xi[jin]
+    yin=yi[jin]
+    xout=xi[jout]
+    yout=yi[jout]
+    jxU=np.where( xout < np.max(xin)+SearchWidth )[0].tolist()
+    jxD=np.where( xout > np.min(xin)-SearchWidth )[0].tolist()
+    jyU=np.where( yout < np.max(yin)+SearchWidth )[0].tolist()
+    jyD=np.where( yout > np.min(yin)-SearchWidth )[0].tolist()
+    j=list( set(jxU) & set(jxD)  & set(jyU)   & set(jyD)  )
+    xout=xout[j]
+    yout=yout[j]
+    din=np.zeros(len(jin))
+    print(len(xin))
+    print(len(xout))
+    for k in range(len(xin)):
+        din[k]=QuickDistance(yin[k],xin[k],yout,xout) # distance from node to closest point not interpolated to
+        if k%10000==0:
+            print("calculating distance to boundary, "+str(k)+":"+ str(nn)+":"+str(k/nn) )
+    dist2bnd[jin]=din
+    return dist2bnd
+
+def QuickDistance(lat1, lon1, lats2, lons2):
+    deg2kmY=111.
+    deg2kmX=np.cos( np.pi * lat1 / 180.)*deg2kmY
+    d= np.min(  np.sqrt( (  (lat1-lats2)*deg2kmY)**2 + ((lon1-lons2)*deg2kmX)**2 )  )
+    return d
+
+def WriteInterpolationWeightsToNetCDF(weights_file,row,col,weights,Nrows,Ncols):
+    #create a esmpy style sparse matrix netcdf file
+    print(Nrows)
+    print(Ncols)
+    print(row)
+    print(col)
+    print(weights)
+    n_s=len(weights)
+
+    if isinstance(row , list):
+        row      = np.array( row  )
+    if isinstance( col, list):
+        col      = np.array( col )
+    if isinstance( weights, list):
+        weights = np.array( weights )
+
+    with nc.Dataset(weights_file, 'w', format='NETCDF4') as ncout:
+            
+        ncout.createDimension('n_s' , n_s)
+        ncout.setncattr("Nrows", Nrows)
+        ncout.setncattr("Ncols", Ncols)
+            
+        r_var=ncout.createVariable('row', 'i4', ('n_s',))
+        r_var.long_name     = 'row index'
+        r_var[:]=row[:]
+            
+        c_var=ncout.createVariable('col', 'i4', ('n_s',))
+        c_var.long_name     = 'column index'
+        c_var[:]=col[:]
+            
+        s_var=ncout.createVariable('S', 'f4', ('n_s',))
+        s_var.long_name     = 'matrix value'
+        s_var[:]=weights[:]
+        
+
+
+def IsInElement(x,y,xp,yp):
+    IsIn=False
+    c1 = (x[1] - x[0]) * (yp - y[0]) - (y[1] - y[0]) * (xp - x[0])
+    c2 = (x[2] - x[1]) * (yp - y[1]) - (y[2] - y[1]) * (xp - x[1])
+    c3 = (x[0] - x[2]) * (yp - y[2]) - (y[0] - y[2]) * (xp - x[2])
+    if ( ( c1 > 0 and c2 > 0 and c3 > 0) or ( c1 < 0 and c2 < 0 and c3 < 0) ):
+        IsIn=True
+    if (c1*c2*c3 == 0): # include points on triangle
+        IsIn=True
+    return IsIn
+
+def FindElement(x,y,e,xi,yi):
+    nd=len(e.shape)
+    e=np.squeeze(e)
+    j=-9999
+    ne=e.shape[0]
+    if nd > 1: 
+        xc = np.squeeze(np.mean(x[e], axis=1))
+        yc = np.squeeze(np.mean(y[e], axis=1))
+        DistanceToElements = np.abs((xi + 1j*yi) - (xc + 1j*yc))
+        n=0
+        while (j < 0 and n < ne)  :
+            n=n+1
+            j = np.argmin( DistanceToElements )
+            xl = np.squeeze(x[e[j,:]])
+            yl = np.squeeze(y[e[j,:]])
+            IsIn=IsInElement(xl,yl,xi,yi)
+            if IsIn :
+                return j
+            else:
+                DistanceToElements[j]=float('inf')
+                j=-9999
+    else:
+        xl = np.squeeze(x[e])
+        yl = np.squeeze(y[e])
+        IsIn=IsInElement(xl,yl,xi,yi)
+        if IsIn:
+            j=0
+            return j
+        else:
+            j=-9999
+            return j
+        
+    return j
+
+
+def compute_mesh_to_mesh_interp_weights(x, y, e, xi, yi):
+    """
+    Compute mesh-to-mesh interpolation weights using barycentric coordinates.
+    
+    Parameters:
+    -----------
+    x : array-like
+        X coordinates of source mesh nodes
+    y : array-like
+        Y coordinates of source mesh nodes
+    e : array-like
+        Element connectivity matrix (n_elements x 3 for triangular elements)
+        (nodes indexed from 1 to nn)
+    xi : array-like
+        X coordinates of target points
+    yi : array-like
+        Y coordinates of target points
+        
+    Returns:
+    --------
+    weights : ndarray
+        Interpolation weights (n_points x 3)
+    nodes : ndarray
+        Node indices for each point (n_points x 3)
+    elenum : ndarray
+        Element number for each point (n_points,)
+    Dist2EleCenter : ndarray
+        Distance to element center for each element
+    """
+    
+    # Convert to numpy arrays
+    
+    deg2kmY=111.
+    UseNearestEle = False #if True just use closest element center
+    N = 12 # search N nearest elements (nearest by element center to target node distance)
+
+    x = np.asarray(x)
+    y = np.asarray(y)
+    e = np.asarray(e)
+    xi = np.asarray(xi)
+    yi = np.asarray(yi)
+
+    # convert node indexs to 0 .. nn-1
+    e0=e-1
+    
+    xc = np.mean(x[e0], axis=1)
+    yc = np.mean(y[e0], axis=1)
+    
+    # Initialize output arrays
+    n_points = len(xi)
+    n_elements = len(e)
+    weights = np.zeros((n_points, 3))
+    nodes = np.zeros((n_points, 3), dtype=int)
+    elenum = np.zeros(n_points, dtype=int)
+    Dist2EleCenter = np.zeros(n_points)
+    
+#    UseNearestEle=True
+    for k in range(n_points):
+        x0=xi[k]
+        y0=yi[k]
+        deg2kmX=np.cos( np.pi * y0 / 180.)*deg2kmY
+        distances = np.abs( deg2kmX*(x0 - xc) + 1j*deg2kmY*(y0 - yc))
+        if UseNearestEle:
+            jg=np.argmin(distances)
+        else:
+            raw_indices = np.argpartition(distances, N)[:N]
+            sorted_sub_indices = np.argsort(distances[raw_indices])
+            jg = raw_indices[sorted_sub_indices]
+        if UseNearestEle :
+            j=jg
+        else:
+            j=FindElement(x,y,e0[jg,:],x0,y0)
+            if j < 0:
+                print("couldn't find element for point : "+ str(k))
+                print(str(x0)+" : "+str(y0))
+                if len(jg)>0:
+                   j = jg[0]
+                else:
+                   j=0
+                print("using closest element at distance : "+ str(distances[j]))
+            else:
+                j=jg[j]
+
+        elenum[k] = j
+        Dist2EleCenter[k] = distances[j]
+        xl = x[e0[j, :]]
+        yl = y[e0[j, :]]
+            
+        # Compute barycentric coordinates using shoelace formula for areas
+        # Area of triangle formed by vertices 1, 2, and point (a3)
+        xt = np.array([xl[0], xl[1], x0, xl[0]])
+        yt = np.array([yl[0], yl[1], y0, yl[0]])
+        a3 = -np.dot(xt[1:4] - xt[0:3], yt[0:3] + yt[1:4]) / 2
+        # Area of triangle formed by vertices 3, 1, and point (a2)
+        xt = np.array([xl[2], xl[0], x0, xl[2]])
+        yt = np.array([yl[2], yl[0], y0, yl[2]])
+        a2 = -np.dot(xt[1:4] - xt[0:3], yt[0:3] + yt[1:4]) / 2
+        # Area of triangle formed by point, vertices 2, 3 (a1)
+        xt = np.array([x0, xl[1], xl[2], x0])
+        yt = np.array([y0, yl[1], yl[2], y0])
+        a1 = -np.dot(xt[1:4] - xt[0:3], yt[0:3] + yt[1:4]) / 2
+
+        # Normalize to get barycentric weights
+        total_area = a1 + a2 + a3
+        weights[k, :] = [a1, a2, a3] / total_area
+
+        nodes[k, :] = e[j, :] #<- indexed 1 .. nn
+        # Progress reporting every 100 iterations
+        if (k + 1) % 100 == 0:
+            print(f"Progress: {k+1}/{n_points}")
+    return weights, nodes, elenum, Dist2EleCenter
+
+
+
+# Move to the directory where the job was submitted
+#-->cd $PBS_O_WORKDIR
+#-->echo "Running on node: $(hostname)"
+#--> python my_script.py --task $PBS_ARRAY_INDEX
+
+import os
+def WriteInterpJobscriptPBS(fl,flin,mshfl,Njobs, ComputeNodes):
+    
+    meshslash=mshfl.rfind('/')+1
+    TmpOutDir="STOFSInterpWeights."+mshfl[meshslash:len(mshfl)-4]
+    WghtFl="STOFS.wght."+mshfl[meshslash:len(mshfl)-4]+".txt"
+    WghtFl="InterpWeights."+mshfl[meshslash:len(mshfl)-4]+".stofs.txt"
+    WghtFlNetCDF="InterpWeights."+mshfl[meshslash:len(mshfl)-4]+".stofs.nc"
+    Njobs=32 #OVERWRITE FOR NOW
+    with open(fl, 'w') as f:
+        f.write("#PBS -N ESMPy\n")
+        f.write("#PBS -j oe\n")
+        f.write("#PBS -S /bin/bash\n")
+        f.write("#PBS -q dev\n")
+        f.write("#PBS -A NWPS-DEV\n")
+        f.write("#PBS -l walltime=01:00:00\n")
+        f.write("#PBS -J 1-"+str(Njobs)+"\n")
+        f.write("#PBS -l select=2:ncpus=32:mem=128gb\n")
+        f.write("#PBS -l place=excl\n")
+        f.write("#PBS -l debug=true\n")
+        f.write("#PBS -r y\n")
+
+        f.write("module reset\n")
+        f.write("module load PrgEnv-intel/8.5.0\n")
+        f.write("module load intel/19.1.3.304\n")
+        f.write("module load craype/2.7.17\n")
+        f.write("module load cray-mpich/8.1.19\n")
+        f.write("module load hdf5-C/1.14.0\n")
+        f.write("module load netcdf-C/4.9.2\n")
+        f.write("module load esmf-C/8.6.0\n")
+        f.write("module load ve/hafs/2.1\n")
+
+        f.write("pip list -v\n")
+
+        current_dir = os.getcwd()
+        f.write("cd "+current_dir+"\n")
+
+        f.write("# calculate interpolation weights in parallel geographically \n")
+        f.write("python compute_unstr_to_rwps_interp_weights.py "+flin+" "+mshfl+" $PBS_ARRAY_INDEX " + str(Njobs)+" > InterpJob.$PBS_ARRAY_INDEX.out \n")
+
+
+def WriteInterpJobscriptSLURM(fl,flin,mshfl,Njobs, ComputeNodes):
+    
+    meshslash=mshfl.rfind('/')+1
+    TmpOutDir="STOFSInterpWeights."+mshfl[meshslash:len(mshfl)-4]
+    WghtFl="STOFS.wght."+mshfl[meshslash:len(mshfl)-4]+".txt"
+    WghtFlNetCDF="STOFS.wght."+mshfl[meshslash:len(mshfl)-4]+".nc"
+    with open(fl, 'w') as f:
+        f.write("#!/bin/bash \n")
+        f.write("#SBATCH --job-name=STOFS_interp_masterscript \n")
+        f.write("#SBATCH --ntasks=1 \n") # ntasks per interpolation
+        f.write("#SBATCH --time=08:00:00 \n") 
+        f.write("#SBATCH --output=mpi_test_%j.log \n")
+        f.write("#SBATCH --error=%j.err \n")
+        f.write("#SBATCH --account=marine-cpu \n")
+        f.write("#SBATCH --nodes="+str(ComputeNodes)+" \n")
+        f.write("#SBATCH --ntasks-per-core=1"+" \n")
+        f.write("#SBATCH --array=0-"+str(Njobs-1)+" \n")
+
+        f.write(" \n")
+
+        f.write("module purge \n")
+        f.write("module use /scratch4/NCEPDEV/marine/Ali.Salimi/Hera_Data/HR4-OPT/FromJessica/Keston/ICunstructuredRuns15km-implicit-450s/global-workflow/sorc/ufs_model.fd/modulefiles \n")
+        f.write("module load ufs_ursa.intel \n")
+        f.write("module load py-scipy/1.14.1 \n")
+        f.write("module load py-netcdf4/1.7.1.post2 \n")
+        f.write("pip list \n")
+        current_dir = os.getcwd()
+        f.write("cd "+current_dir+"\n")
+
+        f.write("# calculate interpolation weights in parallel geographically \n")
+        f.write("python compute_unstr_to_rwps_interp_weights.py "+flin+" "+mshfl+" $SLURM_ARRAY_TASK_ID " + str(Njobs)+" > InterpJob.$SLURM_ARRAY_TASK_ID.out \n")
diff --git a/dev/compute_interpolation_weights/ush/compute_unstr_to_rwps_interp_weights.py b/dev/compute_interpolation_weights/ush/compute_unstr_to_rwps_interp_weights.py
new file mode 100644
index 0000000..356fccd
--- /dev/null
+++ b/dev/compute_interpolation_weights/ush/compute_unstr_to_rwps_interp_weights.py
@@ -0,0 +1,146 @@
+import numpy as np
+import netCDF4 as nc
+import sys
+import compute_interpolation_weights_utilities as  iutil
+import os
+
+PBS=True
+
+nargin = len(sys.argv) - 1
+flin=sys.argv[1]
+mshfl=sys.argv[2]
+
+meshslash=mshfl.rfind('/')+1
+TmpOutDir="STOFSInterpWeights."+mshfl[meshslash:len(mshfl)-4]
+
+#if nargin==3 then write the parallel jobcard that carries out the interpolation weight calculation
+if nargin ==3 :
+    Njobs=int(sys.argv[3])
+    os.makedirs(TmpOutDir, exist_ok=True)
+    try:
+        os.remove(TmpOutDir+"/*.txt")
+    except:
+        print("directory "+TmpOutDir+" is alread empty")
+    iutil.WriteInterpJobscriptPBS("jobcardcompute_unstr_to_rwps_interp_weightsPBS",flin,mshfl,Njobs, 1)
+    iutil.WriteInterpJobscriptSLURM("jobcardcompute_unstr_to_rwps_interp_weightsSLURM",flin,mshfl,Njobs, 1)
+    print("Made parallel jobcards to create interpolation weights with "+str(Njobs)+" processes. Next step:")
+    print("sbatch jobcardcompute_unstr_to_rwps_interp_weightsSLURM")
+    sys.exit()
+
+#if nargin==4 then do the part of the mesh for jobID of Njobs
+if nargin == 4:
+    jobID=int(sys.argv[3])
+    Njobs=int(sys.argv[4])
+
+#if nargin==6 then do the interpolation weight calculation for a specified window
+else:
+    lonW=int(sys.argv[3])
+    lonE=int(sys.argv[4])
+    latS=int(sys.argv[5])
+    latN=int(sys.argv[6])
+
+meshslash=mshfl.rfind('/')+1
+if nargin <  5:
+    weights_file = TmpOutDir+"/Part.IntrpWghts."+str(jobID)+".txt"
+else:
+    weights_file = TmpOutDir+"/Part.IntrpWghts.W"+str(lonW)+".E"+str(lonE)+".S"+str(latS)+".N"+str(latN)+"."+mshfl[meshslash:len(mshfl)-3]+"txt"
+
+
+print("saving output to:"+weights_file)
+
+xi, yi, ei, zi = iutil.loadWW3Mesh(mshfl)
+nni=len(xi)
+j=np.where(xi>90.) # should be empty
+xi[j]=xi[j]-360.
+
+# divide domain from east west based on job id=0 ... Njobs
+# and make North-South window contain the full target domain
+# NOTE: for efficiency the sections of the domain should have approximately
+# the same number of destination nodes and/or the same number of source elements
+
+if nargin <  5:
+# balance node load 
+# on PBS system jobID runs from 1 to Njobs
+# on SLURM system jobID runs from 0 to Njobs-1 so this partition scheam needs to be modified
+    xis=np.sort(xi)
+    mm=round(nni/Njobs)
+    xil=xis[range(0,nni,mm)] # geographic partitions for parallelization
+
+    if PBS:
+        lonW=xil[jobID-1]
+        if jobID==Njobs:
+            lonE=xis[-1]+1.
+        else:
+            lonE=xil[jobID]
+        if jobID==1:
+            lonW=xis[0]-1.
+    else:
+        lonW=xil[jobID]
+        if jobID==Njobs-1:
+            lonE=xis[-1]+1.
+        else:
+            lonE=xil[jobID+1]
+        if jobID==0:
+            lonW=xis[0]-1.
+
+    latS=np.min(yi)-1.
+    latN=np.max(yi)+1.
+
+data = nc.Dataset(flin,"r")
+#read spaital dimensions and determine if input mesh is curvilinear or regular
+x=np.asarray(data["x"][:])
+j=np.where(x>90.) #elements broken here should not effect RWPS interpolation
+x[j]=x[j]-360.
+
+y=np.asarray(data["y"][:])
+e=np.asarray(data["element"][:,:])
+
+e0=e-1
+xc = np.mean(x[e0], axis=1)
+yc = np.mean(y[e0], axis=1)
+
+nn=len(x)
+
+SearchWidth=.5 # search rectangle extended beyond window to this number (degrees lat lon)
+
+print(e.shape)
+ne=e.shape[0]
+print(ne)
+
+xwin=[lonW,lonE]
+ywin=[latS,latN]
+print("calculating interpolation weights for [W"+str(lonW)+": E "+str(lonE)+": S "+str(latS)+": N "+str(latN)+"]")
+#Find target mesh nodes in window
+jxUi=np.where( xi < np.max(xwin))[0].tolist()
+jxDi=np.where( xi >= np.min(xwin))[0].tolist()
+jyUi=np.where( yi <= np.max(ywin))[0].tolist()
+jyDi=np.where( yi >= np.min(ywin))[0].tolist()
+ji=list( set(jxUi) & set(jxDi)  & set(jyUi)   & set(jyDi)  )
+
+#Find source elements in and near window
+SearchWidth=.5
+jxU=np.where( xc < np.max(xwin)+SearchWidth )[0].tolist()
+jxD=np.where( xc > np.min(xwin)-SearchWidth )[0].tolist()
+jyU=np.where( yc < np.max(ywin)+SearchWidth )[0].tolist()
+jyD=np.where( yc > np.min(ywin)-SearchWidth )[0].tolist()
+je=list( set(jxU) & set(jxD)  & set(jyU)   & set(jyD)  )
+
+weights, nodes, elenum, Dist2EleCenter=iutil.compute_mesh_to_mesh_interp_weights(x, y, e[je,:], xi[ji], yi[ji])
+
+#Move back to global index's
+for k in range(len(elenum)):
+    elenum[k]=je[elenum[k]]
+    ji[k]=ji[k]
+
+#return to convention numbering from 1 ...
+elenum = list(np.array(elenum) + 1)
+ji = list(np.array(ji) + 1)
+#nodes = np.array(nodes) + 1 # as writen using the correct 1 .. nn node numbering 
+
+# File format 
+# ["taget node number" "source element number" "source n1" "source n2 "source n3" "distance from target node to element center" "weight1" "weight2" "weight3"]
+
+Fout=np.vstack((np.array(ji),elenum,nodes[:,0],nodes[:,1],nodes[:,2],Dist2EleCenter,weights[:,0],weights[:,1],weights[:,2]))
+np.savetxt(weights_file,Fout.T, fmt='%d %d %d %d %d %.6f %.6f %.6f %.6f')
+
+
diff --git a/dev/compute_interpolation_weights/ush/convert_weights_to_netcdf.py b/dev/compute_interpolation_weights/ush/convert_weights_to_netcdf.py
new file mode 100644
index 0000000..6d346e6
--- /dev/null
+++ b/dev/compute_interpolation_weights/ush/convert_weights_to_netcdf.py
@@ -0,0 +1,134 @@
+import numpy as np
+import netCDF4 as nc
+import sys
+import compute_interpolation_weights_utilities as  iutil
+import scipy.sparse as sp
+
+nargin = len(sys.argv) - 1
+
+flin=sys.argv[1]
+mshfl=sys.argv[2]
+meshslash=mshfl.rfind('/')+1
+
+AddExtrapolationSupport=True
+
+TextWeightFl=sys.argv[3]
+
+flout=TextWeightFl[0:len(TextWeightFl)-4]+".nc"
+
+xi, yi, ei, zi = iutil.loadWW3Mesh(mshfl)
+nn_dst=len(xi)
+data = nc.Dataset(flin,"r")
+#read spaital dimensions and determine if input mesh is curvilinear or regular
+x=np.asarray(data["x"][:])
+#convert to RWPS coordinates
+
+jEast=np.where(x>90)
+x[jEast]=x[jEast]-360.
+
+nn_src=len(x)
+
+#Read in weights from cat of text output(rows are not in order)
+F=np.loadtxt(TextWeightFl)
+
+node_dst=F[:,0] # row number (destination node)
+ele_num_src=F[:,1] # Source element number (not used)
+node_src=F[:,2:5] # Source node
+DistToCenterKM=F[:,5] # distance of destination node from source element center (not used) 
+weights=F[:,6:9] # interpolation weights from source nodes to destination node
+
+node_dst=node_dst.astype(int)
+node_src=node_src.astype(int)
+ele_num_src=ele_num_src.astype(int)
+
+nnz=F.shape[0]
+n_s=3*nnz
+row=np.zeros(n_s,dtype=int)
+col=np.zeros(n_s,dtype=int)
+val=np.zeros(n_s)
+
+n=0
+for k in range(nnz):
+    for j in range(3):
+        row[n]=node_dst[k]
+        col[n]=node_src[k,j]
+        val[n]=weights[k,j]
+        n=n+1
+print("n="+str(n))
+print("n_s="+str(n_s))
+n_s=n
+
+with nc.Dataset(flout, 'w', format='NETCDF4') as ncout:
+    ncout.createDimension('n_s' , n_s)
+    ncout.setncattr("Nrows", nn_dst)
+    ncout.setncattr("Ncols", nn_src)
+    ncout.setncattr("SrcFieldType", "unstructured") 
+    
+    r_var=ncout.createVariable('row', 'i4', ('n_s',))
+    r_var.long_name     = 'row index'
+    r_var[:]=row[:]
+    
+    c_var=ncout.createVariable('col', 'i4', ('n_s',))
+    c_var.long_name     = 'column index'
+    c_var[:]=col[:]
+    
+    s_var=ncout.createVariable('S', 'f8', ('n_s',))
+    s_var.long_name     = 'matrix value'
+    s_var[:]=val[:]
+
+#Consider adding (x,y) destination and (x,y) for source. This is usefull for extrapolation
+    if AddExtrapolationSupport:
+        x=np.asarray(data["x"][:])
+        y=np.asarray(data["y"][:])
+        ncout.createDimension('nn_src' , nn_src)
+        ncout.createDimension('nn_dst' , nn_dst)
+        
+        xsrc_var=ncout.createVariable('x_src', 'f8', ('nn_src',))
+        xsrc_var.long_name     = 'interpolation source node longitude'
+        xsrc_var[:]=x[:]
+        
+        ysrc_var=ncout.createVariable('y_src', 'f8', ('nn_src',))
+        ysrc_var.long_name     = 'interpolation source node latitude'
+        ysrc_var[:]=y[:]
+        
+        xdst_var=ncout.createVariable('x_dst', 'f8', ('nn_dst',))
+        xdst_var.long_name     = 'interpolation destination node longitude'
+        xdst_var[:]=xi[:]
+    
+        ydst_var=ncout.createVariable('y_dst', 'f8', ('nn_dst',))
+        ydst_var.long_name     = 'interpolation destination node latitude'
+        ydst_var[:]=yi[:]
+
+        zdst_var=ncout.createVariable('depth_dst', 'f8', ('nn_dst',))
+        zdst_var.long_name     = 'interpolation destination node latitude'
+        zdst_var[:]=zi[:]
+
+Extrapolate=sys.argv[4]
+
+if Extrapolate:
+    dist2bnd=np.full(len(xi), np.inf) #all points are inside boundary- No boundary with this type of extrapolation
+else:
+    matrix = sp.coo_matrix((weights, (row-1, col-1)), shape=(nn_src,nn_dst)).tocsr()
+    row_sum = matrix.sum(axis=1)
+    j0=np.where( row_sum==0 ) # destination nodes with no coverage from interpolation matrix
+    j0=np.array(j0[0]).tolist()
+    u0=np.ones(xi.shape)
+    nan=float("nan")
+    u0[j0]=nan
+    dist2bnd=iutil.CalculateDistanceToInterpEnvelope(xi,yi,u0, 1.)
+
+dist2bnd_file = "DistToBndy."+mshfl[meshslash:len(mshfl)-4]+".stofs.nc"
+with nc.Dataset(dist2bnd_file, 'w', format='NETCDF4') as ncout:
+    ncout.createDimension('node' , nn_dst)
+    d_var=ncout.createVariable('dist2bnd', 'f4', ('node',))
+    d_var.long_name     = 'distance to boundary'
+    d_var.units         = 'km'
+    d_var.standard_name = 'distance to boundary'
+    d_var[:]=dist2bnd[:]
+    
+    z_var=ncout.createVariable('depth', 'f4', ('node',))
+    z_var.long_name     = 'mesh depth'
+    z_var.units         = 'm'
+    z_var.standard_name = 'depth'
+    z_var[:]=zi[:]
+
diff --git a/sorc/get_externals.sh b/sorc/get_externals.sh
index 47b172e..b48d8f6 100755
--- a/sorc/get_externals.sh
+++ b/sorc/get_externals.sh
@@ -7,7 +7,18 @@
 # --------------------------------------------------------------------------- #
 
 echo 'Fetching externals...'
-cp -p /lfs/h2/emc/couple/noscrub/keston.smith/meshes/RWPS.v0.msh ../fix/
+
+
+# copy mesh to local fix directory
+cp -p $RWPSfix/fix/$meshname/20260722/rwps.$meshname.msh ../fix/
+# copy Interpoplation weights for nbm, rrfs, rtofs and stofs to local fix directory
+cp -p $RWPSfix/fix/$meshname/20260722/InterpolationWeights*$meshname*.nc ../fix/
+# copy distance to boundary for nbm, rrfs, rtofs and stofs to local fix directory
+cp -p $RWPSfix/fix/$meshname/20260722/DistToBndy*$meshname*.nc ../fix/
+
+
+#cp -p /lfs/h2/emc/couple/noscrub/keston.smith/meshes/RWPS.v0.msh ../fix/
+
 cp -p /lfs/h2/emc/couple/noscrub/saeideh.banihashemi/git/fv3gfs/fix/fix_glwu/mesh.glwu ../fix/mesh.rwps
 cp -p /lfs/h2/emc/couple/noscrub/saeideh.banihashemi/git/fv3gfs/fix/fix_glwu/grint_weights.grlc_2p5km ../fix/
 cp -p /lfs/h2/emc/couple/noscrub/saeideh.banihashemi/git/fv3gfs/fix/fix_glwu/grint_weights.grlr ../fix/
diff --git a/versions/preproc.wcoss2.ver b/versions/preproc.wcoss2.ver
new file mode 100755
index 0000000..a53605b
--- /dev/null
+++ b/versions/preproc.wcoss2.ver
@@ -0,0 +1,12 @@
+#! /usr/bin/env bash
+export PrgEnv_intel_ver=8.5.0
+export interl_ver=19.1.3.304
+export craype_ver=2.7.17
+export cray_mpich_ver=8.1.19
+export hdf5_C_ver=1.14.0
+export netcdf_C_ver=4.9.2
+export esmf_C_ver=8.6.0
+export ve_hafs_ver=2.1
+
+export intel_oneapi_ver=2022.2.0.262
+export wgrib2_ver=2.0.8