From 4f05c93fe5962fb1a16f28ec5944554c1bf12d75 Mon Sep 17 00:00:00 2001 From: domingom Date: Wed, 19 Mar 2025 15:45:10 -0600 Subject: [PATCH 1/7] Small change to the title units of a plot --- scripts/python_utilities/coupler/GeoSpec.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/python_utilities/coupler/GeoSpec.py b/scripts/python_utilities/coupler/GeoSpec.py index 44172b5..d8e5406 100644 --- a/scripts/python_utilities/coupler/GeoSpec.py +++ b/scripts/python_utilities/coupler/GeoSpec.py @@ -194,7 +194,7 @@ ax.set_aspect('equal', 'box') ax.set_ylabel(r'$y$ $[\mathrm{km}]$',fontsize=fntSize_labels) ax.set_xlabel(r'$x$ $[\mathrm{km}]$',fontsize=fntSize_labels) - ax.set_title('terrain elevation [m]',fontsize=fntSize_title) + ax.set_title('terrain elevation [m ASL]',fontsize=fntSize_title) ### land cover ### From 39ce90c8acbb960f8b6565f616571ec5be9fd68a Mon Sep 17 00:00:00 2001 From: domingom Date: Tue, 1 Apr 2025 10:49:35 -0600 Subject: [PATCH 2/7] Modifications to GenICBCs.py to handle second increments. --- scripts/python_utilities/coupler/GenICBCs.py | 97 +++++++------------ .../python_utilities/coupler/genicbcs.json | 5 +- 2 files changed, 37 insertions(+), 65 deletions(-) diff --git a/scripts/python_utilities/coupler/GenICBCs.py b/scripts/python_utilities/coupler/GenICBCs.py index b54fd8c..c571642 100644 --- a/scripts/python_utilities/coupler/GenICBCs.py +++ b/scripts/python_utilities/coupler/GenICBCs.py @@ -39,23 +39,23 @@ with open(args.file) as file: params = json.loads(file.read()) -if True: - ICBC_dir = params["ICBC_dir"] - FE_simGrid = params["FE_simGrid"] - WRF_PrntDir = params["WRF_PrntDir"] - WRF_PrntRefOut = params["WRF_PrntRefOut"] - WRF_PrntOutPrefix = params["WRF_PrntOutPrefix"] - dateString = params["dateString"] - timeHour0 = params["timeHour0"] - timeMinute0 = params["timeMinute0"] - itMax = params["itMax"] - itInc = params["itInc"] - - print(f"{mpi_rank}/{mpi_size}: Writing coupler outputs to {ICBC_dir}") - print(f"{mpi_rank}/{mpi_size}: Interpolating to FE-domain from {FE_simGrid}") - print(f"{mpi_rank}/{mpi_size}: Using WRF-reference from {WRF_PrntDir}{WRF_PrntRefOut}") - print(f"{mpi_rank}/{mpi_size}: Processing of WRF-files {WRF_PrntDir}{WRF_PrntOutPrefix}*") - print(f"{mpi_rank}/{mpi_size}: Date and times of WRF-files to process: {dateString}_{timeHour0:02}:{timeMinute0:02}:*, every {itInc} minutes for {itMax} total minutes.") +ICBC_dir = params["ICBC_dir"] +FE_simGrid = params["FE_simGrid"] +WRF_PrntDir = params["WRF_PrntDir"] +WRF_PrntRefOut = params["WRF_PrntRefOut"] +WRF_PrntOutPrefix = params["WRF_PrntOutPrefix"] +dateString = params["dateString"] +timeHour0 = params["timeHour0"] +timeMinute0 = params["timeMinute0"] +timeSecond0 = params["timeSecond0"] +secMax = params["secMax"] +secInc = params["secInc"] + +print(f"{mpi_rank}/{mpi_size}: Writing coupler outputs to {ICBC_dir}") +print(f"{mpi_rank}/{mpi_size}: Interpolating to FE-domain from {FE_simGrid}") +print(f"{mpi_rank}/{mpi_size}: Using WRF-reference from {WRF_PrntDir}{WRF_PrntRefOut}") +print(f"{mpi_rank}/{mpi_size}: Processing of WRF-files {WRF_PrntDir}{WRF_PrntOutPrefix}*") +print(f"{mpi_rank}/{mpi_size}: Date and times of WRF-files to process: {dateString}_{timeHour0:02}:{timeMinute0:02}:*, every {secInc} s for {secMax} total seconds.") ################################################################################################ ### Create a coupler output directory for initial and boundary conditions if necessary @@ -69,48 +69,20 @@ ################################################################################################ files_list=[] times=[] -if False: - timeHour = timeHour0 - timeMinute = timeMinute0 - hourIncFlag = True - for it in range(0,itMax,itInc): - if (it > 0) and ((it+timeMinute0)%60 < itInc) and hourIncFlag : - timeHour+=1 - hourIncFlag = False - else: - hourIncFlag = True - thistime = "{:s}{:02d}:{:02d}:{:s}".format(dateString,timeHour,(timeMinute+it)%60,"00") - times.append(thistime) - for eachtime in times: - file_tmp = f'{WRF_PrntDir}/{WRF_PrntOutPrefix}{eachtime}' - files_list.append(file_tmp) - if(mpi_rank == 0): - print(file_tmp) - #print(files_list) - -if True: # DME - year0 = int(dateString[0:4]) - month0 = int(dateString[5:7]) - day0 = int(dateString[8:10]) - - date_it = dt.datetime(year0,month0,day0,timeHour0,timeMinute0) - for it in range(0,itMax,itInc): - dateString_it = str(date_it.year) + '-' + "{:02d}".format(date_it.month) + '-' + "{:02d}".format(date_it.day) + '_' - thistime = "{:s}{:02d}:{:02d}:{:s}".format(dateString_it,date_it.hour,date_it.minute,"00") - file_tmp = f'{WRF_PrntDir}/{WRF_PrntOutPrefix}{thistime}' - files_list.append(file_tmp) - if(mpi_rank == 0): - print(file_tmp) - date_it = date_it + dt.timedelta(minutes=itInc) - -it00=0 -for Bdy_file_num in range(it00,len(files_list)): - if Bdy_file_num%(60/itInc) == 0: - timeHour=int(Bdy_file_num/(60/itInc))+timeHour0 - timeLabel="{:02d}{:02d}UTC".format(timeHour,timeMinute0) - print('timeLabel for IC file:',timeLabel) - if(mpi_rank == 0): - print("{:s} at Bdy_file_num = {:d}".format(timeLabel,Bdy_file_num)) + +year0 = int(dateString[0:4]) +month0 = int(dateString[5:7]) +day0 = int(dateString[8:10]) + +date_it = dt.datetime(year0,month0,day0,timeHour0,timeMinute0,timeSecond0) +for it in range(0,secMax,secInc): + dateString_it = str(date_it.year) + '-' + "{:02d}".format(date_it.month) + '-' + "{:02d}".format(date_it.day) + '_' + thistime = "{:s}{:02d}:{:02d}:{:02d}".format(dateString_it,date_it.hour,date_it.minute,date_it.second) + file_tmp = f'{WRF_PrntDir}/{WRF_PrntOutPrefix}{thistime}' + files_list.append(file_tmp) + if(mpi_rank == 0): + print(file_tmp) + date_it = date_it + dt.timedelta(seconds=secInc) ################################################################################################ ### Setup mpi task decomposition over the set of files to process @@ -383,10 +355,9 @@ verticalInterpFinal(ds_FEGrid,dsFENew,dsFEFinal,zRect) t3e = time.perf_counter() print('{:d}/{:d}: t3_elapsed = {:f} (s)'.format(mpi_rank, mpi_size, t3e-t3s)) - addTimeDim_FEfinal(dsFEFinal) # DME - if Bdy_file_num%(60/itInc) == 0: - timeHour=int(Bdy_file_num/(60/itInc))+timeHour0 - timeLabel="{:02d}{:02d}UTC".format(timeHour,timeMinute0) + addTimeDim_FEfinal(dsFEFinal) + if Bdy_file_num == 0: + timeLabel="{:02d}{:02d}{:02d}UTC".format(timeHour0,timeMinute0,timeSecond0) dsFEFinal.to_netcdf(ICBC_dir+'FE_interp_{:s}.{:d}'.format(timeLabel,0),format='NETCDF4', encoding={'xIndex': {'dtype': 'i4'},'yIndex': {'dtype': 'i4'},'zIndex': {'dtype': 'i4'}}) t4s = time.perf_counter() diff --git a/scripts/python_utilities/coupler/genicbcs.json b/scripts/python_utilities/coupler/genicbcs.json index e2b2429..f58bc0f 100644 --- a/scripts/python_utilities/coupler/genicbcs.json +++ b/scripts/python_utilities/coupler/genicbcs.json @@ -7,6 +7,7 @@ "dateString": "2024-03-20", "timeHour0": 1, "timeMinute0": 30, - "itMax": 120, - "itInc": 5 + "timeSecond0": 0, + "secMax": 7200, + "secInc": 300 } From 948407fff3078d3b8e229a79771be002b0b6c055 Mon Sep 17 00:00:00 2001 From: domingom Date: Tue, 1 Apr 2025 15:57:48 -0600 Subject: [PATCH 3/7] Removal of wrfRef file and a small fix to the input GIS file name in geospec.json. --- scripts/python_utilities/coupler/GenICBCs.py | 18 ++++++++---------- ...tadata.csv => LandCoverMetadata_NLCD16.csv} | 0 scripts/python_utilities/coupler/genicbcs.json | 3 +-- scripts/python_utilities/coupler/geospec.json | 4 ++-- 4 files changed, 11 insertions(+), 14 deletions(-) rename scripts/python_utilities/coupler/{LandCoverMetadata.csv => LandCoverMetadata_NLCD16.csv} (100%) diff --git a/scripts/python_utilities/coupler/GenICBCs.py b/scripts/python_utilities/coupler/GenICBCs.py index c571642..884a53f 100644 --- a/scripts/python_utilities/coupler/GenICBCs.py +++ b/scripts/python_utilities/coupler/GenICBCs.py @@ -42,9 +42,8 @@ ICBC_dir = params["ICBC_dir"] FE_simGrid = params["FE_simGrid"] WRF_PrntDir = params["WRF_PrntDir"] -WRF_PrntRefOut = params["WRF_PrntRefOut"] WRF_PrntOutPrefix = params["WRF_PrntOutPrefix"] -dateString = params["dateString"] +dateString = params["date0"] timeHour0 = params["timeHour0"] timeMinute0 = params["timeMinute0"] timeSecond0 = params["timeSecond0"] @@ -53,7 +52,6 @@ print(f"{mpi_rank}/{mpi_size}: Writing coupler outputs to {ICBC_dir}") print(f"{mpi_rank}/{mpi_size}: Interpolating to FE-domain from {FE_simGrid}") -print(f"{mpi_rank}/{mpi_size}: Using WRF-reference from {WRF_PrntDir}{WRF_PrntRefOut}") print(f"{mpi_rank}/{mpi_size}: Processing of WRF-files {WRF_PrntDir}{WRF_PrntOutPrefix}*") print(f"{mpi_rank}/{mpi_size}: Date and times of WRF-files to process: {dateString}_{timeHour0:02}:{timeMinute0:02}:*, every {secInc} s for {secMax} total seconds.") @@ -78,7 +76,7 @@ for it in range(0,secMax,secInc): dateString_it = str(date_it.year) + '-' + "{:02d}".format(date_it.month) + '-' + "{:02d}".format(date_it.day) + '_' thistime = "{:s}{:02d}:{:02d}:{:02d}".format(dateString_it,date_it.hour,date_it.minute,date_it.second) - file_tmp = f'{WRF_PrntDir}/{WRF_PrntOutPrefix}{thistime}' + file_tmp = f'{WRF_PrntDir}{WRF_PrntOutPrefix}{thistime}' files_list.append(file_tmp) if(mpi_rank == 0): print(file_tmp) @@ -93,7 +91,7 @@ listCntr=0 list_StartOffset=0 while (listCntr < list_len) and currentExists: - bdyFileName = "{:s}/FE_Bndys.{:d}".format(ICBC_dir,listCntr) + bdyFileName = "{:s}FE_Bndys.{:d}".format(ICBC_dir,listCntr) if (mpi_rank ==0): print("{:d}/{:d}: Checking for {:s}".format(mpi_rank, mpi_size,bdyFileName)) currentExists = os.path.isfile(bdyFileName) @@ -138,7 +136,7 @@ ######################################## ### Load the reference WRF data file ### ######################################## -ds_WRFRef=xr.open_dataset(f'{WRF_PrntDir}/{WRF_PrntRefOut}') +ds_WRFRef=xr.open_dataset(files_list[0]) ############################################################################################ ### Locate the FE-grid-bounding corners as index pairs (j,i) in the WRF-reference domain ### @@ -305,8 +303,8 @@ ### Define a rectilinear vertical coordinate basis for vertical interpolation between WRF and FE ### #################################################################################################### zBottom = 0.0 -zTop = ds_FEGrid['zPos'][-1,0,0].values+90.0 #500 -NzWRFInterp = 275 #125 #100 +zTop = ds_FEGrid['zPos'][-1,0,0].values+90.0 +NzWRFInterp = 275 zRect = np.linspace(zBottom,zTop,NzWRFInterp) if(mpi_rank == 0) and DEBUG_COUPLER: print(zRect[0],zRect[-1]) @@ -329,9 +327,9 @@ for Bdy_file_num in range(it00,it11): bdyFileName = "{:s}/FE_Bndys.{:d}".format(ICBC_dir,Bdy_file_num) if not(os.path.isfile(bdyFileName)): - print('{:d}/{:d}: {:s} does not exist, creating it...'.format(mpi_rank, mpi_size, bdyFileName)) + print('{:d}{:d}: {:s} does not exist, creating it...'.format(mpi_rank, mpi_size, bdyFileName)) if True: - print("{:d}/{:d}: Working on file {:s}".format(mpi_rank, mpi_size, files_list[Bdy_file_num])) + print("{:d}{:d}: Working on file {:s}".format(mpi_rank, mpi_size, files_list[Bdy_file_num])) ds_ref = xr.open_mfdataset(files_list[Bdy_file_num],combine='nested',concat_dim='Time') t0s = time.perf_counter() diff --git a/scripts/python_utilities/coupler/LandCoverMetadata.csv b/scripts/python_utilities/coupler/LandCoverMetadata_NLCD16.csv similarity index 100% rename from scripts/python_utilities/coupler/LandCoverMetadata.csv rename to scripts/python_utilities/coupler/LandCoverMetadata_NLCD16.csv diff --git a/scripts/python_utilities/coupler/genicbcs.json b/scripts/python_utilities/coupler/genicbcs.json index f58bc0f..c41c2bf 100644 --- a/scripts/python_utilities/coupler/genicbcs.json +++ b/scripts/python_utilities/coupler/genicbcs.json @@ -2,9 +2,8 @@ "ICBC_dir": "/path_FEsim/ICBC/", "FE_simGrid": "/path_simgrid/simgrid.0", "WRF_PrntDir": "/path_wrf_data/", - "WRF_PrntRefOut": "wrfout_d01_2024-03-20_01:00:00", "WRF_PrntOutPrefix": "wrf_fasteddy_d01_", - "dateString": "2024-03-20", + "date0": "2024-03-20", "timeHour0": 1, "timeMinute0": 30, "timeSecond0": 0, diff --git a/scripts/python_utilities/coupler/geospec.json b/scripts/python_utilities/coupler/geospec.json index cd97386..3615b88 100644 --- a/scripts/python_utilities/coupler/geospec.json +++ b/scripts/python_utilities/coupler/geospec.json @@ -1,8 +1,8 @@ { "name_dom": "DomainName", "gis_root": "/path_gis/", - "gis_file": "gis.nc", - "nlcd_name": "LandCoverMetadata.csv", + "gis_file": "inputs_gis.nc", + "nlcd_name": "LandCoverMetadata_NLCD16.csv", "water_cats": [11], "FE_dataset_path": "/path_geospec/", "name_dom_add": "", From a589ee7371332436c0ff77de32af5bcffc6b4abe Mon Sep 17 00:00:00 2001 From: domingom Date: Tue, 8 Apr 2025 10:33:40 -0600 Subject: [PATCH 4/7] Change input GIS name of terrain from topoPos to elevation --- scripts/python_utilities/coupler/GeoSpec.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/python_utilities/coupler/GeoSpec.py b/scripts/python_utilities/coupler/GeoSpec.py index d8e5406..5f60a96 100644 --- a/scripts/python_utilities/coupler/GeoSpec.py +++ b/scripts/python_utilities/coupler/GeoSpec.py @@ -68,7 +68,7 @@ vars_gis_ref_v = ['lat','lon','data_topo0','data_land'] if (gis_opt==0): - vars_gis_v = ['lat','lon','topoPos','LandCover'] + vars_gis_v = ['lat','lon','elevation','LandCover'] elif (gis_opt==1): vars_gis_v = ['XLAT','XLONG','HGT','LU_INDEX'] From bda2ebc2053f1efe0322b159c93036dfbbe37447 Mon Sep 17 00:00:00 2001 From: domingom Date: Tue, 8 Apr 2025 10:44:51 -0600 Subject: [PATCH 5/7] Add sea mask plot to GeoSpec --- scripts/python_utilities/coupler/GeoSpec.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/scripts/python_utilities/coupler/GeoSpec.py b/scripts/python_utilities/coupler/GeoSpec.py index 5f60a96..93a6bd9 100644 --- a/scripts/python_utilities/coupler/GeoSpec.py +++ b/scripts/python_utilities/coupler/GeoSpec.py @@ -35,7 +35,7 @@ file_nlcd = gis_root + nlcd_name FE_new_nc = FE_dataset_path + name_dom + name_dom_add + '.nc' -FE_plot = FE_dataset_path + name_dom + name_dom_add + '.png' +FE_plot = FE_dataset_path + name_dom + name_dom_add + '_geospec.png' print('FE_new_nc:', FE_new_nc) # Calculate xPos2d, yPos2d @@ -167,7 +167,6 @@ [0.43921569, 0.63921569, 0.72941176]] nlcd_color_map = ListedColormap(_nlcd_colors_1, name='nlcd_color_map') -# plt.register_cmap(name='nlcd_color_map', cmap=nlcd_color_map) matplotlib.colormaps.register(nlcd_color_map, name='nlcd_color_map', force=False) fntSize = 14.0 @@ -221,8 +220,15 @@ ax.set_xlabel(r'$x$ $[\mathrm{km}]$',fontsize=fntSize_labels) ax.set_title('roughness length, z0m [m]',fontsize=fntSize_title) - ## make the 4th panel invisible ### - axs[1,1].set_visible(False) + ## SeaMask ### + ax=axs[1,1] + + im = ax.pcolormesh(xarr/1e3,yarr/1e3,SeaMask_tmp,cmap='bwr_r',linewidth=0,rasterized=True,shading='nearest') + cbar=fig.colorbar(im, ax=ax, orientation='vertical') + ax.set_aspect('equal', 'box') + ax.set_ylabel(r'$y$ $[\mathrm{km}]$',fontsize=fntSize_labels) + ax.set_xlabel(r'$x$ $[\mathrm{km}]$',fontsize=fntSize_labels) + ax.set_title('sea mask [-]',fontsize=fntSize_title) # save figure CCC = FE_plot From d67a48ea740cd88a5ab48f04e2323d16bdb2f40a Mon Sep 17 00:00:00 2001 From: domingom Date: Tue, 8 Apr 2025 10:48:43 -0600 Subject: [PATCH 6/7] Adjust simgrid plot naeme --- scripts/python_utilities/coupler/SimGrid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/python_utilities/coupler/SimGrid.py b/scripts/python_utilities/coupler/SimGrid.py index 464bddf..765bfea 100644 --- a/scripts/python_utilities/coupler/SimGrid.py +++ b/scripts/python_utilities/coupler/SimGrid.py @@ -41,7 +41,7 @@ start_code_all = time.perf_counter() FE_new_nc = FE_new_nc_path + name_dom + name_dom_add + '.0' -FE_plot = FE_new_nc_path + name_dom + name_dom_add + '.png' +FE_plot = FE_new_nc_path + name_dom + name_dom_add + '_simgrid.png' print('FE_ref_GIS_nc: ',FE_ref_GIS_nc) print('FE_new_nc: ',FE_new_nc) From feedd004c8d1e39d8789110d20f964c6f351f084 Mon Sep 17 00:00:00 2001 From: domingom Date: Tue, 8 Apr 2025 10:59:18 -0600 Subject: [PATCH 7/7] Remove if True/Fasle statements and some comments --- scripts/python_utilities/coupler/GenICBCs.py | 36 +++++------- .../python_utilities/coupler/couplingUtils.py | 58 +++---------------- 2 files changed, 23 insertions(+), 71 deletions(-) diff --git a/scripts/python_utilities/coupler/GenICBCs.py b/scripts/python_utilities/coupler/GenICBCs.py index 884a53f..6e5a33a 100644 --- a/scripts/python_utilities/coupler/GenICBCs.py +++ b/scripts/python_utilities/coupler/GenICBCs.py @@ -262,8 +262,7 @@ print(XvWRF.shape,XvWRF.shape) -if True: - print(ds_WRFRef['HGT'][0,ll_jindx:ll_jindx+j_extent,ll_iindx:ll_iindx+i_extent].values) +print(ds_WRFRef['HGT'][0,ll_jindx:ll_jindx+j_extent,ll_iindx:ll_iindx+i_extent].values) #################################################################################### ### Map the target FE domain into the WRF bounding-grid relative x,y coordinates ### @@ -328,25 +327,20 @@ bdyFileName = "{:s}/FE_Bndys.{:d}".format(ICBC_dir,Bdy_file_num) if not(os.path.isfile(bdyFileName)): print('{:d}{:d}: {:s} does not exist, creating it...'.format(mpi_rank, mpi_size, bdyFileName)) - if True: - print("{:d}{:d}: Working on file {:s}".format(mpi_rank, mpi_size, files_list[Bdy_file_num])) - ds_ref = xr.open_mfdataset(files_list[Bdy_file_num],combine='nested',concat_dim='Time') - - t0s = time.perf_counter() - dsWRF=interpWRFToGrids(ds_ref,it0,varsList,surfVarsList,zRect,ll_iindx,i_extent,ll_jindx,j_extent) - t0e = time.perf_counter() - print('{:d}/{:d}: t0_elapsed = {:f} (s)'.format(mpi_rank, mpi_size, t0e-t0s)) - t1s = time.perf_counter() - ds=copyAndTranspose(dsWRF) - t1e = time.perf_counter() - print('{:d}/{:d}: t1_elapsed = {:f} (s)'.format(mpi_rank, mpi_size, t1e-t1s)) - t2s = time.perf_counter() - dsFENew=interp2DForFE(ds,ds_FEGrid,XvWRF,YvWRF,xVec,yVec) - t2e = time.perf_counter() - else: - t2s = time.perf_counter() - dsFENew=xr.open_dataset(FE_simGrid) - t2e = time.perf_counter() + print("{:d}{:d}: Working on file {:s}".format(mpi_rank, mpi_size, files_list[Bdy_file_num])) + ds_ref = xr.open_mfdataset(files_list[Bdy_file_num],combine='nested',concat_dim='Time') + + t0s = time.perf_counter() + dsWRF=interpWRFToGrids(ds_ref,it0,varsList,surfVarsList,zRect,ll_iindx,i_extent,ll_jindx,j_extent) + t0e = time.perf_counter() + print('{:d}/{:d}: t0_elapsed = {:f} (s)'.format(mpi_rank, mpi_size, t0e-t0s)) + t1s = time.perf_counter() + ds=copyAndTranspose(dsWRF) + t1e = time.perf_counter() + print('{:d}/{:d}: t1_elapsed = {:f} (s)'.format(mpi_rank, mpi_size, t1e-t1s)) + t2s = time.perf_counter() + dsFENew=interp2DForFE(ds,ds_FEGrid,XvWRF,YvWRF,xVec,yVec) + t2e = time.perf_counter() print('{:d}/{:d}: t2_elapsed = {:f} (s)'.format(mpi_rank, mpi_size, t2e-t2s)) t3s = time.perf_counter() dsFEFinal=create_dsFEFinal(ds_FEGrid) diff --git a/scripts/python_utilities/coupler/couplingUtils.py b/scripts/python_utilities/coupler/couplingUtils.py index e12b605..5d6bd2a 100644 --- a/scripts/python_utilities/coupler/couplingUtils.py +++ b/scripts/python_utilities/coupler/couplingUtils.py @@ -127,18 +127,11 @@ def getFEProfileDS(dsWrf,varsList,surfVarsList,zFE): ##### Map (Interp/Extrap-ol surfVarDict = {'TSK':'tskin','Q2':'qskin','HGT':'topoWRF','T2':'t2','PSFC':'psfc'} #Note using Q2 instead of QVG since QVG not default in wrfout files for var in varsList: if var != 'Z': - if True: #JAS -testing check on lack of topoPos in FE for ATAC DFW... - f1=interpolate.interp1d(dsWrf['Z'],dsWrf[var],kind='linear',fill_value='extrapolate') - else: - f1=interpolate.interp1d(dsWrf['Z']-dsWrf['HGT'],dsWrf[var],kind='linear',fill_value='extrapolate') + f1=interpolate.interp1d(dsWrf['Z'],dsWrf[var],kind='linear',fill_value='extrapolate') if var != 'ALT' and var != 'ALB': ds_ret[varDict[var]] = xr.DataArray(f1(zFE),dims=['zIndex']) #,coords={'zIndex':np.array0:zFE.size} else: ds_ret[varDict[var]] = xr.DataArray(1.0/f1(zFE),dims=['zIndex']) #,coords={'zIndex':np.array0:zFE.size} - #Having done all of the 3-d vars, fix rho from modified-moist to dry if moist was included... - ###### THIS TESTED IF ALT WAS LIKE A 1/RHO_M, but according to Klemp 2007 it wouldn't be, rather just 1/rho_d despite using theta_m in rho = P/R(T_m) - #####if 'QVAPOR' in list(dsWrf.keys()): - #####ds_ret['rho'] = ds_ret['rho']/(1.0 + ds_ret['qv']) # Leaving out cloud-water hydrometeor class for now... +ds_ret['ql']) for surfVar in surfVarsList: ds_ret[surfVarDict[surfVar]] = xr.DataArray(dsWrf[surfVar]) return ds_ret @@ -156,18 +149,12 @@ def getWRFProfileDS(it,j,i,dsWrf,varsList,surfVarsList): ##### Destagger and col ds_ret[var] = xr.DataArray((0.5*(dsWrf.PH[it,0:-1,j,i]+dsWrf.PH[it,1:,j,i]) +0.5*(dsWrf.PHB[it,0:-1,j,i]+dsWrf.PHB[it,1:,j,i]))/9.81, dims=(['bottom_top'])) - #ds_ret[var] = xr.DataArray((0.5*(dsWrf.PH[it,0:-1,j,i]+dsWrf.PH[it,1:,j,i]) - # +0.5*(dsWrf.PHB[it,0:-1,j,i]+dsWrf.PHB[it,1:,j,i]))/9.81-dsWrf.HGT[it,j,i], - # dims=(['bottom_top'])) - #elif var == 'U': elif 'west_east_stag' in dsWrf[var].dims: ds_ret[var] = xr.DataArray(0.5*(dsWrf[var][it,:,j,i]+dsWrf[var][it,:,j,i+1]), dims=(['bottom_top'])) - #elif var == 'V': elif 'south_north_stag' in dsWrf[var].dims: ds_ret[var] = xr.DataArray(0.5*(dsWrf[var][it,:,j,i]+dsWrf[var][it,:,j+1,i]), dims=(['bottom_top'])) - #elif var == 'W': elif 'bottom_top_stag' in dsWrf[var].dims: ds_ret[var] = xr.DataArray(0.5*(dsWrf[var][it,0:-1,j,i]+dsWrf[var][it,1:,j,i]), dims=(['bottom_top'])) @@ -177,30 +164,17 @@ def getWRFProfileDS(it,j,i,dsWrf,varsList,surfVarsList): ##### Destagger and col if var == 'T': ds_ret[var] = 300.0+ds_ret[var] - #elif var == 'QVAPOR': - # ds_ret[var] = xr.DataArray(dsWrf[var][it,:,j,i], - # dims=(['bottom_top'])) - #elif var == 'QCLOUD': - # ds_ret[var] = xr.DataArray(dsWrf[var][it,:,j,i], - # dims=(['bottom_top'])) for surfVar in surfVarsList: ds_ret[surfVar] = xr.DataArray(dsWrf[surfVar][it,j,i])#, - #dims=(['surface'])) - #print('-skin[{:d},{:d},{:d}] = {:f}'.format(it,i,j,ds_ret[surfVar].values)) return ds_ret def copyAndTranspose(dsWRF): ds=xr.Dataset() for var in dsWRF.variables: if len(dsWRF[var].dims) == 3: - #ds[var]=xr.DataArray(np.transpose(dsWRF[var].values,(2,1,0)),dims=('zIndex','yIndex','xIndex')) - # !!!!!! Hey !!!!! - #Looking closer at resulting theta fields, it looked a bit transposed so trying this new counter intuitive (2,0,1) ordering of the 3-D transpose... ds[var]=xr.DataArray(np.transpose(dsWRF[var].values,(2,0,1)),dims=('zIndex','yIndex','xIndex')) - #testing do not transpose 2-d fields elif len(dsWRF[var].dims) == 2: ds[var]=xr.DataArray(dsWRF[var].values,dims=('yIndex','xIndex')) #Note no transpose needed here... - # ds[var]=xr.DataArray(np.transpose(dsWRF[var].values,(1,0)),dims=('yIndex','xIndex')) return ds def interp2DForFE(ds,ds_FE,XvWRF,YvWRF,xVec,yVec): @@ -244,12 +218,9 @@ def create_dsFEFinal(ds_FE): def verticalInterpFinal(ds_FE,dsFENew,dsFEFinal,zRect): z3d=ds_FE['zPos'][:,:,:].values.squeeze() for var in dsFENew.variables: - #if (var in FEvarsList or FEsurfVarsList) and var in dsFEFinal.variables: - #if var not in ['ql', 'XLAT','XLONG','topoWRF']: print(var) t1s = time.perf_counter() print(var,list(dsFENew[var].dims)) - #if len(dsFENew[var].dims) == 3: if 'zIndex' in list(dsFENew[var].dims): if 'time' in list(dsFENew[var].dims): fld3dRect=dsFENew[var][0,:,:,:].values.squeeze() @@ -257,12 +228,10 @@ def verticalInterpFinal(ds_FE,dsFENew,dsFEFinal,zRect): fld3dRect=dsFENew[var].values print(z3d.shape,fld3dRect.shape) tmpVar3d=interpolateIrregularVertical(zRect,fld3dRect,z3d) - #print(tmpVar3d.shape) if var in list(dsFEFinal.variables): dsFEFinal[var][:,:,:]=tmpVar3d.astype(np.float32) else: print('Omitting {:s} from dsFEFinal...'.format(var)) - #elif len(dsFENew[var].dims) == 2: else: if 'time' in list(dsFENew[var].dims): dsFEFinal[var][0,:,:]=dsFENew[var][0,:,:].values.astype(np.float32) @@ -287,12 +256,8 @@ def interpolateIrregularVertical(zRect,fld3dRect,z3d): NzR,NyR,NxR = fld3dRect.shape Nz3,Ny3,Nx3 = z3d.shape tmpVar3d=np.zeros(z3d.shape) - #if not(NyR is Ny3) or not(NxR is Nx3): - # print('Uh-oh, NyR,NxR = {:d},{:d} but Ny3,Nx3 = {:d},{:d}'.format(NyR,NxR,Ny3,Nx3)) for j in range(NyR): for i in range(NxR): - #fInterp = interpolate.interp1d(zRect,fld3dRect[:,j,i], kind='linear') - #tmp=fInterp(z3d[:,j,i]) tmp=np.interp(z3d[:,j,i],zRect,fld3dRect[:,j,i]) tmpVar3d[:,j,i]=tmp return tmpVar3d @@ -302,18 +267,13 @@ def create_dsBdy(ds_FE,FEvarsList,FEsurfVarsList): notit=0 for var in FEvarsList: print(var) - ### JAS 12-13-2023 Removing the zeroing out ov qv and ql in the boundary planes - if False: - #if var == "qv" or var == "ql": - ds_Bdy[var+'_YZL']=xr.DataArray(0.0*ds_FE['rho'][:,:,:,0]) - else: - if var in list(ds_FE.variables): - ds_Bdy[var+'_YZL']=xr.DataArray(ds_FE[var][:,:,:,0]) - ds_Bdy[var+'_YZH']=xr.DataArray(ds_FE[var][:,:,:,ds_FE.sizes['xIndex']-1]) - ds_Bdy[var+'_XZL']=xr.DataArray(ds_FE[var][:,:,0,:]) - ds_Bdy[var+'_XZH']=xr.DataArray(ds_FE[var][:,:,ds_FE.sizes['yIndex']-1,:]) - ds_Bdy[var+'_XYL']=xr.DataArray(ds_FE[var][:,0,:,:]) - ds_Bdy[var+'_XYH']=xr.DataArray(ds_FE[var][:,ds_FE.sizes['zIndex']-1,:,:]) + if var in list(ds_FE.variables): + ds_Bdy[var+'_YZL']=xr.DataArray(ds_FE[var][:,:,:,0]) + ds_Bdy[var+'_YZH']=xr.DataArray(ds_FE[var][:,:,:,ds_FE.sizes['xIndex']-1]) + ds_Bdy[var+'_XZL']=xr.DataArray(ds_FE[var][:,:,0,:]) + ds_Bdy[var+'_XZH']=xr.DataArray(ds_FE[var][:,:,ds_FE.sizes['yIndex']-1,:]) + ds_Bdy[var+'_XYL']=xr.DataArray(ds_FE[var][:,0,:,:]) + ds_Bdy[var+'_XYH']=xr.DataArray(ds_FE[var][:,ds_FE.sizes['zIndex']-1,:,:]) for surfVar in FEsurfVarsList: print('{:s}: notit={:d}'.format(surfVar,notit)) if surfVar in ['topoWRF','t2','psfc','tskin','qskin']: @@ -323,7 +283,6 @@ def create_dsBdy(ds_FE,FEvarsList,FEsurfVarsList): return ds_Bdy def createBdysFrom3D(ds_Bdy,ds3D,FEvarsList,FEsurfVarsList): - #ds_Bdy=xr.Dataset() for var in FEvarsList: if var in ds3D.variables: ds_Bdy[var+'_YZL']=ds3D[var][0,:,:,0].expand_dims(dim={'time':ds3D.sizes['time']},axis=0) @@ -347,7 +306,6 @@ def createBdysFrom3D(ds_Bdy,ds3D,FEvarsList,FEsurfVarsList): return ds_Bdy def writeBdyFile(path_out_analysis,fileName,ds_Bdy): - #Bdy_filename = path_out_analysis+"FE_Bndys_WRF_thetaDryRho_HSBal.{:d}".format(fileCounter) Bdy_filename = path_out_analysis+fileName #"FE_Bndys_WRF_PDG_w-rhoALT-moist_qv-gpkg.{:d}".format(fileCounter) print(Bdy_filename) ds_Bdy.isel(time = [0]).to_netcdf(Bdy_filename,format='NETCDF4') #note the isel(time =[#]) preserves the time dimension