[Dart-dev] [9015] DART/trunk: This has an interpolation scheme that honors what the user
nancy at ucar.edu
nancy at ucar.edu
Fri Nov 6 15:20:29 MST 2015
Revision: 9015
Author: thoar
Date: 2015-11-06 15:20:29 -0700 (Fri, 06 Nov 2015)
Log Message:
-----------
This has an interpolation scheme that honors what the user
specifies in the clm_variables namelist. previous version
had hardwired variable names for the KINDS, which may not have
been in the DART state to begin with.
Modified Paths:
--------------
DART/trunk/models/clm/model_mod.f90
DART/trunk/models/clm/model_mod.html
DART/trunk/models/clm/model_mod_check.f90
DART/trunk/models/clm/work/input.nml
DART/trunk/obs_kind/DEFAULT_obs_kind_mod.F90
-------------- next part --------------
Modified: DART/trunk/models/clm/model_mod.f90
===================================================================
--- DART/trunk/models/clm/model_mod.f90 2015-11-06 21:01:09 UTC (rev 9014)
+++ DART/trunk/models/clm/model_mod.f90 2015-11-06 22:20:29 UTC (rev 9015)
@@ -51,17 +51,27 @@
file_exist, find_textfile_dims, file_to_text, &
open_file, close_file
-use obs_kind_mod, only : KIND_SOIL_TEMPERATURE, &
- KIND_SOIL_MOISTURE, &
- KIND_LIQUID_WATER, &
- KIND_ICE, &
- KIND_SNOWCOVER_FRAC, &
- KIND_SNOW_THICKNESS, &
- KIND_LEAF_CARBON, &
- KIND_WATER_TABLE_DEPTH, &
- KIND_GEOPOTENTIAL_HEIGHT,&
- paramname_length, &
- get_raw_obs_kind_index
+use obs_kind_mod, only : KIND_SOIL_TEMPERATURE, &
+ KIND_SOIL_MOISTURE, &
+ KIND_LIQUID_WATER, &
+ KIND_ICE, &
+ KIND_SNOWCOVER_FRAC, &
+ KIND_SNOW_THICKNESS, &
+ KIND_LEAF_CARBON, &
+ KIND_LEAF_AREA_INDEX, &
+ KIND_WATER_TABLE_DEPTH, &
+ KIND_GEOPOTENTIAL_HEIGHT, &
+ KIND_VEGETATION_TEMPERATURE, &
+ KIND_FPAR, &
+ KIND_FPAR_SUNLIT_DIRECT, &
+ KIND_FPAR_SUNLIT_DIFFUSE, &
+ KIND_FPAR_SHADED_DIRECT, &
+ KIND_FPAR_SHADED_DIFFUSE, &
+ KIND_FPAR_SHADED_DIRECT, &
+ KIND_FPAR_SHADED_DIFFUSE, &
+ paramname_length, &
+ get_raw_obs_kind_index, &
+ get_raw_obs_kind_name
use mpi_utilities_mod, only: my_task_id
@@ -214,7 +224,6 @@
type(progvartype), dimension(max_state_variables) :: progvar
-
!----------------------------------------------------------------------
! how many and which columns are in each gridcell
!----------------------------------------------------------------------
@@ -269,17 +278,18 @@
! scan along longitudes, then
! move to next latitude.
-integer :: Ngridcell = -1 ! Number of gridcells containing land
-integer :: Nlandunit = -1 ! Number of land units
-integer :: Ncolumn = -1 ! Number of columns
-integer :: Npft = -1 ! Number of plant functional types
-integer :: Nlevlak = -1 ! Number of
-integer :: Nlevsno = -1 ! Number of snow levels
-integer :: Nlevsno1 = -1 ! Number of snow level ... interfaces?
-integer :: Nlevtot = -1 ! Number of total levels
-integer :: Nnumrad = -1 ! Number of
-integer :: Nrtmlon = -1 ! Number of river transport model longitudes
-integer :: Nrtmlat = -1 ! Number of river transport model latitudes
+integer :: ngridcell = -1 ! Number of gridcells containing land
+integer :: nlandunit = -1 ! Number of land units
+integer :: ncolumn = -1 ! Number of columns
+integer :: npft = -1 ! Number of plant functional types
+integer :: nlevlak = -1 ! Number of
+integer :: nlevsno = -1 ! Number of snow levels
+integer :: nlevsno1 = -1 ! Number of snow level ... interfaces?
+integer :: nlevtot = -1 ! Number of total levels
+integer :: nnumrad = -1 ! Number of
+integer :: nrtmlon = -1 ! Number of river transport model longitudes
+integer :: nrtmlat = -1 ! Number of river transport model latitudes
+integer :: nlevcan = -1 ! Number of canopy layers (*XY*)
integer, allocatable, dimension(:) :: grid1d_ixy, grid1d_jxy ! 2D lon/lat index of corresponding gridcell
integer, allocatable, dimension(:) :: land1d_ixy, land1d_jxy ! 2D lon/lat index of corresponding gridcell
@@ -532,13 +542,13 @@
call get_sparse_dims(ncid, clm_restart_filename, 'open')
-allocate(grid1d_ixy(Ngridcell), grid1d_jxy(Ngridcell))
-allocate(land1d_ixy(Nlandunit), land1d_jxy(Nlandunit), land1d_wtxy(Nlandunit))
-allocate(cols1d_ixy(Ncolumn), cols1d_jxy(Ncolumn))
-allocate(cols1d_wtxy(Ncolumn), cols1d_ityplun(Ncolumn))
-allocate(pfts1d_ixy(Npft), pfts1d_jxy(Npft) , pfts1d_wtxy(Npft))
-allocate(levtot(Nlevtot))
-if (Nlevsno > 0) allocate(zsno(Nlevsno,Ncolumn))
+allocate(grid1d_ixy(ngridcell), grid1d_jxy(ngridcell))
+allocate(land1d_ixy(nlandunit), land1d_jxy(nlandunit), land1d_wtxy(nlandunit))
+allocate(cols1d_ixy(ncolumn), cols1d_jxy(ncolumn))
+allocate(cols1d_wtxy(ncolumn), cols1d_ityplun(ncolumn))
+allocate(pfts1d_ixy(npft), pfts1d_jxy(npft) , pfts1d_wtxy(npft))
+allocate(levtot(nlevtot))
+if (nlevsno > 0) allocate(zsno(nlevsno,ncolumn))
call get_sparse_geog(ncid, clm_restart_filename, 'close')
@@ -753,6 +763,9 @@
! 1D variables are usually from the restart file
! FIXME this same logic is used for 2D variables from the 'vector' file which
! has a singleton dimension of 'time'
+ ! What if I check on the rank of the variable instead of the number of dimensions ...
+ ! If I require 'time' to be the unlimited dimension, it will always be 'last',
+ ! so the first N dimensions and the first N ranks are identical ...
if (progvar(ivar)%numdims == 1) then
@@ -763,10 +776,10 @@
endif
SELECT CASE ( trim(progvar(ivar)%dimnames(1)) )
- CASE ("gridcell")
+ CASE ("gridcell","lndgrid")
do i = 1, progvar(ivar)%dimlens(1)
xi = grid1d_ixy(i)
- xj = grid1d_jxy(i) ! always unity if unstructured
+ xj = grid1d_jxy(i) ! nnnnn_jxy(:) always 1 if unstructured
if (unstructured) then
lonixy( indx) = xi
latjxy( indx) = xi
@@ -782,7 +795,7 @@
CASE ("landunit")
do i = 1, progvar(ivar)%dimlens(1)
xi = land1d_ixy(i)
- xj = land1d_jxy(i) ! always unity if unstructured
+ xj = land1d_jxy(i) ! nnnnn_jxy(:) always 1 if unstructured
if (unstructured) then
lonixy( indx) = xi
latjxy( indx) = xi
@@ -798,7 +811,7 @@
CASE ("column")
do i = 1, progvar(ivar)%dimlens(1)
xi = cols1d_ixy(i)
- xj = cols1d_jxy(i) ! always unity if unstructured
+ xj = cols1d_jxy(i) ! nnnnn_jxy(:) always 1 if unstructured
if (unstructured) then
lonixy( indx) = xi
latjxy( indx) = xi
@@ -814,7 +827,7 @@
CASE ("pft")
do i = 1, progvar(ivar)%dimlens(1)
xi = pfts1d_ixy(i)
- xj = pfts1d_jxy(i) ! always unity if unstructured
+ xj = pfts1d_jxy(i) ! nnnnn_jxy(:) always 1 if unstructured
if (unstructured) then
lonixy( indx) = xi
latjxy( indx) = xi
@@ -855,7 +868,7 @@
if ((debug > 8) .and. do_output()) write(*,*)'length grid1d_ixy ',size(grid1d_ixy)
do j = 1, progvar(ivar)%dimlens(2)
xi = grid1d_ixy(j)
- xj = grid1d_jxy(j) ! always unity if unstructured
+ xj = grid1d_jxy(j) ! nnnnn_jxy(:) always 1 if unstructured
do i = 1, progvar(ivar)%dimlens(1)
if (unstructured) then
lonixy( indx) = xi
@@ -874,7 +887,7 @@
if ((debug > 8) .and. do_output()) write(*,*)'length land1d_ixy ',size(land1d_ixy)
do j = 1, progvar(ivar)%dimlens(2)
xi = land1d_ixy(j)
- xj = land1d_jxy(j) ! always unity if unstructured
+ xj = land1d_jxy(j) ! nnnnn_jxy(:) always 1 if unstructured
do i = 1, progvar(ivar)%dimlens(1)
if (unstructured) then
lonixy( indx) = xi
@@ -904,7 +917,7 @@
call fill_levels(progvar(ivar)%dimnames(1),j,progvar(ivar)%dimlens(1),levtot)
xi = cols1d_ixy(j)
- xj = cols1d_jxy(j) ! always unity if unstructured
+ xj = cols1d_jxy(j) ! nnnnn_jxy(:) always 1 if unstructured
VERTICAL : do i = 1, progvar(ivar)%dimlens(1)
levels( indx) = levtot(i)
if (unstructured) then
@@ -924,7 +937,7 @@
if ((debug > 8) .and. do_output()) write(*,*)'length pfts1d_ixy ',size(pfts1d_ixy)
do j = 1, progvar(ivar)%dimlens(2)
xi = pfts1d_ixy(j)
- xj = pfts1d_jxy(j) ! always unity if unstructured
+ xj = pfts1d_jxy(j) ! nnnnn_jxy(:) always 1 if unstructured
do i = 1, progvar(ivar)%dimlens(1)
if (unstructured) then
lonixy( indx) = xi
@@ -942,12 +955,13 @@
CASE ("time")
! The vector history files can have things 'pft,time' or 'column,time'
+ ! The single-column history files can have things 'lndgrid,time', 'pft,time' or 'column,time'
SELECT CASE ( trim(progvar(ivar)%dimnames(1)) )
CASE ("pft")
do i = 1, progvar(ivar)%dimlens(1)
xi = pfts1d_ixy(i)
- xj = pfts1d_jxy(i) ! always unity if unstructured
+ xj = pfts1d_jxy(i) ! nnnnn_jxy(:) always 1 if unstructured
if (unstructured) then
lonixy( indx) = xi
latjxy( indx) = xi
@@ -963,7 +977,7 @@
CASE ("column")
do i = 1, progvar(ivar)%dimlens(1)
xi = cols1d_ixy(i)
- xj = cols1d_jxy(i) ! always unity if unstructured
+ xj = cols1d_jxy(i) ! nnnnn_jxy(:) always 1 if unstructured
if (unstructured) then
lonixy( indx) = xi
latjxy( indx) = xi
@@ -976,6 +990,22 @@
indx = indx + 1
enddo
+ CASE ("lndgrid")
+ do i = 1, progvar(ivar)%dimlens(1)
+ xi = cols1d_ixy(i)
+ xj = cols1d_jxy(i) ! nnnnn_jxy(:) always 1 if unstructured
+ if (unstructured) then
+ lonixy( indx) = xi
+ latjxy( indx) = xi
+ landarea(indx) = AREA1D(xi) * LANDFRAC1D(xi)
+ else
+ lonixy( indx) = xi
+ latjxy( indx) = xj
+ landarea(indx) = AREA2D(xi,xj) * LANDFRAC2D(xi,xj)
+ endif
+ indx = indx + 1
+ enddo
+
CASE DEFAULT
write(string1,*)'(2d) unknown first Dimension name '//trim(progvar(ivar)%dimnames(1))// &
' while trying to create metadata for '//trim(progvar(ivar)%varname)
@@ -1176,6 +1206,7 @@
! for the dimensions and coordinate variables
integer :: nlonDimID
integer :: nlatDimID
+integer :: lndgridDimID
integer :: gridcellDimID
integer :: landunitDimID
integer :: columnDimID
@@ -1186,6 +1217,7 @@
integer :: levsno1DimID
integer :: levtotDimID
integer :: numradDimID
+integer :: levcanDimID
! for the prognostic variables
integer :: ivar, VarID
@@ -1325,7 +1357,11 @@
call nc_check(nf90_def_dim(ncid=ncFileID, name='lon', len = nlon, &
dimid = nlonDimID),'nc_write_model_atts', 'lon def_dim '//trim(filename))
call nc_check(nf90_def_dim(ncid=ncFileID, name='lat', len = nlat, &
- dimid = NlatDimID),'nc_write_model_atts', 'lat def_dim '//trim(filename))
+ dimid = nlatDimID),'nc_write_model_atts', 'lat def_dim '//trim(filename))
+ if (unstructured) then
+ call nc_check(nf90_def_dim(ncid=ncFileID, name='lndgrid', len = ngridcell, &
+ dimid = lndgridDimID),'nc_write_model_atts', 'lndgrid def_dim '//trim(filename))
+ endif
call nc_check(nf90_def_dim(ncid=ncFileID, name='gridcell', len = ngridcell, &
dimid = gridcellDimID),'nc_write_model_atts', 'gridcell def_dim '//trim(filename))
@@ -1347,6 +1383,9 @@
dimid = levtotDimID),'nc_write_model_atts', 'levtot def_dim '//trim(filename))
call nc_check(nf90_def_dim(ncid=ncFileID, name='numrad', len = nnumrad, &
dimid = numradDimID),'nc_write_model_atts', 'numrad def_dim '//trim(filename))
+ if (nlevcan > 0) &
+ call nc_check(nf90_def_dim(ncid=ncFileID, name='levcan', len = nlevcan, &
+ dimid = levcanDimID),'nc_write_model_atts', 'levcan def_dim'//trim(filename))
!----------------------------------------------------------------------------
! Create the (empty) Coordinate Variables and the Attributes
@@ -1354,7 +1393,7 @@
! Grid Longitudes
call nc_check(nf90_def_var(ncFileID,name='lon', xtype=nf90_real, &
- dimids=(/ NlonDimID /), varid=VarID),&
+ dimids=(/ nlonDimID /), varid=VarID),&
'nc_write_model_atts', 'lon def_var '//trim(filename))
call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', 'coordinate longitude'), &
'nc_write_model_atts', 'lon long_name '//trim(filename))
@@ -1392,11 +1431,11 @@
! grid cell areas
if (unstructured) then
call nc_check(nf90_def_var(ncFileID,name='area', xtype=nf90_real, &
- dimids=(/ NlonDimID /), varid=VarID),&
+ dimids=(/ nlonDimID /), varid=VarID),&
'nc_write_model_atts', 'area def_var '//trim(filename))
else
call nc_check(nf90_def_var(ncFileID,name='area', xtype=nf90_real, &
- dimids=(/ NlonDimID,nlatDimID /), varid=VarID),&
+ dimids=(/ nlonDimID,nlatDimID /), varid=VarID),&
'nc_write_model_atts', 'area def_var '//trim(filename))
endif
call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', 'grid cell areas'), &
@@ -1407,11 +1446,11 @@
! grid cell land fractions
if (unstructured) then
call nc_check(nf90_def_var(ncFileID,name='landfrac', xtype=nf90_real, &
- dimids=(/ NlonDimID /), varid=VarID),&
+ dimids=(/ nlonDimID /), varid=VarID),&
'nc_write_model_atts', 'landfrac def_var '//trim(filename))
else
call nc_check(nf90_def_var(ncFileID,name='landfrac', xtype=nf90_real, &
- dimids=(/ NlonDimID,nlatDimID /), varid=VarID),&
+ dimids=(/ nlonDimID,nlatDimID /), varid=VarID),&
'nc_write_model_atts', 'landfrac def_var '//trim(filename))
endif
call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', 'land fraction'), &
@@ -1952,7 +1991,7 @@
! number of snow layers
! time of restart file
-allocate(snlsno(Ncolumn))
+allocate(snlsno(ncolumn))
call nc_check(nf90_open(trim(clm_restart_filename), NF90_NOWRITE, ncid), &
'clm_to_dart_state_vector', 'open SNLSNO'//clm_restart_filename)
call nc_check(nf90_inq_varid(ncid,'SNLSNO', VarID), &
@@ -2005,7 +2044,7 @@
call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), len=dimlen), &
'clm_to_dart_state_vector', string1)
- ! Time dimension will be unity in progvar, but not necessarily
+ ! Time dimension will be 1 in progvar, but not necessarily
! in origin file. We only want a single matching time.
! static_init_model() only reserves space for a single time.
@@ -2061,7 +2100,7 @@
do j = 1, nj ! loop over columns
numsnowlevels = abs(snlsno(j))
- do i = 1, Nlevsno - numsnowlevels ! loop over layers
+ do i = 1, nlevsno - numsnowlevels ! loop over layers
data_2d_array(i,j) = MISSING_R8
enddo
enddo
@@ -2071,7 +2110,7 @@
do j = 1, nj ! loop over columns
numsnowlevels = abs(snlsno(j))
- do i = 1, Nlevsno - numsnowlevels ! loop over layers
+ do i = 1, nlevsno - numsnowlevels ! loop over layers
data_2d_array(i,j) = MISSING_R8
enddo
enddo
@@ -2362,6 +2401,7 @@
real(r8) :: llon, llat, lheight
real(r8) :: interp_val_2
integer :: istatus_2
+character(len=paramname_length) :: kind_string
if ( .not. module_initialized ) call static_init_model
@@ -2384,40 +2424,11 @@
if ((debug > 6) .and. do_output()) print *, 'requesting interpolation at ', llon, llat, lheight
-! FIXME may be better to check the %maxlevels and kick the interpolation to the
-! appropriate routine based on that ... or check the dimnames for the
-! vertical coordinate ...
+! Some applications just need to know the number of vertical levels.
+! This is done by trying to 'interpolate' height on a large number of levels.
+! When the interpolation fails, you've gone one level too far.
-if (obs_kind == KIND_SOIL_TEMPERATURE) then
- call get_grid_vertval(x, location, 'T_SOISNO', interp_val, istatus )
-
-elseif (obs_kind == KIND_SOIL_MOISTURE) then
- ! TJH FIXME - actually ROLAND FIXME
- ! This is terrible ... the COSMOS operator wants m3/m3 ... CLM is kg/m2
- call get_grid_vertval(x, location, 'H2OSOI_LIQ',interp_val , istatus )
- call get_grid_vertval(x, location, 'H2OSOI_ICE',interp_val_2, istatus_2 )
- if ((istatus == 0) .and. (istatus_2 == 0)) then
- interp_val = interp_val + interp_val_2
- else
- interp_val = MISSING_R8
- istatus = 6
- endif
-
-elseif (obs_kind == KIND_LIQUID_WATER ) then
- call get_grid_vertval(x, location, 'H2OSOI_LIQ',interp_val, istatus )
-elseif (obs_kind == KIND_ICE ) then
- call get_grid_vertval(x, location, 'H2OSOI_ICE',interp_val, istatus )
-elseif (obs_kind == KIND_SNOWCOVER_FRAC ) then
- call compute_gridcell_value(x, location, 'frac_sno', interp_val, istatus)
-elseif (obs_kind == KIND_LEAF_CARBON ) then
- call compute_gridcell_value(x, location, 'leafc', interp_val, istatus)
-elseif (obs_kind == KIND_WATER_TABLE_DEPTH ) then
- call compute_gridcell_value(x, location, 'ZWT', interp_val, istatus)
-elseif (obs_kind == KIND_SNOW_THICKNESS ) then
- write(string1,*)'model_interpolate for DZSNO not written yet.'
- call error_handler(E_ERR,'model_interpolate',string1,source,revision,revdate)
- istatus = 5
-elseif ((obs_kind == KIND_GEOPOTENTIAL_HEIGHT) .and. vert_is_level(location)) then
+if ((obs_kind == KIND_GEOPOTENTIAL_HEIGHT) .and. vert_is_level(location)) then
if (nint(lheight) > nlevgrnd) then
interp_val = MISSING_R8
istatus = 1
@@ -2425,12 +2436,51 @@
interp_val = LEVGRND(nint(lheight))
istatus = 0
endif
-else
- write(string1,*)'model_interpolate not written for (integer) kind ',obs_kind
- call error_handler(E_ERR,'model_interpolate',string1,source,revision,revdate)
- istatus = 5
+ return ! Early Return
endif
+! Standard method of interpolation.
+! get_grid_vertval() for quantities that have a vertical profile.
+! compute_gridcell_value() for quantities computed for the entire gridcell.
+
+select case( obs_kind )
+
+ case ( KIND_SOIL_MOISTURE )
+
+ ! TJH FIXME the COSMOS operator wants m3/m3 ... CLM is kg/m2
+ ! TJH FIXME ... what units do other operators want/need ... ugly
+ call get_grid_vertval(x, location, KIND_LIQUID_WATER, interp_val, istatus)
+ call get_grid_vertval(x, location, KIND_ICE, interp_val_2, istatus_2)
+ if ((istatus == 0) .and. (istatus_2 == 0)) then
+ interp_val = interp_val + interp_val_2
+ else
+ interp_val = MISSING_R8
+ istatus = 6
+ endif
+
+ case ( KIND_SOIL_TEMPERATURE, KIND_LIQUID_WATER, KIND_ICE )
+
+ call get_grid_vertval(x, location, obs_kind, interp_val, istatus)
+
+ case ( KIND_SNOWCOVER_FRAC, KIND_LEAF_AREA_INDEX, KIND_LEAF_CARBON, &
+ KIND_WATER_TABLE_DEPTH, KIND_VEGETATION_TEMPERATURE, KIND_FPAR, &
+ KIND_FPAR_SUNLIT_DIRECT, KIND_FPAR_SUNLIT_DIFFUSE, &
+ KIND_FPAR_SHADED_DIRECT, KIND_FPAR_SHADED_DIFFUSE)
+
+ call compute_gridcell_value(x, location, obs_kind, interp_val, istatus)
+
+ case default
+
+ kind_string = get_raw_obs_kind_name(obs_kind)
+
+ write(string1,*)'not written for (integer) kind ',obs_kind
+ write(string2,*)'AKA '//trim(kind_string)
+ call error_handler(E_ERR, 'model_interpolate', string1, &
+ source, revision, revdate, text2=string2)
+ istatus = 5
+
+end select
+
if ((debug > 6) .and. do_output()) write(*,*)'interp_val ',interp_val
end subroutine model_interpolate
@@ -2439,7 +2489,7 @@
!------------------------------------------------------------------
-subroutine compute_gridcell_value(x, location, varstring, interp_val, istatus)
+subroutine compute_gridcell_value(x, location, kind_index, interp_val, istatus)
!
! Each gridcell may contain values for several land units, each land unit may contain
! several columns, each column may contain several pft's. BUT this routine never
@@ -2450,7 +2500,7 @@
real(r8), intent(in) :: x(:) ! state vector
type(location_type), intent(in) :: location ! location somewhere in a grid cell
-character(len=*), intent(in) :: varstring ! frac_sno, leafc
+integer, intent(in) :: kind_index ! KIND in DART state needed for interpolation
real(r8), intent(out) :: interp_val ! area-weighted result
integer, intent(out) :: istatus ! error code (0 == good)
@@ -2479,13 +2529,13 @@
loc_lat = loc(2)
! determine the portion of interest of the state vector
-ivar = findVarIndex(varstring, 'compute_gridcell_value')
+ivar = findKindIndex(kind_index, 'compute_gridcell_value')
index1 = progvar(ivar)%index1 ! in the DART state vector, start looking here
indexN = progvar(ivar)%indexN ! in the DART state vector, stop looking here
! BOMBPROOFING - check for a vertical dimension for this variable
if (progvar(ivar)%maxlevels > 1) then
- write(string1, *)'Variable '//trim(varstring)//' cannot use this routine.'
+ write(string1, *)'Variable '//trim(progvar(ivar)%varname)//' cannot use this routine.'
write(string2, *)'use get_grid_vertval() instead.'
call error_handler(E_ERR,'compute_gridcell_value', string1, &
source, revision, revdate, text2=string2)
@@ -2543,7 +2593,7 @@
istatus = 0
else
if ((debug > 4) .and. do_output()) then
- write(string1, *)'Variable '//trim(varstring)//' had no viable data'
+ write(string1, *)'Variable '//trim(progvar(ivar)%varname)//' had no viable data'
write(string2, *)'at gridcell ilon/jlat = (',gridloni,',',gridlatj,')'
write(string3, *)'obs lon/lat = (',loc_lon,',',loc_lat,')'
call error_handler(E_MSG,'compute_gridcell_value', string1, &
@@ -2564,7 +2614,7 @@
!------------------------------------------------------------------
-subroutine get_grid_vertval(x, location, varstring, interp_val, istatus)
+subroutine get_grid_vertval(x, location, kind_index, interp_val, istatus)
!
! Calculate the expected vertical value for the gridcell.
! Each gridcell value is an area-weighted value of an unknown number of
@@ -2574,7 +2624,7 @@
real(r8), intent(in) :: x(:) ! state vector
type(location_type), intent(in) :: location ! location somewhere in a grid cell
-character(len=*), intent(in) :: varstring ! T_SOISNO, H2OSOI_LIQ, H2OSOI_ICE
+integer, intent(in) :: kind_index
real(r8), intent(out) :: interp_val ! area-weighted result
integer, intent(out) :: istatus ! error code (0 == good)
@@ -2592,6 +2642,8 @@
real(r8), allocatable, dimension(:) :: above, below
real(r8), allocatable, dimension(:,:) :: myarea
+character(len=paramname_length) :: varstring
+
if ( .not. module_initialized ) call static_init_model
! Let's assume failure. Set return val to missing, then the code can
@@ -2617,10 +2669,12 @@
endif
! determine the portion of interest of the state vector
-ivar = findVarIndex(varstring, 'get_grid_vertval')
+ivar = findKindIndex(kind_index, 'get_grid_vertval')
index1 = progvar(ivar)%index1 ! in the DART state vector, start looking here
indexN = progvar(ivar)%indexN ! in the DART state vector, stop looking here
+varstring = progvar(ivar)%varname ! used in a lot of error messages
+
! BOMBPROOFING - check for a vertical dimension for this variable
if (progvar(ivar)%maxlevels < 2) then
write(string1, *)'Variable '//trim(varstring)//' should not use this routine.'
@@ -2713,34 +2767,28 @@
if ( latjxy(indexi) /= gridlatj ) cycle ELEMENTS
if ( x(indexi) == MISSING_R8) cycle ELEMENTS
-! write(*,*)'level ',indexi,' is ',levels(indexi),' location depth is ',loc_lev
-
- if (levels(indexi) == depthabove) then
+ if (levels(indexi) == depthabove) then
counter1 = counter1 + 1
above( counter1) = x(indexi)
myarea(counter1,1) = landarea(indexi)
endif
- if (levels(indexi) == depthbelow) then
+
+ if (levels(indexi) == depthbelow) then
counter2 = counter2 + 1
below( counter2) = x(indexi)
myarea(counter2,2) = landarea(indexi)
endif
- if ((levels(indexi) /= depthabove) .and. &
- (levels(indexi) /= depthbelow)) then
- cycle ELEMENTS
- endif
-
if ((debug > 4) .and. do_output()) then
- write(*,*)
- write(*,*)'gridcell location match at statevector index',indexi
- write(*,*)'statevector value is (',x(indexi),')'
- write(*,*)'area is (',landarea(indexi),')'
- write(*,*)'LON index is (',lonixy(indexi),')'
- write(*,*)'LAT index is (',latjxy(indexi),')'
- write(*,*)'gridcell LON is (',LON(gridloni),')'
- write(*,*)'gridcell LAT is (',LAT(gridlatj),')'
- write(*,*)'depth is (',levels(indexi),')'
+ write(*,*)
+ write(*,*)'gridcell location match at statevector index',indexi
+ write(*,*)'statevector value is (',x(indexi),')'
+ write(*,*)'area is (',landarea(indexi),')'
+ write(*,*)'LON index is (',lonixy(indexi),')'
+ write(*,*)'LAT index is (',latjxy(indexi),')'
+ write(*,*)'gridcell LON is (',LON(gridloni),')'
+ write(*,*)'gridcell LAT is (',LAT(gridlatj),')'
+ write(*,*)'depth is (',levels(indexi),')'
endif
enddo ELEMENTS
@@ -2766,7 +2814,7 @@
value_above = sum(above(1:counter1) * myarea(1:counter1,1))
else
write(string1, *)'Variable '//trim(varstring)//' had no viable data above'
- write(string2, *)'at gridcell lon/lat/lev = (',gridloni,',',gridlatj,',',depthabove,')'
+ write(string2, *)'at gridcell lon/lat/lev = (',LON(gridloni),',',LAT(gridlatj),',',depthabove,')'
write(string3, *)'obs lon/lat/lev (',loc_lon,',',loc_lat,',',loc_lev,')'
call error_handler(E_ERR,'get_grid_vertval', string1, &
source, revision, revdate, text2=string2,text3=string3)
@@ -2782,7 +2830,7 @@
value_below = sum(below(1:counter2) * myarea(1:counter2,2))
else
write(string1, *)'Variable '//trim(varstring)//' had no viable data below'
- write(string2, *)'at gridcell lon/lat/lev = (',gridloni,',',gridlatj,',',depthbelow,')'
+ write(string2, *)'at gridcell lon/lat/lev = (',LON(gridloni),',',LAT(gridlatj),',',depthbelow,')'
write(string3, *)'obs lon/lat/lev (',loc_lon,',',loc_lat,',',loc_lev,')'
call error_handler(E_ERR,'get_grid_vertval', string1, &
source, revision, revdate, text2=string2,text3=string3)
@@ -3246,22 +3294,22 @@
call nc_check(nf90_inq_dimid(ncid, 'gridcell', dimid), &
'get_sparse_dims','inq_dimid gridcell '//trim(fname))
-call nc_check(nf90_inquire_dimension(ncid, dimid, len=Ngridcell), &
+call nc_check(nf90_inquire_dimension(ncid, dimid, len=ngridcell), &
'get_sparse_dims','inquire_dimension gridcell '//trim(fname))
call nc_check(nf90_inq_dimid(ncid, 'landunit', dimid), &
'get_sparse_dims','inq_dimid landunit '//trim(fname))
-call nc_check(nf90_inquire_dimension(ncid, dimid, len=Nlandunit), &
+call nc_check(nf90_inquire_dimension(ncid, dimid, len=nlandunit), &
'get_sparse_dims','inquire_dimension landunit '//trim(fname))
call nc_check(nf90_inq_dimid(ncid, 'column', dimid), &
'get_sparse_dims','inq_dimid column '//trim(fname))
-call nc_check(nf90_inquire_dimension(ncid, dimid, len=Ncolumn), &
+call nc_check(nf90_inquire_dimension(ncid, dimid, len=ncolumn), &
'get_sparse_dims','inquire_dimension column '//trim(fname))
call nc_check(nf90_inq_dimid(ncid, 'pft', dimid), &
'get_sparse_dims','inq_dimid pft '//trim(fname))
-call nc_check(nf90_inquire_dimension(ncid, dimid, len=Npft), &
+call nc_check(nf90_inquire_dimension(ncid, dimid, len=npft), &
'get_sparse_dims','inquire_dimension pft '//trim(fname))
call nc_check(nf90_inq_dimid(ncid, 'levgrnd', dimid), &
@@ -3277,26 +3325,34 @@
call nc_check(nf90_inq_dimid(ncid, 'levlak', dimid), &
'get_sparse_dims','inq_dimid levlak '//trim(fname))
-call nc_check(nf90_inquire_dimension(ncid, dimid, len=Nlevlak), &
+call nc_check(nf90_inquire_dimension(ncid, dimid, len=nlevlak), &
'get_sparse_dims','inquire_dimension levlak '//trim(fname))
call nc_check(nf90_inq_dimid(ncid, 'levtot', dimid), &
'get_sparse_dims','inq_dimid levtot '//trim(fname))
-call nc_check(nf90_inquire_dimension(ncid, dimid, len=Nlevtot), &
+call nc_check(nf90_inquire_dimension(ncid, dimid, len=nlevtot), &
'get_sparse_dims','inquire_dimension levtot '//trim(fname))
call nc_check(nf90_inq_dimid(ncid, 'numrad', dimid), &
'get_sparse_dims','inq_dimid numrad '//trim(fname))
-call nc_check(nf90_inquire_dimension(ncid, dimid, len=Nnumrad), &
+call nc_check(nf90_inquire_dimension(ncid, dimid, len=nnumrad), &
'get_sparse_dims','inquire_dimension numrad '//trim(fname))
+! CLM4 does not have a multi-level canopy.
+! CLM4.5 has a multi-level canopy.
+istatus = nf90_inq_dimid(ncid, 'levcan', dimid)
+if (istatus == nf90_noerr) then
+ call nc_check(nf90_inquire_dimension(ncid, dimid, len=nlevcan), &
+ 'get_sparse_dims','inquire_dimension levcan '//trim(fname))
+endif
+
! levsno is presently required, but I can envision a domain/experiment that
! will not have snow levels. How this relates to variables dimensioned 'levtot'
! is unclear. For that reason, levsno is presently required.
istatus = nf90_inq_dimid(ncid, 'levsno', dimid)
if (istatus == nf90_noerr) then
- call nc_check(nf90_inquire_dimension(ncid, dimid, len=Nlevsno), &
+ call nc_check(nf90_inquire_dimension(ncid, dimid, len=nlevsno), &
'get_sparse_dims','inquire_dimension levsno '//trim(fname))
endif
@@ -3304,19 +3360,19 @@
istatus = nf90_inq_dimid(ncid, 'levsno1', dimid)
if (istatus == nf90_noerr) then
- call nc_check(nf90_inquire_dimension(ncid, dimid, len=Nlevsno1), &
+ call nc_check(nf90_inquire_dimension(ncid, dimid, len=nlevsno1), &
'get_sparse_dims','inquire_dimension levsno1 '//trim(fname))
endif
istatus = nf90_inq_dimid(ncid, 'rtmlon', dimid)
if (istatus == nf90_noerr) then
- call nc_check(nf90_inquire_dimension(ncid, dimid, len=Nrtmlon), &
+ call nc_check(nf90_inquire_dimension(ncid, dimid, len=nrtmlon), &
'get_sparse_dims','inquire_dimension rtmlon '//trim(fname))
endif
istatus = nf90_inq_dimid(ncid, 'rtmlat', dimid)
if (istatus == nf90_noerr) then
- call nc_check(nf90_inquire_dimension(ncid, dimid, len=Nrtmlat), &
+ call nc_check(nf90_inquire_dimension(ncid, dimid, len=nrtmlat), &
'get_sparse_dims','inquire_dimension rtmlat '//trim(fname))
endif
@@ -3329,32 +3385,34 @@
if ((debug > 7) .and. do_output()) then
write(logfileunit,*)
write(logfileunit,*)'get_sparse_dims output follows:'
- write(logfileunit,*)'Ngridcell = ',Ngridcell
- write(logfileunit,*)'Nlandunit = ',Nlandunit
- write(logfileunit,*)'Ncolumn = ',Ncolumn
- write(logfileunit,*)'Npft = ',Npft
- write(logfileunit,*)'Nlevgrnd = ',Nlevgrnd
- write(logfileunit,*)'Nlevlak = ',Nlevlak
- write(logfileunit,*)'Nlevsno = ',Nlevsno
- write(logfileunit,*)'Nlevsno1 = ',Nlevsno1
- write(logfileunit,*)'Nlevtot = ',Nlevtot
- write(logfileunit,*)'Nnumrad = ',Nnumrad
- write(logfileunit,*)'Nrtmlon = ',Nrtmlon
- write(logfileunit,*)'Nrtmlat = ',Nrtmlat
+ write(logfileunit,*)'ngridcell = ',ngridcell
+ write(logfileunit,*)'nlandunit = ',nlandunit
+ write(logfileunit,*)'ncolumn = ',ncolumn
+ write(logfileunit,*)'npft = ',npft
+ write(logfileunit,*)'nlevgrnd = ',nlevgrnd
+ write(logfileunit,*)'nlevlak = ',nlevlak
+ write(logfileunit,*)'nlevsno = ',nlevsno
+ write(logfileunit,*)'nlevsno1 = ',nlevsno1
+ write(logfileunit,*)'nlevtot = ',nlevtot
+ write(logfileunit,*)'nnumrad = ',nnumrad
+ write(logfileunit,*)'nlevcan = ',nlevcan
+ write(logfileunit,*)'nrtmlon = ',nrtmlon
+ write(logfileunit,*)'nrtmlat = ',nrtmlat
write( * ,*)
write( * ,*)'get_sparse_dims output follows:'
- write( * ,*)'Ngridcell = ',Ngridcell
- write( * ,*)'Nlandunit = ',Nlandunit
- write( * ,*)'Ncolumn = ',Ncolumn
- write( * ,*)'Npft = ',Npft
- write( * ,*)'Nlevgrnd = ',Nlevgrnd
- write( * ,*)'Nlevlak = ',Nlevlak
- write( * ,*)'Nlevsno = ',Nlevsno
- write( * ,*)'Nlevsno1 = ',Nlevsno1
- write( * ,*)'Nlevtot = ',Nlevtot
- write( * ,*)'Nnumrad = ',Nnumrad
- write( * ,*)'Nrtmlon = ',Nrtmlon
- write( * ,*)'Nrtmlat = ',Nrtmlat
+ write( * ,*)'ngridcell = ',ngridcell
+ write( * ,*)'nlandunit = ',nlandunit
+ write( * ,*)'ncolumn = ',ncolumn
+ write( * ,*)'npft = ',npft
+ write( * ,*)'nlevgrnd = ',nlevgrnd
+ write( * ,*)'nlevlak = ',nlevlak
+ write( * ,*)'nlevsno = ',nlevsno
+ write( * ,*)'nlevsno1 = ',nlevsno1
+ write( * ,*)'nlevtot = ',nlevtot
+ write( * ,*)'nnumrad = ',nnumrad
+ write( * ,*)'nlevcan = ',nlevcan
+ write( * ,*)'nrtmlon = ',nrtmlon
+ write( * ,*)'nrtmlat = ',nrtmlat
endif
end subroutine get_sparse_dims
@@ -3383,22 +3441,22 @@
! Make sure the variables are the right size ...
! by comparing agains the size of the variable ...
-if ( Ngridcell < 0 ) then
+if ( ngridcell < 0 ) then
write(string1,*)'Unable to read the number of gridcells.'
call error_handler(E_ERR,'get_sparse_geog',string1,source,revision,revdate)
endif
-if ( Nlandunit < 0 ) then
+if ( nlandunit < 0 ) then
write(string1,*)'Unable to read the number of land units.'
call error_handler(E_ERR,'get_sparse_geog',string1,source,revision,revdate)
endif
-if ( Ncolumn < 0 ) then
+if ( ncolumn < 0 ) then
write(string1,*)'Unable to read the number of columns.'
call error_handler(E_ERR,'get_sparse_geog',string1,source,revision,revdate)
endif
-if ( Npft < 0 ) then
+if ( npft < 0 ) then
write(string1,*)'Unable to read the number of pfts.'
call error_handler(E_ERR,'get_sparse_geog',string1,source,revision,revdate)
endif
@@ -3670,7 +3728,7 @@
integer, intent(out) :: ngood
character(len=*), dimension(:,:), intent(out) :: table
-integer :: nrows, ncols, i
+integer :: nrows, ncols, i, ivar
character(len=NF90_MAX_NAME) :: varname ! column 1
character(len=NF90_MAX_NAME) :: dartstr ! column 2
character(len=NF90_MAX_NAME) :: minvalstring ! column 3
@@ -3736,6 +3794,23 @@
endif
ngood = ngood + 1
+
+ ! Issue warning if DART kind is already in use by another variable.
+ ! The first variable specified with the DART kind is the one used
+ ! by the forward observation operators. All subsequent variables will
+ ! be updated, just not used for the direct forward operator.
+ EXISTLOOP : do ivar = 1,ngood-1
+
+ if (trim(table(i,VT_KINDINDX)) == trim(table(ivar,VT_KINDINDX)) ) then
+ write(string1,*)'.. WARNING: '//trim(table(i,VT_KINDINDX))//' already in DART from '//trim(table(ivar,VT_VARNAMEINDX))
+ write(string2,*)'WARNING: '//trim(table(ivar,VT_VARNAMEINDX))//' will be used for forward observation operator.'
+ write(string3,*)'WARNING: '//trim(table(i,VT_VARNAMEINDX))//' will still be updated, but not used.'
+
+ call error_handler(E_MSG,'parse_variable_table',string1,text2=string2,text3=string3)
+ endif
+
+ enddo EXISTLOOP
+
enddo MyLoop
if (ngood == nrows) then
@@ -4624,29 +4699,32 @@
-function findVarIndex(varstring, caller)
-character(len=*), intent(in) :: varstring
+function findKindIndex(kind_index, caller)
+integer, intent(in) :: kind_index
character(len=*), intent(in) :: caller
-integer :: findVarIndex
+integer :: findKindIndex
integer :: i
+character(len=paramname_length) :: kind_string
-findVarIndex = -1
+findKindIndex = -1
! Skip to the right variable
-VARTYPES : do i = 1,nfields
- if ( trim(progvar(i)%varname) == varstring) then
- findVarIndex = i
- exit VARTYPES
+KINDLOOP : do i = 1,nfields
+ if (progvar(i)%dart_kind == kind_index) then
+ findKindIndex = i
+ exit KINDLOOP
endif
-enddo VARTYPES
+enddo KINDLOOP
-if (findVarIndex < 1) then
- write(string1,*) trim(caller)//' cannot find "'//trim(varstring)//'" in list of DART state variables.'
- call error_handler(E_ERR,'findVarIndex',string1,source,revision,revdate)
+if (findKindIndex < 1) then
+ kind_string = get_raw_obs_kind_name( kind_index )
+ write(string1,*) trim(caller)//' cannot find "'//trim(kind_string)//'" in list of DART state variables.'
+ write(string2,*) trim(caller)//' looking for DART KIND (index) ',kind_index
+ call error_handler(E_ERR, 'findKindIndex', string1, source, revision, revdate, text2=string2)
endif
-end function findVarIndex
+end function findKindIndex
@@ -4668,7 +4746,7 @@
! Count up how many columns are in each gridcell
-do ij = 1,Ncolumn
+do ij = 1,ncolumn
ilon = cols1d_ixy(ij)
ilat = cols1d_jxy(ij)
gridCellInfo(ilon,ilat)%ncols = gridCellInfo(ilon,ilat)%ncols + 1
@@ -4676,7 +4754,7 @@
! Count up how many pfts are in each gridcell
-do ij = 1,Npft
+do ij = 1,npft
ilon = pfts1d_ixy(ij)
ilat = pfts1d_jxy(ij)
gridCellInfo(ilon,ilat)%npfts = gridCellInfo(ilon,ilat)%npfts + 1
@@ -4698,7 +4776,7 @@
! Fill the column pointer arrays
currenticol(:,:) = 0
-do ij = 1,Ncolumn
+do ij = 1,ncolumn
ilon = cols1d_ixy(ij)
ilat = cols1d_jxy(ij)
@@ -4720,7 +4798,7 @@
! Fill the pft pointer arrays
currentipft(:,:) = 0
-do ij = 1,Npft
+do ij = 1,npft
ilon = pfts1d_ixy(ij)
ilat = pfts1d_jxy(ij)
Modified: DART/trunk/models/clm/model_mod.html
===================================================================
--- DART/trunk/models/clm/model_mod.html 2015-11-06 21:01:09 UTC (rev 9014)
+++ DART/trunk/models/clm/model_mod.html 2015-11-06 22:20:29 UTC (rev 9015)
@@ -284,6 +284,58 @@
</TABLE>
</div>
+<P>The following are only meant to be examples - they are not scientifically validated.
+ Some of these that are UPDATED are probably diagnostic quantities, Some of these that
+ should be updated may be marked NO_COPY_BACK. There are multiple choices for some
+ DART kinds. This list is by no means complete.<P>
+<pre>
+ 'livecrootc', 'KIND_ROOT_CARBON', 'NA', 'NA', 'restart', 'UPDATE',
+ 'deadcrootc', 'KIND_ROOT_CARBON', 'NA', 'NA', 'restart', 'UPDATE',
+ 'livestemc', 'KIND_STEM_CARBON', 'NA', 'NA', 'restart', 'UPDATE',
+ 'deadstemc', 'KIND_STEM_CARBON', 'NA', 'NA', 'restart', 'UPDATE',
+ 'livecrootn', 'KIND_ROOT_NITROGEN', 'NA', 'NA', 'restart', 'UPDATE',
+ 'deadcrootn', 'KIND_ROOT_NITROGEN', 'NA', 'NA', 'restart', 'UPDATE',
+ 'livestemn', 'KIND_STEM_NITROGEN', 'NA', 'NA', 'restart', 'UPDATE',
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list