diff --git a/ccpp/config/ccpp_prebuild_config.py b/ccpp/config/ccpp_prebuild_config.py index 31406939..d2f75268 100755 --- a/ccpp/config/ccpp_prebuild_config.py +++ b/ccpp/config/ccpp_prebuild_config.py @@ -22,6 +22,7 @@ 'scm/src/gmtb_scm_kinds.F90', 'scm/src/gmtb_scm_type_defs.F90', 'scm/src/gmtb_scm_physical_constants.F90', + 'scm/src/gmtb_scm_utils.F90', #no definitions, but gmtb_scm_type_defs.F90 uses a module from this file 'ccpp/physics/physics/rte-rrtmgp/rrtmgp/mo_gas_optics_rrtmgp.F90', 'ccpp/physics/physics/rte-rrtmgp/rrtmgp/mo_gas_concentrations.F90', 'ccpp/physics/physics/rte-rrtmgp/rte/mo_optical_props.F90', @@ -110,6 +111,8 @@ # 'ccpp/physics/physics/bl_acm.F90' , 'ccpp/physics/physics/cires_ugwp.F90' , 'ccpp/physics/physics/cires_ugwp_post.F90' , + 'ccpp/physics/physics/unified_ugwp.F90' , + 'ccpp/physics/physics/unified_ugwp_post.F90' , 'ccpp/physics/physics/cnvc90.f' , 'ccpp/physics/physics/cs_conv.F90' , 'ccpp/physics/physics/cs_conv_aw_adj.F90' , @@ -155,6 +158,7 @@ 'ccpp/physics/physics/ozphys_2015.f' , 'ccpp/physics/physics/precpd.f' , 'ccpp/physics/physics/phys_tend.F90' , + 'ccpp/physics/physics/tracer_sanitizer.F90' , 'ccpp/physics/physics/radlw_main.F90' , 'ccpp/physics/physics/radsw_main.F90' , 'ccpp/physics/physics/rascnv.F90' , diff --git a/ccpp/framework b/ccpp/framework index dca1240e..f1dc8d6f 160000 --- a/ccpp/framework +++ b/ccpp/framework @@ -1 +1 @@ -Subproject commit dca1240e6f19a5bbcfa0b14aa8526f36e99ed135 +Subproject commit f1dc8d6f038e590508c272070f673d1fd7ea566f diff --git a/ccpp/physics b/ccpp/physics index 8ef88ca4..176ea9ad 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 8ef88ca46c11535fc7984d39ec38d1582f9db5ff +Subproject commit 176ea9ad3e564bdae968ea764d592dece50a4754 diff --git a/scm/data/processed_case_input/fv3_model_point_noah.nc b/scm/data/processed_case_input/fv3_model_point_noah.nc index 685e999d..b9564898 100644 Binary files a/scm/data/processed_case_input/fv3_model_point_noah.nc and b/scm/data/processed_case_input/fv3_model_point_noah.nc differ diff --git a/scm/etc/scripts/UFS_IC_generator.py b/scm/etc/scripts/UFS_IC_generator.py index 62e46111..bd13d75c 100755 --- a/scm/etc/scripts/UFS_IC_generator.py +++ b/scm/etc/scripts/UFS_IC_generator.py @@ -23,7 +23,11 @@ zvir = rvgas/rdgas - 1. grav = 9.80665 -missing_value = 9.99e20 +missing_value = -9999.0 #9.99e20 + +missing_variable_snow_layers = 3 +missing_variable_soil_layers = 4 +missing_variable_ice_layers = 2 # Path to the directory containing processed case input files PROCESSED_CASE_DIR = '../../data/processed_case_input' @@ -278,14 +282,57 @@ def sph2cart(az, el, r): return (x, y, z) -def get_UFS_IC_data(dir, tile, i, j, old_chgres): +def read_NetCDF_var(nc_file, var_name, i, j): + try: + var = nc_file[var_name][j,i] + except (KeyError, IndexError): + message = "Variable {0} is not found in {1}. Filling with missing value.".format(var_name,nc_file.filepath()) + logging.debug(message) + var = missing_value + return var + +def read_NetCDF_surface_var(nc_file, var_name, i, j, old_chgres, vert_dim): + if old_chgres: + if vert_dim > 0: + try: + var = nc_file[var_name][:,j,i] + except (KeyError, IndexError): + message = "Variable {0} is not found in {1}. Filling with array of size {2} with missing value.".format(var_name,nc_file.filepath(),vert_dim) + logging.debug(message) + var = missing_value*np.ma.ones(vert_dim) + else: + try: + var = nc_file[var_name][j,i] + except (KeyError, IndexError): + message = "Variable {0} is not found in {1}. Filling with missing value.".format(var_name,nc_file.filepath()) + logging.debug(message) + var = missing_value + else: + #the sfc_data.tileX.nc files created from chgres_cube have an extra time dimension in front compared to those created from global_chgres + if vert_dim > 0: + try: + var = nc_file[var_name][0,:,j,i] + except (KeyError, IndexError): + message = "Variable {0} is not found in {1}. Filling with array of size {2} with missing value.".format(var_name,nc_file.filepath(),vert_dim) + logging.debug(message) + var = missing_value*np.ma.ones(vert_dim) + else: + try: + var = nc_file[var_name][0,j,i] + except (KeyError, IndexError): + message = "Variable {0} is not found in {1}. Filling with missing value.".format(var_name,nc_file.filepath()) + logging.debug(message) + var = missing_value + return var + +def get_UFS_IC_data(dir, grid_dir, tile, i, j, old_chgres): """Get the state, surface, and orographic data for the given tile and indices""" #returns dictionaries with the data state_data = get_UFS_state_data(dir, tile, i, j, old_chgres) surface_data = get_UFS_surface_data(dir, tile, i, j, old_chgres) oro_data = get_UFS_oro_data(dir, tile, i, j) - vgrid_data = get_UFS_vgrid_data(dir) #only needed for ak, bk to calculate pressure + vgrid_data = get_UFS_vgrid_data(grid_dir) #only needed for ak, bk to calculate pressure #calculate derived quantities if old_chgres: @@ -379,80 +426,142 @@ def get_UFS_surface_data(dir, tile, i, j, old_chgres): nc_file = Dataset('{0}/{1}'.format(dir,'sfc_data.tile{0}.nc'.format(tile))) - if old_chgres: - ts_in=nc_file['tsea'][j,i] + #FV3/io/FV3GFS_io.F90/sfc_prop_restart_read was used as reference for variables that can be read in + + #read in scalars (would be 2D variables in a 3D model) + + # surface properties (assuming Noah LSM; may contain variables needed for fractional land fraction) + tsfco_in = read_NetCDF_surface_var(nc_file, 'tsea', i, j, old_chgres, 0) + tg3_in = read_NetCDF_surface_var(nc_file, 'tg3', i, j, old_chgres, 0) + uustar_in = read_NetCDF_surface_var(nc_file, 'uustar', i, j, old_chgres, 0) + alvsf_in = read_NetCDF_surface_var(nc_file, 'alvsf', i, j, old_chgres, 0) + alvwf_in = read_NetCDF_surface_var(nc_file, 'alvwf', i, j, old_chgres, 0) + alnsf_in = read_NetCDF_surface_var(nc_file, 'alnsf', i, j, old_chgres, 0) + alnwf_in = read_NetCDF_surface_var(nc_file, 'alnwf', i, j, old_chgres, 0) + facsf_in = read_NetCDF_surface_var(nc_file, 'facsf', i, j, old_chgres, 0) + facwf_in = read_NetCDF_surface_var(nc_file, 'facwf', i, j, old_chgres, 0) + styp_in = read_NetCDF_surface_var(nc_file, 'stype', i, j, old_chgres, 0) + slope_in = read_NetCDF_surface_var(nc_file, 'slope', i, j, old_chgres, 0) + vtyp_in = read_NetCDF_surface_var(nc_file, 'vtype', i, j, old_chgres, 0) + vfrac_in = read_NetCDF_surface_var(nc_file, 'vfrac', i, j, old_chgres, 0) + shdmin_in = read_NetCDF_surface_var(nc_file, 'shdmin', i, j, old_chgres, 0) + shdmax_in = read_NetCDF_surface_var(nc_file, 'shdmax', i, j, old_chgres, 0) + zorlo_in = read_NetCDF_surface_var(nc_file, 'zorl', i, j, old_chgres, 0) + slmsk_in = read_NetCDF_surface_var(nc_file, 'slmsk', i, j, old_chgres, 0) + canopy_in = read_NetCDF_surface_var(nc_file, 'canopy', i, j, old_chgres, 0) + hice_in = read_NetCDF_surface_var(nc_file, 'hice', i, j, old_chgres, 0) + fice_in = read_NetCDF_surface_var(nc_file, 'fice', i, j, old_chgres, 0) + tisfc_in = read_NetCDF_surface_var(nc_file, 'tisfc', i, j, old_chgres, 0) + snwdph_in = read_NetCDF_surface_var(nc_file, 'snwdph', i, j, old_chgres, 0) + snoalb_in = read_NetCDF_surface_var(nc_file, 'snoalb', i, j, old_chgres, 0) + sheleg_in = read_NetCDF_surface_var(nc_file, 'sheleg', i, j, old_chgres, 0) + f10m_in = read_NetCDF_surface_var(nc_file, 'f10m', i, j, old_chgres, 0) + t2m_in = read_NetCDF_surface_var(nc_file, 't2m', i, j, old_chgres, 0) + q2m_in = read_NetCDF_surface_var(nc_file, 'q2m', i, j, old_chgres, 0) + ffmm_in = read_NetCDF_surface_var(nc_file, 'ffmm', i, j, old_chgres, 0) + ffhh_in = read_NetCDF_surface_var(nc_file, 'ffhh', i, j, old_chgres, 0) + 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) + + #present when cplwav = T + zorlw_in = read_NetCDF_surface_var(nc_file, 'zorlw', i, j, old_chgres, 0) + + #NSST variables that may be in the surface file + tref_in = read_NetCDF_surface_var(nc_file, 'tref', i, j, old_chgres, 0) + z_c_in = read_NetCDF_surface_var(nc_file, 'z_c', i, j, old_chgres, 0) + c_0_in = read_NetCDF_surface_var(nc_file, 'c_0', i, j, old_chgres, 0) + c_d_in = read_NetCDF_surface_var(nc_file, 'c_d', i, j, old_chgres, 0) + w_0_in = read_NetCDF_surface_var(nc_file, 'w_0', i, j, old_chgres, 0) + w_d_in = read_NetCDF_surface_var(nc_file, 'w_d', i, j, old_chgres, 0) + xt_in = read_NetCDF_surface_var(nc_file, 'xt', i, j, old_chgres, 0) + xs_in = read_NetCDF_surface_var(nc_file, 'xs', i, j, old_chgres, 0) + xu_in = read_NetCDF_surface_var(nc_file, 'xu', i, j, old_chgres, 0) + xv_in = read_NetCDF_surface_var(nc_file, 'xv', i, j, old_chgres, 0) + xz_in = read_NetCDF_surface_var(nc_file, 'xz', i, j, old_chgres, 0) + zm_in = read_NetCDF_surface_var(nc_file, 'zm', i, j, old_chgres, 0) + xtts_in = read_NetCDF_surface_var(nc_file, 'xtts', i, j, old_chgres, 0) + xzts_in = read_NetCDF_surface_var(nc_file, 'xzts', i, j, old_chgres, 0) + d_conv_in = read_NetCDF_surface_var(nc_file, 'd_conv', i, j, old_chgres, 0) + ifd_in = read_NetCDF_surface_var(nc_file, 'ifd', i, j, old_chgres, 0) + dt_cool_in = read_NetCDF_surface_var(nc_file, 'dt_cool', i, j, old_chgres, 0) + qrain_in = read_NetCDF_surface_var(nc_file, 'qrain', i, j, old_chgres, 0) - # land state - stc_in=nc_file['stc'][:,j,i] - smc_in=nc_file['smc'][:,j,i] - slc_in=nc_file['slc'][:,j,i] - tg3_in=nc_file['tg3'][j,i] + #NoahMP variables that may be in the surface file + snowxy_in = read_NetCDF_surface_var(nc_file, 'snowxy', i, j, old_chgres, 0) + tvxy_in = read_NetCDF_surface_var(nc_file, 'tvxy', i, j, old_chgres, 0) + tgxy_in = read_NetCDF_surface_var(nc_file, 'tgxy', i, j, old_chgres, 0) + canicexy_in = read_NetCDF_surface_var(nc_file, 'canicexy', i, j, old_chgres, 0) + canliqxy_in = read_NetCDF_surface_var(nc_file, 'canliqxy', i, j, old_chgres, 0) + eahxy_in = read_NetCDF_surface_var(nc_file, 'eahxy', i, j, old_chgres, 0) + tahxy_in = read_NetCDF_surface_var(nc_file, 'tahxy', i, j, old_chgres, 0) + cmxy_in = read_NetCDF_surface_var(nc_file, 'cmxy', i, j, old_chgres, 0) + chxy_in = read_NetCDF_surface_var(nc_file, 'chxy', i, j, old_chgres, 0) + fwetxy_in = read_NetCDF_surface_var(nc_file, 'fwetxy', i, j, old_chgres, 0) + sneqvoxy_in = read_NetCDF_surface_var(nc_file, 'sneqvoxy', i, j, old_chgres, 0) + alboldxy_in = read_NetCDF_surface_var(nc_file, 'alboldxy', i, j, old_chgres, 0) + qsnowxy_in = read_NetCDF_surface_var(nc_file, 'qsnowxy', i, j, old_chgres, 0) + wslakexy_in = read_NetCDF_surface_var(nc_file, 'wslakexy', i, j, old_chgres, 0) + zwtxy_in = read_NetCDF_surface_var(nc_file, 'zwtxy', i, j, old_chgres, 0) + waxy_in = read_NetCDF_surface_var(nc_file, 'waxy', i, j, old_chgres, 0) + wtxy_in = read_NetCDF_surface_var(nc_file, 'wtxy', i, j, old_chgres, 0) + lfmassxy_in = read_NetCDF_surface_var(nc_file, 'lfmassxy', i, j, old_chgres, 0) + rtmassxy_in = read_NetCDF_surface_var(nc_file, 'rtmassxy', i, j, old_chgres, 0) + stmassxy_in = read_NetCDF_surface_var(nc_file, 'stmassxy', i, j, old_chgres, 0) + woodxy_in = read_NetCDF_surface_var(nc_file, 'woodxy', i, j, old_chgres, 0) + stblcpxy_in = read_NetCDF_surface_var(nc_file, 'stblcpxy', i, j, old_chgres, 0) + fastcpxy_in = read_NetCDF_surface_var(nc_file, 'fastcpxy', i, j, old_chgres, 0) + xsaixy_in = read_NetCDF_surface_var(nc_file, 'xsaixy', i, j, old_chgres, 0) + xlaixy_in = read_NetCDF_surface_var(nc_file, 'xlaixy', i, j, old_chgres, 0) + taussxy_in = read_NetCDF_surface_var(nc_file, 'taussxy', i, j, old_chgres, 0) + smcwtdxy_in = read_NetCDF_surface_var(nc_file, 'smcwtdxy', i, j, old_chgres, 0) + deeprechxy_in = read_NetCDF_surface_var(nc_file, 'deeprechxy', i, j, old_chgres, 0) + rechxy_in = read_NetCDF_surface_var(nc_file, 'rechxy', i, j, old_chgres, 0) + + #RUC LSM variables that may be in the surface file + wetness_in = read_NetCDF_surface_var(nc_file, 'wetness', i, j, old_chgres, 0) + clw_surf_in = read_NetCDF_surface_var(nc_file, 'clw_surf', i, j, old_chgres, 0) + qwv_surf_in = read_NetCDF_surface_var(nc_file, 'qwv_surf', i, j, old_chgres, 0) + tsnow_in = read_NetCDF_surface_var(nc_file, 'tsnow', i, j, old_chgres, 0) + snowfall_acc_in = read_NetCDF_surface_var(nc_file, 'snowfall_acc', i, j, old_chgres, 0) + swe_snowfall_acc_in = read_NetCDF_surface_var(nc_file, 'swe_snowfall_acc', i, j, old_chgres, 0) + lai_in = read_NetCDF_surface_var(nc_file, 'lai', i, j, old_chgres, 0) + + #read in profiles (would be 3D variables in a 3D model) + + #land_state + stc_in = read_NetCDF_surface_var(nc_file, 'stc', i, j, old_chgres, missing_variable_soil_layers) + smc_in = read_NetCDF_surface_var(nc_file, 'smc', i, j, old_chgres, missing_variable_soil_layers) + slc_in = read_NetCDF_surface_var(nc_file, 'slc', i, j, old_chgres, missing_variable_soil_layers) + + #NoahMP 3D variables + snicexy_in = read_NetCDF_surface_var(nc_file, 'snicexy', i, j, old_chgres, missing_variable_snow_layers) + snliqxy_in = read_NetCDF_surface_var(nc_file, 'snliqxy', i, j, old_chgres, missing_variable_snow_layers) + tsnoxy_in = read_NetCDF_surface_var(nc_file, 'tsnoxy', i, j, old_chgres, missing_variable_snow_layers) + smoiseq_in = read_NetCDF_surface_var(nc_file, 'smoiseq', i, j, old_chgres, missing_variable_soil_layers) + zsnsoxy_in = read_NetCDF_surface_var(nc_file, 'zsnsoxy', i, j, old_chgres, missing_variable_soil_layers + missing_variable_snow_layers) + + #RUC LSM 3D variables + tslb_in = read_NetCDF_surface_var(nc_file, 'tslb', i, j, old_chgres, missing_variable_soil_layers) + smois_in = read_NetCDF_surface_var(nc_file, 'smois', i, j, old_chgres, missing_variable_soil_layers) + sh2o_in = read_NetCDF_surface_var(nc_file, 'sh2o', i, j, old_chgres, missing_variable_soil_layers) + smfr_in = read_NetCDF_surface_var(nc_file, 'smfr', i, j, old_chgres, missing_variable_soil_layers) + flfr_in = read_NetCDF_surface_var(nc_file, 'flfr', i, j, old_chgres, missing_variable_soil_layers) + + #fractional grid 3D variables + tiice_in = read_NetCDF_surface_var(nc_file, 'tiice', i, j, old_chgres, missing_variable_ice_layers) - # surface properties - uustar_in=nc_file['uustar'][j,i] - alvsf_in=nc_file['alvsf'][j,i] - alvwf_in=nc_file['alvwf'][j,i] - alnsf_in=nc_file['alnsf'][j,i] - alnwf_in=nc_file['alnwf'][j,i] - facsf_in=nc_file['facsf'][j,i] - facwf_in=nc_file['facwf'][j,i] - styp_in=nc_file['stype'][j,i] - slope_in=nc_file['slope'][j,i] - vtyp_in=nc_file['vtype'][j,i] - vfrac_in=nc_file['vfrac'][j,i] - shdmin_in=nc_file['shdmin'][j,i] - shdmax_in=nc_file['shdmax'][j,i] - zorl_in=nc_file['zorl'][j,i] - slmsk_in=nc_file['slmsk'][j,i] - canopy_in=nc_file['canopy'][j,i] - hice_in=nc_file['hice'][j,i] - fice_in=nc_file['fice'][j,i] - tisfc_in=nc_file['tisfc'][j,i] - snwdph_in=nc_file['snwdph'][j,i] - snoalb_in=nc_file['snoalb'][j,i] - sheleg_in=nc_file['sheleg'][j,i] - else: - #the sfc_data.tileX.nc files created from chgres_cube have an extra time dimension in front compared to those created from global_chgres - ts_in=nc_file['tsea'][0,j,i] - - # land state - stc_in=nc_file['stc'][0,:,j,i] - smc_in=nc_file['smc'][0,:,j,i] - slc_in=nc_file['slc'][0,:,j,i] - tg3_in=nc_file['tg3'][0,j,i] - - # surface properties - uustar_in=nc_file['uustar'][0,j,i] - alvsf_in=nc_file['alvsf'][0,j,i] - alvwf_in=nc_file['alvwf'][0,j,i] - alnsf_in=nc_file['alnsf'][0,j,i] - alnwf_in=nc_file['alnwf'][0,j,i] - facsf_in=nc_file['facsf'][0,j,i] - facwf_in=nc_file['facwf'][0,j,i] - styp_in=nc_file['stype'][0,j,i] - slope_in=nc_file['slope'][0,j,i] - vtyp_in=nc_file['vtype'][0,j,i] - vfrac_in=nc_file['vfrac'][0,j,i] - shdmin_in=nc_file['shdmin'][0,j,i] - shdmax_in=nc_file['shdmax'][0,j,i] - zorl_in=nc_file['zorl'][0,j,i] - slmsk_in=nc_file['slmsk'][0,j,i] - canopy_in=nc_file['canopy'][0,j,i] - hice_in=nc_file['hice'][0,j,i] - fice_in=nc_file['fice'][0,j,i] - tisfc_in=nc_file['tisfc'][0,j,i] - snwdph_in=nc_file['snwdph'][0,j,i] - snoalb_in=nc_file['snoalb'][0,j,i] - sheleg_in=nc_file['sheleg'][0,j,i] + #print("zorlw_in = {}".format(zorlw_in)) nc_file.close() #put data in a dictionary surface = { - "T_surf": ts_in, - "stc": stc_in, - "smc": smc_in, - "slc": slc_in, + #Noah LSM + "tsfco": tsfco_in, "tg3": tg3_in, "uustar": uustar_in, "alvsf": alvsf_in, @@ -467,7 +576,7 @@ def get_UFS_surface_data(dir, tile, i, j, old_chgres): "vfrac": vfrac_in, "shdmin": shdmin_in, "shdmax": shdmax_in, - "zorl": zorl_in, + "zorlo": zorlo_in, "slmsk": slmsk_in, "canopy": canopy_in, "hice": hice_in, @@ -475,7 +584,95 @@ def get_UFS_surface_data(dir, tile, i, j, old_chgres): "tisfc": tisfc_in, "snwdph": snwdph_in, "snoalb": snoalb_in, - "sheleg": sheleg_in + "sheleg": sheleg_in, + "f10m": f10m_in, + "t2m": t2m_in, + "q2m": q2m_in, + "ffmm": ffmm_in, + "ffhh": ffhh_in, + "tprcp": tprcp_in, + "srflag": srflag_in, + "sncovr": sncovr_in, + "tsfcl": tsfcl_in, + "zorll": zorll_in, + "zorli": zorli_in, + #cplwav + "zorlw": zorlw_in, + #NSST + "tref": tref_in, + "z_c": z_c_in, + "c_0": c_0_in, + "c_d": c_d_in, + "w_0": w_0_in, + "w_d": w_d_in, + "xt": xt_in, + "xs": xs_in, + "xu": xu_in, + "xv": xv_in, + "xz": xz_in, + "zm": zm_in, + "xtts": xtts_in, + "xzts": xzts_in, + "d_conv": d_conv_in, + "ifd": ifd_in, + "dt_cool": dt_cool_in, + "qrain": qrain_in, + #NoahMP + "snowxy": snowxy_in, + "tvxy": tvxy_in, + "tgxy": tgxy_in, + "canicexy": canicexy_in, + "canliqxy": canliqxy_in, + "eahxy": eahxy_in, + "tahxy": tahxy_in, + "cmxy": cmxy_in, + "chxy": chxy_in, + "fwetxy": fwetxy_in, + "sneqvoxy": sneqvoxy_in, + "alboldxy": alboldxy_in, + "qsnowxy": qsnowxy_in, + "wslakexy": wslakexy_in, + "zwtxy": zwtxy_in, + "waxy": waxy_in, + "wtxy": wtxy_in, + "lfmassxy": lfmassxy_in, + "rtmassxy": rtmassxy_in, + "stmassxy": stmassxy_in, + "woodxy": woodxy_in, + "stblcpxy": stblcpxy_in, + "fastcpxy": fastcpxy_in, + "xsaixy": xsaixy_in, + "xlaixy": xlaixy_in, + "taussxy": taussxy_in, + "smcwtdxy": smcwtdxy_in, + "deeprechxy": deeprechxy_in, + "rechxy": rechxy_in, + #RUC LSM + "wetness": wetness_in, + "clw_surf": clw_surf_in, + "qwv_surf": qwv_surf_in, + "tsnow": tsnow_in, + "snowfall_acc": snowfall_acc_in, + "swe_snowfall_acc": swe_snowfall_acc_in, + "lai": lai_in, + #Noah LSM 3D + "stc": stc_in, + "smc": smc_in, + "slc": slc_in, + #NoahMP LSM 3D + "snicexy": snicexy_in, + "snliqxy": snliqxy_in, + "tsnoxy": tsnoxy_in, + "smoiseq": smoiseq_in, + "zsnsoxy": zsnsoxy_in, + #RUC LSM 3D variables + "tslb": tslb_in, + "smois": smois_in, + "sh2o": sh2o_in, + "smfr": smfr_in, + "flfr": flfr_in, + #fractional grid 3D variables + "tiice": tiice_in, } return surface @@ -490,20 +687,27 @@ def get_UFS_oro_data(dir, tile, i, j): nc_file = Dataset('{0}/{1}'.format(dir,filename)) # orographyic properties - stddev_in=nc_file['stddev'][j,i] - convexity_in=nc_file['convexity'][j,i] - oa1_in=nc_file['oa1'][j,i] - oa2_in=nc_file['oa2'][j,i] - oa3_in=nc_file['oa3'][j,i] - oa4_in=nc_file['oa4'][j,i] - ol1_in=nc_file['ol1'][j,i] - ol2_in=nc_file['ol2'][j,i] - ol3_in=nc_file['ol3'][j,i] - ol4_in=nc_file['ol4'][j,i] - theta_in=nc_file['theta'][j,i] - gamma_in=nc_file['gamma'][j,i] - sigma_in=nc_file['sigma'][j,i] - elvmax_in=nc_file['elvmax'][j,i] + stddev_in = read_NetCDF_var(nc_file, "stddev", i, j) + convexity_in = read_NetCDF_var(nc_file, "convexity", i, j) + oa1_in = read_NetCDF_var(nc_file, "oa1", i, j) + oa2_in = read_NetCDF_var(nc_file, "oa2", i, j) + oa3_in = read_NetCDF_var(nc_file, "oa3", i, j) + oa4_in = read_NetCDF_var(nc_file, "oa4", i, j) + ol1_in = read_NetCDF_var(nc_file, "ol1", i, j) + ol2_in = read_NetCDF_var(nc_file, "ol2", i, j) + ol3_in = read_NetCDF_var(nc_file, "ol3", i, j) + ol4_in = read_NetCDF_var(nc_file, "ol4", i, j) + theta_in = read_NetCDF_var(nc_file, "theta", i, j) + gamma_in = read_NetCDF_var(nc_file, "gamma", i, j) + sigma_in = read_NetCDF_var(nc_file, "sigma", i, j) + elvmax_in = read_NetCDF_var(nc_file, "elvmax", i, j) + orog_filt_in = read_NetCDF_var(nc_file, "orog_filt", i, j) + orog_raw_in = read_NetCDF_var(nc_file, "orog_raw", i, j) + #fractional landmask variables + land_frac_in = read_NetCDF_var(nc_file, "land_frac", i, j) + #lake variables + lake_frac_in = read_NetCDF_var(nc_file, "lake_frac", i, j) + lake_depth_in = read_NetCDF_var(nc_file, "lake_depth", i, j) nc_file.close() @@ -522,7 +726,12 @@ def get_UFS_oro_data(dir, tile, i, j): "theta": theta_in, "gamma": gamma_in, "sigma": sigma_in, - "elvmax": elvmax_in + "elvmax": elvmax_in, + "orog_filt": orog_filt_in, + "orog_raw": orog_raw_in, + "land_frac": land_frac_in, + "lake_frac": lake_frac_in, + "lake_depth": lake_depth_in } return oro @@ -673,11 +882,11 @@ def add_noahmp_coldstart(surface, date): surface["zsnsoxy"] = np.ones(n_snow_layers + n_soil_layers)*missing_value if surface["slmsk"] > 0.01: - surface["tvxy"] = surface["T_surf"] - surface["tgxy"] = surface["T_surf"] - surface["tahxy"] = surface["T_surf"] + surface["tvxy"] = surface["tsfco"] + surface["tgxy"] = surface["tsfco"] + surface["tahxy"] = surface["tsfco"] - if (surface["snwdph"] > 0.01 and surface["T_surf"] > 273.15 ): + if (surface["snwdph"] > 0.01 and surface["tsfco"] > 273.15 ): surface["tvxy"] = 273.15 surface["tgxy"] = 273.15 surface["tahxy"]= 273.15 @@ -931,15 +1140,12 @@ def write_SCM_case_file(state, surface, oro, forcing, case, date): nlevs = state["nlevs"] nsoil = len(surface["stc"]) - - #find out if noahmp ICs have been generated - noahmp = False - if "snicexy" in surface.keys(): - noahmp = True - nsnow = len(surface["snicexy"]) + nsnow = len(surface["snicexy"]) + nice = len(surface["tiice"]) nc_file = Dataset(os.path.join(PROCESSED_CASE_DIR, case + '.nc'), 'w', format='NETCDF4') nc_file.description = "FV3GFS model profile input (no forcing)" + nc_file.missing_value = missing_value #create groups for scalars, intitialization, and forcing @@ -967,9 +1173,9 @@ def write_SCM_case_file(state, surface, oro, forcing, case, date): soil_depth_var.units = 'm' soil_depth_var.description = 'depth of bottom of soil layers' - if noahmp: - snow_dim = nc_file.createDimension('nsnow',None) - soil_plus_snow_dim = nc_file.createDimension('nsoil_plus_nsnow',None) + snow_dim = nc_file.createDimension('nsnow',None) + soil_plus_snow_dim = nc_file.createDimension('nsoil_plus_nsnow',None) + ice_dim = nc_file.createDimension('nice',None) #initial group @@ -1014,45 +1220,74 @@ def write_SCM_case_file(state, surface, oro, forcing, case, date): ozone_var.description = 'initial profile of ozone mass mixing ratio' stc_var = initial_grp.createVariable('stc',real_type,('nsoil',)) - stc_var[:] = surface['stc'][0:nsoil] + stc_var[:] = surface["stc"][0:nsoil] stc_var.units = "K" stc_var.description = "initial profile of soil temperature" smc_var = initial_grp.createVariable('smc',real_type,('nsoil',)) - smc_var[:] = surface['smc'][0:nsoil] + smc_var[:] = surface["smc"][0:nsoil] smc_var.units = "kg" smc_var.description = "initial profile of soil moisture" slc_var = initial_grp.createVariable('slc',real_type,('nsoil',)) - slc_var[:] = surface['slc'][0:nsoil] + slc_var[:] = surface["slc"][0:nsoil] slc_var.units = "kg" slc_var.description = "initial profile of soil liquid moisture" - if noahmp: - snicexy_var = initial_grp.createVariable('snicexy',real_type,('nsnow',)) - snicexy_var[:] = surface['snicexy'][0:nsnow] - snicexy_var.units = "mm" - snicexy_var.description = "initial profile of snow layer ice" + snicexy_var = initial_grp.createVariable('snicexy',real_type,('nsnow',)) + snicexy_var[:] = surface["snicexy"][0:nsnow] + snicexy_var.units = "mm" + snicexy_var.description = "initial profile of snow layer ice" - snliqxy_var = initial_grp.createVariable('snliqxy',real_type,('nsnow',)) - snliqxy_var[:] = surface['snliqxy'][0:nsnow] - snliqxy_var.units = "mm" - snliqxy_var.description = "initial profile of snow layer ice" + snliqxy_var = initial_grp.createVariable('snliqxy',real_type,('nsnow',)) + snliqxy_var[:] = surface["snliqxy"][0:nsnow] + snliqxy_var.units = "mm" + snliqxy_var.description = "initial profile of snow layer liquid" - tsnoxy_var = initial_grp.createVariable('tsnoxy',real_type,('nsnow',)) - tsnoxy_var[:] = surface['tsnoxy'][0:nsnow] - tsnoxy_var.units = "K" - tsnoxy_var.description = "initial profile of snow layer temperature" + tsnoxy_var = initial_grp.createVariable('tsnoxy',real_type,('nsnow',)) + tsnoxy_var[:] = surface["tsnoxy"][0:nsnow] + tsnoxy_var.units = "K" + tsnoxy_var.description = "initial profile of snow layer temperature" - smoiseq_var = initial_grp.createVariable('smoiseq',real_type,('nsoil',)) - smoiseq_var[:] = surface['smoiseq'][0:nsoil] - smoiseq_var.units = "m3 m-3" - smoiseq_var.description = "initial profile of equilibrium soil water content" + smoiseq_var = initial_grp.createVariable('smoiseq',real_type,('nsoil',)) + smoiseq_var[:] = surface["smoiseq"][0:nsoil] + smoiseq_var.units = "m3 m-3" + smoiseq_var.description = "initial profile of equilibrium soil water content" - zsnsoxy_var = initial_grp.createVariable('zsnsoxy',real_type,('nsoil_plus_nsnow',)) - zsnsoxy_var[:] = surface['zsnsoxy'][0:nsoil + nsnow] - zsnsoxy_var.units = "m" - zsnsoxy_var.description = "layer bottom depth from snow surface" + zsnsoxy_var = initial_grp.createVariable('zsnsoxy',real_type,('nsoil_plus_nsnow',)) + zsnsoxy_var[:] = surface["zsnsoxy"][0:nsoil + nsnow] + zsnsoxy_var.units = "m" + zsnsoxy_var.description = "layer bottom depth from snow surface" + + tiice_var = initial_grp.createVariable('tiice',real_type,('nice',)) + tiice_var[:] = surface["tiice"][0:nice] + tiice_var.units = "K" + tiice_var.description = "sea ice internal temperature" + + tslb_var = initial_grp.createVariable('tslb',real_type,('nsoil',)) + tslb_var[:] = surface["tslb"][0:nsoil] + tslb_var.units = "K" + tslb_var.description = "soil temperature for RUC LSM" + + smois_var = initial_grp.createVariable('smois',real_type,('nsoil',)) + smois_var[:] = surface["smois"][0:nsoil] + smois_var.units = "None" + smois_var.description = "volume fraction of soil moisture for RUC LSM" + + sh2o_var = initial_grp.createVariable('sh2o',real_type,('nsoil',)) + sh2o_var[:] = surface["sh2o"][0:nsoil] + sh2o_var.units = "None" + sh2o_var.description = "volume fraction of unfrozen soil moisture for RUC LSM" + + smfr_var = initial_grp.createVariable('smfr',real_type,('nsoil',)) + smfr_var[:] = surface["smfr"][0:nsoil] + smfr_var.units = "None" + smfr_var.description = "volume fraction of frozen soil moisture for RUC LSM" + + flfr_var = initial_grp.createVariable('flfr',real_type,('nsoil',)) + flfr_var[:] = surface["flfr"][0:nsoil] + flfr_var.units = "None" + flfr_var.description = "flag for frozen soil physics for RUC LSM" #forcing group @@ -1062,9 +1297,9 @@ def write_SCM_case_file(state, surface, oro, forcing, case, date): p_surf_var.description = 'surface pressure' T_surf_var = forcing_grp.createVariable('T_surf', real_type, ('time',)) - T_surf_var[:] = surface["T_surf"] + T_surf_var[:] = missing_value T_surf_var.units = 'K' - T_surf_var.description = 'surface absolute temperature' + T_surf_var.description = 'surface absolute temperature forcing' w_ls_var = forcing_grp.createVariable('w_ls', real_type, ('levels','time',)) w_ls_var[:] = forcing["w_ls"] @@ -1184,6 +1419,11 @@ def write_SCM_case_file(state, surface, oro, forcing, case, date): #Noah initial parameters + tsfco = scalar_grp.createVariable('tsfco',real_type) + tsfco[:] = surface["tsfco"] + tsfco.units = "K" + tsfco.description = "sea surface temperature OR surface skin temperature over land OR sea ice surface skin temperature (depends on value of slmsk)" + vegsrc = scalar_grp.createVariable('vegsrc',int_type) vegsrc[:] = 1 #when would this be 2? vegsrc.description = "vegetation soure (1-2)" @@ -1212,10 +1452,10 @@ def write_SCM_case_file(state, surface, oro, forcing, case, date): shdmax[:] = surface["shdmax"] shdmax.description = "maximum vegetation fraction" - zorl = scalar_grp.createVariable('zorl',real_type) - zorl[:] = surface["zorl"] - zorl.units = "cm" - zorl.description = "surface roughness length" + zorlo = scalar_grp.createVariable('zorlo',real_type) + zorlo[:] = surface["zorlo"] + zorlo.units = "cm" + zorlo.description = "surface roughness length over ocean" islmsk = scalar_grp.createVariable('slmsk',real_type) islmsk[:] = surface["slmsk"] @@ -1248,10 +1488,6 @@ def write_SCM_case_file(state, surface, oro, forcing, case, date): snoalb = scalar_grp.createVariable('snoalb',real_type) snoalb[:] = surface["snoalb"] snoalb.description = "maximum snow albedo" - - sncovr = scalar_grp.createVariable('sncovr',real_type) - sncovr[:] = 0.0 - sncovr.description = "snow area fraction" tg3 = scalar_grp.createVariable('tg3',real_type) tg3[:] = surface["tg3"] @@ -1293,6 +1529,71 @@ def write_SCM_case_file(state, surface, oro, forcing, case, date): facwf.units = "None" facwf.description = "fractional coverage with weak cosz dependency" + weasd = scalar_grp.createVariable('weasd',real_type) + weasd[:] = surface["sheleg"] + weasd.units = "mm" + weasd.description = "water equivalent accumulated snow depth" + + f10m = scalar_grp.createVariable('f10m',real_type) + f10m[:] = surface["f10m"] + f10m.units = "None" + f10m.description = "ratio of sigma level 1 wind and 10m wind" + + t2m = scalar_grp.createVariable('t2m',real_type) + t2m[:] = surface["t2m"] + t2m.units = "K" + t2m.description = "2-meter absolute temperature" + + q2m = scalar_grp.createVariable('q2m',real_type) + q2m[:] = surface["q2m"] + q2m.units = "kg kg-1" + q2m.description = "2-meter specific humidity" + + ffmm = scalar_grp.createVariable('ffmm',real_type) + ffmm[:] = surface["ffmm"] + ffmm.units = "None" + ffmm.description = "Monin-Obukhov similarity function for momentum" + + ffhh = scalar_grp.createVariable('ffhh',real_type) + ffhh[:] = surface["ffhh"] + ffhh.units = "None" + ffhh.description = "Monin-Obukhov similarity function for heat" + + tprcp = scalar_grp.createVariable('tprcp',real_type) + tprcp[:] = surface["tprcp"] + tprcp.units = "m" + tprcp.description = "instantaneous total precipitation amount" + + srflag = scalar_grp.createVariable('srflag',real_type) + srflag[:] = surface["srflag"] + srflag.units = "None" + srflag.description = "snow/rain flag for precipitation" + + sncovr = scalar_grp.createVariable('sncovr',real_type) + sncovr[:] = surface["sncovr"] + sncovr.units = "None" + sncovr.description = "surface snow area fraction" + + tsfcl = scalar_grp.createVariable('tsfcl',real_type) + tsfcl[:] = surface["tsfcl"] + tsfcl.units = "K" + tsfcl.description = "surface skin temperature over land" + + zorll = scalar_grp.createVariable('zorll',real_type) + zorll[:] = surface["zorll"] + zorll.units = "cm" + zorll.description = "surface roughness length over land" + + zorli = scalar_grp.createVariable('zorli',real_type) + zorli[:] = surface["zorli"] + zorli.units = "cm" + zorli.description = "surface roughness length over ice" + + zorlw = scalar_grp.createVariable('zorlw',real_type) + zorlw[:] = surface["zorlw"] + zorlw.units = "cm" + zorlw.description = "surface roughness length from wave model" + #Orography initial parameters stddev = scalar_grp.createVariable('stddev',real_type) @@ -1365,152 +1666,303 @@ def write_SCM_case_file(state, surface, oro, forcing, case, date): elvmax.units = "m" elvmax.description = "maximum of subgrid orography" + orog_filt = scalar_grp.createVariable('oro',real_type) + orog_filt[:] = oro["orog_filt"] + orog_filt.units = "m" + orog_filt.description = "orography" + + orog_raw = scalar_grp.createVariable('oro_uf',real_type) + orog_raw[:] = oro["orog_raw"] + orog_raw.units = "m" + orog_raw.description = "unfiltered orography" + + land_frac = scalar_grp.createVariable('landfrac',real_type) + land_frac[:] = oro["land_frac"] + land_frac.units = "None" + land_frac.description = "fraction of horizontal grid area occupied by land" + + lake_frac = scalar_grp.createVariable('lakefrac',real_type) + lake_frac[:] = oro["lake_frac"] + lake_frac.units = "None" + lake_frac.description = "fraction of horizontal grid area occupied by lake" + + lake_depth = scalar_grp.createVariable('lakedepth',real_type) + lake_depth[:] = oro["lake_depth"] + lake_depth.units = "m" + lake_depth.description = "lake depth" + #NoahMP initial scalar parameters - if noahmp: - tvxy = scalar_grp.createVariable('tvxy',real_type) - tvxy[:] = surface["tvxy"] - tvxy.units = "K" - tvxy.description = "vegetation temperature" - - tgxy = scalar_grp.createVariable('tgxy',real_type) - tgxy[:] = surface["tgxy"] - tgxy.units = "K" - tgxy.description = "ground temperature for NoahMP" - - tahxy = scalar_grp.createVariable('tahxy',real_type) - tahxy[:] = surface["tahxy"] - tahxy.units = "K" - tahxy.description = "canopy air temperature" - - canicexy = scalar_grp.createVariable('canicexy',real_type) - canicexy[:] = surface["canicexy"] - canicexy.units = "mm" - canicexy.description = "canopy intercepted ice mass" - - canliqxy = scalar_grp.createVariable('canliqxy',real_type) - canliqxy[:] = surface["canliqxy"] - canliqxy.units = "mm" - canliqxy.description = "canopy intercepted liquid water" - - eahxy = scalar_grp.createVariable('eahxy',real_type) - eahxy[:] = surface["eahxy"] - eahxy.units = "Pa" - eahxy.description = "canopy air vapor pressure" - - cmxy = scalar_grp.createVariable('cmxy',real_type) - cmxy[:] = surface["cmxy"] - cmxy.units = "" - cmxy.description = "surface drag coefficient for momentum for NoahMP" - - chxy = scalar_grp.createVariable('chxy',real_type) - chxy[:] = surface["chxy"] - chxy.units = "" - chxy.description = "surface exchange coeff heat & moisture for NoahMP" - - fwetxy = scalar_grp.createVariable('fwetxy',real_type) - fwetxy[:] = surface["fwetxy"] - fwetxy.units = "" - fwetxy.description = "area fraction of canopy that is wetted/snowed" - - sneqvoxy = scalar_grp.createVariable('sneqvoxy',real_type) - sneqvoxy[:] = surface["sneqvoxy"] - sneqvoxy.units = "mm" - sneqvoxy.description = "snow mass at previous time step" - - alboldxy = scalar_grp.createVariable('alboldxy',real_type) - alboldxy[:] = surface["alboldxy"] - alboldxy.units = "" - alboldxy.description = "snow albedo at previous time step" - - qsnowxy = scalar_grp.createVariable('qsnowxy',real_type) - qsnowxy[:] = surface["qsnowxy"] - qsnowxy.units = "mm s-1" - qsnowxy.description = "snow precipitation rate at surface" - - wslakexy = scalar_grp.createVariable('wslakexy',real_type) - wslakexy[:] = surface["wslakexy"] - wslakexy.units = "mm" - wslakexy.description = "lake water storage" - - taussxy = scalar_grp.createVariable('taussxy',real_type) - taussxy[:] = surface["taussxy"] - taussxy.units = "" - taussxy.description = "non-dimensional snow age" - - waxy = scalar_grp.createVariable('waxy',real_type) - waxy[:] = surface["waxy"] - waxy.units = "mm" - waxy.description = "water storage in aquifer" - - wtxy = scalar_grp.createVariable('wtxy',real_type) - wtxy[:] = surface["wtxy"] - wtxy.units = "mm" - wtxy.description = "water storage in aquifer and saturated soil" - - zwtxy = scalar_grp.createVariable('zwtxy',real_type) - zwtxy[:] = surface["zwtxy"] - zwtxy.units = "m" - zwtxy.description = "water table depth" - - xlaixy = scalar_grp.createVariable('xlaixy',real_type) - xlaixy[:] = surface["xlaixy"] - xlaixy.units = "" - xlaixy.description = "leaf area index" - - xsaixy = scalar_grp.createVariable('xsaixy',real_type) - xsaixy[:] = surface["xsaixy"] - xsaixy.units = "" - xsaixy.description = "stem area index" + tvxy = scalar_grp.createVariable('tvxy',real_type) + tvxy[:] = surface["tvxy"] + tvxy.units = "K" + tvxy.description = "vegetation temperature" + + tgxy = scalar_grp.createVariable('tgxy',real_type) + tgxy[:] = surface["tgxy"] + tgxy.units = "K" + tgxy.description = "ground temperature for NoahMP" + + tahxy = scalar_grp.createVariable('tahxy',real_type) + tahxy[:] = surface["tahxy"] + tahxy.units = "K" + tahxy.description = "canopy air temperature" + + canicexy = scalar_grp.createVariable('canicexy',real_type) + canicexy[:] = surface["canicexy"] + canicexy.units = "mm" + canicexy.description = "canopy intercepted ice mass" + + canliqxy = scalar_grp.createVariable('canliqxy',real_type) + canliqxy[:] = surface["canliqxy"] + canliqxy.units = "mm" + canliqxy.description = "canopy intercepted liquid water" + + eahxy = scalar_grp.createVariable('eahxy',real_type) + eahxy[:] = surface["eahxy"] + eahxy.units = "Pa" + eahxy.description = "canopy air vapor pressure" + + cmxy = scalar_grp.createVariable('cmxy',real_type) + cmxy[:] = surface["cmxy"] + cmxy.units = "" + cmxy.description = "surface drag coefficient for momentum for NoahMP" + + chxy = scalar_grp.createVariable('chxy',real_type) + chxy[:] = surface["chxy"] + chxy.units = "" + chxy.description = "surface exchange coeff heat & moisture for NoahMP" + + fwetxy = scalar_grp.createVariable('fwetxy',real_type) + fwetxy[:] = surface["fwetxy"] + fwetxy.units = "" + fwetxy.description = "area fraction of canopy that is wetted/snowed" + + sneqvoxy = scalar_grp.createVariable('sneqvoxy',real_type) + sneqvoxy[:] = surface["sneqvoxy"] + sneqvoxy.units = "mm" + sneqvoxy.description = "snow mass at previous time step" + + alboldxy = scalar_grp.createVariable('alboldxy',real_type) + alboldxy[:] = surface["alboldxy"] + alboldxy.units = "" + alboldxy.description = "snow albedo at previous time step" + + qsnowxy = scalar_grp.createVariable('qsnowxy',real_type) + qsnowxy[:] = surface["qsnowxy"] + qsnowxy.units = "mm s-1" + qsnowxy.description = "snow precipitation rate at surface" + + wslakexy = scalar_grp.createVariable('wslakexy',real_type) + wslakexy[:] = surface["wslakexy"] + wslakexy.units = "mm" + wslakexy.description = "lake water storage" + + taussxy = scalar_grp.createVariable('taussxy',real_type) + taussxy[:] = surface["taussxy"] + taussxy.units = "" + taussxy.description = "non-dimensional snow age" + + waxy = scalar_grp.createVariable('waxy',real_type) + waxy[:] = surface["waxy"] + waxy.units = "mm" + waxy.description = "water storage in aquifer" + + wtxy = scalar_grp.createVariable('wtxy',real_type) + wtxy[:] = surface["wtxy"] + wtxy.units = "mm" + wtxy.description = "water storage in aquifer and saturated soil" + + zwtxy = scalar_grp.createVariable('zwtxy',real_type) + zwtxy[:] = surface["zwtxy"] + zwtxy.units = "m" + zwtxy.description = "water table depth" + + xlaixy = scalar_grp.createVariable('xlaixy',real_type) + xlaixy[:] = surface["xlaixy"] + xlaixy.units = "" + xlaixy.description = "leaf area index" + + xsaixy = scalar_grp.createVariable('xsaixy',real_type) + xsaixy[:] = surface["xsaixy"] + xsaixy.units = "" + xsaixy.description = "stem area index" - lfmassxy = scalar_grp.createVariable('lfmassxy',real_type) - lfmassxy[:] = surface["lfmassxy"] - lfmassxy.units = "g m-2" - lfmassxy.description = "leaf mass" - - stmassxy = scalar_grp.createVariable('stmassxy',real_type) - stmassxy[:] = surface["stmassxy"] - stmassxy.units = "g m-2" - stmassxy.description = "stem mass" - - rtmassxy = scalar_grp.createVariable('rtmassxy',real_type) - rtmassxy[:] = surface["rtmassxy"] - rtmassxy.units = "g m-2" - rtmassxy.description = "fine root mass" - - woodxy = scalar_grp.createVariable('woodxy',real_type) - woodxy[:] = surface["woodxy"] - woodxy.units = "g m-2" - woodxy.description = "wood mass including woody roots" - - stblcpxy = scalar_grp.createVariable('stblcpxy',real_type) - stblcpxy[:] = surface["stblcpxy"] - stblcpxy.units = "g m-2" - stblcpxy.description = "stable carbon in deep soil" - - fastcpxy = scalar_grp.createVariable('fastcpxy',real_type) - fastcpxy[:] = surface["fastcpxy"] - fastcpxy.units = "g m-2" - fastcpxy.description = "short-lived carbon in shallow soil" - - smcwtdxy = scalar_grp.createVariable('smcwtdxy',real_type) - smcwtdxy[:] = surface["smcwtdxy"] - smcwtdxy.units = "m3 m-3" - smcwtdxy.description = "soil water content between the bottom of the soil and the water table" - - deeprechxy = scalar_grp.createVariable('deeprechxy',real_type) - deeprechxy[:] = surface["deeprechxy"] - deeprechxy.units = "m" - deeprechxy.description = "recharge to or from the water table when deep" - - rechxy = scalar_grp.createVariable('rechxy',real_type) - rechxy[:] = surface["rechxy"] - rechxy.units = "m" - rechxy.description = "recharge to or from the water table when shallow" - - snowxy = scalar_grp.createVariable('snowxy',real_type) - snowxy[:] = surface["snowxy"] - snowxy.units = "" - snowxy.description = "number of snow layers" + lfmassxy = scalar_grp.createVariable('lfmassxy',real_type) + lfmassxy[:] = surface["lfmassxy"] + lfmassxy.units = "g m-2" + lfmassxy.description = "leaf mass" + + stmassxy = scalar_grp.createVariable('stmassxy',real_type) + stmassxy[:] = surface["stmassxy"] + stmassxy.units = "g m-2" + stmassxy.description = "stem mass" + + rtmassxy = scalar_grp.createVariable('rtmassxy',real_type) + rtmassxy[:] = surface["rtmassxy"] + rtmassxy.units = "g m-2" + rtmassxy.description = "fine root mass" + + woodxy = scalar_grp.createVariable('woodxy',real_type) + woodxy[:] = surface["woodxy"] + woodxy.units = "g m-2" + woodxy.description = "wood mass including woody roots" + + stblcpxy = scalar_grp.createVariable('stblcpxy',real_type) + stblcpxy[:] = surface["stblcpxy"] + stblcpxy.units = "g m-2" + stblcpxy.description = "stable carbon in deep soil" + + fastcpxy = scalar_grp.createVariable('fastcpxy',real_type) + fastcpxy[:] = surface["fastcpxy"] + fastcpxy.units = "g m-2" + fastcpxy.description = "short-lived carbon in shallow soil" + + smcwtdxy = scalar_grp.createVariable('smcwtdxy',real_type) + smcwtdxy[:] = surface["smcwtdxy"] + smcwtdxy.units = "m3 m-3" + smcwtdxy.description = "soil water content between the bottom of the soil and the water table" + + deeprechxy = scalar_grp.createVariable('deeprechxy',real_type) + deeprechxy[:] = surface["deeprechxy"] + deeprechxy.units = "m" + deeprechxy.description = "recharge to or from the water table when deep" + + rechxy = scalar_grp.createVariable('rechxy',real_type) + rechxy[:] = surface["rechxy"] + rechxy.units = "m" + rechxy.description = "recharge to or from the water table when shallow" + + snowxy = scalar_grp.createVariable('snowxy',real_type) + snowxy[:] = surface["snowxy"] + snowxy.units = "" + snowxy.description = "number of snow layers" + + #NSST initial scalar parameters + tref = scalar_grp.createVariable('tref',real_type) + tref[:] = surface["tref"] + tref.units = "K" + tref.description = "sea surface reference temperature for NSST" + + z_c = scalar_grp.createVariable('z_c',real_type) + z_c[:] = surface["z_c"] + z_c.units = "m" + z_c.description = "sub-layer cooling thickness for NSST" + + c_0 = scalar_grp.createVariable('c_0',real_type) + c_0[:] = surface["c_0"] + c_0.units = "" + c_0.description = "coefficient 1 to calculate d(Tz)/d(Ts) for NSST" + + c_d = scalar_grp.createVariable('c_d',real_type) + c_d[:] = surface["c_d"] + c_d.units = "" + c_d.description = "coefficient 2 to calculate d(Tz)/d(Ts) for NSST" + + w_0 = scalar_grp.createVariable('w_0',real_type) + w_0[:] = surface["w_0"] + w_0.units = "" + w_0.description = "coefficient 3 to calculate d(Tz)/d(Ts) for NSST" + + w_d = scalar_grp.createVariable('w_d',real_type) + w_d[:] = surface["w_d"] + w_d.units = "" + w_d.description = "coefficient 4 to calculate d(Tz)/d(Ts) for NSST" + + xt = scalar_grp.createVariable('xt',real_type) + xt[:] = surface["xt"] + xt.units = "K m" + xt.description = "heat content in diurnal thermocline layer for NSST" + + xs = scalar_grp.createVariable('xs',real_type) + xs[:] = surface["xs"] + xs.units = "ppt m" + xs.description = "salinity content in diurnal thermocline layer for NSST" + + xu = scalar_grp.createVariable('xu',real_type) + xu[:] = surface["xu"] + xu.units = "m2 s-1" + xu.description = "u-current in diurnal thermocline layer for NSST" + + xv = scalar_grp.createVariable('xv',real_type) + xv[:] = surface["xv"] + xv.units = "m2 s-1" + xv.description = "v-current in diurnal thermocline layer for NSST" + + xz = scalar_grp.createVariable('xz',real_type) + xz[:] = surface["xz"] + xz.units = "m" + xz.description = "thickness of diurnal thermocline layer for NSST" + + zm = scalar_grp.createVariable('zm',real_type) + zm[:] = surface["zm"] + zm.units = "m" + zm.description = "thickness of ocean mixed layer for NSST" + + xtts = scalar_grp.createVariable('xtts',real_type) + xtts[:] = surface["xtts"] + xtts.units = "m" + xtts.description = "sensitivity of diurnal thermocline layer heat content to surface temperature [d(xt)/d(ts)] for NSST" + + xzts = scalar_grp.createVariable('xzts',real_type) + xzts[:] = surface["xzts"] + xzts.units = "m K-1" + xzts.description = "sensitivity of diurnal thermocline layer thickness to surface temperature [d(xz)/d(ts)] for NSST" + + d_conv = scalar_grp.createVariable('d_conv',real_type) + d_conv[:] = surface["d_conv"] + d_conv.units = "m" + d_conv.description = "thickness of free convection layer for NSST" + + ifd = scalar_grp.createVariable('ifd',real_type) + ifd[:] = surface["ifd"] + ifd.units = "" + ifd.description = "index to start DTM run for NSST" + + dt_cool = scalar_grp.createVariable('dt_cool',real_type) + dt_cool[:] = surface["dt_cool"] + dt_cool.units = "K" + dt_cool.description = "sub-layer cooling amount for NSST" + + qrain = scalar_grp.createVariable('qrain',real_type) + qrain[:] = surface["qrain"] + qrain.units = "W" + qrain.description = "sensible heat due to rainfall for NSST" + + #RUC LSM + wetness = scalar_grp.createVariable('wetness',real_type) + wetness[:] = surface["wetness"] + wetness.units = "" + wetness.description = "normalized soil wetness for RUC LSM" + + clw_surf = scalar_grp.createVariable('clw_surf',real_type) + clw_surf[:] = surface["clw_surf"] + clw_surf.units = "kg kg-1" + clw_surf.description = "cloud condensed water mixing ratio at surface for RUC LSM" + + qwv_surf = scalar_grp.createVariable('qwv_surf',real_type) + qwv_surf[:] = surface["qwv_surf"] + qwv_surf.units = "kg kg-1" + qwv_surf.description = "water vapor mixing ratio at surface for RUC LSM" + + tsnow = scalar_grp.createVariable('tsnow',real_type) + tsnow[:] = surface["tsnow"] + tsnow.units = "K" + tsnow.description = "snow temperature at the bottom of the first snow layer for RUC LSM" + + snowfall_acc = scalar_grp.createVariable('snowfall_acc',real_type) + snowfall_acc[:] = surface["snowfall_acc"] + snowfall_acc.units = "kg m-2" + snowfall_acc.description = "run-total snow accumulation on the ground for RUC LSM" + + swe_snowfall_acc = scalar_grp.createVariable('swe_snowfall_acc',real_type) + swe_snowfall_acc[:] = surface["swe_snowfall_acc"] + swe_snowfall_acc.units = "kg m-2" + swe_snowfall_acc.description = "snow water equivalent of run-total frozen precip for RUC LSM" + + lai = scalar_grp.createVariable('lai',real_type) + lai[:] = surface["lai"] + lai.units = "" + lai.description = "leaf area index for RUC LSM" nc_file.close() @@ -1531,7 +1983,7 @@ def main(): #find index of closest point in the tile if indices are not specified if not indices: - (tile_j, tile_i, point_lon, point_lat, dist_min) = find_loc_indices(location, in_dir, tile) + (tile_j, tile_i, point_lon, point_lat, dist_min) = find_loc_indices(location, grid_dir, tile) print 'The closest point in tile {0} has indices [{1},{2}]'.format(tile,tile_i,tile_j) print 'This index has a central longitude/latitude of [{0},{1}]'.format(point_lon,point_lat) print 'This grid cell is approximately {0} km away from the desired location of {1} {2}'.format(dist_min/1.0E3,location[0],location[1]) @@ -1539,12 +1991,12 @@ def main(): tile_i = indices[0] tile_j = indices[1] #still need to grab the lon/lat if the tile and indices are supplied - (point_lon, point_lat) = find_lon_lat_of_indices(indices, in_dir, tile) + (point_lon, point_lat) = find_lon_lat_of_indices(indices, grid_dir, tile) print 'This index has a central longitude/latitude of [{0},{1}]'.format(point_lon,point_lat) #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, tile, tile_i, tile_j, old_chgres) + (state_data, surface_data, oro_data) = get_UFS_IC_data(in_dir, grid_dir, tile, tile_i, tile_j, old_chgres) #cold start NoahMP variables if (noahmp): diff --git a/scm/etc/scripts/plot_configs/arm_sgp_summer_1997_A_all_suites.ini b/scm/etc/scripts/plot_configs/arm_sgp_summer_1997_A_all_suites.ini new file mode 100644 index 00000000..02ca8bca --- /dev/null +++ b/scm/etc/scripts/plot_configs/arm_sgp_summer_1997_A_all_suites.ini @@ -0,0 +1,57 @@ +gmtb_scm_datasets = output_arm_sgp_summer_1997_A_SCM_GFS_v15p2/output.nc, output_arm_sgp_summer_1997_A_SCM_GFS_v16beta/output.nc, output_arm_sgp_summer_1997_A_SCM_csawmg/output.nc, output_arm_sgp_summer_1997_A_SCM_GSD_v1/output.nc, +gmtb_scm_datasets_labels = GFSv15p2, GFSv16beta, csawmg, GSDv1 +plot_dir = plots_arm_sgp_summer_1997_A/ +obs_file = ../data/raw_case_input/sgp3hIOPsndgBasedV2.0_ConstrVarAnaX1.c1.19970618.000000.cdf +obs_compare = True +plot_ind_datasets = False +time_series_resample = True + +[time_slices] + [[A]] + start = 1997, 6, 26, 23 + end = 1997, 6, 30, 23 + +[time_snapshots] + +[plots] + [[profiles_mean]] + vars = qc, qv, T, dT_dt_PBL, dT_dt_conv, dT_dt_micro, dT_dt_lwrad, dT_dt_swrad + vars_labels = 'cloud water mixing ratio ($g$ $kg^{-1}$)', 'specific humidity ($g$ $kg^{-1}$)', 'T (K)', 'PBL tendency (K/day)', 'conv. tendency (K/day)', 'microphysics tendency (K/day)', 'LW tendency (K/day)', 'SW tendency (K/day)' + vert_axis = pres_l + vert_axis_label = 'average pressure (Pa)' + y_inverted = True + y_log = False + y_min_option = min #min, max, val (if val, add y_min = float value) + y_max_option = max #min, max, val (if val, add y_max = float value) + conversion_factor = 1000.0, 1000.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 + + [[profiles_mean_multi]] + [[[T_forcing]]] + vars = T_force_tend, dT_dt_PBL, dT_dt_conv, dT_dt_micro, dT_dt_lwrad, dT_dt_swrad + vars_labels = 'force', 'PBL', 'Conv', 'MP', 'LW', 'SW' + x_label = 'K/day' + [[[conv_tendencies]]] + vars = dT_dt_deepconv, dT_dt_shalconv + vars_labels = 'deep', 'shallow' + x_label = 'K/day' + + [[profiles_instant]] + + [[time_series]] + vars = 'pres_s','lhf','shf','rain' + vars_labels = 'surface pressure (Pa)','latent heat flux ($W$ $m^{-2}$)','sensible heat flux ($W$ $m^{-2}$)','surface rainfall rate ($mm$ $hr{-1}$)' + + [[contours]] + vars = qv, + vars_labels = 'Water Vapor ($g$ $kg^{-1}$)', + vert_axis = pres_l + vert_axis_label = 'p (Pa)' + y_inverted = True + y_log = False + y_min_option = val #min, max, val (if val, add y_min = float value) + y_min = 10000.0 + y_max_option = val #min, max, val (if val, add y_max = float value) + y_max = 100000.0 + x_ticks_num = 10 + y_ticks_num = 10 + conversion_factor = 1000.0, diff --git a/scm/etc/scripts/plot_configs/arm_sgp_summer_1997_A_physics_test.ini b/scm/etc/scripts/plot_configs/arm_sgp_summer_1997_A_physics_test.ini deleted file mode 100644 index 9eab0037..00000000 --- a/scm/etc/scripts/plot_configs/arm_sgp_summer_1997_A_physics_test.ini +++ /dev/null @@ -1,56 +0,0 @@ -gmtb_scm_datasets = output_arm_sgp_summer_1997_A_SCM_GFS_v15/output.nc, output_arm_sgp_summer_1997_A_SCM_GFS_v15plus/output.nc, output_arm_sgp_summer_1997_A_SCM_csawmg/output.nc, output_arm_sgp_summer_1997_A_SCM_GSD_v0/output.nc -gmtb_scm_datasets_labels = GFSv15, GFSv15+, csawmg, GSDv0 -plot_dir = plots_arm_sgp_summer_1997_A/ -obs_file = ../data/raw_case_input/sgp3hIOPsndgBasedV2.0_ConstrVarAnaX1.c1.19970618.000000.cdf -obs_compare = True -plot_ind_datasets = False -time_series_resample = True - -[time_slices] - [[A]] - start = 1997, 6, 26, 23 - end = 1997, 6, 30, 23 - -[time_snapshots] - -[plots] - [[profiles_mean]] - vars = T, qv, cld, - vars_labels = 'T (k)', 'q (kg/kg)', 'cloud fraction', - vert_axis = pres_l - vert_axis_label = 'p (Pa)' - y_inverted = True - y_log = False - y_min_option = min #min, max, val (if val, add y_min = float value) - y_max_option = max #min, max, val (if val, add y_max = float value) - - [[profiles_mean_multi]] - [[[T_forcing]]] - vars = T_force_tend, dT_dt_PBL, dT_dt_deepconv, dT_dt_shalconv, dT_dt_micro, dT_dt_lwrad, dT_dt_swrad - vars_labels = 'force', 'PBL', 'Deep Con', 'Shal Con', 'MP', 'LW', 'SW' - x_label = 'K/day' - [[[q_forcing]]] - vars = qv_force_tend, dq_dt_PBL, dq_dt_deepconv, dq_dt_shalconv, dq_dt_micro - vars_labels = 'force', 'PBL', 'Deep Con', 'Shal Con', 'MP', 'LW' - x_label = 'g/kg/day' - - - [[profiles_instant]] - - [[time_series]] - vars = rain, - vars_labels = 'rain rate (m/s)', - - [[contours]] - vars = T, qv - vars_labels = 'Temperature','Water Vapor' - vert_axis = pres_l - vert_axis_label = 'p (Pa)' - y_inverted = True - y_log = False - y_min_option = val #min, max, val (if val, add y_min = float value) - y_min = 10000.0 - y_max_option = val #min, max, val (if val, add y_max = float value) - y_max = 100000.0 - x_ticks_num = 10 - y_ticks_num = 10 diff --git a/scm/src/GFS_typedefs.F90 b/scm/src/GFS_typedefs.F90 index 52bd6d43..ef4be32e 100644 --- a/scm/src/GFS_typedefs.F90 +++ b/scm/src/GFS_typedefs.F90 @@ -343,7 +343,6 @@ module GFS_typedefs real (kind=kind_phys), pointer :: tslb(:,:) => null() !< soil temperature for land surface model real (kind=kind_phys), pointer :: flag_frsoil(:,:) => null() !< RUC LSM: flag for frozen soil physics ! - real (kind=kind_phys), pointer :: zs(:) => null() !< depth of soil levels for land surface model real (kind=kind_phys), pointer :: clw_surf(:) => null() !< RUC LSM: moist cloud water mixing ratio at surface real (kind=kind_phys), pointer :: qwv_surf(:) => null() !< RUC LSM: water vapor mixing ratio at surface real (kind=kind_phys), pointer :: cndm_surf(:) => null() !< RUC LSM: surface condensation mass @@ -763,6 +762,7 @@ module GFS_typedefs integer :: lsoil_lsm !< number of soil layers internal to land surface model integer :: lsnow_lsm !< maximum number of snow layers internal to land surface model integer :: lsnow_lsm_lbound!< lower bound for snow arrays, depending on lsnow_lsm + real(kind=kind_phys), pointer :: zs(:) => null() !< depth of soil levels for land surface model logical :: rdlai !< read LAI from input file (for RUC LSM or NOAH LSM WRFv4) logical :: ua_phys !< flag for using University of Arizona? extension to NOAH LSM WRFv4 logical :: usemonalb !< flag to read surface diffused shortwave albedo from input file for NOAH LSM WRFv4 @@ -808,9 +808,17 @@ module GFS_typedefs logical :: old_monin !< flag for diff monin schemes logical :: cnvgwd !< flag for conv gravity wave drag integer :: gwd_opt !< gwd_opt = 1 => original GFS gwd (gwdps.f) - !< gwd_opt = 2 => unified GWD (placeholder) - !< gwd_opt = 3 => GSD drag suite - !< gwd_opt = 33 => GSD drag suite with extra output + !< gwd_opt = 2 => unified ugwp GWD + !< gwd_opt = 22 => unified ugwp GWD with extra output + !< gwd_opt = 3 => GSL drag suite + !< gwd_opt = 33 => GSL drag suite with extra output + logical :: do_ugwp_v0 !< flag for version 0 ugwp GWD + logical :: do_ugwp_v0_orog_only !< flag for version 0 ugwp GWD (orographic drag only) + logical :: do_gsl_drag_ls_bl !< flag for GSL drag (large-scale GWD and blocking only) + logical :: do_gsl_drag_ss !< flag for GSL drag (small-scale GWD only) + logical :: do_gsl_drag_tofd !< flag for GSL drag (turbulent orog form drag only) + logical :: do_ugwp_v1 !< flag for version 1 ugwp GWD + logical :: do_ugwp_v1_orog_only !< flag for version 1 ugwp GWD (orographic drag only) logical :: mstrat !< flag for moorthi approach for stratus logical :: moist_adj !< flag for moist convective adjustment logical :: cscnv !< flag for Chikira-Sugiyama convection @@ -877,7 +885,7 @@ module GFS_typedefs integer :: isatmedmf_vdifq = 1 !< flag for updated version of satmedmf (as of May 2019) integer :: nmtvr !< number of topographic variables such as variance etc !< used in the GWD parameterization - 10 more added if - !< GSD orographic drag scheme is used + !< GSL orographic drag scheme is used integer :: jcap !< number of spectral wave trancation used only by sascnv shalcnv real(kind=kind_phys) :: cs_parm(10) !< tunable parameters for Chikira-Sugiyama convection real(kind=kind_phys) :: flgmin(2) !< [in] ice fraction bounds @@ -2043,18 +2051,18 @@ module GFS_typedefs type(ty_source_func_lw) :: sources !< RRTMGP DDT !-- HWRF physics: dry mixing ratios - real (kind=kind_phys), pointer :: qv_r(:,:) => null() !< - real (kind=kind_phys), pointer :: qc_r(:,:) => null() !< - real (kind=kind_phys), pointer :: qi_r(:,:) => null() !< - real (kind=kind_phys), pointer :: qr_r(:,:) => null() !< - real (kind=kind_phys), pointer :: qs_r(:,:) => null() !< - real (kind=kind_phys), pointer :: qg_r(:,:) => null() !< - - !-- GSD drag suite - real (kind=kind_phys), pointer :: varss(:) => null() !< - real (kind=kind_phys), pointer :: ocss(:) => null() !< - real (kind=kind_phys), pointer :: oa4ss(:,:) => null() !< - real (kind=kind_phys), pointer :: clxss(:,:) => null() !< + real (kind=kind_phys), pointer :: qv_r(:,:) => null() !< + real (kind=kind_phys), pointer :: qc_r(:,:) => null() !< + real (kind=kind_phys), pointer :: qi_r(:,:) => null() !< + real (kind=kind_phys), pointer :: qr_r(:,:) => null() !< + real (kind=kind_phys), pointer :: qs_r(:,:) => null() !< + real (kind=kind_phys), pointer :: qg_r(:,:) => null() !< + + !-- GSL drag suite + real (kind=kind_phys), pointer :: varss(:) => null() !< + real (kind=kind_phys), pointer :: ocss(:) => null() !< + real (kind=kind_phys), pointer :: oa4ss(:,:) => null() !< + real (kind=kind_phys), pointer :: clxss(:,:) => null() !< !-- Ferrier-Aligo MP scheme real (kind=kind_phys), pointer :: f_rain (:,:) => null() !< @@ -2475,7 +2483,6 @@ subroutine sfcprop_create (Sfcprop, IM, Model) allocate (Sfcprop%smois (IM,Model%lsoil_lsm)) allocate (Sfcprop%tslb (IM,Model%lsoil_lsm)) allocate (Sfcprop%flag_frsoil (IM,Model%lsoil_lsm)) - allocate (Sfcprop%zs (Model%lsoil_lsm)) allocate (Sfcprop%clw_surf (IM)) allocate (Sfcprop%qwv_surf (IM)) allocate (Sfcprop%cndm_surf (IM)) @@ -2489,7 +2496,6 @@ subroutine sfcprop_create (Sfcprop, IM, Model) Sfcprop%keepsmfr = clear_val Sfcprop%smois = clear_val Sfcprop%tslb = clear_val - Sfcprop%zs = clear_val Sfcprop%clw_surf = clear_val Sfcprop%qwv_surf = clear_val Sfcprop%cndm_surf = clear_val @@ -3038,7 +3044,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & integer :: iopt_tbot = 2 !lower boundary of soil temperature (1->zero-flux; 2->noah) integer :: iopt_stc = 1 !snow/soil temperature time scheme (only layer 1) - logical :: use_ufo = .false. !< flag for gcycle surface option + logical :: use_ufo = .false. !< flag for gcycle surface option logical :: lcurr_sf = .false. !< flag for taking ocean currents into account in GFDL surface layer logical :: pert_cd = .false. !< flag for perturbing the surface drag coefficient for momentum in surface layer scheme @@ -3056,8 +3062,17 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & logical :: old_monin = .false. !< flag for diff monin schemes logical :: cnvgwd = .false. !< flag for conv gravity wave drag integer :: gwd_opt = 1 !< flag for configuring gwd scheme - !< gwd_opt = 3 : GSDdrag suite - !< gwd_opt = 33: GSDdrag suite with extra output + !< gwd_opt = 2 => unified ugwp GWD + !< gwd_opt = 22 => unified ugwp GWD with extra output + !< gwd_opt = 3 : GSL drag suite + !< gwd_opt = 33: GSL drag suite with extra output + logical :: do_ugwp_v0 = .true. !< flag for version 0 ugwp GWD + logical :: do_ugwp_v0_orog_only = .false. !< flag for version 0 ugwp GWD (orographic drag only) + logical :: do_gsl_drag_ls_bl = .false. !< flag for GSL drag (large-scale GWD and blocking only) + logical :: do_gsl_drag_ss = .false. !< flag for GSL drag (small-scale GWD only) + logical :: do_gsl_drag_tofd = .false. !< flag for GSL drag (turbulent orog form drag only) + logical :: do_ugwp_v1 = .false. !< flag for version 1 ugwp GWD + logical :: do_ugwp_v1_orog_only = .false. !< flag for version 1 ugwp GWD (orographic drag only) !--- vay-2018 logical :: ldiag_ugwp = .false. !< flag for UGWP diag fields logical :: do_ugwp = .false. !< flag do UGWP+RF @@ -3342,7 +3357,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & bl_mynn_cloudpdf, bl_mynn_edmf, bl_mynn_edmf_mom, & bl_mynn_edmf_tke, bl_mynn_edmf_part, bl_mynn_cloudmix, & bl_mynn_mixqt, bl_mynn_output, icloud_bl, bl_mynn_tkeadvect, & - gwd_opt, & + gwd_opt, do_ugwp_v0, do_ugwp_v0_orog_only, & + do_gsl_drag_ls_bl, do_gsl_drag_ss, do_gsl_drag_tofd, & + do_ugwp_v1, do_ugwp_v1_orog_only, & var_ric, coef_ric_l, coef_ric_s, hurr_pbl, & do_myjsfc, do_myjpbl, & hwrf_samfdeep, hwrf_samfshal, & @@ -3743,9 +3760,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !--- land/surface model parameters Model%lsm = lsm Model%lsoil = lsoil - ! Consistency check for RUC LSM - if ((Model%lsm == Model%lsm_ruc .or. Model%lsm == Model%lsm_noah_wrfv4) .and. Model%nscyc>0) then - write(0,*) 'Logic error: RUC LSM and NOAH WRFv4 LSM cannot be used with surface data cycling at this point (fhcyc>0)' + ! Consistency check for HWRF Noah LSM + if (Model%lsm == Model%lsm_noah_wrfv4 .and. Model%nscyc>0) then + write(0,*) 'Logic error: NOAH WRFv4 LSM cannot be used with surface data cycling at this point (fhcyc>0)' stop end if ! Flag to read leaf area index from input files (initial conditions) @@ -3760,6 +3777,12 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & else Model%lsoil_lsm = lsoil_lsm end if + ! Allocate variable to store depth of soil layers + if (Model%lsm==Model%lsm_ruc) then + allocate (Model%zs(Model%lsoil_lsm)) + Model%zs = clear_val + end if + ! if (lsnow_lsm /= 3) then write(0,*) 'Logic error: NoahMP expects the maximum number of snow layers to be exactly 3 (see sfc_noahmp_drv.f)' stop @@ -3905,12 +3928,20 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%coef_ric_s = coef_ric_s Model%gwd_opt = gwd_opt - if (Model%gwd_opt==3 .or. Model%gwd_opt==33) then - ! Add 10 more orographic static fields for GSD drag scheme + if (Model%gwd_opt==3 .or. Model%gwd_opt==33 .or. & + Model%gwd_opt==2 .or. Model%gwd_opt==22) then + ! Add 10 more orographic static fields for GSL drag scheme Model%nmtvr = 24 end if - Model%do_myjsfc = do_myjsfc - Model%do_myjpbl = do_myjpbl + Model%do_ugwp_v0 = do_ugwp_v0 + Model%do_ugwp_v0_orog_only = do_ugwp_v0_orog_only + Model%do_gsl_drag_ls_bl = do_gsl_drag_ls_bl + Model%do_gsl_drag_ss = do_gsl_drag_ss + Model%do_gsl_drag_tofd = do_gsl_drag_tofd + Model%do_ugwp_v1 = do_ugwp_v1 + Model%do_ugwp_v1_orog_only = do_ugwp_v1_orog_only + Model%do_myjsfc = do_myjsfc + Model%do_myjpbl = do_myjpbl !--- Rayleigh friction Model%prslrd0 = prslrd0 @@ -4472,8 +4503,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%nieffr = 2 Model%nseffr = 3 if (.not. Model%effr_in) then - print *,' Thompson MP requires effr_in to be set to .true. - job aborted' - stop + print *,' Thompson MP requires effr_in to be set to .true. - job aborted' + stop end if if (Model%me == Model%master) print *,' Using Thompson double moment microphysics', & ' ltaerosol = ',Model%ltaerosol, & @@ -4810,8 +4841,6 @@ subroutine control_print(Model) print *, ' rdlai : ', Model%rdlai print *, ' lsoil_lsm : ', Model%lsoil_lsm print *, ' lsnow_lsm : ', Model%lsnow_lsm - print *, ' ivegsrc : ', Model%ivegsrc - print *, ' isot : ', Model%isot print *, ' iopt_thcnd : ', Model%iopt_thcnd print *, ' ua_phys : ', Model%ua_phys print *, ' usemonalb : ', Model%usemonalb @@ -4903,7 +4932,15 @@ subroutine control_print(Model) print *, ' do_mynnsfclay : ', Model%do_mynnsfclay print *, ' do_myjsfc : ', Model%do_myjsfc print *, ' do_myjpbl : ', Model%do_myjpbl + print *, ' do_ugwp : ', Model%do_ugwp print *, ' gwd_opt : ', Model%gwd_opt + print *, ' do_ugwp_v0 : ', Model%do_ugwp_v0 + print *, ' do_ugwp_v0_orog_only : ', Model%do_ugwp_v0_orog_only + print *, ' do_gsl_drag_ls_bl : ', Model%do_gsl_drag_ls_bl + print *, ' do_gsl_drag_ss : ', Model%do_gsl_drag_ss + print *, ' do_gsl_drag_tofd : ', Model%do_gsl_drag_tofd + print *, ' do_ugwp_v1 : ', Model%do_ugwp_v1 + print *, ' do_ugwp_v1_orog_only : ', Model%do_ugwp_v1_orog_only print *, ' hurr_pbl : ', Model%hurr_pbl print *, ' var_ric : ', Model%var_ric print *, ' coef_ric_l : ', Model%coef_ric_l @@ -5619,7 +5656,7 @@ subroutine diag_create (Diag, IM, Model) endif !--- Drag Suite variables: - if (Model%gwd_opt == 33) then + if (Model%gwd_opt == 33 .or. Model%gwd_opt == 22) then !print*,"Allocating all Drag Suite variables:" allocate (Diag%dtaux2d_ls (IM,Model%levs)) allocate (Diag%dtauy2d_ls (IM,Model%levs)) @@ -5789,12 +5826,33 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center) Diag%TRAIN = zero end if Diag%cldfra = zero - + Diag%totprcpb = zero Diag%cnvprcpb = zero Diag%toticeb = zero Diag%totsnwb = zero Diag%totgrpb = zero +! +!--- MYNN variables: + if (Model%do_mynnedmf) then + if (Model%bl_mynn_output .ne. 0) then + Diag%edmf_a = clear_val + Diag%edmf_w = clear_val + Diag%edmf_qt = clear_val + Diag%edmf_thl = clear_val + Diag%edmf_ent = clear_val + Diag%edmf_qc = clear_val + Diag%sub_thl = clear_val + Diag%sub_sqv = clear_val + Diag%det_thl = clear_val + Diag%det_sqv = clear_val + endif + Diag%nupdraft = 0 + Diag%maxmf = clear_val + Diag%ktop_plume = 0 + Diag%exch_h = clear_val + Diag%exch_m = clear_val + endif ! if (Model%do_ca) then Diag%ca1 = zero @@ -6317,8 +6375,9 @@ subroutine interstitial_create (Interstitial, IM, Model) allocate (Interstitial%dudt_mtb (IM,Model%levs)) allocate (Interstitial%dudt_ogw (IM,Model%levs)) allocate (Interstitial%dudt_tms (IM,Model%levs)) -!-- GSD drag suite - if (Model%gwd_opt==3 .or. Model%gwd_opt==33) then +!-- GSL drag suite + if (Model%gwd_opt==3 .or. Model%gwd_opt==33 .or. & + Model%gwd_opt==2 .or. Model%gwd_opt==22 ) then allocate (Interstitial%varss (IM)) allocate (Interstitial%ocss (IM)) allocate (Interstitial%oa4ss (IM,4)) @@ -6628,7 +6687,7 @@ subroutine interstitial_rad_reset (Interstitial, Model) ! class(GFS_interstitial_type) :: Interstitial type(GFS_control_type), intent(in) :: Model - + ! Interstitial%aerodp = clear_val Interstitial%alb1d = clear_val if (.not. Model%do_RRTMGP) then @@ -6954,8 +7013,9 @@ subroutine interstitial_phys_reset (Interstitial, Model) Interstitial%dudt_mtb = clear_val Interstitial%dudt_ogw = clear_val Interstitial%dudt_tms = clear_val -!-- GSD drag suite - if (Model%gwd_opt==3 .or. Model%gwd_opt==33) then +!-- GSL drag suite + if (Model%gwd_opt==3 .or. Model%gwd_opt==33 .or. & + Model%gwd_opt==2 .or. Model%gwd_opt==22) then Interstitial%varss = clear_val Interstitial%ocss = clear_val Interstitial%oa4ss = clear_val @@ -7341,8 +7401,9 @@ subroutine interstitial_print(Interstitial, Model, mpirank, omprank, blkno) write (0,*) 'sum(Interstitial%dudt_mtb ) = ', sum(Interstitial%dudt_mtb ) write (0,*) 'sum(Interstitial%dudt_ogw ) = ', sum(Interstitial%dudt_ogw ) write (0,*) 'sum(Interstitial%dudt_tms ) = ', sum(Interstitial%dudt_tms ) -!-- GSD drag suite - if (Model%gwd_opt==3 .or. Model%gwd_opt==33) then +!-- GSL drag suite + if (Model%gwd_opt==3 .or. Model%gwd_opt==33 .or. & + Model%gwd_opt==2 .or. Model%gwd_opt==22) then write (0,*) 'sum(Interstitial%varss ) = ', sum(Interstitial%varss) write (0,*) 'sum(Interstitial%ocss ) = ', sum(Interstitial%ocss) write (0,*) 'sum(Interstitial%oa4ss ) = ', sum(Interstitial%oa4ss) diff --git a/scm/src/GFS_typedefs.meta b/scm/src/GFS_typedefs.meta index 338f52a1..18c61a34 100644 --- a/scm/src/GFS_typedefs.meta +++ b/scm/src/GFS_typedefs.meta @@ -1282,14 +1282,6 @@ type = real kind = kind_phys active = (flag_for_land_surface_scheme == flag_for_ruc_land_surface_scheme) -[zs] - standard_name = depth_of_soil_levels_for_land_surface_model - long_name = depth of soil levels for land surface model - units = m - dimensions = (soil_vertical_dimension_for_land_surface_model) - type = real - kind = kind_phys - active = (flag_for_land_surface_scheme == flag_for_ruc_land_surface_scheme) [clw_surf] standard_name = cloud_condensed_water_mixing_ratio_at_surface long_name = moist cloud water mixing ratio at surface @@ -3269,9 +3261,17 @@ units = count dimensions = () type = integer +[zs] + standard_name = depth_of_soil_levels_for_land_surface_model + long_name = depth of soil levels for land surface model + units = m + dimensions = (soil_vertical_dimension_for_land_surface_model) + type = real + kind = kind_phys + active = (flag_for_land_surface_scheme == flag_for_ruc_land_surface_scheme) [rdlai] standard_name = flag_for_reading_leaf_area_index_from_input - long_name = flag for reading leaf area index from initial conditions for RUC LSM + long_name = flag for reading leaf area index from initial conditions units = flag dimensions = () type = logical @@ -4833,6 +4833,62 @@ units = flag dimensions = () type = logical +[do_ugwp_v0] + standard_name = do_ugwp_v0 + long_name = flag to activate ver 0 CIRES UGWP + units = flag + dimensions = () + type = logical + intent = in + optional = F +[do_ugwp_v0_orog_only] + standard_name = do_ugwp_v0_orog_only + long_name = flag to activate ver 0 CIRES UGWP - orographic GWD only + units = flag + dimensions = () + type = logical + intent = in + optional = F +[do_gsl_drag_ls_bl] + standard_name = do_gsl_drag_ls_bl + long_name = flag to activate GSL drag suite - large-scale GWD and blocking + units = flag + dimensions = () + type = logical + intent = in + optional = F +[do_gsl_drag_ss] + standard_name = do_gsl_drag_ss + long_name = flag to activate GSL drag suite - small-scale GWD + units = flag + dimensions = () + type = logical + intent = in + optional = F +[do_gsl_drag_tofd] + standard_name = do_gsl_drag_tofd + long_name = flag to activate GSL drag suite - turb orog form drag + units = flag + dimensions = () + type = logical + intent = in + optional = F +[do_ugwp_v1] + standard_name = do_ugwp_v1 + long_name = flag to activate ver 1 CIRES UGWP + units = flag + dimensions = () + type = logical + intent = in + optional = F +[do_ugwp_v1_orog_only] + standard_name = do_ugwp_v1_orog_only + long_name = flag to activate ver 1 CIRES UGWP - orographic GWD only + units = flag + dimensions = () + type = logical + intent = in + optional = F [lmfdeep2] standard_name = flag_for_scale_aware_mass_flux_convection long_name = flag for some scale-aware mass-flux convection scheme active diff --git a/scm/src/gmtb_scm.F90 b/scm/src/gmtb_scm.F90 index ca7fec7f..4660e923 100644 --- a/scm/src/gmtb_scm.F90 +++ b/scm/src/gmtb_scm.F90 @@ -129,6 +129,7 @@ subroutine gmtb_scm_main_sub() cdata%thrd_no = 1 call physics%associate(scm_state) + call physics%set(scm_input, scm_state) ! When asked to calculate 3-dim. tendencies, set Stateout variables to ! Statein variables here in order to capture the first call to dycore diff --git a/scm/src/gmtb_scm_input.F90 b/scm/src/gmtb_scm_input.F90 index 75a227ca..d5c983a4 100644 --- a/scm/src/gmtb_scm_input.F90 +++ b/scm/src/gmtb_scm_input.F90 @@ -10,6 +10,10 @@ module gmtb_scm_input implicit none +integer :: missing_snow_layers = 3 +integer :: missing_soil_layers = 4 +integer :: missing_ice_layers = 2 + contains !> \ingroup SCM @@ -187,15 +191,15 @@ end subroutine get_config_nml !! "processed_case_input" directory. subroutine get_case_init(scm_state, scm_input) use gmtb_scm_type_defs, only : scm_state_type, scm_input_type + use NetCDF_read, only: NetCDF_read_var, check, missing_value type(scm_state_type), intent(in) :: scm_state type(scm_input_type), target, intent(inout) :: scm_input - logical :: noahmp - integer :: input_nlev !< number of levels in the input file integer :: input_nsoil !< number of soil levels in the input file integer :: input_ntimes !< number of times represented in the input file integer :: input_nsnow !< number of snow levels in the input file + integer :: input_nice !< number of sea ice levels in the input file integer :: input_nsoil_plus_nsnow !< number of combined snow and soil levels in the input file ! dimension variables @@ -216,23 +220,31 @@ subroutine get_case_init(scm_state, scm_input) real(kind=dp), allocatable :: input_stc(:) !< soil temperature (K) real(kind=dp), allocatable :: input_smc(:) !< total soil moisture content (fraction) real(kind=dp), allocatable :: input_slc(:) !< liquid soil moisture content (fraction) - real(kind=dp), allocatable :: input_pres_i(:) !< interface pressures - real(kind=dp), allocatable :: input_pres_l(:) !< layer pressures + !real(kind=dp), allocatable :: input_pres_i(:) !< interface pressures + !real(kind=dp), allocatable :: input_pres_l(:) !< layer pressures real(kind=dp), allocatable :: input_snicexy(:) !< snow layer ice (mm) real(kind=dp), allocatable :: input_snliqxy(:) !< snow layer liquid (mm) real(kind=dp), allocatable :: input_tsnoxy(:) !< snow temperature (K) real(kind=dp), allocatable :: input_smoiseq(:) !< equilibrium soil water content (m3 m-3) real(kind=dp), allocatable :: input_zsnsoxy(:) !< layer bottom depth from snow surface (m) + real(kind=dp), allocatable :: input_tiice(:) !< sea ice internal temperature (K) + real(kind=dp), allocatable :: input_tslb(:) !< soil temperature for RUC LSM (K) + real(kind=dp), allocatable :: input_smois(:) !< volume fraction of soil moisture for RUC LSM (frac) + real(kind=dp), allocatable :: input_sh2o(:) !< volume fraction of unfrozen soil moisture for RUC LSM (frac) + real(kind=dp), allocatable :: input_smfr(:) !< volume fraction of frozen soil moisture for RUC LSM (frac) + real(kind=dp), allocatable :: input_flfr(:) !< flag for frozen soil physics + integer :: input_vegsrc !< vegetation source integer :: input_vegtyp !< vegetation type integer :: input_soiltyp!< soil type integer :: input_slopetype !< slope type real(kind=dp) :: input_lat !< column latitude (deg) real(kind=dp) :: input_lon !< column longitude (deg) + real(kind=dp) :: input_tsfco !< input sea surface temperature OR surface skin temperature over land OR surface skin temperature over ice (depending on slmsk) (K) real(kind=dp) :: input_vegfrac !< vegetation fraction real(kind=dp) :: input_shdmin !< minimun vegetation fraction real(kind=dp) :: input_shdmax !< maximun vegetation fraction - real(kind=dp) :: input_zorl !< surfce roughness length [cm] + real(kind=dp) :: input_zorlo !< surfce roughness length over ocean [cm] real(kind=dp) :: input_slmsk !< sea land ice mask [0,1,2] real(kind=dp) :: input_canopy !< amount of water stored in canopy (kg m-2) real(kind=dp) :: input_hice !< sea ice thickness (m) @@ -241,13 +253,28 @@ subroutine get_case_init(scm_state, scm_input) real(kind=dp) :: input_snwdph !< water equivalent snow depth (mm) real(kind=dp) :: input_snoalb !< maximum snow albedo (frac) real(kind=dp) :: input_sncovr !< snow area fraction (frac) - real(kind=dp) :: input_area !< surfce area [m^2] + real(kind=dp) :: input_area !< surface area [m^2] real(kind=dp) :: input_tg3 !< deep soil temperature (K) real(kind=dp) :: input_uustar !< surface friction velocity (m s-1) real(kind=dp) :: input_alvsf !< 60 degree vis albedo with strong cosz dependency real(kind=dp) :: input_alnsf !< 60 degree nir albedo with strong cosz dependency real(kind=dp) :: input_alvwf !< 60 degree vis albedo with weak cosz dependency real(kind=dp) :: input_alnwf !< 60 degree nir albedo with weak cosz dependency + real(kind=dp) :: input_facsf !< fractional coverage with strong cosz dependency + real(kind=dp) :: input_facwf !< fractional coverage with weak cosz dependency + real(kind=dp) :: input_weasd !< water equivalent accumulated snow depth (mm) + real(kind=dp) :: input_f10m !< ratio of sigma level 1 wind and 10m wind + real(kind=dp) :: input_t2m !< 2-meter absolute temperature (K) + real(kind=dp) :: input_q2m !< 2-meter specific humidity (kg kg-1) + real(kind=dp) :: input_ffmm !< Monin-Obukhov similarity function for momentum + real(kind=dp) :: input_ffhh !< Monin-Obukhov similarity function for heat + real(kind=dp) :: input_tprcp !< instantaneous total precipitation amount (m) + real(kind=dp) :: input_srflag !< snow/rain flag for precipitation + real(kind=dp) :: input_tsfcl !< surface skin temperature over land (K) + real(kind=dp) :: input_zorll !< surface roughness length over land (cm) + real(kind=dp) :: input_zorli !< surface roughness length over ice (cm) + real(kind=dp) :: input_zorlw !< surface roughness length from wave model (cm) + real(kind=dp) :: input_stddev !< standard deviation of subgrid orography (m) real(kind=dp) :: input_convexity !< convexity of subgrid orography real(kind=dp) :: input_ol1 !< fraction of grid box with subgrid orography higher than critical height 1 @@ -262,8 +289,12 @@ subroutine get_case_init(scm_state, scm_input) real(kind=dp) :: input_theta !< angle with respect to east of maximum subgrid orographic variations (deg) real(kind=dp) :: input_gamma !< anisotropy of subgrid orography real(kind=dp) :: input_elvmax!< maximum of subgrid orography (m) - real(kind=dp) :: input_facsf !< fractional coverage with strong cosz dependency - real(kind=dp) :: input_facwf !< fractional coverage with weak cosz dependency + real(kind=dp) :: input_oro !< orography (m) + real(kind=dp) :: input_oro_uf !< unfiltered orography (m) + real(kind=dp) :: input_landfrac !< fraction of horizontal grid area occupied by land + real(kind=dp) :: input_lakefrac !< fraction of horizontal grid area occupied by lake + real(kind=dp) :: input_lakedepth !< lake depth (m) + real(kind=dp) :: input_tvxy !< vegetation temperature (K) real(kind=dp) :: input_tgxy !< ground temperature for Noahmp (K) real(kind=dp) :: input_tahxy !< canopy air temperature (K) @@ -294,6 +325,33 @@ subroutine get_case_init(scm_state, scm_input) real(kind=dp) :: input_rechxy !< recharge to or from the water table when shallow (m) real(kind=dp) :: input_snowxy !< number of snow layers + real(kind=dp) :: input_tref !< sea surface reference temperature for NSST (K) + real(kind=dp) :: input_z_c !< sub-layer cooling thickness for NSST (m) + real(kind=dp) :: input_c_0 !< coefficient 1 to calculate d(Tz)/d(Ts) for NSST + real(kind=dp) :: input_c_d !< coefficient 2 to calculate d(Tz)/d(Ts) for NSST + real(kind=dp) :: input_w_0 !< coefficient 3 to calculate d(Tz)/d(Ts) for NSST + real(kind=dp) :: input_w_d !< coefficient 4 to calculate d(Tz)/d(Ts) for NSST + real(kind=dp) :: input_xt !< heat content in diurnal thermocline layer for NSST (K m) + real(kind=dp) :: input_xs !< salinity content in diurnal thermocline layer for NSST (ppt m) + real(kind=dp) :: input_xu !< u-current in diurnal thermocline layer for NSST (m2 s-1) + real(kind=dp) :: input_xv !< v-current in diurnal thermocline layer for NSST (m2 s-1) + real(kind=dp) :: input_xz !< thickness of diurnal thermocline layer for NSST (m) + real(kind=dp) :: input_zm !< thickness of ocean mixed layer for NSST (m) + real(kind=dp) :: input_xtts !< sensitivity of diurnal thermocline layer heat content to surface temperature [d(xt)/d(ts)] for NSST (m) + real(kind=dp) :: input_xzts !< sensitivity of diurnal thermocline layer thickness to surface temperature [d(xz)/d(ts)] for NSST (m K-1) + real(kind=dp) :: input_d_conv !< thickness of free convection layer for NSST (m) + real(kind=dp) :: input_ifd !< index to start DTM run for NSST + real(kind=dp) :: input_dt_cool !< sub-layer cooling amount for NSST (K) + real(kind=dp) :: input_qrain !< sensible heat due to rainfall for NSST (W) + + real(kind=dp) :: input_wetness !< normalized soil wetness for RUC LSM + real(kind=dp) :: input_clw_surf !< cloud condensed water mixing ratio at surface for RUC LSM (kg kg-1) + real(kind=dp) :: input_qwv_surf !< water vapor mixing ratio at surface for RUC LSM (kg kg-1) + real(kind=dp) :: input_tsnow !< snow temperature at the bottom of the first snow layer for RUC LSM (K) + real(kind=dp) :: input_snowfallac !< run-total snow accumulation on the ground for RUC LSM (kg m-2) + real(kind=dp) :: input_acsnow !< snow water equivalent of run-total frozen precip for RUC LSM (kg m-2) + real(kind=dp) :: input_lai !< leaf area index for RUC LSM + !surface time-series variables real(kind=dp), allocatable :: input_pres_surf(:) !< time-series of surface pressure (Pa) real(kind=dp), allocatable :: input_T_surf(:) !< time-series of surface temperature (K) @@ -318,261 +376,241 @@ subroutine get_case_init(scm_state, scm_input) CHARACTER(LEN=nf90_max_name) :: tmpName integer :: ncid, varID, grp_ncid, allocate_status,ierr + real(kind=dp) :: nc_missing_value !> \section get_case_init_alg Algorithm !! @{ !> - Open the case input file found in the processed_case_input dir corresponding to the experiment name. call check(NF90_OPEN(trim(adjustl(scm_state%case_data_dir))//'/'//trim(adjustl(scm_state%case_name))//'.nc',nf90_nowrite,ncid)) - + + !> - Read in missing value from file (replace module variable if present) + ierr = NF90_GET_ATT(ncid, NF90_GLOBAL, 'missing_value', nc_missing_value) + if(ierr == NF90_NOERR) then + missing_value = nc_missing_value + end if + !> - Get the dimensions (global group). + !required dimensions call check(NF90_INQ_DIMID(ncid,"levels",varID)) call check(NF90_INQUIRE_DIMENSION(ncid, varID, tmpName, input_nlev)) - if (scm_state%model_ics) then - call check(NF90_INQ_DIMID(ncid,"nsoil",varID)) - call check(NF90_INQUIRE_DIMENSION(ncid, varID, tmpName, input_nsoil)) - ierr = NF90_INQ_DIMID(ncid,"nsnow",varID) - noahmp = .false. - if (ierr == 0) then - call check(NF90_INQUIRE_DIMENSION(ncid, varID, tmpName, input_nsnow)) - noahmp = .true. !if the nsnow variable is present, NoahMP ICs should be present - endif - endif call check(NF90_INQ_DIMID(ncid,"time",varID)) call check(NF90_INQUIRE_DIMENSION(ncid, varID, tmpName, input_ntimes)) - + + !possible dimensions (if using model ICs) + ierr = NF90_INQ_DIMID(ncid,"nsoil",varID) + if(ierr /= NF90_NOERR) then + input_nsoil = missing_soil_layers + else + call check(NF90_INQUIRE_DIMENSION(ncid, varID, tmpName, input_nsoil)) + end if + ierr = NF90_INQ_DIMID(ncid,"nsnow",varID) + if(ierr /= NF90_NOERR) then + input_nsnow = missing_snow_layers + else + call check(NF90_INQUIRE_DIMENSION(ncid, varID, tmpName, input_nsnow)) + end if + ierr = NF90_INQ_DIMID(ncid,"nice",varID) + if(ierr /= NF90_NOERR) then + input_nice = missing_ice_layers + else + call check(NF90_INQUIRE_DIMENSION(ncid, varID, tmpName, input_nice)) + end if + !> - Allocate the dimension variables. allocate(input_pres(input_nlev),input_time(input_ntimes), stat=allocate_status) - !> - Read in the dimension variables. - call check(NF90_INQ_VARID(ncid,"levels",varID)) - call check(NF90_GET_VAR(ncid,varID,input_pres)) - call check(NF90_INQ_VARID(ncid,"time",varID)) - call check(NF90_GET_VAR(ncid,varID,input_time)) + !> - Read in the dimension variables (required). + call NetCDF_read_var(ncid, "levels", .True., input_pres) + call NetCDF_read_var(ncid, "time", .True., input_time) !> - Read in the initial conditions. !> - Find group ncid for initial group. call check(NF90_INQ_GRP_NCID(ncid,"initial",grp_ncid)) - !> - Allocate the initial profiles. - allocate(input_thetail(input_nlev), input_qt(input_nlev), input_ql(input_nlev), input_qi(input_nlev), & + !> - Allocate the initial profiles (required). One of thetail or temp is required. + allocate(input_thetail(input_nlev), input_temp(input_nlev), input_qt(input_nlev), input_ql(input_nlev), input_qi(input_nlev), & input_u(input_nlev), input_v(input_nlev), input_tke(input_nlev), input_ozone(input_nlev), stat=allocate_status) - if (scm_state%model_ics) then - allocate(input_stc(input_nsoil), input_temp(input_nlev),input_smc(input_nsoil), input_slc(input_nsoil), & - input_pres_i(input_nlev+1),input_pres_l(input_nlev), stat=allocate_status) - input_pres_i(:) = -999.9 - input_pres_l(:) = -999.9 - if (noahmp) then - allocate(input_snicexy(input_nsnow), input_snliqxy(input_nsnow), input_tsnoxy(input_nsnow), & - input_smoiseq(input_nsoil), input_zsnsoxy(input_nsnow + input_nsoil)) - endif - endif !> - Read in the initial profiles. The variable names in all input files are expected to be identical. - if (.NOT. scm_state%model_ics) then - call check(NF90_INQ_VARID(grp_ncid,"thetail",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_thetail)) - endif - call check(NF90_INQ_VARID(grp_ncid,"qt",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_qt)) - call check(NF90_INQ_VARID(grp_ncid,"ql",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_ql)) - call check(NF90_INQ_VARID(grp_ncid,"qi",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_qi)) - call check(NF90_INQ_VARID(grp_ncid,"u",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_u)) - call check(NF90_INQ_VARID(grp_ncid,"v",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_v)) - call check(NF90_INQ_VARID(grp_ncid,"tke",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_tke)) - call check(NF90_INQ_VARID(grp_ncid,"ozone",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_ozone)) - if (scm_state%model_ics) then - call check(NF90_INQ_VARID(grp_ncid,"temp",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_temp)) - call check(NF90_INQ_VARID(grp_ncid,"stc",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_stc)) - call check(NF90_INQ_VARID(grp_ncid,"smc",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_smc)) - call check(NF90_INQ_VARID(grp_ncid,"slc",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_slc)) - ierr = NF90_INQ_VARID(grp_ncid,"pres_i",varID) - if (ierr.EQ.0) then ! input file should have pres_i and pres_l - call check(NF90_GET_VAR(grp_ncid,varID,input_pres_i)) - call check(NF90_INQ_VARID(grp_ncid,"pres_l",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_pres_l)) - endif - if (noahmp) then - call check(NF90_INQ_VARID(grp_ncid,"snicexy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_snicexy)) - call check(NF90_INQ_VARID(grp_ncid,"snliqxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_snliqxy)) - call check(NF90_INQ_VARID(grp_ncid,"tsnoxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_tsnoxy)) - call check(NF90_INQ_VARID(grp_ncid,"smoiseq",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_smoiseq)) - call check(NF90_INQ_VARID(grp_ncid,"zsnsoxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_zsnsoxy)) - endif - endif + + !Either thetail or T must be present + call NetCDF_read_var(grp_ncid, "thetail", .False., input_thetail) + call NetCDF_read_var(grp_ncid, "temp", .False., input_temp) + if (maxval(input_thetail) < 0 .and. maxval(input_temp) < 0) then + write(*,*) "One of thetail or temp variables must be present in ",trim(adjustl(scm_state%case_name))//'.nc',". Stopping..." + STOP + end if + call NetCDF_read_var(grp_ncid, "qt", .True., input_qt ) + call NetCDF_read_var(grp_ncid, "ql", .True., input_ql ) + call NetCDF_read_var(grp_ncid, "qi", .True., input_qi ) + call NetCDF_read_var(grp_ncid, "u", .True., input_u ) + call NetCDF_read_var(grp_ncid, "v", .True., input_v ) + call NetCDF_read_var(grp_ncid, "tke", .True., input_tke ) + call NetCDF_read_var(grp_ncid, "ozone", .True., input_ozone) + + !possible initial profiles + !needed for Noah LSM and others (when running with model ICs) + allocate(input_stc(input_nsoil), input_smc(input_nsoil), input_slc(input_nsoil), & + stat=allocate_status) + call NetCDF_read_var(grp_ncid, "stc", .False., input_stc) + call NetCDF_read_var(grp_ncid, "smc", .False., input_smc) + call NetCDF_read_var(grp_ncid, "slc", .False., input_slc) + + !needed for NoahMP LSM (when running with model ICs) + allocate(input_snicexy(input_nsnow), input_snliqxy(input_nsnow), input_tsnoxy(input_nsnow), & + input_smoiseq(input_nsoil), input_zsnsoxy(input_nsnow + input_nsoil)) + call NetCDF_read_var(grp_ncid, "snicexy", .False., input_snicexy) + call NetCDF_read_var(grp_ncid, "snliqxy", .False., input_snliqxy) + call NetCDF_read_var(grp_ncid, "tsnoxy", .False., input_tsnoxy ) + call NetCDF_read_var(grp_ncid, "smoiseq", .False., input_smoiseq) + call NetCDF_read_var(grp_ncid, "zsnsoxy", .False., input_zsnsoxy) + + !needed for fractional grid (when running with model ICs) + allocate(input_tiice(input_nice)) + call NetCDF_read_var(grp_ncid, "tiice", .False., input_tiice) + + !needed for RUC LSM (when running with model ICs) + allocate(input_tslb(input_nsoil), input_smois(input_nsoil), input_sh2o(input_nsoil), & + input_smfr(input_nsoil), input_flfr(input_nsoil)) + call NetCDF_read_var(grp_ncid, "tslb", .False., input_tslb ) + call NetCDF_read_var(grp_ncid, "smois", .False., input_smois) + call NetCDF_read_var(grp_ncid, "sh2o", .False., input_sh2o ) + call NetCDF_read_var(grp_ncid, "smfr", .False., input_smfr ) + call NetCDF_read_var(grp_ncid, "flfr", .False., input_flfr ) !> - Find group ncid for scalar group. call check(NF90_INQ_GRP_NCID(ncid,"scalars",grp_ncid)) - ! ierr = NF90_INQ_VARID(grp_ncid,"lat",varID) - ! if (ierr .EQ. 0) then - ! call check(NF90_GET_VAR(grp_ncid,varID,input_lat)) - ! call check(NF90_INQ_VARID(grp_ncid,"lon",varID)) - ! call check(NF90_GET_VAR(grp_ncid,varID,input_lon)) - ! endif - call check(NF90_INQ_VARID(grp_ncid,"lat",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_lat)) - call check(NF90_INQ_VARID(grp_ncid,"lon",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_lon)) - - if (scm_state%model_ics) then - call check(NF90_INQ_VARID(grp_ncid,"vegsrc",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_vegsrc)) - call check(NF90_INQ_VARID(grp_ncid,"vegtyp",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_vegtyp)) - call check(NF90_INQ_VARID(grp_ncid,"soiltyp",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_soiltyp)) - call check(NF90_INQ_VARID(grp_ncid,"slopetyp",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_slopetype)) - call check(NF90_INQ_VARID(grp_ncid,"vegfrac",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_vegfrac)) - call check(NF90_INQ_VARID(grp_ncid,"shdmin",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_shdmin)) - call check(NF90_INQ_VARID(grp_ncid,"shdmax",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_shdmax)) - call check(NF90_INQ_VARID(grp_ncid,"zorl",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_zorl)) - call check(NF90_INQ_VARID(grp_ncid,"slmsk",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_slmsk)) - call check(NF90_INQ_VARID(grp_ncid,"canopy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_canopy)) - call check(NF90_INQ_VARID(grp_ncid,"hice",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_hice)) - call check(NF90_INQ_VARID(grp_ncid,"fice",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_fice)) - call check(NF90_INQ_VARID(grp_ncid,"tisfc",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_tisfc)) - call check(NF90_INQ_VARID(grp_ncid,"snwdph",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_snwdph)) - call check(NF90_INQ_VARID(grp_ncid,"snoalb",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_snoalb)) - call check(NF90_INQ_VARID(grp_ncid,"sncovr",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_sncovr)) - call check(NF90_INQ_VARID(grp_ncid,"area",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_area)) - call check(NF90_INQ_VARID(grp_ncid,"tg3",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_tg3)) - call check(NF90_INQ_VARID(grp_ncid,"uustar",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_uustar)) - call check(NF90_INQ_VARID(grp_ncid,"alvsf",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_alvsf)) - call check(NF90_INQ_VARID(grp_ncid,"alnsf",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_alnsf)) - call check(NF90_INQ_VARID(grp_ncid,"alvwf",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_alvwf)) - call check(NF90_INQ_VARID(grp_ncid,"alnwf",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_alnwf)) - call check(NF90_INQ_VARID(grp_ncid,"stddev",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_stddev)) - call check(NF90_INQ_VARID(grp_ncid,"convexity",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_convexity)) - call check(NF90_INQ_VARID(grp_ncid,"oa1",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_oa1)) - call check(NF90_INQ_VARID(grp_ncid,"oa2",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_oa2)) - call check(NF90_INQ_VARID(grp_ncid,"oa3",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_oa3)) - call check(NF90_INQ_VARID(grp_ncid,"oa4",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_oa4)) - call check(NF90_INQ_VARID(grp_ncid,"ol1",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_ol1)) - call check(NF90_INQ_VARID(grp_ncid,"ol2",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_ol2)) - call check(NF90_INQ_VARID(grp_ncid,"ol3",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_ol3)) - call check(NF90_INQ_VARID(grp_ncid,"ol4",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_ol4)) - call check(NF90_INQ_VARID(grp_ncid,"sigma",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_sigma)) - call check(NF90_INQ_VARID(grp_ncid,"gamma",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_gamma)) - call check(NF90_INQ_VARID(grp_ncid,"theta",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_theta)) - call check(NF90_INQ_VARID(grp_ncid,"elvmax",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_elvmax)) - call check(NF90_INQ_VARID(grp_ncid,"facsf",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_facsf)) - call check(NF90_INQ_VARID(grp_ncid,"facwf",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_facwf)) - if (noahmp) then - call check(NF90_INQ_VARID(grp_ncid,"tvxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_tvxy)) - call check(NF90_INQ_VARID(grp_ncid,"tgxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_tgxy)) - call check(NF90_INQ_VARID(grp_ncid,"tahxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_tahxy)) - call check(NF90_INQ_VARID(grp_ncid,"canicexy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_canicexy)) - call check(NF90_INQ_VARID(grp_ncid,"canliqxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_canliqxy)) - call check(NF90_INQ_VARID(grp_ncid,"eahxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_eahxy)) - call check(NF90_INQ_VARID(grp_ncid,"cmxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_cmxy)) - call check(NF90_INQ_VARID(grp_ncid,"chxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_chxy)) - call check(NF90_INQ_VARID(grp_ncid,"fwetxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_fwetxy)) - call check(NF90_INQ_VARID(grp_ncid,"sneqvoxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_sneqvoxy)) - call check(NF90_INQ_VARID(grp_ncid,"alboldxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_alboldxy)) - call check(NF90_INQ_VARID(grp_ncid,"qsnowxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_qsnowxy)) - call check(NF90_INQ_VARID(grp_ncid,"wslakexy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_wslakexy)) - call check(NF90_INQ_VARID(grp_ncid,"taussxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_taussxy)) - call check(NF90_INQ_VARID(grp_ncid,"waxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_waxy)) - call check(NF90_INQ_VARID(grp_ncid,"wtxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_wtxy)) - call check(NF90_INQ_VARID(grp_ncid,"zwtxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_zwtxy)) - call check(NF90_INQ_VARID(grp_ncid,"xlaixy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_xlaixy)) - call check(NF90_INQ_VARID(grp_ncid,"xsaixy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_xsaixy)) - call check(NF90_INQ_VARID(grp_ncid,"lfmassxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_lfmassxy)) - call check(NF90_INQ_VARID(grp_ncid,"stmassxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_stmassxy)) - call check(NF90_INQ_VARID(grp_ncid,"rtmassxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_rtmassxy)) - call check(NF90_INQ_VARID(grp_ncid,"woodxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_woodxy)) - call check(NF90_INQ_VARID(grp_ncid,"stblcpxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_stblcpxy)) - call check(NF90_INQ_VARID(grp_ncid,"fastcpxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_fastcpxy)) - call check(NF90_INQ_VARID(grp_ncid,"smcwtdxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_smcwtdxy)) - call check(NF90_INQ_VARID(grp_ncid,"deeprechxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_deeprechxy)) - call check(NF90_INQ_VARID(grp_ncid,"rechxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_rechxy)) - call check(NF90_INQ_VARID(grp_ncid,"snowxy",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_snowxy)) - endif - endif - + !required + call NetCDF_read_var(grp_ncid, "lat", .True., input_lat) + call NetCDF_read_var(grp_ncid, "lon", .True., input_lon) + !time data and area in file ignored? + call NetCDF_read_var(grp_ncid, "area", .False., input_area) + + !possible scalars + !Noah LSM parameters (when running with model ICs) + call NetCDF_read_var(grp_ncid, "vegsrc", .False., input_vegsrc ) + call NetCDF_read_var(grp_ncid, "vegtyp", .False., input_vegtyp ) + call NetCDF_read_var(grp_ncid, "soiltyp", .False., input_soiltyp ) + call NetCDF_read_var(grp_ncid, "slopetyp", .False., input_slopetype) + call NetCDF_read_var(grp_ncid, "tsfco", .False., input_tsfco) + call NetCDF_read_var(grp_ncid, "vegfrac", .False., input_vegfrac) + call NetCDF_read_var(grp_ncid, "shdmin", .False., input_shdmin) + call NetCDF_read_var(grp_ncid, "shdmax", .False., input_shdmax) + call NetCDF_read_var(grp_ncid, "zorlo", .False., input_zorlo) + call NetCDF_read_var(grp_ncid, "slmsk", .False., input_slmsk) + call NetCDF_read_var(grp_ncid, "canopy", .False., input_canopy) + call NetCDF_read_var(grp_ncid, "hice", .False., input_hice) + call NetCDF_read_var(grp_ncid, "fice", .False., input_fice) + call NetCDF_read_var(grp_ncid, "tisfc", .False., input_tisfc) + call NetCDF_read_var(grp_ncid, "snwdph", .False., input_snwdph) + call NetCDF_read_var(grp_ncid, "snoalb", .False., input_snoalb) + call NetCDF_read_var(grp_ncid, "tg3", .False., input_tg3) + call NetCDF_read_var(grp_ncid, "uustar", .False., input_uustar) + call NetCDF_read_var(grp_ncid, "alvsf", .False., input_alvsf) + call NetCDF_read_var(grp_ncid, "alnsf", .False., input_alnsf) + call NetCDF_read_var(grp_ncid, "alvwf", .False., input_alvwf) + call NetCDF_read_var(grp_ncid, "alnwf", .False., input_alnwf) + call NetCDF_read_var(grp_ncid, "facsf", .False., input_facsf) + call NetCDF_read_var(grp_ncid, "facwf", .False., input_facwf) + call NetCDF_read_var(grp_ncid, "weasd", .False., input_weasd) + call NetCDF_read_var(grp_ncid, "f10m", .False., input_f10m) + call NetCDF_read_var(grp_ncid, "t2m", .False., input_t2m) + call NetCDF_read_var(grp_ncid, "q2m", .False., input_q2m) + call NetCDF_read_var(grp_ncid, "ffmm", .False., input_ffmm) + call NetCDF_read_var(grp_ncid, "ffhh", .False., input_ffhh) + call NetCDF_read_var(grp_ncid, "tprcp", .False., input_tprcp) + call NetCDF_read_var(grp_ncid, "srflag", .False., input_srflag) + call NetCDF_read_var(grp_ncid, "sncovr", .False., input_sncovr) + call NetCDF_read_var(grp_ncid, "tsfcl", .False., input_tsfcl) + call NetCDF_read_var(grp_ncid, "zorll", .False., input_zorll) + call NetCDF_read_var(grp_ncid, "zorli", .False., input_zorli) + call NetCDF_read_var(grp_ncid, "zorlw", .False., input_zorlw) + + !orographic parameters + call NetCDF_read_var(grp_ncid, "stddev", .False., input_stddev) + call NetCDF_read_var(grp_ncid, "convexity", .False., input_convexity) + call NetCDF_read_var(grp_ncid, "oa1", .False., input_oa1) + call NetCDF_read_var(grp_ncid, "oa2", .False., input_oa2) + call NetCDF_read_var(grp_ncid, "oa3", .False., input_oa3) + call NetCDF_read_var(grp_ncid, "oa4", .False., input_oa4) + call NetCDF_read_var(grp_ncid, "ol1", .False., input_ol1) + call NetCDF_read_var(grp_ncid, "ol2", .False., input_ol2) + call NetCDF_read_var(grp_ncid, "ol3", .False., input_ol3) + call NetCDF_read_var(grp_ncid, "ol4", .False., input_ol4) + call NetCDF_read_var(grp_ncid, "theta", .False., input_theta) + call NetCDF_read_var(grp_ncid, "gamma", .False., input_gamma) + call NetCDF_read_var(grp_ncid, "sigma", .False., input_sigma) + call NetCDF_read_var(grp_ncid, "elvmax", .False., input_elvmax) + call NetCDF_read_var(grp_ncid, "oro", .False., input_oro) + call NetCDF_read_var(grp_ncid, "oro_uf", .False., input_oro_uf) + call NetCDF_read_var(grp_ncid, "landfrac", .False., input_landfrac) + call NetCDF_read_var(grp_ncid, "lakefrac", .False., input_lakefrac) + call NetCDF_read_var(grp_ncid, "lakedepth", .False., input_lakedepth) + + !NoahMP parameters + call NetCDF_read_var(grp_ncid, "tvxy", .False., input_tvxy) + call NetCDF_read_var(grp_ncid, "tgxy", .False., input_tgxy) + call NetCDF_read_var(grp_ncid, "tahxy", .False., input_tahxy) + call NetCDF_read_var(grp_ncid, "canicexy", .False., input_canicexy) + call NetCDF_read_var(grp_ncid, "canliqxy", .False., input_canliqxy) + call NetCDF_read_var(grp_ncid, "eahxy", .False., input_eahxy) + call NetCDF_read_var(grp_ncid, "cmxy", .False., input_cmxy) + call NetCDF_read_var(grp_ncid, "chxy", .False., input_chxy) + call NetCDF_read_var(grp_ncid, "fwetxy", .False., input_fwetxy) + call NetCDF_read_var(grp_ncid, "sneqvoxy", .False., input_sneqvoxy) + call NetCDF_read_var(grp_ncid, "alboldxy", .False., input_alboldxy) + call NetCDF_read_var(grp_ncid, "qsnowxy", .False., input_qsnowxy) + call NetCDF_read_var(grp_ncid, "wslakexy", .False., input_wslakexy) + call NetCDF_read_var(grp_ncid, "taussxy", .False., input_taussxy) + call NetCDF_read_var(grp_ncid, "waxy", .False., input_waxy) + call NetCDF_read_var(grp_ncid, "wtxy", .False., input_wtxy) + call NetCDF_read_var(grp_ncid, "zwtxy", .False., input_zwtxy) + call NetCDF_read_var(grp_ncid, "xlaixy", .False., input_xlaixy) + call NetCDF_read_var(grp_ncid, "xsaixy", .False., input_xsaixy) + call NetCDF_read_var(grp_ncid, "lfmassxy", .False., input_lfmassxy) + call NetCDF_read_var(grp_ncid, "stmassxy", .False., input_stmassxy) + call NetCDF_read_var(grp_ncid, "rtmassxy", .False., input_rtmassxy) + call NetCDF_read_var(grp_ncid, "woodxy", .False., input_woodxy) + call NetCDF_read_var(grp_ncid, "stblcpxy", .False., input_stblcpxy) + call NetCDF_read_var(grp_ncid, "fastcpxy", .False., input_fastcpxy) + call NetCDF_read_var(grp_ncid, "smcwtdxy", .False., input_smcwtdxy) + call NetCDF_read_var(grp_ncid, "deeprechxy",.False., input_deeprechxy) + call NetCDF_read_var(grp_ncid, "rechxy", .False., input_rechxy) + call NetCDF_read_var(grp_ncid, "snowxy", .False., input_snowxy) + + !NSST variables + call NetCDF_read_var(grp_ncid, "tref", .False., input_tref) + call NetCDF_read_var(grp_ncid, "z_c", .False., input_z_c) + call NetCDF_read_var(grp_ncid, "c_0", .False., input_c_0) + call NetCDF_read_var(grp_ncid, "c_d", .False., input_c_d) + call NetCDF_read_var(grp_ncid, "w_0", .False., input_w_0) + call NetCDF_read_var(grp_ncid, "w_d", .False., input_w_d) + call NetCDF_read_var(grp_ncid, "xt", .False., input_xt) + call NetCDF_read_var(grp_ncid, "xs", .False., input_xs) + call NetCDF_read_var(grp_ncid, "xu", .False., input_xu) + call NetCDF_read_var(grp_ncid, "xv", .False., input_xv) + call NetCDF_read_var(grp_ncid, "xz", .False., input_xz) + call NetCDF_read_var(grp_ncid, "zm", .False., input_zm) + call NetCDF_read_var(grp_ncid, "xtts", .False., input_xtts) + call NetCDF_read_var(grp_ncid, "xzts", .False., input_xzts) + call NetCDF_read_var(grp_ncid, "d_conv", .False., input_d_conv) + call NetCDF_read_var(grp_ncid, "ifd", .False., input_ifd) + call NetCDF_read_var(grp_ncid, "dt_cool", .False., input_dt_cool) + call NetCDF_read_var(grp_ncid, "qrain", .False., input_qrain) + + !RUC LSM variables + call NetCDF_read_var(grp_ncid, "wetness", .False., input_wetness) + call NetCDF_read_var(grp_ncid, "clw_surf", .False., input_clw_surf) + call NetCDF_read_var(grp_ncid, "qwv_surf", .False., input_qwv_surf) + call NetCDF_read_var(grp_ncid, "tsnow", .False., input_tsnow) + call NetCDF_read_var(grp_ncid, "snowfall_acc", .False., input_snowfallac) + call NetCDF_read_var(grp_ncid, "swe_snowfall_acc", .False., input_acsnow) + call NetCDF_read_var(grp_ncid, "lai", .False., input_lai) + !> - Read in the forcing data. !> - Find group ncid for forcing group. @@ -591,65 +629,38 @@ subroutine get_case_init(scm_state, scm_input) stat=allocate_status) !> - Read in the time-series and 2D forcing data. - call check(NF90_INQ_VARID(grp_ncid,"p_surf",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_pres_surf)) - call check(NF90_INQ_VARID(grp_ncid,"T_surf",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_T_surf)) - if(scm_state%sfc_flux_spec) then - call check(NF90_INQ_VARID(grp_ncid,"sh_flux_sfc",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_sh_flux_sfc)) - call check(NF90_INQ_VARID(grp_ncid,"lh_flux_sfc",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_lh_flux_sfc)) - else - input_sh_flux_sfc = 0.0 - input_lh_flux_sfc = 0.0 - end if - call check(NF90_INQ_VARID(grp_ncid,"w_ls",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_w_ls)) - call check(NF90_INQ_VARID(grp_ncid,"omega",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_omega)) - call check(NF90_INQ_VARID(grp_ncid,"u_g",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_u_g)) - call check(NF90_INQ_VARID(grp_ncid,"v_g",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_v_g)) - call check(NF90_INQ_VARID(grp_ncid,"u_nudge",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_u_nudge)) - call check(NF90_INQ_VARID(grp_ncid,"v_nudge",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_v_nudge)) - call check(NF90_INQ_VARID(grp_ncid,"T_nudge",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_T_nudge)) - call check(NF90_INQ_VARID(grp_ncid,"thil_nudge",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_thil_nudge)) - call check(NF90_INQ_VARID(grp_ncid,"qt_nudge",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_qt_nudge)) - call check(NF90_INQ_VARID(grp_ncid,"dT_dt_rad",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_dT_dt_rad)) - call check(NF90_INQ_VARID(grp_ncid,"h_advec_thetail",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_h_advec_thetail)) - call check(NF90_INQ_VARID(grp_ncid,"h_advec_qt",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_h_advec_qt)) - call check(NF90_INQ_VARID(grp_ncid,"v_advec_thetail",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_v_advec_thetail)) - call check(NF90_INQ_VARID(grp_ncid,"v_advec_qt",varID)) - call check(NF90_GET_VAR(grp_ncid,varID,input_v_advec_qt)) + call NetCDF_read_var(grp_ncid, "p_surf", .True., input_pres_surf) + call NetCDF_read_var(grp_ncid, "T_surf", .True., input_T_surf) + call NetCDF_read_var(grp_ncid, "sh_flux_sfc", .False., input_sh_flux_sfc) + call NetCDF_read_var(grp_ncid, "lh_flux_sfc", .False., input_lh_flux_sfc) + + call NetCDF_read_var(grp_ncid, "w_ls", .True., input_w_ls) + call NetCDF_read_var(grp_ncid, "omega", .True., input_omega) + call NetCDF_read_var(grp_ncid, "u_g", .True., input_u_g) + call NetCDF_read_var(grp_ncid, "v_g", .True., input_v_g) + call NetCDF_read_var(grp_ncid, "u_nudge", .True., input_u_nudge) + call NetCDF_read_var(grp_ncid, "v_nudge", .True., input_v_nudge) + call NetCDF_read_var(grp_ncid, "T_nudge", .True., input_T_nudge) + call NetCDF_read_var(grp_ncid, "thil_nudge", .True., input_thil_nudge) + call NetCDF_read_var(grp_ncid, "qt_nudge", .True., input_qt_nudge) + call NetCDF_read_var(grp_ncid, "dT_dt_rad", .True., input_dT_dt_rad) + call NetCDF_read_var(grp_ncid, "h_advec_thetail", .True., input_h_advec_thetail) + call NetCDF_read_var(grp_ncid, "h_advec_qt", .True., input_h_advec_qt) + call NetCDF_read_var(grp_ncid, "v_advec_thetail", .True., input_v_advec_thetail) + call NetCDF_read_var(grp_ncid, "v_advec_qt", .True., input_v_advec_qt) + call check(NF90_CLOSE(NCID=ncid)) - call scm_input%create(input_ntimes, input_nlev) - if (scm_state%model_ics) then - call scm_input%create_modelics(input_nsoil,input_nsnow,input_nlev,noahmp) - endif - + call scm_input%create(input_ntimes, input_nlev, input_nsoil, input_nsnow, input_nice) + ! GJF already done in scm_input%create routine !scm_input%input_nlev = input_nlev !scm_input%input_ntimes = input_ntimes scm_input%input_pres = input_pres scm_input%input_time = input_time - if (scm_state%model_ics) then - scm_input%input_temp = input_temp - else - scm_input%input_thetail = input_thetail - endif + scm_input%input_temp = input_temp + scm_input%input_thetail = input_thetail scm_input%input_qt = input_qt scm_input%input_ql = input_ql scm_input%input_qi = input_qi @@ -677,96 +688,148 @@ subroutine get_case_init(scm_state, scm_input) scm_input%input_T_nudge = input_T_nudge scm_input%input_thil_nudge = input_thil_nudge scm_input%input_qt_nudge = input_qt_nudge - if (scm_state%model_ics) then - scm_input%input_stc = input_stc - scm_input%input_smc = input_smc - scm_input%input_slc = input_slc - scm_input%input_vegsrc = input_vegsrc - scm_input%input_vegtyp = REAL(input_vegtyp, kind=dp) - scm_input%input_soiltyp = REAL(input_soiltyp, kind=dp) - scm_input%input_slopetype = REAL(input_slopetype, kind=dp) - scm_input%input_vegfrac = input_vegfrac - scm_input%input_shdmin = input_shdmin - scm_input%input_shdmax = input_shdmax - scm_input%input_zorl = input_zorl - scm_input%input_slmsk = input_slmsk - scm_input%input_canopy = input_canopy - scm_input%input_hice = input_hice - scm_input%input_fice = input_fice - scm_input%input_tisfc = input_tisfc - scm_input%input_snwdph = input_snwdph - scm_input%input_snoalb = input_snoalb - scm_input%input_sncovr = input_sncovr - scm_input%input_area = input_area - scm_input%input_tg3 = input_tg3 - scm_input%input_uustar = input_uustar - scm_input%input_alvsf = input_alvsf - scm_input%input_alnsf = input_alnsf - scm_input%input_alvwf = input_alvwf - scm_input%input_alnwf = input_alnwf - scm_input%input_convexity= input_convexity - scm_input%input_stddev = input_stddev - scm_input%input_oa1 = input_oa1 - scm_input%input_oa2 = input_oa2 - scm_input%input_oa3 = input_oa3 - scm_input%input_oa4 = input_oa4 - scm_input%input_ol1 = input_ol1 - scm_input%input_ol2 = input_ol2 - scm_input%input_ol3 = input_ol3 - scm_input%input_ol4 = input_ol4 - scm_input%input_sigma = input_sigma - scm_input%input_theta = input_theta - scm_input%input_gamma = input_gamma - scm_input%input_elvmax = input_elvmax - scm_input%input_facsf = input_facsf - scm_input%input_facwf = input_facwf - scm_input%input_pres_i = input_pres_i - scm_input%input_pres_l = input_pres_l - if (noahmp) then - scm_input%input_snicexy = input_snicexy - scm_input%input_snliqxy = input_snliqxy - scm_input%input_tsnoxy = input_tsnoxy - scm_input%input_smoiseq = input_smoiseq - scm_input%input_zsnsoxy = input_zsnsoxy - scm_input%input_tvxy = input_tvxy - scm_input%input_tgxy = input_tgxy - scm_input%input_tahxy = input_tahxy - scm_input%input_canicexy = input_canicexy - scm_input%input_canliqxy = input_canliqxy - scm_input%input_eahxy = input_eahxy - scm_input%input_cmxy = input_cmxy - scm_input%input_chxy = input_chxy - scm_input%input_fwetxy = input_fwetxy - scm_input%input_sneqvoxy = input_sneqvoxy - scm_input%input_alboldxy = input_alboldxy - scm_input%input_qsnowxy = input_qsnowxy - scm_input%input_wslakexy = input_wslakexy - scm_input%input_taussxy = input_taussxy - scm_input%input_waxy = input_waxy - scm_input%input_wtxy = input_wtxy - scm_input%input_zwtxy = input_zwtxy - scm_input%input_xlaixy = input_xlaixy - scm_input%input_xsaixy = input_xsaixy - scm_input%input_lfmassxy = input_lfmassxy - scm_input%input_stmassxy = input_stmassxy - scm_input%input_rtmassxy = input_rtmassxy - scm_input%input_woodxy = input_woodxy - scm_input%input_stblcpxy = input_stblcpxy - scm_input%input_fastcpxy = input_fastcpxy - scm_input%input_smcwtdxy = input_smcwtdxy - scm_input%input_deeprechxy = input_deeprechxy - scm_input%input_rechxy = input_rechxy - scm_input%input_snowxy = input_snowxy - endif - endif - + + scm_input%input_stc = input_stc + scm_input%input_smc = input_smc + scm_input%input_slc = input_slc + + scm_input%input_snicexy = input_snicexy + scm_input%input_snliqxy = input_snliqxy + scm_input%input_tsnoxy = input_tsnoxy + scm_input%input_smoiseq = input_smoiseq + scm_input%input_zsnsoxy = input_zsnsoxy + + scm_input%input_tiice = input_tiice + scm_input%input_tslb = input_tslb + scm_input%input_smois = input_smois + scm_input%input_sh2o = input_sh2o + scm_input%input_smfr = input_smfr + scm_input%input_flfr = input_flfr + + scm_input%input_vegsrc = input_vegsrc + scm_input%input_vegtyp = REAL(input_vegtyp, kind=dp) + scm_input%input_soiltyp = REAL(input_soiltyp, kind=dp) + scm_input%input_slopetype = REAL(input_slopetype, kind=dp) + scm_input%input_tsfco = input_tsfco + scm_input%input_vegfrac = input_vegfrac + scm_input%input_shdmin = input_shdmin + scm_input%input_shdmax = input_shdmax + scm_input%input_zorlo = input_zorlo + scm_input%input_slmsk = input_slmsk + scm_input%input_canopy = input_canopy + scm_input%input_hice = input_hice + scm_input%input_fice = input_fice + scm_input%input_tisfc = input_tisfc + scm_input%input_snwdph = input_snwdph + scm_input%input_snoalb = input_snoalb + scm_input%input_sncovr = input_sncovr + scm_input%input_area = input_area + scm_input%input_tg3 = input_tg3 + scm_input%input_uustar = input_uustar + scm_input%input_alvsf = input_alvsf + scm_input%input_alnsf = input_alnsf + scm_input%input_alvwf = input_alvwf + scm_input%input_alnwf = input_alnwf + scm_input%input_facsf = input_facsf + scm_input%input_facwf = input_facwf + scm_input%input_weasd = input_weasd + scm_input%input_f10m = input_f10m + scm_input%input_t2m = input_t2m + scm_input%input_q2m = input_q2m + scm_input%input_ffmm = input_ffmm + scm_input%input_ffhh = input_ffhh + scm_input%input_tprcp = input_tprcp + scm_input%input_srflag = input_srflag + scm_input%input_tsfcl = input_tsfcl + scm_input%input_zorll = input_zorll + scm_input%input_zorli = input_zorli + scm_input%input_zorlw = input_zorlw + + scm_input%input_stddev = input_stddev + scm_input%input_convexity= input_convexity + scm_input%input_oa1 = input_oa1 + scm_input%input_oa2 = input_oa2 + scm_input%input_oa3 = input_oa3 + scm_input%input_oa4 = input_oa4 + scm_input%input_ol1 = input_ol1 + scm_input%input_ol2 = input_ol2 + scm_input%input_ol3 = input_ol3 + scm_input%input_ol4 = input_ol4 + scm_input%input_sigma = input_sigma + scm_input%input_theta = input_theta + scm_input%input_gamma = input_gamma + scm_input%input_elvmax = input_elvmax + scm_input%input_oro = input_oro + scm_input%input_oro_uf = input_oro_uf + scm_input%input_landfrac = input_landfrac + scm_input%input_lakefrac = input_lakefrac + scm_input%input_lakedepth= input_lakedepth + + scm_input%input_tvxy = input_tvxy + scm_input%input_tgxy = input_tgxy + scm_input%input_tahxy = input_tahxy + scm_input%input_canicexy = input_canicexy + scm_input%input_canliqxy = input_canliqxy + scm_input%input_eahxy = input_eahxy + scm_input%input_cmxy = input_cmxy + scm_input%input_chxy = input_chxy + scm_input%input_fwetxy = input_fwetxy + scm_input%input_sneqvoxy = input_sneqvoxy + scm_input%input_alboldxy = input_alboldxy + scm_input%input_qsnowxy = input_qsnowxy + scm_input%input_wslakexy = input_wslakexy + scm_input%input_taussxy = input_taussxy + scm_input%input_waxy = input_waxy + scm_input%input_wtxy = input_wtxy + scm_input%input_zwtxy = input_zwtxy + scm_input%input_xlaixy = input_xlaixy + scm_input%input_xsaixy = input_xsaixy + scm_input%input_lfmassxy = input_lfmassxy + scm_input%input_stmassxy = input_stmassxy + scm_input%input_rtmassxy = input_rtmassxy + scm_input%input_woodxy = input_woodxy + scm_input%input_stblcpxy = input_stblcpxy + scm_input%input_fastcpxy = input_fastcpxy + scm_input%input_smcwtdxy = input_smcwtdxy + scm_input%input_deeprechxy = input_deeprechxy + scm_input%input_rechxy = input_rechxy + scm_input%input_snowxy = input_snowxy + + scm_input%input_tref = input_tref + scm_input%input_z_c = input_z_c + scm_input%input_c_0 = input_c_0 + scm_input%input_c_d = input_c_d + scm_input%input_w_0 = input_w_0 + scm_input%input_w_d = input_w_d + scm_input%input_xt = input_xt + scm_input%input_xs = input_xs + scm_input%input_xu = input_xu + scm_input%input_xv = input_xv + scm_input%input_xz = input_xz + scm_input%input_zm = input_zm + scm_input%input_xtts = input_xtts + scm_input%input_xzts = input_xzts + scm_input%input_d_conv = input_d_conv + scm_input%input_ifd = input_ifd + scm_input%input_dt_cool = input_dt_cool + scm_input%input_qrain = input_qrain + + scm_input%input_wetness = input_wetness + scm_input%input_clw_surf = input_clw_surf + scm_input%input_qwv_surf = input_qwv_surf + scm_input%input_tsnow = input_tsnow + scm_input%input_snowfallac = input_snowfallac + scm_input%input_acsnow = input_acsnow + scm_input%input_lai = input_lai + !> @} end subroutine get_case_init !> Subroutine to get reference profile to use above the case data (temporarily hard-coded profile) subroutine get_reference_profile(scm_state, scm_reference) use gmtb_scm_type_defs, only : scm_state_type, scm_reference_type - + use NetCDF_read, only: check + type(scm_state_type), target, intent(in) :: scm_state type(scm_reference_type), target, intent(inout) :: scm_reference @@ -904,15 +967,6 @@ subroutine get_tracers(tracer_names) close (fu) end subroutine get_tracers -!> Generic subroutine to check for netCDF I/O errors -subroutine check(status) - integer, intent ( in) :: status - - if(status /= nf90_noerr) then - print *, trim(nf90_strerror(status)) - stop "stopped" - end if -end subroutine check !> @} !> @} end module gmtb_scm_input diff --git a/scm/src/gmtb_scm_output.F90 b/scm/src/gmtb_scm_output.F90 index e2bee7a8..616adb8c 100644 --- a/scm/src/gmtb_scm_output.F90 +++ b/scm/src/gmtb_scm_output.F90 @@ -4,7 +4,7 @@ module gmtb_scm_output use netcdf -use gmtb_scm_input, only: check +use NetCDF_read, only: check use gmtb_scm_kinds, only: sp, dp, qp implicit none diff --git a/scm/src/gmtb_scm_physical_constants.meta b/scm/src/gmtb_scm_physical_constants.meta index 04fa0b7c..6bca2ab3 100644 --- a/scm/src/gmtb_scm_physical_constants.meta +++ b/scm/src/gmtb_scm_physical_constants.meta @@ -202,3 +202,10 @@ dimensions = () type = real kind = kind_phys +[con_omega] + standard_name = angular_velocity_of_earth + long_name = angular velocity of earth + units = s-1 + dimensions = () + type = real + kind = kind_phys diff --git a/scm/src/gmtb_scm_setup.F90 b/scm/src/gmtb_scm_setup.F90 index 605667f4..f21bfca1 100644 --- a/scm/src/gmtb_scm_setup.F90 +++ b/scm/src/gmtb_scm_setup.F90 @@ -125,92 +125,14 @@ subroutine set_state(scm_input, scm_reference, scm_state) scm_state%state_T(i,:,1) = scm_input%input_temp(:) scm_state%state_tracer(i,:,scm_state%water_vapor_index,1)=scm_input%input_qt scm_state%state_tracer(i,:,scm_state%ozone_index,1)=scm_input%input_ozone - scm_state%veg_type(i) = scm_input%input_vegtyp - scm_state%soil_type(i) = scm_input%input_soiltyp - scm_state%slope_type(i) = scm_input%input_slopetype - scm_state%veg_frac(i) = scm_input%input_vegfrac - scm_state%shdmin(i) = scm_input%input_shdmin - scm_state%shdmax(i) = scm_input%input_shdmax - scm_state%sfc_roughness_length_cm = scm_input%input_zorl - scm_state%sfc_type(i) = scm_input%input_slmsk !< this "overwrites" what is in the SCM case namelist if model ICs are present - scm_state%canopy(i) = scm_input%input_canopy - scm_state%hice(i) = scm_input%input_hice - scm_state%fice(i) = scm_input%input_fice - scm_state%tisfc(i) = scm_input%input_tisfc - scm_state%snwdph(i) = scm_input%input_snwdph - scm_state%snoalb(i) = scm_input%input_snoalb - scm_state%sncovr(i) = scm_input%input_sncovr scm_state%area(i) = scm_input%input_area - scm_state%tg3(i) = scm_input%input_tg3 - scm_state%uustar(i) = scm_input%input_uustar - scm_state%stc(i,:,1)=scm_input%input_stc - scm_state%smc(i,:,1)=scm_input%input_smc - scm_state%slc(i,:,1)=scm_input%input_slc + if (scm_input%input_pres_i(1).GT. 0.0) then ! pressure are read in, overwrite values scm_state%pres_i(i,:)=scm_input%input_pres_i scm_state%pres_l(i,:)=scm_input%input_pres_l endif - scm_state%alvsf(i)=scm_input%input_alvsf - scm_state%alnsf(i)=scm_input%input_alnsf - scm_state%alvwf(i)=scm_input%input_alvwf - scm_state%alnwf(i)=scm_input%input_alnwf - scm_state%hprime(i,1)=scm_input%input_stddev - scm_state%hprime(i,2)=scm_input%input_convexity - scm_state%hprime(i,3)=scm_input%input_oa1 - scm_state%hprime(i,4)=scm_input%input_oa2 - scm_state%hprime(i,5)=scm_input%input_oa3 - scm_state%hprime(i,6)=scm_input%input_oa4 - scm_state%hprime(i,7)=scm_input%input_ol1 - scm_state%hprime(i,8)=scm_input%input_ol2 - scm_state%hprime(i,9)=scm_input%input_ol3 - scm_state%hprime(i,10)=scm_input%input_ol4 - scm_state%hprime(i,11)=scm_input%input_theta - scm_state%hprime(i,12)=scm_input%input_gamma - scm_state%hprime(i,13)=scm_input%input_sigma - scm_state%hprime(i,14)=scm_input%input_elvmax - scm_state%facsf(i)=scm_input%input_facsf - scm_state%facwf(i)=scm_input%input_facwf enddo - !check for nonzero NoahMP input variable and fill in the scm_state with values from scm_input if found - if (scm_input%input_tvxy /= 0.0) then - do i=1, scm_state%n_cols - scm_state%tvxy(i) = scm_input%input_tvxy - scm_state%tgxy(i) = scm_input%input_tgxy - scm_state%tahxy(i) = scm_input%input_tahxy - scm_state%canicexy(i) = scm_input%input_canicexy - scm_state%canliqxy(i) = scm_input%input_canliqxy - scm_state%eahxy(i) = scm_input%input_eahxy - scm_state%cmxy(i) = scm_input%input_cmxy - scm_state%chxy(i) = scm_input%input_chxy - scm_state%fwetxy(i) = scm_input%input_fwetxy - scm_state%sneqvoxy(i) = scm_input%input_sneqvoxy - scm_state%alboldxy(i) = scm_input%input_alboldxy - scm_state%qsnowxy(i) = scm_input%input_qsnowxy - scm_state%wslakexy(i) = scm_input%input_wslakexy - scm_state%taussxy(i) = scm_input%input_taussxy - scm_state%waxy(i) = scm_input%input_waxy - scm_state%wtxy(i) = scm_input%input_wtxy - scm_state%zwtxy(i) = scm_input%input_zwtxy - scm_state%xlaixy(i) = scm_input%input_xlaixy - scm_state%xsaixy(i) = scm_input%input_xsaixy - scm_state%lfmassxy(i) = scm_input%input_lfmassxy - scm_state%stmassxy(i) = scm_input%input_stmassxy - scm_state%rtmassxy(i) = scm_input%input_rtmassxy - scm_state%woodxy(i) = scm_input%input_woodxy - scm_state%stblcpxy(i) = scm_input%input_stblcpxy - scm_state%fastcpxy(i) = scm_input%input_fastcpxy - scm_state%smcwtdxy(i) = scm_input%input_smcwtdxy - scm_state%deeprechxy(i) = scm_input%input_deeprechxy - scm_state%rechxy(i) = scm_input%input_rechxy - scm_state%snowxy(i) = scm_input%input_snowxy - - scm_state%snicexy(i,:) = scm_input%input_snicexy(:) - scm_state%snliqxy(i,:) = scm_input%input_snliqxy(:) - scm_state%tsnoxy(i,:) = scm_input%input_tsnoxy(:) - scm_state%smoiseq(i,:) = scm_input%input_smoiseq(:) - scm_state%zsnsoxy(i,:) = scm_input%input_zsnsoxy(:) - end do - endif + endif !> @} end subroutine set_state diff --git a/scm/src/gmtb_scm_type_defs.F90 b/scm/src/gmtb_scm_type_defs.F90 index f794de33..8a83ea24 100644 --- a/scm/src/gmtb_scm_type_defs.F90 +++ b/scm/src/gmtb_scm_type_defs.F90 @@ -122,66 +122,9 @@ module gmtb_scm_type_defs real(kind=dp), allocatable :: T_surf(:), pres_surf(:) !< surface temperature and pressure interpolated to the model time real(kind=dp), allocatable :: sh_flux(:), lh_flux(:) !< surface sensible and latent heat fluxes interpolated to the model time real(kind=dp), allocatable :: sfc_roughness_length_cm(:) !< surface roughness length used for calculating surface layer parameters from specified fluxes - real(kind=dp), allocatable :: alvsf(:), alnsf(:),alvwf(:),alnwf(:) !< surface albedos - real(kind=dp), allocatable :: facsf(:), facwf(:), hprime(:,:) !< other surface stuff - real(kind=dp), allocatable :: stc(:,:,:) !< soil temperature - real(kind=dp), allocatable :: smc(:,:,:) !< soil moisture - real(kind=dp), allocatable :: slc(:,:,:) !< soil liquid content real(kind=dp), allocatable :: sfc_type(:) !< 0: sea surface, 1: land surface, 2: sea-ice surface - real(kind=dp), allocatable :: veg_type(:) !< vegetation type classification - real(kind=dp), allocatable :: slope_type(:) !< surface slope classification - real(kind=dp), allocatable :: soil_type(:) !< soil type classification - real(kind=dp), allocatable :: veg_frac(:) !< vegetation area fraction - real(kind=dp), allocatable :: shdmin(:) !< minimun vegetation fraction - real(kind=dp), allocatable :: shdmax(:) !< maximun vegetation fraction - real(kind=dp), allocatable :: tg3(:) !< deep soil temperature (K) - real(kind=dp), allocatable :: slmsk(:) !< sea land ice mask [0,1,2] - real(kind=dp), allocatable :: canopy(:) !< amount of water stored in canopy (kg m-2) - real(kind=dp), allocatable :: hice(:) !< sea ice thickness (m) - real(kind=dp), allocatable :: fice(:) !< ice fraction (frac) - real(kind=dp), allocatable :: tisfc(:) !< ice surface temperature (K) - real(kind=dp), allocatable :: snwdph(:) !< water equivalent snow depth (mm) - real(kind=dp), allocatable :: snoalb(:) !< maximum snow albedo (frac) - real(kind=dp), allocatable :: sncovr(:) !< snow area fraction (frac) - real(kind=dp), allocatable :: uustar(:) !< surface friction velocity (m s-1) - - real(kind=dp), allocatable :: tvxy(:) !< vegetation temperature (K) - real(kind=dp), allocatable :: tgxy(:) !< ground temperature for Noahmp (K) - real(kind=dp), allocatable :: tahxy(:) !< canopy air temperature (K) - real(kind=dp), allocatable :: canicexy(:) !< canopy intercepted ice mass (mm) - real(kind=dp), allocatable :: canliqxy(:) !< canopy intercepted liquid water (mm) - real(kind=dp), allocatable :: eahxy(:) !< canopy air vapor pressure (Pa) - real(kind=dp), allocatable :: cmxy(:) !< surface drag coefficient for momentum for noahmp - real(kind=dp), allocatable :: chxy(:) !< surface exchange coeff heat & moisture for noahmp - real(kind=dp), allocatable :: fwetxy(:) !< area fraction of canopy that is wetted/snowed - real(kind=dp), allocatable :: sneqvoxy(:) !< snow mass at previous time step (mm) - real(kind=dp), allocatable :: alboldxy(:) !< snow albedo at previous time step (frac) - real(kind=dp), allocatable :: qsnowxy(:) !< snow precipitation rate at surface (mm s-1) - real(kind=dp), allocatable :: wslakexy(:) !< lake water storage (mm) - real(kind=dp), allocatable :: taussxy(:) !< non-dimensional snow age - real(kind=dp), allocatable :: waxy(:) !< water storage in aquifer (mm) - real(kind=dp), allocatable :: wtxy(:) !< water storage in aquifer and saturated soil (mm) - real(kind=dp), allocatable :: zwtxy(:) !< water table depth (m) - real(kind=dp), allocatable :: xlaixy(:) !< leaf area index - real(kind=dp), allocatable :: xsaixy(:) !< stem area index - real(kind=dp), allocatable :: lfmassxy(:) !< leaf mass (g m-2) - real(kind=dp), allocatable :: stmassxy(:) !< stem mass (g m-2) - real(kind=dp), allocatable :: rtmassxy(:) !< fine root mass (g m-2) - real(kind=dp), allocatable :: woodxy(:) !< wood mass including woody roots (g m-2) - real(kind=dp), allocatable :: stblcpxy(:) !< stable carbon in deep soil (g m-2) - real(kind=dp), allocatable :: fastcpxy(:) !< short-lived carbon in shallow soil (g m-2) - real(kind=dp), allocatable :: smcwtdxy(:) !< soil water content between the bottom of the soil and the water table (m3 m-3) - real(kind=dp), allocatable :: deeprechxy(:) !< recharge to or from the water table when deep (m) - real(kind=dp), allocatable :: rechxy(:) !< recharge to or from the water table when shallow (m) - real(kind=dp), allocatable :: snowxy(:) !< number of snow layers - - real(kind=dp), allocatable :: snicexy(:,:) !< snow layer ice (mm) - real(kind=dp), allocatable :: snliqxy(:,:) !< snow layer liquid (mm) - real(kind=dp), allocatable :: tsnoxy(:,:) !< snow temperature (K) - real(kind=dp), allocatable :: smoiseq(:,:) !< equilibrium soil water content (m3 m-3) - real(kind=dp), allocatable :: zsnsoxy(:,:) !< layer bottom depth from snow surface (m) - + contains procedure :: create => scm_state_create @@ -192,77 +135,126 @@ module gmtb_scm_type_defs integer :: input_nlev !< number of levels in the input file integer :: input_nsoil !< number of soil levels in the input file integer :: input_nsnow !< number of snow layers in the input file + integer :: input_nice !< number of sea ice layers in the input file integer :: input_ntimes !< number of times in the input file where forcing is available - integer :: input_vegsrc !< - real(kind=dp) :: input_vegtyp !< + integer :: input_vegsrc !< vegetation source + real(kind=dp) :: input_vegtyp !< vegetation type classification real(kind=dp) :: input_soiltyp !< - real(kind=dp) :: input_slopetype !< + real(kind=dp) :: input_slopetype !< surface slope classification real(kind=dp) :: input_lat !< latitude of column center real(kind=dp) :: input_lon !< longitude of column center - real(kind=dp) :: input_vegfrac !< - real(kind=dp) :: input_shdmin !< - real(kind=dp) :: input_shdmax !< - real(kind=dp) :: input_zorl !< + real(kind=dp) :: input_tsfco !< input sea surface temperature OR surface skin temperature over land OR surface skin temperature over ice (depending on slmsk) (K) + real(kind=dp) :: input_vegfrac !< vegetation area fraction + real(kind=dp) :: input_shdmin !< minimun vegetation fraction + real(kind=dp) :: input_shdmax !< maximun vegetation fraction + real(kind=dp) :: input_zorlo !< surfce roughness length over ocean [cm] real(kind=dp) :: input_slmsk !< sea land ice mask [0,1,2] - real(kind=dp) :: input_canopy !< amount of water stored in canopy - real(kind=dp) :: input_hice !< ice thickness - real(kind=dp) :: input_fice !< ice fraction - real(kind=dp) :: input_tisfc !< ice surface temperature - real(kind=dp) :: input_snwdph !< snow depth - real(kind=dp) :: input_snoalb !< snow albedo - real(kind=dp) :: input_sncovr !< snow cover - real(kind=dp) :: input_area !< - real(kind=dp) :: input_tg3 !< - real(kind=dp) :: input_uustar !< - real(kind=dp) :: input_alvsf !< - real(kind=dp) :: input_alnsf !< - real(kind=dp) :: input_alvwf !< - real(kind=dp) :: input_alnwf !< - real(kind=dp) :: input_convexity !< - real(kind=dp) :: input_stddev !< - real(kind=dp) :: input_oa1 !< - real(kind=dp) :: input_oa2 !< - real(kind=dp) :: input_oa3 !< - real(kind=dp) :: input_oa4 !< - real(kind=dp) :: input_ol1 !< - real(kind=dp) :: input_ol2 !< - real(kind=dp) :: input_ol3 !< - real(kind=dp) :: input_ol4 !< - real(kind=dp) :: input_theta !< - real(kind=dp) :: input_gamma !< - real(kind=dp) :: input_sigma !< - real(kind=dp) :: input_elvmax !< - real(kind=dp) :: input_facsf !< fraction of surface depedent on sun angle - real(kind=dp) :: input_facwf !< fraction of surface not depedent on sun angle - real(kind=dp) :: input_tvxy !< - real(kind=dp) :: input_tgxy !< - real(kind=dp) :: input_tahxy !< - real(kind=dp) :: input_canicexy !< - real(kind=dp) :: input_canliqxy !< - real(kind=dp) :: input_eahxy !< - real(kind=dp) :: input_cmxy !< - real(kind=dp) :: input_chxy !< - real(kind=dp) :: input_fwetxy !< - real(kind=dp) :: input_sneqvoxy !< - real(kind=dp) :: input_alboldxy !< - real(kind=dp) :: input_qsnowxy !< - real(kind=dp) :: input_wslakexy !< - real(kind=dp) :: input_taussxy !< - real(kind=dp) :: input_waxy !< - real(kind=dp) :: input_wtxy !< - real(kind=dp) :: input_zwtxy !< - real(kind=dp) :: input_xlaixy !< - real(kind=dp) :: input_xsaixy !< - real(kind=dp) :: input_lfmassxy !< - real(kind=dp) :: input_stmassxy !< - real(kind=dp) :: input_rtmassxy !< - real(kind=dp) :: input_woodxy !< - real(kind=dp) :: input_stblcpxy !< - real(kind=dp) :: input_fastcpxy !< - real(kind=dp) :: input_smcwtdxy !< - real(kind=dp) :: input_deeprechxy !< - real(kind=dp) :: input_rechxy !< - real(kind=dp) :: input_snowxy !< + real(kind=dp) :: input_canopy !< amount of water stored in canopy (kg m-2) + real(kind=dp) :: input_hice !< sea ice thickness (m) + real(kind=dp) :: input_fice !< ice fraction (frac) + real(kind=dp) :: input_tisfc !< ice surface temperature (K) + real(kind=dp) :: input_snwdph !< water equivalent snow depth (mm) + real(kind=dp) :: input_snoalb !< maximum snow albedo (frac) + real(kind=dp) :: input_sncovr !< snow area fraction (frac) + real(kind=dp) :: input_area !< surface area [m^2] + real(kind=dp) :: input_tg3 !< deep soil temperature (K) + real(kind=dp) :: input_uustar !< surface friction velocity (m s-1) + real(kind=dp) :: input_alvsf !< 60 degree vis albedo with strong cosz dependency + real(kind=dp) :: input_alnsf !< 60 degree nir albedo with strong cosz dependency + real(kind=dp) :: input_alvwf !< 60 degree vis albedo with weak cosz dependency + real(kind=dp) :: input_alnwf !< 60 degree nir albedo with weak cosz dependency + real(kind=dp) :: input_facsf !< fractional coverage with strong cosz dependency + real(kind=dp) :: input_facwf !< fractional coverage with weak cosz dependency + real(kind=dp) :: input_weasd !< water equivalent accumulated snow depth (mm) + real(kind=dp) :: input_f10m !< ratio of sigma level 1 wind and 10m wind + real(kind=dp) :: input_t2m !< 2-meter absolute temperature (K) + real(kind=dp) :: input_q2m !< 2-meter specific humidity (kg kg-1) + real(kind=dp) :: input_ffmm !< Monin-Obukhov similarity function for momentum + real(kind=dp) :: input_ffhh !< Monin-Obukhov similarity function for heat + real(kind=dp) :: input_tprcp !< instantaneous total precipitation amount (m) + real(kind=dp) :: input_srflag !< snow/rain flag for precipitation + real(kind=dp) :: input_tsfcl !< surface skin temperature over land (K) + real(kind=dp) :: input_zorll !< surface roughness length over land (cm) + real(kind=dp) :: input_zorli !< surface roughness length over ice (cm) + real(kind=dp) :: input_zorlw !< surface roughness length from wave model (cm) + + real(kind=dp) :: input_stddev !< standard deviation of subgrid orography (m) + real(kind=dp) :: input_convexity !< convexity of subgrid orography + real(kind=dp) :: input_ol1 !< fraction of grid box with subgrid orography higher than critical height 1 + real(kind=dp) :: input_ol2 !< fraction of grid box with subgrid orography higher than critical height 2 + real(kind=dp) :: input_ol3 !< fraction of grid box with subgrid orography higher than critical height 3 + real(kind=dp) :: input_ol4 !< fraction of grid box with subgrid orography higher than critical height 4 + real(kind=dp) :: input_oa1 !< assymetry of subgrid orography 1 + real(kind=dp) :: input_oa2 !< assymetry of subgrid orography 2 + real(kind=dp) :: input_oa3 !< assymetry of subgrid orography 3 + real(kind=dp) :: input_oa4 !< assymetry of subgrid orography 4 + real(kind=dp) :: input_sigma !< slope of subgrid orography + real(kind=dp) :: input_theta !< angle with respect to east of maximum subgrid orographic variations (deg) + real(kind=dp) :: input_gamma !< anisotropy of subgrid orography + real(kind=dp) :: input_elvmax!< maximum of subgrid orography (m) + real(kind=dp) :: input_oro !< orography (m) + real(kind=dp) :: input_oro_uf !< unfiltered orography (m) + real(kind=dp) :: input_landfrac !< fraction of horizontal grid area occupied by land + real(kind=dp) :: input_lakefrac !< fraction of horizontal grid area occupied by lake + real(kind=dp) :: input_lakedepth !< lake depth (m) + + real(kind=dp) :: input_tvxy !< vegetation temperature (K) + real(kind=dp) :: input_tgxy !< ground temperature for Noahmp (K) + real(kind=dp) :: input_tahxy !< canopy air temperature (K) + real(kind=dp) :: input_canicexy !< canopy intercepted ice mass (mm) + real(kind=dp) :: input_canliqxy !< canopy intercepted liquid water (mm) + real(kind=dp) :: input_eahxy !< canopy air vapor pressure (Pa) + real(kind=dp) :: input_cmxy !< surface drag coefficient for momentum for noahmp + real(kind=dp) :: input_chxy !< surface exchange coeff heat & moisture for noahmp + real(kind=dp) :: input_fwetxy !< area fraction of canopy that is wetted/snowed + real(kind=dp) :: input_sneqvoxy !< snow mass at previous time step (mm) + real(kind=dp) :: input_alboldxy !< snow albedo at previous time step (frac) + real(kind=dp) :: input_qsnowxy !< snow precipitation rate at surface (mm s-1) + real(kind=dp) :: input_wslakexy !< lake water storage (mm) + real(kind=dp) :: input_taussxy !< non-dimensional snow age + real(kind=dp) :: input_waxy !< water storage in aquifer (mm) + real(kind=dp) :: input_wtxy !< water storage in aquifer and saturated soil (mm) + real(kind=dp) :: input_zwtxy !< water table depth (m) + real(kind=dp) :: input_xlaixy !< leaf area index + real(kind=dp) :: input_xsaixy !< stem area index + real(kind=dp) :: input_lfmassxy !< leaf mass (g m-2) + real(kind=dp) :: input_stmassxy !< stem mass (g m-2) + real(kind=dp) :: input_rtmassxy !< fine root mass (g m-2) + real(kind=dp) :: input_woodxy !< wood mass including woody roots (g m-2) + real(kind=dp) :: input_stblcpxy !< stable carbon in deep soil (g m-2) + real(kind=dp) :: input_fastcpxy !< short-lived carbon in shallow soil (g m-2) + real(kind=dp) :: input_smcwtdxy !< soil water content between the bottom of the soil and the water table (m3 m-3) + real(kind=dp) :: input_deeprechxy !< recharge to or from the water table when deep (m) + real(kind=dp) :: input_rechxy !< recharge to or from the water table when shallow (m) + real(kind=dp) :: input_snowxy !< number of snow layers + + real(kind=dp) :: input_tref !< sea surface reference temperature for NSST (K) + real(kind=dp) :: input_z_c !< sub-layer cooling thickness for NSST (m) + real(kind=dp) :: input_c_0 !< coefficient 1 to calculate d(Tz)/d(Ts) for NSST + real(kind=dp) :: input_c_d !< coefficient 2 to calculate d(Tz)/d(Ts) for NSST + real(kind=dp) :: input_w_0 !< coefficient 3 to calculate d(Tz)/d(Ts) for NSST + real(kind=dp) :: input_w_d !< coefficient 4 to calculate d(Tz)/d(Ts) for NSST + real(kind=dp) :: input_xt !< heat content in diurnal thermocline layer for NSST (K m) + real(kind=dp) :: input_xs !< salinity content in diurnal thermocline layer for NSST (ppt m) + real(kind=dp) :: input_xu !< u-current in diurnal thermocline layer for NSST (m2 s-1) + real(kind=dp) :: input_xv !< v-current in diurnal thermocline layer for NSST (m2 s-1) + real(kind=dp) :: input_xz !< thickness of diurnal thermocline layer for NSST (m) + real(kind=dp) :: input_zm !< thickness of ocean mixed layer for NSST (m) + real(kind=dp) :: input_xtts !< sensitivity of diurnal thermocline layer heat content to surface temperature [d(xt)/d(ts)] for NSST (m) + real(kind=dp) :: input_xzts !< sensitivity of diurnal thermocline layer thickness to surface temperature [d(xz)/d(ts)] for NSST (m K-1) + real(kind=dp) :: input_d_conv !< thickness of free convection layer for NSST (m) + real(kind=dp) :: input_ifd !< index to start DTM run for NSST + real(kind=dp) :: input_dt_cool !< sub-layer cooling amount for NSST (K) + real(kind=dp) :: input_qrain !< sensible heat due to rainfall for NSST (W) + + real(kind=dp) :: input_wetness !< normalized soil wetness for RUC LSM + real(kind=dp) :: input_clw_surf !< cloud condensed water mixing ratio at surface for RUC LSM (kg kg-1) + real(kind=dp) :: input_qwv_surf !< water vapor mixing ratio at surface for RUC LSM (kg kg-1) + real(kind=dp) :: input_tsnow !< snow temperature at the bottom of the first snow layer for RUC LSM (K) + real(kind=dp) :: input_snowfallac !< run-total snow accumulation on the ground for RUC LSM (kg m-2) + real(kind=dp) :: input_acsnow !< snow water equivalent of run-total frozen precip for RUC LSM (kg m-2) + real(kind=dp) :: input_lai !< leaf area index for RUC LSM + real(kind=dp), allocatable :: input_pres_i(:) !< pressure (Pa) of input interface real(kind=dp), allocatable :: input_pres_l(:) !< pressure (Pa) of input levels real(kind=dp), allocatable :: input_pres(:) !< pressure (Pa) of input levels @@ -280,11 +272,17 @@ module gmtb_scm_type_defs real(kind=dp), allocatable :: input_stc(:) !< soil temperature (k) (initial) real(kind=dp), allocatable :: input_smc(:) !< soil moisture conteng (g/g) (initial) real(kind=dp), allocatable :: input_slc(:) !< soil liquid content (g/g) (initial) - real(kind=dp), allocatable :: input_snicexy(:) !< - real(kind=dp), allocatable :: input_snliqxy(:) !< - real(kind=dp), allocatable :: input_tsnoxy(:) !< - real(kind=dp), allocatable :: input_smoiseq(:) !< - real(kind=dp), allocatable :: input_zsnsoxy(:) !< + real(kind=dp), allocatable :: input_snicexy(:) !< snow layer ice (mm) + real(kind=dp), allocatable :: input_snliqxy(:) !< snow layer liquid (mm) + real(kind=dp), allocatable :: input_tsnoxy(:) !< snow temperature (K) + real(kind=dp), allocatable :: input_smoiseq(:) !< equilibrium soil water content (m3 m-3) + real(kind=dp), allocatable :: input_zsnsoxy(:) !< layer bottom depth from snow surface (m) + real(kind=dp), allocatable :: input_tiice(:) !< sea ice internal temperature (K) + real(kind=dp), allocatable :: input_tslb(:) !< soil temperature for RUC LSM (K) + real(kind=dp), allocatable :: input_smois(:) !< volume fraction of soil moisture for RUC LSM (frac) + real(kind=dp), allocatable :: input_sh2o(:) !< volume fraction of unfrozen soil moisture for RUC LSM (frac) + real(kind=dp), allocatable :: input_smfr(:) !< volume fraction of frozen soil moisture for RUC LSM (frac) + real(kind=dp), allocatable :: input_flfr(:) !< flag for frozen soil physics real(kind=dp), allocatable :: input_pres_surf(:) !< time-series of surface pressure (Pa) real(kind=dp), allocatable :: input_T_surf(:) !< time-series of surface temperture real(kind=dp), allocatable :: input_w_ls(:,:) !< large-scale vertical velocity (m/s) (time, levels) @@ -306,7 +304,6 @@ module gmtb_scm_type_defs contains procedure :: create => scm_input_create - procedure :: create_modelics => scm_input_create_model_ics end type scm_input_type @@ -344,6 +341,7 @@ module gmtb_scm_type_defs contains procedure :: create => physics_create procedure :: associate => physics_associate + procedure :: set => physics_set end type physics_type type(physics_type), target :: physics @@ -459,9 +457,6 @@ subroutine scm_state_create(scm_state, n_columns, n_levels, n_soil, n_snow, n_ti allocate(scm_state%temp_tracer(n_columns, n_levels, scm_state%n_tracers, n_time_levels), & scm_state%temp_T(n_columns, n_levels, n_time_levels), & scm_state%temp_u(n_columns, n_levels, n_time_levels), scm_state%temp_v(n_columns, n_levels, n_time_levels)) - allocate(scm_state%stc(n_columns, n_soil, n_time_levels)) - allocate(scm_state%smc(n_columns, n_soil, n_time_levels)) - allocate(scm_state%slc(n_columns, n_soil, n_time_levels)) scm_state%temp_tracer = real_zero scm_state%temp_T = real_zero scm_state%temp_u = real_zero @@ -501,123 +496,21 @@ subroutine scm_state_create(scm_state, n_columns, n_levels, n_soil, n_snow, n_ti scm_state%qv_force_tend = real_zero scm_state%sfc_roughness_length_cm = real_one - allocate(scm_state%alvsf(n_columns), scm_state%alnsf(n_columns), scm_state%alvwf(n_columns), scm_state%alnwf(n_columns), & - scm_state%facsf(n_columns), scm_state%facwf(n_columns), scm_state%hprime(n_columns,14), & - scm_state%sfc_type(n_columns), scm_state%veg_type(n_columns), & - scm_state%veg_frac(n_columns), scm_state%slope_type(n_columns), scm_state%shdmin(n_columns), scm_state%shdmax(n_columns), & - scm_state%tg3(n_columns), scm_state%slmsk(n_columns), scm_state%canopy(n_columns), scm_state%hice(n_columns), scm_state%fice(n_columns), & - scm_state%tisfc(n_columns), scm_state%snwdph(n_columns), scm_state%snoalb(n_columns), scm_state%sncovr(n_columns), scm_state%uustar(n_columns), & - scm_state%soil_type(n_columns)) - - scm_state%alvsf = real_zero - scm_state%alnsf = real_zero - scm_state%alvwf = real_zero - scm_state%alnwf = real_zero - scm_state%facsf = real_zero - scm_state%facwf = real_zero - scm_state%hprime = real_zero + allocate(scm_state%sfc_type(n_columns)) + scm_state%sfc_type = real_zero - scm_state%veg_type = real_zero - scm_state%veg_frac = real_zero - scm_state%slope_type = real_zero - scm_state%shdmin = real_zero - scm_state%shdmax = real_zero - scm_state%tg3 = real_zero - scm_state%slmsk = real_zero - scm_state%canopy = real_zero - scm_state%hice = real_zero - scm_state%fice = real_zero - scm_state%tisfc = real_zero - scm_state%snwdph = real_zero - scm_state%snoalb = real_zero - scm_state%sncovr = real_zero - scm_state%uustar = real_zero - scm_state%soil_type = real_zero - - allocate(scm_state%tvxy(n_columns), scm_state%tgxy(n_columns), scm_state%tahxy(n_columns), scm_state%canicexy(n_columns), & - scm_state%canliqxy(n_columns), scm_state%eahxy(n_columns), scm_state%cmxy(n_columns), scm_state%chxy(n_columns), & - scm_state%fwetxy(n_columns), scm_state%sneqvoxy(n_columns), scm_state%alboldxy(n_columns), scm_state%qsnowxy(n_columns), & - scm_state%wslakexy(n_columns), scm_state%taussxy(n_columns), scm_state%waxy(n_columns), scm_state%wtxy(n_columns), & - scm_state%zwtxy(n_columns), scm_state%xlaixy(n_columns), scm_state%xsaixy(n_columns), scm_state%lfmassxy(n_columns), & - scm_state%stmassxy(n_columns), scm_state%rtmassxy(n_columns), scm_state%woodxy(n_columns), scm_state%stblcpxy(n_columns), & - scm_state%fastcpxy(n_columns), scm_state%smcwtdxy(n_columns), scm_state%deeprechxy(n_columns), scm_state%rechxy(n_columns), & - scm_state%snowxy(n_columns)) - - scm_state%tvxy = real_zero - scm_state%tgxy = real_zero - scm_state%tahxy = real_zero - scm_state%canicexy = real_zero - scm_state%canliqxy = real_zero - scm_state%eahxy = real_zero - scm_state%cmxy = real_zero - scm_state%chxy = real_zero - scm_state%fwetxy = real_zero - scm_state%sneqvoxy = real_zero - scm_state%alboldxy = real_zero - scm_state%qsnowxy = real_zero - scm_state%wslakexy = real_zero - scm_state%taussxy = real_zero - scm_state%waxy = real_zero - scm_state%wtxy = real_zero - scm_state%zwtxy = real_zero - scm_state%xlaixy = real_zero - scm_state%xsaixy = real_zero - scm_state%lfmassxy = real_zero - scm_state%stmassxy = real_zero - scm_state%rtmassxy = real_zero - scm_state%woodxy = real_zero - scm_state%stblcpxy = real_zero - scm_state%fastcpxy = real_zero - scm_state%smcwtdxy = real_zero - scm_state%deeprechxy = real_zero - scm_state%rechxy = real_zero - scm_state%snowxy = real_zero - - allocate(scm_state%snicexy(n_columns,n_snow), scm_state%snliqxy(n_columns,n_snow), scm_state%tsnoxy(n_columns,n_snow), & - scm_state%smoiseq(n_columns,n_soil), scm_state%zsnsoxy(n_columns,n_snow + n_soil)) - - scm_state%snicexy = real_zero - scm_state%snliqxy = real_zero - scm_state%tsnoxy = real_zero - scm_state%smoiseq = real_zero - scm_state%zsnsoxy = real_zero end subroutine scm_state_create - subroutine scm_input_create_model_ics(scm_input,nlev_soil,nlev_snow,nlev,noahmp) - class(scm_input_type) :: scm_input - integer, intent(in) :: nlev,nlev_soil,nlev_snow - logical, intent(in) :: noahmp - - scm_input%input_nsoil= nlev_soil - allocate(scm_input%input_stc(nlev_soil), scm_input%input_smc(nlev_soil), scm_input%input_slc(nlev_soil)) - allocate(scm_input%input_temp(nlev), scm_input%input_pres_i(nlev+1),scm_input%input_pres_l(nlev)) - scm_input%input_stc = real_zero - scm_input%input_smc = real_zero - scm_input%input_slc = real_zero - scm_input%input_temp = real_zero - scm_input%input_pres_i = real_zero - scm_input%input_pres_l = real_zero - - if (noahmp) then - scm_input%input_nsnow = nlev_snow - allocate(scm_input%input_snicexy(nlev_snow), scm_input%input_snliqxy(nlev_snow), scm_input%input_tsnoxy(nlev_snow), & - scm_input%input_smoiseq(nlev_soil), scm_input%input_zsnsoxy(nlev_snow + nlev_soil)) - scm_input%input_snicexy = real_zero - scm_input%input_snliqxy = real_zero - scm_input%input_tsnoxy = real_zero - scm_input%input_smoiseq = real_zero - scm_input%input_zsnsoxy = real_zero - endif - - end subroutine scm_input_create_model_ics - - subroutine scm_input_create(scm_input, ntimes, nlev) + subroutine scm_input_create(scm_input, ntimes, nlev, nsoil, nsnow, nice) class(scm_input_type) :: scm_input - integer, intent(in) :: ntimes, nlev + integer, intent(in) :: ntimes, nlev, nsoil, nsnow, nice scm_input%input_nlev = nlev scm_input%input_ntimes = ntimes + scm_input%input_nsoil = nsoil + scm_input%input_nsnow = nsnow + scm_input%input_nice = nice allocate(scm_input%input_pres(nlev),scm_input%input_time(ntimes)) scm_input%input_pres = real_zero @@ -625,9 +518,10 @@ subroutine scm_input_create(scm_input, ntimes, nlev) allocate(scm_input%input_height(nlev), scm_input%input_thetail(nlev), scm_input%input_qt(nlev), scm_input%input_ql(nlev), & scm_input%input_qi(nlev), scm_input%input_u(nlev), scm_input%input_v(nlev), scm_input%input_tke(nlev), & - scm_input%input_ozone(nlev)) + scm_input%input_ozone(nlev), scm_input%input_temp(nlev)) scm_input%input_height = real_zero scm_input%input_thetail = real_zero + scm_input%input_temp = real_zero scm_input%input_qt = real_zero scm_input%input_ql = real_zero scm_input%input_qi = real_zero @@ -635,6 +529,10 @@ subroutine scm_input_create(scm_input, ntimes, nlev) scm_input%input_v = real_zero scm_input%input_tke = real_zero scm_input%input_ozone = real_zero + + allocate(scm_input%input_pres_i(nlev+1),scm_input%input_pres_l(nlev)) + scm_input%input_pres_i = real_zero + scm_input%input_pres_l = real_zero allocate(scm_input%input_pres_surf(ntimes), & scm_input%input_T_surf(ntimes), scm_input%input_sh_flux_sfc(ntimes), scm_input%input_lh_flux_sfc(ntimes), & @@ -647,6 +545,32 @@ subroutine scm_input_create(scm_input, ntimes, nlev) scm_input%input_T_surf = real_zero scm_input%input_lat = real_zero scm_input%input_lon = real_zero + + allocate(scm_input%input_stc(nsoil), scm_input%input_smc(nsoil), scm_input%input_slc(nsoil)) + scm_input%input_stc = real_zero + scm_input%input_smc = real_zero + scm_input%input_slc = real_zero + + allocate(scm_input%input_snicexy(nsnow),scm_input%input_snliqxy(nsnow), scm_input%input_tsnoxy(nsnow), & + scm_input%input_smoiseq(nsoil), scm_input%input_zsnsoxy(nsnow + nsoil)) + scm_input%input_snicexy = real_zero + scm_input%input_snliqxy = real_zero + scm_input%input_tsnoxy = real_zero + scm_input%input_smoiseq = real_zero + scm_input%input_zsnsoxy = real_zero + + allocate(scm_input%input_tiice(nice)) + scm_input%input_tiice = real_zero + + allocate(scm_input%input_tslb(nsoil), scm_input%input_smois(nsoil), scm_input%input_sh2o(nsoil), & + scm_input%input_smfr(nsoil), scm_input%input_flfr(nsoil)) + scm_input%input_tslb = real_zero + scm_input%input_smois = real_zero + scm_input%input_sh2o = real_zero + scm_input%input_smfr = real_zero + scm_input%input_flfr = real_zero + + scm_input%input_tsfco = real_zero scm_input%input_vegsrc = int_zero scm_input%input_vegtyp = real_zero scm_input%input_soiltyp = real_zero @@ -654,7 +578,7 @@ subroutine scm_input_create(scm_input, ntimes, nlev) scm_input%input_vegfrac = real_zero scm_input%input_shdmin = real_zero scm_input%input_shdmax = real_zero - scm_input%input_zorl = real_zero + scm_input%input_zorlo = real_zero scm_input%input_slmsk = real_zero scm_input%input_canopy = real_zero scm_input%input_hice = real_zero @@ -666,6 +590,45 @@ subroutine scm_input_create(scm_input, ntimes, nlev) scm_input%input_area = real_zero scm_input%input_tg3 = real_zero scm_input%input_uustar = real_zero + scm_input%input_alvsf = real_zero + scm_input%input_alnsf = real_zero + scm_input%input_alvwf = real_zero + scm_input%input_alnwf = real_zero + scm_input%input_facsf = real_zero + scm_input%input_facwf = real_zero + scm_input%input_weasd = real_zero + scm_input%input_f10m = real_zero + scm_input%input_t2m = real_zero + scm_input%input_q2m = real_zero + scm_input%input_ffmm = real_zero + scm_input%input_ffhh = real_zero + scm_input%input_tprcp = real_zero + scm_input%input_srflag = real_zero + scm_input%input_tsfcl = real_zero + scm_input%input_zorll = real_zero + scm_input%input_zorli = real_zero + scm_input%input_zorlw = real_zero + + scm_input%input_stddev = real_zero + scm_input%input_convexity = real_zero + scm_input%input_oa1 = real_zero + scm_input%input_oa2 = real_zero + scm_input%input_oa3 = real_zero + scm_input%input_oa4 = real_zero + scm_input%input_ol1 = real_zero + scm_input%input_ol2 = real_zero + scm_input%input_ol3 = real_zero + scm_input%input_ol4 = real_zero + scm_input%input_theta = real_zero + scm_input%input_gamma = real_zero + scm_input%input_sigma = real_zero + scm_input%input_elvmax = real_zero + scm_input%input_oro = real_zero + scm_input%input_oro_uf = real_zero + scm_input%input_landfrac = real_zero + scm_input%input_lakefrac = real_zero + scm_input%input_lakedepth = real_zero + scm_input%input_tvxy = real_zero scm_input%input_tgxy = real_zero scm_input%input_tahxy = real_zero @@ -695,37 +658,34 @@ subroutine scm_input_create(scm_input, ntimes, nlev) scm_input%input_deeprechxy = real_zero scm_input%input_rechxy = real_zero scm_input%input_snowxy = real_zero - scm_input%input_alvsf = real_zero - scm_input%input_alnsf = real_zero - scm_input%input_alvwf = real_zero - scm_input%input_alnwf = real_zero - scm_input%input_stddev = real_zero - scm_input%input_convexity = real_zero - scm_input%input_oa1 = real_zero - scm_input%input_oa2 = real_zero - scm_input%input_oa3 = real_zero - scm_input%input_oa4 = real_zero - scm_input%input_ol1 = real_zero - scm_input%input_ol2 = real_zero - scm_input%input_ol3 = real_zero - scm_input%input_ol4 = real_zero - scm_input%input_theta = real_zero - scm_input%input_gamma = real_zero - scm_input%input_sigma = real_zero - scm_input%input_elvmax = real_zero - scm_input%input_facsf = real_zero - scm_input%input_facwf = real_zero - scm_input%input_sh_flux_sfc = real_zero - scm_input%input_lh_flux_sfc = real_zero - scm_input%input_w_ls = real_zero - scm_input%input_omega = real_zero - scm_input%input_u_g = real_zero - scm_input%input_v_g = real_zero - scm_input%input_dT_dt_rad = real_zero - scm_input%input_h_advec_thetail = real_zero - scm_input%input_h_advec_qt = real_zero - scm_input%input_v_advec_thetail = real_zero - scm_input%input_v_advec_qt = real_zero + + scm_input%input_tref = real_zero !< sea surface reference temperature for NSST (K) + scm_input%input_z_c = real_zero !< sub-layer cooling thickness for NSST (m) + scm_input%input_c_0 = real_zero !< coefficient 1 to calculate d(Tz)/d(Ts) for NSST + scm_input%input_c_d = real_zero !< coefficient 2 to calculate d(Tz)/d(Ts) for NSST + scm_input%input_w_0 = real_zero !< coefficient 3 to calculate d(Tz)/d(Ts) for NSST + scm_input%input_w_d = real_zero !< coefficient 4 to calculate d(Tz)/d(Ts) for NSST + scm_input%input_xt = real_zero !< heat content in diurnal thermocline layer for NSST (K m) + scm_input%input_xs = real_zero !< salinity content in diurnal thermocline layer for NSST (ppt m) + scm_input%input_xu = real_zero !< u-current in diurnal thermocline layer for NSST (m2 s-1) + scm_input%input_xv = real_zero !< v-current in diurnal thermocline layer for NSST (m2 s-1) + scm_input%input_xz = real_zero !< thickness of diurnal thermocline layer for NSST (m) + scm_input%input_zm = real_zero !< thickness of ocean mixed layer for NSST (m) + scm_input%input_xtts = real_zero !< sensitivity of diurnal thermocline layer heat content to surface temperature [d(xt)/d(ts)] for NSST (m) + scm_input%input_xzts = real_zero !< sensitivity of diurnal thermocline layer thickness to surface temperature [d(xz)/d(ts)] for NSST (m K-1) + scm_input%input_d_conv = real_zero !< thickness of free convection layer for NSST (m) + scm_input%input_ifd = real_zero !< index to start DTM run for NSST + scm_input%input_dt_cool = real_zero !< sub-layer cooling amount for NSST (K) + scm_input%input_qrain = real_zero !< sensible heat due to rainfall for NSST (W) + + scm_input%input_wetness = real_zero !< normalized soil wetness for RUC LSM + scm_input%input_clw_surf = real_zero !< cloud condensed water mixing ratio at surface for RUC LSM (kg kg-1) + scm_input%input_qwv_surf = real_zero !< water vapor mixing ratio at surface for RUC LSM (kg kg-1) + scm_input%input_tsnow = real_zero !< snow temperature at the bottom of the first snow layer for RUC LSM (K) + scm_input%input_snowfallac = real_zero !< run-total snow accumulation on the ground for RUC LSM (kg m-2) + scm_input%input_acsnow = real_zero !< snow water equivalent of run-total frozen precip for RUC LSM (kg m-2) + scm_input%input_lai = real_zero !< leaf area index for RUC LSM + scm_input%input_sh_flux_sfc = real_zero scm_input%input_lh_flux_sfc = real_zero scm_input%input_w_ls = real_zero @@ -821,53 +781,25 @@ subroutine physics_associate(physics, scm_state) physics%Statein%vvl => scm_state%omega physics%Statein%tgrs => scm_state%state_T(:,:,1) physics%Statein%qgrs => scm_state%state_tracer(:,:,:,1) - - physics%Sfcprop%tsfc => scm_state%T_surf - physics%Sfcprop%tref => scm_state%T_surf - physics%Sfcprop%tsfco => scm_state%T_surf - physics%Sfcprop%tisfc => scm_state%T_surf !sea ice temperature is the same as SST for now? - physics%Sfcprop%slmsk => scm_state%sfc_type - physics%Sfcprop%vtype => scm_state%veg_type - physics%Sfcprop%vfrac => scm_state%veg_frac - physics%Sfcprop%slope => scm_state%slope_type - physics%Interstitial%sigmaf = min(scm_state%veg_frac,0.01) - physics%Sfcprop%shdmax => scm_state%shdmax - physics%Sfcprop%shdmin => scm_state%shdmin - physics%Sfcprop%tg3 => scm_state%tg3 - physics%Sfcprop%uustar => scm_state%uustar - physics%Sfcprop%stype => scm_state%soil_type - physics%Sfcprop%alvsf => scm_state%alvsf - physics%Sfcprop%alnsf => scm_state%alnsf - physics%Sfcprop%hprime=> scm_state%hprime - physics%Sfcprop%alvwf => scm_state%alvwf - physics%Sfcprop%alnwf => scm_state%alnwf - physics%Sfcprop%facsf => scm_state%facsf - physics%Sfcprop%facwf => scm_state%facwf - - physics%Sfcprop%stc => scm_state%stc(:,:,1) - physics%Sfcprop%smc => scm_state%smc(:,:,1) - physics%Sfcprop%slc => scm_state%slc(:,:,1) - - !GJF : the following logic was introduced into FV3GFS_io.F90 as part of the fractional landmask update (additional logic exists in the same file if the fractional landmask is actually used!) - do i = 1, scm_state%n_cols - if (physics%Sfcprop%slmsk(i) < 0.1 .or. physics%Sfcprop%slmsk(i) > 1.9) then - physics%Sfcprop%landfrac(i) = 0.0 - if (physics%Sfcprop%oro_uf(i) > 0.01) then - physics%Sfcprop%lakefrac(i) = 1.0 ! lake - else - physics%Sfcprop%lakefrac(i) = 0.0 ! ocean - endif - else - physics%Sfcprop%landfrac(i) = 1.0 - end if - if (physics%Sfcprop%lakefrac(i) > 0.0) then - physics%Sfcprop%oceanfrac(i) = 0.0 ! lake & ocean don't coexist in a cell, lake dominates - else - physics%Sfcprop%oceanfrac(i) = 1.0 - physics%Sfcprop%landfrac(i) !LHS:ocean frac [0:1] - end if - end do - !GJF - + + if (.not. scm_state%model_ics .and. .not. scm_state%sfc_flux_spec) then + do i =1, physics%Model%ncols + if (scm_state%sfc_type(i) == 0) then + physics%Sfcprop%tsfco => scm_state%T_surf + elseif (scm_state%sfc_type(i) == 1) then + physics%Sfcprop%tsfcl => scm_state%T_surf + elseif (scm_state%sfc_type(i) == 2) then + physics%Sfcprop%tisfc => scm_state%T_surf + end if + end do + end if + + if (scm_state%sfc_flux_spec) then + do i =1, physics%Model%ncols + physics%Sfcprop%tsfc => scm_state%T_surf + end do + end if + if(scm_state%time_scheme == 2) then physics%Stateout%gu0 => scm_state%state_u(:,:,2) physics%Stateout%gv0 => scm_state%state_v(:,:,2) @@ -880,61 +812,743 @@ subroutine physics_associate(physics, scm_state) physics%Stateout%gq0 => scm_state%state_tracer(:,:,:,1) endif - physics%Sfcprop%zorl => scm_state%sfc_roughness_length_cm if(scm_state%sfc_flux_spec) then physics%Sfcprop%spec_sh_flux => scm_state%sh_flux physics%Sfcprop%spec_lh_flux => scm_state%lh_flux endif - if(scm_state%model_ics) then - !physics%Sfcprop%zorl => scm_state%sfc_roughness_length_cm - physics%Sfcprop%canopy => scm_state%canopy - physics%Sfcprop%hice => scm_state%hice - physics%Sfcprop%fice => scm_state%fice - physics%Sfcprop%tisfc => scm_state%tisfc - physics%Sfcprop%snowd => scm_state%snwdph - physics%Sfcprop%snoalb => scm_state%snoalb - physics%Sfcprop%sncovr => scm_state%sncovr - if (physics%Model%lsm == physics%Model%lsm_noahmp) then - physics%Sfcprop%snowxy => scm_state%snowxy - physics%Sfcprop%tvxy => scm_state%tvxy - physics%Sfcprop%tgxy => scm_state%tgxy - physics%Sfcprop%canicexy => scm_state%canicexy - physics%Sfcprop%canliqxy => scm_state%canliqxy - physics%Sfcprop%eahxy => scm_state%eahxy - physics%Sfcprop%tahxy => scm_state%tahxy - physics%Sfcprop%cmxy => scm_state%cmxy - physics%Sfcprop%chxy => scm_state%chxy - physics%Sfcprop%fwetxy => scm_state%fwetxy - physics%Sfcprop%sneqvoxy => scm_state%sneqvoxy - physics%Sfcprop%alboldxy => scm_state%alboldxy - physics%Sfcprop%qsnowxy => scm_state%qsnowxy - physics%Sfcprop%wslakexy => scm_state%wslakexy - physics%Sfcprop%zwtxy => scm_state%zwtxy - physics%Sfcprop%waxy => scm_state%waxy - physics%Sfcprop%wtxy => scm_state%wtxy - physics%Sfcprop%lfmassxy => scm_state%lfmassxy - physics%Sfcprop%rtmassxy => scm_state%rtmassxy - physics%Sfcprop%stmassxy => scm_state%stmassxy - physics%Sfcprop%woodxy => scm_state%woodxy - physics%Sfcprop%stblcpxy => scm_state%stblcpxy - physics%Sfcprop%fastcpxy => scm_state%fastcpxy - physics%Sfcprop%xsaixy => scm_state%xsaixy - physics%Sfcprop%xlaixy => scm_state%xlaixy - physics%Sfcprop%taussxy => scm_state%taussxy - physics%Sfcprop%smcwtdxy => scm_state%smcwtdxy - physics%Sfcprop%deeprechxy => scm_state%deeprechxy - physics%Sfcprop%rechxy => scm_state%rechxy - - physics%Sfcprop%snicexy => scm_state%snicexy - physics%Sfcprop%snliqxy => scm_state%snliqxy - physics%Sfcprop%tsnoxy => scm_state%tsnoxy - physics%Sfcprop%smoiseq => scm_state%smoiseq - physics%Sfcprop%zsnsoxy => scm_state%zsnsoxy + + end subroutine physics_associate + + subroutine physics_set(physics, scm_input, scm_state) + !used for initializing variables in the physics DDT (but not pointer association); + !this should be utilized for variables that cannot be modified or "forced" by the SCM; + !most of this routine follows what is in FV3/io/FV3GFS_io.F90/sfc_prop_restart_read + use gmtb_scm_physical_constants, only: con_tice + use data_qc, only: conditionally_set_var + use NetCDF_read, only: missing_value + use noahmp_tables, only: laim_table,saim_table,sla_table, & + bexp_table,smcmax_table,smcwlt_table, & + dwsat_table,dksat_table,psisat_table, & + isurban_table,isbarren_table, & + isice_table,iswater_table + + class(physics_type) :: physics + type(scm_input_type), intent(in) :: scm_input + type(scm_state_type), intent(in) :: scm_state + + integer :: i,j,n,vegtyp,imn,lsoil,isnow,ns,soiltyp,iter + real(kind=dp) :: tem1, tem, masslai, masssai, snd + logical :: missing_var(100) + logical :: cold_start_noahmp + real(kind=dp), parameter:: drythresh = 1.e-4 + + real(kind=kind_phys) :: ddz,expon,aa,bb,smc,func,dfunc,dx + real(kind=kind_phys) :: bexp, smcmax, smcwlt,dwsat,dksat,psisat + + real(kind=kind_phys), dimension(-2:0) :: dzsno + real(kind=kind_phys), dimension(-2:4) :: dzsnso + + real(kind=kind_phys), dimension(4), save :: zsoil,dzs + data dzs / 0.1_dp, 0.3_dp, 0.6_dp, 1.0_dp/ + data zsoil /-0.1_dp,-0.4_dp,-1.0_dp,-2.0_dp/ + + cold_start_noahmp = .false. + + !double check under what circumstances these should actually be set from input!!! (these overwrite the initialzation in GFS_typedefs) + missing_var = .false. + do i = 1, physics%Model%ncols + if (scm_state%model_ics) then + write(0,'(a)') "Setting internal physics variables from the orographic section of the case input file (scalars)..." + call conditionally_set_var(scm_input%input_stddev, physics%Sfcprop%hprime(i,1), "stddev", .true., missing_var(1)) + call conditionally_set_var(scm_input%input_convexity, physics%Sfcprop%hprime(i,2), "convexity", .true., missing_var(2)) + call conditionally_set_var(scm_input%input_oa1, physics%Sfcprop%hprime(i,3), "oa1", .true., missing_var(3)) + call conditionally_set_var(scm_input%input_oa2, physics%Sfcprop%hprime(i,4), "oa2", .true., missing_var(4)) + call conditionally_set_var(scm_input%input_oa3, physics%Sfcprop%hprime(i,5), "oa3", .true., missing_var(5)) + call conditionally_set_var(scm_input%input_oa4, physics%Sfcprop%hprime(i,6), "oa4", .true., missing_var(6)) + call conditionally_set_var(scm_input%input_ol1, physics%Sfcprop%hprime(i,7), "ol1", .true., missing_var(7)) + call conditionally_set_var(scm_input%input_ol2, physics%Sfcprop%hprime(i,8), "ol2", .true., missing_var(8)) + call conditionally_set_var(scm_input%input_ol3, physics%Sfcprop%hprime(i,9), "ol3", .true., missing_var(9)) + call conditionally_set_var(scm_input%input_ol4, physics%Sfcprop%hprime(i,10), "ol4", .true., missing_var(10)) + call conditionally_set_var(scm_input%input_theta, physics%Sfcprop%hprime(i,11), "theta", .true., missing_var(11)) + call conditionally_set_var(scm_input%input_gamma, physics%Sfcprop%hprime(i,12), "gamma", .true., missing_var(12)) + call conditionally_set_var(scm_input%input_sigma, physics%Sfcprop%hprime(i,13), "sigma", .true., missing_var(13)) + call conditionally_set_var(scm_input%input_elvmax, physics%Sfcprop%hprime(i,14), "elvmax", .true., missing_var(14)) + call conditionally_set_var(scm_input%input_oro, physics%Sfcprop%oro(i), "oro", .true., missing_var(15)) + call conditionally_set_var(scm_input%input_oro_uf, physics%Sfcprop%oro_uf(i), "oro_uf", (physics%Model%do_ugwp .and. physics%Model%nmtvr == 14), missing_var(16)) + call conditionally_set_var(scm_input%input_landfrac, physics%Sfcprop%landfrac(i), "landfrac", physics%Model%frac_grid, missing_var(17)) + call conditionally_set_var(scm_input%input_lakefrac, physics%Sfcprop%lakefrac(i), "lakefrac", (physics%Model%lkm == 1), missing_var(18)) + call conditionally_set_var(scm_input%input_lakedepth, physics%Sfcprop%lakedepth(i), "lakedepth", (physics%Model%lkm == 1), missing_var(19)) + + !write out warning if missing data for non-required variables + n = 19 + if ( i==1 .and. ANY( missing_var(1:n) ) ) then + write(0,'(a)') "INPUT CHECK: Some missing input data was found related to (potentially non-required) orography and gravity wave drag parameters. This may lead to crashes or other strange behavior." + write(0,'(a)') "Check gmtb_scm_type_defs.F90/physics_set to see the names of variables that are missing, corresponding to the following indices:" + do j=1, n + if (missing_var(j)) write(0,'(a,i0)') "variable index ",j + end do + end if + missing_var = .false. + + write(0,'(a)') "Setting internal physics variables from the surface section of the case input file (scalars)..." + call conditionally_set_var(scm_input%input_slmsk, physics%Sfcprop%slmsk(i), "slmsk", (.not. physics%Model%frac_grid), missing_var(1)) + call conditionally_set_var(scm_input%input_tsfco, physics%Sfcprop%tsfco(i), "tsfco", .true., missing_var(2)) + call conditionally_set_var(scm_input%input_weasd, physics%Sfcprop%weasd(i), "weasd", .false., missing_var(3)) + call conditionally_set_var(scm_input%input_tg3, physics%Sfcprop%tg3(i), "tg3", .true., missing_var(4)) + call conditionally_set_var(scm_input%input_zorlo, physics%Sfcprop%zorlo(i), "zorlo", .true., missing_var(5)) + call conditionally_set_var(scm_input%input_alvsf, physics%Sfcprop%alvsf(i), "alvsf", .true., missing_var(6)) + call conditionally_set_var(scm_input%input_alnsf, physics%Sfcprop%alnsf(i), "alnsf", .true., missing_var(7)) + call conditionally_set_var(scm_input%input_alvwf, physics%Sfcprop%alvwf(i), "alvwf", .true., missing_var(8)) + call conditionally_set_var(scm_input%input_alnwf, physics%Sfcprop%alnwf(i), "alnwf", .true., missing_var(9)) + call conditionally_set_var(scm_input%input_facsf, physics%Sfcprop%facsf(i), "facsf", .true., missing_var(10)) + call conditionally_set_var(scm_input%input_facwf, physics%Sfcprop%facwf(i), "facwf", .true., missing_var(11)) + call conditionally_set_var(scm_input%input_vegfrac, physics%Sfcprop%vfrac(i), "vegfrac", .true., missing_var(12)) + !GJF: is this needed anymore (not in FV3GFS_io)? + physics%Interstitial%sigmaf(i) = min(physics%Sfcprop%vfrac(i),0.01) + call conditionally_set_var(scm_input%input_canopy, physics%Sfcprop%canopy(i), "canopy", .true., missing_var(13)) + call conditionally_set_var(scm_input%input_f10m, physics%Sfcprop%f10m(i), "f10m", .false., missing_var(14)) + call conditionally_set_var(scm_input%input_t2m, physics%Sfcprop%t2m(i), "t2m", physics%Model%cplflx, missing_var(15)) + call conditionally_set_var(scm_input%input_q2m, physics%Sfcprop%q2m(i), "q2m", physics%Model%cplflx, missing_var(16)) + call conditionally_set_var(scm_input%input_vegtyp, physics%Sfcprop%vtype(i), "vegtyp", .true., missing_var(17)) + call conditionally_set_var(scm_input%input_soiltyp, physics%Sfcprop%stype(i), "soiltyp", .true., missing_var(18)) + call conditionally_set_var(scm_input%input_uustar, physics%Sfcprop%uustar(i), "uustar", .true., missing_var(19)) + call conditionally_set_var(scm_input%input_ffmm, physics%Sfcprop%ffmm(i), "ffmm", .false., missing_var(20)) + call conditionally_set_var(scm_input%input_ffhh, physics%Sfcprop%ffhh(i), "ffhh", .false., missing_var(21)) + call conditionally_set_var(scm_input%input_hice, physics%Sfcprop%hice(i), "hice", .true., missing_var(22)) + call conditionally_set_var(scm_input%input_fice, physics%Sfcprop%fice(i), "fice", .true., missing_var(23)) + call conditionally_set_var(scm_input%input_tisfc, physics%Sfcprop%tisfc(i), "tisfc", .true., missing_var(24)) + call conditionally_set_var(scm_input%input_tprcp, physics%Sfcprop%tprcp(i), "tprcp", .false., missing_var(25)) + call conditionally_set_var(scm_input%input_srflag, physics%Sfcprop%srflag(i), "srflag", .false., missing_var(26)) + call conditionally_set_var(scm_input%input_snwdph, physics%Sfcprop%snowd(i), "snwdph", .false., missing_var(27)) + call conditionally_set_var(scm_input%input_shdmin, physics%Sfcprop%shdmin(i), "shdmin", .true., missing_var(28)) + call conditionally_set_var(scm_input%input_shdmax, physics%Sfcprop%shdmax(i), "shdmax", .true., missing_var(29)) + call conditionally_set_var(scm_input%input_slopetype, physics%Sfcprop%slope(i), "slopetyp", .true., missing_var(30)) + call conditionally_set_var(scm_input%input_snoalb, physics%Sfcprop%snoalb(i), "snoalb", .true., missing_var(31)) + call conditionally_set_var(scm_input%input_sncovr, physics%Sfcprop%sncovr(i), "sncovr", .false., missing_var(32)) + call conditionally_set_var(scm_input%input_tsfcl, physics%Sfcprop%tsfcl(i), "tsfcl", .false., missing_var(33)) + call conditionally_set_var(scm_input%input_zorll, physics%Sfcprop%zorll(i), "zorll", .false., missing_var(34)) + call conditionally_set_var(scm_input%input_zorli, physics%Sfcprop%zorli(i), "zorli", .false., missing_var(35)) + if (physics%Model%cplwav) then + call conditionally_set_var(scm_input%input_zorlw, physics%Sfcprop%zorlw(i), "zorlw", .true., missing_var(36)) + else + physics%Sfcprop%zorlw(i) = physics%Sfcprop%zorlo(i) + end if + + !GJF: computing sncovr if model_ics and sncovr is missing: (this requires data from namelist_soilveg module, which is not set until after the CCPP physics_init stage) + !this should happen inside of the init stage of a CCPP scheme (after the LSM is initialized -- or at least after set_soilveg is called) + !For now, just set sncovr = 0 is missing + if (missing_var(32)) physics%Sfcprop%sncovr(i) = real_zero + ! physics%Sfcprop%sncovr(i) = zero + ! if (physics%Sfcprop%landfrac(i) >= drythresh .or. physics%Sfcprop%fice(i) >= physics%Model%min_seaice) then + ! vegtyp = physics%Sfcprop%vtype(i) + ! if (vegtyp == 0) vegtyp = 7 + ! rsnow = 0.001_dp*physics%Sfcprop%weasd(i)/snupx(vegtyp) + ! if (0.001_dp*physics%Sfcprop%weasd(i) < snupx(vegtyp)) then + ! physics%Sfcprop%sncovr(i) = real_one - (exp(-salp_data*rsnow) - rsnow*exp(-salp_data)) + ! else + ! physics%Sfcprop%sncovr(i) = real_one + ! endif + ! endif + + !write out warning if missing data for non-required variables + n = 36 + if ( i==1 .and. ANY( missing_var(1:n) ) ) then + write(0,'(a)') "INPUT CHECK: Some missing input data was found related to (potentially non-required) surface variables. This may lead to crashes or other strange behavior." + write(0,'(a)') "Check gmtb_scm_type_defs.F90/physics_set to see the names of variables that are missing, corresponding to the following indices:" + do j=1, n + if (missing_var(j)) write(0,'(a,i0)') "variable index ",j + end do + end if + missing_var = .false. + else + physics%Sfcprop%slmsk(i) = scm_state%sfc_type(i) + ! tsfco is already pointing to T_surf forcing in physics_associate + ! physics%Sfcprop%tsfco(i) => scm_state%T_surf + physics%Sfcprop%zorlo(i) = scm_state%sfc_roughness_length_cm(i) + ! tisfc is already pointing to T_surf forcing in physics_associate + ! physics%Sfcprop%tisfc(i) => scm_state%T_surf + ! tsfcl is already pointing to T_surf forcing in physics_associate + ! physics%Sfcprop%tsfcl(i) => scm_state%T_surf + if (physics%Sfcprop%slmsk(i) > 1.9_dp) physics%Sfcprop%fice(i) = 1.0 !needed to calculate tsfc and zorl below when model_ics == .false. + end if + + !this overwrites what is in the suite namelist file -- is that desirable? + !if (scm_state%model_ics) then + ! physics%Model%ivegsrc = scm_input%input_vegsrc + !end if + + if(scm_state%model_ics .and. physics%Model%frac_grid) then ! obtain slmsk from landfrac + !! next 5 lines are temporary till lake model is available + !if (physics%Sfcprop%lakefrac(i) > real_zero) then + ! !physics%Sfcprop%lakefrac(i) = nint(physics%Sfcprop%lakefrac(i)) + ! physics%Sfcprop%landfrac(i) = real_one - physics%Sfcprop%lakefrac(i) + ! if (physics%Sfcprop%lakefrac(i) == real_zero) physics%Sfcprop%fice(i) = real_zero + !endif + !physics%Sfcprop%slmsk(i) = ceiling(physics%Sfcprop%landfrac(i)) + !if (physics%Sfcprop%fice(i) > physics%Model%min_lakeice .and. physics%Sfcprop%landfrac(i) == real_zero) physics%Sfcprop%slmsk(i) = 2 ! land dominates ice if co-exist + physics%Sfcprop%slmsk(i) = ceiling(physics%Sfcprop%landfrac(i)) !nint/floor are options + else ! obtain landfrac from slmsk + if (physics%Sfcprop%slmsk(i) > 1.9_dp) then + physics%Sfcprop%landfrac(i) = real_zero + else + physics%Sfcprop%landfrac(i) = physics%Sfcprop%slmsk(i) + endif endif - endif + + if (physics%Sfcprop%lakefrac(i) > real_zero) then + physics%Sfcprop%oceanfrac(i) = real_zero ! lake & ocean don't coexist in a cell + !GJF: the following code was taken from fv3atm/develop branch commit d10e450 from 2020/12/2, but the if statements seem like they are in error (reproduced here) + if (physics%Sfcprop%slmsk(i) /= real_one) then + if (physics%Sfcprop%fice(i) >= physics%Model%min_lakeice) then + if (physics%Sfcprop%slmsk(i) < 1.9_dp) & + write(*,'(a,2i3,3f6.2)') 'reset lake slmsk=2 at i=' & + ,i,physics%Sfcprop%fice(i),physics%Sfcprop%slmsk(i),physics%Sfcprop%lakefrac(i) + physics%Sfcprop%slmsk(i) = 2. + else if (physics%Sfcprop%slmsk(i) > 1.e-7) then + write(*,'(a,2i3,3f6.2)') 'reset lake slmsk=0 at i=' & + ,i,physics%Sfcprop%fice(i),physics%Sfcprop%slmsk(i),physics%Sfcprop%lakefrac(i) + physics%Sfcprop%slmsk(i) = real_zero + end if + end if + !if (physics%Sfcprop%fice(i) < physics%Model%min_lakeice) then + ! physics%Sfcprop%fice(i) = real_zero + ! if (physics%Sfcprop%slmsk(i) == 2) physics%Sfcprop%slmsk(i) = 0 + !endif + else + physics%Sfcprop%oceanfrac(i) = real_one - physics%Sfcprop%landfrac(i) + if (physics%Sfcprop%slmsk(i) /= real_one) then + if (physics%Sfcprop%fice(i) >= physics%Model%min_seaice) then + if (physics%Sfcprop%slmsk(i) < 1.9_dp) & + write(*,'(a,2i3,3f6.2)') 'reset sea slmsk=2 at i,=' & + ,i,physics%Sfcprop%fice(i),physics%Sfcprop%slmsk(i),physics%Sfcprop%landfrac(i) + physics%Sfcprop%slmsk(i) = 2. + else if (physics%Sfcprop%slmsk(i) > 1.e-7) then + write(*,'(a,2i3,4f6.2)') 'reset sea slmsk=0 at i=' & + ,i,physics%Sfcprop%fice(i),physics%Sfcprop%slmsk(i),physics%Sfcprop%landfrac(i) + physics%Sfcprop%slmsk(i) = real_zero + end if + end if + ! if (physics%Sfcprop%fice(i) < physics%Model%min_seaice) then + ! physics%Sfcprop%fice(i) = real_zero + ! if (physics%Sfcprop%slmsk(i) == 2) physics%Sfcprop%slmsk(i) = 0 + ! endif + endif + + !--- NSSTM variables + if (physics%Model%nstf_name(1) > 0) then + if (physics%Model%nstf_name(2) == 1 .or. .not. scm_state%model_ics) then ! nsst spinup + !--- nsstm tref + physics%Sfcprop%tref(i) = physics%Sfcprop%tsfco(i) + physics%Sfcprop%z_c(i) = real_zero + physics%Sfcprop%c_0(i) = real_zero + physics%Sfcprop%c_d(i) = real_zero + physics%Sfcprop%w_0(i) = real_zero + physics%Sfcprop%w_d(i) = real_zero + physics%Sfcprop%xt(i) = real_zero + physics%Sfcprop%xs(i) = real_zero + physics%Sfcprop%xu(i) = real_zero + physics%Sfcprop%xv(i) = real_zero + physics%Sfcprop%xz(i) = 30.0_dp + physics%Sfcprop%zm(i) = real_zero + physics%Sfcprop%xtts(i) = real_zero + physics%Sfcprop%xzts(i) = real_zero + physics%Sfcprop%d_conv(i) = real_zero + physics%Sfcprop%ifd(i) = real_zero + physics%Sfcprop%dt_cool(i) = real_zero + physics%Sfcprop%qrain(i) = real_zero + elseif (physics%Model%nstf_name(2) == 0) then ! nsst restart + write(0,'(a)') "Setting internal physics variables from the NSST section of the case input file (scalars)..." + call conditionally_set_var(scm_input%input_tref, physics%Sfcprop%tref(i), "tref", .true., missing_var(1)) + call conditionally_set_var(scm_input%input_z_c, physics%Sfcprop%z_c(i), "z_c", .true., missing_var(2)) + call conditionally_set_var(scm_input%input_c_0, physics%Sfcprop%c_0(i), "c_0", .true., missing_var(3)) + call conditionally_set_var(scm_input%input_c_d, physics%Sfcprop%c_d(i), "c_d", .true., missing_var(4)) + call conditionally_set_var(scm_input%input_w_0, physics%Sfcprop%w_0(i), "w_0", .true., missing_var(5)) + call conditionally_set_var(scm_input%input_w_d, physics%Sfcprop%w_d(i), "w_d", .true., missing_var(6)) + call conditionally_set_var(scm_input%input_xt, physics%Sfcprop%xt(i), "xt", .true., missing_var(7)) + call conditionally_set_var(scm_input%input_xs, physics%Sfcprop%xs(i), "xs", .true., missing_var(8)) + call conditionally_set_var(scm_input%input_xu, physics%Sfcprop%xu(i), "xu", .true., missing_var(9)) + call conditionally_set_var(scm_input%input_xv, physics%Sfcprop%xv(i), "xv", .true., missing_var(10)) + call conditionally_set_var(scm_input%input_xz, physics%Sfcprop%xz(i), "xz", .true., missing_var(11)) + call conditionally_set_var(scm_input%input_zm, physics%Sfcprop%zm(i), "zm", .true., missing_var(12)) + call conditionally_set_var(scm_input%input_xtts, physics%Sfcprop%xtts(i), "xtts", .true., missing_var(13)) + call conditionally_set_var(scm_input%input_xzts, physics%Sfcprop%xzts(i), "xzts", .true., missing_var(14)) + call conditionally_set_var(scm_input%input_d_conv, physics%Sfcprop%d_conv(i), "d_conv", .true., missing_var(15)) + call conditionally_set_var(scm_input%input_ifd, physics%Sfcprop%ifd(i), "ifd", .true., missing_var(16)) + call conditionally_set_var(scm_input%input_dt_cool, physics%Sfcprop%dt_cool(i), "dt_cool", .true., missing_var(17)) + call conditionally_set_var(scm_input%input_qrain, physics%Sfcprop%qrain(i), "qrain", .true., missing_var(18)) + ! all NNST variables are required when NNST spin-up is off, so no need to write out warning for missing data (the model would have already stopped) + missing_var = .false. + endif + endif + + if (scm_state%model_ics .and. physics%Model%lsm == physics%Model%lsm_ruc) then !.and. warm_start (not implemented here -- assuming that RUC LSM has warm start data from file) + !--- Extra RUC LSM variables + write(0,'(a)') "Setting internal physics variables from the RUC LSM section of the case input file (scalars)..." + call conditionally_set_var(scm_input%input_wetness, physics%Sfcprop%wetness(i), "wetness", .false., missing_var(1)) + call conditionally_set_var(scm_input%input_clw_surf, physics%Sfcprop%clw_surf(i), "clw_surf", .false., missing_var(2)) + call conditionally_set_var(scm_input%input_qwv_surf, physics%Sfcprop%qwv_surf(i), "qwv_surf", .false., missing_var(3)) + call conditionally_set_var(scm_input%input_tsnow, physics%Sfcprop%tsnow(i), "tsnow", .false., missing_var(4)) + call conditionally_set_var(scm_input%input_snowfallac, physics%Sfcprop%snowfallac(i), "snowfallac", .false., missing_var(5)) + call conditionally_set_var(scm_input%input_acsnow, physics%Sfcprop%acsnow(i), "acsnow", .false., missing_var(6)) + if (physics%Model%lsm == physics%Model%lsm_ruc .and. physics%Model%rdlai) then + !when rdlai = T, RUC LSM expects the LAI to be read in, hence the required variable attribute below + call conditionally_set_var(scm_input%input_lai, physics%Sfcprop%xlaixy(i), "lai", .true., missing_var(7)) + end if + + !write out warning if missing data for non-required variables + n = 7 + if ( i==1 .and. ANY( missing_var(1:n) ) ) then + write(0,'(a)') "INPUT CHECK: Some missing input data was found related to (potentially non-required) surface variables for RUC LSM. This may lead to crashes or other strange behavior." + write(0,'(a)') "Check gmtb_scm_type_defs.F90/physics_set to see the names of variables that are missing, corresponding to the following indices:" + do j=1, n + if (missing_var(j)) write(0,'(a,i0)') "variable index ",j + end do + end if + missing_var = .false. + elseif (scm_state%model_ics .and. physics%Model%lsm == physics%Model%lsm_noahmp) then + write(0,'(a)') "Setting internal physics variables from the NoahMP section of the case input file (scalars)..." + !all of these can be missing, since a method exists to "cold start" these variables + call conditionally_set_var(scm_input%input_snowxy, physics%Sfcprop%snowxy(i), "snowxy", .false., missing_var(1)) + call conditionally_set_var(scm_input%input_tvxy, physics%Sfcprop%tvxy(i), "tvxy", .false., missing_var(2)) + call conditionally_set_var(scm_input%input_tgxy, physics%Sfcprop%tgxy(i), "tgxy", .false., missing_var(3)) + call conditionally_set_var(scm_input%input_canicexy, physics%Sfcprop%canicexy(i), "canicexy", .false., missing_var(4)) + call conditionally_set_var(scm_input%input_canliqxy, physics%Sfcprop%canliqxy(i), "canliqxy", .false., missing_var(5)) + call conditionally_set_var(scm_input%input_eahxy, physics%Sfcprop%eahxy(i), "eahxy", .false., missing_var(6)) + call conditionally_set_var(scm_input%input_tahxy, physics%Sfcprop%tahxy(i), "tahxy", .false., missing_var(7)) + call conditionally_set_var(scm_input%input_cmxy, physics%Sfcprop%cmxy(i), "cmxy", .false., missing_var(8)) + call conditionally_set_var(scm_input%input_chxy, physics%Sfcprop%chxy(i), "chxy", .false., missing_var(9)) + call conditionally_set_var(scm_input%input_fwetxy, physics%Sfcprop%fwetxy(i), "fwetxy", .false., missing_var(10)) + call conditionally_set_var(scm_input%input_sneqvoxy, physics%Sfcprop%sneqvoxy(i), "sneqvoxy", .false., missing_var(11)) + call conditionally_set_var(scm_input%input_alboldxy, physics%Sfcprop%alboldxy(i), "alboldxy", .false., missing_var(12)) + call conditionally_set_var(scm_input%input_qsnowxy, physics%Sfcprop%qsnowxy(i), "qsnowxy", .false., missing_var(13)) + call conditionally_set_var(scm_input%input_wslakexy, physics%Sfcprop%wslakexy(i), "wslakexy", .false., missing_var(14)) + call conditionally_set_var(scm_input%input_zwtxy, physics%Sfcprop%zwtxy(i), "zwtxy", .false., missing_var(15)) + call conditionally_set_var(scm_input%input_waxy, physics%Sfcprop%waxy(i), "waxy", .false., missing_var(16)) + call conditionally_set_var(scm_input%input_wtxy, physics%Sfcprop%wtxy(i), "wtxy", .false., missing_var(17)) + call conditionally_set_var(scm_input%input_lfmassxy, physics%Sfcprop%lfmassxy(i), "lfmassxy", .false., missing_var(18)) + call conditionally_set_var(scm_input%input_rtmassxy, physics%Sfcprop%rtmassxy(i), "rtmassxy", .false., missing_var(19)) + call conditionally_set_var(scm_input%input_stmassxy, physics%Sfcprop%stmassxy(i), "stmassxy", .false., missing_var(20)) + call conditionally_set_var(scm_input%input_woodxy, physics%Sfcprop%woodxy(i), "woodxy", .false., missing_var(21)) + call conditionally_set_var(scm_input%input_stblcpxy, physics%Sfcprop%stblcpxy(i), "stblcpxy", .false., missing_var(22)) + call conditionally_set_var(scm_input%input_fastcpxy, physics%Sfcprop%fastcpxy(i), "fastcpxy", .false., missing_var(23)) + call conditionally_set_var(scm_input%input_xsaixy, physics%Sfcprop%xsaixy(i), "xsaixy", .false., missing_var(24)) + call conditionally_set_var(scm_input%input_xlaixy, physics%Sfcprop%xlaixy(i), "xlaixy", .false., missing_var(25)) + call conditionally_set_var(scm_input%input_taussxy, physics%Sfcprop%taussxy(i), "taussxy", .false., missing_var(26)) + call conditionally_set_var(scm_input%input_smcwtdxy, physics%Sfcprop%smcwtdxy(i), "smcwtdxy", .false., missing_var(27)) + call conditionally_set_var(scm_input%input_deeprechxy, physics%Sfcprop%deeprechxy(i), "deeprechxy", .false., missing_var(28)) + call conditionally_set_var(scm_input%input_rechxy, physics%Sfcprop%rechxy(i), "rechxy", .false., missing_var(29)) + + !write out warning if missing data for non-required variables + n = 29 + if ( i==1 .and. ANY( missing_var(1:n) ) ) then + cold_start_noahmp = .true. + write(0,'(a)') "INPUT CHECK: Some missing input data was found related to surface variables for NoahMP LSM. Due to this, a cold-start algorithm to initialize variables will be used." + write(0,'(a)') "Check gmtb_scm_type_defs.F90/physics_set to see the names of variables that are missing, corresponding to the following indices:" + do j=1, n + if (missing_var(j)) write(0,'(a,i0)') "variable index ",j + end do + end if + missing_var = .false. + end if + + if (scm_state%model_ics .and. (physics%Model%lsm == physics%Model%lsm_noah .or. & + physics%Model%lsm == physics%Model%lsm_noahmp)) then !.or. (.not. warm_start) from FV3GFS_io is not implemented + !check for nonmissing values + !--- 3D variables + ! do lsoil = 1,physics%Model%lsoil + ! physics%Sfcprop%stc(i,lsoil) = scm_input%input_stc(lsoil) !--- stc + ! physics%Sfcprop%smc(i,lsoil) = scm_input%input_smc(lsoil) !--- smc + ! physics%Sfcprop%slc(i,lsoil) = scm_input%input_slc(lsoil) !--- slc + ! enddo + call conditionally_set_var(scm_input%input_stc(:), physics%Sfcprop%stc(i,:), "stc", .true., missing_var(1)) + call conditionally_set_var(scm_input%input_smc(:), physics%Sfcprop%smc(i,:), "smc", .true., missing_var(2)) + call conditionally_set_var(scm_input%input_slc(:), physics%Sfcprop%slc(i,:), "slc", .true., missing_var(3)) + + if (physics%Model%lsm == physics%Model%lsm_noahmp) then + ! do lsoil = -2, 0 + ! physics%Sfcprop%snicexy(i,lsoil) = scm_input%input_snicexy(lsoil) + ! physics%Sfcprop%snliqxy(i,lsoil) = scm_input%input_snliqxy(lsoil) + ! physics%Sfcprop%tsnoxy(i,lsoil) = scm_input%input_tsnoxy(lsoil) + ! enddo + call conditionally_set_var(scm_input%input_snicexy(:), physics%Sfcprop%snicexy(i,:), "snicexy", .false., missing_var(4)) + call conditionally_set_var(scm_input%input_snliqxy(:), physics%Sfcprop%snliqxy(i,:), "sliqexy", .false., missing_var(5)) + call conditionally_set_var(scm_input%input_tsnoxy(:), physics%Sfcprop%tsnoxy(i,:), "tsnoxy", .false., missing_var(6)) + + ! do lsoil = 1, 4 + ! physics%Sfcprop%smoiseq(i,lsoil) = scm_input%input_smoiseq(lsoil) + ! enddo + call conditionally_set_var(scm_input%input_smoiseq(:), physics%Sfcprop%smoiseq(i,:), "smoiseq", .false., missing_var(7)) + + ! do lsoil = -2, 4 + ! physics%Sfcprop%zsnsoxy(i,lsoil) = scm_input%input_zsnsoxy(lsoil) + ! enddo + call conditionally_set_var(scm_input%input_zsnsoxy(:), physics%Sfcprop%zsnsoxy(i,:), "zsnzoxy", .false., missing_var(8)) + + !write out warning if missing data for non-required variables + n = 8 + if ( i==1 .and. ANY( missing_var(1:n) ) ) then + cold_start_noahmp = .true. + write(0,'(a)') "INPUT CHECK: Some missing input data was found related to surface variables for NoahMP LSM. Due to this, a cold-start algorithm to initialize variables will be used." + write(0,'(a)') "Check gmtb_scm_type_defs.F90/physics_set to see the names of variables that are missing, corresponding to the following indices:" + do j=1, n + if (missing_var(j)) write(0,'(a,i0)') "variable index ",j + end do + end if + missing_var = .false. + endif + else if (scm_state%model_ics .and. physics%Model%lsm == physics%Model%lsm_ruc) then + !--- 3D variables + ! do lsoil = 1,Model%lsoil_lsm + ! physics%Sfcprop%tslb(i,lsoil) = scm_input%input_tslb(lsoil) + ! physics%Sfcprop%smois(i,lsoil) = scm_input%input_smois(lsoil) + ! physics%Sfcprop%sh2o(i,lsoil) = scm_input%input_sh2o(lsoil) + ! physics%Sfcprop%keepsmfr(i,lsoil) = scm_input%input_smfr(lsoil) + ! physics%Sfcprop%flag_frsoil(i,lsoil) = scm_input%input_flfr(lsoil) + ! enddo + call conditionally_set_var(scm_input%input_tslb(:), physics%Sfcprop%tslb(i,:), "tslb", .false., missing_var(1)) + call conditionally_set_var(scm_input%input_smois(:), physics%Sfcprop%smois(i,:), "smois", .false., missing_var(2)) + call conditionally_set_var(scm_input%input_sh2o(:), physics%Sfcprop%sh2o(i,:), "sh2o", .false., missing_var(3)) + call conditionally_set_var(scm_input%input_smfr(:), physics%Sfcprop%keepsmfr(i,:), "smfr", .false., missing_var(4)) + call conditionally_set_var(scm_input%input_flfr(:), physics%Sfcprop%flag_frsoil(i,:), "flfr", .false., missing_var(5)) + + !write out warning if missing data for non-required variables + n = 5 + if ( i==1 .and. ANY( missing_var(1:n) ) ) then + write(0,'(a)') "INPUT CHECK: Some missing input data was found related to (potentially non-required) surface variables for RUC LSM. This may lead to crashes or other strange behavior." + write(0,'(a)') "Check gmtb_scm_type_defs.F90/physics_set to see the names of variables that are missing, corresponding to the following indices:" + do j=1, n + if (missing_var(j)) write(0,'(a,i0)') "variable index ",j + end do + end if + missing_var = .false. + end if - end subroutine physics_associate + if (scm_state%model_ics) then + !check for nonmissing values + ! do k = 1,Model%kice + ! physics%Sfcprop%tiice(i,k) = scm_input%input_tiice(k) !--- internal ice temp + ! enddo + call conditionally_set_var(scm_input%input_tiice(:), physics%Sfcprop%tiice(i,:), "tiice", .false., missing_var(1)) + if (missing_var(1)) then + write(0,'(a)') "INPUT CHECK: Some missing input data was found related to the internal sea ice temperature. These will be set from the internal soil temperature variable (stc)." + end if + end if + + if (scm_input%input_tsfcl < real_zero) then + physics%Sfcprop%tsfcl(i) = physics%Sfcprop%tsfco(i) !--- compute tsfcl from existing variables + end if + + if (scm_input%input_zorll < real_zero) then + physics%Sfcprop%zorll(i) = physics%Sfcprop%zorlo(i) !--- compute zorll from existing variables + end if + + if (scm_input%input_zorli < real_zero) then + physics%Sfcprop%zorli(i) = physics%Sfcprop%zorlo(i) !--- compute zorli from existing variables + end if + + if (scm_input%input_zorlw < real_zero) then + physics%Sfcprop%zorlw(i) = physics%Sfcprop%zorlo(i) !--- compute zorlw from existing variables + end if + + if(physics%Model%frac_grid .and. scm_state%model_ics) then ! 3-way composite + if( physics%Model%phour < 1.e-7) physics%Sfcprop%tsfco(i) = max(con_tice, physics%Sfcprop%tsfco(i)) + tem1 = real_one - physics%Sfcprop%landfrac(i) + tem = tem1 * physics%Sfcprop%fice(i) ! tem = ice fraction wrt whole cell + physics%Sfcprop%zorl(i) = physics%Sfcprop%zorll(i) * physics%Sfcprop%landfrac(i) & + + physics%Sfcprop%zorli(i) * tem & + + physics%Sfcprop%zorlo(i) * (tem1-tem) + + physics%Sfcprop%tsfc(i) = physics%Sfcprop%tsfcl(i) * physics%Sfcprop%landfrac(i) & + + physics%Sfcprop%tisfc(i) * tem & + + physics%Sfcprop%tsfco(i) * (tem1-tem) + else + !--- specify tsfcl/zorll/zorli from existing variable tsfco/zorlo + ! physics%Sfcprop%tsfcl(i) = physics%Sfcprop%tsfco(i) + ! physics%Sfcprop%zorll(i) = physics%Sfcprop%zorlo(i) + ! physics%Sfcprop%zorli(i) = physics%Sfcprop%zorlo(i) + ! physics%Sfcprop%zorl(i) = physics%Sfcprop%zorlo(i) + ! physics%Sfcprop%tsfc(i) = physics%Sfcprop%tsfco(i) + if (physics%Sfcprop%slmsk(i) == 1) then + physics%Sfcprop%zorl(i) = physics%Sfcprop%zorll(i) + physics%Sfcprop%tsfc(i) = physics%Sfcprop%tsfcl(i) + else + tem = real_one - physics%Sfcprop%fice(i) + physics%Sfcprop%zorl(i) = physics%Sfcprop%zorli(i) * physics%Sfcprop%fice(i) & + + physics%Sfcprop%zorlo(i) * tem + + physics%Sfcprop%tsfc(i) = physics%Sfcprop%tisfc(i) * physics%Sfcprop%fice(i) & + + physics%Sfcprop%tsfco(i) * tem + endif + endif ! if (Model%frac_grid) + + if (scm_state%model_ics .and. MAXVAL(scm_input%input_tiice) < real_zero) then + physics%Sfcprop%tiice(i,1) = physics%Sfcprop%stc(i,1) !--- initialize internal ice temp from soil temp at layer 1 + physics%Sfcprop%tiice(i,2) = physics%Sfcprop%stc(i,2) !--- initialize internal ice temp from soil temp at layer 2 + end if + + if (cold_start_noahmp) then + !NoahMP cold start algorithm from FV3/io/FV3GFS_io + physics%Sfcprop%tvxy(i) = missing_value + physics%Sfcprop%tgxy(i) = missing_value + physics%Sfcprop%tahxy(i) = missing_value + physics%Sfcprop%canicexy(i) = missing_value + physics%Sfcprop%canliqxy(i) = missing_value + physics%Sfcprop%eahxy(i) = missing_value + physics%Sfcprop%cmxy(i) = missing_value + physics%Sfcprop%chxy(i) = missing_value + physics%Sfcprop%fwetxy(i) = missing_value + physics%Sfcprop%sneqvoxy(i) = missing_value + physics%Sfcprop%alboldxy(i) = missing_value + physics%Sfcprop%qsnowxy(i) = missing_value + physics%Sfcprop%wslakexy = missing_value + physics%Sfcprop%taussxy = missing_value + physics%Sfcprop%waxy(i) = missing_value + physics%Sfcprop%wtxy(i) = missing_value + physics%Sfcprop%zwtxy(i) = missing_value + physics%Sfcprop%xlaixy(i) = missing_value + physics%Sfcprop%xsaixy(i) = missing_value + + physics%Sfcprop%lfmassxy(i) = missing_value + physics%Sfcprop%stmassxy(i) = missing_value + physics%Sfcprop%rtmassxy(i) = missing_value + physics%Sfcprop%woodxy(i) = missing_value + physics%Sfcprop%stblcpxy(i) = missing_value + physics%Sfcprop%fastcpxy(i) = missing_value + physics%Sfcprop%smcwtdxy(i) = missing_value + physics%Sfcprop%deeprechxy(i) = missing_value + physics%Sfcprop%rechxy(i) = missing_value + + physics%Sfcprop%snowxy (i) = missing_value + physics%Sfcprop%snicexy(i, physics%Model%lsnow_lsm_lbound:0) = missing_value + physics%Sfcprop%snliqxy(i, physics%Model%lsnow_lsm_lbound:0) = missing_value + physics%Sfcprop%tsnoxy (i, physics%Model%lsnow_lsm_lbound:0) = missing_value + physics%Sfcprop%smoiseq(i, 1:physics%Model%lsoil_lsm) = missing_value + physics%Sfcprop%zsnsoxy(i, physics%Model%lsnow_lsm_lbound:physics%Model%lsoil_lsm) = missing_value + + if (physics%Sfcprop%landfrac(i) >= drythresh) then + + physics%Sfcprop%tvxy(i) = physics%Sfcprop%tsfcl(i) + physics%Sfcprop%tgxy(i) = physics%Sfcprop%tsfcl(i) + physics%Sfcprop%tahxy(i) = physics%Sfcprop%tsfcl(i) + + if (physics%Sfcprop%snowd(i) > 0.01 .and. physics%Sfcprop%tsfcl(i) > 273.15 ) then + physics%Sfcprop%tvxy(i) = 273.15 + physics%Sfcprop%tgxy(i) = 273.15 + physics%Sfcprop%tahxy(i) = 273.15 + end if + + physics%Sfcprop%canicexy(i) = 0.0 + physics%Sfcprop%canliqxy(i) = physics%Sfcprop%canopy(i) + + physics%Sfcprop%eahxy(i) = 2000.0 + +! eahxy = psfc*qv/(0.622+qv); qv is mixing ratio, converted from sepcific +! humidity specific humidity /(1.0 - specific humidity) + + physics%Sfcprop%cmxy(i) = 0.0 + physics%Sfcprop%chxy(i) = 0.0 + physics%Sfcprop%fwetxy(i) = 0.0 + physics%Sfcprop%sneqvoxy(i) = physics%Sfcprop%weasd(i) ! mm + physics%Sfcprop%alboldxy(i) = 0.65 + physics%Sfcprop%qsnowxy(i) = 0.0 + +! if (Sfcprop(nb)%srflag(ix) > 0.001) Sfcprop(nb)%qsnowxy(ix) = Sfcprop(nb)%tprcp(ix)/Model%dtp +! already set to 0.0 + physics%Sfcprop%wslakexy(i) = 0.0 + physics%Sfcprop%taussxy(i) = 0.0 + + + physics%Sfcprop%waxy(i) = 4900.0 + physics%Sfcprop%wtxy(i) = physics%Sfcprop%waxy(i) + physics%Sfcprop%zwtxy(i) = (25.0 + 2.0) - physics%Sfcprop%waxy(i) / 1000.0 /0.2 +! + vegtyp = physics%Sfcprop%vtype(i) + if (vegtyp == 0) vegtyp = 7 + imn = physics%Model%idate(2) + + if ((vegtyp == isbarren_table) .or. (vegtyp == isice_table) .or. (vegtyp == isurban_table) .or. (vegtyp == iswater_table)) then + + physics%Sfcprop%xlaixy(i) = 0.0 + physics%Sfcprop%xsaixy(i) = 0.0 + + physics%Sfcprop%lfmassxy(i) = 0.0 + physics%Sfcprop%stmassxy(i) = 0.0 + physics%Sfcprop%rtmassxy(i) = 0.0 + + physics%Sfcprop%woodxy (i) = 0.0 + physics%Sfcprop%stblcpxy (i) = 0.0 + physics%Sfcprop%fastcpxy (i) = 0.0 + + else + +! print *, 'vegtyp', vegtyp +! print *, 'imn', imn +! print *, 'xlaixy', Sfcprop(nb)%xlaixy(ix) + + physics%Sfcprop%xlaixy(i) = max(laim_table(vegtyp, imn),0.05) +! Sfcprop(nb)%xsaixy(ix) = max(saim_table(vegtyp, imn),0.05) + physics%Sfcprop%xsaixy(i) = max(physics%Sfcprop%xlaixy(i)*0.1,0.05) + + masslai = 1000.0 / max(sla_table(vegtyp),1.0) + physics%Sfcprop%lfmassxy(i) = physics%Sfcprop%xlaixy(i)*masslai + masssai = 1000.0 / 3.0 + physics%Sfcprop%stmassxy(i) = physics%Sfcprop%xsaixy(i)* masssai + + physics%Sfcprop%rtmassxy(i) = 500.0 + + physics%Sfcprop%woodxy (i) = 500.0 + physics%Sfcprop%stblcpxy(i) = 1000.0 + physics%Sfcprop%fastcpxy(i) = 1000.0 + + endif ! non urban ... + + if ( vegtyp == isice_table ) then + do lsoil = 1,physics%Model%lsoil + physics%Sfcprop%stc(i,lsoil) = min(physics%Sfcprop%stc(i,lsoil),min(physics%Sfcprop%tg3(i),263.15)) + physics%Sfcprop%smc(i,lsoil) = 1 + physics%Sfcprop%slc(i,lsoil) = 0 + enddo + endif + + snd = physics%Sfcprop%snowd(i)/1000.0 ! go to m from snwdph + + if (physics%Sfcprop%weasd(i) /= 0.0 .and. snd == 0.0 ) then + snd = physics%Sfcprop%weasd(i)/1000.0 + endif + + if (vegtyp == 15) then ! land ice in MODIS/IGBP + if ( physics%Sfcprop%weasd(i) < 0.1) then + physics%Sfcprop%weasd(i) = 0.1 + snd = 0.01 + endif + endif + + if (snd < 0.025 ) then + physics%Sfcprop%snowxy(i) = 0.0 + dzsno(-2:0) = 0.0 + elseif (snd >= 0.025 .and. snd <= 0.05 ) then + physics%Sfcprop%snowxy(i) = -1.0 + dzsno(0) = snd + elseif (snd > 0.05 .and. snd <= 0.10 ) then + physics%Sfcprop%snowxy(i) = -2.0 + dzsno(-1) = 0.5*snd + dzsno(0) = 0.5*snd + elseif (snd > 0.10 .and. snd <= 0.25 ) then + physics%Sfcprop%snowxy(i) = -2.0 + dzsno(-1) = 0.05 + dzsno(0) = snd - 0.05 + elseif (snd > 0.25 .and. snd <= 0.45 ) then + physics%Sfcprop%snowxy(i) = -3.0 + dzsno(-2) = 0.05 + dzsno(-1) = 0.5*(snd-0.05) + dzsno(0) = 0.5*(snd-0.05) + elseif (snd > 0.45) then + physics%Sfcprop%snowxy(i) = -3.0 + dzsno(-2) = 0.05 + dzsno(-1) = 0.20 + dzsno(0) = snd - 0.05 - 0.20 + else + write(0,'(a)') 'problem with the logic assigning snow layers in gmtb_scm_type_defs.F90/physics_set' + STOP + endif + +! Now we have the snowxy field +! snice + snliq + tsno allocation and compute them from what we have + +! + physics%Sfcprop%tsnoxy(i,physics%Model%lsnow_lsm_lbound:0) = 0.0 + physics%Sfcprop%snicexy(i,physics%Model%lsnow_lsm_lbound:0) = 0.0 + physics%Sfcprop%snliqxy(i,physics%Model%lsnow_lsm_lbound:0) = 0.0 + physics%Sfcprop%zsnsoxy(i,physics%Model%lsnow_lsm_lbound:physics%Model%lsoil_lsm) = 0.0 + + isnow = nint(physics%Sfcprop%snowxy(i))+1 ! snowxy <=0.0, dzsno >= 0.0 + + do ns = isnow , 0 + physics%Sfcprop%tsnoxy(i,ns) = physics%Sfcprop%tgxy(i) + physics%Sfcprop%snliqxy(i,ns) = 0.0 + physics%Sfcprop%snicexy(i,ns) = 1.00 * dzsno(ns) * physics%Sfcprop%weasd(i)/snd + enddo +! +!zsnsoxy, all negative ? +! + do ns = isnow, 0 + dzsnso(ns) = -dzsno(ns) + enddo + + do ns = 1 , physics%Model%lsoil_lsm + dzsnso(ns) = -dzs(ns) + enddo +! +! Assign to zsnsoxy +! + physics%Sfcprop%zsnsoxy(i,isnow) = dzsnso(isnow) + do ns = isnow+1,physics%Model%lsoil_lsm + physics%Sfcprop%zsnsoxy(i,ns) = physics%Sfcprop%zsnsoxy(i,ns-1) + dzsnso(ns) + enddo + +! +! smoiseq +! Init water table related quantities here +! + soiltyp = physics%Sfcprop%stype(i) + + if (soiltyp /= 0) then + bexp = bexp_table(soiltyp) + smcmax = smcmax_table(soiltyp) + smcwlt = smcwlt_table(soiltyp) + dwsat = dwsat_table(soiltyp) + dksat = dksat_table(soiltyp) + psisat = -psisat_table(soiltyp) + endif + + if (vegtyp == isurban_table) then + smcmax = 0.45 + smcwlt = 0.40 + endif + + if ((bexp > 0.0) .and. (smcmax > 0.0) .and. (-psisat > 0.0 )) then + do ns = 1, physics%Model%lsoil + if ( ns == 1 )then + ddz = -zsoil(ns+1) * 0.5 + elseif ( ns < physics%Model%lsoil ) then + ddz = ( zsoil(ns-1) - zsoil(ns+1) ) * 0.5 + else + ddz = zsoil(ns-1) - zsoil(ns) + endif +! +! Use newton-raphson method to find eq soil moisture +! + expon = bexp + 1. + aa = dwsat / ddz + bb = dksat / smcmax ** expon + + smc = 0.5 * smcmax + + do iter = 1, 100 + func = (smc - smcmax) * aa + bb * smc ** expon + dfunc = aa + bb * expon * smc ** bexp + dx = func / dfunc + smc = smc - dx + if ( abs (dx) < 1.e-6) exit + enddo ! iteration + physics%Sfcprop%smoiseq(i,ns) = min(max(smc,1.e-4),smcmax*0.99) + enddo ! ddz soil layer + else ! bexp <= 0.0 + physics%Sfcprop%smoiseq(i,1:physics%Model%lsoil) = smcmax + endif ! end the bexp condition +! + physics%Sfcprop%smcwtdxy(i) = smcmax + physics%Sfcprop%deeprechxy(i) = 0.0 + physics%Sfcprop%rechxy(i) = 0.0 + + endif !end if slmsk>0.01 (land only) + + end if !noahmp_cold_start + + end do + + end subroutine physics_set function get_tracer_index (tracer_names, name) diff --git a/scm/src/gmtb_scm_utils.F90 b/scm/src/gmtb_scm_utils.F90 index 3bc40b1b..442237c7 100644 --- a/scm/src/gmtb_scm_utils.F90 +++ b/scm/src/gmtb_scm_utils.F90 @@ -95,3 +95,199 @@ end subroutine w_to_omega !> @} !> @} end module gmtb_scm_utils + +module NetCDF_read + use gmtb_scm_kinds, only : sp, dp, qp + use netcdf + + implicit none + + real(kind=dp) :: missing_value = -9999.0 + + interface NetCDF_read_var + module procedure NetCDF_read_var_0d_int + module procedure NetCDF_read_var_0d + module procedure NetCDF_read_var_1d + module procedure NetCDF_read_var_2d + end interface + + contains + + subroutine NetCDF_read_var_0d_int(ncid, var_name, req, var_data) + + integer, intent(in) :: ncid + character (*), intent(in) :: var_name + logical, intent(in) :: req + integer, intent(out) :: var_data + + integer :: varID, ierr + + if (req) then + call check(NF90_INQ_VARID(ncid,var_name,varID)) + call check(NF90_GET_VAR(ncid,varID,var_data)) + else + ierr = NF90_INQ_VARID(ncid,var_name,varID) + if (ierr /= NF90_NOERR) then + var_data = missing_value + else + call check(NF90_GET_VAR(ncid,varID,var_data)) + end if + end if + + end subroutine NetCDF_read_var_0d_int + + subroutine NetCDF_read_var_0d(ncid, var_name, req, var_data) + + integer, intent(in) :: ncid + character (*), intent(in) :: var_name + logical, intent(in) :: req + real(kind=dp), intent(out) :: var_data + + integer :: varID, ierr + + if (req) then + call check(NF90_INQ_VARID(ncid,var_name,varID)) + call check(NF90_GET_VAR(ncid,varID,var_data)) + else + ierr = NF90_INQ_VARID(ncid,var_name,varID) + if (ierr /= NF90_NOERR) then + var_data = missing_value + else + call check(NF90_GET_VAR(ncid,varID,var_data)) + end if + end if + + end subroutine NetCDF_read_var_0d + + subroutine NetCDF_read_var_1d(ncid, var_name, req, var_data) + + integer, intent(in) :: ncid + character (*), intent(in) :: var_name + logical, intent(in) :: req + real(kind=dp), dimension(:), intent(out) :: var_data + + integer :: varID, ierr + + if (req) then + call check(NF90_INQ_VARID(ncid,var_name,varID)) + call check(NF90_GET_VAR(ncid,varID,var_data)) + else + ierr = NF90_INQ_VARID(ncid,var_name,varID) + if (ierr /= NF90_NOERR) then + var_data = missing_value + else + call check(NF90_GET_VAR(ncid,varID,var_data)) + end if + end if + + end subroutine NetCDF_read_var_1d + + subroutine NetCDF_read_var_2d(ncid, var_name, req, var_data) + + integer, intent(in) :: ncid + character (*), intent(in) :: var_name + logical, intent(in) :: req + real(kind=dp), dimension(:,:), intent(out) :: var_data + + integer :: varID, ierr + + if (req) then + call check(NF90_INQ_VARID(ncid,var_name,varID)) + call check(NF90_GET_VAR(ncid,varID,var_data)) + else + ierr = NF90_INQ_VARID(ncid,var_name,varID) + if (ierr /= NF90_NOERR) then + var_data = missing_value + else + call check(NF90_GET_VAR(ncid,varID,var_data)) + end if + end if + + end subroutine NetCDF_read_var_2d + + !Generic subroutine to check for netCDF I/O errors + subroutine check(status) + integer, intent ( in) :: status + + if(status /= NF90_NOERR) then + print *, trim(nf90_strerror(status)) + stop "stopped" + end if + end subroutine check + +end module NetCDF_read + +module data_qc + use gmtb_scm_kinds, only : sp, dp, qp + use NetCDF_read, only: missing_value + + implicit none + + interface check_missing + module procedure check_missing_0d + module procedure check_missing_1d + end interface + + interface conditionally_set_var + module procedure conditionally_set_var_0d + module procedure conditionally_set_var_1d + end interface + + contains + + subroutine check_missing_0d(var, missing) + real(kind=dp), intent(in) :: var + logical, intent(out) :: missing + + missing = .false. + if (var == missing_value) missing = .true. + end subroutine check_missing_0d + + subroutine check_missing_1d(var, missing) + real(kind=dp), dimension(:), intent(in) :: var + logical, intent(out) :: missing + + missing = .false. + if ( ANY(var == missing_value)) missing = .true. + end subroutine check_missing_1d + + subroutine conditionally_set_var_0d(input, set_var, input_name, req, missing) + real(kind=dp), intent(in) :: input + real(kind=dp), intent(inout) :: set_var + character (*), intent(in) :: input_name + logical, intent(in) :: req + logical, intent(out) :: missing + + call check_missing(input, missing) + + if (.not. missing) then + set_var = input + else + if (req) then + write(0,'(a,i0,a)') "The variable '" // input_name // "' in the case data file had missing data, but it is required for the given physics configuration. Stopping..." + STOP + end if + end if + + end subroutine conditionally_set_var_0d + + subroutine conditionally_set_var_1d(input, set_var, input_name, req, missing) + real(kind=dp), dimension(:), intent(in) :: input + real(kind=dp), dimension(:), intent(inout) :: set_var + character (*), intent(in) :: input_name + logical, intent(in) :: req + logical, intent(out) :: missing + + call check_missing(input, missing) + + if (.not. missing) then + set_var = input + else + if (req) then + write(0,'(a,i0,a)') "The variable '" // input_name // "' in the case data file had missing data, but it is required for the given physics configuration. Stopping..." + STOP + end if + end if + + end subroutine conditionally_set_var_1d +end module data_qc