When selecting single-point execution (SOIL_DATA_OPTION=2), how should I set SOILCL1, SOILCL2, SOILCL3, and SOILCL4?

[yxjdbly]

When I selected SOIL_DATA_OPTION=2 for single-point execution, I attempted to modify the bondville_xx.dat file (my data file) based on my understanding. The specific modifications are as follows:
!----------------------------------------------------------------------------------------------------------------------
!Location Information
!----------------------------------------------------------------------------------------------------------------------
&location
latitude = 42.4025 ! latitude [degrees]
longitude = 128.0958 ! longitude [degrees]
vegetation_category = 5 ! category based on land_cover_source
soil_category = 4 ! soil texture class
deep_soil_temperature = 276.537 ! deep soil temperature [K]
elevation = 731.0 ! elevation [m]
sea_ice = 0.0 ! sea ice fraction [-] (should be 0.0 for a land point)
maximum_vegetation_pct = 96.0 ! maximum annual vegetation percentage [0-100]
minimum_vegetation_pct = 1.0 ! minimum annual vegetation percentage [0-100]
land_mask = 1 ! land mask [-] (should be 1 for a land point)
nsoil = 4 ! number of soil layers
/
!----------------------------------------------------------------------------------------------------------------------
!Initial States
!----------------------------------------------------------------------------------------------------------------------
&initial
snow_depth = 0.00 ! snow depth [m]
snow_water_equivalent = 0.0 ! snow water equivalent [mm]
canopy_water = 0.00 ! canopy water [mm]
skin_temperature = 271.551 ! skin temperature [K]
soil_layer_thickness = 0.05, 0.15, 0.30, 0.50 ! layer thicknesses [m]
soil_layer_nodes = 0.25, 0.125, 0.35, 0.75 ! layer centroids from surface [m]
soil_temperature = 271.551, 272.019, 273.01, 273.01 ! layer soil temperature [K]
soil_moisture = 0.133, 0.19, 0.337, 0.337 ! layer volumetric total water content [m3/m3]
leaf_area_index = 0.2 ! leaf area index [m2/m2] used to initialize carbon budgets
soilcl1 = 0.5:1.0,0.0:0.3
In addition, I also modified create_point_data.f90. The specific modifications are as follows, but I'm not sure if these changes are correct or necessary:
use netcdf
implicit none

! Location Information

real :: latitude ! latitude [degrees]
real :: longitude ! longitude [degrees]
integer :: vegetation_category ! category based on land_cover_source
integer :: soil_category ! soil texture class
real :: deep_soil_temperature ! deep soil temperature [K]
real :: elevation ! elevation [m]
real :: sea_ice ! sea ice fraction [-] (should be 0.0 for a land point)
real :: maximum_vegetation_pct ! maximum annual vegetation percentage [0-100]
real :: minimum_vegetation_pct ! minimum annual vegetation percentage [0-100]
integer :: land_mask ! land mask [-] (should be 1 for a land point)
integer :: nsoil ! number of soil layers
! Initial States

real :: snow_depth ! snow depth [m]
real :: snow_water_equivalent ! snow water equivalent [m]
real :: canopy_water ! canopy water [mm]
real :: skin_temperature ! skin temperature [K]
real :: leaf_area_index ! leaf area index [m2/m2]
real, allocatable, dimension(:) :: soil_layer_thickness ! layer thicknesses [m]
real, allocatable, dimension(:) :: soil_layer_nodes ! layer centroids from surface [m]
real, allocatable, dimension(:) :: soil_temperature ! layer soil temperature [K]
real, allocatable, dimension(:) :: soil_moisture ! layer volumetric total water content [m3/m3]
real, allocatable, dimension(:) :: soil_liquid ! layer volumetric liquid content [m3/m3]
real, allocatable, dimension(:) :: soilcl1
real, allocatable, dimension(:) :: soilcl2
real, allocatable, dimension(:) :: soilcl3
real, allocatable, dimension(:) :: soilcl4

