Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix UFS replay when using active LSM #434

Merged
merged 18 commits into from
Apr 10, 2024
Merged

Fix UFS replay when using active LSM #434

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
18 commits
Select commit Hold shift + click to select a range
c3e01ea
If using model initial conditions and forcing from the UFS (e.g. Repl…
dustinswales Feb 22, 2024
53548ea
Some changes for Man to test
dustinswales Feb 28, 2024
f22b07e
Replay with RUC and NoahMP LSM working
dustinswales Feb 28, 2024
8e21d36
Small cahnge from physics needed for replay
dustinswales Feb 29, 2024
3f06191
Synced physics
dustinswales Feb 29, 2024
5a7cab0
Merge branch 'main' of https://github.com/NCAR/ccpp-scm into replay_f…
dustinswales Feb 29, 2024
82dd9ca
Cleanup/Housekeeping. Remove unnecessary interpolation step
dustinswales Mar 11, 2024
7f7c80b
Add capability to use input file with longitude and latitude lists
dustinswales Mar 21, 2024
650e207
Omission from previous commit
dustinswales Mar 21, 2024
dac377e
Added utility script to create input file containing replay points
dustinswales Mar 21, 2024
b93bc51
Added step to skip gridpoint in ensemble if error during case generation
dustinswales Mar 21, 2024
9e3663a
Synced physics
dustinswales Apr 3, 2024
0ddfd34
Merge branch 'main' of https://github.com/NCAR/ccpp-scm into replay_f…
dustinswales Apr 3, 2024
528aefe
Merge branch 'main' of https://github.com/NCAR/ccpp-scm into replay_f…
dustinswales Apr 10, 2024
c25071b
Address reviewers comments
dustinswales Apr 10, 2024
88b3df8
Update physics
dustinswales Apr 10, 2024
1c12612
Merge branch 'main' of https://github.com/NCAR/ccpp-scm into replay_f…
dustinswales Apr 10, 2024
51c5420
Revert submodules
dustinswales Apr 10, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion ccpp/physics
Open in desktop
Submodule physics updated 1 files
+7 −2 physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.F90
124 changes: 68 additions & 56 deletions scm/etc/scripts/UFS_IC_generator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -464,18 +464,20 @@ def get_UFS_IC_data(dir, grid_dir, tile, i, j, old_chgres, lam):
#returns dictionaries with the data

vgrid_data = get_UFS_vgrid_data(grid_dir) #only needed for ak, bk to calculate pressure
state_data = get_UFS_state_data(vgrid_data, dir, tile, i, j, old_chgres, lam)
(state_data, error_msg) = get_UFS_state_data(vgrid_data, dir, tile, i, j, old_chgres, lam)
surface_data = get_UFS_surface_data(dir, tile, i, j, old_chgres, lam)
oro_data = get_UFS_oro_data(dir, tile, i, j, lam)

return (state_data, surface_data, oro_data)
return (state_data, surface_data, oro_data, error_msg)

########################################################################################
#
########################################################################################
def get_UFS_state_data(vgrid, dir, tile, i, j, old_chgres, lam):
"""Get the state data for the given tile and indices"""


state = {}
error_msg=None
if lam:
nc_file_data = Dataset('{0}/{1}'.format(dir,'gfs_data.nc'))
else:
Expand Down Expand Up @@ -600,7 +602,10 @@ def get_UFS_state_data(vgrid, dir, tile, i, j, old_chgres, lam):
icewat_model_rev[k] = cloud_water - liqwat_model_rev[0,k]
(liqwat_model_rev[0,k], dummy_rain, icewat_model_rev[k], dummy_snow) = fv3_remap.mp_auto_conversion(liqwat_model_rev[0,k], icewat_model_rev[k])

[u_s, u_n, v_w, v_e] = get_zonal_and_meridional_winds_on_cd_grid(tile, dir, i, j, nc_file_data, lam)
[u_s, u_n, v_w, v_e, unknown_grid] = get_zonal_and_meridional_winds_on_cd_grid(tile, dir, i, j, nc_file_data, lam)
if unknown_grid:
error_msg='unknown grid orientation'
return(state,error_msg)

#put C/D grid zonal/meridional winds on model pressure levels
u_s_model_rev = fv3_remap.mappm(levp_data, pressure_from_data_rev[np.newaxis, :], u_s[np.newaxis, :], nlevs_model, pressure_model_interfaces_rev[np.newaxis, :], 1, 1, -1, 8, ptop_data)
Expand Down Expand Up @@ -638,7 +643,7 @@ def get_UFS_state_data(vgrid, dir, tile, i, j, old_chgres, lam):
"pa_i": pressure_model_interfaces
}

