More ocean and ice shelf figures for intercomparison

circumpolar_plot.py

@@ -89,7 +89,9 @@ def circumpolar_plot (file_path, var_name, tstep, depth_key, depth, depth_bounds
                 u_data = id.variables[var_name.replace('v','u')][tstep-1,k,:-15,:]
                 u_data_lonlat, v_data_lonlat = rotate_vector_roms(u_data, v_data, angle)
                 data_full[k,:,:] = v_data_lonlat        
-
+
+    id.close()
+    id = Dataset(grid_path, 'r')
     # Read grid variables
     h = id.variables['h'][:-15,:]
     zice = id.variables['zice'][:-15,:]
@@ -359,7 +361,7 @@ def average_btw_depths (data_3d, z_3d, dz_3d, z_bounds):
         # Will need the grid file to get the angle
         grid_path = raw_input("Path to ROMS grid file: ")
     else:
-        grid_path = None
+        grid_path = raw_input("Path to ROMS grid file: ") #grid_path = None
 
     # Get index of time axis in ROMS history/averages file
     tstep = int(raw_input("Timestep number (starting at 1): "))

file_guide.txt

@@ -1165,15 +1165,16 @@ mip_timeseries.py: Plot MetROMS and FESOM timeseries together, for Drake
 			   plotted). The script will create a bunch of png
 			   files.
 
-mip_massloss_error_map.py: Make a 2x1 circumpolar Antarctic map of percentage
+mip_massloss_error_map.py: Make a 3x1 circumpolar Antarctic map of percentage
                            error in mass loss for each ice shelf, outside the
-			   bounds given by Rignot et al. (2013), in MetROMS
-			   (left) and FESOM (right), for the 2003-2008 average.
+			   bounds given by Rignot et al. (2013), in MetROMS,
+			   FESOM low-res, and FESOM high-res, for the 2003-2008
+			   average. Also print the values to the screen.
 			   To run: First make sure you have run
 			           timeseries_massloss.py for ROMS through to
 				   at least the end of 2008, and the equivalent
 				   timeseries_massloss.py in the "fesomtools"
-				   repository for the FESOM simulation. Then
+				   repository for each FESOM simulation. Then
 				   open python or ipython and type
 				   "run mip_massloss_error_map.py". The script
 				   will prompt you for the path to the ROMS
@@ -1184,23 +1185,22 @@ mip_massloss_error_map.py: Make a 2x1 circumpolar Antarctic map of percentage
 
 mip_massloss_difference_map.py: Make a circumpolar Antarctic figure showing the
                                 percentage change in mass loss for each ice
-				shelf in the FESOM simulation with respect to
-				the MetROMS simulation, for the 2003-2008
-				average.
+				shelf in the high-resolution FESOM simulation
+				with respect to the low-resolution FESOM
+				simulation, discarding the first 10 years (i.e.
+				2002-2016 average). Also print the values to
+				the screen.
 				To run: First make sure you have run
-			                timeseries_massloss.py for ROMS through
-					to at least the end of 2008, and the
-					equivalent timeseries_massloss.py in
-					the "fesomtools" repository for the
-					FESOM simulation. Then open python or
-					ipython and type
+			                timeseries_massloss.py in the
+					"fesomtools" repository for each FESOM
+					simulation. Then open python or ipython
+					and type
 					"run mip_massloss_error_map.py".
 					The script will prompt you for the path
-					to the ROMS grid file, the ROMS and
-					FESOM massloss logfiles, and whether
-					you want to save the figure (and if so,
-					what filename) or display it on the
-					screen.
+					to the ROMS grid file, the FESOM
+					massloss logfiles, and whether you want
+					to save the figure (and if so, what
+					filename) or display it on the screen.
 
 mip_cavity_fields.py: For each major ice shelf, make a 2x1 plot of the given
                       field for MetROMS (left) and FESOM (right), zoomed into
@@ -1347,23 +1347,34 @@ mip_ts_distribution.py: Make a 2x1 plot of T/S distributions south of 65S,
 				mesh directory, and the FESOM time-averaged
 				file. It will output a png.
 
