From febf9fcbd2e5821fe596eee3d669c8235e28f50a Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Mon, 6 May 2024 14:19:14 -0700 Subject: [PATCH 01/39] add interpolate_mpasOcean_to_mali Creates new script that interpolates MPAS-Ocean thermal forcing and sub ice shelf melt rates to a MALI mesh for one way coupling. Still a work in progress at this point --- .../interpolate_mpasOcean_to_mali.py | 142 ++++++++++++++++++ 1 file changed, 142 insertions(+) create mode 100644 landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py new file mode 100644 index 000000000..2b1ad930d --- /dev/null +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -0,0 +1,142 @@ +#!/usr/bin/env python + +import os +import shutil +import xarray as xr +import subprocess + +from mpas_tools.logging import check_call +from mpas_tools.scrip.from_mpas import scrip_from_mpas +from optparse import OptionParser + +parser = OptionParser() +parser.add_option("--oceanMesh", dest="mpasMeshFile", help="MPAS-Ocean mesh to be interpolated from", metavar="FILENAME") +parser.add_option("--oceanDiag", dest="mpasDiagFile", help="MPAS-Ocean diagnostic output file", metavar="FILENAME") +parser.add_option("--ice", dest="maliFile", help="MALI mesh to be interpolated onto", metavar="FILENAME") +parser.add_option("-n", "--ntasks", dest="ntasks", help="Number of processors to use with ESMF_RegridWeightGen") +parser.add_option("-m", "--method", dest="method", help="Remapping method, either 'bilinear' or 'neareststod'") +parser.add_option("-o", "--outFile",dest="outputFile", help="Desired name of output file", metavar="FILENAME", default="mpas_to_mali_remapped.nc") +options, args = parser.parse_args() + +#make copy of mpasMeshFile +tmp_mpasMeshFile = "tmp_mpasMeshFile.nc" +subprocess.run(["ncks", options.mpasMeshFile, tmp_mpasMeshFile]) + +#Preprocessing: +#TO DO: Create capability to read multiple files and time average over relavant fields + +#TO DO: Vertically interpolate between MPAS-O to MALI grids + +#main: +#calculate thermal forcing +calc_ocean_thermal_forcing(options.mpasMeshFile, options.mpasDiagFile) + +#TO DO: copy landIceFreshwaterFlux to mpas mesh file copy for interpolation. Make name relevant to MALI + +#remap mpas to mali +remap_mpas_to_mali(options.mpasMeshFile, options.maliMeshFile, options.ntasks, options.method, options.outputFile) + +def calc_ocean_thermal_forcing(tmp_mpasMeshFile, mpasDiagFile) + print("Calculating Thermal Forcing") + print("... gathering data ...") + fM = xr.open_dataset(options.mpasMeshFile, decode_times=False, decode_cf=False) + minLevelCell = fM['minLevelCell'][:].data + maxLevelCell = fM['maxLevelCell'][:].data + layerThickness = fM['layerThickness'][:,:,:].data + landIceFloatingMask = fM['landIceFloatingMask'][:,:].data + nCells = fM.sizes['nCells'] + + fD = xr.open_dataset(options.mpasDiagFile, decode_times=False, decode_cf=False) + temperature = fD['timeMonthly_avg_activeTracers_temperature'][:,:,:].data + salinity = fD['timeMonthly_avg_activeTracers_salinity'][:,:,:].data + density = fD['time_monthy_avg_density'][:,:,:].data + atmPressure = fD['time_monthly_avg_atmosphericPressure'][:,:].data + coeff_0_openOcean = fD.attrs['config_open_ocean_freezing_temperature_coeff_0'] + coeff_S_openOcean = fD.attrs['config_open_ocean_freezing_temperature_coeff_S'] + coeff_p_openOcean = fD.attrs['config_open_ocean_freezing_temperature_coeff_p'] + coeff_pS_openOcean = fD.attrs['config_open_ocean_freezing_temperature_coeff_pS'] + coeff_mushy_az1_liq = fD.attrs['config_open_ocean_freezing_temperature_coeff_mushy_az1_liq'] + coeff_mushy_openOcean = 1/coeff_mushy_az1_liq + coeff_0_cavity = fD.attrs['config_land_ice_cavity_freezing_temperature_coeff_0'] + coeff_S_cavity = fD.attrs['config_land_ice_cavity_freezing_temperature_coeff_S'] + coeff_p_cavity = fD.attrs['config_land_ice_cavity_temperature_coeff_p'] + coeff_pS_cavity = fD.attrs['config_land_ice_cavity_freezing_temperature_coeff_pS'] + coeff_mushy_az1_liq = fD.attrs['config_land_ice_cavity_freezing_temperature_coeff_mushy_az1_liq'] + coeff_mushy_cavity = 1/coeff_mushy_az1_liq + Time = fD.sizes['Time'] + + gravity = 9.81 + + for iTime in range(Time-1): + + #calculate pressure: NEED TO CHECK KMIN AND K INDEXING + for iCell in range(nCells-1): + kmin = minLevelCell[iCell] + kmax = maxLevelCell[iCell] + pressure[iTime,kmin,iCell] = atmPressure[iTime,iCell] + density[iTime,iCell,kmin]*gravity*0.5*layerThickness[iTime,iCell,kmin] + + for k in np.arange(kmin, kmax): + pressure[iTime,k,iCell] = pressure[iTime,iCell,k] + 0.5*gravity*(density[iTime,iCell,k]*layerThickness[iTime,iCell,k] + density[iTime,iCell,k-1]*layerThickness[iTime,iCell,k-1]) + + if (landIceFloatingMask[iTime,iCell] == 1): + ocn_freezing_temperature[iTime,iCell,:] = coeff_0_openOcean + coeff_S_openOcean * salinity[iTime,iCell,:] + coeff_p_openOcean * pressure[iTime,iCell,:] \ + + coeff_pS_openOcean * pressure[iTime,iCell,:] * salinity[iTime,iCell,:] + coeff_mushy_openOcean * salinity[iTime,iCell,:] / (1.0 - salinity[iTime,iCell,:] / 1e3) + + elif (landIceFloatingMask[iTime,iCell] == 0): + ocn_freezing_temperature[iTime,iCell,:] = coeff_0_cavity + coeff_S_cavity * salinity[iTime,iCell,:] + coeff_p_cavity * pressure[iTime,iCell,:] \ + + coeff_pS_cavity * pressure[iTime,iCell,:] * salinity[iTime,iCell,:] + coeff_mushy_cavity * salinity[iTime,iCell,:] / (1.0 - salinity[iTime,iCell,:] / 1e3) + + # Calculate thermal forcing + oceanThermalForcing = temperature - ocn_freezing_temperature + + # Save thermal forcing to mesh file + #TO DO: Change name to variable that MALI will recognize + tf = xr.DataArray(oceanThermalForcing, dims=('Time','nCells','nVertLevels') + fM['oceanThermalForcing'] = tf + fM.to_netcdf(tmp_mpasMeshFile) + fM.close() + +def remap_mpas_to_mali(mpas_meshFile, mali_meshFile, ntasks, method, outputFile) + #create scrip files + print("Creating Scrip Files") + mali_scripfile = "tmp_mali_scrip.nc" + mpas_scripfile = "tmp_mpas_scrip.nc" + + scrip_from_mpas(mali_meshFile, mali_scripfile) + scrip_from_mpas(mpas_meshFile, mpas_scripfile) + + #create mapping file + print("Creating Mapping File") + args_esmf = ['srun', + '-n', ntasks, 'ESMF_RegridWeightGen', + '--source', mpas_scripfile, + '--destination', mali_scripfile, + '--weight', "tmp_mapping_file.nc", + '--method', method, + '-i', '-64bit_offset', + "--dst_regional", "--src_regional", '--ignore_unmapped'] + check_call(args_esmf) + + #TO DO: Create mask indentifying valid overlapping ocean cells + + #permute to work with ncremap + shutil.copy(mpas_meshFile, "tmp_" + mpas_meshFile) + tmp_mpasFile = "tmp_" + mpas_meshFile + + subprocess.run(["ncpdq", "-O", "-a", "Time,nVertLevels,nCells", tmp_mpasFile, tmp_mpasFile]) + subprocess.run(["ncpdq", "-O", "-a", "Time,nVertLevelsP1,nCells", tmp_mpasFile, tmp_mpasFile]) + subprocess.run(["ncpdq", "-O", "-a", "Time,nVertInterfaces,nCells", tmp_mpasFile, tmp_mpasFile]) + + # remap the input data + args_remap = ["ncremap", + "-i", tmp_mpasFile, + "-o", outputFile, + "-m", "tmp_mapping_file.nc"] + check_call(args_remap) + + #remove temporary files + os.remove(mpas_scripfile) + os.remove(mali_scripfile) + os.remove(tmp_mpasFile) + os.remove("tmp_mapping_file.nc") + From 55a569ba0e8f73b9a40cee0ffa518ed570170623 Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Tue, 7 May 2024 17:03:39 -0700 Subject: [PATCH 02/39] thermal forcing calc. and reorganize Implements the thermal forcing calculation, as well as reorganizes script to object oriented. Further work is labeled with "<< TO DO >>:" --- .../interpolate_mpasOcean_to_mali.py | 304 ++++++++++-------- 1 file changed, 174 insertions(+), 130 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index 2b1ad930d..9d1ac011e 100644 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -4,139 +4,183 @@ import shutil import xarray as xr import subprocess - +import numpy as np +from scipy import stats as st from mpas_tools.logging import check_call from mpas_tools.scrip.from_mpas import scrip_from_mpas from optparse import OptionParser -parser = OptionParser() -parser.add_option("--oceanMesh", dest="mpasMeshFile", help="MPAS-Ocean mesh to be interpolated from", metavar="FILENAME") -parser.add_option("--oceanDiag", dest="mpasDiagFile", help="MPAS-Ocean diagnostic output file", metavar="FILENAME") -parser.add_option("--ice", dest="maliFile", help="MALI mesh to be interpolated onto", metavar="FILENAME") -parser.add_option("-n", "--ntasks", dest="ntasks", help="Number of processors to use with ESMF_RegridWeightGen") -parser.add_option("-m", "--method", dest="method", help="Remapping method, either 'bilinear' or 'neareststod'") -parser.add_option("-o", "--outFile",dest="outputFile", help="Desired name of output file", metavar="FILENAME", default="mpas_to_mali_remapped.nc") -options, args = parser.parse_args() - -#make copy of mpasMeshFile -tmp_mpasMeshFile = "tmp_mpasMeshFile.nc" -subprocess.run(["ncks", options.mpasMeshFile, tmp_mpasMeshFile]) - -#Preprocessing: -#TO DO: Create capability to read multiple files and time average over relavant fields - -#TO DO: Vertically interpolate between MPAS-O to MALI grids - -#main: -#calculate thermal forcing -calc_ocean_thermal_forcing(options.mpasMeshFile, options.mpasDiagFile) - -#TO DO: copy landIceFreshwaterFlux to mpas mesh file copy for interpolation. Make name relevant to MALI - -#remap mpas to mali -remap_mpas_to_mali(options.mpasMeshFile, options.maliMeshFile, options.ntasks, options.method, options.outputFile) - -def calc_ocean_thermal_forcing(tmp_mpasMeshFile, mpasDiagFile) - print("Calculating Thermal Forcing") - print("... gathering data ...") - fM = xr.open_dataset(options.mpasMeshFile, decode_times=False, decode_cf=False) - minLevelCell = fM['minLevelCell'][:].data - maxLevelCell = fM['maxLevelCell'][:].data - layerThickness = fM['layerThickness'][:,:,:].data - landIceFloatingMask = fM['landIceFloatingMask'][:,:].data - nCells = fM.sizes['nCells'] - - fD = xr.open_dataset(options.mpasDiagFile, decode_times=False, decode_cf=False) - temperature = fD['timeMonthly_avg_activeTracers_temperature'][:,:,:].data - salinity = fD['timeMonthly_avg_activeTracers_salinity'][:,:,:].data - density = fD['time_monthy_avg_density'][:,:,:].data - atmPressure = fD['time_monthly_avg_atmosphericPressure'][:,:].data - coeff_0_openOcean = fD.attrs['config_open_ocean_freezing_temperature_coeff_0'] - coeff_S_openOcean = fD.attrs['config_open_ocean_freezing_temperature_coeff_S'] - coeff_p_openOcean = fD.attrs['config_open_ocean_freezing_temperature_coeff_p'] - coeff_pS_openOcean = fD.attrs['config_open_ocean_freezing_temperature_coeff_pS'] - coeff_mushy_az1_liq = fD.attrs['config_open_ocean_freezing_temperature_coeff_mushy_az1_liq'] - coeff_mushy_openOcean = 1/coeff_mushy_az1_liq - coeff_0_cavity = fD.attrs['config_land_ice_cavity_freezing_temperature_coeff_0'] - coeff_S_cavity = fD.attrs['config_land_ice_cavity_freezing_temperature_coeff_S'] - coeff_p_cavity = fD.attrs['config_land_ice_cavity_temperature_coeff_p'] - coeff_pS_cavity = fD.attrs['config_land_ice_cavity_freezing_temperature_coeff_pS'] - coeff_mushy_az1_liq = fD.attrs['config_land_ice_cavity_freezing_temperature_coeff_mushy_az1_liq'] - coeff_mushy_cavity = 1/coeff_mushy_az1_liq - Time = fD.sizes['Time'] - - gravity = 9.81 - - for iTime in range(Time-1): - - #calculate pressure: NEED TO CHECK KMIN AND K INDEXING - for iCell in range(nCells-1): - kmin = minLevelCell[iCell] - kmax = maxLevelCell[iCell] - pressure[iTime,kmin,iCell] = atmPressure[iTime,iCell] + density[iTime,iCell,kmin]*gravity*0.5*layerThickness[iTime,iCell,kmin] - - for k in np.arange(kmin, kmax): - pressure[iTime,k,iCell] = pressure[iTime,iCell,k] + 0.5*gravity*(density[iTime,iCell,k]*layerThickness[iTime,iCell,k] + density[iTime,iCell,k-1]*layerThickness[iTime,iCell,k-1]) - - if (landIceFloatingMask[iTime,iCell] == 1): - ocn_freezing_temperature[iTime,iCell,:] = coeff_0_openOcean + coeff_S_openOcean * salinity[iTime,iCell,:] + coeff_p_openOcean * pressure[iTime,iCell,:] \ - + coeff_pS_openOcean * pressure[iTime,iCell,:] * salinity[iTime,iCell,:] + coeff_mushy_openOcean * salinity[iTime,iCell,:] / (1.0 - salinity[iTime,iCell,:] / 1e3) - - elif (landIceFloatingMask[iTime,iCell] == 0): - ocn_freezing_temperature[iTime,iCell,:] = coeff_0_cavity + coeff_S_cavity * salinity[iTime,iCell,:] + coeff_p_cavity * pressure[iTime,iCell,:] \ - + coeff_pS_cavity * pressure[iTime,iCell,:] * salinity[iTime,iCell,:] + coeff_mushy_cavity * salinity[iTime,iCell,:] / (1.0 - salinity[iTime,iCell,:] / 1e3) - - # Calculate thermal forcing - oceanThermalForcing = temperature - ocn_freezing_temperature +class mpasToMaliInterp: + + def __init__(self): + parser = OptionParser() + parser.add_option("--oceanMesh", dest="mpasMeshFile", help="MPAS-Ocean mesh to be interpolated from", metavar="FILENAME") + parser.add_option("--oceanDiags", dest="mpasDiagsDir", help="Directory path where MPAS-Ocean diagnostic files are stored. Should be only netcdf files in that directory", metavar="FILENAME") + parser.add_option("--ice", dest="maliFile", help="MALI mesh to be interpolated onto", metavar="FILENAME") + parser.add_option("-n", "--ntasks", dest="ntasks", help="Number of processors to use with ESMF_RegridWeightGen") + parser.add_option("-m", "--method", dest="method", help="Remapping method, either 'bilinear' or 'neareststod'") + parser.add_option("-o", "--outFile",dest="outputFile", help="Desired name of output file", metavar="FILENAME", default="mpas_to_mali_remapped.nc") + parser.add_option("-a","--yearlyAvg", dest="yearlyAvg", help="true or false option to average monthly output to yearly intervals", default="true") + self.options, args = parser.parse_args() + # << TO DO >>: Get rid of options sub-class, define options directly in self class + + # open and concatenate diagnostic dataset + DS = xr.open_mfdataset(self.options.mpasDiagsDir + '/' + '*.nc', combine='nested', concat_dim='Time', decode_timedelta=False) + #open mpas ocean mesh + OM = xr.open_dataset(tmp_mpasMeshFile, decode_times=False, decode_cf=False) + #open MALI mesh + IM = xr.open_dataset(options.maliMeshFile, decode_times=False, decode_cf=False) + + # variables for time averaging + # << TO DO >>: Adjust for variables with no time dimension after open_mfdataset (e.g., Matt's areaCell example) + self.temperature = DS['timeMonthly_avg_activeTracers_temperature'][:,:,:].data + self.salinity = DS['timeMonthly_avg_activeTracers_temperature'][:,:,:].data + self.density = DS['timeMonthly_avg_density'][:,:,:].data + self.daysSinceStart = DS['daysSinceStart'][:].data + self.stTime = DS.attrs['config_start_time'] + self.landIceFloatingMask = OM['landIceFloatingMiask'][:,:,:].data + self.layerThickness = OM['layerThickness'][:,:].data + # << TO DO >> : add landIceFreshwaterFlux field + # << TO DO >>: Figure out how to import xtime with xarray + + # additional variables for computing thermal forcing + self.minLevelCell = OM['minLevelCell'][:].data + self.maxLevelCell = OM['maxLevelCell'][:].data + self.nCells = OM.sizes['nCells'] + self.coeff_0_openOcean = DS.attrs['config_open_ocean_freezing_temperature_coeff_0'] + self.coeff_S_openOcean = DS.attrs['config_open_ocean_freezing_temperature_coeff_S'] + self.coeff_p_openOcean = DS.attrs['config_open_ocean_freezing_temperature_coeff_p'] + self.coeff_pS_openOcean = DS.attrs['config_open_ocean_freezing_temperature_coeff_pS'] + self.coeff_mushy_az1_liq = DS.attrs['config_open_ocean_freezing_temperature_coeff_mushy_az1_liq'] + self.coeff_mushy_openOcean = 1/coeff_mushy_az1_liq + self.coeff_0_cavity = DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_0'] + self.coeff_S_cavity = DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_S'] + self.coeff_p_cavity = DS.attrs['config_land_ice_cavity_temperature_coeff_p'] + self.coeff_pS_cavity = DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_pS'] + self.coeff_mushy_az1_liq = DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_mushy_az1_liq'] + self.coeff_mushy_cavity = 1/coeff_mushy_az1_liq + self.Time = DS.sizes['Time'] + + def time_average_output(self) + if (self.options.yearlyAvg == 'true'): + yearsSinceStart = self.daysSinceStart / 365.0 + finalYear = np.floor(np.max(yearsSinceStart)) + timeStride = 1 # 1 year average + + years = np.arange(0, finalYear, timeStride) + for i in range(len(years)-1): + ind = np.where(yearsSinceStart >= years(i) and yearsSinceStart < years(i+1)) + self.newTemp[i,:,:] = np.mean(self.temperature[ind,:,:], axis=0) + self.newSal[i,:,:] = np.mean(self.salinity[ind,:,:], axis=0) + self.newDens[i,:,:] = np.mean(self.density[ind,:,:], axis=0) + self.newLThick[i,:,:] = np.mean(self.layerThickness[ind,:,:], axis=0) + self.newAtmPr[i,:] = np.mean(self.atmPressure[ind,:], axis=0) + self.newFloatIceMask[i,:] = st.mode(self.landIceFloatingMask[ind,:], axis=0) + + # Build new xtime string. Defined on first day of each year + stTime = np.datetime64(self.stTime[0:4]) + stYear = np.datetime_as_string(stTime,unit='Y').astype('float64') + yearVec = yearsSinceStart + stYear + self.newXtime = np.array([f"{int(year)}-01-01_00:00:00" for year in yearVec],dtype=np.str_) + else #do nothing if not time averaging + self.newTemp = self.temperature + self.newSal = self.salinity + self.newDens= self.density + self.newLThick = self.layerThickness + self.newAtmPr = self.atmPressure + self.newFloatIceMask = self.landIceFloatingMask + # << TO DO >>: Figure out what to do with xtime if no change. How do we get xarray to load original xtime variable? - # Save thermal forcing to mesh file - #TO DO: Change name to variable that MALI will recognize - tf = xr.DataArray(oceanThermalForcing, dims=('Time','nCells','nVertLevels') - fM['oceanThermalForcing'] = tf - fM.to_netcdf(tmp_mpasMeshFile) - fM.close() - -def remap_mpas_to_mali(mpas_meshFile, mali_meshFile, ntasks, method, outputFile) - #create scrip files - print("Creating Scrip Files") - mali_scripfile = "tmp_mali_scrip.nc" - mpas_scripfile = "tmp_mpas_scrip.nc" - - scrip_from_mpas(mali_meshFile, mali_scripfile) - scrip_from_mpas(mpas_meshFile, mpas_scripfile) - - #create mapping file - print("Creating Mapping File") - args_esmf = ['srun', - '-n', ntasks, 'ESMF_RegridWeightGen', - '--source', mpas_scripfile, - '--destination', mali_scripfile, - '--weight', "tmp_mapping_file.nc", - '--method', method, - '-i', '-64bit_offset', - "--dst_regional", "--src_regional", '--ignore_unmapped'] - check_call(args_esmf) - - #TO DO: Create mask indentifying valid overlapping ocean cells - - #permute to work with ncremap - shutil.copy(mpas_meshFile, "tmp_" + mpas_meshFile) - tmp_mpasFile = "tmp_" + mpas_meshFile - - subprocess.run(["ncpdq", "-O", "-a", "Time,nVertLevels,nCells", tmp_mpasFile, tmp_mpasFile]) - subprocess.run(["ncpdq", "-O", "-a", "Time,nVertLevelsP1,nCells", tmp_mpasFile, tmp_mpasFile]) - subprocess.run(["ncpdq", "-O", "-a", "Time,nVertInterfaces,nCells", tmp_mpasFile, tmp_mpasFile]) - - # remap the input data - args_remap = ["ncremap", - "-i", tmp_mpasFile, - "-o", outputFile, - "-m", "tmp_mapping_file.nc"] - check_call(args_remap) - - #remove temporary files - os.remove(mpas_scripfile) - os.remove(mali_scripfile) - os.remove(tmp_mpasFile) - os.remove("tmp_mapping_file.nc") + def calc_ocean_thermal_forcing(self) + + gravity = 9.81 + + for iTime in range(self.Time-1): # << TO DO >>: "Time" here needs to change to be consistent with time averaging + + #calculate pressure: + # << TO DO >>: NEED TO CHECK KMIN AND K INDEXING + ## << TO DO >>: pre-allocate ocn_freezing_temperature + for iCell in range(self.nCells-1): + kmin = self.minLevelCell[iCell] + kmax = self.maxLevelCell[iCell] + self.newPr[iTime,kmin,iCell] = self.newAtmPr[iTime,iCell] + self.newDens[iTime,iCell,kmin]*gravity*0.5*self.newLThick[iTime,iCell,kmin] + + for k in np.arange(kmin, kmax): + self.newPr[iTime,k,iCell] = self.newPr[iTime,iCell,k] + 0.5*gravity*(self.newDens[iTime,iCell,k]*self.newLThick[iTime,iCell,k] \ + + self.newDens[iTime,iCell,k-1]*self.newLThick[iTime,iCell,k-1]) + + if (self.newFloatIceMask[iTime,iCell] == 1): + ocn_freezing_temperature[iTime,iCell,:] = self.coeff_0_openOcean + self.coeff_S_openOcean * self.newSal[iTime,iCell,:] + self.coeff_p_openOcean * self.newPr[iTime,iCell,:] \ + + self.coeff_pS_openOcean * self.newPr[iTime,iCell,:] * self.newSal[iTime,iCell,:] + self.coeff_mushy_openOcean * self.newSal[iTime,iCell,:] / (1.0 - self.newSal[iTime,iCell,:] / 1e3) + + elif (self.newFloatingIceMask[iTime,iCell] == 0): + ocn_freezing_temperature[iTime,iCell,:] = self.coeff_0_cavity + self.coeff_S_cavity * self.newSal[iTime,iCell,:] + self.coeff_p_cavity * self.newPr[iTime,iCell,:] \ + + self.coeff_pS_cavity * self.newPr[iTime,iCell,:] * self.newSal[iTime,iCell,:] + self.coeff_mushy_cavity * self.newSal[iTime,iCell,:] / (1.0 - self.newSal[iTime,iCell,:] / 1e3) + + # Calculate thermal forcing + self.oceanThermalForcing = self.newTemp - ocn_freezing_temperature + + def remap_mpas_to_mali(self) + + # make copy of mpasMeshFile to save interpolated variable to + tmp_mpasFile = "tmp_" + self.options.mpasMeshFile + shutil.copy(self.options.mpasMeshFile, tmp_mpasMeshFile) + # << TO DO >>: populate tmp_mpasMeshFile with variables to be interpolated + + #create scrip files + print("Creating Scrip Files") + mali_scripfile = "tmp_mali_scrip.nc" + mpas_scripfile = "tmp_mpas_scrip.nc" + + scrip_from_mpas(self.options.maliFile, mali_scripfile) + scrip_from_mpas(tmp_mpasMeshFile, mpas_scripfile) + + #create mapping file + print("Creating Mapping File") + args_esmf = ['srun', + '-n', self.options.ntasks, 'ESMF_RegridWeightGen', + '--source', mpas_scripfile, + '--destination', mali_scripfile, + '--weight', "tmp_mapping_file.nc", + '--method', self.options.method, + '-i', '-64bit_offset', + "--dst_regional", "--src_regional", '--ignore_unmapped'] + check_call(args_esmf) + + #<< TO DO >>: Create mask indentifying valid overlapping ocean cells + + subprocess.run(["ncpdq", "-O", "-a", "Time,nVertLevels,nCells", tmp_mpasMeshFile, tmp_mpasMeshFile]) + subprocess.run(["ncpdq", "-O", "-a", "Time,nVertLevelsP1,nCells", tmp_mpasMeshFile, tmp_mpasMeshFile]) + subprocess.run(["ncpdq", "-O", "-a", "Time,nVertInterfaces,nCells", tmp_mpasMeshFile, tmp_mpasMeshFile]) + + # remap the input data + args_remap = ["ncremap", + "-i", tmp_mpasMeshFile, + "-o", self.options.outputFile, + "-m", "tmp_mapping_file.nc"] + check_call(args_remap) + + #remove temporary files + os.remove("tmp_*.nc") + +def main(): + run = mpasToMaliInterp() + + #compute yearly time average if necessary + run.time_average_output(self) + + #calculate thermal forcing + run.calc_ocean_thermal_forcing(self) + + # << TO DO >>: copy landIceFreshwaterFlux to mpas mesh file copy for interpolation. Make name relevant to MALI + + # << TO DO >>: Vertically interpolate between MPAS-O to MALI grids + + #remap mpas to mali + run.remap_mpas_to_mali(self) + +if __name__ == "__main__": + main() + + From c0df6c69d472e76a817fccbbfaa4d3c54021654b Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Fri, 10 May 2024 11:22:43 -0700 Subject: [PATCH 03/39] Finish interpolate_mpasOcean_to_mali Finishes to-do list of interpolation script. Script capabilities are now all complete, but script still needs to be tested and likely debugged --- .../interpolate_mpasOcean_to_mali.py | 105 +++++++++++++----- 1 file changed, 75 insertions(+), 30 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index 9d1ac011e..2ce84bc91 100644 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -8,6 +8,7 @@ from scipy import stats as st from mpas_tools.logging import check_call from mpas_tools.scrip.from_mpas import scrip_from_mpas +from mpas_tools.ocean.depth import compute_zmid from optparse import OptionParser class mpasToMaliInterp: @@ -21,27 +22,32 @@ def __init__(self): parser.add_option("-m", "--method", dest="method", help="Remapping method, either 'bilinear' or 'neareststod'") parser.add_option("-o", "--outFile",dest="outputFile", help="Desired name of output file", metavar="FILENAME", default="mpas_to_mali_remapped.nc") parser.add_option("-a","--yearlyAvg", dest="yearlyAvg", help="true or false option to average monthly output to yearly intervals", default="true") - self.options, args = parser.parse_args() - # << TO DO >>: Get rid of options sub-class, define options directly in self class + options, args = parser.parse_args() + self.options = options # open and concatenate diagnostic dataset DS = xr.open_mfdataset(self.options.mpasDiagsDir + '/' + '*.nc', combine='nested', concat_dim='Time', decode_timedelta=False) #open mpas ocean mesh - OM = xr.open_dataset(tmp_mpasMeshFile, decode_times=False, decode_cf=False) + OM = xr.open_dataset(self.options.mpasMeshFile, decode_times=False, decode_cf=False) #open MALI mesh - IM = xr.open_dataset(options.maliMeshFile, decode_times=False, decode_cf=False) + IM = xr.open_dataset(self.options.maliMeshFile, decode_times=False, decode_cf=False) # variables for time averaging - # << TO DO >>: Adjust for variables with no time dimension after open_mfdataset (e.g., Matt's areaCell example) self.temperature = DS['timeMonthly_avg_activeTracers_temperature'][:,:,:].data self.salinity = DS['timeMonthly_avg_activeTracers_temperature'][:,:,:].data self.density = DS['timeMonthly_avg_density'][:,:,:].data - self.daysSinceStart = DS['daysSinceStart'][:].data + self.daysSinceStart = DS['timeMonthly_avg_daysSinceStartOfSim'][:].data self.stTime = DS.attrs['config_start_time'] - self.landIceFloatingMask = OM['landIceFloatingMiask'][:,:,:].data - self.layerThickness = OM['layerThickness'][:,:].data - # << TO DO >> : add landIceFreshwaterFlux field - # << TO DO >>: Figure out how to import xtime with xarray + self.landIceFloatingMask = OM['landIceFloatingMask'][:,:,:].data + self.layerThickness = DS['timeMonthly_avg_layerThickness'][:,:].data + xt = DS['xtime'][:,:].data + xtime = np.array([xt],dtype=('S',64)) + # << NOTE >>: may need to eventually use: "xtime.data.tobytes().decode()" but not sure yet + try: + self.landIceFreshwaterFlux = DS['timeMonthly_avg_landIceFreshwaterFlux'][:,:].data + finally: + print("No landIceFreshwaterFlux variable") + self.landIceFreshwaterFlux = 0 # additional variables for computing thermal forcing self.minLevelCell = OM['minLevelCell'][:].data @@ -59,8 +65,12 @@ def __init__(self): self.coeff_pS_cavity = DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_pS'] self.coeff_mushy_az1_liq = DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_mushy_az1_liq'] self.coeff_mushy_cavity = 1/coeff_mushy_az1_liq - self.Time = DS.sizes['Time'] + #variables for interpolation << NOTE >>: a couple of these fields are redundant but we want them in xarray form for compute_zmid + self.bottomDepth_dataArray = OM['bottomDepth'][:] + self.maxLevelCell_dataArray = OM['maxLevelCell'][:] + self.layerThickness_dataArray = DS['timeMonthly_avg_layerThickness'] + def time_average_output(self) if (self.options.yearlyAvg == 'true'): yearsSinceStart = self.daysSinceStart / 365.0 @@ -76,6 +86,10 @@ def time_average_output(self) self.newLThick[i,:,:] = np.mean(self.layerThickness[ind,:,:], axis=0) self.newAtmPr[i,:] = np.mean(self.atmPressure[ind,:], axis=0) self.newFloatIceMask[i,:] = st.mode(self.landIceFloatingMask[ind,:], axis=0) + if (self.landIceFreshwaterFlux != 0): + self.newLandIceFWFlux[i,:] = np.mean(self.landIceFreshwaterFlux[ind,:], axis=0) + elif + self.newLandIceFWFlux = 0 # Build new xtime string. Defined on first day of each year stTime = np.datetime64(self.stTime[0:4]) @@ -89,27 +103,28 @@ def time_average_output(self) self.newLThick = self.layerThickness self.newAtmPr = self.atmPressure self.newFloatIceMask = self.landIceFloatingMask - # << TO DO >>: Figure out what to do with xtime if no change. How do we get xarray to load original xtime variable? + self.newLandIceFWFlux = self.landIceFreshwaterFlux + self.newXtime = xtime def calc_ocean_thermal_forcing(self) gravity = 9.81 + nt,nc,nz = self.newTemp.shape - for iTime in range(self.Time-1): # << TO DO >>: "Time" here needs to change to be consistent with time averaging + ocn_freezing_temperature = np.zeros((nt,nc,nz)) + for iTime in range(nt): #calculate pressure: - # << TO DO >>: NEED TO CHECK KMIN AND K INDEXING - ## << TO DO >>: pre-allocate ocn_freezing_temperature - for iCell in range(self.nCells-1): - kmin = self.minLevelCell[iCell] - kmax = self.maxLevelCell[iCell] + for iCell in range(self.nCells): + kmin = self.minLevelCell[iCell] - 1 + kmax = self.maxLevelCell[iCell] - 1 self.newPr[iTime,kmin,iCell] = self.newAtmPr[iTime,iCell] + self.newDens[iTime,iCell,kmin]*gravity*0.5*self.newLThick[iTime,iCell,kmin] - for k in np.arange(kmin, kmax): - self.newPr[iTime,k,iCell] = self.newPr[iTime,iCell,k] + 0.5*gravity*(self.newDens[iTime,iCell,k]*self.newLThick[iTime,iCell,k] \ - + self.newDens[iTime,iCell,k-1]*self.newLThick[iTime,iCell,k-1]) + for k in np.arange(kmin + 1, kmax): + self.newPr[iTime,k,iCell] = self.newPr[iTime,iCell,k-1] + 0.5*gravity*(self.newDens[iTime,iCell,k-1]*self.newLThick[iTime,iCell,k-1] \ + + self.newDens[iTime,iCell,k]*self.newLThick[iTime,iCell,k]) - if (self.newFloatIceMask[iTime,iCell] == 1): + if (self.newFloatingIceMask[iTime,iCell] == 1): ocn_freezing_temperature[iTime,iCell,:] = self.coeff_0_openOcean + self.coeff_S_openOcean * self.newSal[iTime,iCell,:] + self.coeff_p_openOcean * self.newPr[iTime,iCell,:] \ + self.coeff_pS_openOcean * self.newPr[iTime,iCell,:] * self.newSal[iTime,iCell,:] + self.coeff_mushy_openOcean * self.newSal[iTime,iCell,:] / (1.0 - self.newSal[iTime,iCell,:] / 1e3) @@ -125,7 +140,19 @@ def remap_mpas_to_mali(self) # make copy of mpasMeshFile to save interpolated variable to tmp_mpasFile = "tmp_" + self.options.mpasMeshFile shutil.copy(self.options.mpasMeshFile, tmp_mpasMeshFile) - # << TO DO >>: populate tmp_mpasMeshFile with variables to be interpolated + + # populate tmp_mpasMeshFile with variables to be interpolated + f = xr.open_dataset(tmp_mpasMeshFile, decode_times=False, decode_cf=False) + tf = xr.DataArray(self.oceanThermalForcing.astype('float64'),dims=('Time','nCells','nVertLevels')) + f['ismip6shelfMelt_3dThermalForcing'] = tf + + if (self.newLandIceFWFlux != 0): + m = xr.DataArray(self.newLandIceFWFlux('float64'),dims=('Time','nCells')) + f['floatingBasalMassBal'] = m + + # save new mesh file + f.to_netcdf(tmp_mpasMeshFile) + f.close() #create scrip files print("Creating Scrip Files") @@ -147,8 +174,6 @@ def remap_mpas_to_mali(self) "--dst_regional", "--src_regional", '--ignore_unmapped'] check_call(args_esmf) - #<< TO DO >>: Create mask indentifying valid overlapping ocean cells - subprocess.run(["ncpdq", "-O", "-a", "Time,nVertLevels,nCells", tmp_mpasMeshFile, tmp_mpasMeshFile]) subprocess.run(["ncpdq", "-O", "-a", "Time,nVertLevelsP1,nCells", tmp_mpasMeshFile, tmp_mpasMeshFile]) subprocess.run(["ncpdq", "-O", "-a", "Time,nVertInterfaces,nCells", tmp_mpasMeshFile, tmp_mpasMeshFile]) @@ -156,10 +181,34 @@ def remap_mpas_to_mali(self) # remap the input data args_remap = ["ncremap", "-i", tmp_mpasMeshFile, - "-o", self.options.outputFile, + "-o", "tmp_" + self.options.outputFile, "-m", "tmp_mapping_file.nc"] check_call(args_remap) + # Create mask indentifying valid overlapping ocean cells. Combine all MALI terms in one input file + interpDS = xr.open_dataset("tmp_" + self.options.outputFile, decode_times=False, decode_cf=False) + interpTF = interpDS['ismip6shelfMelt_3dThermalForcing'][:,:,:] + if (self.newLandIceFWFlux != 0): + interpFWF = interpDS['floatingBasalMassBal'][:,:] + + validOceanMask = np.zeros(self.interpTF.data.shape, np.dtype = 'int') + ind = np.where(interpTF != 0) + validOceanMask[ind] = 1 + + mask = xr.dataArray(validOceanMask,dims=("Time","nCells"),long_name="Mask of MALI cells overlapping interpolated MPAS-Oceans cells") + IM = open_dataset(self.options.maliFile,decode_times=False,decode_cf=False) + IM['validOceanMask'] = mask + if (self.newLandIceFWFlux != 0): + IM['floatingBasalMassBal'] = interpFWF + IM['ismip6shelfMelt_3dThermalForcing'] = interpTF + + # define zmid coordinate in MPAS-O and translate to ISMIP6 vertical coordinate + zmid = compute_zmid(self.bottomDepth_dataArray, self.maxLevelCell_dataArray, self.layerThickness_dataArray) + IM['ismip6shelfMelt_zOcean'] = zmid + + IM.to_netcdf(self.options.outputFile) + IM.close() + #remove temporary files os.remove("tmp_*.nc") @@ -172,10 +221,6 @@ def main(): #calculate thermal forcing run.calc_ocean_thermal_forcing(self) - # << TO DO >>: copy landIceFreshwaterFlux to mpas mesh file copy for interpolation. Make name relevant to MALI - - # << TO DO >>: Vertically interpolate between MPAS-O to MALI grids - #remap mpas to mali run.remap_mpas_to_mali(self) From 5d7651215d0f6966431f7e3dedd5703dd6722351 Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Fri, 17 May 2024 15:57:42 -0700 Subject: [PATCH 04/39] Debugging through time averaging step --- .../interpolate_mpasOcean_to_mali.py | 125 ++++++++++++------ 1 file changed, 87 insertions(+), 38 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index 2ce84bc91..a9dcb7f76 100644 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -14,6 +14,7 @@ class mpasToMaliInterp: def __init__(self): + print("Gathering Information ... ") parser = OptionParser() parser.add_option("--oceanMesh", dest="mpasMeshFile", help="MPAS-Ocean mesh to be interpolated from", metavar="FILENAME") parser.add_option("--oceanDiags", dest="mpasDiagsDir", help="Directory path where MPAS-Ocean diagnostic files are stored. Should be only netcdf files in that directory", metavar="FILENAME") @@ -30,22 +31,23 @@ def __init__(self): #open mpas ocean mesh OM = xr.open_dataset(self.options.mpasMeshFile, decode_times=False, decode_cf=False) #open MALI mesh - IM = xr.open_dataset(self.options.maliMeshFile, decode_times=False, decode_cf=False) + IM = xr.open_dataset(self.options.maliFile, decode_times=False, decode_cf=False) # variables for time averaging - self.temperature = DS['timeMonthly_avg_activeTracers_temperature'][:,:,:].data - self.salinity = DS['timeMonthly_avg_activeTracers_temperature'][:,:,:].data - self.density = DS['timeMonthly_avg_density'][:,:,:].data - self.daysSinceStart = DS['timeMonthly_avg_daysSinceStartOfSim'][:].data + self.temperature = DS['timeMonthly_avg_activeTracers_temperature'][:,:,:].compute() + self.salinity = DS['timeMonthly_avg_activeTracers_temperature'][:,:,:].compute() + self.density = DS['timeMonthly_avg_density'][:,:,:].compute() + self.atmPressure = DS['timeMonthly_avg_atmosphericPressure'][:,:].compute() + self.daysSinceStart = DS['timeMonthly_avg_daysSinceStartOfSim'][:].compute() self.stTime = DS.attrs['config_start_time'] - self.landIceFloatingMask = OM['landIceFloatingMask'][:,:,:].data + self.landIceFloatingMask = OM['landIceFloatingMask'][0,:].data #not letting floating ice mask evolve for now because it's only in the mesh file self.layerThickness = DS['timeMonthly_avg_layerThickness'][:,:].data - xt = DS['xtime'][:,:].data + xt = DS['xtime_startMonthly'] xtime = np.array([xt],dtype=('S',64)) # << NOTE >>: may need to eventually use: "xtime.data.tobytes().decode()" but not sure yet try: self.landIceFreshwaterFlux = DS['timeMonthly_avg_landIceFreshwaterFlux'][:,:].data - finally: + except KeyError: print("No landIceFreshwaterFlux variable") self.landIceFreshwaterFlux = 0 @@ -57,61 +59,104 @@ def __init__(self): self.coeff_S_openOcean = DS.attrs['config_open_ocean_freezing_temperature_coeff_S'] self.coeff_p_openOcean = DS.attrs['config_open_ocean_freezing_temperature_coeff_p'] self.coeff_pS_openOcean = DS.attrs['config_open_ocean_freezing_temperature_coeff_pS'] - self.coeff_mushy_az1_liq = DS.attrs['config_open_ocean_freezing_temperature_coeff_mushy_az1_liq'] - self.coeff_mushy_openOcean = 1/coeff_mushy_az1_liq + az1_liq = DS.attrs['config_open_ocean_freezing_temperature_coeff_mushy_az1_liq'] + self.coeff_mushy_openOcean = 1/az1_liq self.coeff_0_cavity = DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_0'] self.coeff_S_cavity = DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_S'] - self.coeff_p_cavity = DS.attrs['config_land_ice_cavity_temperature_coeff_p'] + self.coeff_p_cavity = DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_p'] self.coeff_pS_cavity = DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_pS'] - self.coeff_mushy_az1_liq = DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_mushy_az1_liq'] - self.coeff_mushy_cavity = 1/coeff_mushy_az1_liq + self.coeff_mushy_cavity = 0 #variables for interpolation << NOTE >>: a couple of these fields are redundant but we want them in xarray form for compute_zmid self.bottomDepth_dataArray = OM['bottomDepth'][:] self.maxLevelCell_dataArray = OM['maxLevelCell'][:] self.layerThickness_dataArray = DS['timeMonthly_avg_layerThickness'] - def time_average_output(self) + def time_average_output(self): + print("Time Averaging ...") if (self.options.yearlyAvg == 'true'): + + print("daysSinceStart: {}".format(self.daysSinceStart)) yearsSinceStart = self.daysSinceStart / 365.0 + print("yearsSinceStart: {}".format(yearsSinceStart)) finalYear = np.floor(np.max(yearsSinceStart)) + startYear = np.floor(np.min(yearsSinceStart)) + print("final year: {}".format(finalYear)) + print("start year: {}".format(startYear)) timeStride = 1 # 1 year average - - years = np.arange(0, finalYear, timeStride) - for i in range(len(years)-1): - ind = np.where(yearsSinceStart >= years(i) and yearsSinceStart < years(i+1)) - self.newTemp[i,:,:] = np.mean(self.temperature[ind,:,:], axis=0) - self.newSal[i,:,:] = np.mean(self.salinity[ind,:,:], axis=0) - self.newDens[i,:,:] = np.mean(self.density[ind,:,:], axis=0) - self.newLThick[i,:,:] = np.mean(self.layerThickness[ind,:,:], axis=0) - self.newAtmPr[i,:] = np.mean(self.atmPressure[ind,:], axis=0) - self.newFloatIceMask[i,:] = st.mode(self.landIceFloatingMask[ind,:], axis=0) + + if (startYear != finalYear): + years = np.arange(startYear, finalYear, timeStride, dtype=int) + nt = len(years) + else : + years = startYear + nt = 1 + + print("years {}: ".format(years)) + + #pre-allocate + _,nc,nz = self.temperature.shape + self.newTemp = np.zeros((nt,nc,nz)) + self.newSal = np.zeros((nt,nc,nz)) + self.newDens = np.zeros((nt,nc,nz)) + self.newLThick= np.zeros((nt,nc,nz)) + self.newAtmPr = np.zeros((nt,nc)) + if (self.landIceFreshwaterFlux != 0): + self.newLandIceFWFlux = np.zeros((nt,nc)) + + print("starting loop ...") + if (years.ndim == 0): + print("Years: {}".format(years)) + log = np.logical_and(yearsSinceStart >= years, yearsSinceStart < years + timeStride) + self.newTemp[0,:,:] = np.mean(self.temperature[log,:,:], axis=0) + self.newSal[0,:,:] = np.mean(self.salinity[log,:,:], axis=0) + self.newDens[0,:,:] = np.mean(self.density[log,:,:], axis=0) + self.newLThick[0,:,:] = np.mean(self.layerThickness[log,:,:], axis=0) + self.newAtmPr[0,:] = np.mean(self.atmPressure[log,:], axis=0) if (self.landIceFreshwaterFlux != 0): - self.newLandIceFWFlux[i,:] = np.mean(self.landIceFreshwaterFlux[ind,:], axis=0) - elif + self.newLandIceFWFlux[0,:] = np.mean(self.landIceFreshwaterFlux[log,:], axis=0) + else : self.newLandIceFWFlux = 0 + else : + ct = 0 + for i in years: + print("Year: {}".format(years[i])) + log = np.logical_and(yearsSinceStart >= years[i], yearsSinceStart < years[i] + timeStride) + self.newTemp[ct,:,:] = np.mean(self.temperature[log,:,:], axis=0) + self.newSal[ct,:,:] = np.mean(self.salinity[log,:,:], axis=0) + self.newDens[ct,:,:] = np.mean(self.density[log,:,:], axis=0) + self.newLThick[ct,:,:] = np.mean(self.layerThickness[log,:,:], axis=0) + self.newAtmPr[ct,:] = np.mean(self.atmPressure[log,:], axis=0) + if (self.landIceFreshwaterFlux != 0): + self.newLandIceFWFlux[ct,:] = np.mean(self.landIceFreshwaterFlux[log,:], axis=0) + else : + self.newLandIceFWFlux = 0 + ct = ct + 1 # Build new xtime string. Defined on first day of each year stTime = np.datetime64(self.stTime[0:4]) stYear = np.datetime_as_string(stTime,unit='Y').astype('float64') yearVec = yearsSinceStart + stYear self.newXtime = np.array([f"{int(year)}-01-01_00:00:00" for year in yearVec],dtype=np.str_) - else #do nothing if not time averaging + else : #do nothing if not time averaging self.newTemp = self.temperature self.newSal = self.salinity self.newDens= self.density self.newLThick = self.layerThickness self.newAtmPr = self.atmPressure - self.newFloatIceMask = self.landIceFloatingMask self.newLandIceFWFlux = self.landIceFreshwaterFlux self.newXtime = xtime - def calc_ocean_thermal_forcing(self) - + def calc_ocean_thermal_forcing(self): + print("Calculating thermal forcing ... ") gravity = 9.81 - nt,nc,nz = self.newTemp.shape + #pre-allocate + nt,nc,nz = self.newTemp.shape ocn_freezing_temperature = np.zeros((nt,nc,nz)) + self.newPr = np.zeros((nt,nc,nz)) + + print("floating ice mask: {}".format(self.landIceFloatingMask.shape)) for iTime in range(nt): #calculate pressure: @@ -124,18 +169,20 @@ def calc_ocean_thermal_forcing(self) self.newPr[iTime,k,iCell] = self.newPr[iTime,iCell,k-1] + 0.5*gravity*(self.newDens[iTime,iCell,k-1]*self.newLThick[iTime,iCell,k-1] \ + self.newDens[iTime,iCell,k]*self.newLThick[iTime,iCell,k]) - if (self.newFloatingIceMask[iTime,iCell] == 1): + print("floating ice mask iCell: {}".format(self.landIceFloatingMask[iCell])) + if (self.landIceFloatingMask[iCell] == 1): ocn_freezing_temperature[iTime,iCell,:] = self.coeff_0_openOcean + self.coeff_S_openOcean * self.newSal[iTime,iCell,:] + self.coeff_p_openOcean * self.newPr[iTime,iCell,:] \ + self.coeff_pS_openOcean * self.newPr[iTime,iCell,:] * self.newSal[iTime,iCell,:] + self.coeff_mushy_openOcean * self.newSal[iTime,iCell,:] / (1.0 - self.newSal[iTime,iCell,:] / 1e3) - elif (self.newFloatingIceMask[iTime,iCell] == 0): + elif (self.landIceFloatingMask[iCell] == 0): ocn_freezing_temperature[iTime,iCell,:] = self.coeff_0_cavity + self.coeff_S_cavity * self.newSal[iTime,iCell,:] + self.coeff_p_cavity * self.newPr[iTime,iCell,:] \ + self.coeff_pS_cavity * self.newPr[iTime,iCell,:] * self.newSal[iTime,iCell,:] + self.coeff_mushy_cavity * self.newSal[iTime,iCell,:] / (1.0 - self.newSal[iTime,iCell,:] / 1e3) # Calculate thermal forcing self.oceanThermalForcing = self.newTemp - ocn_freezing_temperature - def remap_mpas_to_mali(self) + def remap_mpas_to_mali(self): + print("Start remapping ... ") # make copy of mpasMeshFile to save interpolated variable to tmp_mpasFile = "tmp_" + self.options.mpasMeshFile @@ -178,6 +225,7 @@ def remap_mpas_to_mali(self) subprocess.run(["ncpdq", "-O", "-a", "Time,nVertLevelsP1,nCells", tmp_mpasMeshFile, tmp_mpasMeshFile]) subprocess.run(["ncpdq", "-O", "-a", "Time,nVertInterfaces,nCells", tmp_mpasMeshFile, tmp_mpasMeshFile]) + print("ncremap ...") # remap the input data args_remap = ["ncremap", "-i", tmp_mpasMeshFile, @@ -185,13 +233,14 @@ def remap_mpas_to_mali(self) "-m", "tmp_mapping_file.nc"] check_call(args_remap) + print("Finishing processing and saving ...") # Create mask indentifying valid overlapping ocean cells. Combine all MALI terms in one input file interpDS = xr.open_dataset("tmp_" + self.options.outputFile, decode_times=False, decode_cf=False) interpTF = interpDS['ismip6shelfMelt_3dThermalForcing'][:,:,:] if (self.newLandIceFWFlux != 0): interpFWF = interpDS['floatingBasalMassBal'][:,:] - validOceanMask = np.zeros(self.interpTF.data.shape, np.dtype = 'int') + validOceanMask = np.zeros(self.interpTF.data.shape, dtype='int') ind = np.where(interpTF != 0) validOceanMask[ind] = 1 @@ -216,13 +265,13 @@ def main(): run = mpasToMaliInterp() #compute yearly time average if necessary - run.time_average_output(self) + run.time_average_output() #calculate thermal forcing - run.calc_ocean_thermal_forcing(self) + run.calc_ocean_thermal_forcing() #remap mpas to mali - run.remap_mpas_to_mali(self) + run.remap_mpas_to_mali() if __name__ == "__main__": main() From 7a3722739851227d7db80ffb54642e441dd8c9a9 Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Fri, 17 May 2024 16:28:02 -0700 Subject: [PATCH 05/39] Debug through calc ocean thermal forcing --- landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index a9dcb7f76..1c6f7cb5a 100644 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -155,7 +155,7 @@ def calc_ocean_thermal_forcing(self): nt,nc,nz = self.newTemp.shape ocn_freezing_temperature = np.zeros((nt,nc,nz)) self.newPr = np.zeros((nt,nc,nz)) - + print("NewPr: {}".format(self.newPr.shape)) print("floating ice mask: {}".format(self.landIceFloatingMask.shape)) for iTime in range(nt): @@ -163,13 +163,15 @@ def calc_ocean_thermal_forcing(self): for iCell in range(self.nCells): kmin = self.minLevelCell[iCell] - 1 kmax = self.maxLevelCell[iCell] - 1 - self.newPr[iTime,kmin,iCell] = self.newAtmPr[iTime,iCell] + self.newDens[iTime,iCell,kmin]*gravity*0.5*self.newLThick[iTime,iCell,kmin] + print("kmin: {}".format(kmin)) + print("iCell: {}".format(iCell)) + + self.newPr[iTime,iCell,kmin] = self.newAtmPr[iTime,iCell] + self.newDens[iTime,iCell,kmin]*gravity*0.5*self.newLThick[iTime,iCell,kmin] for k in np.arange(kmin + 1, kmax): - self.newPr[iTime,k,iCell] = self.newPr[iTime,iCell,k-1] + 0.5*gravity*(self.newDens[iTime,iCell,k-1]*self.newLThick[iTime,iCell,k-1] \ + self.newPr[iTime,iCell,k] = self.newPr[iTime,iCell,k-1] + 0.5*gravity*(self.newDens[iTime,iCell,k-1]*self.newLThick[iTime,iCell,k-1] \ + self.newDens[iTime,iCell,k]*self.newLThick[iTime,iCell,k]) - print("floating ice mask iCell: {}".format(self.landIceFloatingMask[iCell])) if (self.landIceFloatingMask[iCell] == 1): ocn_freezing_temperature[iTime,iCell,:] = self.coeff_0_openOcean + self.coeff_S_openOcean * self.newSal[iTime,iCell,:] + self.coeff_p_openOcean * self.newPr[iTime,iCell,:] \ + self.coeff_pS_openOcean * self.newPr[iTime,iCell,:] * self.newSal[iTime,iCell,:] + self.coeff_mushy_openOcean * self.newSal[iTime,iCell,:] / (1.0 - self.newSal[iTime,iCell,:] / 1e3) From 0581fa5d954730a7a8b987f2be42c009821cbb08 Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Thu, 23 May 2024 13:24:48 -0700 Subject: [PATCH 06/39] Debugging until remapping --- .../interpolate_mpasOcean_to_mali.py | 36 ++++++++++--------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index 1c6f7cb5a..7bf986308 100644 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -163,22 +163,26 @@ def calc_ocean_thermal_forcing(self): for iCell in range(self.nCells): kmin = self.minLevelCell[iCell] - 1 kmax = self.maxLevelCell[iCell] - 1 - print("kmin: {}".format(kmin)) - print("iCell: {}".format(iCell)) - self.newPr[iTime,iCell,kmin] = self.newAtmPr[iTime,iCell] + self.newDens[iTime,iCell,kmin]*gravity*0.5*self.newLThick[iTime,iCell,kmin] + self.newPr[iTime,iCell,kmin] = self.newAtmPr[iTime,iCell] + \ + self.newDens[iTime,iCell,kmin]*gravity*0.5*self.newLThick[iTime,iCell,kmin] for k in np.arange(kmin + 1, kmax): - self.newPr[iTime,iCell,k] = self.newPr[iTime,iCell,k-1] + 0.5*gravity*(self.newDens[iTime,iCell,k-1]*self.newLThick[iTime,iCell,k-1] \ + self.newPr[iTime,iCell,k] = self.newPr[iTime,iCell,k-1] + \ + 0.5*gravity*(self.newDens[iTime,iCell,k-1]*self.newLThick[iTime,iCell,k-1] \ + self.newDens[iTime,iCell,k]*self.newLThick[iTime,iCell,k]) if (self.landIceFloatingMask[iCell] == 1): - ocn_freezing_temperature[iTime,iCell,:] = self.coeff_0_openOcean + self.coeff_S_openOcean * self.newSal[iTime,iCell,:] + self.coeff_p_openOcean * self.newPr[iTime,iCell,:] \ - + self.coeff_pS_openOcean * self.newPr[iTime,iCell,:] * self.newSal[iTime,iCell,:] + self.coeff_mushy_openOcean * self.newSal[iTime,iCell,:] / (1.0 - self.newSal[iTime,iCell,:] / 1e3) + ocn_freezing_temperature[iTime,iCell,:] = self.coeff_0_openOcean + \ + self.coeff_S_openOcean * self.newSal[iTime,iCell,:] + self.coeff_p_openOcean * self.newPr[iTime,iCell,:] \ + + self.coeff_pS_openOcean * self.newPr[iTime,iCell,:] * self.newSal[iTime,iCell,:] \ + + self.coeff_mushy_openOcean * self.newSal[iTime,iCell,:] / (1.0 - self.newSal[iTime,iCell,:] / 1e3) elif (self.landIceFloatingMask[iCell] == 0): - ocn_freezing_temperature[iTime,iCell,:] = self.coeff_0_cavity + self.coeff_S_cavity * self.newSal[iTime,iCell,:] + self.coeff_p_cavity * self.newPr[iTime,iCell,:] \ - + self.coeff_pS_cavity * self.newPr[iTime,iCell,:] * self.newSal[iTime,iCell,:] + self.coeff_mushy_cavity * self.newSal[iTime,iCell,:] / (1.0 - self.newSal[iTime,iCell,:] / 1e3) + ocn_freezing_temperature[iTime,iCell,:] = self.coeff_0_cavity + \ + self.coeff_S_cavity * self.newSal[iTime,iCell,:] + self.coeff_p_cavity * self.newPr[iTime,iCell,:] \ + + self.coeff_pS_cavity * self.newPr[iTime,iCell,:] * self.newSal[iTime,iCell,:] \ + + self.coeff_mushy_cavity * self.newSal[iTime,iCell,:] / (1.0 - self.newSal[iTime,iCell,:] / 1e3) # Calculate thermal forcing self.oceanThermalForcing = self.newTemp - ocn_freezing_temperature @@ -187,7 +191,7 @@ def remap_mpas_to_mali(self): print("Start remapping ... ") # make copy of mpasMeshFile to save interpolated variable to - tmp_mpasFile = "tmp_" + self.options.mpasMeshFile + tmp_mpasMeshFile = self.options.mpasMeshFile + ".tmp" shutil.copy(self.options.mpasMeshFile, tmp_mpasMeshFile) # populate tmp_mpasMeshFile with variables to be interpolated @@ -205,8 +209,8 @@ def remap_mpas_to_mali(self): #create scrip files print("Creating Scrip Files") - mali_scripfile = "tmp_mali_scrip.nc" - mpas_scripfile = "tmp_mpas_scrip.nc" + mali_scripfile = "mali_scrip.nc.tmp" + mpas_scripfile = "mpas_scrip.tmp" scrip_from_mpas(self.options.maliFile, mali_scripfile) scrip_from_mpas(tmp_mpasMeshFile, mpas_scripfile) @@ -217,7 +221,7 @@ def remap_mpas_to_mali(self): '-n', self.options.ntasks, 'ESMF_RegridWeightGen', '--source', mpas_scripfile, '--destination', mali_scripfile, - '--weight', "tmp_mapping_file.nc", + '--weight', "mapping_file.nc.tmp", '--method', self.options.method, '-i', '-64bit_offset', "--dst_regional", "--src_regional", '--ignore_unmapped'] @@ -231,13 +235,13 @@ def remap_mpas_to_mali(self): # remap the input data args_remap = ["ncremap", "-i", tmp_mpasMeshFile, - "-o", "tmp_" + self.options.outputFile, - "-m", "tmp_mapping_file.nc"] + "-o", self.options.outputFile + ".tmp", + "-m", "mapping_file.nc.tmp"] check_call(args_remap) print("Finishing processing and saving ...") # Create mask indentifying valid overlapping ocean cells. Combine all MALI terms in one input file - interpDS = xr.open_dataset("tmp_" + self.options.outputFile, decode_times=False, decode_cf=False) + interpDS = xr.open_dataset(self.options.outputFile + ".tmp", decode_times=False, decode_cf=False) interpTF = interpDS['ismip6shelfMelt_3dThermalForcing'][:,:,:] if (self.newLandIceFWFlux != 0): interpFWF = interpDS['floatingBasalMassBal'][:,:] @@ -261,7 +265,7 @@ def remap_mpas_to_mali(self): IM.close() #remove temporary files - os.remove("tmp_*.nc") + os.remove("*.tmp") def main(): run = mpasToMaliInterp() From 3a2339b680ed4c029115df5e9dad0e6702327953 Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Fri, 14 Jun 2024 08:34:21 -0700 Subject: [PATCH 07/39] Add vertical interpolation --- .../interpolate_mpasOcean_to_mali.py | 260 +++++++++++++++--- 1 file changed, 215 insertions(+), 45 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index 7bf986308..fc4b20318 100644 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -10,6 +10,7 @@ from mpas_tools.scrip.from_mpas import scrip_from_mpas from mpas_tools.ocean.depth import compute_zmid from optparse import OptionParser +import time class mpasToMaliInterp: @@ -41,7 +42,12 @@ def __init__(self): self.daysSinceStart = DS['timeMonthly_avg_daysSinceStartOfSim'][:].compute() self.stTime = DS.attrs['config_start_time'] self.landIceFloatingMask = OM['landIceFloatingMask'][0,:].data #not letting floating ice mask evolve for now because it's only in the mesh file - self.layerThickness = DS['timeMonthly_avg_layerThickness'][:,:].data + avg_layerThickness = DS['timeMonthly_avg_layerThickness'] + self.layerThickness = avg_layerThickness.data + alT = avg_layerThickness.compute().values + print("avg_layerThickness: {}".format(alT.shape)) + print("avg_layerTHickness NaNs: {}".format(np.sum(np.isnan(alT)))) + xt = DS['xtime_startMonthly'] xtime = np.array([xt],dtype=('S',64)) # << NOTE >>: may need to eventually use: "xtime.data.tobytes().decode()" but not sure yet @@ -51,9 +57,24 @@ def __init__(self): print("No landIceFreshwaterFlux variable") self.landIceFreshwaterFlux = 0 - # additional variables for computing thermal forcing + #check for nans in original data + if np.isnan(self.temperature).any(): + print("Original temperature contains NaNs") + else: + print("Original temperature does not contain NaNs") + + # additional variables for computing thermal forcing + bottomDepth = OM['bottomDepth'] + bD = bottomDepth.data + print("bottomDepth: {}".format(bD.shape)) + print("bottomDepth NaNs: {}".format(np.sum(np.isnan(bD)))) + self.minLevelCell = OM['minLevelCell'][:].data - self.maxLevelCell = OM['maxLevelCell'][:].data + maxLevelCell = OM['maxLevelCell'] + self.maxLevelCell = maxLevelCell.data + print("maxLevelCell: {}".format(self.maxLevelCell.shape)) + print("maxLevelCell NaNs: {}".format(np.sum(np.isnan(self.maxLevelCell)))) + self.nCells = OM.sizes['nCells'] self.coeff_0_openOcean = DS.attrs['config_open_ocean_freezing_temperature_coeff_0'] self.coeff_S_openOcean = DS.attrs['config_open_ocean_freezing_temperature_coeff_S'] @@ -67,22 +88,20 @@ def __init__(self): self.coeff_pS_cavity = DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_pS'] self.coeff_mushy_cavity = 0 - #variables for interpolation << NOTE >>: a couple of these fields are redundant but we want them in xarray form for compute_zmid - self.bottomDepth_dataArray = OM['bottomDepth'][:] - self.maxLevelCell_dataArray = OM['maxLevelCell'][:] - self.layerThickness_dataArray = DS['timeMonthly_avg_layerThickness'] - + # Define vertical coordinates of mpas output + mpas_cellCenterElev = compute_zmid(bottomDepth, maxLevelCell, avg_layerThickness) + self.mpas_cellCenterElev = mpas_cellCenterElev.data + print("mpas_cellCenterElev {}".format(self.mpas_cellCenterElev.shape)) + mpas_cellCenterElev = mpas_cellCenterElev.compute().values + print("mpas_cellCenterElev number of NaNs: {}".format(np.sum(np.isnan(mpas_cellCenterElev)))) + def time_average_output(self): print("Time Averaging ...") if (self.options.yearlyAvg == 'true'): - print("daysSinceStart: {}".format(self.daysSinceStart)) yearsSinceStart = self.daysSinceStart / 365.0 - print("yearsSinceStart: {}".format(yearsSinceStart)) finalYear = np.floor(np.max(yearsSinceStart)) startYear = np.floor(np.min(yearsSinceStart)) - print("final year: {}".format(finalYear)) - print("start year: {}".format(startYear)) timeStride = 1 # 1 year average if (startYear != finalYear): @@ -92,8 +111,6 @@ def time_average_output(self): years = startYear nt = 1 - print("years {}: ".format(years)) - #pre-allocate _,nc,nz = self.temperature.shape self.newTemp = np.zeros((nt,nc,nz)) @@ -101,17 +118,19 @@ def time_average_output(self): self.newDens = np.zeros((nt,nc,nz)) self.newLThick= np.zeros((nt,nc,nz)) self.newAtmPr = np.zeros((nt,nc)) + self.newMpasCCE = np.zeros((nt,nc,nz)) if (self.landIceFreshwaterFlux != 0): self.newLandIceFWFlux = np.zeros((nt,nc)) print("starting loop ...") + st = time.time() if (years.ndim == 0): - print("Years: {}".format(years)) log = np.logical_and(yearsSinceStart >= years, yearsSinceStart < years + timeStride) self.newTemp[0,:,:] = np.mean(self.temperature[log,:,:], axis=0) self.newSal[0,:,:] = np.mean(self.salinity[log,:,:], axis=0) self.newDens[0,:,:] = np.mean(self.density[log,:,:], axis=0) self.newLThick[0,:,:] = np.mean(self.layerThickness[log,:,:], axis=0) + self.newMpasCCE[0,:,:] = np.mean(self.mpas_cellCenterElev[log,:,:], axis=0) self.newAtmPr[0,:] = np.mean(self.atmPressure[log,:], axis=0) if (self.landIceFreshwaterFlux != 0): self.newLandIceFWFlux[0,:] = np.mean(self.landIceFreshwaterFlux[log,:], axis=0) @@ -120,18 +139,21 @@ def time_average_output(self): else : ct = 0 for i in years: - print("Year: {}".format(years[i])) log = np.logical_and(yearsSinceStart >= years[i], yearsSinceStart < years[i] + timeStride) self.newTemp[ct,:,:] = np.mean(self.temperature[log,:,:], axis=0) self.newSal[ct,:,:] = np.mean(self.salinity[log,:,:], axis=0) self.newDens[ct,:,:] = np.mean(self.density[log,:,:], axis=0) self.newLThick[ct,:,:] = np.mean(self.layerThickness[log,:,:], axis=0) + self.newMpasCCE[ct,:,:] = np.mean(self.mpas_cellCenterElev[log,:,:], axis=0) self.newAtmPr[ct,:] = np.mean(self.atmPressure[log,:], axis=0) if (self.landIceFreshwaterFlux != 0): self.newLandIceFWFlux[ct,:] = np.mean(self.landIceFreshwaterFlux[log,:], axis=0) else : self.newLandIceFWFlux = 0 ct = ct + 1 + nd = time.time() + tm = nd - st + print("Time averaging loop", tm, "seconds") # Build new xtime string. Defined on first day of each year stTime = np.datetime64(self.stTime[0:4]) @@ -143,9 +165,34 @@ def time_average_output(self): self.newSal = self.salinity self.newDens= self.density self.newLThick = self.layerThickness + self.newMpasCCE = self.mpas_cellCenterElev self.newAtmPr = self.atmPressure self.newLandIceFWFlux = self.landIceFreshwaterFlux self.newXtime = xtime + + print("neMPASCCE: {}".format(self.newMpasCCE.shape)) + print("NaNs in newMpasCCE: {}".format(np.sum(np.isnan(self.newMpasCCE)))) + + #check for nans + if np.isnan(self.mpas_cellCenterElev).any(): + print("mpas_cellCenterElev contains NaNs") + else: + print("mpas_cellCenterElev does not contain NaNs") + + if np.isnan(self.newMpasCCE).any(): + print("newMpasCCE contains NaNs") + else: + print("newMpasCCE does not contain NaNs") + + if np.isnan(self.temperature).any(): + print("temperature contains NaNs") + else: + print("temperature does not contain NaNs") + + if np.isnan(self.newTemp).any(): + print("newTemp contains NaNs") + else: + print("newTemp does not contain NaNs") def calc_ocean_thermal_forcing(self): print("Calculating thermal forcing ... ") @@ -155,8 +202,7 @@ def calc_ocean_thermal_forcing(self): nt,nc,nz = self.newTemp.shape ocn_freezing_temperature = np.zeros((nt,nc,nz)) self.newPr = np.zeros((nt,nc,nz)) - print("NewPr: {}".format(self.newPr.shape)) - print("floating ice mask: {}".format(self.landIceFloatingMask.shape)) + st = time.time() for iTime in range(nt): #calculate pressure: @@ -183,34 +229,68 @@ def calc_ocean_thermal_forcing(self): self.coeff_S_cavity * self.newSal[iTime,iCell,:] + self.coeff_p_cavity * self.newPr[iTime,iCell,:] \ + self.coeff_pS_cavity * self.newPr[iTime,iCell,:] * self.newSal[iTime,iCell,:] \ + self.coeff_mushy_cavity * self.newSal[iTime,iCell,:] / (1.0 - self.newSal[iTime,iCell,:] / 1e3) - + nd = time.time() + tm = nd - st + print("Ocean thermal forcing loop", tm, "seconds") # Calculate thermal forcing self.oceanThermalForcing = self.newTemp - ocn_freezing_temperature + if np.isnan(self.oceanThermalForcing).any(): + print("oceanThermalForcing contains NaNs") + else: + print("oceanThermalForcing does not contain NaNs") + def remap_mpas_to_mali(self): print("Start remapping ... ") - # make copy of mpasMeshFile to save interpolated variable to - tmp_mpasMeshFile = self.options.mpasMeshFile + ".tmp" - shutil.copy(self.options.mpasMeshFile, tmp_mpasMeshFile) - # populate tmp_mpasMeshFile with variables to be interpolated - f = xr.open_dataset(tmp_mpasMeshFile, decode_times=False, decode_cf=False) + f = xr.open_dataset(self.options.mpasMeshFile, decode_times=False, decode_cf=False) tf = xr.DataArray(self.oceanThermalForcing.astype('float64'),dims=('Time','nCells','nVertLevels')) + f['ismip6shelfMelt_3dThermalForcing'] = tf + mcce = xr.DataArray(self.newMpasCCE.astype('float64'),dims=('Time','nCells','nVertLevels')) + f['mpas_cellCenterElev'] = mcce + if (self.newLandIceFWFlux != 0): - m = xr.DataArray(self.newLandIceFWFlux('float64'),dims=('Time','nCells')) + m = xr.DataArray(self.newLandIceFWFlux.astype('float64'),dims=('Time','nCells')) f['floatingBasalMassBal'] = m + + #delete unncessary variables with 'string1' dimension or will cause errors. + try: + f = f.drop_vars('simulationStartTime') + except ValueError: + f = f + try: + f = f.drop_vars('forcingGroupNames') + except ValueError: + f = f + + try: + f = f.drop_vars('forcingGroupRestartTimes') + except ValueError: + f = f + + st = time.time() # save new mesh file + tmp_mpasMeshFile = "tmp_mpasMeshFile.nc" f.to_netcdf(tmp_mpasMeshFile) f.close() + nd = time.time() + tm = nd - st + print("saving updates mpas mesh file", tm, "seconds") + + st = time.time() + subprocess.run(["ncatted", "-a", "_FillValue,,d,,", tmp_mpasMeshFile]) + nd = time.time() + tm = nd - st + print("Removing fill value:", tm, "seconds") #create scrip files print("Creating Scrip Files") - mali_scripfile = "mali_scrip.nc.tmp" - mpas_scripfile = "mpas_scrip.tmp" + mali_scripfile = "tmp_mali_scrip.nc" + mpas_scripfile = "tmp_mpas_scrip.nc" scrip_from_mpas(self.options.maliFile, mali_scripfile) scrip_from_mpas(tmp_mpasMeshFile, mpas_scripfile) @@ -221,7 +301,7 @@ def remap_mpas_to_mali(self): '-n', self.options.ntasks, 'ESMF_RegridWeightGen', '--source', mpas_scripfile, '--destination', mali_scripfile, - '--weight', "mapping_file.nc.tmp", + '--weight', "tmp_mapping_file.nc", '--method', self.options.method, '-i', '-64bit_offset', "--dst_regional", "--src_regional", '--ignore_unmapped'] @@ -235,37 +315,126 @@ def remap_mpas_to_mali(self): # remap the input data args_remap = ["ncremap", "-i", tmp_mpasMeshFile, - "-o", self.options.outputFile + ".tmp", - "-m", "mapping_file.nc.tmp"] + "-o", "tmp_outputFile.nc", + "-m", "tmp_mapping_file.nc"] check_call(args_remap) - print("Finishing processing and saving ...") + #Vertical interpolation + print("Vertically interpolating onto mali grid") + interpDS = xr.open_dataset("tmp_outputFile.nc", decode_times=False, decode_cf=False) + TF = interpDS['ismip6shelfMelt_3dThermalForcing'][:,:,:].data + mpas_cellCenterElev = interpDS['mpas_cellCenterElev'][:,:,:].data + + #Define vertical coordinates for mali grid + IM = xr.open_dataset(self.options.maliFile,decode_times=False,decode_cf=False) + layerThicknessFractions = IM['layerThicknessFractions'][:].data + bedTopography = IM['bedTopography'][:,:].data + thickness = IM['thickness'][:,:].data + + nz, = layerThicknessFractions.shape + nt,nc = thickness.shape + + layerThicknessFractions = layerThicknessFractions.reshape((1,1,nz)) + thickness = thickness.reshape((nt,nc,1)) + bedTopography = bedTopography.reshape((nt,nc,1)) + + layerThickness = thickness * layerThicknessFractions + mali_cellCenterElev = np.cumsum(layerThickness,axis=2) - layerThickness/2 + bedTopography + + # Reshape to (nt*nc, nz) before interpolation to avoid looping + mali_cellCenterElev = mali_cellCenterElev.reshape(nt*nc, -1) + TF = np.transpose(TF,(0,2,1)) + TF = TF.reshape(nt*nc, -1) + mpas_cellCenterElev = np.transpose(mpas_cellCenterElev,(0,2,1)) + mpas_cellCenterElev = mpas_cellCenterElev.reshape(nt*nc, -1) + + print("Before Interpolation") + print("TF zeroth dim: {}".format(TF[0:5,0])) + print("TF first dim: {}".format(TF[0,0:5])) + + print("mpas_cellCenterElev zeroth dim: {}".format(mpas_cellCenterElev[0:5,0])) + print("mpas_cellCenterElev first dim: {}".format(mpas_cellCenterElev[0,0:5])) + + st = time.time() + # Linear interpolation + + if np.isnan(mali_cellCenterElev).any(): + print("mali_cellCenterElev contains NaNs") + else: + print("mali_cellCenterElev does not contain NaNs") + + if np.isnan(mpas_cellCenterElev).any(): + print("mpas_cellCenterElev contains NaNs") + else: + print("mpas_cellCenterElev does not contain NaNs") + + if np.isnan(TF).any(): + print("TF contains NaNs") + else: + print("TF does not contain NaNs") + + nan_count = 0 + for i in range(nt*nc): + ind = np.where(~np.isnan(TF[i,:].flatten() * mpas_cellCenterElev[i,:].flatten())) + + #if (ind.size != 0): #only interpolation where valid data + ct = 0 + if len(ind[0]) != 0: + if (ct == 0): + print("mali_cellCenterElev[i,:]".format(mali_cellCenterElev[i,:].flatten().shape)) + print("mpas_cellCenterElev[i,ind]".format(mpas_cellCenterElev[i,ind].flatten().shape)) + print("TF[i,ind]".format(mali_cellCenterElev[i,ind].flatten().shape)) + ct = ct + 1 + interpTF = np.array(np.interp(mali_cellCenterElev[i,:].flatten(), mpas_cellCenterElev[i,ind].flatten(), TF[i,ind].flatten())) + else: + nan_count = nan_count + 1 + print("NaN Count: {}; nt*nc = {}".format(nan_count, nt*nc)) + nd = time.time() + tm = nd-st + print("vertical interpolation time: {}",tm, "seconds") + + print("After interpolation:") + print("interpTF: {}".format(interpTF[0:5,0])) + + # Reshape interpTF back to (nt, nc, nz) + interpTF = interpTF.reshape(nt, nc, -1) + + print("After Reshape:") + print("interpTF: {}".format(interpTF[0,0:5,0])) + # Create mask indentifying valid overlapping ocean cells. Combine all MALI terms in one input file - interpDS = xr.open_dataset(self.options.outputFile + ".tmp", decode_times=False, decode_cf=False) - interpTF = interpDS['ismip6shelfMelt_3dThermalForcing'][:,:,:] + #Making this a 3-D variable for now, but may want to make 2-D eventually + validOpenOceanMask = np.zeros((nt,nc), dtype='int') + ind = np.where(interpTF[:,:,0].data != 0) + validOpenOceanMask[ind] = 1 + if (self.newLandIceFWFlux != 0): interpFWF = interpDS['floatingBasalMassBal'][:,:] + print("interpFWF: {}".format(interpFWF.data.shape)) - validOceanMask = np.zeros(self.interpTF.data.shape, dtype='int') - ind = np.where(interpTF != 0) - validOceanMask[ind] = 1 + print("Saving") + + print("interpTF: {}".format(interpTF.shape)) + print("validOpenOceanMask: {}".format(validOpenOceanMask.shape)) + + mask = xr.DataArray(validOpenOceanMask,dims=("Time","nCells"),attrs={'long_name':'Mask of MALI cells overlapping interpolated MPAS-Oceans cells'}) + interptf = xr.DataArray(interpTF, dims=("Time","nCells","nVertLevels")) - mask = xr.dataArray(validOceanMask,dims=("Time","nCells"),long_name="Mask of MALI cells overlapping interpolated MPAS-Oceans cells") - IM = open_dataset(self.options.maliFile,decode_times=False,decode_cf=False) IM['validOceanMask'] = mask if (self.newLandIceFWFlux != 0): IM['floatingBasalMassBal'] = interpFWF - IM['ismip6shelfMelt_3dThermalForcing'] = interpTF - - # define zmid coordinate in MPAS-O and translate to ISMIP6 vertical coordinate - zmid = compute_zmid(self.bottomDepth_dataArray, self.maxLevelCell_dataArray, self.layerThickness_dataArray) - IM['ismip6shelfMelt_zOcean'] = zmid + IM['ismip6shelfMelt_3dThermalForcing'] = interptf IM.to_netcdf(self.options.outputFile) IM.close() +# subprocess.run(["ncatted", "-a", "_FillValue,,d,,", self.options.outputFile]) + #remove temporary files - os.remove("*.tmp") + files = os.listdir() + for file in files: + if file.startswith("tmp"): + os.remove(file) def main(): run = mpasToMaliInterp() @@ -278,7 +447,8 @@ def main(): #remap mpas to mali run.remap_mpas_to_mali() - + + print("Finished.") if __name__ == "__main__": main() From 274513b7d0547326b63066f12e9fd2a86ced2c34 Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Tue, 18 Jun 2024 15:45:53 -0700 Subject: [PATCH 08/39] Misc. Debugging --- .../interpolate_mpasOcean_to_mali.py | 100 +++--------------- 1 file changed, 12 insertions(+), 88 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index fc4b20318..c05d9f679 100644 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -44,9 +44,6 @@ def __init__(self): self.landIceFloatingMask = OM['landIceFloatingMask'][0,:].data #not letting floating ice mask evolve for now because it's only in the mesh file avg_layerThickness = DS['timeMonthly_avg_layerThickness'] self.layerThickness = avg_layerThickness.data - alT = avg_layerThickness.compute().values - print("avg_layerThickness: {}".format(alT.shape)) - print("avg_layerTHickness NaNs: {}".format(np.sum(np.isnan(alT)))) xt = DS['xtime_startMonthly'] xtime = np.array([xt],dtype=('S',64)) @@ -57,23 +54,13 @@ def __init__(self): print("No landIceFreshwaterFlux variable") self.landIceFreshwaterFlux = 0 - #check for nans in original data - if np.isnan(self.temperature).any(): - print("Original temperature contains NaNs") - else: - print("Original temperature does not contain NaNs") - # additional variables for computing thermal forcing bottomDepth = OM['bottomDepth'] bD = bottomDepth.data - print("bottomDepth: {}".format(bD.shape)) - print("bottomDepth NaNs: {}".format(np.sum(np.isnan(bD)))) self.minLevelCell = OM['minLevelCell'][:].data maxLevelCell = OM['maxLevelCell'] self.maxLevelCell = maxLevelCell.data - print("maxLevelCell: {}".format(self.maxLevelCell.shape)) - print("maxLevelCell NaNs: {}".format(np.sum(np.isnan(self.maxLevelCell)))) self.nCells = OM.sizes['nCells'] self.coeff_0_openOcean = DS.attrs['config_open_ocean_freezing_temperature_coeff_0'] @@ -91,9 +78,6 @@ def __init__(self): # Define vertical coordinates of mpas output mpas_cellCenterElev = compute_zmid(bottomDepth, maxLevelCell, avg_layerThickness) self.mpas_cellCenterElev = mpas_cellCenterElev.data - print("mpas_cellCenterElev {}".format(self.mpas_cellCenterElev.shape)) - mpas_cellCenterElev = mpas_cellCenterElev.compute().values - print("mpas_cellCenterElev number of NaNs: {}".format(np.sum(np.isnan(mpas_cellCenterElev)))) def time_average_output(self): print("Time Averaging ...") @@ -170,30 +154,6 @@ def time_average_output(self): self.newLandIceFWFlux = self.landIceFreshwaterFlux self.newXtime = xtime - print("neMPASCCE: {}".format(self.newMpasCCE.shape)) - print("NaNs in newMpasCCE: {}".format(np.sum(np.isnan(self.newMpasCCE)))) - - #check for nans - if np.isnan(self.mpas_cellCenterElev).any(): - print("mpas_cellCenterElev contains NaNs") - else: - print("mpas_cellCenterElev does not contain NaNs") - - if np.isnan(self.newMpasCCE).any(): - print("newMpasCCE contains NaNs") - else: - print("newMpasCCE does not contain NaNs") - - if np.isnan(self.temperature).any(): - print("temperature contains NaNs") - else: - print("temperature does not contain NaNs") - - if np.isnan(self.newTemp).any(): - print("newTemp contains NaNs") - else: - print("newTemp does not contain NaNs") - def calc_ocean_thermal_forcing(self): print("Calculating thermal forcing ... ") gravity = 9.81 @@ -234,11 +194,6 @@ def calc_ocean_thermal_forcing(self): print("Ocean thermal forcing loop", tm, "seconds") # Calculate thermal forcing self.oceanThermalForcing = self.newTemp - ocn_freezing_temperature - - if np.isnan(self.oceanThermalForcing).any(): - print("oceanThermalForcing contains NaNs") - else: - print("oceanThermalForcing does not contain NaNs") def remap_mpas_to_mali(self): print("Start remapping ... ") @@ -330,7 +285,6 @@ def remap_mpas_to_mali(self): layerThicknessFractions = IM['layerThicknessFractions'][:].data bedTopography = IM['bedTopography'][:,:].data thickness = IM['thickness'][:,:].data - nz, = layerThicknessFractions.shape nt,nc = thickness.shape @@ -348,64 +302,31 @@ def remap_mpas_to_mali(self): mpas_cellCenterElev = np.transpose(mpas_cellCenterElev,(0,2,1)) mpas_cellCenterElev = mpas_cellCenterElev.reshape(nt*nc, -1) - print("Before Interpolation") - print("TF zeroth dim: {}".format(TF[0:5,0])) - print("TF first dim: {}".format(TF[0,0:5])) - - print("mpas_cellCenterElev zeroth dim: {}".format(mpas_cellCenterElev[0:5,0])) - print("mpas_cellCenterElev first dim: {}".format(mpas_cellCenterElev[0,0:5])) - - st = time.time() # Linear interpolation - - if np.isnan(mali_cellCenterElev).any(): - print("mali_cellCenterElev contains NaNs") - else: - print("mali_cellCenterElev does not contain NaNs") - - if np.isnan(mpas_cellCenterElev).any(): - print("mpas_cellCenterElev contains NaNs") - else: - print("mpas_cellCenterElev does not contain NaNs") - - if np.isnan(TF).any(): - print("TF contains NaNs") - else: - print("TF does not contain NaNs") - + st = time.time() nan_count = 0 + interpTF = np.zeros((mali_cellCenterElev.shape)) for i in range(nt*nc): ind = np.where(~np.isnan(TF[i,:].flatten() * mpas_cellCenterElev[i,:].flatten())) - #if (ind.size != 0): #only interpolation where valid data ct = 0 if len(ind[0]) != 0: - if (ct == 0): - print("mali_cellCenterElev[i,:]".format(mali_cellCenterElev[i,:].flatten().shape)) - print("mpas_cellCenterElev[i,ind]".format(mpas_cellCenterElev[i,ind].flatten().shape)) - print("TF[i,ind]".format(mali_cellCenterElev[i,ind].flatten().shape)) - ct = ct + 1 - interpTF = np.array(np.interp(mali_cellCenterElev[i,:].flatten(), mpas_cellCenterElev[i,ind].flatten(), TF[i,ind].flatten())) + interpTF[i,:] = np.array(np.interp(mali_cellCenterElev[i,:].flatten(), mpas_cellCenterElev[i,ind].flatten(), TF[i,ind].flatten())) else: nan_count = nan_count + 1 - print("NaN Count: {}; nt*nc = {}".format(nan_count, nt*nc)) nd = time.time() tm = nd-st print("vertical interpolation time: {}",tm, "seconds") - print("After interpolation:") - print("interpTF: {}".format(interpTF[0:5,0])) - # Reshape interpTF back to (nt, nc, nz) interpTF = interpTF.reshape(nt, nc, -1) - print("After Reshape:") - print("interpTF: {}".format(interpTF[0,0:5,0])) - # Create mask indentifying valid overlapping ocean cells. Combine all MALI terms in one input file #Making this a 3-D variable for now, but may want to make 2-D eventually validOpenOceanMask = np.zeros((nt,nc), dtype='int') - ind = np.where(interpTF[:,:,0].data != 0) + print("interpTF: {}".format(interpTF.data.shape)) + surfaceTF = interpTF[:,:,0] + ind = np.where(surfaceTF != 0) validOpenOceanMask[ind] = 1 if (self.newLandIceFWFlux != 0): @@ -418,17 +339,21 @@ def remap_mpas_to_mali(self): print("validOpenOceanMask: {}".format(validOpenOceanMask.shape)) mask = xr.DataArray(validOpenOceanMask,dims=("Time","nCells"),attrs={'long_name':'Mask of MALI cells overlapping interpolated MPAS-Oceans cells'}) - interptf = xr.DataArray(interpTF, dims=("Time","nCells","nVertLevels")) + + IM = IM.expand_dims(nISMIP6OceanLayers=mali_cellCenterElev) + interptf = xr.DataArray(interpTF, dims=("Time","nCells","nISMIP6OceanLayers")) + zlayers = xr.DataArray(mali_cellCenterElev, dims=("nISMIP6OceanLayers")) IM['validOceanMask'] = mask if (self.newLandIceFWFlux != 0): IM['floatingBasalMassBal'] = interpFWF IM['ismip6shelfMelt_3dThermalForcing'] = interptf + IM['ismip6shelfMelt_zOcean'] = zlayers IM.to_netcdf(self.options.outputFile) IM.close() -# subprocess.run(["ncatted", "-a", "_FillValue,,d,,", self.options.outputFile]) + subprocess.run(["ncatted", "-a", "_FillValue,,d,,", self.options.outputFile]) #remove temporary files files = os.listdir() @@ -453,4 +378,3 @@ def main(): main() - From 739a7d06aa791d5bebb146769a42d5451abf6733 Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Mon, 24 Jun 2024 15:50:27 -0700 Subject: [PATCH 09/39] Save xtime --- .../interpolate_mpasOcean_to_mali.py | 157 +++++++++++++----- 1 file changed, 112 insertions(+), 45 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index c05d9f679..c4ac27c1d 100644 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -11,6 +11,7 @@ from mpas_tools.ocean.depth import compute_zmid from optparse import OptionParser import time +import cftime class mpasToMaliInterp: @@ -45,9 +46,11 @@ def __init__(self): avg_layerThickness = DS['timeMonthly_avg_layerThickness'] self.layerThickness = avg_layerThickness.data - xt = DS['xtime_startMonthly'] - xtime = np.array([xt],dtype=('S',64)) + #xt = DS['xtime_startMonthly'] + + #xtime = np.array([xt],dtype=('S',64)) # << NOTE >>: may need to eventually use: "xtime.data.tobytes().decode()" but not sure yet + try: self.landIceFreshwaterFlux = DS['timeMonthly_avg_landIceFreshwaterFlux'][:,:].data except KeyError: @@ -83,11 +86,11 @@ def time_average_output(self): print("Time Averaging ...") if (self.options.yearlyAvg == 'true'): - yearsSinceStart = self.daysSinceStart / 365.0 + yearsSinceStart = self.daysSinceStart.data / 365.0 finalYear = np.floor(np.max(yearsSinceStart)) startYear = np.floor(np.min(yearsSinceStart)) timeStride = 1 # 1 year average - + if (startYear != finalYear): years = np.arange(startYear, finalYear, timeStride, dtype=int) nt = len(years) @@ -105,8 +108,15 @@ def time_average_output(self): self.newMpasCCE = np.zeros((nt,nc,nz)) if (self.landIceFreshwaterFlux != 0): self.newLandIceFWFlux = np.zeros((nt,nc)) - + yearVec = np.zeros((nt,)) + #prepare time vector + stTime = np.datetime64(self.stTime[0:4]) + print("stTime = {}".format(stTime)) + stYear = np.datetime_as_string(stTime,unit='Y').astype('float64') + print("stYear = {}".format(stYear)) + print("starting loop ...") + st = time.time() if (years.ndim == 0): log = np.logical_and(yearsSinceStart >= years, yearsSinceStart < years + timeStride) @@ -120,6 +130,12 @@ def time_average_output(self): self.newLandIceFWFlux[0,:] = np.mean(self.landIceFreshwaterFlux[log,:], axis=0) else : self.newLandIceFWFlux = 0 + + #Define time at the first of each year + + print("Test A: {}".format(np.floor(np.min(yearsSinceStart[log])) + stYear)) + yearVec[0] = np.floor(np.min(yearsSinceStart[log])) + stYear + print("yearVec = {}".format(yearVec)) else : ct = 0 for i in years: @@ -134,16 +150,21 @@ def time_average_output(self): self.newLandIceFWFlux[ct,:] = np.mean(self.landIceFreshwaterFlux[log,:], axis=0) else : self.newLandIceFWFlux = 0 + + yearVec[ct] = np.floor(np.min(yearsSinceStart[log])) + stYear + print("yearVec = {}".format(yearVec)) ct = ct + 1 nd = time.time() tm = nd - st print("Time averaging loop", tm, "seconds") + + #establish xtime + dates = [] + xt = cftime.num2date(yearVec*365, units="days since 0000-01-01_00:00:00",calendar='noleap') + xt_reformat = [dt.strftime("%Y-%m-%d_%H:%M:%S") for dt in xt] + dates.append(xt_reformat) + self.newXtime = dates - # Build new xtime string. Defined on first day of each year - stTime = np.datetime64(self.stTime[0:4]) - stYear = np.datetime_as_string(stTime,unit='Y').astype('float64') - yearVec = yearsSinceStart + stYear - self.newXtime = np.array([f"{int(year)}-01-01_00:00:00" for year in yearVec],dtype=np.str_) else : #do nothing if not time averaging self.newTemp = self.temperature self.newSal = self.salinity @@ -152,8 +173,9 @@ def time_average_output(self): self.newMpasCCE = self.mpas_cellCenterElev self.newAtmPr = self.atmPressure self.newLandIceFWFlux = self.landIceFreshwaterFlux - self.newXtime = xtime + self.newXtime = np.nan #use as placeholder for now + print("New xtime: {}".format(self.newXtime)) def calc_ocean_thermal_forcing(self): print("Calculating thermal forcing ... ") gravity = 9.81 @@ -281,37 +303,33 @@ def remap_mpas_to_mali(self): mpas_cellCenterElev = interpDS['mpas_cellCenterElev'][:,:,:].data #Define vertical coordinates for mali grid - IM = xr.open_dataset(self.options.maliFile,decode_times=False,decode_cf=False) - layerThicknessFractions = IM['layerThicknessFractions'][:].data - bedTopography = IM['bedTopography'][:,:].data - thickness = IM['thickness'][:,:].data - nz, = layerThicknessFractions.shape - nt,nc = thickness.shape - - layerThicknessFractions = layerThicknessFractions.reshape((1,1,nz)) - thickness = thickness.reshape((nt,nc,1)) - bedTopography = bedTopography.reshape((nt,nc,1)) - - layerThickness = thickness * layerThicknessFractions - mali_cellCenterElev = np.cumsum(layerThickness,axis=2) - layerThickness/2 + bedTopography + #For now, hardcode in depths from ISMIP6 AIS. Should work fine for our purposes + ismip6shelfMelt_zOcean = np.array((-30, -90, -150, -210, -270, -330, -390, -450, -510,\ + -570, -630, -690, -750, -810, -870, -930, -990, -1050, -1110, -1170,\ + -1230, -1290, -1350, -1410, -1470, -1530, -1590, -1650, -1710, -1770),dtype='float64') + + nz, = ismip6shelfMelt_zOcean.shape # Reshape to (nt*nc, nz) before interpolation to avoid looping - mali_cellCenterElev = mali_cellCenterElev.reshape(nt*nc, -1) TF = np.transpose(TF,(0,2,1)) - TF = TF.reshape(nt*nc, -1) mpas_cellCenterElev = np.transpose(mpas_cellCenterElev,(0,2,1)) + nt,nc,_ = mpas_cellCenterElev.shape + print("mpas_cellCenterElev: {}".format(mpas_cellCenterElev.shape)) + print("TF: {}".format(TF.shape)) + print("nz: {}:".format(nz)) + TF = TF.reshape(nt*nc, -1) mpas_cellCenterElev = mpas_cellCenterElev.reshape(nt*nc, -1) # Linear interpolation st = time.time() nan_count = 0 - interpTF = np.zeros((mali_cellCenterElev.shape)) + interpTF = np.zeros((nt*nc,nz)) for i in range(nt*nc): ind = np.where(~np.isnan(TF[i,:].flatten() * mpas_cellCenterElev[i,:].flatten())) ct = 0 if len(ind[0]) != 0: - interpTF[i,:] = np.array(np.interp(mali_cellCenterElev[i,:].flatten(), mpas_cellCenterElev[i,ind].flatten(), TF[i,ind].flatten())) + interpTF[i,:] = np.array(np.interp(ismip6shelfMelt_zOcean, mpas_cellCenterElev[i,ind].flatten(), TF[i,ind].flatten())) else: nan_count = nan_count + 1 nd = time.time() @@ -329,29 +347,78 @@ def remap_mpas_to_mali(self): ind = np.where(surfaceTF != 0) validOpenOceanMask[ind] = 1 - if (self.newLandIceFWFlux != 0): - interpFWF = interpDS['floatingBasalMassBal'][:,:] - print("interpFWF: {}".format(interpFWF.data.shape)) print("Saving") + #open original mali file + IM = xr.open_dataset(self.options.maliFile,decode_times=False,decode_cf=False) - print("interpTF: {}".format(interpTF.shape)) - print("validOpenOceanMask: {}".format(validOpenOceanMask.shape)) - - mask = xr.DataArray(validOpenOceanMask,dims=("Time","nCells"),attrs={'long_name':'Mask of MALI cells overlapping interpolated MPAS-Oceans cells'}) + #create new mali forcing file that only desired forcing variables + ds_out = IM.copy(deep=False) + ds_out = ds_out.drop_vars(ds_out.data_vars) - IM = IM.expand_dims(nISMIP6OceanLayers=mali_cellCenterElev) - interptf = xr.DataArray(interpTF, dims=("Time","nCells","nISMIP6OceanLayers")) - zlayers = xr.DataArray(mali_cellCenterElev, dims=("nISMIP6OceanLayers")) + # introduce new ismip6 depth dimension + ds_out = ds_out.expand_dims(nISMIP6OceanLayers=ismip6shelfMelt_zOcean) # introduce new ismip6 depth dimension - IM['validOceanMask'] = mask + # Save variables + mask = xr.DataArray(validOpenOceanMask,dims=("Time","nCells"),attrs={'long_name':'Mask of MALI cells overlapping interpolated MPAS-Oceans cells'}) + interptf = xr.DataArray(interpTF, dims=("Time","nCells","nISMIP6OceanLayers")) + zlayers = xr.DataArray(ismip6shelfMelt_zOcean, dims=("nISMIP6OceanLayers")) + + ds_out['angleEdge'] = IM['angleEdge'] + ds_out['areaCell'] = IM['areaCell'] + ds_out['areaTriangle'] = IM['areaTriangle'] + ds_out['cellsOnCell'] = IM['cellsOnCell'] + ds_out['cellsOnEdge'] = IM['cellsOnEdge'] + ds_out['cellsOnVertex'] = IM['cellsOnVertex'] + ds_out['dcEdge'] = IM['dcEdge'] + ds_out['dvEdge'] = IM['dvEdge'] + ds_out['edgesOnCell'] = IM['edgesOnCell'] + ds_out['edgesOnEdge'] = IM['edgesOnEdge'] + ds_out['edgesOnVertex'] = IM['edgesOnVertex'] + ds_out['indexToCellID'] = IM['indexToCellID'] + ds_out['indexToEdgeID'] = IM['indexToEdgeID'] + ds_out['indexToVertexID'] = IM['indexToVertexID'] + ds_out['kiteAreasOnVertex'] = IM['kiteAreasOnVertex'] + ds_out['latCell'] = IM['latCell'] + ds_out['latEdge'] = IM['latEdge'] + ds_out['latVertex'] = IM['latVertex'] + ds_out['lonCell'] = IM['lonCell'] + ds_out['lonEdge'] = IM['lonEdge'] + ds_out['lonVertex'] = IM['lonVertex'] + ds_out['meshDensity'] = IM['meshDensity'] + ds_out['nEdgesOnCell'] = IM['nEdgesOnCell'] + ds_out['nEdgesOnEdge'] = IM['nEdgesOnEdge'] + ds_out['verticesOnCell'] = IM['verticesOnCell'] + ds_out['verticesOnEdge'] = IM['verticesOnEdge'] + ds_out['weightsOnEdge'] = IM['weightsOnEdge'] + ds_out['xCell'] = IM['xCell'] + ds_out['xEdge'] = IM['xEdge'] + ds_out['xVertex'] = IM['xVertex'] + ds_out['yCell'] = IM['yCell'] + ds_out['yEdge'] = IM['yEdge'] + ds_out['yVertex'] = IM['yVertex'] + ds_out['zCell'] = IM['zCell'] + ds_out['zEdge'] = IM['zEdge'] + ds_out['zVertex'] = IM['zVertex'] + + ds_out['validOceanMask'] = mask + ds_out['ismip6shelfMelt_3dThermalForcing'] = interptf + ds_out['ismip6shelfMelt_zOcean'] = zlayers if (self.newLandIceFWFlux != 0): - IM['floatingBasalMassBal'] = interpFWF - IM['ismip6shelfMelt_3dThermalForcing'] = interptf - IM['ismip6shelfMelt_zOcean'] = zlayers - - IM.to_netcdf(self.options.outputFile) - IM.close() + ds_out['floatingBasalMassBal'] = interpDS['floatingBasalMassBal'][:,:] + + # Save xtime + ds_out = ds_out.expand_dims({'StrLen': self.newXtime}, axis=1) + xtime = xr.DataArray(self.newXtime, dims=("Time", "StrLen")) + xtime.encoding.update({"char_dim_name":"StrLen"}) + ds_out['xtime'] = xtime + + #clean history for new file + if 'history' in ds_out.attrs: + del ds_out.attrs['history'] + + ds_out.to_netcdf(self.options.outputFile) + ds_out.close() subprocess.run(["ncatted", "-a", "_FillValue,,d,,", self.options.outputFile]) From 262acb644508bd2fbeea722d44ee840ea7af104c Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Tue, 25 Jun 2024 17:14:05 -0700 Subject: [PATCH 10/39] xtime debugging --- .../interpolate_mpasOcean_to_mali.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index c4ac27c1d..fa639884c 100644 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -163,8 +163,8 @@ def time_average_output(self): xt = cftime.num2date(yearVec*365, units="days since 0000-01-01_00:00:00",calendar='noleap') xt_reformat = [dt.strftime("%Y-%m-%d_%H:%M:%S") for dt in xt] dates.append(xt_reformat) - self.newXtime = dates - + #self.newXtime = dates + self.newXtime = xt_reformat else : #do nothing if not time averaging self.newTemp = self.temperature self.newSal = self.salinity @@ -408,16 +408,19 @@ def remap_mpas_to_mali(self): ds_out['floatingBasalMassBal'] = interpDS['floatingBasalMassBal'][:,:] # Save xtime - ds_out = ds_out.expand_dims({'StrLen': self.newXtime}, axis=1) - xtime = xr.DataArray(self.newXtime, dims=("Time", "StrLen")) - xtime.encoding.update({"char_dim_name":"StrLen"}) + #ds_out = ds_out.expand_dims({'StrLen': self.newXtime}, axis=1) + + print("self.newXtime: {}".format(self.newXtime)) + #ds_out = ds_out.expand_dims({'StrLen': self.newXtime}) + xtime = xr.DataArray(np.array(self.newXtime, dtype = np.dtype('S64')), dims=["Time"]) + xtime = xtime.encoding.update({"char_dim_name":"StrLen"}) ds_out['xtime'] = xtime #clean history for new file if 'history' in ds_out.attrs: del ds_out.attrs['history'] - ds_out.to_netcdf(self.options.outputFile) + ds_out.to_netcdf(self.options.outputFile, mode='w', unlimited_dims=['Time']) ds_out.close() subprocess.run(["ncatted", "-a", "_FillValue,,d,,", self.options.outputFile]) From c106ddbee3d66846ca6e65fc600a98b757f8cacd Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 26 Jun 2024 16:28:08 -0600 Subject: [PATCH 11/39] Fix xtime variable Two changes to get xtime written correctly: * when the char array is generated, pad it to 64 characters * fix a bug in that the return code to encoding.update was being assigned to xtime --- landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index fa639884c..1e0a94c2f 100644 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -161,7 +161,7 @@ def time_average_output(self): #establish xtime dates = [] xt = cftime.num2date(yearVec*365, units="days since 0000-01-01_00:00:00",calendar='noleap') - xt_reformat = [dt.strftime("%Y-%m-%d_%H:%M:%S") for dt in xt] + xt_reformat = [dt.strftime("%Y-%m-%d_%H:%M:%S").ljust(64) for dt in xt] dates.append(xt_reformat) #self.newXtime = dates self.newXtime = xt_reformat @@ -413,7 +413,7 @@ def remap_mpas_to_mali(self): print("self.newXtime: {}".format(self.newXtime)) #ds_out = ds_out.expand_dims({'StrLen': self.newXtime}) xtime = xr.DataArray(np.array(self.newXtime, dtype = np.dtype('S64')), dims=["Time"]) - xtime = xtime.encoding.update({"char_dim_name":"StrLen"}) + xtime.encoding.update({"char_dim_name":"StrLen"}) ds_out['xtime'] = xtime #clean history for new file From a825898392d725ba7a8a689f729e4ae29ea02c41 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Thu, 27 Jun 2024 16:26:59 -0500 Subject: [PATCH 12/39] Fix handling of if landIceFreshwaterFlux is available Previous code was giving an error when landIceFreshwaterFlux exists. --- landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index 1e0a94c2f..9c0710b5e 100644 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -51,11 +51,12 @@ def __init__(self): #xtime = np.array([xt],dtype=('S',64)) # << NOTE >>: may need to eventually use: "xtime.data.tobytes().decode()" but not sure yet + self.have_landIceFreshwaterFlux = False try: self.landIceFreshwaterFlux = DS['timeMonthly_avg_landIceFreshwaterFlux'][:,:].data + self.have_landIceFreshwaterFlux = True except KeyError: print("No landIceFreshwaterFlux variable") - self.landIceFreshwaterFlux = 0 # additional variables for computing thermal forcing bottomDepth = OM['bottomDepth'] @@ -106,7 +107,7 @@ def time_average_output(self): self.newLThick= np.zeros((nt,nc,nz)) self.newAtmPr = np.zeros((nt,nc)) self.newMpasCCE = np.zeros((nt,nc,nz)) - if (self.landIceFreshwaterFlux != 0): + if self.have_landIceFreshwaterFlux: self.newLandIceFWFlux = np.zeros((nt,nc)) yearVec = np.zeros((nt,)) #prepare time vector @@ -126,7 +127,7 @@ def time_average_output(self): self.newLThick[0,:,:] = np.mean(self.layerThickness[log,:,:], axis=0) self.newMpasCCE[0,:,:] = np.mean(self.mpas_cellCenterElev[log,:,:], axis=0) self.newAtmPr[0,:] = np.mean(self.atmPressure[log,:], axis=0) - if (self.landIceFreshwaterFlux != 0): + if self.have_landIceFreshwaterFlux: self.newLandIceFWFlux[0,:] = np.mean(self.landIceFreshwaterFlux[log,:], axis=0) else : self.newLandIceFWFlux = 0 @@ -146,7 +147,7 @@ def time_average_output(self): self.newLThick[ct,:,:] = np.mean(self.layerThickness[log,:,:], axis=0) self.newMpasCCE[ct,:,:] = np.mean(self.mpas_cellCenterElev[log,:,:], axis=0) self.newAtmPr[ct,:] = np.mean(self.atmPressure[log,:], axis=0) - if (self.landIceFreshwaterFlux != 0): + if self.have_landIceFreshwaterFlux: self.newLandIceFWFlux[ct,:] = np.mean(self.landIceFreshwaterFlux[log,:], axis=0) else : self.newLandIceFWFlux = 0 From 56bce681faf62bbb905b87ef5319085910fac647 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Thu, 27 Jun 2024 16:28:10 -0500 Subject: [PATCH 13/39] Set minLevelCell to 1 everywhere if variable is not present minLevelCell does not exist in older meshes, which assumed it was 1. --- landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index 9c0710b5e..7e862f6db 100644 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -61,10 +61,12 @@ def __init__(self): # additional variables for computing thermal forcing bottomDepth = OM['bottomDepth'] bD = bottomDepth.data - - self.minLevelCell = OM['minLevelCell'][:].data maxLevelCell = OM['maxLevelCell'] self.maxLevelCell = maxLevelCell.data + if 'minLevelCell' in OM: + self.minLevelCell = OM['minLevelCell'][:].data + else: + self.minLevelCell = self.maxLevelCell * 0 + 1 self.nCells = OM.sizes['nCells'] self.coeff_0_openOcean = DS.attrs['config_open_ocean_freezing_temperature_coeff_0'] From 94303474846583e207cfb0549d461e798a08eb63 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Thu, 27 Jun 2024 16:52:14 -0500 Subject: [PATCH 14/39] Switch to getting simulation start time from variable instead of global attr The global attribute version could sometimes have the value 'file'. --- landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index 7e862f6db..c3e6c95b6 100644 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -41,7 +41,8 @@ def __init__(self): self.density = DS['timeMonthly_avg_density'][:,:,:].compute() self.atmPressure = DS['timeMonthly_avg_atmosphericPressure'][:,:].compute() self.daysSinceStart = DS['timeMonthly_avg_daysSinceStartOfSim'][:].compute() - self.stTime = DS.attrs['config_start_time'] + self.stTime = OM['simulationStartTime'].data.tobytes().decode() + print(f'Using simulation start time of: {self.stTime}') self.landIceFloatingMask = OM['landIceFloatingMask'][0,:].data #not letting floating ice mask evolve for now because it's only in the mesh file avg_layerThickness = DS['timeMonthly_avg_layerThickness'] self.layerThickness = avg_layerThickness.data From 473393ca3275391346a684ca55c65ebb5110ef3a Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Thu, 27 Jun 2024 16:56:02 -0500 Subject: [PATCH 15/39] Further fix to handling of logic for if landiceFreshwaterFlux is present --- landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index c3e6c95b6..7dab84883 100644 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -132,8 +132,6 @@ def time_average_output(self): self.newAtmPr[0,:] = np.mean(self.atmPressure[log,:], axis=0) if self.have_landIceFreshwaterFlux: self.newLandIceFWFlux[0,:] = np.mean(self.landIceFreshwaterFlux[log,:], axis=0) - else : - self.newLandIceFWFlux = 0 #Define time at the first of each year @@ -152,8 +150,6 @@ def time_average_output(self): self.newAtmPr[ct,:] = np.mean(self.atmPressure[log,:], axis=0) if self.have_landIceFreshwaterFlux: self.newLandIceFWFlux[ct,:] = np.mean(self.landIceFreshwaterFlux[log,:], axis=0) - else : - self.newLandIceFWFlux = 0 yearVec[ct] = np.floor(np.min(yearsSinceStart[log])) + stYear print("yearVec = {}".format(yearVec)) @@ -176,7 +172,8 @@ def time_average_output(self): self.newLThick = self.layerThickness self.newMpasCCE = self.mpas_cellCenterElev self.newAtmPr = self.atmPressure - self.newLandIceFWFlux = self.landIceFreshwaterFlux + if self.have_landIceFreshwaterFlux: + self.newLandIceFWFlux = self.landIceFreshwaterFlux self.newXtime = np.nan #use as placeholder for now print("New xtime: {}".format(self.newXtime)) @@ -233,7 +230,7 @@ def remap_mpas_to_mali(self): mcce = xr.DataArray(self.newMpasCCE.astype('float64'),dims=('Time','nCells','nVertLevels')) f['mpas_cellCenterElev'] = mcce - if (self.newLandIceFWFlux != 0): + if self.have_landIceFreshwaterFlux: m = xr.DataArray(self.newLandIceFWFlux.astype('float64'),dims=('Time','nCells')) f['floatingBasalMassBal'] = m From f4dfa53461b4ce60adf1a8b3d59283096be0755c Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Thu, 27 Jun 2024 17:05:56 -0500 Subject: [PATCH 16/39] Fix integer type for mask variable --- landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index 7dab84883..4751ef173 100644 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -342,7 +342,7 @@ def remap_mpas_to_mali(self): # Create mask indentifying valid overlapping ocean cells. Combine all MALI terms in one input file #Making this a 3-D variable for now, but may want to make 2-D eventually - validOpenOceanMask = np.zeros((nt,nc), dtype='int') + validOpenOceanMask = np.zeros((nt,nc), dtype='int32') print("interpTF: {}".format(interpTF.data.shape)) surfaceTF = interpTF[:,:,0] ind = np.where(surfaceTF != 0) From 1e12de1b74ddfccb7c0c5d66d9ff3eb2e393e73f Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Thu, 27 Jun 2024 16:09:25 -0700 Subject: [PATCH 17/39] Fix remaining landIceFWFlux != 0 argument --- landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index 4751ef173..20645ba0f 100644 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -405,7 +405,8 @@ def remap_mpas_to_mali(self): ds_out['validOceanMask'] = mask ds_out['ismip6shelfMelt_3dThermalForcing'] = interptf ds_out['ismip6shelfMelt_zOcean'] = zlayers - if (self.newLandIceFWFlux != 0): + + if self.have_landIceFreshwaterFlux: ds_out['floatingBasalMassBal'] = interpDS['floatingBasalMassBal'][:,:] # Save xtime From c815154e27ad334b9eb8053fe909377e29e70b84 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 3 Jul 2024 07:45:11 -0500 Subject: [PATCH 18/39] correct sign of floatingBasalMassBal --- landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index 20645ba0f..b8cd46ccd 100644 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -231,8 +231,8 @@ def remap_mpas_to_mali(self): f['mpas_cellCenterElev'] = mcce if self.have_landIceFreshwaterFlux: - m = xr.DataArray(self.newLandIceFWFlux.astype('float64'),dims=('Time','nCells')) - f['floatingBasalMassBal'] = m + melt = xr.DataArray(self.newLandIceFWFlux.astype('float64'),dims=('Time','nCells')) + f['floatingBasalMassBal'] = -1.0 * melt #delete unncessary variables with 'string1' dimension or will cause errors. try: From 03afbfe48d8438c593d28b33929c53dbd991d085 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 3 Jul 2024 14:42:26 -0500 Subject: [PATCH 19/39] Restrict analysis to only *.mpaso.hist.am.timeSeriesStatsMonthly.*.nc files There will likely be other files in output directories, so ensure we are only working with the relevant ones. --- landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index b8cd46ccd..17450722c 100644 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -29,7 +29,7 @@ def __init__(self): self.options = options # open and concatenate diagnostic dataset - DS = xr.open_mfdataset(self.options.mpasDiagsDir + '/' + '*.nc', combine='nested', concat_dim='Time', decode_timedelta=False) + DS = xr.open_mfdataset(self.options.mpasDiagsDir + '/' + '*.mpaso.hist.am.timeSeriesStatsMonthly.*.nc', combine='nested', concat_dim='Time', decode_timedelta=False) #open mpas ocean mesh OM = xr.open_dataset(self.options.mpasMeshFile, decode_times=False, decode_cf=False) #open MALI mesh From 8fe90dd56fc823c19dbc1162293ef353980d2f90 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 3 Jul 2024 14:59:03 -0500 Subject: [PATCH 20/39] make script executable for user --- landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100644 => 100755 landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py old mode 100644 new mode 100755 From a5bb344efd10bb662c975ff53ff290f1c1b986e3 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 3 Jul 2024 16:12:22 -0500 Subject: [PATCH 21/39] correct bug in retrieving salinity --- landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index 17450722c..1276c0dae 100755 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -37,7 +37,7 @@ def __init__(self): # variables for time averaging self.temperature = DS['timeMonthly_avg_activeTracers_temperature'][:,:,:].compute() - self.salinity = DS['timeMonthly_avg_activeTracers_temperature'][:,:,:].compute() + self.salinity = DS['timeMonthly_avg_activeTracers_salinity'][:,:,:].compute() self.density = DS['timeMonthly_avg_density'][:,:,:].compute() self.atmPressure = DS['timeMonthly_avg_atmosphericPressure'][:,:].compute() self.daysSinceStart = DS['timeMonthly_avg_daysSinceStartOfSim'][:].compute() From c0dab1662532218695a864daa7e1493c5ab2b9fa Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Thu, 25 Jul 2024 13:26:45 -0500 Subject: [PATCH 22/39] Allow script to run multiple years MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * introduce start and end year options * add time loop to main function * create get_data function to retrieve needed data for each year * adjust remap_mpas_to_mali to handle a year argument A few comments: * I’ve commented out writing of the mesh fields to keep file size down. This is unrelated to the handling of years and perhaps should be a command line option * I had to adjust handling of various details of xarray data types to get this to work * creation of the mapping file should be moved out of the time loop so it only occurs once! --- .../interpolate_mpasOcean_to_mali.py | 194 ++++++++++-------- 1 file changed, 108 insertions(+), 86 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index 1276c0dae..cd1a44f77 100755 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -25,26 +25,42 @@ def __init__(self): parser.add_option("-m", "--method", dest="method", help="Remapping method, either 'bilinear' or 'neareststod'") parser.add_option("-o", "--outFile",dest="outputFile", help="Desired name of output file", metavar="FILENAME", default="mpas_to_mali_remapped.nc") parser.add_option("-a","--yearlyAvg", dest="yearlyAvg", help="true or false option to average monthly output to yearly intervals", default="true") + parser.add_option("-s","--startYr", dest="startYr", type="int", help="starting year to process") + parser.add_option("-e","--endYr", dest="endYr", type="int", help="ending year to process (inclusive)") options, args = parser.parse_args() self.options = options - # open and concatenate diagnostic dataset - DS = xr.open_mfdataset(self.options.mpasDiagsDir + '/' + '*.mpaso.hist.am.timeSeriesStatsMonthly.*.nc', combine='nested', concat_dim='Time', decode_timedelta=False) #open mpas ocean mesh OM = xr.open_dataset(self.options.mpasMeshFile, decode_times=False, decode_cf=False) + self.stTime = OM['simulationStartTime'].data.tobytes().decode() + print(f'Using simulation start time of: {self.stTime}') + self.landIceFloatingMask = OM['landIceFloatingMask'][0,:].data #not letting floating ice mask evolve for now because it's only in the mesh file #open MALI mesh IM = xr.open_dataset(self.options.maliFile, decode_times=False, decode_cf=False) + + # variables needed from MPAS-Ocean mesh file + self.bottomDepth = OM['bottomDepth'] + self.maxLevelCell = OM['maxLevelCell'] + if 'minLevelCell' in OM: + self.minLevelCell = OM['minLevelCell'][:].values + else: + self.minLevelCell = self.maxLevelCell * 0 + 1 + + self.nCells = OM.sizes['nCells'] + + def get_data(self, year): + + # open and concatenate diagnostic dataset + self.DS = xr.open_mfdataset(os.path.join(self.options.mpasDiagsDir, + f'*.mpaso.hist.am.timeSeriesStatsMonthly.{year}-*-01.nc'), combine='nested', concat_dim='Time', decode_timedelta=False) # variables for time averaging - self.temperature = DS['timeMonthly_avg_activeTracers_temperature'][:,:,:].compute() - self.salinity = DS['timeMonthly_avg_activeTracers_salinity'][:,:,:].compute() - self.density = DS['timeMonthly_avg_density'][:,:,:].compute() - self.atmPressure = DS['timeMonthly_avg_atmosphericPressure'][:,:].compute() - self.daysSinceStart = DS['timeMonthly_avg_daysSinceStartOfSim'][:].compute() - self.stTime = OM['simulationStartTime'].data.tobytes().decode() - print(f'Using simulation start time of: {self.stTime}') - self.landIceFloatingMask = OM['landIceFloatingMask'][0,:].data #not letting floating ice mask evolve for now because it's only in the mesh file - avg_layerThickness = DS['timeMonthly_avg_layerThickness'] + self.temperature = self.DS['timeMonthly_avg_activeTracers_temperature'] + self.salinity = self.DS['timeMonthly_avg_activeTracers_salinity'] + self.density = self.DS['timeMonthly_avg_density'] + self.atmPressure =self. DS['timeMonthly_avg_atmosphericPressure'] + self.daysSinceStart = self.DS['timeMonthly_avg_daysSinceStartOfSim'] + avg_layerThickness = self.DS['timeMonthly_avg_layerThickness'] self.layerThickness = avg_layerThickness.data #xt = DS['xtime_startMonthly'] @@ -54,46 +70,37 @@ def __init__(self): self.have_landIceFreshwaterFlux = False try: - self.landIceFreshwaterFlux = DS['timeMonthly_avg_landIceFreshwaterFlux'][:,:].data + self.landIceFreshwaterFlux = self.DS['timeMonthly_avg_landIceFreshwaterFlux'] self.have_landIceFreshwaterFlux = True except KeyError: print("No landIceFreshwaterFlux variable") - # additional variables for computing thermal forcing - bottomDepth = OM['bottomDepth'] - bD = bottomDepth.data - maxLevelCell = OM['maxLevelCell'] - self.maxLevelCell = maxLevelCell.data - if 'minLevelCell' in OM: - self.minLevelCell = OM['minLevelCell'][:].data - else: - self.minLevelCell = self.maxLevelCell * 0 + 1 - - self.nCells = OM.sizes['nCells'] - self.coeff_0_openOcean = DS.attrs['config_open_ocean_freezing_temperature_coeff_0'] - self.coeff_S_openOcean = DS.attrs['config_open_ocean_freezing_temperature_coeff_S'] - self.coeff_p_openOcean = DS.attrs['config_open_ocean_freezing_temperature_coeff_p'] - self.coeff_pS_openOcean = DS.attrs['config_open_ocean_freezing_temperature_coeff_pS'] - az1_liq = DS.attrs['config_open_ocean_freezing_temperature_coeff_mushy_az1_liq'] + self.coeff_0_openOcean = self.DS.attrs['config_open_ocean_freezing_temperature_coeff_0'] + self.coeff_S_openOcean = self.DS.attrs['config_open_ocean_freezing_temperature_coeff_S'] + self.coeff_p_openOcean = self.DS.attrs['config_open_ocean_freezing_temperature_coeff_p'] + self.coeff_pS_openOcean = self.DS.attrs['config_open_ocean_freezing_temperature_coeff_pS'] + az1_liq = self.DS.attrs['config_open_ocean_freezing_temperature_coeff_mushy_az1_liq'] self.coeff_mushy_openOcean = 1/az1_liq - self.coeff_0_cavity = DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_0'] - self.coeff_S_cavity = DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_S'] - self.coeff_p_cavity = DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_p'] - self.coeff_pS_cavity = DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_pS'] + self.coeff_0_cavity = self.DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_0'] + self.coeff_S_cavity = self.DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_S'] + self.coeff_p_cavity = self.DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_p'] + self.coeff_pS_cavity = self.DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_pS'] self.coeff_mushy_cavity = 0 # Define vertical coordinates of mpas output - mpas_cellCenterElev = compute_zmid(bottomDepth, maxLevelCell, avg_layerThickness) + mpas_cellCenterElev = compute_zmid(self.bottomDepth, self.maxLevelCell, avg_layerThickness) self.mpas_cellCenterElev = mpas_cellCenterElev.data def time_average_output(self): print("Time Averaging ...") if (self.options.yearlyAvg == 'true'): - yearsSinceStart = self.daysSinceStart.data / 365.0 + yearsSinceStart = self.daysSinceStart.values / 365.0 finalYear = np.floor(np.max(yearsSinceStart)) startYear = np.floor(np.min(yearsSinceStart)) timeStride = 1 # 1 year average + + print(f'startYear={startYear}, finalYear={finalYear}') if (startYear != finalYear): years = np.arange(startYear, finalYear, timeStride, dtype=int) @@ -102,6 +109,8 @@ def time_average_output(self): years = startYear nt = 1 + print(f'years={years}', nt) + #pre-allocate _,nc,nz = self.temperature.shape self.newTemp = np.zeros((nt,nc,nz)) @@ -124,14 +133,19 @@ def time_average_output(self): st = time.time() if (years.ndim == 0): log = np.logical_and(yearsSinceStart >= years, yearsSinceStart < years + timeStride) - self.newTemp[0,:,:] = np.mean(self.temperature[log,:,:], axis=0) - self.newSal[0,:,:] = np.mean(self.salinity[log,:,:], axis=0) - self.newDens[0,:,:] = np.mean(self.density[log,:,:], axis=0) + #ind1 = np.nonzero(yearsSinceStart >= years)[0][0] + #ind2 = np.nonzero(yearsSinceStart < years + timeStride)[0][-1] + #print(ind1, ind2) + #self.newTemp[0,:,:] = np.mean(self.temperature[ind1:ind2+1,:,:], axis=0) + #self.newTemp[0,:,:] = np.mean(self.temperature[log,:,:], axis=0) + self.newTemp[0,:,:] = np.mean(self.temperature[log,:,:].values, axis=0) + self.newSal[0,:,:] = np.mean(self.salinity[log,:,:].values, axis=0) + self.newDens[0,:,:] = np.mean(self.density[log,:,:].values, axis=0) self.newLThick[0,:,:] = np.mean(self.layerThickness[log,:,:], axis=0) self.newMpasCCE[0,:,:] = np.mean(self.mpas_cellCenterElev[log,:,:], axis=0) - self.newAtmPr[0,:] = np.mean(self.atmPressure[log,:], axis=0) + self.newAtmPr[0,:] = np.mean(self.atmPressure[log,:].values, axis=0) if self.have_landIceFreshwaterFlux: - self.newLandIceFWFlux[0,:] = np.mean(self.landIceFreshwaterFlux[log,:], axis=0) + self.newLandIceFWFlux[0,:] = np.mean(self.landIceFreshwaterFlux[log,:].values, axis=0) #Define time at the first of each year @@ -191,7 +205,7 @@ def calc_ocean_thermal_forcing(self): #calculate pressure: for iCell in range(self.nCells): kmin = self.minLevelCell[iCell] - 1 - kmax = self.maxLevelCell[iCell] - 1 + kmax = self.maxLevelCell[iCell].values - 1 self.newPr[iTime,iCell,kmin] = self.newAtmPr[iTime,iCell] + \ self.newDens[iTime,iCell,kmin]*gravity*0.5*self.newLThick[iTime,iCell,kmin] @@ -218,7 +232,7 @@ def calc_ocean_thermal_forcing(self): # Calculate thermal forcing self.oceanThermalForcing = self.newTemp - ocn_freezing_temperature - def remap_mpas_to_mali(self): + def remap_mpas_to_mali(self, year): print("Start remapping ... ") # populate tmp_mpasMeshFile with variables to be interpolated @@ -365,42 +379,42 @@ def remap_mpas_to_mali(self): interptf = xr.DataArray(interpTF, dims=("Time","nCells","nISMIP6OceanLayers")) zlayers = xr.DataArray(ismip6shelfMelt_zOcean, dims=("nISMIP6OceanLayers")) - ds_out['angleEdge'] = IM['angleEdge'] - ds_out['areaCell'] = IM['areaCell'] - ds_out['areaTriangle'] = IM['areaTriangle'] - ds_out['cellsOnCell'] = IM['cellsOnCell'] - ds_out['cellsOnEdge'] = IM['cellsOnEdge'] - ds_out['cellsOnVertex'] = IM['cellsOnVertex'] - ds_out['dcEdge'] = IM['dcEdge'] - ds_out['dvEdge'] = IM['dvEdge'] - ds_out['edgesOnCell'] = IM['edgesOnCell'] - ds_out['edgesOnEdge'] = IM['edgesOnEdge'] - ds_out['edgesOnVertex'] = IM['edgesOnVertex'] - ds_out['indexToCellID'] = IM['indexToCellID'] - ds_out['indexToEdgeID'] = IM['indexToEdgeID'] - ds_out['indexToVertexID'] = IM['indexToVertexID'] - ds_out['kiteAreasOnVertex'] = IM['kiteAreasOnVertex'] - ds_out['latCell'] = IM['latCell'] - ds_out['latEdge'] = IM['latEdge'] - ds_out['latVertex'] = IM['latVertex'] - ds_out['lonCell'] = IM['lonCell'] - ds_out['lonEdge'] = IM['lonEdge'] - ds_out['lonVertex'] = IM['lonVertex'] - ds_out['meshDensity'] = IM['meshDensity'] - ds_out['nEdgesOnCell'] = IM['nEdgesOnCell'] - ds_out['nEdgesOnEdge'] = IM['nEdgesOnEdge'] - ds_out['verticesOnCell'] = IM['verticesOnCell'] - ds_out['verticesOnEdge'] = IM['verticesOnEdge'] - ds_out['weightsOnEdge'] = IM['weightsOnEdge'] - ds_out['xCell'] = IM['xCell'] - ds_out['xEdge'] = IM['xEdge'] - ds_out['xVertex'] = IM['xVertex'] - ds_out['yCell'] = IM['yCell'] - ds_out['yEdge'] = IM['yEdge'] - ds_out['yVertex'] = IM['yVertex'] - ds_out['zCell'] = IM['zCell'] - ds_out['zEdge'] = IM['zEdge'] - ds_out['zVertex'] = IM['zVertex'] + #ds_out['angleEdge'] = IM['angleEdge'] + #ds_out['areaCell'] = IM['areaCell'] + #ds_out['areaTriangle'] = IM['areaTriangle'] + #ds_out['cellsOnCell'] = IM['cellsOnCell'] + #ds_out['cellsOnEdge'] = IM['cellsOnEdge'] + #ds_out['cellsOnVertex'] = IM['cellsOnVertex'] + #ds_out['dcEdge'] = IM['dcEdge'] + #ds_out['dvEdge'] = IM['dvEdge'] + #ds_out['edgesOnCell'] = IM['edgesOnCell'] + #ds_out['edgesOnEdge'] = IM['edgesOnEdge'] + #ds_out['edgesOnVertex'] = IM['edgesOnVertex'] + #ds_out['indexToCellID'] = IM['indexToCellID'] + #ds_out['indexToEdgeID'] = IM['indexToEdgeID'] + #ds_out['indexToVertexID'] = IM['indexToVertexID'] + #ds_out['kiteAreasOnVertex'] = IM['kiteAreasOnVertex'] + #ds_out['latCell'] = IM['latCell'] + #ds_out['latEdge'] = IM['latEdge'] + #ds_out['latVertex'] = IM['latVertex'] + #ds_out['lonCell'] = IM['lonCell'] + #ds_out['lonEdge'] = IM['lonEdge'] + #ds_out['lonVertex'] = IM['lonVertex'] + #ds_out['meshDensity'] = IM['meshDensity'] + #ds_out['nEdgesOnCell'] = IM['nEdgesOnCell'] + #ds_out['nEdgesOnEdge'] = IM['nEdgesOnEdge'] + #ds_out['verticesOnCell'] = IM['verticesOnCell'] + #ds_out['verticesOnEdge'] = IM['verticesOnEdge'] + #ds_out['weightsOnEdge'] = IM['weightsOnEdge'] + #ds_out['xCell'] = IM['xCell'] + #ds_out['xEdge'] = IM['xEdge'] + #ds_out['xVertex'] = IM['xVertex'] + #ds_out['yCell'] = IM['yCell'] + #ds_out['yEdge'] = IM['yEdge'] + #ds_out['yVertex'] = IM['yVertex'] + #ds_out['zCell'] = IM['zCell'] + #ds_out['zEdge'] = IM['zEdge'] + #ds_out['zVertex'] = IM['zVertex'] ds_out['validOceanMask'] = mask ds_out['ismip6shelfMelt_3dThermalForcing'] = interptf @@ -422,10 +436,11 @@ def remap_mpas_to_mali(self): if 'history' in ds_out.attrs: del ds_out.attrs['history'] - ds_out.to_netcdf(self.options.outputFile, mode='w', unlimited_dims=['Time']) + out_name = f'{self.options.outputFile}_{year}.nc' + ds_out.to_netcdf(out_name, mode='w', unlimited_dims=['Time']) ds_out.close() - subprocess.run(["ncatted", "-a", "_FillValue,,d,,", self.options.outputFile]) + subprocess.run(["ncatted", "-a", "_FillValue,,d,,", out_name]) #remove temporary files files = os.listdir() @@ -435,16 +450,23 @@ def remap_mpas_to_mali(self): def main(): run = mpasToMaliInterp() + + for yr in range(run.options.startYr, run.options.endYr+1): + print(f'### Processing year {yr}') + + run.get_data(yr) - #compute yearly time average if necessary - run.time_average_output() + #compute yearly time average if necessary + run.time_average_output() - #calculate thermal forcing - run.calc_ocean_thermal_forcing() + #calculate thermal forcing + run.calc_ocean_thermal_forcing() + + #remap mpas to mali + run.remap_mpas_to_mali(yr) + + run.DS.close() - #remap mpas to mali - run.remap_mpas_to_mali() - print("Finished.") if __name__ == "__main__": main() From 4ae6f86c6d917e835bcfd067b2df2555827b5fb5 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Fri, 26 Jul 2024 12:25:23 -0500 Subject: [PATCH 23/39] Switch to argparse and add option for including MALI mesh vars or not --- .../interpolate_mpasOcean_to_mali.py | 102 +++++++++--------- 1 file changed, 53 insertions(+), 49 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index cd1a44f77..dd283779c 100755 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -9,7 +9,7 @@ from mpas_tools.logging import check_call from mpas_tools.scrip.from_mpas import scrip_from_mpas from mpas_tools.ocean.depth import compute_zmid -from optparse import OptionParser +from argparse import ArgumentParser import time import cftime @@ -17,18 +17,21 @@ class mpasToMaliInterp: def __init__(self): print("Gathering Information ... ") - parser = OptionParser() - parser.add_option("--oceanMesh", dest="mpasMeshFile", help="MPAS-Ocean mesh to be interpolated from", metavar="FILENAME") - parser.add_option("--oceanDiags", dest="mpasDiagsDir", help="Directory path where MPAS-Ocean diagnostic files are stored. Should be only netcdf files in that directory", metavar="FILENAME") - parser.add_option("--ice", dest="maliFile", help="MALI mesh to be interpolated onto", metavar="FILENAME") - parser.add_option("-n", "--ntasks", dest="ntasks", help="Number of processors to use with ESMF_RegridWeightGen") - parser.add_option("-m", "--method", dest="method", help="Remapping method, either 'bilinear' or 'neareststod'") - parser.add_option("-o", "--outFile",dest="outputFile", help="Desired name of output file", metavar="FILENAME", default="mpas_to_mali_remapped.nc") - parser.add_option("-a","--yearlyAvg", dest="yearlyAvg", help="true or false option to average monthly output to yearly intervals", default="true") - parser.add_option("-s","--startYr", dest="startYr", type="int", help="starting year to process") - parser.add_option("-e","--endYr", dest="endYr", type="int", help="ending year to process (inclusive)") - options, args = parser.parse_args() - self.options = options + parser = ArgumentParser( + prog='interpolate_mpasOcean_to_mali.py', + description='Interpolates between MPAS-Ocean data to melt rates and thermal forcing on a MALI mesh') + parser.add_argument("--oceanMesh", dest="mpasMeshFile", help="MPAS-Ocean mesh to be interpolated from", metavar="FILENAME") + parser.add_argument("--oceanDiags", dest="mpasDiagsDir", help="Directory path where MPAS-Ocean diagnostic files are stored. Should be only netcdf files in that directory", metavar="FILENAME") + parser.add_argument("--ice", dest="maliFile", help="MALI mesh to be interpolated onto", metavar="FILENAME") + parser.add_argument("-n", "--ntasks", dest="ntasks", help="Number of processors to use with ESMF_RegridWeightGen") + parser.add_argument("-m", "--method", dest="method", help="Remapping method, either 'bilinear' or 'neareststod'") + parser.add_argument("-o", "--outFile", dest="outputFile", help="Desired name of output file", metavar="FILENAME", default="mpas_to_mali_remapped.nc") + parser.add_argument("-a","--yearlyAvg", dest="yearlyAvg", help="true or false option to average monthly output to yearly intervals", default="true") + parser.add_argument("-s","--startYr", dest="startYr", type=int, help="starting year to process") + parser.add_argument("-e","--endYr", dest="endYr", type=int, help="ending year to process (inclusive)") + parser.add_argument("--meshVars", dest="includeMeshVars", help="whether to include mesh variables in resulting MALI forcing file", action='store_true') + args = parser.parse_args() + self.options = args #open mpas ocean mesh OM = xr.open_dataset(self.options.mpasMeshFile, decode_times=False, decode_cf=False) @@ -379,42 +382,43 @@ def remap_mpas_to_mali(self, year): interptf = xr.DataArray(interpTF, dims=("Time","nCells","nISMIP6OceanLayers")) zlayers = xr.DataArray(ismip6shelfMelt_zOcean, dims=("nISMIP6OceanLayers")) - #ds_out['angleEdge'] = IM['angleEdge'] - #ds_out['areaCell'] = IM['areaCell'] - #ds_out['areaTriangle'] = IM['areaTriangle'] - #ds_out['cellsOnCell'] = IM['cellsOnCell'] - #ds_out['cellsOnEdge'] = IM['cellsOnEdge'] - #ds_out['cellsOnVertex'] = IM['cellsOnVertex'] - #ds_out['dcEdge'] = IM['dcEdge'] - #ds_out['dvEdge'] = IM['dvEdge'] - #ds_out['edgesOnCell'] = IM['edgesOnCell'] - #ds_out['edgesOnEdge'] = IM['edgesOnEdge'] - #ds_out['edgesOnVertex'] = IM['edgesOnVertex'] - #ds_out['indexToCellID'] = IM['indexToCellID'] - #ds_out['indexToEdgeID'] = IM['indexToEdgeID'] - #ds_out['indexToVertexID'] = IM['indexToVertexID'] - #ds_out['kiteAreasOnVertex'] = IM['kiteAreasOnVertex'] - #ds_out['latCell'] = IM['latCell'] - #ds_out['latEdge'] = IM['latEdge'] - #ds_out['latVertex'] = IM['latVertex'] - #ds_out['lonCell'] = IM['lonCell'] - #ds_out['lonEdge'] = IM['lonEdge'] - #ds_out['lonVertex'] = IM['lonVertex'] - #ds_out['meshDensity'] = IM['meshDensity'] - #ds_out['nEdgesOnCell'] = IM['nEdgesOnCell'] - #ds_out['nEdgesOnEdge'] = IM['nEdgesOnEdge'] - #ds_out['verticesOnCell'] = IM['verticesOnCell'] - #ds_out['verticesOnEdge'] = IM['verticesOnEdge'] - #ds_out['weightsOnEdge'] = IM['weightsOnEdge'] - #ds_out['xCell'] = IM['xCell'] - #ds_out['xEdge'] = IM['xEdge'] - #ds_out['xVertex'] = IM['xVertex'] - #ds_out['yCell'] = IM['yCell'] - #ds_out['yEdge'] = IM['yEdge'] - #ds_out['yVertex'] = IM['yVertex'] - #ds_out['zCell'] = IM['zCell'] - #ds_out['zEdge'] = IM['zEdge'] - #ds_out['zVertex'] = IM['zVertex'] + if self.options.includeMeshVars: + ds_out['angleEdge'] = IM['angleEdge'] + ds_out['areaCell'] = IM['areaCell'] + ds_out['areaTriangle'] = IM['areaTriangle'] + ds_out['cellsOnCell'] = IM['cellsOnCell'] + ds_out['cellsOnEdge'] = IM['cellsOnEdge'] + ds_out['cellsOnVertex'] = IM['cellsOnVertex'] + ds_out['dcEdge'] = IM['dcEdge'] + ds_out['dvEdge'] = IM['dvEdge'] + ds_out['edgesOnCell'] = IM['edgesOnCell'] + ds_out['edgesOnEdge'] = IM['edgesOnEdge'] + ds_out['edgesOnVertex'] = IM['edgesOnVertex'] + ds_out['indexToCellID'] = IM['indexToCellID'] + ds_out['indexToEdgeID'] = IM['indexToEdgeID'] + ds_out['indexToVertexID'] = IM['indexToVertexID'] + ds_out['kiteAreasOnVertex'] = IM['kiteAreasOnVertex'] + ds_out['latCell'] = IM['latCell'] + ds_out['latEdge'] = IM['latEdge'] + ds_out['latVertex'] = IM['latVertex'] + ds_out['lonCell'] = IM['lonCell'] + ds_out['lonEdge'] = IM['lonEdge'] + ds_out['lonVertex'] = IM['lonVertex'] + ds_out['meshDensity'] = IM['meshDensity'] + ds_out['nEdgesOnCell'] = IM['nEdgesOnCell'] + ds_out['nEdgesOnEdge'] = IM['nEdgesOnEdge'] + ds_out['verticesOnCell'] = IM['verticesOnCell'] + ds_out['verticesOnEdge'] = IM['verticesOnEdge'] + ds_out['weightsOnEdge'] = IM['weightsOnEdge'] + ds_out['xCell'] = IM['xCell'] + ds_out['xEdge'] = IM['xEdge'] + ds_out['xVertex'] = IM['xVertex'] + ds_out['yCell'] = IM['yCell'] + ds_out['yEdge'] = IM['yEdge'] + ds_out['yVertex'] = IM['yVertex'] + ds_out['zCell'] = IM['zCell'] + ds_out['zEdge'] = IM['zEdge'] + ds_out['zVertex'] = IM['zVertex'] ds_out['validOceanMask'] = mask ds_out['ismip6shelfMelt_3dThermalForcing'] = interptf From ce2d473f2b10887b150aae9eba9cfa9b0659a7ab Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Fri, 26 Jul 2024 12:36:41 -0500 Subject: [PATCH 24/39] Move creation of mapping file outside of time loop This only needs to happen once. This commit also changes some variable names to be more precise. --- .../interpolate_mpasOcean_to_mali.py | 75 ++++++++++--------- 1 file changed, 41 insertions(+), 34 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index dd283779c..c5df11be7 100755 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -20,8 +20,8 @@ def __init__(self): parser = ArgumentParser( prog='interpolate_mpasOcean_to_mali.py', description='Interpolates between MPAS-Ocean data to melt rates and thermal forcing on a MALI mesh') - parser.add_argument("--oceanMesh", dest="mpasMeshFile", help="MPAS-Ocean mesh to be interpolated from", metavar="FILENAME") - parser.add_argument("--oceanDiags", dest="mpasDiagsDir", help="Directory path where MPAS-Ocean diagnostic files are stored. Should be only netcdf files in that directory", metavar="FILENAME") + parser.add_argument("--oceanMesh", dest="mpasoMeshFile", help="MPAS-Ocean mesh to be interpolated from", metavar="FILENAME") + parser.add_argument("--oceanDiags", dest="mpasoHistDir", help="Directory where MPAS-Ocean history files are stored. Should be only netcdf files in that directory", metavar="FILENAME") parser.add_argument("--ice", dest="maliFile", help="MALI mesh to be interpolated onto", metavar="FILENAME") parser.add_argument("-n", "--ntasks", dest="ntasks", help="Number of processors to use with ESMF_RegridWeightGen") parser.add_argument("-m", "--method", dest="method", help="Remapping method, either 'bilinear' or 'neareststod'") @@ -34,7 +34,7 @@ def __init__(self): self.options = args #open mpas ocean mesh - OM = xr.open_dataset(self.options.mpasMeshFile, decode_times=False, decode_cf=False) + OM = xr.open_dataset(self.options.mpasoMeshFile, decode_times=False, decode_cf=False) self.stTime = OM['simulationStartTime'].data.tobytes().decode() print(f'Using simulation start time of: {self.stTime}') self.landIceFloatingMask = OM['landIceFloatingMask'][0,:].data #not letting floating ice mask evolve for now because it's only in the mesh file @@ -51,10 +51,12 @@ def __init__(self): self.nCells = OM.sizes['nCells'] + self.mapping_file_name = f'mpaso_to_mali_mapping_{self.options.method}.nc' + def get_data(self, year): # open and concatenate diagnostic dataset - self.DS = xr.open_mfdataset(os.path.join(self.options.mpasDiagsDir, + self.DS = xr.open_mfdataset(os.path.join(self.options.mpasoHistDir, f'*.mpaso.hist.am.timeSeriesStatsMonthly.{year}-*-01.nc'), combine='nested', concat_dim='Time', decode_timedelta=False) # variables for time averaging @@ -194,6 +196,7 @@ def time_average_output(self): self.newXtime = np.nan #use as placeholder for now print("New xtime: {}".format(self.newXtime)) + def calc_ocean_thermal_forcing(self): print("Calculating thermal forcing ... ") gravity = 9.81 @@ -234,12 +237,34 @@ def calc_ocean_thermal_forcing(self): print("Ocean thermal forcing loop", tm, "seconds") # Calculate thermal forcing self.oceanThermalForcing = self.newTemp - ocn_freezing_temperature + + def create_mapping_file(self): + + #create scrip files + print("Creating Scrip Files") + mpaso_scripfile = "tmp_mpaso_scrip.nc" + mali_scripfile = "tmp_mali_scrip.nc" + + scrip_from_mpas(self.options.mpasoMeshFile, mpaso_scripfile) + scrip_from_mpas(self.options.maliFile, mali_scripfile) + + #create mapping file + print("Creating Mapping File") + args_esmf = ['srun', + '-n', self.options.ntasks, 'ESMF_RegridWeightGen', + '--source', mpaso_scripfile, + '--destination', mali_scripfile, + '--weight', self.mapping_file_name, + '--method', self.options.method, + '-i', '-64bit_offset', + "--dst_regional", "--src_regional", '--ignore_unmapped'] + check_call(args_esmf) def remap_mpas_to_mali(self, year): print("Start remapping ... ") - # populate tmp_mpasMeshFile with variables to be interpolated - f = xr.open_dataset(self.options.mpasMeshFile, decode_times=False, decode_cf=False) + # populate tmp_mpasoMeshFile with variables to be interpolated + f = xr.open_dataset(self.options.mpasoMeshFile, decode_times=False, decode_cf=False) tf = xr.DataArray(self.oceanThermalForcing.astype('float64'),dims=('Time','nCells','nVertLevels')) f['ismip6shelfMelt_3dThermalForcing'] = tf @@ -269,49 +294,29 @@ def remap_mpas_to_mali(self, year): st = time.time() # save new mesh file - tmp_mpasMeshFile = "tmp_mpasMeshFile.nc" - f.to_netcdf(tmp_mpasMeshFile) + tmp_mpasoMeshFile = "tmp_mpasoMeshFile.nc" + f.to_netcdf(tmp_mpasoMeshFile) f.close() nd = time.time() tm = nd - st print("saving updates mpas mesh file", tm, "seconds") st = time.time() - subprocess.run(["ncatted", "-a", "_FillValue,,d,,", tmp_mpasMeshFile]) + subprocess.run(["ncatted", "-a", "_FillValue,,d,,", tmp_mpasoMeshFile]) nd = time.time() tm = nd - st print("Removing fill value:", tm, "seconds") - #create scrip files - print("Creating Scrip Files") - mali_scripfile = "tmp_mali_scrip.nc" - mpas_scripfile = "tmp_mpas_scrip.nc" - - scrip_from_mpas(self.options.maliFile, mali_scripfile) - scrip_from_mpas(tmp_mpasMeshFile, mpas_scripfile) - - #create mapping file - print("Creating Mapping File") - args_esmf = ['srun', - '-n', self.options.ntasks, 'ESMF_RegridWeightGen', - '--source', mpas_scripfile, - '--destination', mali_scripfile, - '--weight', "tmp_mapping_file.nc", - '--method', self.options.method, - '-i', '-64bit_offset', - "--dst_regional", "--src_regional", '--ignore_unmapped'] - check_call(args_esmf) - - subprocess.run(["ncpdq", "-O", "-a", "Time,nVertLevels,nCells", tmp_mpasMeshFile, tmp_mpasMeshFile]) - subprocess.run(["ncpdq", "-O", "-a", "Time,nVertLevelsP1,nCells", tmp_mpasMeshFile, tmp_mpasMeshFile]) - subprocess.run(["ncpdq", "-O", "-a", "Time,nVertInterfaces,nCells", tmp_mpasMeshFile, tmp_mpasMeshFile]) + subprocess.run(["ncpdq", "-O", "-a", "Time,nVertLevels,nCells", tmp_mpasoMeshFile, tmp_mpasoMeshFile]) + subprocess.run(["ncpdq", "-O", "-a", "Time,nVertLevelsP1,nCells", tmp_mpasoMeshFile, tmp_mpasoMeshFile]) + subprocess.run(["ncpdq", "-O", "-a", "Time,nVertInterfaces,nCells", tmp_mpasoMeshFile, tmp_mpasoMeshFile]) print("ncremap ...") # remap the input data args_remap = ["ncremap", - "-i", tmp_mpasMeshFile, + "-i", tmp_mpasoMeshFile, "-o", "tmp_outputFile.nc", - "-m", "tmp_mapping_file.nc"] + "-m", self.mapping_file_name] check_call(args_remap) #Vertical interpolation @@ -455,6 +460,8 @@ def remap_mpas_to_mali(self, year): def main(): run = mpasToMaliInterp() + run.create_mapping_file() + for yr in range(run.options.startYr, run.options.endYr+1): print(f'### Processing year {yr}') From 80fc72e876f07ed4dbc4ef6fa875867a1dc46010 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 29 Jul 2024 11:34:42 -0500 Subject: [PATCH 25/39] Move mask creation to before interpolation in a separate output file This commit creates 2 masks about the extent of ocean data and moves them to before interpolation so that MALI has access to the interpolated mask values, rather than trying to construct them after the fact. The masks get interpolated in a separate file from the TF/melt data. It is not yet clear if there will be a single mask definition that is appropriate for all applications, so this approach forces the MALI user to make a conscious decision about how to define the mask. The two masks are: * the entire ocean domain * the open ocean --- .../interpolate_mpasOcean_to_mali.py | 44 ++++++++++++------- 1 file changed, 28 insertions(+), 16 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index c5df11be7..af851e9f0 100755 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -33,13 +33,15 @@ def __init__(self): args = parser.parse_args() self.options = args + self.mapping_file_name = f'mpaso_to_mali_mapping_{self.options.method}.nc' + + def prepare_mpaso_mesh_data(self): + #open mpas ocean mesh OM = xr.open_dataset(self.options.mpasoMeshFile, decode_times=False, decode_cf=False) self.stTime = OM['simulationStartTime'].data.tobytes().decode() print(f'Using simulation start time of: {self.stTime}') self.landIceFloatingMask = OM['landIceFloatingMask'][0,:].data #not letting floating ice mask evolve for now because it's only in the mesh file - #open MALI mesh - IM = xr.open_dataset(self.options.maliFile, decode_times=False, decode_cf=False) # variables needed from MPAS-Ocean mesh file self.bottomDepth = OM['bottomDepth'] @@ -51,7 +53,24 @@ def __init__(self): self.nCells = OM.sizes['nCells'] - self.mapping_file_name = f'mpaso_to_mali_mapping_{self.options.method}.nc' + # create needed masks and remap them in separate file + mpasoDomainMask = np.ones((self.nCells,), dtype=np.double) + mpasoDomainMaskDA = xr.DataArray(mpasoDomainMask, name='mpasoDomainMask', dims=("nCells"), attrs={'long_name':'MPAS-Ocean domain mask'}) + mpasoOpenOceanMaskDA = xr.DataArray(np.logical_not(self.landIceFloatingMask).astype(np.double), name='mpasoOpenOceanMask', dims=("nCells"), attrs={'long_name':'MPAS-Ocean open ocean mask'}) + + # prepare file to be remapped + out_data_vars = xr.merge([mpasoDomainMaskDA, mpasoOpenOceanMaskDA]) + dataOut = xr.Dataset(data_vars=out_data_vars) + mpaso_mask_file = 'mpaso_mask_file.nc' + dataOut.to_netcdf(mpaso_mask_file, mode='w') + + print("Calling ncremap for masks...") + # remap the input data + args_remap = ["ncremap", + "-i", mpaso_mask_file, + "-o", "mpaso_masks_on_mali_mesh.nc", + "-m", self.mapping_file_name] + check_call(args_remap) def get_data(self, year): @@ -261,6 +280,9 @@ def create_mapping_file(self): check_call(args_esmf) def remap_mpas_to_mali(self, year): + + out_name = f'{self.options.outputFile}_{year}.nc' + print("Start remapping ... ") # populate tmp_mpasoMeshFile with variables to be interpolated @@ -362,16 +384,7 @@ def remap_mpas_to_mali(self, year): # Reshape interpTF back to (nt, nc, nz) interpTF = interpTF.reshape(nt, nc, -1) - # Create mask indentifying valid overlapping ocean cells. Combine all MALI terms in one input file - #Making this a 3-D variable for now, but may want to make 2-D eventually - validOpenOceanMask = np.zeros((nt,nc), dtype='int32') - print("interpTF: {}".format(interpTF.data.shape)) - surfaceTF = interpTF[:,:,0] - ind = np.where(surfaceTF != 0) - validOpenOceanMask[ind] = 1 - - - print("Saving") + print(f"Saving data to {out_name}") #open original mali file IM = xr.open_dataset(self.options.maliFile,decode_times=False,decode_cf=False) @@ -383,7 +396,6 @@ def remap_mpas_to_mali(self, year): ds_out = ds_out.expand_dims(nISMIP6OceanLayers=ismip6shelfMelt_zOcean) # introduce new ismip6 depth dimension # Save variables - mask = xr.DataArray(validOpenOceanMask,dims=("Time","nCells"),attrs={'long_name':'Mask of MALI cells overlapping interpolated MPAS-Oceans cells'}) interptf = xr.DataArray(interpTF, dims=("Time","nCells","nISMIP6OceanLayers")) zlayers = xr.DataArray(ismip6shelfMelt_zOcean, dims=("nISMIP6OceanLayers")) @@ -425,7 +437,6 @@ def remap_mpas_to_mali(self, year): ds_out['zEdge'] = IM['zEdge'] ds_out['zVertex'] = IM['zVertex'] - ds_out['validOceanMask'] = mask ds_out['ismip6shelfMelt_3dThermalForcing'] = interptf ds_out['ismip6shelfMelt_zOcean'] = zlayers @@ -445,7 +456,6 @@ def remap_mpas_to_mali(self, year): if 'history' in ds_out.attrs: del ds_out.attrs['history'] - out_name = f'{self.options.outputFile}_{year}.nc' ds_out.to_netcdf(out_name, mode='w', unlimited_dims=['Time']) ds_out.close() @@ -462,6 +472,8 @@ def main(): run.create_mapping_file() + run.prepare_mpaso_mesh_data() + for yr in range(run.options.startYr, run.options.endYr+1): print(f'### Processing year {yr}') From 74fe08cce50378870830ea13098cd4b6cddfcaa9 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 29 Jul 2024 18:26:02 -0500 Subject: [PATCH 26/39] Move vert. interp. before horiz. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit We want to do ver. before horiz interp to avoid any horiz/vert "mixing" due to potential mpas-ocean hybrid coordinate. Also, this should speed up ncremap because the number of ismip6 vert levels is less than most ocean meshes. After this change, “scalloping” in the open ocean is gone. As this required substantial refactoring, there are a lot of related changes to file handling and xarray objects. --- .../interpolate_mpasOcean_to_mali.py | 184 ++++++++---------- 1 file changed, 80 insertions(+), 104 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index af851e9f0..f748c4382 100755 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -35,6 +35,12 @@ def __init__(self): self.mapping_file_name = f'mpaso_to_mali_mapping_{self.options.method}.nc' + # Define ocean vertical coordinates for mali grid + # For now, hardcode in depths from ISMIP6 AIS. Should work fine for our purposes + self.ismip6shelfMelt_zOcean = np.array((-30, -90, -150, -210, -270, -330, -390, -450, -510,\ + -570, -630, -690, -750, -810, -870, -930, -990, -1050, -1110, -1170,\ + -1230, -1290, -1350, -1410, -1470, -1530, -1590, -1650, -1710, -1770), dtype='float64') + def prepare_mpaso_mesh_data(self): #open mpas ocean mesh @@ -282,124 +288,68 @@ def create_mapping_file(self): def remap_mpas_to_mali(self, year): out_name = f'{self.options.outputFile}_{year}.nc' + tmp_mpasoSourceFile = "tmp_mpasoSourceFile.nc" + tmp_maliDestFile = "tmp_maliDestFile.nc" print("Start remapping ... ") - # populate tmp_mpasoMeshFile with variables to be interpolated + # Get data needed for remapping f = xr.open_dataset(self.options.mpasoMeshFile, decode_times=False, decode_cf=False) tf = xr.DataArray(self.oceanThermalForcing.astype('float64'),dims=('Time','nCells','nVertLevels')) - - f['ismip6shelfMelt_3dThermalForcing'] = tf - mcce = xr.DataArray(self.newMpasCCE.astype('float64'),dims=('Time','nCells','nVertLevels')) - f['mpas_cellCenterElev'] = mcce - - if self.have_landIceFreshwaterFlux: - melt = xr.DataArray(self.newLandIceFWFlux.astype('float64'),dims=('Time','nCells')) - f['floatingBasalMassBal'] = -1.0 * melt - - #delete unncessary variables with 'string1' dimension or will cause errors. - try: - f = f.drop_vars('simulationStartTime') - except ValueError: - f = f - try: - f = f.drop_vars('forcingGroupNames') - except ValueError: - f = f - - try: - f = f.drop_vars('forcingGroupRestartTimes') - except ValueError: - f = f - + # perform vertical interpolation + # do this before horiz interp to avoid any horiz/vert "mixing" due to + # potential mpas-ocean hybrid coordinate. Also should speed up ncremap + # because the number of ismip6 vert levels is less than most ocean meshes st = time.time() - # save new mesh file - tmp_mpasoMeshFile = "tmp_mpasoMeshFile.nc" - f.to_netcdf(tmp_mpasoMeshFile) - f.close() + vertInterpTF = _vertical_interpolate(self, tf.data, mcce.data) nd = time.time() tm = nd - st - print("saving updates mpas mesh file", tm, "seconds") + print("vertical interpolation:", tm, "seconds") st = time.time() - subprocess.run(["ncatted", "-a", "_FillValue,,d,,", tmp_mpasoMeshFile]) + # create new file of ocean data to be remapped + tf_DA = xr.DataArray(vertInterpTF.astype('float64'), + dims=('Time','nCells','nISMIP6OceanLayers')) + dataOut = xr.Dataset(data_vars={'ismip6shelfMelt_3dThermalForcing': tf_DA}) + if self.have_landIceFreshwaterFlux: + melt = xr.DataArray(self.newLandIceFWFlux.astype('float64'),dims=('Time','nCells')) + dataOut['floatingBasalMassBal'] = -1.0 * melt + dataOut.to_netcdf(tmp_mpasoSourceFile, mode='w') + f.close() + + subprocess.run(["ncatted", "-a", "_FillValue,,d,,", tmp_mpasoSourceFile]) + subprocess.run(["ncpdq", "-O", "-a", "Time,nISMIP6OceanLayers,nCells", tmp_mpasoSourceFile, tmp_mpasoSourceFile]) + nd = time.time() tm = nd - st - print("Removing fill value:", tm, "seconds") - - subprocess.run(["ncpdq", "-O", "-a", "Time,nVertLevels,nCells", tmp_mpasoMeshFile, tmp_mpasoMeshFile]) - subprocess.run(["ncpdq", "-O", "-a", "Time,nVertLevelsP1,nCells", tmp_mpasoMeshFile, tmp_mpasoMeshFile]) - subprocess.run(["ncpdq", "-O", "-a", "Time,nVertInterfaces,nCells", tmp_mpasoMeshFile, tmp_mpasoMeshFile]) + print(f"saved modified mpaso data file for remapping and file prep: {tm} seconds") print("ncremap ...") + st = time.time() # remap the input data args_remap = ["ncremap", - "-i", tmp_mpasoMeshFile, - "-o", "tmp_outputFile.nc", + "-i", tmp_mpasoSourceFile, + "-o", tmp_maliDestFile, "-m", self.mapping_file_name] check_call(args_remap) - - #Vertical interpolation - print("Vertically interpolating onto mali grid") - interpDS = xr.open_dataset("tmp_outputFile.nc", decode_times=False, decode_cf=False) - TF = interpDS['ismip6shelfMelt_3dThermalForcing'][:,:,:].data - mpas_cellCenterElev = interpDS['mpas_cellCenterElev'][:,:,:].data - - #Define vertical coordinates for mali grid - #For now, hardcode in depths from ISMIP6 AIS. Should work fine for our purposes - ismip6shelfMelt_zOcean = np.array((-30, -90, -150, -210, -270, -330, -390, -450, -510,\ - -570, -630, -690, -750, -810, -870, -930, -990, -1050, -1110, -1170,\ - -1230, -1290, -1350, -1410, -1470, -1530, -1590, -1650, -1710, -1770),dtype='float64') - - nz, = ismip6shelfMelt_zOcean.shape - - # Reshape to (nt*nc, nz) before interpolation to avoid looping - TF = np.transpose(TF,(0,2,1)) - mpas_cellCenterElev = np.transpose(mpas_cellCenterElev,(0,2,1)) - nt,nc,_ = mpas_cellCenterElev.shape - print("mpas_cellCenterElev: {}".format(mpas_cellCenterElev.shape)) - print("TF: {}".format(TF.shape)) - print("nz: {}:".format(nz)) - TF = TF.reshape(nt*nc, -1) - mpas_cellCenterElev = mpas_cellCenterElev.reshape(nt*nc, -1) - - # Linear interpolation - st = time.time() - nan_count = 0 - interpTF = np.zeros((nt*nc,nz)) - for i in range(nt*nc): - ind = np.where(~np.isnan(TF[i,:].flatten() * mpas_cellCenterElev[i,:].flatten())) - - ct = 0 - if len(ind[0]) != 0: - interpTF[i,:] = np.array(np.interp(ismip6shelfMelt_zOcean, mpas_cellCenterElev[i,ind].flatten(), TF[i,ind].flatten())) - else: - nan_count = nan_count + 1 + subprocess.run(["ncpdq", "-O", "-a", "Time,nCells,nISMIP6OceanLayers", tmp_maliDestFile, tmp_maliDestFile]) nd = time.time() - tm = nd-st - print("vertical interpolation time: {}",tm, "seconds") - - # Reshape interpTF back to (nt, nc, nz) - interpTF = interpTF.reshape(nt, nc, -1) + tm = nd - st + print(f"ncremap completed: {tm} seconds") print(f"Saving data to {out_name}") - #open original mali file - IM = xr.open_dataset(self.options.maliFile,decode_times=False,decode_cf=False) - - #create new mali forcing file that only desired forcing variables - ds_out = IM.copy(deep=False) - ds_out = ds_out.drop_vars(ds_out.data_vars) - - # introduce new ismip6 depth dimension - ds_out = ds_out.expand_dims(nISMIP6OceanLayers=ismip6shelfMelt_zOcean) # introduce new ismip6 depth dimension - # Save variables - interptf = xr.DataArray(interpTF, dims=("Time","nCells","nISMIP6OceanLayers")) - zlayers = xr.DataArray(ismip6shelfMelt_zOcean, dims=("nISMIP6OceanLayers")) + #create new mali forcing file that only contains desired forcing variables + ds_remapped = xr.open_dataset(tmp_maliDestFile, decode_times=False, decode_cf=False) + ds_out = xr.Dataset(data_vars={'ismip6shelfMelt_3dThermalForcing': ds_remapped['ismip6shelfMelt_3dThermalForcing']}) + ds_out['ismip6shelfMelt_zOcean'] = xr.DataArray(self.ismip6shelfMelt_zOcean, dims=("nISMIP6OceanLayers")) + if self.have_landIceFreshwaterFlux: + ds_out['floatingBasalMassBal'] = ds_remapped['floatingBasalMassBal'] if self.options.includeMeshVars: + IM = xr.open_dataset(self.options.maliFile,decode_times=False,decode_cf=False) ds_out['angleEdge'] = IM['angleEdge'] ds_out['areaCell'] = IM['areaCell'] ds_out['areaTriangle'] = IM['areaTriangle'] @@ -436,13 +386,8 @@ def remap_mpas_to_mali(self, year): ds_out['zCell'] = IM['zCell'] ds_out['zEdge'] = IM['zEdge'] ds_out['zVertex'] = IM['zVertex'] + IM.close() - ds_out['ismip6shelfMelt_3dThermalForcing'] = interptf - ds_out['ismip6shelfMelt_zOcean'] = zlayers - - if self.have_landIceFreshwaterFlux: - ds_out['floatingBasalMassBal'] = interpDS['floatingBasalMassBal'][:,:] - # Save xtime #ds_out = ds_out.expand_dims({'StrLen': self.newXtime}, axis=1) @@ -452,10 +397,6 @@ def remap_mpas_to_mali(self, year): xtime.encoding.update({"char_dim_name":"StrLen"}) ds_out['xtime'] = xtime - #clean history for new file - if 'history' in ds_out.attrs: - del ds_out.attrs['history'] - ds_out.to_netcdf(out_name, mode='w', unlimited_dims=['Time']) ds_out.close() @@ -467,6 +408,36 @@ def remap_mpas_to_mali(self, year): if file.startswith("tmp"): os.remove(file) +def _vertical_interpolate(self, TF, mpas_cellCenterElev): + #Vertical interpolation + print("Vertically interpolating onto standardized z-level grid") + + nz2, = self.ismip6shelfMelt_zOcean.shape + nt,nc,nz1 = mpas_cellCenterElev.shape + + # Reshape to (nt*nc, nz) before interpolation to avoid looping + #print("mpas_cellCenterElev: {}".format(mpas_cellCenterElev.shape)) + #print("TF: {}".format(TF.shape)) + #print(f"MPAS-Ocean nz: {nz1}; MALI ocean nz: {nz2}") + TF = TF.reshape(nt*nc, -1) + mpas_cellCenterElev = mpas_cellCenterElev.reshape(nt*nc, -1) + + # Linear interpolation + nan_count = 0 + vertInterpTF = np.zeros((nt*nc, nz2)) + for i in range(nt*nc): + ind = np.where(~np.isnan(TF[i,:].flatten() * mpas_cellCenterElev[i,:].flatten())) + + if len(ind[0]) != 0: + vertInterpTF[i,:] = np.interp(self.ismip6shelfMelt_zOcean, mpas_cellCenterElev[i,ind].flatten(), TF[i,ind].flatten()) + else: + nan_count = nan_count + 1 + + # Reshape vertInterpTF back to (nt, nc, nz) + vertInterpTF = vertInterpTF.reshape(nt, nc, nz2) + + return vertInterpTF + def main(): run = mpasToMaliInterp() @@ -475,7 +446,7 @@ def main(): run.prepare_mpaso_mesh_data() for yr in range(run.options.startYr, run.options.endYr+1): - print(f'### Processing year {yr}') + st = time.time() run.get_data(yr) @@ -490,7 +461,12 @@ def main(): run.DS.close() - print("Finished.") + + nd = time.time() + tm = nd-st + print(f"\n### Processed year {yr} in {tm} seconds\n") + + print("Finished all years.") if __name__ == "__main__": main() From e085dced22de8c8c8861d5492190fc44cace63e4 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Tue, 30 Jul 2024 10:57:27 -0500 Subject: [PATCH 27/39] Corrections to vertical interpolation Previous results had the vertical coordinate upside down. These changes correct that and are more careful about indexing conventions. --- .../interpolate_mpasOcean_to_mali.py | 26 ++++++++++++++----- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index f748c4382..27a211f32 100755 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -412,15 +412,25 @@ def _vertical_interpolate(self, TF, mpas_cellCenterElev): #Vertical interpolation print("Vertically interpolating onto standardized z-level grid") - nz2, = self.ismip6shelfMelt_zOcean.shape nt,nc,nz1 = mpas_cellCenterElev.shape + nz2, = self.ismip6shelfMelt_zOcean.shape + + # numpy.interp expects the x-coordinate sequence to be increasing + # This is not explicitly enforced, and it's unclear if it is + # necessary, but we will follow that convention. + # So we need to flip the ordering of the zOcean array + z_target = np.flip(self.ismip6shelfMelt_zOcean, 0) + assert(np.all(np.diff(z_target) > 0)) # Reshape to (nt*nc, nz) before interpolation to avoid looping - #print("mpas_cellCenterElev: {}".format(mpas_cellCenterElev.shape)) - #print("TF: {}".format(TF.shape)) - #print(f"MPAS-Ocean nz: {nz1}; MALI ocean nz: {nz2}") - TF = TF.reshape(nt*nc, -1) + # Also, flip the ordering of the z-dimension + # numpy.interp requires the independent variable to be increasing + # but cellCenterElev will be decreasing (it is indexed from ocean + # surface down with positive up) + TF = TF.reshape(nt*nc, -nz1) + TF = np.flip(TF, 1) mpas_cellCenterElev = mpas_cellCenterElev.reshape(nt*nc, -1) + mpas_cellCenterElev = np.flip(mpas_cellCenterElev, 1) # Linear interpolation nan_count = 0 @@ -429,12 +439,14 @@ def _vertical_interpolate(self, TF, mpas_cellCenterElev): ind = np.where(~np.isnan(TF[i,:].flatten() * mpas_cellCenterElev[i,:].flatten())) if len(ind[0]) != 0: - vertInterpTF[i,:] = np.interp(self.ismip6shelfMelt_zOcean, mpas_cellCenterElev[i,ind].flatten(), TF[i,ind].flatten()) + assert(np.all(np.diff(mpas_cellCenterElev[i,ind]) > 0)) # confirm correct ordering + vertInterpTF[i,:] = np.interp(z_target, mpas_cellCenterElev[i,ind].flatten(), TF[i,ind].flatten()) else: nan_count = nan_count + 1 - # Reshape vertInterpTF back to (nt, nc, nz) + # Reshape vertInterpTF back to (nt, nc, nz), flip back vertical coordinate vertInterpTF = vertInterpTF.reshape(nt, nc, nz2) + vertInterpTF = np.flip(vertInterpTF, 2) return vertInterpTF From 6fd2405742189ce32b1ea848d9f1daf5fc01aac9 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Tue, 13 Aug 2024 23:41:04 -0500 Subject: [PATCH 28/39] Correct logic for open ocean and cavity coefficients They had been reversed. Co-authored-by: Alexander Hager --- landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index 27a211f32..262ce0ef2 100755 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -246,13 +246,13 @@ def calc_ocean_thermal_forcing(self): 0.5*gravity*(self.newDens[iTime,iCell,k-1]*self.newLThick[iTime,iCell,k-1] \ + self.newDens[iTime,iCell,k]*self.newLThick[iTime,iCell,k]) - if (self.landIceFloatingMask[iCell] == 1): + if (self.landIceFloatingMask[iCell] == 0): ocn_freezing_temperature[iTime,iCell,:] = self.coeff_0_openOcean + \ self.coeff_S_openOcean * self.newSal[iTime,iCell,:] + self.coeff_p_openOcean * self.newPr[iTime,iCell,:] \ + self.coeff_pS_openOcean * self.newPr[iTime,iCell,:] * self.newSal[iTime,iCell,:] \ + self.coeff_mushy_openOcean * self.newSal[iTime,iCell,:] / (1.0 - self.newSal[iTime,iCell,:] / 1e3) - elif (self.landIceFloatingMask[iCell] == 0): + elif (self.landIceFloatingMask[iCell] == 1): ocn_freezing_temperature[iTime,iCell,:] = self.coeff_0_cavity + \ self.coeff_S_cavity * self.newSal[iTime,iCell,:] + self.coeff_p_cavity * self.newPr[iTime,iCell,:] \ + self.coeff_pS_cavity * self.newPr[iTime,iCell,:] * self.newSal[iTime,iCell,:] \ From 5be2a0ff0c5f37befdcbb6f65e9784f0cdeac2d8 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Fri, 16 Aug 2024 11:11:05 -0500 Subject: [PATCH 29/39] Mark vertical levels outside the source data range as nan This allows locations inside the ice shelf or below the bathymetry to be identified on the MALI mesh. --- landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index 262ce0ef2..bafce8e7e 100755 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -440,7 +440,14 @@ def _vertical_interpolate(self, TF, mpas_cellCenterElev): if len(ind[0]) != 0: assert(np.all(np.diff(mpas_cellCenterElev[i,ind]) > 0)) # confirm correct ordering - vertInterpTF[i,:] = np.interp(z_target, mpas_cellCenterElev[i,ind].flatten(), TF[i,ind].flatten()) + # Note: Setting locations outside the extent of valid ocean data to nan + # so that those locations can be avoided + # This linear interpolation approach means that any vertical levels on the destination + # vertical grid that are not bounded on both sides by valid vertical levels on the + # source grid will be marked as nan. We may want to replace this method with one + # that uses and source data overlapping the range of the destination vertical level, + # as well as doing an area weighted remapping of the overlap. + vertInterpTF[i,:] = np.interp(z_target, mpas_cellCenterElev[i,ind].flatten(), TF[i,ind].flatten(), right=np.nan, left=np.nan) else: nan_count = nan_count + 1 From 77e1f980d38aca626085481a4f875209267025af Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 20 Nov 2024 17:20:34 -0600 Subject: [PATCH 30/39] Fix two bugs in Tf calculation * Use cavity Tf coefficients everywhere because we want TF appropriate for cavities (even if extrapolated from the open ocean) * include landIcePressure in the Tf calculation --- .../interpolate_mpasOcean_to_mali.py | 32 +++++++------------ 1 file changed, 11 insertions(+), 21 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index bafce8e7e..9d7e6415f 100755 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -56,6 +56,7 @@ def prepare_mpaso_mesh_data(self): self.minLevelCell = OM['minLevelCell'][:].values else: self.minLevelCell = self.maxLevelCell * 0 + 1 + self.landIcePressure = OM['landIcePressure'][0,:].values self.nCells = OM.sizes['nCells'] @@ -105,17 +106,11 @@ def get_data(self, year): except KeyError: print("No landIceFreshwaterFlux variable") - self.coeff_0_openOcean = self.DS.attrs['config_open_ocean_freezing_temperature_coeff_0'] - self.coeff_S_openOcean = self.DS.attrs['config_open_ocean_freezing_temperature_coeff_S'] - self.coeff_p_openOcean = self.DS.attrs['config_open_ocean_freezing_temperature_coeff_p'] - self.coeff_pS_openOcean = self.DS.attrs['config_open_ocean_freezing_temperature_coeff_pS'] - az1_liq = self.DS.attrs['config_open_ocean_freezing_temperature_coeff_mushy_az1_liq'] - self.coeff_mushy_openOcean = 1/az1_liq self.coeff_0_cavity = self.DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_0'] self.coeff_S_cavity = self.DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_S'] self.coeff_p_cavity = self.DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_p'] self.coeff_pS_cavity = self.DS.attrs['config_land_ice_cavity_freezing_temperature_coeff_pS'] - self.coeff_mushy_cavity = 0 + self.coeff_mushy_cavity = 0.0 # Define vertical coordinates of mpas output mpas_cellCenterElev = compute_zmid(self.bottomDepth, self.maxLevelCell, avg_layerThickness) @@ -239,24 +234,19 @@ def calc_ocean_thermal_forcing(self): kmax = self.maxLevelCell[iCell].values - 1 self.newPr[iTime,iCell,kmin] = self.newAtmPr[iTime,iCell] + \ + self.landIcePressure[iCell] + \ self.newDens[iTime,iCell,kmin]*gravity*0.5*self.newLThick[iTime,iCell,kmin] for k in np.arange(kmin + 1, kmax): self.newPr[iTime,iCell,k] = self.newPr[iTime,iCell,k-1] + \ - 0.5*gravity*(self.newDens[iTime,iCell,k-1]*self.newLThick[iTime,iCell,k-1] \ - + self.newDens[iTime,iCell,k]*self.newLThick[iTime,iCell,k]) - - if (self.landIceFloatingMask[iCell] == 0): - ocn_freezing_temperature[iTime,iCell,:] = self.coeff_0_openOcean + \ - self.coeff_S_openOcean * self.newSal[iTime,iCell,:] + self.coeff_p_openOcean * self.newPr[iTime,iCell,:] \ - + self.coeff_pS_openOcean * self.newPr[iTime,iCell,:] * self.newSal[iTime,iCell,:] \ - + self.coeff_mushy_openOcean * self.newSal[iTime,iCell,:] / (1.0 - self.newSal[iTime,iCell,:] / 1e3) - - elif (self.landIceFloatingMask[iCell] == 1): - ocn_freezing_temperature[iTime,iCell,:] = self.coeff_0_cavity + \ - self.coeff_S_cavity * self.newSal[iTime,iCell,:] + self.coeff_p_cavity * self.newPr[iTime,iCell,:] \ - + self.coeff_pS_cavity * self.newPr[iTime,iCell,:] * self.newSal[iTime,iCell,:] \ - + self.coeff_mushy_cavity * self.newSal[iTime,iCell,:] / (1.0 - self.newSal[iTime,iCell,:] / 1e3) + 0.5 * gravity * (self.newDens[iTime,iCell,k-1] * self.newLThick[iTime,iCell,k-1] \ + + self.newDens[iTime,iCell,k] * self.newLThick[iTime,iCell,k]) + + ocn_freezing_temperature[iTime,iCell,:] = self.coeff_0_cavity + \ + self.coeff_S_cavity * self.newSal[iTime,iCell,:] + \ + self.coeff_p_cavity * self.newPr[iTime,iCell,:] \ + + self.coeff_pS_cavity * self.newPr[iTime,iCell,:] * self.newSal[iTime,iCell,:] \ + + self.coeff_mushy_cavity * self.newSal[iTime,iCell,:] / (1.0 - self.newSal[iTime,iCell,:] / 1e3) nd = time.time() tm = nd - st print("Ocean thermal forcing loop", tm, "seconds") From 3a8f805dc103bfea709a4414c06ac04c785ff3b9 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 20 Nov 2024 23:02:47 -0600 Subject: [PATCH 31/39] Add Bnds array --- .../mesh_tools_li/interpolate_mpasOcean_to_mali.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index 9d7e6415f..7d703151d 100755 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -41,6 +41,18 @@ def __init__(self): -570, -630, -690, -750, -810, -870, -930, -990, -1050, -1110, -1170,\ -1230, -1290, -1350, -1410, -1470, -1530, -1590, -1650, -1710, -1770), dtype='float64') + # create bnds array + self.ismip6shelfMelt_zBndsOcean = np.zeros((2, len(self.ismip6shelfMelt_zOcean))) + self.ismip6shelfMelt_zBndsOcean[0, 0] = 0.0 + for z in range(1, len(self.ismip6shelfMelt_zOcean)): + self.ismip6shelfMelt_zBndsOcean[0, z] = 0.5 * (self.ismip6shelfMelt_zOcean[z - 1] + + self.ismip6shelfMelt_zOcean[z]) + for z in range(0, len(self.ismip6shelfMelt_zOcean) - 1): + self.ismip6shelfMelt_zBndsOcean[1, z] = 0.5 * (self.ismip6shelfMelt_zOcean[z] + + self.ismip6shelfMelt_zOcean[z + 1]) + self.ismip6shelfMelt_zBndsOcean[1, -1] = self.ismip6shelfMelt_zOcean[-1] + \ + (self.ismip6shelfMelt_zOcean[-1] - self.ismip6shelfMelt_zBndsOcean[0, -1]) + def prepare_mpaso_mesh_data(self): #open mpas ocean mesh @@ -335,6 +347,7 @@ def remap_mpas_to_mali(self, year): ds_remapped = xr.open_dataset(tmp_maliDestFile, decode_times=False, decode_cf=False) ds_out = xr.Dataset(data_vars={'ismip6shelfMelt_3dThermalForcing': ds_remapped['ismip6shelfMelt_3dThermalForcing']}) ds_out['ismip6shelfMelt_zOcean'] = xr.DataArray(self.ismip6shelfMelt_zOcean, dims=("nISMIP6OceanLayers")) + ds_out['ismip6shelfMelt_zBndsOcean'] = xr.DataArray(self.ismip6shelfMelt_zBndsOcean, dims=("TWO", "nISMIP6OceanLayers")) if self.have_landIceFreshwaterFlux: ds_out['floatingBasalMassBal'] = ds_remapped['floatingBasalMassBal'] From adce94745b5deecce9570e798c49d685b6e3030b Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 20 Nov 2024 23:59:19 -0600 Subject: [PATCH 32/39] Generate open ocean mask in orig3dOceanMask format --- .../interpolate_mpasOcean_to_mali.py | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index 7d703151d..e37f01f0b 100755 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -75,7 +75,11 @@ def prepare_mpaso_mesh_data(self): # create needed masks and remap them in separate file mpasoDomainMask = np.ones((self.nCells,), dtype=np.double) mpasoDomainMaskDA = xr.DataArray(mpasoDomainMask, name='mpasoDomainMask', dims=("nCells"), attrs={'long_name':'MPAS-Ocean domain mask'}) - mpasoOpenOceanMaskDA = xr.DataArray(np.logical_not(self.landIceFloatingMask).astype(np.double), name='mpasoOpenOceanMask', dims=("nCells"), attrs={'long_name':'MPAS-Ocean open ocean mask'}) + mask1d = np.logical_not(self.landIceFloatingMask).astype(np.double) + mpasoOpenOceanMaskDA = xr.DataArray(mask1d, + name='mpasoOpenOceanMask', + dims=("nCells",), + attrs={'long_name':'MPAS-Ocean open ocean mask'}) # prepare file to be remapped out_data_vars = xr.merge([mpasoDomainMaskDA, mpasoOpenOceanMaskDA]) @@ -91,6 +95,19 @@ def prepare_mpaso_mesh_data(self): "-m", self.mapping_file_name] check_call(args_remap) + # Now on MALI mesh, clean up open ocean mask so it can be used by MALI extrap code + ds = xr.open_dataset("mpaso_masks_on_mali_mesh.nc", decode_times=False, decode_cf=False) + mask1d = ds.mpasoOpenOceanMask[:].values + mask1d = np.where(mask1d > 0.99, 1.0, 0.0).astype(np.int32) # only keep values close to 1 + mask2d = np.tile(mask1d.reshape(-1, 1), (1, len(self.ismip6shelfMelt_zOcean))) # tile mask across all vertical layers + mask3d = mask2d[np.newaxis, :, :] # add time dimension + mpasoOpenOceanMaskDA = xr.DataArray(mask3d, + name='orig3dOceanMask', + dims=("Time", "nCells", "nISMIP6OceanLayers"), + attrs={'long_name':'MPAS-Ocean open ocean mask'}) + mpasoOpenOceanMaskDA.to_netcdf('orig3dOceanMask_open_ocean_mask.nc', mode='w') + + def get_data(self, year): # open and concatenate diagnostic dataset From 78a7b2f6efdc07bee73dab68481c5603c75fb00d Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Thu, 21 Nov 2024 00:53:46 -0600 Subject: [PATCH 33/39] Add 3d ocean mask generation This is a mask on the MALI mesh with the ISMIP6-ocean vertical coordinate that indicates where valid MPAS-Ocean data exists (including inside cavities). To support this, a notable change is that nan is no longer inserted in the TF field outside of the valid range. Instead data is vertically extrapolated and one should rely on orig3dOceanMask to identify where the TF data can be used. --- .../interpolate_mpasOcean_to_mali.py | 45 +++++++++++++++++-- 1 file changed, 41 insertions(+), 4 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index e37f01f0b..76c2ef6c9 100755 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -81,8 +81,32 @@ def prepare_mpaso_mesh_data(self): dims=("nCells",), attrs={'long_name':'MPAS-Ocean open ocean mask'}) + # create mask of 3d ocean data (including cavities) + # Note: Technically, the 3d ocean mask will vary with time because fluctuations in + # the layerThickess and surface pressure will cause the depths that are valid to + # change slightly over time. Eventually, we may want orig3dOceanMask to vary every + # time step, but for now it is assumed to be constant in time. To support the mask being + # constant in time while the depths of the TF data is not, we have to ensure there is + # valid TF data if any mismatch occurs (i.e. where the mask says valid data exists but + # due to changes in layerThickness, it is outside the depth range of valid TF data). + # This is handled by changing the vertical interpolation function to extrapolate instead + # of insert nan when TF is calculated. But here where we calculate the mask we mark + # regions outside the valid depth range with nan so they can be masked. + layerThickness = OM['layerThickness'] + mpas_cellCenterElev = compute_zmid(self.bottomDepth, self.maxLevelCell, layerThickness) + vertInterpResult = _vertical_interpolate(self, np.zeros(layerThickness.shape), mpas_cellCenterElev.data, markExtrap=True) + # the above commands follow the operations used for actual TF data below + # the result will have zero where there is valid data and nan where there is not + maskOcean3d = np.logical_not(np.isnan(vertInterpResult)).astype(np.double) + maskTmp = np.swapaxes(maskOcean3d[0,:,:], 0, 1) # necessary for ncremap + mpaso3dOceanMaskDA = xr.DataArray(maskTmp, + name='mpaso3dOceanMask', + #dims=("nCells", "nISMIP6OceanLayers"), + dims=("nISMIP6OceanLayers", "nCells"), + attrs={'long_name':'MPAS-Ocean 3d ocean mask'}) + # prepare file to be remapped - out_data_vars = xr.merge([mpasoDomainMaskDA, mpasoOpenOceanMaskDA]) + out_data_vars = xr.merge([mpasoDomainMaskDA, mpasoOpenOceanMaskDA, mpaso3dOceanMaskDA]) dataOut = xr.Dataset(data_vars=out_data_vars) mpaso_mask_file = 'mpaso_mask_file.nc' dataOut.to_netcdf(mpaso_mask_file, mode='w') @@ -95,7 +119,7 @@ def prepare_mpaso_mesh_data(self): "-m", self.mapping_file_name] check_call(args_remap) - # Now on MALI mesh, clean up open ocean mask so it can be used by MALI extrap code + # Now on MALI mesh, clean up open ocean masks so they can be used by MALI extrap code ds = xr.open_dataset("mpaso_masks_on_mali_mesh.nc", decode_times=False, decode_cf=False) mask1d = ds.mpasoOpenOceanMask[:].values mask1d = np.where(mask1d > 0.99, 1.0, 0.0).astype(np.int32) # only keep values close to 1 @@ -107,6 +131,16 @@ def prepare_mpaso_mesh_data(self): attrs={'long_name':'MPAS-Ocean open ocean mask'}) mpasoOpenOceanMaskDA.to_netcdf('orig3dOceanMask_open_ocean_mask.nc', mode='w') + mask = ds.mpaso3dOceanMask[:].values + mask = np.swapaxes(mask, 0, 1) # undo swap prior to remapping + mask = np.where(mask > 0.99, 1.0, 0.0).astype(np.int32) # only keep values close to 1 + mask3d = mask[np.newaxis, :, :] # add time dimension + mpaso3dOceanMaskDA = xr.DataArray(mask3d, + name='orig3dOceanMask', + dims=("Time", "nCells", "nISMIP6OceanLayers"), + attrs={'long_name':'MPAS-Ocean full 3d ocean mask'}) + mpaso3dOceanMaskDA.to_netcdf('orig3dOceanMask_3d_ocean_mask_with_cavities.nc', mode='w') + def get_data(self, year): @@ -428,7 +462,7 @@ def remap_mpas_to_mali(self, year): if file.startswith("tmp"): os.remove(file) -def _vertical_interpolate(self, TF, mpas_cellCenterElev): +def _vertical_interpolate(self, TF, mpas_cellCenterElev, markExtrap=False): #Vertical interpolation print("Vertically interpolating onto standardized z-level grid") @@ -467,7 +501,10 @@ def _vertical_interpolate(self, TF, mpas_cellCenterElev): # source grid will be marked as nan. We may want to replace this method with one # that uses and source data overlapping the range of the destination vertical level, # as well as doing an area weighted remapping of the overlap. - vertInterpTF[i,:] = np.interp(z_target, mpas_cellCenterElev[i,ind].flatten(), TF[i,ind].flatten(), right=np.nan, left=np.nan) + if markExtrap: + vertInterpTF[i,:] = np.interp(z_target, mpas_cellCenterElev[i,ind].flatten(), TF[i,ind].flatten(), right=np.nan, left=np.nan) + else: + vertInterpTF[i,:] = np.interp(z_target, mpas_cellCenterElev[i,ind].flatten(), TF[i,ind].flatten()) else: nan_count = nan_count + 1 From 3520ddea526272280499fd430ed02f767d08e038 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Tue, 26 Nov 2024 12:56:24 -0600 Subject: [PATCH 34/39] Correct indexing for ismip6shelfMelt_zBndsOcean --- .../mesh_tools_li/interpolate_mpasOcean_to_mali.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index 76c2ef6c9..fad6e5f88 100755 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -42,16 +42,16 @@ def __init__(self): -1230, -1290, -1350, -1410, -1470, -1530, -1590, -1650, -1710, -1770), dtype='float64') # create bnds array - self.ismip6shelfMelt_zBndsOcean = np.zeros((2, len(self.ismip6shelfMelt_zOcean))) + self.ismip6shelfMelt_zBndsOcean = np.zeros((len(self.ismip6shelfMelt_zOcean), 2)) self.ismip6shelfMelt_zBndsOcean[0, 0] = 0.0 for z in range(1, len(self.ismip6shelfMelt_zOcean)): - self.ismip6shelfMelt_zBndsOcean[0, z] = 0.5 * (self.ismip6shelfMelt_zOcean[z - 1] + + self.ismip6shelfMelt_zBndsOcean[z, 0] = 0.5 * (self.ismip6shelfMelt_zOcean[z - 1] + self.ismip6shelfMelt_zOcean[z]) for z in range(0, len(self.ismip6shelfMelt_zOcean) - 1): - self.ismip6shelfMelt_zBndsOcean[1, z] = 0.5 * (self.ismip6shelfMelt_zOcean[z] + + self.ismip6shelfMelt_zBndsOcean[z, 1] = 0.5 * (self.ismip6shelfMelt_zOcean[z] + self.ismip6shelfMelt_zOcean[z + 1]) - self.ismip6shelfMelt_zBndsOcean[1, -1] = self.ismip6shelfMelt_zOcean[-1] + \ - (self.ismip6shelfMelt_zOcean[-1] - self.ismip6shelfMelt_zBndsOcean[0, -1]) + self.ismip6shelfMelt_zBndsOcean[-1, 1] = self.ismip6shelfMelt_zOcean[-1] + \ + (self.ismip6shelfMelt_zOcean[-1] - self.ismip6shelfMelt_zBndsOcean[-1, 0]) def prepare_mpaso_mesh_data(self): @@ -398,7 +398,7 @@ def remap_mpas_to_mali(self, year): ds_remapped = xr.open_dataset(tmp_maliDestFile, decode_times=False, decode_cf=False) ds_out = xr.Dataset(data_vars={'ismip6shelfMelt_3dThermalForcing': ds_remapped['ismip6shelfMelt_3dThermalForcing']}) ds_out['ismip6shelfMelt_zOcean'] = xr.DataArray(self.ismip6shelfMelt_zOcean, dims=("nISMIP6OceanLayers")) - ds_out['ismip6shelfMelt_zBndsOcean'] = xr.DataArray(self.ismip6shelfMelt_zBndsOcean, dims=("TWO", "nISMIP6OceanLayers")) + ds_out['ismip6shelfMelt_zBndsOcean'] = xr.DataArray(self.ismip6shelfMelt_zBndsOcean, dims=("nISMIP6OceanLayers", "TWO")) if self.have_landIceFreshwaterFlux: ds_out['floatingBasalMassBal'] = ds_remapped['floatingBasalMassBal'] From 93650c626a0a5a4290b2e6ad1a6f43229a16e980 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 27 Nov 2024 14:56:35 -0600 Subject: [PATCH 35/39] Simplify time-averaging function Now that the year loop has been moved to the main driver, the time averaging function can be simplified to only handle the single case of time-averaging one year. The year loop had to be moved to the main driver to avoid running out of memory when processing many years. Note there still is code in the function for the case of not time-averaging. The non-time-averaged case cannot currently be handled by the script, and this is left as a placeholder until later when the driver is restructured to support non-time-averaging. --- .../interpolate_mpasOcean_to_mali.py | 77 +++++-------------- 1 file changed, 19 insertions(+), 58 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index fad6e5f88..adb0adcad 100755 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -181,25 +181,16 @@ def get_data(self, year): def time_average_output(self): print("Time Averaging ...") - if (self.options.yearlyAvg == 'true'): - - yearsSinceStart = self.daysSinceStart.values / 365.0 - finalYear = np.floor(np.max(yearsSinceStart)) - startYear = np.floor(np.min(yearsSinceStart)) - timeStride = 1 # 1 year average - print(f'startYear={startYear}, finalYear={finalYear}') - - if (startYear != finalYear): - years = np.arange(startYear, finalYear, timeStride, dtype=int) - nt = len(years) - else : - years = startYear - nt = 1 + #prepare time vector + yearsSinceStart = self.daysSinceStart.values / 365.0 + stTime = np.datetime64(self.stTime[0:4]) + stYear = np.datetime_as_string(stTime,unit='Y').astype('float64') - print(f'years={years}', nt) + if (self.options.yearlyAvg == 'true'): #pre-allocate + nt = 1 _,nc,nz = self.temperature.shape self.newTemp = np.zeros((nt,nc,nz)) self.newSal = np.zeros((nt,nc,nz)) @@ -210,52 +201,21 @@ def time_average_output(self): if self.have_landIceFreshwaterFlux: self.newLandIceFWFlux = np.zeros((nt,nc)) yearVec = np.zeros((nt,)) - #prepare time vector - stTime = np.datetime64(self.stTime[0:4]) - print("stTime = {}".format(stTime)) - stYear = np.datetime_as_string(stTime,unit='Y').astype('float64') - print("stYear = {}".format(stYear)) - - print("starting loop ...") st = time.time() - if (years.ndim == 0): - log = np.logical_and(yearsSinceStart >= years, yearsSinceStart < years + timeStride) - #ind1 = np.nonzero(yearsSinceStart >= years)[0][0] - #ind2 = np.nonzero(yearsSinceStart < years + timeStride)[0][-1] - #print(ind1, ind2) - #self.newTemp[0,:,:] = np.mean(self.temperature[ind1:ind2+1,:,:], axis=0) - #self.newTemp[0,:,:] = np.mean(self.temperature[log,:,:], axis=0) - self.newTemp[0,:,:] = np.mean(self.temperature[log,:,:].values, axis=0) - self.newSal[0,:,:] = np.mean(self.salinity[log,:,:].values, axis=0) - self.newDens[0,:,:] = np.mean(self.density[log,:,:].values, axis=0) - self.newLThick[0,:,:] = np.mean(self.layerThickness[log,:,:], axis=0) - self.newMpasCCE[0,:,:] = np.mean(self.mpas_cellCenterElev[log,:,:], axis=0) - self.newAtmPr[0,:] = np.mean(self.atmPressure[log,:].values, axis=0) - if self.have_landIceFreshwaterFlux: - self.newLandIceFWFlux[0,:] = np.mean(self.landIceFreshwaterFlux[log,:].values, axis=0) + + self.newTemp[0,:,:] = np.mean(self.temperature.values, axis=0) + self.newSal[0,:,:] = np.mean(self.salinity.values, axis=0) + self.newDens[0,:,:] = np.mean(self.density.values, axis=0) + self.newLThick[0,:,:] = np.mean(self.layerThickness, axis=0) + self.newMpasCCE[0,:,:] = np.mean(self.mpas_cellCenterElev, axis=0) + self.newAtmPr[0,:] = np.mean(self.atmPressure.values, axis=0) + if self.have_landIceFreshwaterFlux: + self.newLandIceFWFlux[0,:] = np.mean(self.landIceFreshwaterFlux.values, axis=0) - #Define time at the first of each year - - print("Test A: {}".format(np.floor(np.min(yearsSinceStart[log])) + stYear)) - yearVec[0] = np.floor(np.min(yearsSinceStart[log])) + stYear - print("yearVec = {}".format(yearVec)) - else : - ct = 0 - for i in years: - log = np.logical_and(yearsSinceStart >= years[i], yearsSinceStart < years[i] + timeStride) - self.newTemp[ct,:,:] = np.mean(self.temperature[log,:,:], axis=0) - self.newSal[ct,:,:] = np.mean(self.salinity[log,:,:], axis=0) - self.newDens[ct,:,:] = np.mean(self.density[log,:,:], axis=0) - self.newLThick[ct,:,:] = np.mean(self.layerThickness[log,:,:], axis=0) - self.newMpasCCE[ct,:,:] = np.mean(self.mpas_cellCenterElev[log,:,:], axis=0) - self.newAtmPr[ct,:] = np.mean(self.atmPressure[log,:], axis=0) - if self.have_landIceFreshwaterFlux: - self.newLandIceFWFlux[ct,:] = np.mean(self.landIceFreshwaterFlux[log,:], axis=0) - - yearVec[ct] = np.floor(np.min(yearsSinceStart[log])) + stYear - print("yearVec = {}".format(yearVec)) - ct = ct + 1 + #Define time at the first of each year + yearVec[0] = np.floor(np.min(yearsSinceStart)) + stYear + nd = time.time() tm = nd - st print("Time averaging loop", tm, "seconds") @@ -523,6 +483,7 @@ def main(): for yr in range(run.options.startYr, run.options.endYr+1): st = time.time() + print(f'\n**** Processing year {yr} ****\n') run.get_data(yr) From c887b219d02549443eaf77bf6da7b186228ba346 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 27 Nov 2024 10:57:28 -0600 Subject: [PATCH 36/39] Change masking threshold and marking of bad data This commit changes the threshold for overlap between an MPAS-Ocean grid cell and a MALI grid cell for retaining interpolated values from 0.99 to 1.0. It also fills all locations that don't have valid data with nan instead of extrapolating, making it easier to see what was intentionally interpolated and what wasn't. --- landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index adb0adcad..ab5ba4132 100755 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -122,7 +122,7 @@ def prepare_mpaso_mesh_data(self): # Now on MALI mesh, clean up open ocean masks so they can be used by MALI extrap code ds = xr.open_dataset("mpaso_masks_on_mali_mesh.nc", decode_times=False, decode_cf=False) mask1d = ds.mpasoOpenOceanMask[:].values - mask1d = np.where(mask1d > 0.99, 1.0, 0.0).astype(np.int32) # only keep values close to 1 + mask1d = np.where(mask1d >= 1.0, 1.0, 0.0).astype(np.int32) # only keep values of 1 mask2d = np.tile(mask1d.reshape(-1, 1), (1, len(self.ismip6shelfMelt_zOcean))) # tile mask across all vertical layers mask3d = mask2d[np.newaxis, :, :] # add time dimension mpasoOpenOceanMaskDA = xr.DataArray(mask3d, @@ -133,7 +133,7 @@ def prepare_mpaso_mesh_data(self): mask = ds.mpaso3dOceanMask[:].values mask = np.swapaxes(mask, 0, 1) # undo swap prior to remapping - mask = np.where(mask > 0.99, 1.0, 0.0).astype(np.int32) # only keep values close to 1 + mask = np.where(mask >= 1.0, 1.0, 0.0).astype(np.int32) # only keep values of 1 mask3d = mask[np.newaxis, :, :] # add time dimension mpaso3dOceanMaskDA = xr.DataArray(mask3d, name='orig3dOceanMask', @@ -422,7 +422,7 @@ def remap_mpas_to_mali(self, year): if file.startswith("tmp"): os.remove(file) -def _vertical_interpolate(self, TF, mpas_cellCenterElev, markExtrap=False): +def _vertical_interpolate(self, TF, mpas_cellCenterElev, markExtrap=True): #Vertical interpolation print("Vertically interpolating onto standardized z-level grid") From fbe08a5e5bc8ac2491edb5b4511cea64dc70c734 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 27 Nov 2024 12:27:12 -0600 Subject: [PATCH 37/39] Refactor TF function With vectorization, this function goes from 243.6 to 2.5 seconds for a single call on the v2.1 SORRM mesh. I confirmed the output looks the same as before but speckling in TF is now gone. --- .../interpolate_mpasOcean_to_mali.py | 53 ++++++++++--------- 1 file changed, 29 insertions(+), 24 deletions(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index ab5ba4132..87497cc2f 100755 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -240,41 +240,46 @@ def time_average_output(self): print("New xtime: {}".format(self.newXtime)) + def calc_ocean_thermal_forcing(self): print("Calculating thermal forcing ... ") + st = time.time() gravity = 9.81 #pre-allocate - nt,nc,nz = self.newTemp.shape - ocn_freezing_temperature = np.zeros((nt,nc,nz)) - self.newPr = np.zeros((nt,nc,nz)) - st = time.time() + nt, nc, nz = self.newTemp.shape + self.oceanThermalForcing = np.zeros((nt, nc, nz)) for iTime in range(nt): - #calculate pressure: - for iCell in range(self.nCells): - kmin = self.minLevelCell[iCell] - 1 - kmax = self.maxLevelCell[iCell].values - 1 - - self.newPr[iTime,iCell,kmin] = self.newAtmPr[iTime,iCell] + \ - self.landIcePressure[iCell] + \ - self.newDens[iTime,iCell,kmin]*gravity*0.5*self.newLThick[iTime,iCell,kmin] - - for k in np.arange(kmin + 1, kmax): - self.newPr[iTime,iCell,k] = self.newPr[iTime,iCell,k-1] + \ - 0.5 * gravity * (self.newDens[iTime,iCell,k-1] * self.newLThick[iTime,iCell,k-1] \ - + self.newDens[iTime,iCell,k] * self.newLThick[iTime,iCell,k]) - - ocn_freezing_temperature[iTime,iCell,:] = self.coeff_0_cavity + \ - self.coeff_S_cavity * self.newSal[iTime,iCell,:] + \ - self.coeff_p_cavity * self.newPr[iTime,iCell,:] \ - + self.coeff_pS_cavity * self.newPr[iTime,iCell,:] * self.newSal[iTime,iCell,:] \ - + self.coeff_mushy_cavity * self.newSal[iTime,iCell,:] / (1.0 - self.newSal[iTime,iCell,:] / 1e3) + minLev = self.minLevelCell - 1 + maxLev = self.maxLevelCell.values - 1 + levMask = np.zeros((nc, nz)) + for k in np.arange(nz): + levMask[np.logical_and(k >= minLev, k < maxLev), k] = 1.0 + + newPr = np.zeros((nc,nz)) + newPr[:, 0] = self.newAtmPr[iTime, :] + self.landIcePressure[:] + \ + levMask[:, 0] * 0.5 * gravity * self.newDens[iTime, :, 0] * self.newLThick[iTime, :, 0] + for k in np.arange(1, nz): + newPr[:, k] = newPr[:, k - 1] + \ + levMask[:, k - 1] * 0.5 * gravity * self.newDens[iTime, :, k - 1] * self.newLThick[iTime, :, k - 1] \ + + levMask[:, k] * 0.5 * gravity * self.newDens[iTime, :, k] * self.newLThick[iTime, :, k] + + ocn_freezing_temperature = np.zeros((nc,nz)) + for k in np.arange(nz): + ocn_freezing_temperature[:, k] = self.coeff_0_cavity + \ + self.coeff_S_cavity * self.newSal[iTime, :, k] + \ + self.coeff_p_cavity * newPr[:, k] \ + + self.coeff_pS_cavity * newPr[:, k] * self.newSal[iTime, :, k] \ + + self.coeff_mushy_cavity * self.newSal[iTime, :, k] / (1.0 - self.newSal[iTime, :, k] / 1.e3) + + self.oceanThermalForcing[iTime, :, :] = self.newTemp[iTime, :, :] - ocn_freezing_temperature[:, :] + nd = time.time() tm = nd - st print("Ocean thermal forcing loop", tm, "seconds") # Calculate thermal forcing - self.oceanThermalForcing = self.newTemp - ocn_freezing_temperature + def create_mapping_file(self): From 672a159d2d8b16b689e1fb276effed6c326d2b18 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 27 Nov 2024 15:53:12 -0600 Subject: [PATCH 38/39] Fix indexing for maxLevelCell --- landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index 87497cc2f..7a80b04de 100755 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -255,7 +255,7 @@ def calc_ocean_thermal_forcing(self): maxLev = self.maxLevelCell.values - 1 levMask = np.zeros((nc, nz)) for k in np.arange(nz): - levMask[np.logical_and(k >= minLev, k < maxLev), k] = 1.0 + levMask[np.logical_and(k >= minLev, k <= maxLev), k] = 1.0 newPr = np.zeros((nc,nz)) newPr[:, 0] = self.newAtmPr[iTime, :] + self.landIcePressure[:] + \ From b5415281912f671173d0e18ffb48e4802cb7052b Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Tue, 21 Jan 2025 21:11:10 -0600 Subject: [PATCH 39/39] Allow vertical extrapolation This avoids issues where small changes in vertical coordinate cause nans to occur in locations that the reference mask indicates are valid --- landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py index 7a80b04de..5a94c6fca 100755 --- a/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py +++ b/landice/mesh_tools_li/interpolate_mpasOcean_to_mali.py @@ -321,7 +321,10 @@ def remap_mpas_to_mali(self, year): # potential mpas-ocean hybrid coordinate. Also should speed up ncremap # because the number of ismip6 vert levels is less than most ocean meshes st = time.time() - vertInterpTF = _vertical_interpolate(self, tf.data, mcce.data) + # allow vertical extrapolation. + # this avoids issues where small changes in vertical coordinate cause nans to + # occur in locations that the reference mask indicates are valid + vertInterpTF = _vertical_interpolate(self, tf.data, mcce.data, markExtrap=False) nd = time.time() tm = nd - st print("vertical interpolation:", tm, "seconds")