Added observations to plots of basal mass loss and area-averaged ice …

…shelf melt rate

aice_animation.py

@@ -14,16 +14,16 @@
 # Directory containing CICE output files
 directory = '/short/y99/kaa561/roms_spinup_newest/cice/'
 # Number of time indices in each file
-num_ts = [73]
+num_ts = [18, 270, 270, 270, 252]
 # File number to start with for the animation (1-based)
-start_file = 1
+start_file = 5
 # Degrees to radians conversion factor
 deg2rad = pi/180
 # Names of each month for making titles
 month_names = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December']
 
 # Read grid from the first file
-id = Dataset(directory + 'iceh1.nc', 'r')
+id = Dataset(directory + 'iceh' + str(start_file) + '.nc', 'r')
 lon = id.variables['TLON'][:,:]
 lat = id.variables['TLAT'][:,:]
 id.close()
@@ -62,10 +62,10 @@ def animate(i):
     img = contourf(x, y, data, lev, cmap='jet', extend='both')
     axis('off')
     # Add a title containing the date
-    title(str(time.day) + ' ' + month_names[time.month-1] + ' ' + str(time.year), fontsize=30)
+    title(str(time.day) + ' ' + month_names[time.month-1], fontsize=30) # + ' ' + str(time.year), fontsize=30)
     return img
 
 # Animate once every time index from start_file to the last file
-anim = FuncAnimation(fig, func=animate, frames=sum(array(num_ts[start_file-1:])))
+anim = FuncAnimation(fig, func=animate, frames=range(179,252)) #sum(array(num_ts[start_file-1:])))
 # Save as an mp4 with one frame per second
 anim.save('aice.mp4', fps=1)

circumpolar_plot.py

