diff --git a/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py b/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py index eacdc0a323..780cec959b 100644 --- a/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py +++ b/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py @@ -132,12 +132,6 @@ def _remap_mali_topo(self): self._partition_scrip_file('mali.scrip.nc') out_mesh_name = self.mesh_name - out_descriptor = MpasCellMeshDescriptor( - filename='base_mesh.nc', - mesh_name=out_mesh_name) - out_descriptor.format = 'NETCDF3_64BIT' - out_descriptor.to_scrip('mpaso.scrip.nc') - self._partition_scrip_file('mpaso.scrip.nc') map_filename = \ f'map_{in_mesh_name}_to_{out_mesh_name}_mbtraave.nc' @@ -168,14 +162,16 @@ def _modify_mali_topo(self): g = constants['SHR_CONST_G'] - draft = - (ice_density / ocean_density) * thickness + draft = - (ice_density / ocean_density) * thickness + sea_level ice_mask = ds_mali.thickness > 0 - floating_mask = np.logical_and( - ice_mask, - ice_density / ocean_density * thickness <= sea_level - bed) + floating_mask = np.logical_and(ice_mask, draft > bed) grounded_mask = np.logical_and(ice_mask, np.logical_not(floating_mask)) + # draft is determined by the bed where the ice is grounded and by + # flotation where the ice is not grounded (floating or no ice) + draft = xr.where(grounded_mask, bed, draft) + if self.ocean_includes_grounded: ocean_mask = bed < sea_level else: @@ -191,10 +187,10 @@ def _modify_mali_topo(self): # we will remap conservatively ds_in = xr.Dataset() - ds_in['bed_elevation'] = bed - ds_in['landIceThkObserved'] = thickness - ds_in['landIcePressureObserved'] = lithop - ds_in['landIceDraftObserved'] = draft + ds_in['bed_elevation'] = ocean_frac * bed + ds_in['landIceThkObserved'] = ocean_frac * thickness + ds_in['landIcePressureObserved'] = ocean_frac * lithop + ds_in['landIceDraftObserved'] = ocean_frac * draft ds_in['landIceFrac'] = ice_frac ds_in['landIceGroundedFrac'] = grounded_frac ds_in['oceanFrac'] = ocean_frac @@ -211,10 +207,11 @@ def _create_mali_to_mpaso_weights(self, map_filename): logger = self.logger logger.info('Create weights file') + # target h5m file was created in parent class args = [ 'mbtempest', '--type', '5', '--load', f'mali.scrip.p{self.ntasks}.h5m', - '--load', f'mpaso.scrip.p{self.ntasks}.h5m', + '--load', f'target.scrip.p{self.ntasks}.h5m', '--file', map_filename, '--weights', '--gnomonic', '--boxeps', '1e-9', @@ -233,6 +230,9 @@ def _remap_mali_to_mpaso(self, map_filename): """ logger = self.logger logger.info('Remap to target') + config = self.config + section = config['remap_topography'] + renorm_threshold = section.getfloat('renorm_threshold') args = [ 'ncremap', @@ -252,6 +252,20 @@ def _remap_mali_to_mpaso(self, map_filename): var = np.minimum(var, 1.) ds_out[var_name] = var + # renormalize topography variables based on ocean fraction + ocean_frac = ds_out['oceanFrac'] + ocean_mask = ocean_frac > renorm_threshold + norm = xr.where(ocean_mask, 1.0 / ocean_frac, 0.0) + for var_name in [ + 'bed_elevation', + 'landIceThkObserved', + 'landIcePressureObserved', + 'landIceDraftObserved' + ]: + attrs = ds_out[var_name].attrs + ds_out[var_name] = ds_out[var_name] * norm + ds_out[var_name].attrs = attrs + write_netcdf(ds_out, 'mali_topography_remapped.nc') logger.info(' Done.') @@ -275,7 +289,10 @@ def _combine_topo(self): ds_out['landIceFloatingFracObserved'] = ( ds_mali['landIceFrac'] - - ds_mali['landIceGroundedFrac']) + ds_mali['landIceGroundedFrac'] + ) + + ds_out['landIceGroundedFracObserved'] = ds_mali['landIceGroundedFrac'] for var in ['maliFrac']: mali_field = ds_mali[var] @@ -289,12 +306,27 @@ def _combine_topo(self): alpha * ds_mali.bed_elevation + (1.0 - alpha) * ds_bedmachine.bed_elevation) + # our topography blending technique can lead to locations where the + # ice draft is slightly below the bed elevation; we correct that here + ds_out['landIceDraftObserved'] = xr.where( + ds_out['landIceDraftObserved'] < ds_out['bed_elevation'], + ds_out['bed_elevation'], + ds_out['landIceDraftObserved']) + alpha = ds_out.maliFrac # NOTE: MALI's ocean fraction is already scaled by the MALI fraction ds_out['oceanFracObserved'] = ( ds_mali.oceanFrac + (1.0 - alpha) * ds_bedmachine.oceanFracObserved) + # limit land ice floatation fraction to ocean fraction -- because of + # machine-precision noise in the combined ocean fraciton, + # land ice floating fraction can exceed ocean fraction by ~1e-15 + ds_out['landIceFloatingFracObserved'] = np.minimum( + ds_out['landIceFloatingFracObserved'], + ds_out['oceanFracObserved'] + ) + ds_out['ssh'] = ds_out.landIceDraftObserved write_netcdf(ds_out, 'topography_remapped.nc')