-mip_drift_slices.py: Make a 4x2 plot showing zonal slices through 0E of
-                     temperature (left) and salinity (right), in MetROMS (top)
-		     and FESOM (bottom). The top 2x2 plots show the average
-		     over the last year of simulation, while the bottom 2x2
-		     plots show anomalies with respect to the first year.
+mip_drift_slices.py: Make a 3x2 plot of temperature (left) and salinity (right)
+                     through 0E. The top row is the initial conditions from
+		     ECCO2. The middle and bottom rows are the end of the
+		     simulation (last 5 day average) from MetROMS and FESOM
+		     respectively.
 		     To run: First clone my "fesomtools" GitHub repository and
 		             replace '/short/y99/kaa561/fesomtools' (near the
 			     top of the file) with the path to the cloned
-			     repository on your system. Next, make
-			     time-averaged files of temperature and salinity
-			     for each model over the first and last year (just
-			     use ncra). Then open python or ipython and type
-			     "run mip_drift_slices.py". The script will prompt
-			     you for the paths to the ROMS grid file, the ROMS
-			     time-averaged files, the FESOM mesh directory, and
-			     the FESOM time-averaged files. It will output a
-			     png.
+			     repository on your system. Also make sure the
+			     paths to the original ECCO2 files (ecco_temp_file,
+			     ecco_salt_file) are correct. Then open python or
+			     ipython and type "run mip_drift_slices.py". The
+			     script will prompt you for the paths to the ROMS
+			     grid file, the last ocean_avg file output by ROMS,
+			     the FESOM mesh directory, and the last oce.mean
+			     file output by FESOM. It will output a png.
+
+mip_calc_massloss.py: Calculate the average mass loss for each ice shelf over
+                      2002-2016 for MetROMS, FESOM low-res, and FESOM high-res.
+		      Print to the screen along with Rignot's estimates for
+		      2003-2008.
+		      To run: First make sure you have run
+		              timeseries_massloss.py for the MetROMS simulation,
+			      as well as the equivalent timeseries_massloss.py
+			      in the fesomtools repository for both FESOM
+			      simulations. Then open python or ipython and type
+			      "run mip_calc_massloss.py". The script will
+			      prompt you for the paths to the three logfiles.
 
 
 ***UTILITY FUNCTIONS***

mip_calc_massloss.py