@@ -137,7 +137,7 @@ def circumpolar_plot (file_path, var_name, tstep, depth_key, depth, depth_bounds
             colour_map = 'jet'
     else:
         # Determine bounds automatically
-        if var_name in ['u', 'v', 'ubar', 'vbar', 'm', 'zeta', 'shflux', 'ssflux', 'sustr', 'svstr', 'bustr', 'bvstr']:
+        if var_name in ['u', 'v', 'ubar', 'vbar', 'm', 'shflux', 'ssflux', 'sustr', 'svstr', 'bustr', 'bvstr']:
             # Center levels on 0 for certain variables, with a blue-to-red
             # colourmap
             max_val = amax(abs(data))

dpt_timeseries.py

@@ -91,6 +91,7 @@ def dpt_timeseries (file_path):
     plot(time, transport)
     xlabel('Years')
     ylabel('Drake Passage Transport (Sv)')
+    grid(True)
     show(block=False)
     #savefig('dpt_timeseries.png')
 

file_guide.txt

@@ -108,18 +108,18 @@ plot_volume.py: Extracts the volume values written in the ocean.log files (you
 spinup_plots.py: Analyse a ROMS spinup by calculating and plotting 9 timeseries:
                  total heat content, total salt content, area-averaged ice shelf
 		 melt rate, ice shelf basal mass loss, total kinetic energy,
-		 maximum velocity, Drake Passage transport, total sea ice
-		 extent, and average bottom water temperature in ice shelf
-		 cavities. Also write the timeseries to a log file so they don't
-		 have to be recomputed if the run is extended; at the beginning
-		 of this script, previous values will be read from the same
-		 log file if it exists. Note that to run this file, density
-		 anomalies (variable name "rho") must be written to the ocean
-		 history/averages file. This option can be activated in the .in
-		 ROMS configuration file with Hout(idDano) and Aout(idDano).
-		 Otherwise, you can use make_density_file.py to calculate
-		 absolute density based on the seawater equation of state;
-		 however, this process is quite slow.
+		 Drake Passage transport, total sea ice extent, net sea ice-to-
+		 ocean freshwater flux, and area-averaged bottom water
+		 temperature in ice shelf cavities. Also write the timeseries to
+		 a log file so they don't have to be recomputed if the run is
+		 extended; at the beginning of this script, previous values will
+		 be read from the same log file if it exists. Note that to run
+		 this file, density anomalies (variable name "rho") must be
+		 written to the ocean history/averages file. This option can be
+		 activated in the .in ROMS configuration file with Hout(idDano)
+		 and Aout(idDano). Otherwise, you can use make_density_file.py
+		 to calculate absolute density based on the seawater equation of
+		 state; however, this process is quite slow.
 		 To run: First edit the i and j coordinates of the Drake Passage
 		         (in the function calc_drakepsgtrans) for your grid, and
 			 the value of rho0 (in the function get_rho) to match

h_circumpolar.py

@@ -32,7 +32,7 @@ def h_circumpolar (grid_path, fig_name):
     x = -(lat+90)*cos(lon*deg2rad+pi/2)
     y = (lat+90)*sin(lon*deg2rad+pi/2)
 
-    lev = arange(0,6000,500)
+    lev = arange(0,max(data),500)
 
     # Plot
     fig = figure(figsize=(16,12))
@@ -43,7 +43,8 @@ def h_circumpolar (grid_path, fig_name):
     title('Bathymetry (m)', fontsize=30)
     axis('off')
 
-    savefig(fig_name)
+    show()
+    #savefig(fig_name)
 
 
 # Command-line interface

massloss.py

@@ -25,7 +25,13 @@ def massloss (file_path, log_path):
     i_max = [350, 872, 975, 1030, 1090, 1155, 1190, 1240, 1369, 1443, 12]
     j_min = [1,   20,  150, 140,  160,  150,  187,  1,    65,   80,   100]
     j_max = [125, 123, 175, 160,  185,  200,  220,  135,  116,  150,  120]
-
+    # Observed mass loss (Rignot 2013) and uncertainty for each ice shelf
+    obs_massloss = [35.5, 47.7, 144.9, 101.2, 51.8, 89, 20.7, 155.4, 9.7, 22.5]
+    obs_massloss_error = [23, 34, 14, 8, 19, 17, 67, 45, 16, 12]
+    # Observed ice shelf melt rates and uncertainty
+    obs_ismr = [0.6, 0.1, 4.3, 16.2, 1.7, 3.8, 0.4, 0.3, 0.1, 0.5]
+    obs_ismr_error = [0.4, 0.1, 0.4, 1, 0.6, 0.7, 1, 0.1, 0.2, 0.2]
+
     # Density of ice in kg/m^3
     rho_ice = 916
 
@@ -82,7 +88,7 @@ def massloss (file_path, log_path):
 
     # Set up array of conversion factors from mass loss to area-averaged melt
     # rate for each ice shelf
-    factors = empty(12)
+    factors = empty(len(names))
 
     # Process each timestep separately to prevent memory overflow
     for t in range(start_t, size(time)):
@@ -120,18 +126,42 @@ def massloss (file_path, log_path):
     # Plot each timeseries
     print 'Plotting'
     for index in range(len(names)):
+        # Calculate the bounds on observed mass loss and melt rate
+        massloss_low = obs_massloss[index] - obs_massloss_error[index]
+        massloss_high = obs_massloss[index] + obs_massloss_error[index]
+        ismr_low = obs_ismr[index] - obs_ismr_error[index]
+        ismr_high = obs_ismr[index] + obs_ismr_error[index]
+        # Set up plot: mass loss and melt rate are directly proportional (with
+        # a different constant of proportionality for each ice shelf depending
+        # on its area) so plot one line with two y-axes
         fig, ax1 = subplots()
-        # Mass loss and melt rate are directly proportional (with a different
-        # constant of proportionality for each ice shelf depending on its area)
-        # so plot one line with two y-axes scales.
-        ax1.plot(time, massloss[index,:])        
+        ax1.plot(time, massloss[index,:], color='black')
+        # In blue, add dashed lines for observed mass loss
+        ax1.axhline(massloss_low, color='b', linestyle='dashed')
+        ax1.axhline(massloss_high, color='b', linestyle='dashed')
+        # Make sure y-limits won't cut off observed melt rate
+        ymin = min(0.95*ismr_low/factors[index], ax1.get_ylim()[0])
+        ymax = max(1.05*ismr_high/factors[index], ax1.get_ylim()[1])
+        ax1.set_ylim([ymin, ymax])
+        # Title and ticks in blue for this side of the plot
+        ax1.set_ylabel('Basal Mass Loss (Gt/y)', color='b')
+        for t1 in ax1.get_yticklabels():
+            t1.set_color('b')
         ax1.set_xlabel('Years')
-        ax1.set_ylabel('Basal Mass Loss (Gt/y)')
+        ax1.grid(True)
+        # Twin axis for melt rates
         ax2 = ax1.twinx()
-        ax2.set_ylabel('Area-Averaged Ice Shelf Melt Rate (m/y)')
         # Make sure the scales line up
-        limits = ax1.get_ylim()        
+        limits = ax1.get_ylim()
         ax2.set_ylim([limits[0]*factors[index], limits[1]*factors[index]])
+        # In red, add dashed lines for observed ice shelf melt rates
+        ax2.axhline(ismr_low, color='r', linestyle='dashed')
+        ax2.axhline(ismr_high, color='r', linestyle='dashed')
+        # Title and ticks in red for this side of the plot
+        ax2.set_ylabel('Area-Averaged Ice Shelf Melt Rate (m/y)', color='r')
+        for t2 in ax2.get_yticklabels():
+            t2.set_color('r')
+        # Name of the ice shelf for the main title
         title(names[index])
         fig.savefig(fig_names[index])
 

plot_kinetic_energy.py

@@ -69,6 +69,7 @@ def plot_kinetic_energy (files, dt, freq):
     plot(time, ke_all)
     xlabel('Years')
     ylabel('Kinetic energy')
+    grid(True)
     show()
 
 

plot_maxspeed.py

@@ -75,6 +75,7 @@ def plot_maxspeed (files, dt, freq):
     plot(time, maxv_all)
     xlabel('Years')
     ylabel('Maximum speed')
+    grid(True)
     show(block=False)
 
 

plot_volume.py

@@ -73,6 +73,7 @@ def plot_volume (files, dt, freq):
     plot(time, vol_all)
     xlabel('Years')
     ylabel('Volume Percent Anomalies')
+    grid(True)
     show()
 
 

repeat_forcing.py

@@ -74,12 +74,12 @@ def process (directory, head, tail, year_start, perday, var_list):
 # User parameters to edit here
 
 # Data where files <an_head>yyyy<tail> and <fc_head>yyyy<tail> exist
-directory = '../data/ERA_Interim/'
+directory = '/short/m68/kaa561/ROMS-CICE-MCT/data/ERA_Interim/subdaily/'
 # First part of filename for AN and FC files
 an_head = 'AN_'
 fc_head = 'FC_'
 # Last part of filename (common to both AN and FC)
-tail = '_unlim.nc'
+tail = '_subdaily.nc'
 # Year to build data from
 year_start = 1995
 # Number of records per day

romscice_atm_subdaily.py

@@ -3,17 +3,19 @@
 from scipy.interpolate import LinearNDInterpolator, RectBivariateSpline
 
 # Convert two ERA-Interim files:
-# AN_yyyy_unlim_orig.nc: one year of 6-hour measurements for surface pressure 
-#                        (sp), 2-metre temperature (t2m) and dew point (d2m),
-#                        total cloud cover (tcc), and 10-metre winds (u10, v10)
-# FC_yyyy_unlim_orig.nc: one year of 12-hour measurements for total 
-#                        precipitation (tp) 
+# AN_yyyy_subdaily_orig.nc: one year of 6-hour measurements for surface pressure
+#                           (sp), 2-metre temperature (t2m) and dew point (d2m),
+#                           total cloud cover (tcc), and 10-metre winds (u10, 
+#                           v10)
+# FC_yyyy_subdaily_orig.nc: one year of 12-hour measurements for total 
+#                           precipitation (tp) and snowfall (sf) 
 # to two ROMS-CICE input forcing files with the correct units and naming 
 # conventions:
-# AN_yyyy_unlim.nc: one year of 6-hour measurements for surface pressure
+# AN_yyyy_subdaily.nc: one year of 6-hour measurements for surface pressure
 #                   (Pair), temperature (Tair), specific humidity (Qair), 
 #                   cloud fraction (cloud), and winds (Uwind, Vwind)
-# FC_yyyy_unlim.nc: one year of 12-hour measurements for rainfall (rain)
+# FC_yyyy_subdaily.nc: one year of 12-hour measurements for rainfall (rain) and
+#                   snowfall (snow)
 # Input: year = integer containing the year to process
 #        count = time record in the given year to start with
 
@@ -30,11 +32,11 @@ def convert_file (year, count):
 
     # Paths of ROMS grid file, input ERA-Interim files, and output ROMS-CICE
     # files; other users will need to change these
-    grid_file = '/short/m68/kaa561/roms_circumpolar/data/caisom001_OneQuartergrd.nc'
-    input_atm_file = '../data/ERA_Interim/AN_' + str(year) + '_unlim_orig.nc'
-    input_ppt_file = '../data/ERA_Interim/FC_' + str(year) + '_unlim_orig.nc'
-    output_atm_file = '../data/ERA_Interim/AN_' + str(year) + '_unlim.nc'
-    output_ppt_file = '../data/ERA_Interim/FC_' + str(year) + '_unlim.nc'
+    grid_file = '/short/m68/kaa561/ROMS-CICE-MCT/apps/common/grid/circ30S_quarterdegree_rp5.nc'
+    input_atm_file = '/short/m68/kaa561/ROMS-CICE-MCT/data/ERA_Interim/subdaily/AN_' + str(year) + '_subdaily_orig.nc'
+    input_ppt_file = '/short/m68/kaa561/ROMS-CICE-MCT/data/ERA_Interim/subdaily/FC_' + str(year) + '_subdaily_orig.nc'
+    output_atm_file = '/short/m68/kaa561/ROMS-CICE-MCT/data/ERA_Interim/subdaily/AN_' + str(year) + '_subdaily.nc'
+    output_ppt_file = '/short/m68/kaa561/ROMS-CICE-MCT/data/ERA_Interim/subdaily/FC_' + str(year) + '_subdaily.nc'
     logfile = str(year) + '.log'
 
     Lv = 2.5e6 # Latent heat of vapourisation, J/kg
@@ -139,6 +141,9 @@ def convert_file (year, count):
         tair = interp_era2roms(t2m, lon_era, lat_era, lon_roms, lat_roms)
         oatm_fid.variables['Tair'][t,:,:] = tair-273.15
         qair = interp_era2roms(rh, lon_era, lat_era, lon_roms, lat_roms)
+        # Constrain humidity values to be between 0 and 1
+        qair[qair < 0] = 0.0
+        qair[qair > 1] = 1.0
         oatm_fid.variables['Qair'][t,:,:] = qair
         cloud = interp_era2roms(tcc, lon_era, lat_era, lon_roms, lat_roms)
         # Constrain cloud fractions to be between 0 and 1
@@ -183,6 +188,9 @@ def convert_file (year, count):
         oppt_fid.createVariable('rain', 'f8', ('time', 'eta_rho', 'xi_rho'))
         oppt_fid.variables['rain'].long_name = 'rain fall rate'
         oppt_fid.variables['rain'].units = 'm_per_12hr'
+        oppt_fid.createVariable('snow', 'f8', ('time', 'eta_rho', 'xi_rho'))
+        oppt_fid.variables['snow'].long_name = 'snow fall rate'
+        oppt_fid.variables['snow'].units = 'm_per_12hr'
         oppt_fid.close()
 
     log = open(logfile, 'a')
@@ -198,15 +206,19 @@ def convert_file (year, count):
         log.close()
         # Write the current time value to output FC file
         oppt_fid.variables['time'][t] = ppt_time[t]
-        # Read rainfall for this timestep
+        # Read data for this timestep
         ippt_fid = Dataset(input_ppt_file, 'r')
         tp = transpose(ippt_fid.variables['tp'][t,:,:])
+        sf = transpose(ippt_fid.variables['sf'][t,:,:])
         ippt_fid.close()
         # Interpolate to ROMS grid and write to output FC file
         rain = interp_era2roms(tp, lon_era, lat_era, lon_roms, lat_roms)
+        snow = interp_era2roms(sf, lon_era, lat_era, lon_roms, lat_roms)
         # Make sure there are no negative values
         rain[rain < 0] = 0.0
         oppt_fid.variables['rain'][t,:,:] = rain
+        snow[snow < 0] = 0.0
+        oppt_fid.variables['snow'][t,:,:] = snow
         oppt_fid.close()
 
     log = open(logfile, 'a')