To generate the grid surface data file, we will load an existing
surface data file, locate the grid cell that is the closest to the point
of interest, abd extract the data. For some of the variables, such as
soil texture and soil organic carbon, there may be local observations.
If this is the case, we replace the gridded information with the local
observations.
Load the gridded surface data
Open the surface data file, and retrieve the grid dimensions.
# Open gridded surface data
if ("nc_gridin" %in% ls()){dummy=nc_close(nc_gridin);rm(nc_gridin)}
nc_gridin = nc_open(filename=input_surfgrid)
Retrieve global attributes from the input file, and add information
about this processing.
# Retrieve original attributes. We only make minor edits.
glob_att = ncatt_get(nc = nc_gridin, varid = 0)
# Define the code developer information (indirect way so the email is not visible).
developer_name = c( 111L, 103L, 110L, 111L, 76L, 32L, 115L, 111L, 99L, 114L, 97L, 77L)
developer_email = c( 118L, 111L, 103L, 46L, 108L, 98L, 108L, 64L, 111L, 103L, 110L, 111L
, 108L, 109L)
# Get current time for the log
create_time = paste0(as.character(now(tzone="UTC"))," UTC")
# Generate information about this file.
glob_new = list( title = paste0( "Domain file for ",site_desc,".")
, Version = dat_version
, Revision_Id = paste0( "$Id: make_fates_domain+surface.Rmd"
, " 0 ",create_time,Sys.info()["user"],"$"
)#end paste0
, History_Log = paste0( "created on: ", create_time)
, Logname = paste0( "code originally developed by "
, intToUtf8(rev(developer_name))
," <"
, intToUtf8(rev(developer_email))
,">"
)#end paste0
, Host = Sys.info()["nodename"]
, File_Author = paste0(author_name," <",author_email,">")
)#end list
# Update attribute list
glob_att = modifyList( glob_att, glob_new)
Load coordinates, and identify the nearest grid point.
# Retrieve longitude and latitude
in_lon = melt(ncvar_get(nc=nc_gridin,"LONGXY" )); names(in_lon) = c("x","y","lon" )
in_lat = melt(ncvar_get(nc=nc_gridin,"LATIXY" )); names(in_lat) = c("x","y","lat" )
in_veg = melt(ncvar_get(nc=nc_gridin,"PCT_NATVEG" )); names(in_veg) = c("x","y","pcveg" )
in_land = melt(ncvar_get(nc=nc_gridin,"LANDFRAC_PFT")); names(in_land) = c("x","y","fland" )
# Create a tibble with coordinates
in_coord = tibble(merge(merge(merge(in_lon,in_lat),in_veg),in_land))
# Define function that takes longitude/latitude and find distances.
rdist_wrap = function(lon,lat,lon0,lat0){
ans = c(fields::rdist.earth(x1=cbind(lon0,lat0),x2=cbind(lon,lat),miles=FALSE))
return(ans)
}#end function
# Find the distances.
in_coord = in_coord %>%
mutate( dist = rdist_wrap(lon=lon,lat=lat,lon0=clon,lat0=clat))
# Select the nearest neighbour
in_near = in_coord %>%
filter((pcveg > pcveg_min) & (fland > landfrac_min)) %>%
slice(which.min(dist))
in_xy = in_near %>% select(c(x,y)) %>% unlist(); names(in_xy) = NULL
Load all the dimensions from the input netCDF file, and create
dimensions for the output file that are essentially copies of the input.
The only exceptions are the longitude and latitude, which will be always
one (single site).
# List of dimension names
nc_dnames = names(nc_gridin$dim)
# Loop through all the dimensions
nc_dims = list()
for (d in seq_along(nc_dnames)){
# Retrieve dimensions
d_dname = nc_dnames[d]
nc_dimnow = nc_gridin$dim[[d_dname]]
# Fix dimensions and values for longitude and latitude.
if (d_dname %in% c("lsmlon","lsmlat")){
# Longitude or latitude, make data a single point
nc_dimnow = modifyList( x = nc_dimnow
, val = list( len = 1L, vals = sequence(1L))
)#end modifyList
# }else if (d_dname %in% "time"){
# # Time is not properly set in some nc maps.
# nc_dimnow = modifyList( x = nc_dimnow
# , val = list(vals=seq_along(nc_dimnow$vals))
# )#end modifyList
}#end if (nc_dimnow)
# Decide whether to append calendar information
if (nc_dimnow$unlim && (! is.null(nc_dimnow$calendar))){
nc_calendar = nc_dimnow$calendar
}else{
nc_calendar = NA_character_
}#end if (nc_dimnow$unlim && (! is.null(nc_dimnow$calendar)))
# Create dimensions. Distinguish unlimited and limited dimensions (if needed).
nc_dims[[d_dname]] = ncdim_def( name = d_dname
, units = nc_dimnow$units
, vals = nc_dimnow$vals
, unlim = nc_dimnow$unlim
, create_dimvar = nc_dimnow$create_dimvar
, calendar = nc_calendar
)#end ncdim_def
}#end for (d in seq_along(nc_dnames))
In case depth to bedrock in known, include it in the list of
variables. We also reset the number of active soil layers
nlevsoi and adjust the variables that depend upon
nlevsoi
# Check whether or not to overwrite the number of active layers
use_depth_geotiff = tolower(site_bedrock) %in% "gridded"
use_depth_constant = (! use_depth_geotiff) && (! is.na(site_bedrock))
if (use_depth_geotiff || use_depth_constant){
# Check how to define the depth to bedrock
if (use_depth_constant){
# Depth to bedrock is from user-defined constant
ud_bedrock = br_scale * site_bedrock
}else{
# Obtain depth to bedrock from the GeoTIFF file
rast_bedrock = rast(x=input_soildepth)
# Append the coordinates into a spatial vector
site_coord = data.frame(lon = site_lon, lat = site_lat)
site_coord = vect(x = site_coord, geom = c("lon","lat"), crs = "epsg:4326")
# Decide how to extract the sample
if (is.na(br_window)){
# Reproject coordinates to the same projection used by the raster
site_coord = project(x = site_coord, y = crs(rast_bedrock))
# No window, simply extract the nearest pixel
br_neighbour = terra::extract(x=rast_bedrock, y = site_coord)[,2L]
}else{
# Add a buffer so we extract multiple pixels
site_buffer = buffer(x=site_coord,width=br_window)
# Reproject the coordinates so they match the raster.
site_coord = project(x = site_coord , y = crs(rast_bedrock))
site_buffer = project(x = site_buffer, y = crs(rast_bedrock))
# Extract the pixels within the buffer
br_neighbour = terra::extract(x=rast_bedrock, y = site_buffer)[,2L]
}#end if (! is.na(br_window))
# Average the sampled points, and fix units
ud_bedrock = br_scale * mean(br_neighbour,na.rm = TRUE)
}#end if (is.na(br_direction))
# Find the layer associated with the depth to bedrock
nlevsoi = max(which(ud_bedrock >= ztop))
nlevsoi_orig = nc_dims$nlevsoi$len
# Overwrite dimensions for nlevsoi
nc_dims$nlevsoi = ncdim_def( name = nc_dims$nlevsoi$name
, units = nc_dims$nlevsoi$units
, vals = sequence(nlevsoi)
, unlim = nc_dims$nlevsoi$unlim
, create_dimvar = nc_dims$nlevsoi$create_dimvar
, calendar = nc_calendar
)#end ncdim_def
# Find the sequence of indices for updating variables that depend on the soil
idx_levsoi = pmin(nlevsoi_orig,sequence(nlevsoi))
# Append information for depth to bedrock
if (use_depth_geotiff){
map_soil_depth_file = basename(input_soildepth)
}else{
map_soil_depth_file = "site-specific information"
}#end if (use_soilgrids)
glob_att = modifyList( x = glob_att
, val = list( map_soil_depth_file = map_soil_depth_file )
)#end modifyList
}else{
# Save the number of active soil layers in case we need for soil interpolation
nlevsoi_orig = nc_dims$nlevsoi$len
idx_levsoi = sequence(nlevsoi)
}#end if (use_geotiff || use_constant)
Define variables that should feature in the output. Every variable in
the input should be in the output too.
# List of variable names
nc_vnames = names(nc_gridin$var)
# This will contain variables that were actually used. We will create dummy variables for
# those not used, to ensure they exist in the output.
nc_dimused = character(0)
# Loop through all the variables, retrieve data and settings
nc_vars = list()
nc_values = list()
for (v in seq_along(nc_vnames)){
# Retrieve variable information
v_vname = nc_vnames[v]
v_vinfo = nc_gridin$var [[v_vname]]
v_vlong = v_vinfo$longname
v_vunit = v_vinfo$unit
v_vprec = if( v_vinfo$prec %in% "int"){"integer"}else{v_vinfo$prec}
v_vmiss = if( v_vinfo$prec %in% c("float","double")){v_vinfo$missval}else{NULL}
# Switch the longitude and latitude dimensions with point
if (is.null(v_vinfo$varsize)){
# Dummy dimensions
v_count = NA
v_start = NA
v_vdims = list()
}else{
# Load count, and assume initially that all dimensions start at 1.
v_count = v_vinfo$varsize
v_start = rep(1L,times=length(v_count))
# Find dimensions
v_dname = c(unlist(do.call(rbind,v_vinfo$dim)[,"name"]))
v_vdims = nc_dims[v_dname]
# Find longitude and latitude dimensions
ixy = match(c("lsmlon","lsmlat"),v_dname)
sel = ! is.na(ixy)
# Make sure we only read the data from the nearest neighbour
v_count[ixy[sel]] = c(1,1)[sel]
v_start[ixy[sel]] = in_xy [sel]
names(v_start) = NULL
}#end if (length(v_start) == 0)
# Load information to define the data set.
nc_vars[[v_vname]] = ncvar_def( name = v_vname
, units = v_vunit
, dim = v_vdims
, missval = v_vmiss
, longname = v_vlong
, prec = v_vprec
)#end ncvar_def
# Load the data to a temporary array. We may need to adjust it if the soil depth
# dimension was updated.
nc_datum = ncvar_get( nc = nc_gridin
, varid = v_vname
, start = v_start
, count = v_count
, collapse_degen = FALSE
)#end ncvar_get
# Update the data with the revised number of dimensions if this is a soil variable.
if ("nlevsoi" %in% names(nc_vars[[v_vname]]$dim)){
nc_values[[v_vname]] = nc_datum[,,idx_levsoi,drop=FALSE]
}else{
nc_values[[v_vname]] = nc_datum
}#end if ("nlevsoi" %in% names(nc_vars[[v_vname]]$dim))
# Update the list of dimensions used.
nc_dimused = sort(unique(c(nc_dimused,names(v_vdims))))
}#end for (v in seq_along(nc_vnames))
Append soil depth.
# Append soil layer information to surface data.
nc_vars$ZSOI = ncvar_def( name = "ZSOI"
, units = "m"
, dim = nc_dims[c("nlevsoi")]
, missval = nc_gridin$var$PCT_SAND$missval
, longname = "mid-layer soil depth. Note: in ELM, this will overwrite the default soil depth."
)#end ncvar_def
# Append value to the list of values
nc_values$ZSOI = zsoi
If soil depth was provided, make sure that the value is appended to
the output. This must be done after loading the other variables and
updating the dimensions. The default name in CLM is
zbedrock, whilst ELM uses aveDTB. To ensure
that the surface file is is agnostic to the host land model, we add it
twice. If ELM accepts the pull request for cross-model compatibility, we
can keep one of these defaults only.
if (use_depth_geotiff || use_depth_constant){
# Append depth to bedrock to surface data.
nc_vars$zbedrock = ncvar_def( name = "zbedrock"
, units = "m"
, dim = nc_dims[c("lsmlon","lsmlat")]
, missval = nc_gridin$var$PCT_SAND$missval
, longname = "soil depth"
)#end ncvar_def
nc_vars$aveDTB = ncvar_def( name = "aveDTB"
, units = "m"
, dim = nc_dims[c("lsmlon","lsmlat")]
, missval = nc_gridin$var$PCT_SAND$missval
, longname = "soil depth"
)#end ncvar_def
# Append value to the list of values
nc_values$zbedrock = array(data=ud_bedrock,dim=c(1L,1L))
nc_values$aveDTB = array(data=ud_bedrock,dim=c(1L,1L))
}#end if (use_depth_geotiff || use_depth_constant)
Create dummy variables for unused dimensions. It may not be
necessary, but this ensures every variable is included.
# Make data set
is_unused = which(! (names(nc_dims) %in% nc_dimused))
for (u in is_unused){
# Retrieve dimension information
u_dname = names(nc_dims)[u]
un_dims = nc_dims[[u_dname]]
u_dlong = paste0("Dummy variable to retain dimension ",un_dims$longname)
u_dunit = un_dims$units
u_dvals = un_dims$vals
u_dprec = typeof(u_dvals)
# Load information to define the data set.
nc_vars[[u_dname]] = ncvar_def( name = u_dname
, units = u_dunit
, dim = un_dims
, longname = u_dlong
, prec = u_dprec
)#end ncvar_def
# Load the data
nc_values[[u_dname]] = u_dvals
}#end for (u in is_unused)
In case force_natural is true, we impose percentage of
natural vegetation to 100%.
# Re-set coordinates
nc_values$LONGXY = 0. * nc_values$LONGXY + clon
nc_values$LATIXY = 0. * nc_values$LATIXY + clat
# Re-define percentage of natural vegetation to 100%, and other classes to 0%
if (force_natural){
nc_values$PCT_NATVEG = 0. * nc_values$PCT_NATVEG + 100.
nc_values$PCT_URBAN = 0. * nc_values$PCT_URBAN
nc_values$PCT_CROP = 0. * nc_values$PCT_CROP
nc_values$PCT_WETLAND = 0. * nc_values$PCT_WETLAND
nc_values$PCT_LAKE = 0. * nc_values$PCT_LAKE
nc_values$PCT_GLACIER = 0. * nc_values$PCT_GLACIER
}#end if (force_natural)
In case texture and soil organic content were provided, we use the
site-specific information instead. We test whether the user provided the
information directly or if we are supposed to retrieve it from
SoilGrids.
# In case texture was defined, set the values.
use_soilgrids = any(tolower(c(site_xsand,site_xclay)) %in% "soilgrids")
use_site = (! use_soilgrids) && all(! is.na(c(site_xsand,site_xclay)))
if ( use_soilgrids || use_site){
# Check how to retrieve the soil layers.
if (use_soilgrids){
# Retrieve soil texture from SoilGrids
ud_texture = fetchSoilGrids( x = data.frame(id=1,lat=site_lat,lon=site_lon)
, variables = c("sand","clay")
)#end fetchSoilGrids
# Save variables to vectors.
n_ud_layer = length(ud_texture$id)
ud_ztop = 0.01 * ud_texture$hzdept
ud_zbot = 0.01 * ud_texture$hzdepb
ud_xsand = ud_texture$sandmean
ud_xclay = ud_texture$claymean
}else{
# Retrieve soil texture from user input.
# Define soil layers from user input. If the user did not provide user layers,
# assume the layers to be the same as the default soil grid in the host land models.
if (is.na(site_zbot)){
# Soil depths not specified. Assume the same as the host land model.
ud_ztop = ztop[sequence(nlevsoi)]
ud_zbot = zbot[sequence(nlevsoi)]
}else{
# Soil depths provided. Define the top of each layer, assuming top depth = 0 for
# the topmost layer.
ud_ztop = c(0,site_zbot[-length(site_zbot)])
ud_zbot = site_zbot
}#end if (is.na(site_zbot))
# Number of user-defined layer
n_ud_layer = length(ud_ztop)
# Make sure that the length of user-defined textures and soil layers match, or
# that the user provided a single value to be used everywhere.
if (length(site_xsand) == 1L ){
ud_xsand = rep(site_xsand,times=n_ud_layer)
}else{
ud_xsand = site_xsand
}#end if (length(site_xsand) == 1L )
if (length(site_xclay) == 1L ){
ud_xclay = rep(site_xclay,times=n_ud_layer)
}else{
ud_xclay = site_xclay
}#end if (length(site_xclay) == 1L )
# Make sure the lengths of the input soil textures are all consistent with the
# number of input soil layers.
if (any(c(length(ud_xsand),length(ud_xclay)) != n_ud_layer)){
cat("---~---\n")
cat(" - Length of input sand texture: ",length(site_xsand),".\n",sep="")
cat(" - Length of input clay texture: ",length(site_xclay),".\n",sep="")
cat(" - Number of soil layers: " ,length(n_ud_layer),".\n",sep="")
cat(" Soil texture must be either a scalar or a vector that matches\n")
cat(" the number of layers.\n")
cat("---~---\n")
stop(" Mismatch between input texture and soil layers.")
}#end if (any(c(length(ud_xsand),length(ud_xclay)) != n_ud_layer))
}# end if (use_soilgrids)
# Keep only layers that are within the model grid.
ud_keep = ud_ztop < zbot[nlevsoi]
ud_ztop = ud_ztop [ud_keep]
ud_zbot = ud_zbot [ud_keep]
ud_xsand = ud_xsand[ud_keep]
ud_xclay = ud_xclay[ud_keep]
# We carry over the information of the deepest SoilGrids layer all the way to the
# bottom of the host land model grid
ud_zbot[n_ud_layer] = zbot[nlevsoi]
# Loop through the surface file soil layers and interpolate texture.
for (z in sequence(nlevsoi)){
# For the interpolation, we find how much of the input layer overlaps with the
# host land model layer. We use the overlapping thickness of each input layer as
# weights for the host land model interpolated texture.
ud_lyr_ztop = ifelse( test = ud_zbot < ztop[z]
, yes = ztop[z]
, no = ifelse( test = ud_ztop > zbot[z]
, yes = zbot[z]
, no = pmax(ud_ztop,ztop[z])
)#end ifelse
)#end ifelse
ud_lyr_zbot = ifelse( test = ud_zbot < ztop[z]
, yes = ztop[z]
, no = ifelse( test = ud_ztop > zbot[z]
, yes = zbot[z]
, no = pmin(ud_zbot,zbot[z])
)#end ifelse
)#end ifelse
ud_lyr_thick = ud_lyr_zbot - ud_lyr_ztop
# Find the interpolated texture.
ud_lyr_xsand = weighted.mean(x=ud_xsand,w=ud_lyr_thick)
ud_lyr_xclay = weighted.mean(x=ud_xclay,w=ud_lyr_thick)
ud_lyr_xsilt = 100. - ud_lyr_xsand - ud_lyr_xclay
# Set a vector with textures to test whether or not to use site-specific data
ud_lyr_xtext = c(ud_lyr_xsand,ud_lyr_xsilt,ud_lyr_xclay)
# If the interpolated texture is valid, update the values
if (all( (ud_lyr_xtext >= 0.) & (ud_lyr_xtext <= 100.) ) ){
nc_values$PCT_SAND[,,z] = ud_lyr_xsand
nc_values$PCT_CLAY[,,z] = ud_lyr_xclay
}#end if (all( (ud_lyr_xtext >= 0.) & (ud_lyr_xtext <= 100.) ) )
}#end for (z in sequence(nlevsoi))
# Update global attribute list to inform that the data are site-specific
if (use_soilgrids){
map_soil_texture_file = "SoilGrids 2.0 (obtained for the site through R package soilDB)"
}else{
map_soil_texture_file = "site-specific information"
}#end if (use_soilgrids)
glob_att = modifyList( x = glob_att
, val = list( map_soil_texture_file = map_soil_texture_file )
)#end modifyList
}#end if (use_soilgrids || use_site)
In case texture and soil organic content were provided, we use the
site-specific information instead. We test whether the user provided the
information directly or if we are supposed to retrieve it from
SoilGrids.
# In case texture was defined, set the values.
use_soilgrids = tolower(site_organic) %in% "soilgrids"
use_site = ( ! use_soilgrids ) && (! is.na(site_organic) )
if ( use_soilgrids || use_site){
# Check how to retrieve the soil layers.
if (use_soilgrids){
# Retrieve soil texture from SoilGrids
ud_extract = fetchSoilGrids( x = data.frame(id=1,lat=site_lat,lon=site_lon)
, variables = "ocd"
)#end fetchSoilGrids
# Retrieve the variables
n_ud_layer = length(ud_extract$id)
ud_ztop = 0.01 * ud_extract$hzdept
ud_zbot = 0.01 * ud_extract$hzdepb
ud_organic = ud_extract$ocdmean / 0.58 # To convert it to kg organic matter
}else{
# Retrieve soil organic matter density from user input.
# Define soil layers from user input. If the user did not provide user layers,
# assume the layers to be the same as the default soil grid in the host land models.
if (is.na(site_zbot)){
# Soil depths not specified. Assume the same as the host land model.
ud_ztop = ztop[sequence(nlevsoi)]
ud_zbot = zbot[sequence(nlevsoi)]
}else{
# Soil depths provided. Define the top of each layer, assuming top depth = 0 for
# the topmost layer.
ud_ztop = c(0,site_zbot[-length(site_zbot)])
ud_zbot = site_zbot
}#end if (is.na(site_zbot))
# Number of user-defined layer
n_ud_layer = length(ud_ztop)
# If the user provided constant soil organic matter density, turn it into a vector.
if (length(site_organic) == 1L ){
ud_organic = rep(site_organic,times=n_ud_layer)
}else{
ud_organic = site_organic
}#end if (length(site_organic) == 1L )
# Make sure the length of the input data is consistent with the number of input
# soil layers.
if (length(ud_organic) != n_ud_layer){
cat("---~---\n")
cat(" - Length of input organic matter: ",length(site_organic),".\n",sep="")
cat(" - Number of soil layers: " ,length(n_ud_layer) ,".\n",sep="")
cat(" Site organic matter density must be either a scalar or a vector\n")
cat(" that matches the number of layers.\n")
cat("---~---\n")
stop(" Mismatch between input texture and soil layers.")
}#end if (length(ud_organic) != n_ud_layer)
}# end if (use_soilgrids)
# Keep only layers that are within the model grid.
ud_keep = ud_ztop < zbot[nlevsoi]
ud_ztop = ud_ztop [ud_keep]
ud_zbot = ud_zbot [ud_keep]
ud_xsand = ud_xsand[ud_keep]
ud_xclay = ud_xclay[ud_keep]
# We carry over the information of the deepest SoilGrids layer all the way to the
# bottom of the host land model grid
ud_zbot[n_ud_layer] = zbot[nlevsoi]
# Loop through the surface file soil layers and interpolate soil organic matter
# density.
for (z in sequence(nlevsoi)){
# For the interpolation, we find how much of the input layer overlaps with the
# host land model layer. We use the overlapping thickness of each input layer as
# weights for the host land model interpolated texture.
ud_lyr_ztop = ifelse( test = ud_zbot < ztop[z]
, yes = ztop[z]
, no = ifelse( test = ud_ztop > zbot[z]
, yes = zbot[z]
, no = pmax(ud_ztop,ztop[z])
)#end ifelse
)#end ifelse
ud_lyr_zbot = ifelse( test = ud_zbot < ztop[z]
, yes = ztop[z]
, no = ifelse( test = ud_ztop > zbot[z]
, yes = zbot[z]
, no = pmin(ud_zbot,zbot[z])
)#end ifelse
)#end ifelse
ud_lyr_thick = ud_lyr_zbot - ud_lyr_ztop
# Find the interpolated soil organic matter density.
nc_values$ORGANIC[,,z] = weighted.mean(x=ud_organic,w=ud_lyr_thick)
}#end for (z in sequence(nlevsoi))
# Update global attribute list to inform that the data are site-specific
if (use_soilgrids){
map_soil_organic_file = "SoilGrids 2.0 (obtained for the site through R package soilDB)"
}else{
map_soil_organic_file = "site-specific information"
}#end if (use_soilgrids)
glob_att = modifyList( x = glob_att
, val = list( map_soil_organic_file = map_soil_organic_file )
)#end modifyList
}#end if (use_soilgrids || use_site)
Variable T_BUILDING_MAX used to be part of every surface
data. Apparently CLM dropped this variable, but ELM continues to require
it. If the reference surface data set is a recent one, this variable
will be missing and the surface data file will not work with ELM. To
address this, we use T_BUILDING_MIN to re-generate
T_BUILDING_MAX, by applying an offset compatible with older
surface data files.
if (! ("T_BUILDING_MAX" %in% names(nc_vars))){
# Create maximum interior building temperature
nc_vars$T_BUILDING_MAX = ncvar_def( name = "T_BUILDING_MAX"
, units = "K"
, dim = nc_dims[c("lsmlon","lsmlat","numurbl")]
, missval = nc_gridin$var$T_BUILDING_MIN$missval
, longname = "maximum interior building temperature"
)#end ncvar_def
# Append value to the list of values
dT_building = c(15,rep(39,times=nc_dims$numurbl$len-1))
nc_values$T_BUILDING_MAX = nc_values$T_BUILDING_MIN + dT_building
}# end if (! ("T_BUILDING_MAX" %in% names(nc_vars)))
Create the output files.
# Create output netCDF
nc_surfout = nc_create(filename = output_surfdata,vars=nc_vars)
# Add global attributes
for (a in seq_along(glob_att)){
# Retrieve attribute information
a_name = names(glob_att)[a]
a_value = glob_att[[a_name]]
# Include global attribute
dummy = ncatt_put(nc = nc_surfout, varid=0,attname=a_name,attval=a_value)
}#end for (a in seq_along(glob_att))
# Add variables
for (v in seq_along(nc_values)){
# Retrieve variable information
v_name = names(nc_values)[v]
v_value = nc_values[[v_name]]
# Add variable values
dummy = ncvar_put(nc = nc_surfout,varid=v_name,vals=v_value)
}#end for (v in seq_along(nc_values))
# Close input and output files
dummy = nc_close(nc = nc_gridin )
dummy = nc_close(nc = nc_surfout)