| 1 | import numpy as np
|
|---|
| 2 | from triangle import triangle
|
|---|
| 3 | from model import *
|
|---|
| 4 | from netCDF4 import Dataset
|
|---|
| 5 | from InterpFromGridToMesh import InterpFromGridToMesh
|
|---|
| 6 | from bamg import bamg
|
|---|
| 7 | from xy2ll import xy2ll
|
|---|
| 8 | from plotmodel import plotmodel
|
|---|
| 9 | from export_netCDF import export_netCDF
|
|---|
| 10 | from loadmodel import loadmodel
|
|---|
| 11 | from setmask import setmask
|
|---|
| 12 | from parameterize import parameterize
|
|---|
| 13 | from setflowequation import setflowequation
|
|---|
| 14 | from socket import gethostname
|
|---|
| 15 | from solve import solve
|
|---|
| 16 | from ll2xy import ll2xy
|
|---|
| 17 | from BamgTriangulate import BamgTriangulate
|
|---|
| 18 | from InterpFromMeshToMesh2d import InterpFromMeshToMesh2d
|
|---|
| 19 | from scipy.interpolate import griddata
|
|---|
| 20 |
|
|---|
| 21 | steps = [1]
|
|---|
| 22 |
|
|---|
| 23 | if 1 in steps:
|
|---|
| 24 | # Step 1: Mesh creation {{{
|
|---|
| 25 | print(' Step 1: Mesh creation')
|
|---|
| 26 |
|
|---|
| 27 | #Generate initial uniform mesh (resolution = 20000 m)
|
|---|
| 28 | md = triangle(model(), './DomainOutline.exp', 20000)
|
|---|
| 29 |
|
|---|
| 30 | ncdata = Dataset('../Data/Greenland_5km_dev1.2.nc', mode='r')
|
|---|
| 31 |
|
|---|
| 32 | # Get velocities (Note: You can use ncprint('file') to see an ncdump)
|
|---|
| 33 | x1 = np.squeeze(ncdata.variables['x1'][:].data)
|
|---|
| 34 | y1 = np.squeeze(ncdata.variables['y1'][:].data)
|
|---|
| 35 | velx = np.squeeze(ncdata.variables['surfvelx'][:].data)
|
|---|
| 36 | vely = np.squeeze(ncdata.variables['surfvely'][:].data)
|
|---|
| 37 | ncdata.close()
|
|---|
| 38 |
|
|---|
| 39 | vx = InterpFromGridToMesh(x1, y1, velx, md.mesh.x, md.mesh.y, 0)
|
|---|
| 40 | vy = InterpFromGridToMesh(x1, y1, vely, md.mesh.x, md.mesh.y, 0)
|
|---|
| 41 | vel = np.sqrt(vx**2 + vy**2)
|
|---|
| 42 |
|
|---|
| 43 | #Mesh Greenland
|
|---|
| 44 | md = bamg(md, 'hmax', 400000, 'hmin', 5000, 'gradation', 1.4, 'field', vel, 'err', 8)
|
|---|
| 45 |
|
|---|
| 46 | #convert x, y coordinates (Polar stereo) to lat / lon
|
|---|
| 47 | [md.mesh.lat, md.mesh.long] = xy2ll(md.mesh.x, md.mesh.y, + 1, 39, 71)
|
|---|
| 48 |
|
|---|
| 49 | export_netCDF(md, './Models/Greenland.Mesh_generation.nc')
|
|---|
| 50 | plotmodel(md, 'data', 'mesh')
|
|---|
| 51 | # }}}
|
|---|
| 52 |
|
|---|
| 53 | if 2 in steps:
|
|---|
| 54 | # Step 2: Parameterization{{{
|
|---|
| 55 | print(' Step 2: Parameterization')
|
|---|
| 56 | md = loadmodel('./Models/Greenland.Mesh_generation.nc')
|
|---|
| 57 |
|
|---|
| 58 | md = setmask(md, '', '')
|
|---|
| 59 | md = parameterize(md, './Greenland.py')
|
|---|
| 60 | md = setflowequation(md, 'SSA', 'all')
|
|---|
| 61 |
|
|---|
| 62 | export_netCDF(md, "./Models/Greenland.Parameterization.nc")
|
|---|
| 63 | # }}}
|
|---|
| 64 |
|
|---|
| 65 | if 3 in steps:
|
|---|
| 66 | # Step 3: Control method friction {{{
|
|---|
| 67 | print(' Step 3: Control method friction')
|
|---|
| 68 | md = loadmodel('./Models/Greenland.Parameterization.nc')
|
|---|
| 69 | #Control general
|
|---|
| 70 | md.inversion.iscontrol = 1
|
|---|
| 71 | md.inversion.nsteps = 30
|
|---|
| 72 | md.inversion.step_threshold = 0.99 * np.ones((md.inversion.nsteps))
|
|---|
| 73 | md.inversion.maxiter_per_step = 5 * np.ones((md.inversion.nsteps))
|
|---|
| 74 |
|
|---|
| 75 | #Cost functions
|
|---|
| 76 | md.inversion.cost_functions = [101, 103, 501]
|
|---|
| 77 | md.inversion.cost_functions_coefficients = np.ones((md.mesh.numberofvertices, 3))
|
|---|
| 78 | md.inversion.cost_functions_coefficients[:, 0] = 350
|
|---|
| 79 | md.inversion.cost_functions_coefficients[:, 1] = 0.2
|
|---|
| 80 | md.inversion.cost_functions_coefficients[:, 2] = 2e-6
|
|---|
| 81 |
|
|---|
| 82 | #Controls
|
|---|
| 83 | md.inversion.control_parameters = ['FrictionCoefficient']
|
|---|
| 84 | md.inversion.gradient_scaling = 50 * np.ones((md.inversion.nsteps, 1))
|
|---|
| 85 | md.inversion.min_parameters = 1 * np.ones((md.mesh.numberofvertices, 1))
|
|---|
| 86 | md.inversion.max_parameters = 200 * np.ones((md.mesh.numberofvertices, 1))
|
|---|
| 87 |
|
|---|
| 88 | #Additional parameters
|
|---|
| 89 | md.stressbalance.restol = 0.01
|
|---|
| 90 | md.stressbalance.reltol = 0.1
|
|---|
| 91 | md.stressbalance.abstol = np.nan
|
|---|
| 92 | md.stressbalance.loadingforce = np.zeros((md.mesh.numberofvertices, 3))
|
|---|
| 93 |
|
|---|
| 94 | #Go solve
|
|---|
| 95 | md.cluster = generic('name', gethostname(), 'np', 2)
|
|---|
| 96 | md.toolkits = toolkits()
|
|---|
| 97 | md.verbose = verbose('solution', True, 'control', True)
|
|---|
| 98 | md = solve(md, 'Stressbalance')
|
|---|
| 99 |
|
|---|
| 100 | #Update model friction fields accordingly
|
|---|
| 101 | md.friction.coefficient = md.results.StressbalanceSolution.FrictionCoefficient
|
|---|
| 102 |
|
|---|
| 103 | export_netCDF(md, "./Models/Greenland.Control_drag.nc")
|
|---|
| 104 | # }}}
|
|---|
| 105 |
|
|---|
| 106 | if 4 in steps:
|
|---|
| 107 | # Step 4: Transient run {{{
|
|---|
| 108 | print(' Step 4: Transient run')
|
|---|
| 109 | md = loadmodel('./Models/Greenland.Control_drag.nc')
|
|---|
| 110 |
|
|---|
| 111 | #Set surface mass balance
|
|---|
| 112 | ncdata = Dataset('../Data/Greenland_5km_dev1.2.nc', mode='r')
|
|---|
| 113 | x1 = np.squeeze(ncdata.variables['x1'][:].data)
|
|---|
| 114 | y1 = np.squeeze(ncdata.variables['y1'][:].data)
|
|---|
| 115 | smb = np.squeeze(ncdata.variables['smb'][:].data)
|
|---|
| 116 | ncdata.close()
|
|---|
| 117 |
|
|---|
| 118 | smb = InterpFromGridToMesh(x1, y1, smb, md.mesh.x, md.mesh.y, 0)
|
|---|
| 119 | smb = smb * md.materials.rho_freshwater / md.materials.rho_ice
|
|---|
| 120 | smb = np.vstack((smb, smb, smb - 1.0)).T
|
|---|
| 121 | md.smb.mass_balance = np.vstack((smb, np.asarray([[1, 10, 20]])))
|
|---|
| 122 | #Set transient options, run for 20 years, saving every timestep
|
|---|
| 123 | md.timestepping.time_step = 0.2
|
|---|
| 124 | md.timestepping.final_time = 20
|
|---|
| 125 | md.settings.output_frequency = 1
|
|---|
| 126 |
|
|---|
| 127 | #Additional options
|
|---|
| 128 | md.inversion.iscontrol = 0
|
|---|
| 129 | md.transient.requested_outputs = ['IceVolume', 'TotalSmb', 'SmbMassBalance']
|
|---|
| 130 | md.verbose = verbose('solution', True, 'module', True, 'convergence', True)
|
|---|
| 131 |
|
|---|
| 132 | #Go solve
|
|---|
| 133 | md.cluster = generic('name', gethostname(), 'np', 2)
|
|---|
| 134 | md = solve(md, 'Transient')
|
|---|
| 135 |
|
|---|
| 136 | export_netCDF(md, './Models/Greenland.Transient.nc')
|
|---|
| 137 | # }}}
|
|---|
| 138 |
|
|---|
| 139 | if 5 in steps:
|
|---|
| 140 | # Step 5: Plotting {{{
|
|---|
| 141 | print(' Step 5: Plotting')
|
|---|
| 142 | md = loadmodel('./Models/Greenland.Transient.nc')
|
|---|
| 143 |
|
|---|
| 144 | #Planview plots
|
|---|
| 145 | plotmodel(md, 'data', md.results.TransientSolution[-1].Vel, 'log#1', 1e-1,
|
|---|
| 146 | 'caxis#1', [1e-1, 1e4], 'title#1', 'Velocity (m / y)',
|
|---|
| 147 | 'data', md.results.TransientSolution[1].SmbMassBalance,
|
|---|
| 148 | 'title#2', 'Surface Mass Balance (m / y)',
|
|---|
| 149 | 'data', md.results.TransientSolution[-1].Thickness,
|
|---|
| 150 | 'title', 'Thickness (m)',
|
|---|
| 151 | 'data', md.results.TransientSolution[-1].Surface,
|
|---|
| 152 | 'title', 'Surface (m)')
|
|---|
| 153 |
|
|---|
| 154 | #Line Plots
|
|---|
| 155 | figure
|
|---|
| 156 |
|
|---|
| 157 | #Plot surface mass balance, velocity and volume
|
|---|
| 158 | surfmb = []
|
|---|
| 159 | vel = []
|
|---|
| 160 | volume = []
|
|---|
| 161 | for i in range(0, 100):
|
|---|
| 162 | surfmb.append(md.results.TransientSolution[i].SmbMassBalance)
|
|---|
| 163 | vel.append(md.results.TransientSolution[i].Vel)
|
|---|
| 164 | volume.append(md.results.TransientSolution[i].IceVolume)
|
|---|
| 165 |
|
|---|
| 166 | layout, ax = plt.subplots(3, 1, sharex=True, sharey=False, figsize=(5, 5))
|
|---|
| 167 | ax[0].plot(np.arange(0.2, 20.2, 0.2), np.nanmean(surfmb, axis=1))
|
|---|
| 168 | ax[0].set_title('Mean Surface mass balance')
|
|---|
| 169 | ax[1].plot(np.arange(0.2, 20.2, 0.2), np.nanmean(vel, axis=1))
|
|---|
| 170 | ax[1].set_title('Mean Velocity')
|
|---|
| 171 | ax[2].plot(np.arange(0.2, 20.2, 0.2), volume)
|
|---|
| 172 | ax[2].set_title('Ice Volume')
|
|---|
| 173 | ax[2].set_xlabel('years')
|
|---|
| 174 | # }}}
|
|---|
| 175 |
|
|---|
| 176 | if 6 in steps:
|
|---|
| 177 | # Step 6: Extract Box SMB{{{
|
|---|
| 178 | print(' Step 6: Extract Box SMB')
|
|---|
| 179 | md = loadmodel('./Models/Greenland.Transient.nc')
|
|---|
| 180 |
|
|---|
| 181 | #Set surface mass balance
|
|---|
| 182 | ncbox = Dataset('../Data/Box_Greenland_SMB_monthly_1840-2012_5km_cal_ver20141007.nc', mode='r')
|
|---|
| 183 | lat = np.squeeze(ncbox.variables['lat'][:].data)
|
|---|
| 184 | lon = np.squeeze(ncbox.variables['lon'][:].data)
|
|---|
| 185 | smbbox = np.squeeze(ncbox.variables['MassFlux'][:].data)
|
|---|
| 186 | ncbox.close()
|
|---|
| 187 | [x1, y1] = ll2xy(lat, lon, + 1, 45, 70)
|
|---|
| 188 |
|
|---|
| 189 | years_of_simulation = np.arange(1840, 2012 + 1)
|
|---|
| 190 | t = np.arange(years_of_simulation[0], years_of_simulation[-1] + 11 / 12, 1 / 12)
|
|---|
| 191 | #Area of grid for 5km box
|
|---|
| 192 | area_of_grid = 5000 * 5000
|
|---|
| 193 | totalsmb = reshape(np.nansum(smbbox / 1000, axis=(-2, -1)), (len(t), 1)) * area_of_grid
|
|---|
| 194 |
|
|---|
| 195 | #save surface mass balance mat dataset
|
|---|
| 196 | smbmean = np.nanmean(smbbox, axis=(0, 1))
|
|---|
| 197 | SMB = {}
|
|---|
| 198 | SMB['smbmean'] = smbmean
|
|---|
| 199 | SMB['totalsmb '] = totalsmb
|
|---|
| 200 | SMB['smbbox'] = smbbox
|
|---|
| 201 | SMB['x1'] = x1
|
|---|
| 202 | SMB['y1'] = y1
|
|---|
| 203 | SMB['t'] = t
|
|---|
| 204 |
|
|---|
| 205 | np.savez('./smbbox.npz', **SMB)
|
|---|
| 206 |
|
|---|
| 207 | #plot a time series of total SMB
|
|---|
| 208 | fig = plt.figure(tight_layout=True)
|
|---|
| 209 | ax = fig.add_subplot(111)
|
|---|
| 210 | ax.plot(t, totalsmb / 1e9)
|
|---|
| 211 | ax.set_title('Total Surface mass balance, Gt')
|
|---|
| 212 | ax.set_xlabel('year')
|
|---|
| 213 | ax.set_ylabel('Gt/yr')
|
|---|
| 214 |
|
|---|
| 215 | del smbbox
|
|---|
| 216 | # }}}
|
|---|
| 217 |
|
|---|
| 218 | if 7 in steps:
|
|---|
| 219 | # Step 7: Historical Relaxation run {{{
|
|---|
| 220 | print(' Step 7: Historical Relaxation run')
|
|---|
| 221 | md = loadmodel('./Models/Greenland.Control_drag.nc')
|
|---|
| 222 |
|
|---|
| 223 | with open('./smbbox.npz', "rb") as smbFile:
|
|---|
| 224 | SMB = np.load(smbFile)
|
|---|
| 225 | x1 = np.squeeze(SMB['x1'])
|
|---|
| 226 | y1 = np.squeeze(SMB['y1'])
|
|---|
| 227 | smbmean = np.squeeze(SMB['smbmean'])
|
|---|
| 228 |
|
|---|
| 229 | #convert mesh x, y into the Box projection
|
|---|
| 230 | [md.mesh.lat, md.mesh.long] = xy2ll(md.mesh.x, md.mesh.y, + 1, 39, 71)
|
|---|
| 231 | [xi, yi] = ll2xy(md.mesh.lat, md.mesh.long, + 1, 45, 70)
|
|---|
| 232 |
|
|---|
| 233 | #Interpolate and set surface mass balance
|
|---|
| 234 | x1 = x1.flatten()
|
|---|
| 235 | y1 = y1.flatten()
|
|---|
| 236 | smbmean = smbmean.flatten()
|
|---|
| 237 | index = BamgTriangulate(x1, y1)
|
|---|
| 238 | smb_mo = InterpFromMeshToMesh2d(index, x1, y1, smbmean, xi, yi)
|
|---|
| 239 | smb = smb_mo * 12 / 1000 * md.materials.rho_freshwater / md.materials.rho_ice
|
|---|
| 240 | md.smb.mass_balance = np.append(smb, 1)
|
|---|
| 241 |
|
|---|
| 242 | #Set transient options, run for 20 years, saving every timestep
|
|---|
| 243 | md.timestepping.time_step = 0.2
|
|---|
| 244 | md.timestepping.final_time = 20
|
|---|
| 245 | md.settings.output_frequency = 1
|
|---|
| 246 |
|
|---|
| 247 | #Additional options
|
|---|
| 248 | md.inversion.iscontrol = 0
|
|---|
| 249 | md.transient.requested_outputs = ['IceVolume', 'TotalSmb', 'SmbMassBalance']
|
|---|
| 250 | md.verbose = verbose('solution', True, 'module', True)
|
|---|
| 251 |
|
|---|
| 252 | #Go solve
|
|---|
| 253 | md.cluster = generic('name', gethostname(), 'np', 2)
|
|---|
| 254 | md = solve(md, 'Transient')
|
|---|
| 255 |
|
|---|
| 256 | export_netCDF(md, './Models/Greenland.HistoricTransient.nc')
|
|---|
| 257 | # }}}
|
|---|
| 258 |
|
|---|
| 259 | if 8 in steps:
|
|---|
| 260 | # Step 8: Plotting exercise {{{
|
|---|
| 261 | print(' Step 8: Plotting exercise')
|
|---|
| 262 | #Load historic transient model
|
|---|
| 263 |
|
|---|
| 264 | #Create Line Plots of relaxation run. Create a figure.
|
|---|
| 265 |
|
|---|
| 266 | #Save surface mass balance, by looping through 200 years (1000 steps)
|
|---|
| 267 | #Note, the first output will always contain output from time step 1
|
|---|
| 268 |
|
|---|
| 269 | #Plot surface mass balance time series in first subplot
|
|---|
| 270 |
|
|---|
| 271 | #Title this plot Mean surface mass balance
|
|---|
| 272 |
|
|---|
| 273 | #Save velocity by looping through 200 years
|
|---|
| 274 |
|
|---|
| 275 | #Plot velocity time series in second subplot
|
|---|
| 276 |
|
|---|
| 277 | #Title this plot Mean Velocity
|
|---|
| 278 |
|
|---|
| 279 | #Save Ice Volume by looping through 200 years
|
|---|
| 280 |
|
|---|
| 281 | #Plot volume time series in third subplot
|
|---|
| 282 |
|
|---|
| 283 | #Title this plot Mean Velocity and add an x label of years
|
|---|
| 284 | # }}}
|
|---|
| 285 |
|
|---|
| 286 | if 9 in steps:
|
|---|
| 287 | # Step 9: Box Transient run{{{
|
|---|
| 288 | print(' Step 9: Box Transient run')
|
|---|
| 289 | md = loadmodel('./Models/Greenland.HistoricTransient.nc')
|
|---|
| 290 |
|
|---|
| 291 | #load past transient results
|
|---|
| 292 | md.geometry.base = md.results.TransientSolution[-1].Base
|
|---|
| 293 | md.geometry.thickness = md.results.TransientSolution[-1].Thickness
|
|---|
| 294 | md.geometry.surface = md.geometry.base + md.geometry.thickness
|
|---|
| 295 | md.initialization.vx = md.results.TransientSolution[-1].Vx
|
|---|
| 296 | md.initialization.vy = md.results.TransientSolution[-1].Vy
|
|---|
| 297 | md.results = []
|
|---|
| 298 |
|
|---|
| 299 | #convert mesh x, y into the Box projection
|
|---|
| 300 | [md.mesh.lat, md.mesh.long] = xy2ll(md.mesh.x, md.mesh.y, + 1, 39, 71)
|
|---|
| 301 | [xi, yi] = ll2xy(md.mesh.lat, md.mesh.long, + 1, 45, 70)
|
|---|
| 302 |
|
|---|
| 303 | #Set surface mass balance
|
|---|
| 304 | with open('./smbbox.npz', "rb") as smbFile:
|
|---|
| 305 | SMB = np.load(smbFile)
|
|---|
| 306 | x1 = np.squeeze(SMB['x1'])
|
|---|
| 307 | y1 = np.squeeze(SMB['y1'])
|
|---|
| 308 | smbbox = np.squeeze(SMB['smbbox'])
|
|---|
| 309 |
|
|---|
| 310 | x1 = x1.flatten()
|
|---|
| 311 | y1 = y1.flatten()
|
|---|
| 312 | index = BamgTriangulate(x1, y1)
|
|---|
| 313 | #Set years to run
|
|---|
| 314 | years_of_simulation = np.arange(2003, 2012 + 1)
|
|---|
| 315 |
|
|---|
| 316 | #initialize surface mass balance matrix
|
|---|
| 317 | smb = np.nan * np.ones((md.mesh.numberofvertices, len(years_of_simulation) * 12))
|
|---|
| 318 |
|
|---|
| 319 | #Interpolate and set surface mass balance
|
|---|
| 320 | for year in years_of_simulation:
|
|---|
| 321 | for month in range(0, 12):
|
|---|
| 322 | smb_mo = griddata((np.double(x1), np.double(y1)), np.double(smbbox[year - 1840, month, :, :].flatten()), (xi, yi), method='nearest')
|
|---|
| 323 | smb[:, (year - years_of_simulation[0]) * 12 + month] = smb_mo
|
|---|
| 324 | smb = smb * 12 / 1000 * md.materials.rho_freshwater / md.materials.rho_ice
|
|---|
| 325 | timer = np.arange(1 / 24, len(years_of_simulation), 1 / 12)
|
|---|
| 326 | md.smb.mass_balance = np.vstack((smb, timer))
|
|---|
| 327 |
|
|---|
| 328 | #Set transient options, monthly timestep, saving every month
|
|---|
| 329 | md.timestepping.time_step = 1 / 12
|
|---|
| 330 | md.timestepping.final_time = len(years_of_simulation)
|
|---|
| 331 | md.settings.output_frequency = 1
|
|---|
| 332 |
|
|---|
| 333 | #Additional options
|
|---|
| 334 | md.inversion.iscontrol = 0
|
|---|
| 335 | md.transient.requested_outputs = ['IceVolume', 'TotalSmb', 'SmbMassBalance']
|
|---|
| 336 | md.verbose = verbose('solution', True, 'module', True)
|
|---|
| 337 |
|
|---|
| 338 | #Go solve
|
|---|
| 339 | md.cluster = generic('name', gethostname(), 'np', 2)
|
|---|
| 340 | md = solve(md, 'Transient')
|
|---|
| 341 |
|
|---|
| 342 | export_netCDF(md, './Models/Greenland.BoxTransient.nc')
|
|---|
| 343 | # }}}
|
|---|
| 344 |
|
|---|
| 345 | if 10 in steps:
|
|---|
| 346 | print(' Step 10: Plot Box Transient')
|
|---|
| 347 | md = loadmodel('./Models/Greenland.BoxTransient.nc')
|
|---|
| 348 |
|
|---|
| 349 | #Set years run
|
|---|
| 350 | years_of_simulation = np.arange(2003, 2012 + 1)
|
|---|
| 351 | t = np.arange(years_of_simulation[0], years_of_simulation[-1] + 11 / 12, 1 / 12)
|
|---|
| 352 |
|
|---|
| 353 | #Line Plots
|
|---|
| 354 | layout, ax = plt.subplots(3, 1, sharex=True, sharey=False, figsize=(5, 5))
|
|---|
| 355 | #Plot surface mass balance
|
|---|
| 356 | surfmb = []
|
|---|
| 357 | vel = []
|
|---|
| 358 | volume = []
|
|---|
| 359 | for i in range(0, len(t)):
|
|---|
| 360 | surfmb.append(md.results.TransientSolution[i].TotalSmb)
|
|---|
| 361 | vel.append(md.results.TransientSolution[i].Vel)
|
|---|
| 362 | volume.append(md.results.TransientSolution[i].IceVolume)
|
|---|
| 363 | ax[0].plot(t, surfmb)
|
|---|
| 364 | ax[0].set_title('Total Surface mass balance')
|
|---|
| 365 | ax[1].plot(t, np.nanmax(vel, axis=1))
|
|---|
| 366 | ax[1].set_title('Max Velocity')
|
|---|
| 367 | ax[2].plot(t, volume)
|
|---|
| 368 | ax[2].set_title('Ice Volume')
|
|---|
| 369 | ax[2].set_xlabel('years')
|
|---|
![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
