Skip to content

Commit

Permalink
Merge pull request #80 from GreyEvenson-NOAA/noah_om_grid
Browse files Browse the repository at this point in the history
Add capacity to read-in gridded soils data from NetCDF
  • Loading branch information
GreyREvenson authored Aug 22, 2023
2 parents 3a2e191 + 7e1dfec commit efbde7d
Show file tree
Hide file tree
Showing 7 changed files with 264 additions and 129 deletions.
Binary file renamed data/netcdf_dummy.nc → data/netcdf_landuse.nc
Binary file not shown.
Binary file added data/netcdf_soils.nc
Binary file not shown.
4 changes: 2 additions & 2 deletions run/namelist.input
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,8 @@
startdate = "199801010630" ! UTC time start of simulation (YYYYMMDDhhmm)
enddate = "199901010630" ! UTC time end of simulation (YYYYMMDDhhmm)
forcing_filename = "../data/bondville.dat" ! file containing forcing data
vegtyp_filename = "../data/netcdf_dummy.nc"
landuse_filename = "../data/netcdf_landuse.nc"
soils_filename = "../data/netcdf_soils.nc"
output_filename = "../data/output.nc"
/

Expand Down Expand Up @@ -47,7 +48,6 @@
/

&structure
isltyp = 1 ! soil texture class
nsoil = 4 ! number of soil levels
nsnow = 3 ! number of snow levels
nveg = 20 ! number of vegetation types
Expand Down
2 changes: 1 addition & 1 deletion src/DomainGridType.f90
Original file line number Diff line number Diff line change
Expand Up @@ -141,7 +141,7 @@ subroutine InitTransfer(this,namelist,gridinfo)
this%ZREF = namelist%ZREF
this%vegtyp(:,:) = gridinfo%vegtyp(:,:)
this%croptype(:,:) = namelist%croptype
this%isltyp(:,:) = namelist%isltyp
this%isltyp(:,:) = gridinfo%isltyp(:,:)
this%IST(:,:) = namelist%sfctyp
this%soilcolor(:,:) = namelist%soilcolor
this%start_datetime = date_to_unix(namelist%startdate) ! returns seconds-since-1970-01-01
Expand Down
240 changes: 162 additions & 78 deletions src/GridInfoType.f90

Large diffs are not rendered by default.

20 changes: 10 additions & 10 deletions src/NamelistRead.f90
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,8 @@ module NamelistRead
character(len=12) :: startdate ! UTC start datetime of the model run ( YYYYMMDDHHmm )
character(len=12) :: enddate ! UTC end datetime of the model run ( YYYYMMDDHHmm )
character(len=256) :: forcing_filename ! directory/name of the input/forcing file
character(len=256) :: vegtyp_filename ! directory/name of NetCDF file holding gridded vegtyp values
character(len=256) :: landuse_filename ! directory/name of NetCDF file holding gridded land use (vegtyp) values
character(len=256) :: soils_filename ! directory/name of NetCDF file holding gridded soils (isltyp) values
character(len=256) :: output_filename ! directory/name of the output file
character(len=256) :: parameter_dir ! name of the directory where parameter TBLs reside
character(len=256) :: noahowp_table ! name of noahowp parameter table
Expand All @@ -26,7 +27,6 @@ module NamelistRead
real :: ZREF ! measurement height for wind speed (m)
real :: rain_snow_thresh ! rain-snow temperature threshold (°C)

integer :: isltyp ! soil type
integer :: nsoil ! number of soil layers
integer :: nsnow ! number of snow layers
integer :: nveg ! number of vegetation types
Expand Down Expand Up @@ -92,7 +92,8 @@ subroutine ReadNamelist(this, namelist_file)
character(len=12) :: startdate
character(len=12) :: enddate
character(len=256) :: forcing_filename
character(len=256) :: vegtyp_filename
character(len=256) :: landuse_filename
character(len=256) :: soils_filename
character(len=256) :: output_filename
character(len=256) :: parameter_dir
character(len=256) :: soil_table
Expand All @@ -105,7 +106,6 @@ subroutine ReadNamelist(this, namelist_file)
real :: ZREF ! measurement height for wind speed (m)
real :: rain_snow_thresh

