[Dart-dev] [7195] DART/trunk: The variable selection mechanism for what goes into the CLM DART vector can now come
nancy at ucar.edu
nancy at ucar.edu
Fri Oct 3 11:01:30 MDT 2014
Revision: 7195
Author: thoar
Date: 2014-10-03 11:01:29 -0600 (Fri, 03 Oct 2014)
Log Message:
-----------
The variable selection mechanism for what goes into the CLM DART vector can now come
from 3 sources ... the restart file, and two additional history files. It is now
demonstrated how to create a history file in the 'vector' format used by the restart files.
This should allow for more accurate forward observation operators. There is also now
a namelist-driven range-restriction capability for the variables that will be updated
in the CLM restart file. It is also now possible to make a run-time specification of
whether or not you want to update the variable in the restart file.
In order to remove confusion about this, a namelist variable had to change from
'clm_state_variables' to simply 'clm_variables' because some of the variables
in the DART vector are not technically 'state' in the strictest sense of the word.
The forward observation operators for brightness temperature may now explicitly grab
what they need from the DART state via the model_interpolate() function as opposed to
having to open,read,close a file (for each ensemble member) for every single observation.
Furthermore, the data can be at the same resolution as the restart file. All good.
Modified Paths:
--------------
DART/trunk/models/clm/clm_to_dart.f90
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/shell_scripts/CESM1_1_1_setup_hybrid
DART/trunk/models/clm/shell_scripts/CESM1_1_1_setup_pmo
DART/trunk/models/clm/shell_scripts/CESM_DART_config
DART/trunk/models/clm/shell_scripts/assimilate.csh
DART/trunk/models/clm/shell_scripts/perfect_model.csh
DART/trunk/models/clm/work/input.nml
DART/trunk/obs_kind/DEFAULT_obs_kind_mod.F90
Added Paths:
-----------
DART/trunk/models/clm/shell_scripts/CESM1_2_1_setup_hybrid
DART/trunk/models/clm/shell_scripts/CESM1_2_1_setup_pmo
-------------- next part --------------
Modified: DART/trunk/models/clm/clm_to_dart.f90
===================================================================
--- DART/trunk/models/clm/clm_to_dart.f90 2014-10-02 15:34:58 UTC (rev 7194)
+++ DART/trunk/models/clm/clm_to_dart.f90 2014-10-03 17:01:29 UTC (rev 7195)
@@ -24,8 +24,7 @@
use types_mod, only : r8
use utilities_mod, only : initialize_utilities, finalize_utilities, &
find_namelist_in_file, check_namelist_read
-use model_mod, only : get_model_size, restart_file_to_sv, &
- get_clm_restart_filename
+use model_mod, only : get_model_size, clm_to_dart_state_vector
use assim_model_mod, only : awrite_state_restart, open_restart_write, close_restart
use time_manager_mod, only : time_type, print_time, print_date
@@ -41,7 +40,7 @@
! namelist parameters with default values.
!-----------------------------------------------------------------------
-character(len=128) :: clm_to_dart_output_file = 'dart_ics'
+character(len=512) :: clm_to_dart_output_file = 'dart_ics'
namelist /clm_to_dart_nml/ clm_to_dart_output_file
@@ -52,7 +51,6 @@
integer :: io, iunit, x_size
type(time_type) :: model_time
real(r8), allocatable :: statevector(:)
-character(len=256) :: clm_restart_filename
!======================================================================
@@ -66,13 +64,6 @@
read(iunit, nml = clm_to_dart_nml, iostat = io)
call check_namelist_read(iunit, io, "clm_to_dart_nml") ! closes, too.
-call get_clm_restart_filename( clm_restart_filename )
-
-write(*,*)
-write(*,'(''clm_to_dart:converting clm restart file '',A, &
- &'' to DART file '',A)') &
- trim(clm_restart_filename), trim(clm_to_dart_output_file)
-
!----------------------------------------------------------------------
! get to work
!----------------------------------------------------------------------
@@ -80,7 +71,8 @@
x_size = get_model_size()
allocate(statevector(x_size))
-call restart_file_to_sv(clm_restart_filename, statevector, model_time)
+! Each variable specifies its own file of origin.
+call clm_to_dart_state_vector(statevector, model_time)
iunit = open_restart_write(clm_to_dart_output_file)
Modified: DART/trunk/models/clm/model_mod.f90
===================================================================
--- DART/trunk/models/clm/model_mod.f90 2014-10-02 15:34:58 UTC (rev 7194)
+++ DART/trunk/models/clm/model_mod.f90 2014-10-03 17:01:29 UTC (rev 7195)
@@ -28,11 +28,11 @@
use types_mod, only : r4, r8, SECPERDAY, MISSING_R8, &
MISSING_I, MISSING_R4, rad2deg, deg2rad, PI, &
obstypelength
-use time_manager_mod, only : set_time, get_time, set_date, get_date, &
+use time_manager_mod, only : time_type, set_time, get_time, set_date, get_date,&
print_time, print_date, set_calendar_type, &
- operator(*), operator(+), operator(-), &
- operator(>), operator(<), operator(/), &
- operator(/=), operator(<=), time_type
+ operator(*), operator(+), operator(-), &
+ operator(>), operator(<), operator(/), &
+ operator(/=), operator(<=), operator(==)
use location_mod, only : location_type, get_dist, query_location, &
get_close_maxdist_init, get_close_type, &
@@ -48,8 +48,8 @@
E_ERR, E_WARN, E_MSG, logfileunit, get_unit, &
nc_check, do_output, to_upper, &
find_namelist_in_file, check_namelist_read, &
- open_file, file_exist, find_textfile_dims, &
- file_to_text
+ file_exist, find_textfile_dims, file_to_text, &
+ open_file, close_file
use obs_kind_mod, only : KIND_SOIL_TEMPERATURE, &
KIND_SOIL_MOISTURE, &
@@ -96,13 +96,13 @@
! the interfaces here can be changed as appropriate.
public :: get_gridsize, &
- restart_file_to_sv, &
+ clm_to_dart_state_vector, &
sv_to_restart_file, &
get_clm_restart_filename, &
get_state_time, &
get_grid_vertval, &
compute_gridcell_value, &
- find_gridcell_Npft, &
+ gridcell_components, &
DART_get_var, &
get_model_time
@@ -138,11 +138,25 @@
integer, parameter :: LAKE = 3
+! Codes for restricting the range of a variable
+integer, parameter :: BOUNDED_NONE = 0 ! ... unlimited range
+integer, parameter :: BOUNDED_BELOW = 1 ! ... minimum, but no maximum
+integer, parameter :: BOUNDED_ABOVE = 2 ! ... maximum, but no minimum
+integer, parameter :: BOUNDED_BOTH = 3 ! ... minimum and maximum
+
integer :: nfields
integer, parameter :: max_state_variables = 40
-integer, parameter :: num_state_table_columns = 2
+integer, parameter :: num_state_table_columns = 6
character(len=obstypelength) :: variable_table(max_state_variables, num_state_table_columns)
+! Codes for interpreting the columns of the variable_table
+integer, parameter :: VT_VARNAMEINDX = 1 ! ... variable name
+integer, parameter :: VT_KINDINDX = 2 ! ... DART kind
+integer, parameter :: VT_MINVALINDX = 3 ! ... minimum value if any
+integer, parameter :: VT_MAXVALINDX = 4 ! ... maximum value if any
+integer, parameter :: VT_ORIGININDX = 5 ! ... file of origin
+integer, parameter :: VT_STATEINDX = 6 ! ... update (state) or not
+
! things which can/should be in the model_nml
integer :: assimilation_period_days = 0
@@ -153,18 +167,21 @@
character(len=32) :: calendar = 'Gregorian'
character(len=256) :: clm_restart_filename = 'clm_restart.nc'
character(len=256) :: clm_history_filename = 'clm_history.nc'
-character(len=obstypelength) :: clm_state_variables(max_state_variables*num_state_table_columns) = ' '
+character(len=256) :: clm_vector_history_filename = 'clm_vector_history.nc'
+character(len=obstypelength) :: clm_variables(max_state_variables*num_state_table_columns) = ' '
+
namelist /model_nml/ &
clm_restart_filename, &
clm_history_filename, &
+ clm_vector_history_filename, &
output_state_vector, &
assimilation_period_days, & ! for now, this is the timestep
assimilation_period_seconds, &
model_perturbation_amplitude,&
calendar, &
debug, &
- clm_state_variables
+ clm_variables
! Everything needed to describe a variable
@@ -175,23 +192,42 @@
character(len=NF90_MAX_NAME) :: units
character(len=obstypelength), dimension(NF90_MAX_VAR_DIMS) :: dimnames
integer, dimension(NF90_MAX_VAR_DIMS) :: dimlens
- integer :: numdims
- integer :: maxlevels
- integer :: xtype
- integer :: varsize ! prod(dimlens(1:numdims))
- integer :: index1 ! location in dart state vector of first occurrence
- integer :: indexN ! location in dart state vector of last occurrence
- integer :: dart_kind
- integer :: spvalINT, missingINT
+ integer :: numdims
+ integer :: maxlevels
+ integer :: xtype
+ integer :: varsize ! prod(dimlens(1:numdims))
+ integer :: index1 ! location in dart state vector of first occurrence
+ integer :: indexN ! location in dart state vector of last occurrence
+ integer :: dart_kind
+ integer :: rangeRestricted
+ real(r8) :: minvalue
+ real(r8) :: maxvalue
+ integer :: spvalINT, missingINT
real(r4) :: spvalR4, missingR4
real(r8) :: spvalR8, missingR8
logical :: has_fill_value ! intended for future use
logical :: has_missing_value ! intended for future use
character(len=paramname_length) :: kind_string
+ character(len=512) :: origin ! the file it came from
+ logical :: update
end type progvartype
type(progvartype), dimension(max_state_variables) :: progvar
+
+!----------------------------------------------------------------------
+! how many and which columns are in each gridcell
+!----------------------------------------------------------------------
+
+type gridcellcolumns ! given a gridcell, which columns contribute
+ private
+ integer :: ncols
+ integer, pointer, dimension(:) :: columnids
+ integer :: npfts
+ integer, pointer, dimension(:) :: pftids
+end type gridcellcolumns
+type(gridcellcolumns), allocatable, dimension(:,:), target :: gridCellInfo
+
!------------------------------------------------------------------------------
! Things that come from the CLM history file.
!
@@ -210,7 +246,7 @@
real(r8), allocatable :: LON(:)
real(r8), allocatable :: LAT(:)
real(r8), allocatable :: AREA1D(:), LANDFRAC1D(:) ! unstructured grid
-real(r8), allocatable :: AREA2D(:,:), LANDFRAC2D(:,:) ! 2D grid
+real(r8), allocatable :: AREA2D(:,:), LANDFRAC2D(:,:) ! 2D grid
logical :: unstructured = .false.
@@ -256,8 +292,6 @@
real(r8), allocatable, dimension(:,:):: zsno ! (column,levsno) ... snow layer midpoint
integer, allocatable, dimension(:) :: cols1d_ityplun ! columntype ... lake, forest, city ...
-
-
!------------------------------------------------------------------------------
! These are the metadata arrays that are the same size as the state vector.
@@ -291,6 +325,7 @@
INTERFACE DART_get_var
MODULE PROCEDURE get_var_1d
MODULE PROCEDURE get_var_2d
+ MODULE PROCEDURE get_var_3d
END INTERFACE
INTERFACE get_state_time
@@ -355,7 +390,7 @@
integer, intent(in) :: indx
type(location_type) :: location
-integer, optional, intent(out) :: var_type
+integer, OPTIONAL, intent(out) :: var_type
! Local variables
@@ -426,10 +461,8 @@
! Local variables - all the important ones have module scope
integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs
-character(len=NF90_MAX_NAME) :: varname
character(len=obstypelength) :: dimname
-character(len=paramname_length) :: kind_string
-integer :: ncid, VarID, dimlen, varsize
+integer :: ncid, TimeDimID, VarID, dimlen, varsize
integer :: iunit, io, ivar
integer :: i, j, xi, xj, index1, indexN, indx
integer :: ss, dd
@@ -469,7 +502,7 @@
call get_time(model_timestep,ss,dd) ! set_time() assures the seconds [0,86400)
write(string1,*)'assimilation period is ',dd,' days ',ss,' seconds'
-call error_handler(E_MSG,'static_init_model',string1,source,revision,revdate)
+call error_handler(E_MSG,'static_init_model',string1)
!---------------------------------------------------------------
! The CLM history file (h0?) has the 'superset' of information.
@@ -507,45 +540,44 @@
allocate(levtot(Nlevtot))
if (Nlevsno > 0) allocate(zsno(Nlevsno,Ncolumn))
-call get_sparse_geog(ncid, clm_restart_filename, 'open')
+call get_sparse_geog(ncid, clm_restart_filename, 'close')
!---------------------------------------------------------------
+! Generate list of columns in each gridcell
+
+allocate(gridCellInfo(nlon,nlat))
+call SetLocatorArrays()
+
+!---------------------------------------------------------------
! Compile the list of clm variables to use in the creation
-! of the DART state vector. Required to determine model_size.
-!
-! Verify all variables are in the clm restart file
-!
+! of the DART state vector.
+
+call parse_variable_table( clm_variables, nfields, variable_table )
+
! Compute the offsets into the state vector for the start of each
! variable type. Requires reading shapes from the clm restart file.
! Record the extent of the data type in the state vector.
-call verify_state_variables( clm_state_variables, ncid, clm_restart_filename, &
- nfields, variable_table )
-
index1 = 1;
indexN = 0;
do ivar = 1, nfields
- varname = trim(variable_table(ivar,1))
- kind_string = trim(variable_table(ivar,2))
- progvar(ivar)%varname = varname
- progvar(ivar)%kind_string = kind_string
- progvar(ivar)%dart_kind = get_raw_obs_kind_index( progvar(ivar)%kind_string )
- progvar(ivar)%dimlens = 0
- progvar(ivar)%dimnames = ' '
- progvar(ivar)%spvalINT = -9999 ! from CESM/CLM clm_varcon.F90
- progvar(ivar)%spvalR4 = 1.e36_r4 ! from CESM/CLM clm_varcon.F90
- progvar(ivar)%spvalR8 = 1.e36_r8 ! from CESM/CLM clm_varcon.F90
- progvar(ivar)%missingINT = MISSING_I
- progvar(ivar)%missingR4 = MISSING_R4
- progvar(ivar)%missingR8 = MISSING_R8
- progvar(ivar)%has_fill_value = .true.
- progvar(ivar)%has_missing_value = .true.
- progvar(ivar)%maxlevels = 0
+ ! convey the information in the variable_table to each progvar.
- string2 = trim(clm_restart_filename)//' '//trim(varname)
+ call SetVariableAttributes(ivar)
- call nc_check(nf90_inq_varid(ncid, trim(varname), VarID), &
+ ! Open the file for each variable and get dimensions, etc.
+
+ call nc_check(nf90_open(trim(progvar(ivar)%origin), NF90_NOWRITE, ncid), &
+ 'static_init_model','open '//trim(progvar(ivar)%origin))
+
+ ! File is not required to have a time dimension
+ io = nf90_inq_dimid(ncid, 'time', TimeDimID)
+ if (io /= NF90_NOERR) TimeDimID = MISSING_I
+
+ string2 = trim(progvar(ivar)%origin)//' '//trim(progvar(ivar)%varname)
+
+ call nc_check(nf90_inq_varid(ncid, trim(progvar(ivar)%varname), VarID), &
'static_init_model', 'inq_varid '//trim(string2))
call nc_check(nf90_inquire_variable(ncid, VarID, dimids=dimIDs, &
@@ -559,7 +591,7 @@
call nc_check( nf90_get_att(ncid, VarID, 'long_name' , progvar(ivar)%long_name), &
'static_init_model', 'get_att long_name '//trim(string2))
else
- progvar(ivar)%long_name = varname
+ progvar(ivar)%long_name = progvar(ivar)%varname
endif
if( nf90_inquire_attribute( ncid, VarID, 'units') == NF90_NOERR ) then
@@ -611,6 +643,10 @@
write(string1,'(''inquire dimension'',i2,A)') i,trim(string2)
call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), name=dimname, len=dimlen), &
'static_init_model', string1)
+
+ ! Only reserve space for a single time slice
+ if (dimIDs(i) == TimeDimID) dimlen = 1
+
progvar(ivar)%dimlens( i) = dimlen
progvar(ivar)%dimnames(i) = dimname
varsize = varsize * dimlen
@@ -625,10 +661,12 @@
if ((debug > 0) .and. do_output()) then
write(logfileunit,*)
write(logfileunit,*) trim(progvar(ivar)%varname),' variable number ',ivar
+ write(logfileunit,*) ' filename ',trim(progvar(ivar)%origin)
+ write(logfileunit,*) ' update ',progvar(ivar)%update
write(logfileunit,*) ' long_name ',trim(progvar(ivar)%long_name)
write(logfileunit,*) ' units ',trim(progvar(ivar)%units)
write(logfileunit,*) ' xtype ',progvar(ivar)%xtype
- write(logfileunit,*) ' dimnames ',progvar(ivar)%dimnames(1:progvar(ivar)%numdims)
+ write(logfileunit,*) ' dimnames ',(/ (trim(progvar(ivar)%dimnames(i))//' ', i=1,progvar(ivar)%numdims ) /)
write(logfileunit,*) ' dimlens ',progvar(ivar)%dimlens( 1:progvar(ivar)%numdims)
write(logfileunit,*) ' numdims ',progvar(ivar)%numdims
write(logfileunit,*) ' varsize ',progvar(ivar)%varsize
@@ -644,13 +682,18 @@
write(logfileunit,*) ' missingR8 ',progvar(ivar)%missingR8
write(logfileunit,*) ' has_fill_value ',progvar(ivar)%has_fill_value
write(logfileunit,*) ' has_missing_value ',progvar(ivar)%has_missing_value
+ write(logfileunit,*)' rangeRestricted ',progvar(ivar)%rangeRestricted
+ write(logfileunit,*)' minvalue ',progvar(ivar)%minvalue
+ write(logfileunit,*)' maxvalue ',progvar(ivar)%maxvalue
write( * ,*)
write( * ,*) trim(progvar(ivar)%varname),' variable number ',ivar
+ write( * ,*) ' filename ',trim(progvar(ivar)%origin)
+ write( * ,*) ' update ',progvar(ivar)%update
write( * ,*) ' long_name ',trim(progvar(ivar)%long_name)
write( * ,*) ' units ',trim(progvar(ivar)%units)
write( * ,*) ' xtype ',progvar(ivar)%xtype
- write( * ,*) ' dimnames ',progvar(ivar)%dimnames(1:progvar(ivar)%numdims)
+ write( * ,*) ' dimnames ',(/ (trim(progvar(ivar)%dimnames(i))//' ', i=1,progvar(ivar)%numdims ) /)
write( * ,*) ' dimlens ',progvar(ivar)%dimlens( 1:progvar(ivar)%numdims)
write( * ,*) ' numdims ',progvar(ivar)%numdims
write( * ,*) ' varsize ',progvar(ivar)%varsize
@@ -666,13 +709,16 @@
write( * ,*) ' missingR8 ',progvar(ivar)%missingR8
write( * ,*) ' has_fill_value ',progvar(ivar)%has_fill_value
write( * ,*) ' has_missing_value ',progvar(ivar)%has_missing_value
+ write( * ,*)' rangeRestricted ',progvar(ivar)%rangeRestricted
+ write( * ,*)' minvalue ',progvar(ivar)%minvalue
+ write( * ,*)' maxvalue ',progvar(ivar)%maxvalue
endif
+ call nc_check(nf90_close(ncid),'static_init_model','close '//trim(string2))
+ ncid = 0
+
enddo
-call nc_check(nf90_close(ncid),'static_init_model','close '//trim(clm_restart_filename))
-ncid = 0
-
model_size = progvar(nfields)%indexN
if ((debug > 0) .and. do_output()) then
@@ -704,6 +750,10 @@
indx = progvar(ivar)%index1
+ ! 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'
+
if (progvar(ivar)%numdims == 1) then
if ((debug > 8) .and. do_output()) then
@@ -779,13 +829,17 @@
CASE DEFAULT
write(string1,*)'(1d) unknown Dimension name '//trim(progvar(ivar)%dimnames(1))// &
- ' while trying to create metadata arrays.'
+ ' while trying to create metadata for '//trim(progvar(ivar)%varname)
call error_handler(E_ERR,'static_init_model', string1, source, revision, revdate)
END SELECT
elseif (progvar(ivar)%numdims == 2) then
+ ! restart file variables may have 2 dimensions, one of which is layers.
+ ! vector_history files may have 2 dimensions, one of which is time
+ ! history file variables always have 3 dimension [time,lat,lon]
+
! In the ncdump output, dimension 2 is the leftmost dimension.
! Only dimension 2 matters for the weights.
@@ -885,15 +939,107 @@
enddo
enddo
+ CASE ("time")
+
+ ! The vector history files can have things '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
+ if (unstructured) then
+ lonixy( indx) = xi
+ latjxy( indx) = xi
+ landarea(indx) = AREA1D(xi) * LANDFRAC1D(xi) * pfts1d_wtxy(i)
+ else
+ lonixy( indx) = xi
+ latjxy( indx) = xj
+ landarea(indx) = AREA2D(xi,xj) * LANDFRAC2D(xi,xj) * pfts1d_wtxy(i)
+ endif
+ indx = indx + 1
+ enddo
+
+ CASE ("column")
+ do i = 1, progvar(ivar)%dimlens(1)
+ xi = cols1d_ixy(i)
+ xj = cols1d_jxy(i) ! always unity if unstructured
+ if (unstructured) then
+ lonixy( indx) = xi
+ latjxy( indx) = xi
+ landarea(indx) = AREA1D(xi) * LANDFRAC1D(xi) * cols1d_wtxy(i)
+ else
+ lonixy( indx) = xi
+ latjxy( indx) = xj
+ landarea(indx) = AREA2D(xi,xj) * LANDFRAC2D(xi,xj) * cols1d_wtxy(i)
+ 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)
+ call error_handler(E_ERR,'static_init_model', string1, source, revision, revdate)
+ END SELECT
+
CASE DEFAULT
write(string1,*)'(2d) unknown Dimension name '//trim(progvar(ivar)%dimnames(2))// &
- ' while trying to create metadata arrays.'
+ ' while trying to create metadata for '//trim(progvar(ivar)%varname)
call error_handler(E_ERR,'static_init_model', string1, source, revision, revdate)
END SELECT
+
+ elseif (progvar(ivar)%numdims == 3) then
+
+ ! restart file variables never have 3 dimensions
+ ! vector_history variables may have 3 dimensions [time, lat, lon]
+ ! history file variables always have 3 dimensions [time, lat, lon]
+ ! exception is float H2OSOI(time, levgrnd, lat, lon) ... but we
+ ! have access to restart file h2osoi_[liq,ice]
+
+ if ((debug > 8) .and. do_output()) then
+ write(*,*)
+ write(*,*)'variable ',trim(progvar(ivar)%varname)
+ write(*,*)'dimension 1 (i) ',progvar(ivar)%dimnames(1),progvar(ivar)%dimlens(1)
+ write(*,*)'dimension 2 (j) ',progvar(ivar)%dimnames(2),progvar(ivar)%dimlens(2)
+ write(*,*)'dimension 3 (k) ',progvar(ivar)%dimnames(3),progvar(ivar)%dimlens(3)
+ endif
+
+ ! Remember the order is reversed from ncdump to fortran
+ ! The messages are in ncdump-order, checking is fortran-order
+ if ((trim(progvar(ivar)%dimnames(1)) .ne. 'lon') .or. &
+ (trim(progvar(ivar)%dimnames(2)) .ne. 'lat' ) .or. &
+ (trim(progvar(ivar)%dimnames(3)) .ne. 'time' )) then
+ write(string1,*)'3D variables must be [time,lat,lon] (as reported by ncdump)'
+ write(string2,*)trim(progvar(ivar)%varname),' is ', &
+ trim(progvar(ivar)%dimnames(3)), ' ', &
+ trim(progvar(ivar)%dimnames(2)), ' ', &
+ trim(progvar(ivar)%dimnames(1))
+ call error_handler(E_ERR,'static_init_model', string1, &
+ source, revision, revdate, text2=string2)
+ endif
+
+ ! The get_var_3d() routine ensures there is only 1 timestep.
+ ! So there is no need to loop over dimlens(3)
+
+ do j = 1, progvar(ivar)%dimlens(2) ! time-invariant
+ do i = 1, progvar(ivar)%dimlens(1) ! time-invariant
+ lonixy( indx) = i
+ latjxy( indx) = j
+ landarea(indx) = AREA2D(i,j) * LANDFRAC2D(i,j)
+ indx = indx + 1
+ enddo
+ enddo
+
else
+ ! Cannot support 4D variables yet.
+ ! These only occurr in 2D history files for variables with levels (and time)
+ ! i.e. H2OSOI(time, levgrnd, lat, lon)
+ !
+ ! You can get the same information from the restart file, usually.
+
write(string1,*)'variables of rank ',progvar(ivar)%numdims,' are unsupported.'
write(string2,*)trim(progvar(ivar)%varname),' is dimensioned ',&
progvar(ivar)%dimlens(1:progvar(ivar)%numdims)
@@ -1413,7 +1559,7 @@
call nc_check(nf90_put_var(ncFileID, VarID, LEVGRND ), &
'nc_write_model_atts', 'levgrnd put_var '//trim(filename))
- ! AREA can be 1D or 2D
+ ! AREA can be 1D or 2D
call nc_check(nf90_inq_varid(ncFileID, 'area', VarID), &
'nc_write_model_atts', 'put_var area '//trim(filename))
if (unstructured) then
@@ -1425,7 +1571,7 @@
endif
- ! LANDFRAC can be 1D or 2D
+ ! LANDFRAC can be 1D or 2D
call nc_check(nf90_inq_varid(ncFileID, 'landfrac', VarID), &
'nc_write_model_atts', 'put_var landfrac '//trim(filename))
if (unstructured) then
@@ -1517,7 +1663,7 @@
integer, intent(in) :: timeindex
integer :: ierr ! return value of function
-integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs, mystart, mycount
+integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs, ncstart, nccount
character(len=NF90_MAX_NAME) :: varname
integer :: i, ivar, VarID, ncNdims, dimlen
integer :: TimeDimID, CopyDimID
@@ -1578,14 +1724,16 @@
call nc_check(nf90_inquire_variable(ncFileID,VarID,dimids=dimIDs,ndims=ncNdims), &
'nc_write_model_vars', 'inquire '//trim(string2))
- mystart = 1 ! These are arrays, actually
- mycount = 1
+ ncstart = 1 ! These are arrays, actually
+ nccount = 1
DimCheck : do i = 1,progvar(ivar)%numdims
write(string1,'(a,i2,A)') 'inquire dimension ',i,trim(string2)
call nc_check(nf90_inquire_dimension(ncFileID, dimIDs(i), len=dimlen), &
'nc_write_model_vars', trim(string1))
+ if (progvar(ivar)%dimnames(i) == 'time') cycle DimCheck
+
if ( dimlen /= progvar(ivar)%dimlens(i) ) then
write(string1,*)trim(string2),' dim/dimlen ',i,dimlen, &
' not ',progvar(ivar)%dimlens(i)
@@ -1594,22 +1742,27 @@
source, revision, revdate, text2=trim(string2))
endif
- mycount(i) = dimlen
+ nccount(i) = dimlen
enddo DimCheck
- where(dimIDs == CopyDimID) mystart = copyindex
- where(dimIDs == CopyDimID) mycount = 1
- where(dimIDs == TimeDimID) mystart = timeindex
- where(dimIDs == TimeDimID) mycount = 1
+ where(dimIDs == CopyDimID) ncstart = copyindex
+ where(dimIDs == CopyDimID) nccount = 1
+ where(dimIDs == TimeDimID) ncstart = timeindex
+ where(dimIDs == TimeDimID) nccount = 1
- if ((debug > 7) .and. do_output()) then
- write(*,*)'nc_write_model_vars '//trim(varname)//' start is ',mystart(1:ncNdims)
- write(*,*)'nc_write_model_vars '//trim(varname)//' count is ',mycount(1:ncNdims)
+ if ((debug > 0) .and. do_output()) then
+ write(*,*)'nc_write_model_vars '//trim(varname)//' start is ',ncstart(1:ncNdims)
+ write(*,*)'nc_write_model_vars '//trim(varname)//' count is ',nccount(1:ncNdims)
endif
- if ( progvar(ivar)%numdims == 1 ) then
+ ! Since 'time' is a singleton dimension, we can use the same logic
+ ! as if it the variable had one less dimension.
+ if ( (progvar(ivar)%numdims == 1) .or. &
+ ((progvar(ivar)%numdims == 2) .and. &
+ (progvar(ivar)%dimnames(2) == 'time')) )then
+
if ( ncNdims /= 3 ) then
write(string1,*)trim(varname),' no room for copy,time dimensions.'
write(string2,*)'netcdf file should have 3 dimensions, has ',ncNdims
@@ -1620,11 +1773,13 @@
allocate(data_1d_array( progvar(ivar)%dimlens(1) ) )
call vector_to_prog_var(state_vec, ivar, data_1d_array)
call nc_check(nf90_put_var(ncFileID, VarID, data_1d_array, &
- start = mystart(1:ncNdims), count=mycount(1:ncNdims)), &
+ start = ncstart(1:ncNdims), count=nccount(1:ncNdims)), &
'nc_write_model_vars', 'put_var '//trim(string2))
deallocate(data_1d_array)
- elseif ( progvar(ivar)%numdims == 2 ) then
+ elseif ( (progvar(ivar)%numdims == 2) .or. &
+ ((progvar(ivar)%numdims == 3) .and. &
+ (progvar(ivar)%dimnames(3) == 'time')) )then
if ( ncNdims /= 4 ) then
write(string1,*)trim(varname),' no room for copy,time dimensions.'
@@ -1637,7 +1792,7 @@
progvar(ivar)%dimlens(2) ))
call vector_to_prog_var(state_vec, ivar, data_2d_array)
call nc_check(nf90_put_var(ncFileID, VarID, data_2d_array, &
- start = mystart(1:ncNdims), count=mycount(1:ncNdims)), &
+ start = ncstart(1:ncNdims), count=nccount(1:ncNdims)), &
'nc_write_model_vars', 'put_var '//trim(string2))
deallocate(data_2d_array)
@@ -1753,25 +1908,25 @@
-subroutine restart_file_to_sv(filename, state_vector, restart_time)
+subroutine clm_to_dart_state_vector(state_vector, restart_time)
!------------------------------------------------------------------
! Reads the current time and state variables from a clm restart
! file and packs them into a dart state vector. This better happen
! in the same fashion as the metadata arrays are built.
-character(len=*), intent(in) :: filename
-real(r8), intent(inout) :: state_vector(:)
-type(time_type), intent(out) :: restart_time
+real(r8), intent(out) :: state_vector(:)
+type(time_type), intent(out) :: restart_time
! temp space to hold data while we are reading it
integer :: i, j, ni, nj, ivar, indx, numsnowlevels
integer, allocatable, dimension(:) :: snlsno
real(r8), allocatable, dimension(:) :: data_1d_array
real(r8), allocatable, dimension(:,:) :: data_2d_array
+real(r8), allocatable, dimension(:,:,:) :: data_3d_array
integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs
character(len=NF90_MAX_NAME) :: varname
-integer :: VarID, ncNdims, dimlen
+integer :: io, TimeDimID, VarID, ncNdims, dimlen
integer :: ncid
character(len=256) :: myerrorstring
@@ -1779,21 +1934,6 @@
state_vector = MISSING_R8
-! Check that the input file exists ...
-
-if ( .not. file_exist(filename) ) then
- write(string1,*) 'cannot open file ', trim(filename),' for reading.'
- call error_handler(E_ERR,'restart_file_to_sv',string1,source,revision,revdate)
-endif
-
-call nc_check(nf90_open(trim(filename), NF90_NOWRITE, ncid), &
- 'restart_file_to_sv','open '//trim(filename))
-
-restart_time = get_state_time(ncid)
-
-if (do_output()) call print_time(restart_time,'time in restart file '//trim(filename))
-if (do_output()) call print_date(restart_time,'date in restart file '//trim(filename))
-
! Must check anything with a dimension of 'levtot' or 'levsno' and manually
! set the values to DART missing. If only it were that easy ...
!
@@ -1808,32 +1948,52 @@
! The snow over lakes is wholly contained in the bulk formulation variables
! as opposed to the snow layer variables.
-! read number of snow layers
+! Some things must be gotten explicitly from the restart file - ONCE.
+! number of snow layers
+! time of restart file
allocate(snlsno(Ncolumn))
-call nc_check(nf90_inq_varid(ncid,'SNLSNO', VarID), 'restart_file_to_sv', 'inq_varid SNLSNO')
-call nc_check(nf90_get_var(ncid, VarID, snlsno), 'restart_file_to_sv', 'get_var SNLSNO')
+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), &
+ 'clm_to_dart_state_vector', 'inq_varid SNLSNO'//clm_restart_filename)
+call nc_check(nf90_get_var(ncid, VarID, snlsno), &
+ 'clm_to_dart_state_vector', 'get_var SNLSNO'//clm_restart_filename)
+restart_time = get_state_time(ncid)
+
+if (do_output()) call print_time(restart_time,'time in restart file '//clm_restart_filename)
+if (do_output()) call print_date(restart_time,'date in restart file '//clm_restart_filename)
+
+call nc_check(nf90_close(ncid),'clm_to_dart_state_vector','close '//clm_restart_filename)
+
! Start counting and filling the state vector one item at a time,
! repacking the Nd arrays into a single 1d list of numbers.
do ivar=1, nfields
varname = trim(progvar(ivar)%varname)
- myerrorstring = trim(filename)//' '//trim(varname)
+ myerrorstring = trim(progvar(ivar)%origin)//' '//trim(progvar(ivar)%varname)
+ call nc_check(nf90_open(trim(progvar(ivar)%origin), NF90_NOWRITE, ncid), &
+ 'clm_to_dart_state_vector','open '//trim(myerrorstring))
+
+ ! File is not required to have a time dimension
+ io = nf90_inq_dimid(ncid, 'time', TimeDimID)
+ if (io /= NF90_NOERR) TimeDimID = MISSING_I
+
call nc_check(nf90_inq_varid(ncid, varname, VarID), &
- 'restart_file_to_sv', 'inq_varid '//trim(myerrorstring))
+ 'clm_to_dart_state_vector', 'inq_varid '//trim(myerrorstring))
call nc_check(nf90_inquire_variable(ncid,VarID,dimids=dimIDs,ndims=ncNdims), &
- 'restart_file_to_sv', 'inquire '//trim(myerrorstring))
+ 'clm_to_dart_state_vector', 'inquire '//trim(myerrorstring))
! Check the rank of the variable
if ( ncNdims /= progvar(ivar)%numdims ) then
write(string1, *) 'netCDF rank of '//trim(varname)//' does not match derived type knowledge'
write(string2, *) 'netCDF rank is ',ncNdims,' expected ',progvar(ivar)%numdims
- call error_handler(E_ERR,'restart_file_to_sv', string1, &
+ call error_handler(E_ERR,'clm_to_dart_state_vector', string1, &
source,revision,revdate,text2=string2)
endif
@@ -1843,11 +2003,18 @@
write(string1,'(''inquire dimension'',i2,A)') i,trim(myerrorstring)
call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), len=dimlen), &
- 'restart_file_to_sv', string1)
+ 'clm_to_dart_state_vector', string1)
+ ! Time dimension will be unity 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.
+
+ if ( dimIDs(i) == TimeDimID ) dimlen = 1
+
if ( dimlen /= progvar(ivar)%dimlens(i) ) then
- write(string1,*) trim(myerrorstring),' dim/dimlen ',i,dimlen,' not ',progvar(ivar)%dimlens(i)
- call error_handler(E_ERR,'restart_file_to_sv',string1,source,revision,revdate)
+ write(string1,*) trim(myerrorstring),' dim/dimlen ',i,dimlen, &
+ ' not ',progvar(ivar)%dimlens(i)
+ call error_handler(E_ERR,'clm_to_dart_state_vector',string1,source,revision,revdate)
endif
enddo
@@ -1946,28 +2113,71 @@
enddo
deallocate(data_2d_array)
+ elseif (ncNdims == 3) then
+
+ ! restart file variables never have 3 dimensions
+ ! vector_history variables may have 3 dimensions [time, lat, lon]
+ ! history file variables always have 3 dimensions [time, lat, lon]
+ ! exception is float H2OSOI(time, levgrnd, lat, lon) ... but we
+ ! have access to restart file h2osoi_[liq,ice]
+
+ if ( (trim(progvar(ivar)%dimnames(1)) == 'lon') .and. &
+ (trim(progvar(ivar)%dimnames(2)) == 'lat') .and. &
+ (trim(progvar(ivar)%dimnames(3)) == 'time') ) then
+
+ ni = progvar(ivar)%dimlens(1)
+ nj = progvar(ivar)%dimlens(2)
+ ! nk = progvar(ivar)%dimlens(3) not needed ... time is always a singleton
+
+ allocate(data_3d_array(ni, nj, 1))
+ call DART_get_var(ncid, varname, data_3d_array)
+
+ ! In the CLM history files, the _missing_value_ flag seems to be
+ ! applied correctly for PBOT, TBOT ... so there is no need for the
+ ! extra processing that is present in the previous loops.
+
+ do j = 1, nj
+ do i = 1, ni
+ state_vector(indx) = data_3d_array(i, j, 1)
+ indx = indx + 1
+ enddo
+ enddo
+ deallocate(data_3d_array)
+ else
+
+ write(string1, *) '3D variable unexpected shape -- only support nlon, nlat, time(=1)'
+ write(string2, *) 'variable [',trim(progvar(ivar)%varname),']'
+ write(string3, *) 'file [',trim(progvar(ivar)%origin),']'
+ call error_handler(E_ERR,'clm_to_dart_state_vector', string1, &
+ source, revision, revdate, text2=string2, text3=string3)
+
+ endif
+
else
+
write(string1, *) 'no support for data array of dimension ', ncNdims
- call error_handler(E_ERR,'restart_file_to_sv', string1, &
- source,revision,revdate)
+ write(string2, *) 'variable [',trim(progvar(ivar)%varname),']'
+ write(string3, *) 'file [',trim(progvar(ivar)%origin),']'
+ call error_handler(E_ERR,'clm_to_dart_state_vector', string1, &
+ source, revision, revdate, text2=string2, text3=string3)
endif
indx = indx - 1
if ( indx /= progvar(ivar)%indexN ) then
write(string1, *)'Variable '//trim(varname)//' filled wrong.'
write(string2, *)'Should have ended at ',progvar(ivar)%indexN,' actually ended at ',indx
- call error_handler(E_ERR,'restart_file_to_sv', string1, &
+ call error_handler(E_ERR,'clm_to_dart_state_vector', string1, &
source,revision,revdate,text2=string2)
endif
+ call nc_check(nf90_close(ncid),'clm_to_dart_state_vector','close '//progvar(ivar)%origin)
+ ncid = 0
+
enddo
-call nc_check(nf90_close(ncid),'restart_file_to_sv','close '//trim(filename))
-ncid = 0
-
deallocate(snlsno)
-end subroutine restart_file_to_sv
+end subroutine clm_to_dart_state_vector
@@ -2038,10 +2248,10 @@
varname = trim(progvar(ivar)%varname)
string2 = trim(filename)//' '//trim(varname)
- if (trim(varname) == 'frac_sno') then
- ! frac_sno (snow cover fraction) is a diagnosed field.
- ! Simply update SNOWDP (snow depth) and CLM recalculates frac_sno.
- ! Actually updating frac_sno causes restart problems for CLM.
+ if ( .not. progvar(ivar)%update ) then
+ write(string1,*)'intentionally not updating '//trim(string2)
+ write(string3,*)'as per namelist control in model_nml:clm_variables'
+ call error_handler(E_MSG, 'sv_to_restart_file', string1, text2=string3)
cycle UPDATE
endif
@@ -2148,7 +2358,7 @@
! Local storage
-real(r8), dimension(3) :: loc_array
+real(r8), dimension(LocationDims) :: loc_array
real(r8) :: llon, llat, lheight
real(r8) :: interp_val_2
integer :: istatus_2
@@ -2269,7 +2479,7 @@
loc_lat = loc(2)
! determine the portion of interest of the state vector
-ivar = findVarindex(varstring, 'compute_gridcell_value')
+ivar = findVarIndex(varstring, '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
@@ -2296,6 +2506,9 @@
! If there is no vertical component, the problem is greatly simplified.
! Simply area-weight an average of all pieces in the grid cell.
+! FIXME ... this is the loop that can exploit the knowledge of what
+! columnids or pftids are needed for any particular gridcell.
+! gridCellInfo%pftids, gridCellInfo%columnids
counter = 0
total = 0.0_r8 ! temp storage for state vector
@@ -2329,12 +2542,12 @@
interp_val = total/total_area
istatus = 0
else
- if ((debug > 0) .and. do_output()) then
+ if ((debug > 4) .and. do_output()) then
write(string1, *)'Variable '//trim(varstring)//' 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, &
- source, revision, revdate, text2=string2,text3=string3)
+ text2=string2,text3=string3)
endif
endif
@@ -2342,8 +2555,7 @@
if ((debug > 5) .and. do_output()) then
write(string1,*)'counter, total, total_area', counter, total, total_area
write(string2,*)'interp_val, istatus', interp_val, istatus
- call error_handler(E_MSG,'compute_gridcell_value', string1, &
- source, revision, revdate, text2=string2)
+ call error_handler(E_MSG,'compute_gridcell_value', string1, text2=string2)
endif
end subroutine compute_gridcell_value
@@ -2405,7 +2617,7 @@
endif
! determine the portion of interest of the state vector
-ivar = findVarindex(varstring, 'get_grid_vertval')
+ivar = findVarIndex(varstring, '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
@@ -2483,7 +2695,7 @@
write(string2, *)'at gridcell lon/lat = (',gridloni,',',gridlatj,')'
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list