Bug fix

calc_ice_prod.py

@@ -0,0 +1,93 @@
+from netCDF4 import Dataset, num2date
+from numpy import *
+
+def calc_ice_prod (in_file, start_year, end_year, out_file):
+
+    # cm/day to m/day conversion
+    cm_to_m = 1e-2
+
+    # Read the grid
+    id = Dataset(in_file, 'r')
+    lon = id.variables['TLON'][:,:]
+    lat = id.variables['TLAT'][:,:]
+    num_lon = size(lon,1)
+    num_lat = size(lat,0)
+    # Set up array to hold output data
+    ice_prod = ma.empty([num_lat, num_lon])
+    ice_prod[:,:] = 0.0
+
+    # Read the time values
+    time_id = id.variables['time']
+    cice_time = num2date(time_id[:], units=time_id.units, calendar='standard')
+    num_time = time_id.size
+    # Find the first timestep and how many days in it we care about
+    if start_year == 1992:
+        # Starting at the beginning of the simulation
+        start_t = 0
+        start_days = 5
+    else:
+        for t in range(num_time):
+            if cice_time[t].year == start_year and cice_time[t].month-1 == 0 and cice_time[t].day in range(2,6+1):
+                start_t = t
+                start_days = cice_time[t].day-1
+                break
+    # Find the last timestep and how many days in it we care about
+    if end_year == 2016:
+        # Ending at the end of the simulation
+        end_t = num_time-1
+        end_days = 5
+    else:
+        for t in range(num_time):
+            if cice_time[t].year == end_year+1 and cice_time[t].month-1 == 0 and cice_time[t].day in range(1,5+1):
+                end_t = t
+                end_days = 6-cice_time[t].day
+                break
+    print 'Starting at index ' + str(start_t) + ' (' + str(cice_time[start_t].year) + '-' + str(cice_time[start_t].month) + '-' + str(cice_time[start_t].day) + '), ' + str(start_days) + ' days'
+    print 'Ending at index ' + str(end_t) + ' (' + str(cice_time[end_t].year) + '-' + str(cice_time[end_t].month) + '-' + str(cice_time[end_t].day) + '), ' + str(end_days) + ' days'
+
+    # Integrate the start days
+    thdgr_start = id.variables['frazil'][start_t,:,:] + id.variables['congel'][start_t,:,:] - id.variables['meltt'][start_t,:,:] - id.variables['meltb'][start_t,:,:] - id.variables['meltl'][start_t,:,:]
+    index = thdgr_start < 0
+    thdgr_start[index] = 0
+    ice_prod += thdgr_start*cm_to_m*start_days
+    # Integrate the middle days
+    thdgr_mid = id.variables['frazil'][start_t+1:end_t,:,:] + id.variables['congel'][start_t+1:end_t,:,:] - id.variables['meltt'][start_t+1:end_t,:,:] - id.variables['meltb'][start_t+1:end_t,:,:] - id.variables['meltl'][start_t+1:end_t,:,:]
+    index = thdgr_mid < 0
+    thdgr_mid[index] = 0
+    ice_prod += sum(thdgr_mid, axis=0)*cm_to_m*5
+    # Integrate the end days
+    thdgr_end = id.variables['frazil'][end_t,:,:] + id.variables['congel'][end_t,:,:] - id.variables['meltt'][end_t,:,:] - id.variables['meltb'][end_t,:,:] - id.variables['meltl'][end_t,:,:]
+    index = thdgr_end < 0
+    thdgr_end[index] = 0
+    ice_prod += thdgr_end*cm_to_m*end_days
+    # Get annual average in m/y
+    ice_prod = ice_prod/(end_year-start_year+1)
+
+    # Write to file
+    id = Dataset(out_file, 'w')
+    id.createDimension('ni', size(lon,1))
+    id.createDimension('nj', size(lon,0))
+    id.createDimension('time', 4)
+    id.createVariable('TLON', 'f8', ('nj', 'ni'))
+    id.variables['TLON'][:,:] = lon
+    id.createVariable('TLAT', 'f8', ('nj', 'ni'))
+    id.variables['TLAT'][:,:] = lat
+    id.createVariable('ice_prod', 'f8', ('nj', 'ni'))
+    id.variables['ice_prod'].units = 'm/y'
+    id.variables['ice_prod'][:,:] = ice_prod
+    id.close()
+
+
+# Command-line interface
+if __name__ == "__main__":
+
+    in_file = raw_input("Path to CICE iceh_tot.nc file for entire simulation: ")
+    start_year = int(raw_input("First year to process: "))
+    end_year = int(raw_input("End year to process: "))
+    out_file = raw_input("Path to desired output file: ")
+    calc_ice_prod (in_file, start_year, end_year, out_file)
+
+
+
+
+