! Metadata

integer :: grid_id ! used for grid labeling
integer :: water_classification ! water type in land classification
integer :: urban_classification ! urban type in land classification
integer :: ice_classification ! snow/ice type in land classification
character*100 :: land_cover_source ! land classification

! Conversion

logical :: have_relative_humidity ! set to true if you need to convert relative to specific
real :: temperature_offset ! 273.15 to convert C to K
real :: temperature_scale ! can be used to convert F to C
real :: pressure_scale ! 100.0 to convert mb to Pa
real :: precipitation_scale ! 0.014111 to convert inches/30min to mm/s

! Text forcing variables

integer :: yyyy
integer :: mm
integer :: dd
integer :: hh
integer :: nn
real :: uwind_speed
real :: vwind_speed = 0.0
real :: temperature
real :: humidity
real :: pressure
real :: solar_radiation
real :: longwave_radiation
real :: precipitation

! NETCDF variables

integer :: ncid
integer :: iret
integer :: dim_id_time, dim_id_sn, dim_id_we, dim_id_soil
integer :: var_id_xlat , var_id_xlong , var_id_hgt , var_id_tmn , &
var_id_veg , var_id_soil , var_id_seaice , var_id_vegmax , &
var_id_vegmin , var_id_mask
integer :: var_id_depth , var_id_swe , var_id_canwat , var_id_tsk , &
var_id_thick , var_id_nodes , var_id_stemp , var_id_smois , &
var_id_sliq , var_id_lai
integer :: var_id_grid , var_id_water , var_id_urban , var_id_ice , &
var_id_source
integer :: var_id_force_temp , var_id_force_humid , var_id_force_pres , &
var_id_force_uwind , var_id_force_vwind , &
var_id_force_sw , var_id_force_lw , var_id_force_precip
! Misc variables

integer :: itime
character(len=27) :: filename
real :: e, svp ! for relative humidity conversion

namelist / location / latitude, longitude, vegetation_category, soil_category, &
deep_soil_temperature, elevation, sea_ice, maximum_vegetation_pct, &
minimum_vegetation_pct, land_mask, nsoil

namelist / initial / snow_depth, snow_water_equivalent, canopy_water, &
skin_temperature, soil_layer_thickness, soil_layer_nodes, &
soil_temperature, soil_moisture, soil_liquid, leaf_area_index

namelist / metadata / grid_id, water_classification, urban_classification, &
ice_classification, land_cover_source

namelist / conversion / have_relative_humidity, temperature_offset, &
temperature_scale, pressure_scale, precipitation_scale

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Start creation of hrldas setup file
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

open(30, file="bondville_2003-2005.dat", form="FORMATTED")
read(30, location)

allocate(soil_layer_thickness(nsoil))
allocate(soil_layer_nodes(nsoil))
allocate(soil_temperature(nsoil))
allocate(soil_moisture(nsoil))
allocate(soil_liquid(nsoil))

read(30, initial)
read(30, metadata)
read(30, conversion)
close(30)

! Create the NetCDF file mimicing WRF input files.

iret = nf90_create("hrldas_setup_single_point.nc", NF90_CLOBBER, ncid)

! Define dimensions in the file.

iret = nf90_def_dim(ncid, "west_east" , 1 , dim_id_we)
iret = nf90_def_dim(ncid, "south_north" , 1 , dim_id_sn)
iret = nf90_def_dim(ncid, "soil_layers_stag", nsoil , dim_id_soil)
iret = nf90_def_dim(ncid, "Time" , NF90_UNLIMITED, dim_id_time)

iret = nf90_def_var(ncid, "XLAT", NF90_FLOAT, (/dim_id_we,dim_id_sn,dim_id_time/), var_id_xlat)
iret = nf90_put_att(ncid, var_id_xlat, "description", "latitude")
iret = nf90_put_att(ncid, var_id_xlat, "units", "degrees_north")

