From bd6d7bec2d81d25efaf6c3543fd8c253be42e031 Mon Sep 17 00:00:00 2001 From: Kaitlin Naughten Date: Fri, 2 Jun 2017 16:32:36 +1000 Subject: [PATCH] Added flux correction calculation --- file_guide.txt | 16 +++++ iceberg_melt.py | 14 +--- romscice_flux_correction.py | 135 ++++++++++++++++++++++++++++++++++++ 3 files changed, 154 insertions(+), 11 deletions(-) create mode 100644 romscice_flux_correction.py diff --git a/file_guide.txt b/file_guide.txt index fcd4643..aed8250 100644 --- a/file_guide.txt +++ b/file_guide.txt @@ -308,6 +308,22 @@ romscice_sss_nudging.py: Interpolate the World Ocean Atlas 2013 monthly file). Then open python or ipython and type "run romscice_sss_nudging.py". +romscice_flux_correction.py: Build a monthly climatology of the extra surface + salt flux due to surface salinity restoring in a + long ROMS run. Save as a NetCDF file that can be + used as a forcing file for another ROMS run + (option SSFLUX_EXTRA) with surface salinity + restoring off. + To run: First make sure you have concatenated the + ROMS averages files from the long + simulation into a single file, or at least + the ssflux_restoring variable. Then open + python or ipython and type + "run romscice_flux_correction.py". The + script will prompt you for the paths to + the concatenated ROMS averages file and + the desired output forcing file. + ***POST-PROCESSING DIAGNOSTIC FILES*** diff --git a/iceberg_melt.py b/iceberg_melt.py index 879ab3f..ef93029 100644 --- a/iceberg_melt.py +++ b/iceberg_melt.py @@ -16,12 +16,6 @@ def iceberg_melt (out_file): iceberg_grid = '/short/m68/kaa561/metroms_iceshelf/data/originals/MartinAdcroft2010_iceberg_meltfluxes/icebergs.static.nc' # ROMS grid file roms_grid ='/short/m68/kaa561/metroms_iceshelf/apps/common/grid/circ30S_quarterdegree.nc' - # Density of freshwater - rho_w = 1e3 - # Seconds in 12 hours - seconds_per_12h = 60.*60.*12. - # Days per month (neglect leap years) - days_per_month = array([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]) # Read ROMS grid id = Dataset(roms_grid, 'r') @@ -61,7 +55,7 @@ def iceberg_melt (out_file): out_id.variables['time'].cycle_length = 365.25 out_id.createVariable('icebergs', 'f8', ('time', 'eta_rho', 'xi_rho')) out_id.variables['icebergs'].long_name = 'freshwater flux from iceberg melt' - out_id.variables['icebergs'].units = 'm_per_12hr' + out_id.variables['icebergs'].units = 'kg/m^2/s' # Loop over months for month in range(12): @@ -78,10 +72,8 @@ def iceberg_melt (out_file): id.close() # Interpolate to ROMS grid melt_roms = interp_iceberg2roms(melt_iceberg, lon_iceberg, lat_iceberg, lon_roms, lat_roms) - # Convert to m per 12 h - melt_roms = melt_roms/rho_w*seconds_per_12h - # Calculate time values: 15th of this month at midnight - out_id.variables['time'][month] = sum(days_per_month[0:month]) + 14 + # Calculate time values centered in the middle of each month + out_id.variables['time'][month] = 365.25/12*(month+0.5) # Save data out_id.variables['icebergs'][month,:,:] = melt_roms out_id.close() diff --git a/romscice_flux_correction.py b/romscice_flux_correction.py new file mode 100644 index 0000000..7b260d9 --- /dev/null +++ b/romscice_flux_correction.py @@ -0,0 +1,135 @@ +from netCDF4 import Dataset, num2date +from numpy import * + +# Build a monthly climatology of the extra surface salt flux due to surface +# salinity restoring in a long ROMS run. Save as a NetCDF file that can be +# used as a forcing file for another ROMS run (option SSFLUX_EXTRA) with +# surface salinity restoring off. +# Input: +# in_file = path to ROMS ocean averages file concatenated for the entire +# simulation (or just the variable ssflux_restoring extracted) +# out_file = desired path to output ROMS forcing file +def romscice_flux_correction (in_file, out_file): + + # Starting and ending days for each month + # Ignore leap years, they will be dealt with later + start_day = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] + end_day = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] + + print 'Reading ' + in_file + id = Dataset(in_file, 'r') + lon = id.variables['lon_rho'][:,:] + lat = id.variables['lat_rho'][:,:] + # Convert time to Date objects + time_id = id.variables['ocean_time'] + time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) + restoring = id.variables['ssflux_restoring'][:,:,:] + id.close() + + print 'Calculating climatology' + num_lat = size(lat,0) + num_lon = size(lon,1) + # Initialise integral arrays + climatology = ma.empty([12, num_lat, num_lon]) + climatology[:,:,:] = 0.0 + num_days = zeros(12) + # Loop over timesteps + for t in range(size(time)): + print '...time index ' + str(t+1) + ' of ' + str(size(time)) + # 5-day averages marked with middle day's date + year = time[t].year + month = time[t].month-1 + day = time[t].day + # Check for leap years + leap_year = False + if mod(year, 4) == 0: + leap_year = True + if mod(year, 100) == 0: + leap_year = False + if mod(year, 400) == 0: + leap_year = True + if leap_year: + end_day[1] = 29 + else: + end_day[1] = 28 + # Integrate ssflux weighted by days for the correct month(s) + if day-2 >= start_day[month] and day+2 <= end_day[month]: + # The 5 days in this index are entirely within 1 month + climatology[month,:,:] += restoring[t,:,:]*5 + # Also integrate the number of days + num_days[month] += 5 + elif day-2 < start_day[month]: + # Partially spills over into the previous month + prev_month = mod(month-1,12) + # How many days does it spill over by? + spill_days = start_day[month]-day+2 + # Should be either 1 or 2 + if spill_days < 1 or spill_days > 2: + print 'Problem: spill_days is ' + str(spill_days) + print 'Timestep ' + str(t+1) + print 'Year ' + str(year) + print 'Month ' + str(month+1) + print 'Day ' + str(day) + return + # Split between previous month and this month appropriately + climatology[prev_month,:,:] += restoring[t,:,:]*spill_days + num_days[prev_month] += spill_days + climatology[month,:,:] += restoring[t,:,:]*(5-spill_days) + num_days[month] += 5-spill_days + elif day+2 > end_day[month]: + # Partially spills over into the next month + next_month = mod(month+1,12) + # How many days does it spill over by? + spill_days = day+2-end_day[month] + # Should be either 1 or 2 + if spill_days < 1 or spill_days > 2: + print 'Problem: spill_days is ' + str(spill_days) + print 'Timestep ' + str(t+1) + print 'Year ' + str(year) + print 'Month ' + str(month+1) + print 'Day ' + str(day) + return + # Split between previous month and this month appropriately + climatology[next_month,:,:] += restoring[t,:,:]*spill_days + num_days[next_month] += spill_days + climatology[month,:,:] += restoring[t,:,:]*(5-spill_days) + num_days[month] += 5-spill_days + # Make sure we counted days correctly + if sum(num_days) != size(time)*5: + print 'Problem with the number of days: found ' + str(num_days) + ' instead of ' + sum(size(time)*5) + # Convert from sum to averages + for month in range(12): + climatology[month,:,:] /= num_days[month] + + print 'Writing ' + out_file + id = Dataset(out_file, 'w') + id.createDimension('xi_rho', num_lon) + id.createDimension('eta_rho', num_lat) + id.createDimension('time', None) + id.createVariable('lon_rho', 'f8', ('eta_rho', 'xi_rho')) + id.variables['lon_rho'].long_name = 'longitude of rho-points' + id.variables['lon_rho'].units = 'degree_east' + id.variables['lon_rho'][:,:] = lon + id.createVariable('lat_rho', 'f8', ('eta_rho', 'xi_rho')) + id.variables['lat_rho'].long_name = 'latitude of rho-points' + id.variables['lat_rho'].units = 'degree_north' + id.variables['lat_rho'][:,:] = lat + id.createVariable('time', 'f8', ('time')) + id.variables['time'].units = 'days since 1992-01-01 00:00:0.0' + # Add a cycle length so it repeats every year + id.variables['time'].cycle_length = 365.25 + # Time values in the middle of each month + id.variables['time'][:] = 365.25/12*(arange(12)+0.5) + id.createVariable('ssflux_extra', 'f8', ('time', 'eta_rho', 'xi_rho')) + id.variables['ssflux_extra'].long_name = 'additional surface salt flux' + id.variables['ssflux_extra'].units = 'kg/m^2/s' + id.variables['ssflux_extra'][:,:,:] = climatology + id.close() + + +# Command-line interface +if __name__ == "__main__": + + in_file = raw_input("Path to ROMS ocean averages file containing ssflux_restoring for the entire simulation: ") + out_file = raw_input("Path to desired output forcing file: ") + romscice_flux_correction(in_file, out_file)