diff --git a/conf/hermes.conf b/conf/hermes.conf index ad9ea14736848678cbfc107aa8a6a3e677db1d05..364be0d4d728c61fec214a8d511abe11a71fd1ec 100755 --- a/conf/hermes.conf +++ b/conf/hermes.conf @@ -9,15 +9,15 @@ start_date = 2016/11/29 00:00:00 # ----- end_date = start_date [DEFAULT] ----- # end_date = 2010/01/01 00:00:00 output_timestep_num = 24 -auxiliary_files_path = /scratch/Earth/HERMESv3_BU_aux/__test +auxiliary_files_path = /scratch/Earth/ctena/HERMESv3_BU_aux/__test first_time = 0 erase_auxiliary_files = 0 [DOMAIN] -# domain_type=[lcc, rotated, mercator, regular] -domain_type = lcc +# domain_type=[lcc, rotated, mercator, regular, rotated_nested] +domain_type = rotated_nested # output_type=[MONARCH, CMAQ, WRF_CHEM, DEFAULT] output_model = DEFAULT output_attributes = /writing/global_attributes_WRF-Chem.csv @@ -39,6 +39,16 @@ vertical_description = /profiles/vertical/MONARCH_Global_48layers_ver # inc_rlat = 0.4 # inc_rlon = 0.4 +# if domain_type == rotated_nested: + parent_grid_path = /scratch/Earth/ctena/temporal_coords.nc + parent_ratio = 4 + # i -> lon + i_parent_start = 200 + # j -> lat + j_parent_start = 100 + n_rlat = 400 + n_rlon = 400 + # if domain_type == lcc: # CALIOPE diff --git a/hermesv3_bu/config/config.py b/hermesv3_bu/config/config.py index f99e3e4072e4067c68e908aff90b853f4f46fcc5..4def4711b29de6b6c46f0adcb4eebcac09e72da0 100755 --- a/hermesv3_bu/config/config.py +++ b/hermesv3_bu/config/config.py @@ -74,12 +74,12 @@ class Config(ArgParser): p.add_argument('--output_attributes', required=False, help='Path to the file that contains the global attributes.') - p.add_argument('--domain_type', required=True, help='Type of domain to simulate.', - choices=['lcc', 'rotated', 'mercator', 'regular']) - p.add_argument('--vertical_description', required=True, help='Path to the file that contains the vertical description of the desired output.') + p.add_argument('--domain_type', required=True, help='Type of domain to simulate.', + choices=['lcc', 'rotated', 'mercator', 'regular', 'rotated_nested']) + # Rotated options p.add_argument('--centre_lat', required=False, type=float, help='Central geographic latitude of grid (non-rotated degrees). Corresponds to the TPH0D ' + @@ -100,6 +100,20 @@ class Config(ArgParser): help='Longitudinal grid resolution (rotated degrees). Corresponds to the DLMD parameter ' + 'in NMMB-MONARCH.') + # Rotated_nested options + p.add_argument('--parent_grid_path', required=False, type=str, + help='Path to the netCDF that contains the grid definition.') + p.add_argument('--parent_ratio', required=False, type=int, + help='Ratio between the parent and the nested domain.') + p.add_argument('--i_parent_start', required=False, type=int, + help='Location of the I to start the nested.') + p.add_argument('--j_parent_start', required=False, type=int, + help='Location of the J to start the nested.') + p.add_argument('--n_rlat', required=False, type=int, + help='Number of rotated latitude points.') + p.add_argument('--n_rlon', required=False, type=int, + help='Number of rotated longitude points.') + # Lambert conformal conic options p.add_argument('--lat_1', required=False, type=float, help='Standard parallel 1 (in deg). Corresponds to the P_ALP parameter of the GRIDDESC file.') @@ -686,6 +700,9 @@ class Config(ArgParser): elif arguments.domain_type == 'rotated': arguments.__dict__[item] = arguments.__dict__[item].replace('', '{1}_{2}'.format( item, arguments.inc_rlat, arguments.inc_rlon)) + elif arguments.domain_type == 'rotated_nested': + arguments.__dict__[item] = arguments.__dict__[item].replace('', '{1}_{2}'.format( + item, arguments.n_rlat, arguments.n_rlon)) elif arguments.domain_type == 'lcc' or arguments.domain_type == 'mercator': arguments.__dict__[item] = arguments.__dict__[item].replace('', '{1}_{2}'.format( item, arguments.inc_x, arguments.inc_y)) diff --git a/hermesv3_bu/grids/grid.py b/hermesv3_bu/grids/grid.py index 1f1e49b5e134bc8fba691c9e91093b92b844c6da..ff91a4acb648a3164fb9386401725d8a20152226 100755 --- a/hermesv3_bu/grids/grid.py +++ b/hermesv3_bu/grids/grid.py @@ -45,6 +45,12 @@ def select_grid(comm, logger, arguments): logger, arguments.auxiliary_files_path, arguments.output_timestep_num, arguments.vertical_description, arguments.centre_lat, arguments.centre_lon, arguments.west_boundary, arguments.south_boundary, arguments.inc_rlat, arguments.inc_rlon) + elif arguments.domain_type == 'rotated_nested': + from hermesv3_bu.grids.grid_rotated_nested import RotatedNestedGrid + grid = RotatedNestedGrid( + logger, arguments.auxiliary_files_path, arguments.output_timestep_num, + arguments.vertical_description, arguments.parent_grid_path, arguments.parent_ratio, + arguments.i_parent_start, arguments.j_parent_start, arguments.n_rlat, arguments.n_rlon) elif arguments.domain_type == 'mercator': from hermesv3_bu.grids.grid_mercator import MercatorGrid diff --git a/hermesv3_bu/grids/grid_rotated.py b/hermesv3_bu/grids/grid_rotated.py index 219570769a3dbd3dc6c1dedaa5bbe5c939bc8dc8..5d1fbec643c8ba30bc114ba7d24fac7eea592c12 100755 --- a/hermesv3_bu/grids/grid_rotated.py +++ b/hermesv3_bu/grids/grid_rotated.py @@ -172,7 +172,7 @@ class RotatedGrid(Grid): boundary_latitudes=self.boundary_latitudes, boundary_longitudes=self.boundary_longitudes, rotated=True, rotated_lats=self.rlat, rotated_lons=self.rlon, - north_pole_lat=90 - self.attributes['new_pole_latitude_degrees'], + north_pole_lat=self.attributes['new_pole_latitude_degrees'], north_pole_lon=self.attributes['new_pole_longitude_degrees']) self.logger.write_log("\tGrid created at '{0}'".format(self.netcdf_path), 3) self.logger.write_time_log('RotatedGrid', 'write_netcdf', timeit.default_timer() - spent_time, 3) diff --git a/hermesv3_bu/grids/grid_rotated_nested.py b/hermesv3_bu/grids/grid_rotated_nested.py new file mode 100755 index 0000000000000000000000000000000000000000..9b5ed392397f893d5608ffa91195fe1de37c18a0 --- /dev/null +++ b/hermesv3_bu/grids/grid_rotated_nested.py @@ -0,0 +1,199 @@ +#!/usr/bin/env python + +import os +import timeit +from hermesv3_bu.grids.grid import Grid +import numpy as np +import math + +from hermesv3_bu.logger.log import Log + + +class RotatedNestedGrid(Grid): + def __init__(self, logger, auxiliary_path, tstep_num, vertical_description_path, parent_grid_path, parent_ratio, + i_parent_start, j_parent_start, n_rlat, n_rlon): + """ + + :param logger: Logger. + :type logger: Log + + :param auxiliary_path: + :param tstep_num: + :param vertical_description_path: + + """ + spent_time = timeit.default_timer() + + self.rlat = None + self.rlon = None + + logger.write_log('Rotated Nested grid selected.') + self.grid_type = 'Rotated_nested' + + attributes = {'parent_grid_path': parent_grid_path, 'parent_ratio': parent_ratio, + 'i_parent_start': i_parent_start, 'j_parent_start': j_parent_start, + 'n_rlat': n_rlat, 'n_rlon': n_rlon, 'crs': {'init': 'epsg:4326'}} + attributes = self.get_parent_attributes(attributes) + + # Initialises with parent class + super(RotatedNestedGrid, self).__init__(logger, attributes, auxiliary_path, vertical_description_path) + + self.shape = (tstep_num, len(self.vertical_desctiption), attributes['n_rlat'], attributes['n_rlon']) + self.logger.write_time_log('RotatedNestedGrid', '__init__', timeit.default_timer() - spent_time, 3) + + @staticmethod + def get_parent_attributes(attributes): + from netCDF4 import Dataset + + netcdf = Dataset(attributes['parent_grid_path'], mode='r') + + rlat = netcdf.variables['rlat'][:] + attributes['inc_rlat'] = (rlat[1] - rlat[0]) / attributes['parent_ratio'] + # j_parent_start start the index at 1 so we must - 1 + attributes['1st_rlat'] = rlat[int(attributes['j_parent_start']) - 1] + + rlon = netcdf.variables['rlon'][:] + attributes['inc_rlon'] = (rlon[1] - rlon[0]) / attributes['parent_ratio'] + # i_parent_start start the index at 1 so we must - 1 + attributes['1st_rlon'] = rlon[attributes['i_parent_start'] - 1] + + rotated_pole = netcdf.variables['rotated_pole'] + attributes['new_pole_longitude_degrees'] = rotated_pole.grid_north_pole_longitude + attributes['new_pole_latitude_degrees'] = 90 - rotated_pole.grid_north_pole_latitude + + netcdf.close() + + return attributes + + def create_regular_rotated(self): + """ + Create a regular grid on the rotated domain. + + :return: center_latitudes, center_longitudes, corner_latitudes, corner_longitudes + :rtype: tuple + """ + spent_time = timeit.default_timer() + + center_latitudes = np.linspace(self.attributes['1st_rlat'], self.attributes['1st_rlat'] + + (self.attributes['inc_rlat'] * (self.attributes['n_rlat'] - 1)), + self.attributes['n_rlat'], dtype=np.float) + center_longitudes = np.linspace(self.attributes['1st_rlon'], self.attributes['1st_rlon'] + + (self.attributes['inc_rlon'] * (self.attributes['n_rlon'] - 1)), + self.attributes['n_rlon'], dtype=np.float) + + corner_latitudes = self.create_bounds(center_latitudes, self.attributes['inc_rlat'], number_vertices=4, + inverse=True) + corner_longitudes = self.create_bounds(center_longitudes, self.attributes['inc_rlon'], number_vertices=4) + + self.logger.write_time_log('RotatedNestedGrid', 'create_regular_rotated', + timeit.default_timer() - spent_time, 3) + return center_latitudes, center_longitudes, corner_latitudes, corner_longitudes + + def create_coords(self): + """ + Create the coordinates for a rotated domain. + """ + spent_time = timeit.default_timer() + # Create rotated coordinates + (self.rlat, self.rlon, br_lats_single, br_lons_single) = self.create_regular_rotated() + + # 1D to 2D + c_lats = np.array([self.rlat] * len(self.rlon)).T + c_lons = np.array([self.rlon] * len(self.rlat)) + + # Create rotated boundary coordinates + b_lats = super(RotatedNestedGrid, self).create_bounds(c_lats, self.attributes['inc_rlat'], number_vertices=4, + inverse=True) + b_lons = super(RotatedNestedGrid, self).create_bounds(c_lons, self.attributes['inc_rlon'], number_vertices=4) + + # Rotated to Lat-Lon + self.boundary_longitudes, self.boundary_latitudes = self.rotated2latlon(b_lons, b_lats) + self.center_longitudes, self.center_latitudes = self.rotated2latlon(c_lons, c_lats) + + self.logger.write_time_log('RotatedNestedGrid', 'create_coords', timeit.default_timer() - spent_time, 3) + return True + + def rotated2latlon(self, lon_deg, lat_deg, lon_min=-180): + """ + Calculate the unrotated coordinates using the rotated ones. + + :param lon_deg: Rotated longitude coordinate. + :type lon_deg: numpy.array + + :param lat_deg: Rotated latitude coordinate. + :type lat_deg: numpy.array + + :param lon_min: Minimum value for the longitudes: -180 (-180 to 180) or 0 (0 to 360) + :type lon_min: float + + :return: Unrotated coordinates. Longitudes, Latitudes + :rtype: tuple(numpy.array, numpy.array) + """ + spent_time = timeit.default_timer() + degrees_to_radians = math.pi / 180. + # radians_to_degrees = 180. / math.pi + + # Positive east to negative east + # self.new_pole_longitude_degrees -= 180 + + tph0 = self.attributes['new_pole_latitude_degrees'] * degrees_to_radians + tlm = lon_deg * degrees_to_radians + tph = lat_deg * degrees_to_radians + tlm0d = self.attributes['new_pole_longitude_degrees'] + ctph0 = np.cos(tph0) + stph0 = np.sin(tph0) + + stlm = np.sin(tlm) + ctlm = np.cos(tlm) + stph = np.sin(tph) + ctph = np.cos(tph) + + # Latitude + sph = (ctph0 * stph) + (stph0 * ctph * ctlm) + # if sph > 1.: + # sph = 1. + # if sph < -1.: + # sph = -1. + sph[sph > 1.] = 1. + sph[sph < -1.] = -1. + + aph = np.arcsin(sph) + aphd = aph / degrees_to_radians + + # Longitude + anum = ctph * stlm + denom = (ctlm * ctph - stph0 * sph) / ctph0 + relm = np.arctan2(anum, denom) - math.pi + almd = relm / degrees_to_radians + tlm0d + + # if almd < min_lon: + # almd += 360 + # elif almd > max_lon: + # almd -= 360 + almd[almd > (lon_min + 360)] -= 360 + almd[almd < lon_min] += 360 + + self.logger.write_time_log('RotatedNestedGrid', 'rotated2latlon', timeit.default_timer() - spent_time, 3) + + return almd, aphd + + def write_netcdf(self): + """ + Write a rotated grid NetCDF with empty data + """ + from hermesv3_bu.io_server.io_netcdf import write_coords_netcdf + spent_time = timeit.default_timer() + if not os.path.exists(self.netcdf_path): + if not os.path.exists(os.path.dirname(self.netcdf_path)): + os.makedirs(os.path.dirname(self.netcdf_path)) + # Writes an auxiliary empty NetCDF only with the coordinates and an empty variable. + write_coords_netcdf(self.netcdf_path, self.center_latitudes, self.center_longitudes, + [{'name': 'var_aux', 'units': '', 'data': 0}], + boundary_latitudes=self.boundary_latitudes, + boundary_longitudes=self.boundary_longitudes, + rotated=True, rotated_lats=self.rlat, rotated_lons=self.rlon, + north_pole_lat=90 - self.attributes['new_pole_latitude_degrees'], + north_pole_lon=self.attributes['new_pole_longitude_degrees']) + self.logger.write_log("\tGrid created at '{0}'".format(self.netcdf_path), 3) + self.logger.write_time_log('RotatedNestedGrid', 'write_netcdf', timeit.default_timer() - spent_time, 3) + return True diff --git a/hermesv3_bu/writer/default_writer.py b/hermesv3_bu/writer/default_writer.py index 5bde50bb8bcab1f262f3445ec05b92eee0f18cf8..1f2546fac076c3b5c94f298c55120cccc5e2e0a6 100755 --- a/hermesv3_bu/writer/default_writer.py +++ b/hermesv3_bu/writer/default_writer.py @@ -110,7 +110,7 @@ class DefaultWriter(Writer): var_dim = ('y', 'x',) lat_dim = lon_dim = var_dim - elif self.grid.grid_type == 'Rotated': + elif self.grid.grid_type in ['Rotated', 'Rotated_nested']: netcdf.createDimension('rlat', len(self.grid.rlat)) netcdf.createDimension('rlon', len(self.grid.rlon)) var_dim = ('rlat', 'rlon') @@ -177,7 +177,7 @@ class DefaultWriter(Writer): y_var.standard_name = "projection_y_coordinate" y_var[:] = self.grid.y - elif self.grid.grid_type == 'Rotated': + elif self.grid.grid_type in ['Rotated', 'Rotated_nested']: self.logger.write_log('\t\tCreating rlat variable', message_level=3) rlat = netcdf.createVariable('rlat', np.float64, ('rlat',)) rlat.long_name = "latitude in rotated pole grid" @@ -224,7 +224,7 @@ class DefaultWriter(Writer): var.grid_mapping = 'Latitude_Longitude' elif self.grid.grid_type == 'Lambert Conformal Conic': var.grid_mapping = 'Lambert_Conformal' - elif self.grid.grid_type == 'Rotated': + elif self.grid.grid_type in ['Rotated', 'Rotated_nested']: var.grid_mapping = 'rotated_pole' elif self.grid.grid_type == 'Mercator': var.grid_mapping = 'mercator' @@ -247,7 +247,7 @@ class DefaultWriter(Writer): mapping.longitude_of_central_meridian = self.grid.attributes['lon_0'] mapping.latitude_of_projection_origin = self.grid.attributes['lat_0'] - elif self.grid.grid_type == 'Rotated': + elif self.grid.grid_type in ['Rotated', 'Rotated_nested']: mapping = netcdf.createVariable('rotated_pole', 'c') mapping.grid_mapping_name = 'rotated_latitude_longitude' mapping.grid_north_pole_latitude = 90 - self.grid.attributes['new_pole_latitude_degrees'] diff --git a/hermesv3_bu/writer/monarch_writer.py b/hermesv3_bu/writer/monarch_writer.py index deb5d4cd24864f448bc6514bcd562877fffb1913..435811b2f1ae05a2864df7b7824f9e9b04b27a90 100755 --- a/hermesv3_bu/writer/monarch_writer.py +++ b/hermesv3_bu/writer/monarch_writer.py @@ -62,8 +62,8 @@ class MonarchWriter(Writer): super(MonarchWriter, self).__init__(comm_world, comm_write, logger, netcdf_path, grid, date_array, pollutant_info, rank_distribution, emission_summary) - if self.grid.grid_type not in ['Rotated']: - error_exit("ERROR: Only Rotated grid is implemented for MONARCH. " + + if self.grid.grid_type not in ['Rotated', 'Rotated_nested']: + error_exit("ERROR: Only Rotated or Rotated-nested grid is implemented for MONARCH. " + "The current grid type is '{0}'".format(self.grid.grid_type)) for i, (pollutant, variable) in enumerate(self.pollutant_info.iterrows()): diff --git a/hermesv3_bu/writer/writer.py b/hermesv3_bu/writer/writer.py index 0a92322b392c2c7906313218d4951cd46fc77242..f7ca8337f5c5b4cf912731f6d8db9dcf8443b050 100755 --- a/hermesv3_bu/writer/writer.py +++ b/hermesv3_bu/writer/writer.py @@ -14,6 +14,8 @@ CHUNKING = True BALANCED = False MPI_TAG_CONSTANT = 10**6 +BUFFER_SIZE = 2**28 + def select_writer(logger, comm_world, arguments, grid, date_array): """ @@ -328,7 +330,7 @@ class Writer(object): for i_rank in range(self.comm_world.Get_size()): self.logger.write_log( '\tFrom {0} to {1}'.format(i_rank, self.comm_world.Get_rank()), message_level=3) - req = self.comm_world.irecv(2**27, source=i_rank, tag=i_rank) + req = self.comm_world.irecv(BUFFER_SIZE, source=i_rank, tag=i_rank) dataframe = req.wait() data_list.append(dataframe.reset_index())