iret = nf90_def_var(ncid, "XLONG", NF90_FLOAT, (/dim_id_we,dim_id_sn,dim_id_time/), var_id_xlong)
iret = nf90_put_att(ncid, var_id_xlong, "description", "longitude")
iret = nf90_put_att(ncid, var_id_xlong, "units", "degrees_east")

iret = nf90_def_var(ncid, "TMN", NF90_FLOAT, (/dim_id_we,dim_id_sn,dim_id_time/), var_id_tmn)
iret = nf90_put_att(ncid, var_id_tmn, "description", "soil temperature lower boundary")
iret = nf90_put_att(ncid, var_id_tmn, "units", "K")

iret = nf90_def_var(ncid, "HGT", NF90_FLOAT, (/dim_id_we,dim_id_sn,dim_id_time/), var_id_hgt)
iret = nf90_put_att(ncid, var_id_hgt, "description", "elevation")
iret = nf90_def_var(ncid, "SEAICE", NF90_FLOAT, (/dim_id_we,dim_id_sn,dim_id_time/), var_id_seaice)
iret = nf90_put_att(ncid, var_id_seaice, "description", "sea ice fraction")
iret = nf90_put_att(ncid, var_id_seaice, "units", "-")

iret = nf90_def_var(ncid, "SHDMAX", NF90_FLOAT, (/dim_id_we,dim_id_sn,dim_id_time/), var_id_vegmax)
iret = nf90_put_att(ncid, var_id_vegmax, "description", "maximum annual vegetation cover")
iret = nf90_put_att(ncid, var_id_vegmax, "units", "%")

iret = nf90_def_var(ncid, "SHDMIN", NF90_FLOAT, (/dim_id_we,dim_id_sn,dim_id_time/), var_id_vegmin)
iret = nf90_put_att(ncid, var_id_vegmin, "description", "minimum annual vegetation cover")
iret = nf90_put_att(ncid, var_id_vegmin, "units", "%")

iret = nf90_def_var(ncid, "LAI", NF90_FLOAT, (/dim_id_we,dim_id_sn,dim_id_time/), var_id_lai)
iret = nf90_put_att(ncid, var_id_lai, "description", "leaf area index")
iret = nf90_put_att(ncid, var_id_lai, "units", "m2/m2")

iret = nf90_def_var(ncid, "XLAND", NF90_INT, (/dim_id_we,dim_id_sn,dim_id_time/), var_id_mask)
iret = nf90_put_att(ncid, var_id_mask, "description", "land mask: 1=land/ice")
iret = nf90_put_att(ncid, var_id_mask, "units", "-")

iret = nf90_def_var(ncid, "IVGTYP", NF90_INT, (/dim_id_we,dim_id_sn,dim_id_time/), var_id_veg)
iret = nf90_put_att(ncid, var_id_veg, "description", "vegetation type")
iret = nf90_put_att(ncid, var_id_veg, "units", "-")

iret = nf90_def_var(ncid, "ISLTYP", NF90_INT, (/dim_id_we,dim_id_sn,dim_id_time/), var_id_soil)
iret = nf90_put_att(ncid, var_id_soil, "description", "soil texture type")
iret = nf90_put_att(ncid, var_id_soil, "units", "-")
iret = nf90_def_var(ncid, "SNOW", NF90_FLOAT, (/dim_id_we,dim_id_sn,dim_id_time/), var_id_swe)
iret = nf90_put_att(ncid, var_id_swe, "description", "snow water equivalent")
iret = nf90_put_att(ncid, var_id_swe, "units", "mm")

iret = nf90_def_var(ncid, "SNODEP", NF90_FLOAT, (/dim_id_we,dim_id_sn,dim_id_time/), var_id_depth)
iret = nf90_put_att(ncid, var_id_depth, "description", "snow depth")
iret = nf90_put_att(ncid, var_id_depth, "units", "m")