integer :: isltyp
integer :: nsoil
integer :: nsnow
integer :: nveg
Expand Down Expand Up @@ -146,7 +146,7 @@ subroutine ReadNamelist(this, namelist_file)
!--------------------------- !
! define namelist groups !
!--------------------------- !
namelist / timing / dt,startdate,enddate,forcing_filename,vegtyp_filename,output_filename
namelist / timing / dt,startdate,enddate,forcing_filename,landuse_filename,soils_filename,output_filename
namelist / parameters / parameter_dir, soil_table, general_table, noahowp_table,&
soil_class_name, veg_class_name
namelist / location / terrain_slope,azimuth
Expand All @@ -157,7 +157,7 @@ subroutine ReadNamelist(this, namelist_file)
crop_model_option,snowsoil_temp_time_option,soil_temp_boundary_option,&
supercooled_water_option,stomatal_resistance_option,&
evap_srfc_resistance_option,subsurface_option
namelist / structure / isltyp,nsoil,nsnow,nveg,croptype,sfctyp,soilcolor
namelist / structure / nsoil,nsnow,nveg,croptype,sfctyp,soilcolor
namelist / initial_values / zsoil,dzsnso,sice,sh2o,zwt

! missing values against which namelist options can be checked
Expand All @@ -177,7 +177,8 @@ subroutine ReadNamelist(this, namelist_file)
startdate = stringMissing
enddate = stringMissing
forcing_filename = stringMissing
vegtyp_filename = stringMissing
landuse_filename = stringMissing
soils_filename = stringMissing
output_filename = stringMissing
parameter_dir = stringMissing
soil_table = stringMissing
Expand All @@ -189,7 +190,6 @@ subroutine ReadNamelist(this, namelist_file)
ZREF = realMissing
rain_snow_thresh = realMissing

isltyp = integerMissing
nsoil = integerMissing
nveg = integerMissing
nsnow = integerMissing
Expand Down Expand Up @@ -297,7 +297,8 @@ subroutine ReadNamelist(this, namelist_file)
if(startdate /= stringMissing) then; this%startdate = startdate; else; write(*,'(A)') 'ERROR: required entry startdate not found in namelist'; stop; end if
if(enddate /= stringMissing) then; this%enddate = enddate; else; write(*,'(A)') 'ERROR: required entry enddate not found in namelist'; stop; end if
if(forcing_filename /= stringMissing) then; this%forcing_filename = forcing_filename; else; write(*,'(A)') 'ERROR: required entry forcing_filename not found in namelist'; stop; end if
if(vegtyp_filename /= stringMissing) then; this%vegtyp_filename = vegtyp_filename; else; write(*,'(A)') 'ERROR: required entry vegtyp_filename not found in namelist'; stop; end if
if(landuse_filename /= stringMissing) then; this%landuse_filename = landuse_filename; else; write(*,'(A)') 'ERROR: required entry landuse_filename not found in namelist'; stop; end if
if(soils_filename /= stringMissing) then; this%soils_filename = soils_filename; else; write(*,'(A)') 'ERROR: required entry soils_filename not found in namelist'; stop; end if
if(output_filename /= stringMissing) then; this%output_filename = output_filename; else; write(*,'(A)') 'ERROR: required entry output_filename not found in namelist'; stop; end if
if(parameter_dir /= stringMissing) then; this%parameter_dir = parameter_dir; else; write(*,'(A)') 'ERROR: required entry parameter_dir not found in namelist'; stop; end if
if(soil_table /= stringMissing) then; this%soil_table = soil_table; else; write(*,'(A)') 'ERROR: required entry soil_table not found in namelist'; stop; end if
Expand All @@ -311,7 +312,6 @@ subroutine ReadNamelist(this, namelist_file)
if(zref /= realMissing) then; this%ZREF = ZREF; else; write(*,'(A)') 'ERROR: required entry ZREF not found in namelist'; stop; end if
if(rain_snow_thresh /= realMissing) then; this%rain_snow_thresh = rain_snow_thresh; else; write(*,'(A)') 'ERROR: required entry rain_snow_thresh not found in namelist'; stop; end if

if(isltyp /= integerMissing) then; this%isltyp = isltyp; else; write(*,'(A)') 'ERROR: required entry isltyp not found in namelist'; stop; end if
if(nsoil /= integerMissing) then; this%nsoil = nsoil; else; write(*,'(A)') 'ERROR: required entry nsoil not found in namelist'; stop; end if
if(nsnow /= integerMissing) then; this%nsnow = nsnow; else; write(*,'(A)') 'ERROR: required entry nsnow not found in namelist'; stop; end if
if(nveg /= integerMissing) then; this%nveg = nveg; else; write(*,'(A)') 'ERROR: required entry nveg not found in namelist'; stop; end if
Expand Down
127 changes: 89 additions & 38 deletions test/analysis/create_dummy_netcdf.py
Original file line number Diff line number Diff line change
@@ -1,51 +1,102 @@
import netCDF4 as nc
import numpy as np
import random
import random, os