circumpolar_plot.py

@@ -264,51 +264,25 @@ def average_btw_depths (data_3d, z_3d, dz_3d, z_bounds):
                 # both depth bounds are too deep (i.e. in the seafloor)
                 data[j,i] = ma.masked
             else:
-                # Initialise integral and depth_range (vertical distance we're integrating over) to 0
-                integral = 0.0
-                depth_range = 0.0
+                # Find the first level above z_deep
                 if any(z_col < z_deep):
                     # There exist ocean cells below z_deep
-                    # Linearly interpolate to z_deep
-                    k_above_deep = nonzero(z_col > z_deep)[0][0]
-                    k_below_deep = k_above_deep - 1
-                    coeff1 = (z_col[k_below_deep] - z_deep)/(z_col[k_below_deep] - z_col[k_above_deep])
-                    coeff2 = 1 - coeff1
-                    data_deep = coeff1*data_col[k_above_deep] + coeff2*data_col[k_below_deep]
-                    # Now integrate between z_deep and z_col[k_above_deep]
-                    dz_curr = z_col[k_above_deep] - z_deep
-                    integral += 0.5*(data_deep + data_col[k_above_deep])*dz_curr # Trapezoidal rule
-                    depth_range += dz_curr
-                    # Save index of k_above_deep; we will start normal (non-interpolated) integration here
-                    k_start = k_above_deep
+                    k_start = nonzero(z_col > z_deep)[0][0]
                 else:
-                    # z_deep is deeper than the seafloor at this location
-                    # Start normal (non-interpolated) integration at the seafloor
-                    k_start = 0
+                    # z_deep is deeper than the seafloor at this location,
+                    # so start at the seafloor
+                    k_end = 0
+                # Find the first level above z_shallow                     
                 if any(z_col > z_shallow):
                     # There exist ocean cells above z_shallow
-                    # Linearly interpolate to z_shallow
-                    k_above_shallow = nonzero(z_col > z_shallow)[0][0]
-                    k_below_shallow = k_above_shallow - 1
-                    coeff1 = (z_col[k_below_shallow] - z_shallow)/(z_col[k_below_shallow] - z_col[k_above_shallow])
-                    coeff2 = 1 - coeff1
-                    data_shallow = coeff1*data_col[k_above_shallow] + coeff2*data_col[k_below_shallow]
-                    # Now integrate between z_col[k_below_shallow] and z_shallow
-                    dz_curr = z_shallow - z_col[k_below_shallow]
-                    integral += 0.5*(data_col[k_below_shallow] + data_shallow)*dz_curr # Trapezoidal rule
-                    depth_range += dz_curr
-                    # Save index of k_above_shallow; we will stop normal (non-interpolated) integration one level below it
-                    k_end = k_above_shallow
+                    k_end = nonzero(z_col > z_shallow)[0][0]
                 else:
-                    # z_shallow is shallower than the sea surface (or ice shelf draft) at this location
-                    # Continue normal (non-interpolated) integration all the way to the surface/ice shelf
+                    # z_shallow is shallower than the ice shelf draft at this location
+                    # Continue normal (non-interpolated) integration all the way to the ice draft
                     k_end = size(z_col)
                 # Now integrate between k_start and k_end
                 if k_start < k_end:
-                    integral += sum(data_col[k_start:k_end]*dz_col[k_start:k_end])
-                    depth_range += sum(dz_col[k_start:k_end])
-                # Divide integral by depth_range to get the vertical average
-                data[j,i] = integral/depth_range
+                    data[j,i] = sum(data_col[k_start:k_end]*dz_col[k_start:k_end])/sum(dz_col[k_start:k_end])
 
     return data    
 

ismr_plot.py

@@ -35,7 +35,7 @@ def ismr_plot (file_path, save=False, fig_name=None):
     zice = id.variables['zice'][:-15,:-1]
     # Read the last year of ice shelf melt rates (assume 5-day averages here),
     # average over time, and convert from m/s to m/y
-    ismr = mean(id.variables['m'][-73:,:-15,:-1], axis=0)*60*60*24*365.25
+    ismr = id.variables['m'][0,:-15,:-1]*60*60*24*365.25 #mean(id.variables['m'][-73:,:-15,:-1], axis=0)*60*60*24*365.25
     id.close()
     # Mask the open ocean and land out of the melt rates
     ismr = ma.masked_where(zice==0, ismr)