iret = nf90_def_var(ncid, "CANWAT", NF90_FLOAT, (/dim_id_we,dim_id_sn,dim_id_time/), var_id_canwat)
iret = nf90_put_att(ncid, var_id_canwat, "description", "canopy surface water")
iret = nf90_put_att(ncid, var_id_canwat, "units", "mm")

iret = nf90_def_var(ncid, "TSK", NF90_FLOAT, (/dim_id_we,dim_id_sn,dim_id_time/), var_id_tsk)
iret = nf90_put_att(ncid, var_id_tsk, "description", "surface skin temperature")
iret = nf90_put_att(ncid, var_id_tsk, "units", "K")

iret = nf90_def_var(ncid, "DZS", NF90_FLOAT, (/dim_id_soil,dim_id_time/) , var_id_thick)
iret = nf90_put_att(ncid, var_id_thick, "description", "soil layer thicknesses")
iret = nf90_put_att(ncid, var_id_thick, "units", "m")

iret = nf90_def_var(ncid, "ZS", NF90_FLOAT, (/dim_id_soil,dim_id_time/) , var_id_nodes)
iret = nf90_put_att(ncid, var_id_nodes, "description", "soil node depths")
iret = nf90_put_att(ncid, var_id_nodes, "units", "m")

iret = nf90_def_var(ncid, "TSLB", NF90_FLOAT, (/dim_id_we,dim_id_sn,dim_id_soil,dim_id_time/), var_id_stemp)
iret = nf90_put_att(ncid, var_id_stemp, "description", "soil temperature")
iret = nf90_put_att(ncid, var_id_stemp, "units", "K")

iret = nf90_def_var(ncid, "SMOIS", NF90_FLOAT, (/dim_id_we,dim_id_sn,dim_id_soil,dim_id_time/), var_id_smois)
iret = nf90_put_att(ncid, var_id_smois, "description", "volumetric soil moisture")
iret = nf90_put_att(ncid, var_id_smois, "units", "m3/m3")
iret = nf90_def_var(ncid, "SOILCL1", NF90_INT, (/dim_id_we,dim_id_sn,dim_id_time/), var_id_soilol1)
iret = nf90_put_att(ncid, var_id_soilcl1, "description", "soilcl1")
iret = nf90_put_att(ncid, var_id_soilol1, "units", "-")
! File metadata:

iret = nf90_put_att(ncid, NF90_GLOBAL, "GRID_ID" , grid_id)
iret = nf90_put_att(ncid, NF90_GLOBAL, "ISWATER" , water_classification)
iret = nf90_put_att(ncid, NF90_GLOBAL, "ISURBAN" , urban_classification)
iret = nf90_put_att(ncid, NF90_GLOBAL, "ISICE" , ice_classification)
iret = nf90_put_att(ncid, NF90_GLOBAL, "MMINLU" , land_cover_source)
iret = nf90_put_att(ncid, NF90_GLOBAL, "DX" , 0.0)
iret = nf90_put_att(ncid, NF90_GLOBAL, "DY" , 0.0)
iret = nf90_put_att(ncid, NF90_GLOBAL, "TRUELAT1" , 0.0)
iret = nf90_put_att(ncid, NF90_GLOBAL, "TRUELAT2" , 0.0)
iret = nf90_put_att(ncid, NF90_GLOBAL, "STAND_LON", 0.0)
iret = nf90_put_att(ncid, NF90_GLOBAL, "MAP_PROJ" , 0)
iret = nf90_put_att(ncid, NF90_GLOBAL, "TITLE" , "Created from HRLDAS create_point_data")

iret = nf90_enddef(ncid)

