[26559] | 1 | if 7 in steps:
|
---|
| 2 | # Step 7: Historical Relaxation run {{{
|
---|
| 3 | print(' Step 7: Historical Relaxation run')
|
---|
| 4 | md = loadmodel('./Models/Greenland.Control_drag.nc')
|
---|
| 5 |
|
---|
| 6 | with open('./smbbox.npz', "rb") as smbFile:
|
---|
| 7 | SMB = np.load(smbFile)
|
---|
| 8 | x1 = np.squeeze(SMB['x1'])
|
---|
| 9 | y1 = np.squeeze(SMB['y1'])
|
---|
| 10 | smbmean = np.squeeze(SMB['smbmean'])
|
---|
| 11 |
|
---|
| 12 | #convert mesh x, y into the Box projection
|
---|
| 13 | [md.mesh.lat, md.mesh.long] = xy2ll(md.mesh.x, md.mesh.y, + 1, 39, 71)
|
---|
| 14 | [xi, yi] = ll2xy(md.mesh.lat, md.mesh.long, + 1, 45, 70)
|
---|
| 15 |
|
---|
| 16 | #Interpolate and set surface mass balance
|
---|
| 17 | x1 = x1.flatten()
|
---|
| 18 | y1 = y1.flatten()
|
---|
| 19 | smbmean = smbmean.flatten()
|
---|
| 20 | index = BamgTriangulate(x1, y1)
|
---|
| 21 | smb_mo = InterpFromMeshToMesh2d(index, x1, y1, smbmean, xi, yi)
|
---|
| 22 | smb = smb_mo * 12 / 1000 * md.materials.rho_freshwater / md.materials.rho_ice
|
---|
| 23 | md.smb.mass_balance = np.append(smb, 1)
|
---|
| 24 |
|
---|
| 25 | #Set transient options, run for 20 years, saving every timestep
|
---|
| 26 | md.timestepping.time_step = 0.2
|
---|
| 27 | md.timestepping.final_time = 200
|
---|
| 28 | md.settings.output_frequency = 5
|
---|
| 29 |
|
---|
| 30 | #Additional options
|
---|
| 31 | md.inversion.iscontrol = 0
|
---|
| 32 | md.transient.requested_outputs = ['IceVolume', 'TotalSmb', 'SmbMassBalance']
|
---|
| 33 | md.verbose = verbose('solution', True, 'module', True)
|
---|
| 34 |
|
---|
| 35 | #Go solve
|
---|
| 36 | md.cluster = generic('name', gethostname(), 'np', 2)
|
---|
| 37 | md = solve(md, 'Transient')
|
---|
| 38 |
|
---|
| 39 | export_netCDF(md, './Models/Greenland.HistoricTransient_200yr.nc')
|
---|
| 40 | # }}}
|
---|
| 41 |
|
---|
| 42 | if 8 in steps:
|
---|
| 43 | # Step 8: Plotting exercise {{{
|
---|
| 44 | print(' Step 8: Plotting exercise')
|
---|
| 45 | md = loadmodel('./Models/Greenland.HistoricTransient_200yr.nc')
|
---|
| 46 |
|
---|
| 47 | #Create Line Plots of relaxation run. Create a figure.
|
---|
| 48 | fig = plt.figure(tight_layout=True)
|
---|
| 49 |
|
---|
| 50 | #Save surface mass balance, by looping through 200 years (1000 steps)
|
---|
| 51 | #Note, the first output will always contain output from time step 1
|
---|
| 52 |
|
---|
| 53 | Timer = np.arange(0.2, 200.2)
|
---|
| 54 | surfmb = []
|
---|
| 55 | for i in range(0, 200):
|
---|
| 56 | surfmb.append(md.results.TransientSolution[i].SmbMassBalance)
|
---|
| 57 |
|
---|
| 58 | #Plot surface mass balance time series in first subplot
|
---|
| 59 | ax = fig.add_subplot(311)
|
---|
| 60 | ax.plot(Timer, np.nanmean(surfmb, axis=1))
|
---|
| 61 |
|
---|
| 62 | #Title this plot Mean surface mass balance
|
---|
| 63 | ax.set_title('Mean Surface mass balance')
|
---|
| 64 |
|
---|
| 65 | #Save velocity by looping through 200 years
|
---|
| 66 | vel = []
|
---|
| 67 | for i in range(0, 200):
|
---|
| 68 | vel.append(md.results.TransientSolution[i].Vel)
|
---|
| 69 |
|
---|
| 70 | #Plot velocity time series in second subplot
|
---|
| 71 | ay = fig.add_subplot(312)
|
---|
| 72 | ay.plot(Timer, np.nanmean(vel, axis=1))
|
---|
| 73 |
|
---|
| 74 | #Title this plot Mean Velocity
|
---|
| 75 | ay.set_title('Mean Velocity')
|
---|
| 76 |
|
---|
| 77 | #Save Ice Volume by looping through 200 years
|
---|
| 78 | volume = []
|
---|
| 79 | for i in range(0, 200):
|
---|
| 80 | volume.append(md.results.TransientSolution[i].IceVolume)
|
---|
| 81 |
|
---|
| 82 | #Plot volume time series in third subplot
|
---|
| 83 | az = fig.add_subplot(313)
|
---|
| 84 | az.plot(Timer, volume)
|
---|
| 85 |
|
---|
| 86 | #Title this plot Mean Velocity and add an x label of years
|
---|
| 87 | az.set_title('Ice Volume')
|
---|
| 88 | az.set_xlabel('years')
|
---|
| 89 | # }}}
|
---|