return state
return (state,error_msg)

########################################################################################
#
Expand All @@ -665,6 +670,10 @@ def get_zonal_and_meridional_winds_on_cd_grid(tile, dir, i, j, nc_file_data, lam
raise Exception(message)

nc_file_grid = Dataset('{0}/{1}'.format(dir,filename))

# Get grid dimension
nz,nx,ny = np.shape(nc_file_data['u_w'])


if (lam):
#strip ghost/halo points and return supergrid
Expand Down Expand Up @@ -696,7 +705,7 @@ def get_zonal_and_meridional_winds_on_cd_grid(tile, dir, i, j, nc_file_data, lam

east_test_point = np.argmax(test_lon_diff)
north_test_point = np.argmax(test_lat_diff)

unknown=False
if east_test_point == 0:
#longitude increases most along the positive i axis
if north_test_point == 2:
Expand Down Expand Up @@ -775,7 +784,11 @@ def get_zonal_and_meridional_winds_on_cd_grid(tile, dir, i, j, nc_file_data, lam
(ex, ey) = fv3_remap.get_latlon_vector(p3)
v_e = nc_file_data['u_w'][:,j,i+1]*fv3_remap.inner_prod(e1, ex) + nc_file_data['v_w'][:,j,i+1]*fv3_remap.inner_prod(e1, ey)
else:
print('unknown grid orientation')
u_s = np.zeros(nz)
u_n = np.zeros(nz)
v_w = np.zeros(nz)
v_e = np.zeros(nz)
unknown=True
elif east_test_point == 1:
#longitude increases most along the negative i axis
if north_test_point == 2:
Expand Down Expand Up @@ -853,7 +866,11 @@ def get_zonal_and_meridional_winds_on_cd_grid(tile, dir, i, j, nc_file_data, lam
(ex, ey) = fv3_remap.get_latlon_vector(p3)
v_e = nc_file_data['u_w'][:,j,i-1]*fv3_remap.inner_prod(e1, ex) + nc_file_data['v_w'][:,j,i-1]*fv3_remap.inner_prod(e1, ey)
else:
print('unknown grid orientation')
u_s = np.zeros(nz)
u_n = np.zeros(nz)
v_w = np.zeros(nz)
v_e = np.zeros(nz)
unknown=True
elif east_test_point == 2:
#longitude increases most along the positive j axis
if north_test_point == 0:
Expand Down Expand Up @@ -929,9 +946,16 @@ def get_zonal_and_meridional_winds_on_cd_grid(tile, dir, i, j, nc_file_data, lam
p3 = fv3_remap.mid_pt_sphere(p1*deg_to_rad, p2*deg_to_rad)
e1 = fv3_remap.get_unit_vect2(p1*deg_to_rad, p2*deg_to_rad)
(ex, ey) = fv3_remap.get_latlon_vector(p3)
v_e = nc_file_data['u_w'][:,j+1,i]*fv3_remap.inner_prod(e1, ex) + nc_file_data['v_w'][:,j+1,i]*fv3_remap.inner_prod(e1, ey)
if (j < nx -1):
v_e = nc_file_data['u_w'][:,j+1,i]*fv3_remap.inner_prod(e1, ex) + nc_file_data['v_w'][:,j+1,i]*fv3_remap.inner_prod(e1, ey)
else:
v_e = nc_file_data['v_w'][:,j,i]
else:
print('unknown grid orientation')
u_s = np.zeros(nz)
u_n = np.zeros(nz)
v_w = np.zeros(nz)
v_e = np.zeros(nz)
unknown=True
elif east_test_point == 3:
#longitude increases most along the negative j axis
if north_test_point == 0:
Expand Down Expand Up @@ -1009,12 +1033,15 @@ def get_zonal_and_meridional_winds_on_cd_grid(tile, dir, i, j, nc_file_data, lam
(ex, ey) = fv3_remap.get_latlon_vector(p3)
v_e = nc_file_data['u_w'][:,j-1,i]*fv3_remap.inner_prod(e1, ex) + nc_file_data['v_w'][:,j-1,i]*fv3_remap.inner_prod(e1, ey)
else:
print('unknown grid orientation')

u_s = np.zeros(nz)
u_n = np.zeros(nz)
v_w = np.zeros(nz)
v_e = np.zeros(nz)
unknown=True

nc_file_grid.close()

return [u_s, u_n, v_w, v_e]
return [u_s, u_n, v_w, v_e, unknown]

########################################################################################
#
Expand Down Expand Up @@ -1066,10 +1093,14 @@ def get_UFS_surface_data(dir, tile, i, j, old_chgres, lam):
tprcp_in = read_NetCDF_surface_var(nc_file, 'tprcp', i, j, old_chgres, 0)
srflag_in = read_NetCDF_surface_var(nc_file, 'srflag', i, j, old_chgres, 0)
sncovr_in = read_NetCDF_surface_var(nc_file, 'sncovr', i, j, old_chgres, 0)
tsfcl_in = read_NetCDF_surface_var(nc_file, 'tsfcl', i, j, old_chgres, 0)
zorll_in = read_NetCDF_surface_var(nc_file, 'zorll', i, j, old_chgres, 0)
zorli_in = read_NetCDF_surface_var(nc_file, 'zorli', i, j, old_chgres, 0)

tsfcl_in = read_NetCDF_surface_var(nc_file, 'tsea', i, j, old_chgres, 0)
zorll_in = read_NetCDF_surface_var(nc_file, 'zorl', i, j, old_chgres, 0)
zorli_in = read_NetCDF_surface_var(nc_file, 'zorl', i, j, old_chgres, 0)
if (snwdph_in > 0):
sncovr_in = 1.0
else:
sncovr_in = 0.0

# present when cplwav = T
zorlw_in = read_NetCDF_surface_var(nc_file, 'zorlw', i, j, old_chgres, 0)

Expand Down Expand Up @@ -1171,6 +1202,12 @@ def get_UFS_surface_data(dir, tile, i, j, old_chgres, lam):
# fractional grid
tiice_in = read_NetCDF_surface_var(nc_file, 'tiice', i, j, old_chgres, missing_variable_ice_layers)

# soil color (From the UFS FV3: io/fv3atm_sfc_io.F90)
if (slmsk_in == 1):
scolor_in = 4
grantfirl marked this conversation as resolved.
Show resolved Hide resolved
else:
scolor_in = 1

#
nc_file.close()

Expand All @@ -1187,6 +1224,7 @@ def get_UFS_surface_data(dir, tile, i, j, old_chgres, lam):
"facsf": facsf_in,
"facwf": facwf_in,
"soiltyp": styp_in,
"scolor": scolor_in,
"slopetyp": slope_in,
"vegtyp": vtyp_in,
"vegfrac": vfrac_in,
Expand Down Expand Up @@ -1906,37 +1944,6 @@ def get_UFS_forcing_data(nlevs, state_IC, location, use_nearest, forcing_dir, gr
pres_adv[t+1,:] = pres_adv[t,:]
pres_i_adv[t+1,:] = pres_i_adv[t,:]

if save_comp_data:
#
t_layr = np.zeros([n_files+1,nlevs])
qv_layr = np.zeros([n_files+1,nlevs])
u_layr = np.zeros([n_files+1,nlevs])
v_layr = np.zeros([n_files+1,nlevs])
p_layr = np.zeros([n_files+1,nlevs])

#
for t in range(0,n_files):
from_p[0,:] = stateNATIVE["p_lev"][t,::-1]
to_p[0,:] = stateNATIVE["p_lev"][1,::-1]
log_from_p[0,:] = np.log(from_p[0,:])
log_to_p[0,:] = np.log(to_p[0,:])
p_layr[t,:] = stateNATIVE["p_lay"][1,::-1]
for k in range(0,nlevs): dp2[0,k] = to_p[0,k+1] - to_p[0,k]
t_layr[t,:] = fv3_remap.map_scalar(nlevs, log_from_p, stateNATIVE["t_lay"][t:t+1,::-1], \
dummy, nlevs, log_to_p, 0, 0, 1, np.abs(kord_tm), t_min)
qv_layr[t,:] = fv3_remap.map1_q2(nlevs, from_p, stateNATIVE["qv_lay"][t:t+1,::-1], \
nlevs, to_p, dp2, 0, 0, 0, kord_tr, q_min)
u_layr[t,:] = fv3_remap.map1_ppm(nlevs, from_p, stateNATIVE["u_lay"][t:t+1,::-1], \
0.0, nlevs, to_p, 0, 0, -1, kord_tm)
v_layr[t,:] = fv3_remap.map1_ppm(nlevs, from_p, stateNATIVE["v_lay"][t:t+1,::-1], \
0.0, nlevs, to_p, 0, 0, -1, kord_tm)

t_layr[t+1,:] = t_layr[t,:]
qv_layr[t+1,:] = qv_layr[t,:]
u_layr[t+1,:] = u_layr[t,:]
v_layr[t+1,:] = v_layr[t,:]
p_layr[t+1,:] = p_layr[t,:]

####################################################################################
#
# if we had atmf,sfcf files at every timestep (and the SCM timestep is made to match
Expand Down Expand Up @@ -2105,11 +2112,11 @@ def get_UFS_forcing_data(nlevs, state_IC, location, use_nearest, forcing_dir, gr
if (save_comp_data):
comp_data = {
"time": stateNATIVE["time"]*sec_in_hr,
"pa" : p_layr[:,::-1],
"ta" : t_layr[:,::-1],
"qv" : qv_layr[:,::-1],
"ua" : u_layr[:,::-1],
"va" : v_layr[:,::-1],
"pa" : stateNATIVE["p_lay"][:,:],
"ta" : stateNATIVE["t_lay"][:,:],
"qv" : stateNATIVE["qv_lay"][:,:],
"ua" : stateNATIVE["u_lay"][:,:],
"va" : stateNATIVE["v_lay"][:,:],
"vars2d":vars2d}
else:
comp_data = {}
Expand Down Expand Up @@ -2392,7 +2399,7 @@ def write_SCM_case_file(state, surface, oro, forcing, case, date, stateREGRID):
{"name": "mrsos_forc", "type":wp, "dimd": ('time' ), "units": "kg m-2", "desc": "forcing_mass_content_of_water_in_soil_layer"}]

#
var_oro = [{"name": "area", "type":wp, "dimd": ('t0'), "units": "m 2-1", "desc": "grid_cell_area"},\
var_oro = [{"name": "area", "type":wp, "dimd": ('t0'), "units": "m2", "desc": "grid_cell_area"},\
{"name": "stddev", "type":wp, "dimd": ('t0'), "units": "m", "desc": "standard deviation of subgrid orography"}, \
{"name": "convexity", "type":wp, "dimd": ('t0'), "units": "none", "desc": "convexity of subgrid orography"}, \
{"name": "oa1", "type":wp, "dimd": ('t0'), "units": "none", "desc": "assymetry of subgrid orography 1"}, \
Expand Down Expand Up @@ -2450,6 +2457,7 @@ def write_SCM_case_file(state, surface, oro, forcing, case, date, stateREGRID):
{"name": "q2m", "type":wp, "dimd": ('t0'), "units": "kg kg-1", "desc": "2-meter specific humidity"}, \
{"name": "vegtyp", "type":wi, "dimd": ('t0'), "units": "none", "desc": "vegetation type (1-12)"}, \
{"name": "soiltyp", "type":wi, "dimd": ('t0'), "units": "none", "desc": "soil type (1-12)"}, \
{"name": "scolor", "type":wp, "dimd": ('t0'), "units": "none", "desc": "soil color"}, \
{"name": "ffmm", "type":wp, "dimd": ('t0'), "units": "none", "desc": "Monin-Obukhov similarity function for momentum"}, \
{"name": "ffhh", "type":wp, "dimd": ('t0'), "units": "none", "desc": "Monin-Obukhov similarity function for heat"}, \
{"name": "hice", "type":wp, "dimd": ('t0'), "units": "m", "desc": "sea ice thickness"}, \
Expand Down Expand Up @@ -2704,8 +2712,12 @@ def main():

# get UFS IC data (TODO: flag to read in RESTART data rather than IC data and implement
# different file reads)
(state_data, surface_data, oro_data) = get_UFS_IC_data(in_dir, grid_dir, tile, tile_i,\
tile_j, old_chgres, lam)
(state_data, surface_data, oro_data, error_msg) = get_UFS_IC_data(in_dir, grid_dir, tile, tile_i,\
tile_j, old_chgres, lam)
if (error_msg):
print(error_msg)
print("STOPPING")
exit()

if not date:
# date was not included on command line; look in atmf* file for initial date
Expand Down
Loading
Loading