# Define dataset contents
n_y = 3
n_x = 2
dx = 30. # arc-seconds (~1 km)
dy = 30.
vegtyp = 1
isltyp = 1
lat = 40.01
lon = -88.37
cnv_arcsec_to_degrees = 0.0002777777777777778
loc = '../../data'

# Create and open dataset
loc = '../../data/netcdf_dummy.nc'
dum = nc.Dataset(loc, 'w', format='NETCDF4')

# Define dimensions
dim_lat = dum.createDimension('Latitude', n_y)
dim_lon = dum.createDimension('Longitude', n_x)

# Create variables
var_lon = dum.createVariable('Longitude', 'f4', 'Longitude')
var_lat = dum.createVariable('Latitude', 'f4', 'Latitude')
var_vegtyp = dum.createVariable('vegtyp', 'i4', ('Latitude','Longitude'))

# Create attributes
var_lon.dx = dx
var_lat.dy = dy

# Fill variables
var_lat[:] = np.array([lat+(i*dy*cnv_arcsec_to_degrees) for i in range(n_y)])
var_lon[:] = np.array([lon+(i*dx*cnv_arcsec_to_degrees) for i in range(n_x)])
var_vegtyp[:,:] = np.full((n_y, n_x), vegtyp)
for x in range(n_x):
for y in range(n_y):
if x != 0 or y != 0:
var_vegtyp[y,x] = random.randint(1,20)

# Print some stuff
print('\n****************DIMENSIONS****************')
print(dim_lon)
print(dim_lat)
print('\n****************VARIABLES****************')
print(var_lon)
print(var_lat)
print(var_vegtyp)

# Close dataset
dum.close()
def create_dummy_soils():

# Create and open dataset
filename = os.path.join(loc,'netcdf_soils.nc')
dum = nc.Dataset(filename, 'w', format='NETCDF4')

# Define dimensions
dim_lat = dum.createDimension('Latitude', n_y)
dim_lon = dum.createDimension('Longitude', n_x)

# Create variables
var_lon = dum.createVariable('Longitude', 'f4', 'Longitude')
var_lat = dum.createVariable('Latitude', 'f4', 'Latitude')
var_isltyp = dum.createVariable('isltyp', 'i4', ('Latitude','Longitude'))

# Create attributes
var_lon.dx = dx
var_lat.dy = dy

# Fill variables
var_lat[:] = np.array([lat+(i*dy*cnv_arcsec_to_degrees) for i in range(n_y)])
var_lon[:] = np.array([lon+(i*dx*cnv_arcsec_to_degrees) for i in range(n_x)])
var_isltyp[:,:] = np.full((n_y, n_x), isltyp)
for x in range(n_x):
for y in range(n_y):
if x != 0 or y != 0:
var_isltyp[y,x] = random.randint(1,19)

# Print some stuff
print('\n****************DIMENSIONS****************')
print(dim_lon)
print(dim_lat)
print('\n****************VARIABLES****************')
print(var_lon)
print(var_lat)
print(var_isltyp)

# Close dataset
dum.close()

def create_dummy_vegtyp():

# Create and open dataset
filename = os.path.join(loc,'netcdf_landuse.nc')
dum = nc.Dataset(filename, 'w', format='NETCDF4')

# Define dimensions
dim_lat = dum.createDimension('Latitude', n_y)
dim_lon = dum.createDimension('Longitude', n_x)

# Create variables
var_lon = dum.createVariable('Longitude', 'f4', 'Longitude')
var_lat = dum.createVariable('Latitude', 'f4', 'Latitude')
var_vegtyp = dum.createVariable('vegtyp', 'i4', ('Latitude','Longitude'))

# Create attributes
var_lon.dx = dx
var_lat.dy = dy

# Fill variables
var_lat[:] = np.array([lat+(i*dy*cnv_arcsec_to_degrees) for i in range(n_y)])
var_lon[:] = np.array([lon+(i*dx*cnv_arcsec_to_degrees) for i in range(n_x)])
var_vegtyp[:,:] = np.full((n_y, n_x), vegtyp)
for x in range(n_x):
for y in range(n_y):
if x != 0 or y != 0:
var_vegtyp[y,x] = random.randint(1,20)

# Print some stuff
print('\n****************DIMENSIONS****************')
print(dim_lon)
print(dim_lat)
print('\n****************VARIABLES****************')
print(var_lon)
print(var_lat)
print(var_vegtyp)

# Close dataset
dum.close()

def main():
create_dummy_vegtyp()
create_dummy_soils()

if __name__ == '__main__':
main()

0 comments on commit efbde7d

Please sign in to comment.