iret = nf90_put_var(ncid, var_id_xlat , latitude )
iret = nf90_put_var(ncid, var_id_xlong , longitude )
iret = nf90_put_var(ncid, var_id_tmn , deep_soil_temperature )
iret = nf90_put_var(ncid, var_id_hgt , elevation )
iret = nf90_put_var(ncid, var_id_seaice, sea_ice )
iret = nf90_put_var(ncid, var_id_vegmax, maximum_vegetation_pct)
iret = nf90_put_var(ncid, var_id_vegmin, minimum_vegetation_pct)
iret = nf90_put_var(ncid, var_id_lai , leaf_area_index )
iret = nf90_put_var(ncid, var_id_mask , land_mask )
iret = nf90_put_var(ncid, var_id_veg , vegetation_category )
iret = nf90_put_var(ncid, var_id_soil , soil_category )
iret = nf90_put_var(ncid, var_id_soilcl1 , soilcl1 )
These are my modifications, not present in the original file, added here for illustration purposes only
create_point_data.zip
error:forrtl: severe (19): invalid reference to variable in NAMELIST input, unit 30, file /home/xxxx/bondville_2003-2005.dat, line 30, position 25
Dear Professor，Please guide me and provide specific modification suggestions. Thank you

[tslin2]

Please see if this helps
#173 (comment)

[tslin2]

I guess you should use comma, not colon for soilcl1.
soilcl1 = 0.5:1.0,0.0:0.3

for each layer soil type should have only one value
SOILCL1
SOILCL2
SOILCL3
SOILCL4

[yxjdbly]

First of all, thank you for your answer. As a beginner, due to my misunderstanding, the input parameters of the data are somewhat biased, which has resulted in the program not running effectively. Following your answer, I still don't have a clear understanding. This is my problem. Can you further explain what specific actions I should take?
Question one: Should I input all parameters in the create_point_data.f90 file, and should I enter the corresponding running code (like "LAI")?
Here is my understanding, please help me confirm if it is correct

I need to confirm the input parameters here.(like real, allocatable, dimension(:) :: soilcl1)
This part I don't know how to modify, and I think var_id_xxx is a predefined variable, but I don't know how to determine the id for SOILCL1.

This part I'm unsure whether I should proceed with.

This part, I'm not sure if my understanding is correct. Based on these two figures, I should write a similar code to define SOILCL1


I've done a similar operation in this part, but I'm not sure if it's correct.

[tslin2]

Hi,
you can try borrowing the code logistics for other variables like soil_category or leaf_area_index to implement SOILCL1, SOILCL2, SOILCL3, SOILCL4, and you had to give each layer soil type variable a value for soil type either in the namelist or in this code directly.

if you do it in the code directly, for example, you give each variable a number to represent soil type for each layer of soil

real :: SOILCL1 = 1
real :: SOILCL2 = 2
real :: SOILCL3 = 2
real :: SOILCL4 = 1

1 means SAND and 2 means LOAMY SAND, you can change to the type you want based on the number you give

Then you will add each variable in the netCDF format, such as

integer ::        var_id_soilcl1

iret = nf90_def_var(ncid,    "SOILCL1", NF90_FLOAT, (/dim_id_we,dim_id_sn,dim_id_time/), var_id_soilcl1)
iret = nf90_put_att(ncid, var_id_soilcl1, "description", "soil texture type layer 1")
iret = nf90_put_att(ncid, var_id_soilcl1, "units", "-")

iret = nf90_put_var(ncid, var_id_soilcl1   , SOILCL1       )

below is the website that includes some code examples to learn how to read/write netCDF files and format in Fortran or other programming languages.
https://www.unidata.ucar.edu/software/netcdf/examples/programs/

[cenlinhe]

Note that the input capability for SOILCL1~4 is already there in HRLDAS/NoahMP if you use soil_option = 2:
https://github.com/NCAR/hrldas/blob/master/hrldas/IO_code/module_NoahMP_hrldas_driver.F#L231
https://github.com/NCAR/hrldas/blob/master/hrldas/IO_code/module_hrldas_netcdf_io.F#L1374
There is no need to implement the I/O code yourself. You just need to set the SOILCL1-4 in your hrldas_setup file or your bondville_xx.dat file