diff --git a/examples/GEOSInput/Antarctica_input_control.py b/examples/GEOSInput/Antarctica_input_control.py
new file mode 100755
index 000000000..f4e7b9b31
--- /dev/null
+++ b/examples/GEOSInput/Antarctica_input_control.py
@@ -0,0 +1,103 @@
+import devpath
+import os
+ISSM_DIR = os.getenv('ISSM_DIR') # for binaries
+
+import numpy as np
+from triangle import triangle
+from model import *
+from netCDF4 import Dataset
+from InterpFromGridToMesh import InterpFromGridToMesh
+from bamg import bamg
+from loadmodel import loadmodel
+from setmask import setmask
+from parameterize import parameterize
+from clusters.discover_geos import export_discover
+from marshall import marshall
+
+if not os.path.exists('./Models'):
+ os.mkdir('./Models')
+
+# Step 1: Mesh generation
+# Generate initial uniform mesh (resolution = 60000 m)
+# project mesh onto new coordinate system
+md = triangle(model(), '../Data/Ais.exp', 60000)
+
+print(' Loading velocities data from NetCDF')
+nsidc_vel = Dataset('../Data/Antarctica_ice_velocity.nc')
+xmin = nsidc_vel.xmin
+xmin = float(xmin.lstrip()[0:10])
+ymax = nsidc_vel.ymax
+ymax = float(ymax.lstrip()[0:10])
+spacing = nsidc_vel.spacing
+spacing = float((spacing.lstrip())[0:4])
+nx = nsidc_vel.nx
+ny = nsidc_vel.ny
+velx = nsidc_vel['vx'][:].data
+vely = nsidc_vel['vy'][:].data
+# Build coordinates
+x2 = xmin + np.arange(nx + 1) * spacing
+y2 = (ymax - ny * spacing) + np.arange(ny + 1) * spacing
+
+# print(' Set observed velocities')
+md.initialization.vx = InterpFromGridToMesh(x2, y2, np.flipud(velx), md.mesh.x, md.mesh.y, 0)
+md.initialization.vy = InterpFromGridToMesh(x2, y2, np.flipud(vely), md.mesh.x, md.mesh.y, 0)
+md.initialization.vz = np.zeros(md.mesh.numberofvertices)
+md.initialization.vel = np.sqrt(md.initialization.vx**2 + md.initialization.vy**2)
+del velx, vely
+
+md = bamg(md, 'hmax', 75000, 'hmin', 5000, 'gradation', 1.4, 'field', md.initialization.vel, 'err', 8)
+
+# export mesh
+export_discover(md, './Models/Antarctica_mesh.nc')
+
+# Step 2: parameterize model
+md = loadmodel('./Models/Antarctica_mesh.nc')
+
+md = setmask(md, '', '')
+md = parameterize(md, './Antarctica_parameterize.py')
+
+
+# export parameterization
+export_discover(md, "./Models/Antarctica_parameterization.nc",delete_rundir=True)
+
+# Step 3: basal friction inversion
+ # Control general
+md.inversion.nsteps = 400
+md.inversion.iscontrol=1
+md.inversion.maxsteps=400
+md.inversion.maxiter=400
+md.inversion.dxmin=0.01
+md.inversion.gttol=1.0e-8
+
+md.inversion.step_threshold = 0.99 * np.ones((md.inversion.nsteps))
+md.inversion.maxiter_per_step = 40 * np.ones((md.inversion.nsteps))
+
+md.inversion.gradient_scaling = 50 * np.ones((md.inversion.nsteps, 1))
+md.inversion.min_parameters = 1 * np.ones((md.mesh.numberofvertices, 1))
+md.inversion.max_parameters = 200 * np.ones((md.mesh.numberofvertices, 1))
+
+#Cost functions
+md.inversion.cost_functions = [101, 103, 501]
+md.inversion.cost_functions_coefficients = np.ones((md.mesh.numberofvertices, 3))
+md.inversion.cost_functions_coefficients[:, 0] = 1
+md.inversion.cost_functions_coefficients[:, 1] = 1
+md.inversion.cost_functions_coefficients[:, 2] = 2e-10
+
+# Controls
+md.inversion.control_parameters = ['FrictionCoefficient']
+md.inversion.min_parameters=1*np.ones(np.shape(md.mesh.x))
+md.inversion.max_parameters=200*np.ones(np.shape(md.mesh.x))
+
+# Additional parameters
+md.stressbalance.restol=0.0000001
+md.stressbalance.reltol=0.0000001
+md.stressbalance.abstol=np.nan
+
+# Solve
+md.private.solution = 'Stressbalance'
+md.settings.waitonlock = 0
+md.toolkits.ToolkitsFile(md.miscellaneous.name + '.toolkits')
+marshall(md) # create .bin file
+
+# export configuration for loading solution in next step
+export_discover(md,'./Models/Antarctica_inversion.nc',delete_rundir=True)
diff --git a/examples/GEOSInput/Antarctica_input_transient.py b/examples/GEOSInput/Antarctica_input_transient.py
new file mode 100755
index 000000000..537058211
--- /dev/null
+++ b/examples/GEOSInput/Antarctica_input_transient.py
@@ -0,0 +1,29 @@
+import devpath
+import os,sys
+ISSM_DIR = os.getenv('ISSM_DIR') # for binaries
+
+from model import *
+from loadmodel import loadmodel
+from clusters.discover_geos import export_discover
+from loadresultsfromdisk import loadresultsfromdisk
+from marshall import marshall
+from verbose import verbose
+
+md = loadmodel('./Models/Antarctica_inversion.nc')
+md = loadresultsfromdisk(md,'AntarcticaGEOS.outbin')
+md.friction.coefficient = md.results.StressbalanceSolution.FrictionCoefficient
+
+
+# Write the binary input file
+# Additional options
+md.inversion.iscontrol = 0
+md.transient.requested_outputs = ['default']
+md.transient.isthermal=0
+md.settings.waitonlock = 0
+md.private.solution = 'Transient'
+md.verbose = verbose('000000000')
+md.toolkits = toolkits()
+marshall(md) # create .bin file
+md.toolkits.ToolkitsFile(md.miscellaneous.name + '.toolkits')
+export_discover(md,'./Models/Antarctica_initialization.nc',delete_rundir=True)
+
diff --git a/examples/GEOSInput/Antarctica_parameterize.py b/examples/GEOSInput/Antarctica_parameterize.py
new file mode 100755
index 000000000..20077a0fa
--- /dev/null
+++ b/examples/GEOSInput/Antarctica_parameterize.py
@@ -0,0 +1,145 @@
+import numpy as np
+from paterson import paterson
+from netCDF4 import Dataset
+from xy2ll import xy2ll
+from InterpFromGridToMesh import InterpFromGridToMesh
+from SetMarineIceSheetBC import SetMarineIceSheetBC
+from m1qn3inversion import m1qn3inversion
+from setflowequation import setflowequation
+import os
+ISSM_DIR = os.getenv('ISSM_DIR')
+
+
+#Name and Coordinate system
+md.miscellaneous.name='AntarcticaGEOS'
+md.mesh.epsg=3031
+
+nsidc_vel = Dataset('../Data/Antarctica_ice_velocity.nc')
+xmin = nsidc_vel.xmin
+xmin = float(xmin.lstrip()[0:10])
+ymax = nsidc_vel.ymax
+ymax = float(ymax.lstrip()[0:10])
+spacing = nsidc_vel.spacing
+spacing = float((spacing.lstrip())[0:4])
+nx = nsidc_vel.nx
+ny = nsidc_vel.ny
+velx = nsidc_vel['vx'][:].data
+vely = nsidc_vel['vy'][:].data
+# Build coordinates
+x2 = xmin + np.arange(nx + 1) * spacing
+y2 = (ymax - ny * spacing) + np.arange(ny + 1) * spacing
+
+# print(' Set observed velocities')
+md.initialization.vx = InterpFromGridToMesh(x2, y2, np.flipud(velx), md.mesh.x, md.mesh.y, 0)
+md.initialization.vy = InterpFromGridToMesh(x2, y2, np.flipud(vely), md.mesh.x, md.mesh.y, 0)
+md.initialization.vz = np.zeros(md.mesh.numberofvertices)
+md.initialization.vel = np.sqrt(md.initialization.vx**2 + md.initialization.vy**2)
+del velx, vely
+
+# Parameters to change/Try
+friction_coefficient = 10 # default [10]
+Temp_change = 0 # default [0 K]
+
+# NetCDF Loading
+print(' Loading SeaRISE data from NetCDF')
+ncdata = Dataset('../Data/Antarctica_5km_withshelves_v0.75.nc')
+x1 = ncdata['x1'][:].data
+y1 = ncdata['y1'][:].data
+usrf = ncdata['usrf'][:].data[0]
+topg = ncdata['topg'][:].data[0]
+temp = ncdata['presartm'][:].data[0]
+smb = ncdata['presprcp'][:].data[0]
+gflux = ncdata['bheatflx_fox'][:].data[0]
+
+# Geometry
+print(' Interpolating surface and ice base')
+md.geometry.base = InterpFromGridToMesh(x1, y1, topg, md.mesh.x, md.mesh.y, 0)
+md.geometry.surface = InterpFromGridToMesh(x1, y1, usrf, md.mesh.x, md.mesh.y, 0)
+del usrf, topg
+
+thkmask=ncdata['thkmask'][:].data[0]
+
+##interpolate onto our mesh vertices
+groundedice= InterpFromGridToMesh(x1,y1,thkmask,md.mesh.x,md.mesh.y,0)
+groundedice[groundedice<=0]=-1
+del thkmask
+
+#fill in the md.mask structure
+md.mask.ocean_levelset = groundedice #ice is grounded for mask equal one
+md.mask.ice_levelset = -1*np.ones(np.shape(md.mesh.x)) #ice is present when negatvie
+
+print(' Constructing thickness')
+md.geometry.thickness = md.geometry.surface - md.geometry.base
+
+# Ensure hydrostatic equilibrium on ice shelf
+di = md.materials.rho_ice / md.materials.rho_water
+
+# Get the node numbers of floating nodes
+pos = np.where(md.mask.ocean_levelset < 0)
+
+# Apply flotation criterion
+md.geometry.thickness[pos] = 1 / (1 - di) * md.geometry.surface[pos]
+md.geometry.base[pos] = md.geometry.surface[pos] - md.geometry.thickness[pos]
+md.geometry.hydrostatic_ratio = np.ones(md.mesh.numberofvertices)
+
+# Set min thickness to 1 meter
+pos0 = np.where(md.geometry.thickness <= 1)
+md.geometry.thickness[pos0] = 1
+md.geometry.surface = md.geometry.thickness + md.geometry.base
+md.geometry.bed = md.geometry.base.copy()
+md.geometry.bed[pos] = md.geometry.base[pos] - 1000
+
+
+# Initialization parameters
+print(' Interpolating temperatures')
+md.initialization.temperature = InterpFromGridToMesh(
+ x1, y1, temp, md.mesh.x, md.mesh.y, 0
+) + 273.15 + Temp_change
+
+[md.mesh.lat, md.mesh.long] = xy2ll(md.mesh.x, md.mesh.y, -1)
+
+
+print(' Set Pressure')
+md.initialization.pressure = md.materials.rho_ice * md.constants.g * md.geometry.thickness
+
+print(' Construct ice rheological properties')
+md.materials.rheology_n = 3 * np.ones(md.mesh.numberofelements)
+md.materials.rheology_B = paterson(md.initialization.temperature)
+
+# Forcings
+print(' Interpolating surface mass balance')
+mass_balance = InterpFromGridToMesh(x1, y1, smb, md.mesh.x, md.mesh.y, 0)
+md.smb.mass_balance = mass_balance * md.materials.rho_water / md.materials.rho_ice
+
+print(' Set geothermal heat flux')
+md.basalforcings.geothermalflux = InterpFromGridToMesh(x1, y1, gflux, md.mesh.x, md.mesh.y, 0)
+
+# Friction and inversion set up
+print(' Construct basal friction parameters')
+md.friction.coefficient = friction_coefficient * np.ones(md.mesh.numberofvertices)
+md.friction.p = np.ones(md.mesh.numberofelements)
+md.friction.q = np.ones(md.mesh.numberofelements)
+
+# No friction applied on floating ice
+pos = np.where(md.mask.ocean_levelset < 0)[0]
+md.friction.coefficient[pos] = 0
+md.groundingline.migration = 'SubelementMigration'
+
+md.inversion = m1qn3inversion()
+md.inversion.vx_obs = md.initialization.vx
+md.inversion.vy_obs = md.initialization.vy
+md.inversion.vel_obs = md.initialization.vel
+
+print(' Set flow equations')
+md = setflowequation(md,'SSA','all')
+
+print(' Set boundary conditions')
+md = SetMarineIceSheetBC(md)
+md.basalforcings.floatingice_melting_rate = np.zeros(md.mesh.numberofvertices)
+md.basalforcings.groundedice_melting_rate = np.zeros(md.mesh.numberofvertices)
+md.thermal.spctemperature = md.initialization.temperature
+md.masstransport.spcthickness = np.full(md.mesh.numberofvertices, np.nan)
+
+
+
+
diff --git a/examples/GEOSInput/Greenland_input_control.py b/examples/GEOSInput/Greenland_input_control.py
index 0d373b756..38ed2f27f 100755
--- a/examples/GEOSInput/Greenland_input_control.py
+++ b/examples/GEOSInput/Greenland_input_control.py
@@ -9,12 +9,10 @@
from InterpFromGridToMesh import InterpFromGridToMesh
from bamg import bamg
from xy2ll import xy2ll
-from plotmodel import plotmodel
from loadmodel import loadmodel
from setmask import setmask
from parameterize import parameterize
from setflowequation import setflowequation
-from solve import solve
from ll2xy import ll2xy
from clusters.discover_geos import export_discover
from marshall import marshall
diff --git a/examples/GEOSInput/Greenland_input_transient.py b/examples/GEOSInput/Greenland_input_transient.py
index 82a7614c7..45a2407f2 100755
--- a/examples/GEOSInput/Greenland_input_transient.py
+++ b/examples/GEOSInput/Greenland_input_transient.py
@@ -2,7 +2,6 @@
import os,sys
ISSM_DIR = os.getenv('ISSM_DIR') # for binaries
-import numpy as np
from model import *
from loadmodel import loadmodel
from clusters.discover_geos import export_discover
@@ -18,7 +17,8 @@
# Write the binary input file
# Additional options
md.inversion.iscontrol = 0
-md.transient.requested_outputs = ['IceVolume', 'TotalSmb', 'SmbMassBalance']
+md.transient.requested_outputs = ['default']
+md.transient.isthermal=0
md.settings.waitonlock = 0
md.private.solution = 'Transient'
md.verbose = verbose('000000000')
diff --git a/examples/GEOSInput/generate_antarctica_input.sh b/examples/GEOSInput/generate_antarctica_input.sh
new file mode 100755
index 000000000..e114167a4
--- /dev/null
+++ b/examples/GEOSInput/generate_antarctica_input.sh
@@ -0,0 +1,5 @@
+source issm_env_withpy
+rm -f AntarcticaGEOS.bin AntarcticaGEOS.outbin AntarcticaGEOS.outlog AntarcticaGEOS.errlog
+(export LD_LIBRARY_PATH=$PYTHON_LIB:$LD_LIBRARY_PATH; python Antarctica_input_control.py)
+${ISSM_DIR}/bin/issm.exe StressbalanceSolution $(pwd) AntarcticaGEOS 2>> AntarcticaGEOS.errlog
+(export LD_LIBRARY_PATH=$PYTHON_LIB:$LD_LIBRARY_PATH; python Antarctica_input_transient.py)
diff --git a/examples/GEOSInput/generate_geos_input.sh b/examples/GEOSInput/generate_greenland_input.sh
similarity index 96%
rename from examples/GEOSInput/generate_geos_input.sh
rename to examples/GEOSInput/generate_greenland_input.sh
index c4fa657ab..be6acca27 100755
--- a/examples/GEOSInput/generate_geos_input.sh
+++ b/examples/GEOSInput/generate_greenland_input.sh
@@ -1,5 +1,4 @@
source issm_env_withpy
-rm -rf Models
rm -f GreenlandGEOS.bin GreenlandGEOS.outbin GreenlandGEOS.outlog GreenlandGEOS.errlog
(export LD_LIBRARY_PATH=$PYTHON_LIB:$LD_LIBRARY_PATH; python Greenland_input_control.py)
${ISSM_DIR}/bin/issm.exe StressbalanceSolution $(pwd) GreenlandGEOS 2>> GreenlandGEOS.errlog
diff --git a/examples/GEOSInput/main_esmf.f90 b/examples/GEOSInput/main_esmf.f90
index 126054c4a..ea30885b4 100755
--- a/examples/GEOSInput/main_esmf.f90
+++ b/examples/GEOSInput/main_esmf.f90
@@ -5,40 +5,40 @@
program main
use iso_fortran_env, only: dp=>real64
- use iso_c_binding, only: c_ptr, c_double, c_f_pointer,c_null_char, c_loc, c_int
+ use iso_c_binding, only: c_ptr, c_double, c_f_pointer,c_null_char, c_loc, c_int, c_char
use ESMF
use netcdf
implicit none
! Define the interface for the ISSM C++ functions
interface
- subroutine InitializeISSM(argc, argv, num_elements, num_nodes, comm) bind(c, name="InitializeISSM")
- import :: c_ptr, c_int
- integer(c_int), value :: argc
- type(c_ptr), dimension(argc) :: argv
- integer(c_int) :: num_elements
- integer(c_int) :: num_nodes
- integer(c_int) :: comm
+ subroutine InitializeISSM(expdir, num_elements, num_nodes, comm) bind(c, name="InitializeISSM")
+ import :: c_char, c_int
+ character(c_char), dimension(*) :: expdir
+ integer(c_int) :: num_elements
+ integer(c_int) :: num_nodes
+ integer(c_int) :: comm
end subroutine InitializeISSM
subroutine RunISSM(dt, gcm_forcings, issm_outputs) bind(C,NAME="RunISSM")
- import :: c_ptr, c_double
+ import :: c_ptr, c_double, c_int
real(c_double), value :: dt
type(c_ptr), value :: gcm_forcings
type(c_ptr), value :: issm_outputs
end subroutine RunISSM
subroutine GetNodesISSM(nodeIds, nodeCoords) bind(C,NAME="GetNodesISSM")
- import :: c_ptr
+ import :: c_ptr, c_int
type(c_ptr), value :: nodeIds
type(c_ptr), value :: nodeCoords
end subroutine GetNodesISSM
- subroutine GetElementsISSM(elementIds, elementConn, elementCoords) bind(C,NAME="GetElementsISSM")
- import :: c_ptr
+ subroutine GetElementsISSM(elementIds, elementConn, elementCoords, glacIds) bind(C,NAME="GetElementsISSM")
+ import :: c_ptr, c_int
type(c_ptr), value :: elementIds
type(c_ptr), value :: elementConn
type(c_ptr), value :: elementCoords
+ type(c_ptr), value :: glacIds
end subroutine GetElementsISSM
subroutine FinalizeISSM() bind(C,NAME="FinalizeISSM")
@@ -87,6 +87,7 @@ end subroutine FinalizeISSM
integer :: rc
type(ESMF_Mesh) :: mesh
integer :: sdim
+ integer, pointer, dimension(:) :: glacIds => null()
integer, pointer, dimension(:) :: elementIds => null()
integer, pointer, dimension(:) :: elementConn => null()
real(dp),pointer, dimension(:) :: elementCoords => null()
@@ -98,41 +99,30 @@ end subroutine FinalizeISSM
integer :: num_outputs = 3
- dt = 0.05 ! timestep in years
-
- ! Get the environment variable ISSM_DIR
- call get_environment_variable("ISSM_DIR", issm_path, length, status)
-
- ! Manually set argc and argv
- argc = 4 ! Example: 3 arguments
- allocate(argv(argc))
- argv(1) = "this arg does not matter"//c_null_char
- argv(2) = "TransientSolution"//c_null_char
- argv(3) = trim(issm_path)//"/examples/GEOSInput"//c_null_char
- argv(4) = "GreenlandGEOS"//c_null_char
+ ! stuff for file counting
+ character(len=256) :: EXPDIR
- ! Convert Fortran strings to C pointers (argv)
- allocate(argv_ptr(argc))
- do i = 1, argc
- ! Ensure that we are only getting the memory address once per string
- argv_ptr(i) = c_loc(argv(i))
- end do
+ integer :: total_elements ! elements for all models (on each PET)
+ integer :: total_nodes ! elements for all models (on each PET)
+
! initialize ESMF to and get vm info / comm for ISSM MPI
call ESMF_Initialize(vm=vm, defaultlogfilename="VMDefaultBasicsEx.Log", logkindflag=ESMF_LOGKIND_MULTI, rc=rc)
call ESMF_VMGet(vm,mpiCommunicator=comm,localPET=localPET,petCount=petCount,rc=rc)
- ! Call the C++ function for initializing ISSM
- ! gets the number of elements and nodes of the mesh
- call ESMF_VMBarrier(vm, rc=rc)
- call InitializeISSM(argc, argv_ptr,num_elements,num_nodes,comm)
- call ESMF_VMBarrier(vm, rc=rc)
- ! ! print out some info if desired:
- !print *, "number of elements on PET ", localPET, ": ", num_elements
- !print *, "number of nodes on PET ", localPET, ": ", num_nodes
+ dt = 0.05 ! timestep in years
+
+ ! Get the environment variable ISSM_DIR
+ call get_environment_variable("ISSM_DIR", issm_path, length, status)
+
+ EXPDIR = trim(issm_path)//"/examples/GEOSInput"//c_null_char
+
+ call InitializeISSM(EXPDIR, num_elements, num_nodes, comm)
+
+ print *, "rank: ", localPET, " num_elements: ", num_elements
! allocate mesh-related pointers
allocate(nodeCoords(2*num_nodes))
@@ -140,6 +130,7 @@ end subroutine FinalizeISSM
allocate(nodeOwners(num_nodes))
allocate(elementTypes(num_elements))
allocate(elementIds(num_elements))
+ allocate(glacIds(num_elements))
allocate(elementConn(3*num_elements))
allocate(elementCoords(2*num_elements))
@@ -156,7 +147,10 @@ end subroutine FinalizeISSM
! create ESMF mesh corresponding to ISSM mesh
! get information about nodes and elements
call GetNodesISSM(c_loc(nodeIds), c_loc(nodeCoords))
- call GetElementsISSM(c_loc(elementIds), c_loc(elementConn), c_loc(elementCoords))
+
+ call GetElementsISSM(c_loc(elementIds), c_loc(elementConn), c_loc(elementCoords),c_loc(glacIds))
+
+ print *, "max elementIds: ", maxval(elementIds)
elementTypes(:) = ESMF_MESHELEMTYPE_TRI
call ESMF_VMBarrier(vm, rc=rc)
@@ -175,6 +169,7 @@ end subroutine FinalizeISSM
end if
end if
+
call ESMF_VMBarrier(vm, rc=rc)
call ESMF_MeshGet(mesh,nodeOwners=nodeOwners)
@@ -195,12 +190,10 @@ end subroutine FinalizeISSM
!call the C++ routine for running a single time step
call RunISSM(dt, c_loc(SMBToISSM), c_loc(issm_outputs))
-
SurfaceOnElements(:) = issm_outputs(1:num_elements)
ICEEL(:) = issm_outputs(1:num_elements)
ICETHICK(:) = issm_outputs(num_elements+1:2*num_elements)
ICEVEL(:) = issm_outputs(2*num_elements+1:3*num_elements)
-
call ESMF_VMBarrier(vm, rc=rc)
srcField = ESMF_FieldCreate(mesh=mesh,farrayPtr=SurfaceOnElements,meshloc=ESMF_MESHLOC_ELEMENT, &
diff --git a/examples/GEOSInput/plotting.ipynb b/examples/GEOSInput/plotting.ipynb
index 491e0428d..520e40ab4 100644
--- a/examples/GEOSInput/plotting.ipynb
+++ b/examples/GEOSInput/plotting.ipynb
@@ -2,10 +2,19 @@
"cells": [
{
"cell_type": "code",
- "execution_count": 9,
+ "execution_count": 1,
"id": "2e7146ec-37eb-46bc-88c2-d91ef3afb2f5",
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/usr/local/other/anaconda/GEOSpyD/24.3.0-0/2024-08-29/envs/py3.12/envs/viz/lib/python3.11/site-packages/pydap/lib.py:5: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.\n",
+ " from pkg_resources import get_distribution\n"
+ ]
+ }
+ ],
"source": [
"import netCDF4\n",
"import numpy as np\n",
@@ -13,10 +22,12 @@
"from matplotlib.tri import Triangulation\n",
"import cartopy.crs as ccrs\n",
"import cartopy.feature as cfeature\n",
+ "from pyproj import Transformer\n",
"import xarray as xr\n",
"import glob\n",
"from matplotlib import ticker\n",
- "import matplotlib.colors as colors"
+ "import matplotlib.colors as colors\n",
+ "from matplotlib import cm"
]
},
{
@@ -392,36 +403,34 @@
" fill: currentColor;\n",
"}\n",
"<xarray.Dataset>\n",
- "Dimensions: (num_nodes: 4064, num_elements: 7028, num_owned_nodes: 4064)\n",
+ "Dimensions: (num_nodes: 1210, num_elements: 2100, num_owned_nodes: 1210)\n",
"Dimensions without coordinates: num_nodes, num_elements, num_owned_nodes\n",
"Data variables: (12/14)\n",
- " node_lon (num_nodes) float32 -37.96 -38.13 -38.44 ... -22.34 -66.87\n",
- " node_lat (num_nodes) float32 66.21 66.06 66.01 ... 65.15 75.1 76.49\n",
- " node_ids (num_nodes) float32 1.0 2.0 3.0 ... 4.063e+03 4.064e+03\n",
+ " node_lon (num_nodes) float32 37.89 37.06 36.29 ... -41.21 -25.48\n",
+ " node_lat (num_nodes) float32 -69.75 -69.69 -69.28 ... 75.94 74.54\n",
+ " node_ids (num_nodes) float32 291.0 292.0 293.0 ... 7.211e+03 7.531e+03\n",
" node_owners (num_nodes) float32 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0\n",
- " element_lon (num_elements) float32 -37.64 -37.68 -37.64 ... -37.57 -37.61\n",
- " element_lat (num_elements) float32 66.57 66.54 66.55 ... 66.47 66.5 66.5\n",
+ " element_lon (num_elements) float32 18.72 18.24 20.48 ... -39.0 -40.11\n",
+ " element_lat (num_elements) float32 -71.9 -72.16 -72.1 ... 75.16 75.7\n",
" ... ...\n",
- " element_conn2 (num_elements) float32 711.0 3.981e+03 ... 712.0 3.979e+03\n",
- " element_conn3 (num_elements) float32 703.0 705.0 ... 3.979e+03 980.0\n",
- " node_surf (num_owned_nodes) float32 0.0 0.0 0.0 ... 1.215e+03 625.3\n",
- " element_surf (num_elements) float32 756.2 867.7 691.3 ... 563.7 569.3\n",
- " ice_thick (num_elements) float32 808.4 503.7 944.4 ... 415.4 794.0\n",
- " ice_vel (num_elements) float32 2.906e-05 1.647e-05 ... 2.979e-05- num_nodes: 4064
- num_elements: 7028
- num_owned_nodes: 4064
node_lon
(num_nodes)
float32
...
array([-37.95821 , -38.13414 , -38.435802, ..., -40.172302, -22.342548,\n",
- " -66.86685 ], dtype=float32)
node_lat
(num_nodes)
float32
...
array([66.20738 , 66.05826 , 66.00676 , ..., 65.15363 , 75.098915, 76.49378 ],\n",
- " dtype=float32)
node_ids
(num_nodes)
float32
...
array([1.000e+00, 2.000e+00, 3.000e+00, ..., 4.062e+03, 4.063e+03, 4.064e+03],\n",
- " dtype=float32)
node_owners
(num_nodes)
float32
...
array([0., 0., 0., ..., 0., 0., 0.], dtype=float32)
element_lon
(num_elements)
float32
...
array([-37.64498 , -37.68052 , -37.63686 , ..., -37.56055 , -37.57282 ,\n",
- " -37.609463], dtype=float32)
element_lat
(num_elements)
float32
...
array([66.57463 , 66.54456 , 66.549576, ..., 66.467606, 66.49856 , 66.49967 ],\n",
- " dtype=float32)
element_ids
(num_elements)
float32
...
array([1.000e+00, 2.000e+00, 3.000e+00, ..., 7.026e+03, 7.027e+03, 7.028e+03],\n",
- " dtype=float32)
element_conn1
(num_elements)
float32
...
array([4040., 980., 3980., ..., 704., 713., 979.], dtype=float32)
element_conn2
(num_elements)
float32
...
array([ 711., 3981., 711., ..., 3861., 712., 3979.], dtype=float32)
element_conn3
(num_elements)
float32
...
array([ 703., 705., 3981., ..., 979., 3979., 980.], dtype=float32)
node_surf
(num_owned_nodes)
float32
...
array([ 0. , 0. , 0. , ..., 389.0679, 1215.3248, 625.2626],\n",
- " dtype=float32)
element_surf
(num_elements)
float32
...
array([756.23944, 867.7225 , 691.3218 , ..., 682.693 , 563.65063, 569.26044],\n",
- " dtype=float32)
ice_thick
(num_elements)
float32
...
array([808.4347 , 503.74957, 944.4406 , ..., 530.3279 , 415.359 , 793.9868 ],\n",
- " dtype=float32)
ice_vel
(num_elements)
float32
...
array([2.906376e-05, 1.647293e-05, 2.837983e-05, ..., 1.840896e-05,\n",
- " 1.622319e-05, 2.978819e-05], dtype=float32)
- num_nodes: 1210
- num_elements: 2100
- num_owned_nodes: 1210
node_lon
(num_nodes)
float32
...
array([ 37.89375 , 37.060917, 36.286823, ..., -25.981052, -41.20778 ,\n",
+ " -25.480812], dtype=float32)
node_lat
(num_nodes)
float32
...
array([-69.74513 , -69.69277 , -69.27664 , ..., 74.688675, 75.93859 ,\n",
+ " 74.53882 ], dtype=float32)
node_ids
(num_nodes)
float32
...
array([ 291., 292., 293., ..., 6900., 7211., 7531.], dtype=float32)
node_owners
(num_nodes)
float32
...
array([0., 0., 0., ..., 0., 0., 0.], dtype=float32)
element_lon
(num_elements)
float32
...
array([ 18.72468 , 18.242746, 20.483683, ..., -30.09065 , -39.00315 ,\n",
+ " -40.11061 ], dtype=float32)
element_lat
(num_elements)
float32
...
array([-71.89896 , -72.156654, -72.10306 , ..., 77.663994, 75.16285 ,\n",
+ " 75.698784], dtype=float32)
element_ids
(num_elements)
float32
...
array([ 1341., 1342., 1343., ..., 13198., 13466., 13467.], dtype=float32)
element_conn1
(num_elements)
float32
...
array([ 425., 251., 357., ..., 908., 1207., 1209.], dtype=float32)
element_conn2
(num_elements)
float32
...
array([ 82., 425., 148., ..., 821., 1197., 1207.], dtype=float32)
element_conn3
(num_elements)
float32
...
array([ 391., 391., 425., ..., 893., 854., 1200.], dtype=float32)
node_surf
(num_owned_nodes)
float32
...
array([ 0. , 0. , 0. , ..., 1839.776 , 2785.441 , 1711.3314],\n",
+ " dtype=float32)
element_surf
(num_elements)
float32
...
array([1627.6075, 1994.5875, 1673.3105, ..., 1822.6097, 2841.8447, 2804.0322],\n",
+ " dtype=float32)
ice_thick
(num_elements)
float32
...
array([ 442.83444, 537.81976, 699.25867, ..., 1650.8022 , 2962.7068 ,\n",
+ " 2919.2612 ], dtype=float32)
ice_vel
(num_elements)
float32
...
array([1.116207e-06, 1.045563e-06, 1.409277e-06, ..., 2.525359e-06,\n",
+ " 3.088083e-07, 2.584707e-07], dtype=float32)