@@ -0,0 +1,144 @@
+from numpy import *
+
+# Calculate the average mass loss for each ice shelf over 2002-2016 for MetROMS,
+# FESOM low-res, and FESOM high-res. Print to the screen along with Rignot's
+# estimates for 2003-2008.
+# Input:
+# roms_logfile = path to ROMS logfile from timeseries_massloss.py
+# fesom_logfile_lr, fesom_logfile_hr = paths to FESOM logfiles from
+#                   timeseries_massloss.py in the fesomtools repository, for
+#                   low-res and high-res respectively
+def mip_calc_massloss (roms_logfile, fesom_logfile_lr, fesom_logfile_hr):
+
+    # Year simulations start
+    year_start = 1992
+    # Years to avearge over
+    calc_start = 2002
+    calc_end = 2016
+    # Number of output steps per year in FESOM
+    peryear = 365/5
+    # Name of each ice shelf
+    names = ['Total Mass Loss', 'Larsen D Ice Shelf', 'Larsen C Ice Shelf', 'Wilkins & George VI & Stange Ice Shelves', 'Ronne-Filchner Ice Shelf', 'Abbot Ice Shelf', 'Pine Island Glacier Ice Shelf', 'Thwaites Ice Shelf', 'Dotson Ice Shelf', 'Getz Ice Shelf', 'Nickerson Ice Shelf', 'Sulzberger Ice Shelf', 'Mertz Ice Shelf', 'Totten & Moscow University Ice Shelves', 'Shackleton Ice Shelf', 'West Ice Shelf', 'Amery Ice Shelf', 'Prince Harald Ice Shelf', 'Baudouin & Borchgrevink Ice Shelves', 'Lazarev Ice Shelf', 'Nivl Ice Shelf', 'Fimbul & Jelbart & Ekstrom Ice Shelves', 'Brunt & Riiser-Larsen Ice Shelves', 'Ross Ice Shelf']
+    # Observed mass loss (Rignot 2013) and uncertainty for each ice shelf, in Gt/y
+    obs_massloss = [1325, 1.4, 20.7, 135.4, 155.4, 51.8, 101.2, 97.5, 45.2, 144.9, 4.2, 18.2, 7.9, 90.6, 72.6, 27.2, 35.5, -2, 21.6, 6.3, 3.9, 26.8, 9.7, 47.7]
+    obs_massloss_error = [235, 14, 67, 40, 45, 19, 8, 7, 4, 14, 2, 3, 3, 8, 15, 10, 23, 3, 18, 2, 2, 14, 16, 34]
+    num_shelves = len(obs_massloss)
+
+    # Read ROMS logfile
+    roms_time = []
+    f = open(roms_logfile, 'r')
+    # Skip the first line (header for time array)
+    f.readline()
+    for line in f:
+        try:
+            roms_time.append(float(line))
+        except(ValueError):
+            # Reached the header for the next variable
+            break       
+    # Set up array for mass loss values at each ice shelf
+    roms_massloss_ts = empty([num_shelves, len(roms_time)])
+    index = 0
+    # Loop over ice shelves
+    while index < num_shelves:
+        t = 0
+        for line in f:
+            try:
+                roms_massloss_ts[index, t] = float(line)
+                t += 1
+            except(ValueError):
+                # Reached the header for the next ice shelf
+                break
+        index +=1
+    f.close()
+    # Add start year to ROMS time array
+    roms_time = array(roms_time) + year_start
+    # Average between given years
+    t_start = nonzero(roms_time >= calc_start)[0][0]
+    if calc_end == 2016:
+        t_end = size(roms_time)
+    else:
+        t_end = nonzero(roms_time >= calc_end+1)[0][0]
+    roms_massloss = mean(roms_massloss_ts[:,t_start:t_end], axis=1)
+
+    # Read FESOM timeseries
+    # Low-res
+    tmp = []
+    f = open(fesom_logfile_lr, 'r')
+    # Skip the first line (header)
+    f.readline()
+    # Read total mass loss
+    num_time = 0
+    for line in f:
+        try:
+            tmp.append(float(line))
+            num_time += 1
+        except(ValueError):
+            # Reached the header for the next variable
+            break
+    # Set up array for mass loss values at each ice shelf
+    fesom_massloss_ts_lr = empty([num_shelves, num_time])
+    # Save the total values in the first index
+    fesom_massloss_ts_lr[0,:] = array(tmp)
+    # Loop over ice shelves
+    index = 1
+    while index < num_shelves:
+        t = 0
+        for line in f:
+            try:
+                fesom_massloss_ts_lr[index,t] = float(line)
+                t += 1
+            except(ValueError):
+                # Reached the header for the next ice shelf
+                break
+        index += 1
+    f.close()
+    # Average between given years
+    fesom_massloss_lr = mean(fesom_massloss_ts_lr[:,peryear*(calc_start-year_start):peryear*(calc_end+1-year_start)], axis=1)
+    # Repeat for high-res
+    tmp = []
+    f = open(fesom_logfile_hr, 'r')
+    f.readline()
+    num_time = 0
+    for line in f:
+        try:
+            tmp.append(float(line))
+            num_time += 1
+        except(ValueError):
+            break
+    fesom_massloss_ts_hr = empty([num_shelves, num_time])
+    fesom_massloss_ts_hr[0,:] = array(tmp)
+    index = 1
+    while index < num_shelves:
+        t = 0
+        for line in f:
+            try:
+                fesom_massloss_ts_hr[index,t] = float(line)
+                t += 1
+            except(ValueError):
+                break
+        index += 1
+    f.close()
+    fesom_massloss_hr = mean(fesom_massloss_ts_hr[:,peryear*(calc_start-year_start):peryear*(calc_end+1-year_start)], axis=1)
+
+    # Loop over ice shelves
+    for index in range(num_shelves):
+        print names[index]
+        print 'MetROMS: ' + str(roms_massloss[index])
+        print 'FESOM low-res: ' + str(fesom_massloss_lr[index])
+        print 'FESOM high-res: ' + str(fesom_massloss_hr[index])
+        # Find the range of observations
+        massloss_low = obs_massloss[index] - obs_massloss_error[index]
+        massloss_high = obs_massloss[index] + obs_massloss_error[index]
+        print 'Rignot: ' + str(massloss_low) + '-' + str(massloss_high)
+
+
+# Command-line interface
+if __name__ == "__main__":
+
+    roms_logfile = raw_input("Path to ROMS logfile from timeseries_massloss.py: ")
+    fesom_logfile_lr = raw_input("Path to FESOM low-res logfile from timeseries_massloss.py: ")
+    fesom_logfile_hr = raw_input("Path to FESOM high-res logfile from timeseries_massloss.py: ")
+    mip_calc_massloss(roms_logfile, fesom_logfile_lr, fesom_logfile_hr)
+
+
+