diff --git a/.gitignore b/.gitignore index b580da07fc..cadf4a4fd4 100644 --- a/.gitignore +++ b/.gitignore @@ -208,6 +208,7 @@ test_beta_dist test_kde_dist test_window test_force_bounds +test_parse_variables # Directories to NOT IGNORE ... same as executable names # as far as I know, these must be listed after the executables diff --git a/CHANGELOG.rst b/CHANGELOG.rst index a29103d032..7b8761d0f8 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -22,6 +22,14 @@ individual files. The changes are now listed with the most recent at the top. +**February 4 2025 :: Generic model_mod subroutine parse_variables. Tag v11.10.3** + +- Creates two generalized subroutines that convert the table of state variables that is + read in from the &model_nml to a state_var_type: parse_variables and parse_variables_clamp +- Alternate versions for this subroutine were replaced with parse_variables in the + models MOM6, wrf_hydro, aether_lat-lon, cam-fv, cam-se, POP, and cice +- New dev test test_parse_variables added + **February 3 2025 :: Inflation documentation. Tag v11.10.2** - Improved inflation documentation diff --git a/assimilation_code/modules/io/state_structure_mod.f90 b/assimilation_code/modules/io/state_structure_mod.f90 index 97de1ea97e..fb3240ac49 100644 --- a/assimilation_code/modules/io/state_structure_mod.f90 +++ b/assimilation_code/modules/io/state_structure_mod.f90 @@ -18,11 +18,15 @@ module state_structure_mod !> The add_domain() call adds a 'domain' to the state. This may be a component in !> the case of XCESM or another coupled model. !> -!> There are three ways to add a domain (these are overloaded as add_domain): +!> There are four ways to add a domain (these are overloaded as add_domain): !> * add_domain_blank. This takes model size as an argument. !> !> * add_domain_from_file. This takes a netcdf file and a list of variables !> +!> * add_domain_from_state_type. This takes a netcdf file and a state_var_type, +!> which includes nvars, netcdf variable names, qtys (kinds), +!> clamp_values (optional), and updates +!> !> * add_domain_from_spec. This makes a skeleton structure for a domain. Dimensions !> for each variable must be added using add_dimension_to_variable(). This is intended !> to be used to create netcdf output for models like bgrid_solo that are spun up. @@ -68,6 +72,8 @@ module state_structure_mod use sort_mod, only : index_sort +use default_model_mod, only : state_var_type + use netcdf implicit none @@ -240,8 +246,8 @@ module state_structure_mod character(len=256) :: info_file = 'NULL' ! string identifying the manner in which the domain was created - ! 'blank', 'file', or 'spec' - character(len=6) :: method = 'none' + ! 'blank', 'file', 'state_type', or 'spec' + character(len=11) :: method = 'none' end type domain_type @@ -281,6 +287,7 @@ module state_structure_mod module procedure add_domain_blank module procedure add_domain_from_file module procedure add_domain_from_spec + module procedure add_domain_from_state_type end interface interface get_index_start @@ -313,8 +320,6 @@ module state_structure_mod !> into the state_strucutre. !> !> Returns a dom_id that can be used to harvest information of a particular domain -!> -!> Does this need to be a function or a subroutine? function add_domain_from_file(info_file, num_vars, var_names, kind_list, clamp_vals, update_list) result(dom_id) @@ -365,6 +370,58 @@ function add_domain_from_file(info_file, num_vars, var_names, kind_list, clamp_v end function add_domain_from_file +!------------------------------------------------------------------------------- +!> Given an info_file, reads in a state_var_type including nvars, netcdf +!> variable names, qtys (kinds), clamp values (optional), and updates into the +!> state_structure +!> +!> Returns a dom_id that can be used to harvest information of a particular domain + + +function add_domain_from_state_type(info_file, vars) result(dom_id) + +character(len=*), intent(in) :: info_file +type(state_var_type), intent(in) :: vars +integer :: dom_id + +integer :: ivar + +! add to domains +call assert_below_max_num_domains('add_domain_from_state_type') +state%num_domains = state%num_domains + 1 +!>@todo dom_id should be a handle. +dom_id = state%num_domains + +! save information about the information file +state%domain(dom_id)%info_file = info_file +state%domain(dom_id)%method = 'state_type' + +! set number of variables in this domain +state%domain(dom_id)%num_variables = vars%nvars + +! load up the variable names +allocate(state%domain(dom_id)%variable(vars%nvars)) + +do ivar = 1, vars%nvars + state%domain(dom_id)%variable(ivar)%varname = vars%netcdf_var_names(ivar) +enddo + +! load up variable id's and sizes +call load_state_variable_info(state%domain(dom_id),dom_id) + +! load up the domain unique dimension info +call load_unique_dim_info(dom_id) + +! load up any cf-conventions if they exist +call load_common_cf_conventions(state%domain(dom_id)) + +call set_dart_kinds(dom_id, vars%nvars, vars%qtys) +if (allocated(vars%clamp_values)) call set_clamping(dom_id, vars%nvars, vars%clamp_values) +call set_update_list(dom_id, vars%nvars, vars%updates) + +end function add_domain_from_state_type + + !------------------------------------------------------------------------------- !> Defines a skeleton structure for the state structure. Dimension can be !> added to variables with add_dimension_to_variable. diff --git a/conf.py b/conf.py index f7ce20b498..67e3a861ce 100644 --- a/conf.py +++ b/conf.py @@ -21,7 +21,7 @@ author = 'Data Assimilation Research Section' # The full version, including alpha/beta/rc tags -release = '11.10.2' +release = '11.10.3' root_doc = 'index' # -- General configuration --------------------------------------------------- diff --git a/developer_tests/namelist/test_parse_variables.f90 b/developer_tests/namelist/test_parse_variables.f90 new file mode 100644 index 0000000000..16e75ceb42 --- /dev/null +++ b/developer_tests/namelist/test_parse_variables.f90 @@ -0,0 +1,85 @@ +! DART software - Copyright UCAR. This open source software is provided +! by UCAR, "as is", without charge, subject to all terms of use at +! http://www.image.ucar.edu/DAReS/DART/DART_download + +program test_parse_variables + +use utilities_mod, only : find_namelist_in_file, check_namelist_read, & + initialize_utilities, finalize_utilities +use types_mod, only : vtablenamelength, MISSING_R8, r8 +use default_model_mod, only : parse_variables_clamp, parse_variables, & + state_var_type, & + MAX_STATE_VARIABLE_FIELDS_CLAMP, & + MAX_STATE_VARIABLE_FIELDS +use obs_kind_mod, only : QTY_SALINITY, QTY_POTENTIAL_TEMPERATURE, & + QTY_U_CURRENT_COMPONENT + +use test ! fortran-testanything + +implicit none + +integer :: iunit, io +type(state_var_type) :: state_vars, state_vars_clamp + +character(len=vtablenamelength) :: state_variables(MAX_STATE_VARIABLE_FIELDS) = ' ' +character(len=vtablenamelength) :: state_variables_clamp(MAX_STATE_VARIABLE_FIELDS_CLAMP) = ' ' + +namelist /model_nml/ & + state_variables + +namelist /model_nml_clamp/ & + state_variables_clamp + +call initialize_utilities('test_parse_variables') + +call plan(28) + +! Using namelist entry WITHOUT clamping values + +call find_namelist_in_file('input.nml', 'model_nml', iunit) +read(iunit, nml = model_nml, iostat = io) +call check_namelist_read(iunit, io, 'model_nml') + +state_vars = parse_variables(state_variables) + +call ok(state_vars%nvars == 3) +call ok(state_vars%netcdf_var_names(1) == 'SALT_CUR') +call ok(state_vars%netcdf_var_names(2) == 'TEMP_CUR') +call ok(state_vars%netcdf_var_names(3) == 'UVEL_CUR') +call ok(state_vars%qtys(1) == QTY_SALINITY) +call ok(state_vars%qtys(2) == QTY_POTENTIAL_TEMPERATURE) +call ok(state_vars%qtys(3) == QTY_U_CURRENT_COMPONENT) +call ok(allocated(state_vars_clamp%clamp_values) .eqv. .false.) +call ok(state_vars%updates(1) .eqv. .true.) +call ok(state_vars%updates(2) .eqv. .false.) +call ok(state_vars%updates(3) .eqv. .true.) + +! Using namelist entry WITH clamping values + +call find_namelist_in_file('input.nml', 'model_nml_clamp', iunit) +read(iunit, nml = model_nml_clamp, iostat = io) +call check_namelist_read(iunit, io, 'model_nml_clamp') + +state_vars_clamp = parse_variables_clamp(state_variables_clamp) + +call ok(state_vars_clamp%nvars == 3) +call ok(state_vars_clamp%netcdf_var_names(1) == 'SALT_CUR') +call ok(state_vars_clamp%netcdf_var_names(2) == 'TEMP_CUR') +call ok(state_vars_clamp%netcdf_var_names(3) == 'UVEL_CUR') +call ok(state_vars_clamp%qtys(1) == QTY_SALINITY) +call ok(state_vars_clamp%qtys(2) == QTY_POTENTIAL_TEMPERATURE) +call ok(state_vars_clamp%qtys(3) == QTY_U_CURRENT_COMPONENT) +call ok(allocated(state_vars_clamp%clamp_values) .eqv. .true.) +call ok(state_vars_clamp%clamp_values(1,1) == 0.0_r8) +call ok(state_vars_clamp%clamp_values(1,2) == 0.0_r8) +call ok(state_vars_clamp%clamp_values(2,1) == 0.0_r8) +call ok(state_vars_clamp%clamp_values(2,2) == 0.0_r8) +call ok(state_vars_clamp%clamp_values(3,1) == 0.0_r8) +call ok(state_vars_clamp%clamp_values(3,2) == 0.0_r8) +call ok(state_vars_clamp%updates(1) .eqv. .true.) +call ok(state_vars_clamp%updates(2) .eqv. .false.) +call ok(state_vars_clamp%updates(3) .eqv. .true.) + +call finalize_utilities() + +end program test_parse_variables diff --git a/developer_tests/namelist/work/input.nml b/developer_tests/namelist/work/input.nml new file mode 100644 index 0000000000..5cb8222ada --- /dev/null +++ b/developer_tests/namelist/work/input.nml @@ -0,0 +1,67 @@ +&model_nml_clamp + state_variables_clamp = 'SALT_CUR ', 'QTY_SALINITY', '0.0', '0.0', 'UPDATE', + 'TEMP_CUR ', 'QTY_POTENTIAL_TEMPERATURE', '0.0', '0.0', 'NO_COPY_BACK', + 'UVEL_CUR ', 'QTY_U_CURRENT_COMPONENT ', '0', '0.0', 'UPDATE', + / + +&model_nml + state_variables = 'SALT_CUR ', 'QTY_SALINITY', 'UPDATE', + 'TEMP_CUR ', 'QTY_POTENTIAL_TEMPERATURE', 'NO_COPY_BACK', + 'UVEL_CUR ', 'QTY_U_CURRENT_COMPONENT ', 'update' + / + +&utilities_nml + module_details = .false. + write_nml = 'none' + / + +&preprocess_nml + input_obs_kind_mod_file = '../../../assimilation_code/modules/observations/DEFAULT_obs_kind_mod.F90' + quantity_files = '../../../assimilation_code/modules/observations/land_quantities_mod.f90', + '../../../assimilation_code/modules/observations/default_quantities_mod.f90' + '../../../assimilation_code/modules/observations/atmosphere_quantities_mod.f90' + '../../../assimilation_code/modules/observations/ocean_quantities_mod.f90' + output_obs_kind_mod_file = '../../../assimilation_code/modules/observations/obs_kind_mod.f90' + input_obs_def_mod_file = '../../../observations/forward_operators/DEFAULT_obs_def_mod.F90' + output_obs_def_mod_file = '../../../observations/forward_operators/obs_def_mod.f90' + input_files = '../../../observations/forward_operators/obs_def_AIRS_mod.f90', + '../../../observations/forward_operators/obs_def_AOD_mod.f90', + '../../../observations/forward_operators/obs_def_AURA_mod.f90', + '../../../observations/forward_operators/obs_def_COSMOS_mod.f90', + '../../../observations/forward_operators/obs_def_CO_Nadir_mod.f90', + '../../../observations/forward_operators/obs_def_GWD_mod.f90', + '../../../observations/forward_operators/obs_def_QuikSCAT_mod.f90', + '../../../observations/forward_operators/obs_def_SABER_mod.f90', + '../../../observations/forward_operators/obs_def_altimeter_mod.f90', + '../../../observations/forward_operators/obs_def_cloud_mod.f90', + '../../../observations/forward_operators/obs_def_dew_point_mod.f90', + '../../../observations/forward_operators/obs_def_dwl_mod.f90', + '../../../observations/forward_operators/obs_def_eval_mod.f90', + '../../../observations/forward_operators/obs_def_gps_mod.f90', + '../../../observations/forward_operators/obs_def_gts_mod.f90', + '../../../observations/forward_operators/obs_def_land_mod.f90', + '../../../observations/forward_operators/obs_def_metar_mod.f90', + '../../../observations/forward_operators/obs_def_ocean_mod.f90', + '../../../observations/forward_operators/obs_def_pe2lyr_mod.f90', + '../../../observations/forward_operators/obs_def_radar_mod.f90', + '../../../observations/forward_operators/obs_def_reanalysis_bufr_mod.f90', + '../../../observations/forward_operators/obs_def_rel_humidity_mod.f90', + '../../../observations/forward_operators/obs_def_sqg_mod.f90', + '../../../observations/forward_operators/obs_def_tower_mod.f90', + '../../../observations/forward_operators/obs_def_tpw_mod.f90', + '../../../observations/forward_operators/obs_def_upper_atm_mod.f90', + '../../../observations/forward_operators/obs_def_vortex_mod.f90', + '../../../observations/forward_operators/obs_def_wind_speed_mod.f90' + / + +&mpi_utilities_nml + / + +&obs_kind_nml + / + +&ensemble_manager_nml + / + +&state_vector_io_nml + / diff --git a/developer_tests/namelist/work/quickbuild.sh b/developer_tests/namelist/work/quickbuild.sh new file mode 100755 index 0000000000..61d4f99ebf --- /dev/null +++ b/developer_tests/namelist/work/quickbuild.sh @@ -0,0 +1,43 @@ +#!/usr/bin/env bash + +# DART software - Copyright UCAR. This open source software is provided +# by UCAR, "as is", without charge, subject to all terms of use at +# http://www.image.ucar.edu/DAReS/DART/DART_download + +main() { + +export DART=$(git rev-parse --show-toplevel) +source "$DART"/build_templates/buildfunctions.sh + +MODEL="none" +EXTRA="$DART"/models/template/threed_model_mod.f90 +dev_test=1 +TEST="namelist" +LOCATION="threed_sphere" + +serial_programs=( +test_parse_variables +) + + +# quickbuild arguments +arguments "$@" + +# clean the directory +\rm -f -- *.o *.mod Makefile .cppdefs + +# build any NetCDF files from .cdl files +cdl_to_netcdf + +# build and run preprocess before making any other DART executables +buildpreprocess + +# build +buildit + +# clean up +\rm -f -- *.o *.mod + +} + +main "$@" diff --git a/models/MOM6/model_mod.f90 b/models/MOM6/model_mod.f90 index e0b88ae288..7fd9c0b642 100644 --- a/models/MOM6/model_mod.f90 +++ b/models/MOM6/model_mod.f90 @@ -21,11 +21,9 @@ module model_mod get_location, query_location, VERTISLEVEL, & VERTISHEIGHT, set_vertical -use utilities_mod, only : error_handler, & - E_ERR, E_MSG, & +use utilities_mod, only : error_handler, E_ERR, E_MSG, & nmlfileunit, do_output, do_nml_file, do_nml_term, & - find_namelist_in_file, check_namelist_read, & - to_upper + find_namelist_in_file, check_namelist_read use netcdf_utilities_mod, only : nc_add_global_attribute, nc_synchronize_file, & nc_add_global_creation_time, & @@ -61,7 +59,9 @@ module model_mod use default_model_mod, only : pert_model_copies, write_model_time, & init_time => fail_init_time, & init_conditions => fail_init_conditions, & - convert_vertical_obs, adv_1step + convert_vertical_obs, adv_1step, & + parse_variables, & + MAX_STATE_VARIABLE_FIELDS implicit none private @@ -106,10 +106,6 @@ module model_mod ! Ocean vs land real(r8), allocatable :: wet(:,:), basin_depth(:,:) -! DART state vector contents are specified in the input.nml:&model_nml namelist. -integer, parameter :: MAX_STATE_VARIABLES = 10 -integer, parameter :: NUM_STATE_TABLE_COLUMNS = 3 - ! model_interpolate failure codes integer, parameter :: NOT_IN_STATE = 12 integer, parameter :: THICKNESS_NOT_IN_STATE = 13 @@ -128,7 +124,7 @@ module model_mod character(len=256) :: ocean_geometry = 'ocean_geometry.nc' integer :: assimilation_period_days = -1 integer :: assimilation_period_seconds = -1 -character(len=vtablenamelength) :: model_state_variables(MAX_STATE_VARIABLES * NUM_STATE_TABLE_COLUMNS ) = ' ' +character(len=vtablenamelength) :: model_state_variables(MAX_STATE_VARIABLE_FIELDS) = ' ' namelist /model_nml/ template_file, static_file, ocean_geometry, assimilation_period_days, & assimilation_period_seconds, model_state_variables @@ -151,15 +147,6 @@ module model_mod subroutine static_init_model() integer :: iunit, io -character(len=vtablenamelength) :: variable_table(MAX_STATE_VARIABLES, NUM_STATE_TABLE_COLUMNS) - -integer :: state_qty_list(MAX_STATE_VARIABLES) -logical :: update_var_list(MAX_STATE_VARIABLES) - -! identifiers for variable_table -integer, parameter :: VAR_NAME_INDEX = 1 -integer, parameter :: VAR_QTY_INDEX = 2 -integer, parameter :: VAR_UPDATE_INDEX = 3 module_initialized = .true. @@ -167,7 +154,7 @@ subroutine static_init_model() read(iunit, nml = model_nml, iostat = io) call check_namelist_read(iunit, io, "model_nml") -! Record the namelist values used for the run +! Record the namelist values used for the run if (do_nml_file()) write(nmlfileunit, nml=model_nml) if (do_nml_term()) write( * , nml=model_nml) @@ -181,15 +168,11 @@ subroutine static_init_model() assimilation_time_step = set_time(assimilation_period_seconds, & assimilation_period_days) -! verify that the model_state_variables namelist was filled in correctly. -! returns variable_table which has variable names, kinds and update strings. -call verify_state_variables(model_state_variables, nfields, variable_table, state_qty_list, update_var_list) - -! Define which variables are in the model state -dom_id = add_domain(template_file, nfields, & - var_names = variable_table(1:nfields, VAR_NAME_INDEX), & - kind_list = state_qty_list(1:nfields), & - update_list = update_var_list(1:nfields)) +! Define which variables are in the model state; +! parse_variables converts the character table that was read in from +! model_nml:model_state_variables to a state_var_type that can be passed +! to add_domain +dom_id = add_domain(template_file, parse_variables(model_state_variables)) model_size = get_domain_size(dom_id) @@ -1018,77 +1001,6 @@ elemental function sensible_temp(potemp, s, lpres) end function sensible_temp -!------------------------------------------------------------------ -! Verify that the namelist was filled in correctly, and check -! that there are valid entries for the dart_kind. -! Returns a table with columns: -! -! netcdf_variable_name ; dart_qty_string ; update_string - -subroutine verify_state_variables(state_variables, ngood, table, qty_list, update_var) - -character(len=*), intent(inout) :: state_variables(:) -integer, intent(out) :: ngood -character(len=*), intent(out) :: table(:,:) -integer, intent(out) :: qty_list(:) ! kind number -logical, intent(out) :: update_var(:) ! logical update - -integer :: nrows, i -character(len=NF90_MAX_NAME) :: varname, dartstr, update -character(len=256) :: string1, string2 - -if ( .not. module_initialized ) call static_init_model - -nrows = size(table,1) - -ngood = 0 - -MyLoop : do i = 1, nrows - - varname = trim(state_variables(3*i -2)) - dartstr = trim(state_variables(3*i -1)) - update = trim(state_variables(3*i )) - - call to_upper(update) - - table(i,1) = trim(varname) - table(i,2) = trim(dartstr) - table(i,3) = trim(update) - - if ( table(i,1) == ' ' .and. table(i,2) == ' ' .and. table(i,3) == ' ') exit MyLoop ! Found end of list. - - if ( table(i,1) == ' ' .or. table(i,2) == ' ' .or. table(i,3) == ' ' ) then - string1 = 'model_nml:model_state_variables not fully specified' - call error_handler(E_ERR,'verify_state_variables',string1) - endif - - ! Make sure DART qty is valid - - qty_list(i) = get_index_for_quantity(dartstr) - if( qty_list(i) < 0 ) then - write(string1,'(''there is no obs_kind <'',a,''> in obs_kind_mod.f90'')') trim(dartstr) - call error_handler(E_ERR,'verify_state_variables',string1) - endif - - ! Make sure the update variable has a valid name - - select case (update) - case ('UPDATE') - update_var(i) = .true. - case ('NO_COPY_BACK') - update_var(i) = .false. - case default - write(string1,'(A)') 'only UPDATE or NO_COPY_BACK supported in model_state_variable namelist' - write(string2,'(6A)') 'you provided : ', trim(varname), ', ', trim(dartstr), ', ', trim(update) - call error_handler(E_ERR,'verify_state_variables',string1, text2=string2) - end select - - ngood = ngood + 1 -enddo MyLoop - - -end subroutine verify_state_variables - !-------------------------------------------------------------------- function read_model_time(filename) diff --git a/models/POP/model_mod.f90 b/models/POP/model_mod.f90 index 23ddfb9a7c..24a66017f3 100644 --- a/models/POP/model_mod.f90 +++ b/models/POP/model_mod.f90 @@ -2,7 +2,6 @@ ! by UCAR, "as is", without charge, subject to all terms of use at ! http://www.image.ucar.edu/DAReS/DART/DART_download ! -! $Id$ module model_mod @@ -49,8 +48,11 @@ module model_mod use state_structure_mod, only : add_domain, get_model_variable_indices, & get_num_variables, get_index_start, & get_num_dims, get_domain_size, & - get_dart_vector_index -use default_model_mod, only : adv_1step, init_time, init_conditions, nc_write_model_vars + get_dart_vector_index, get_varid_from_kind, & + get_kind_index +use default_model_mod, only : adv_1step, init_time, init_conditions, & + nc_write_model_vars, parse_variables, & + MAX_STATE_VARIABLE_FIELDS use typesizes use netcdf @@ -106,24 +108,12 @@ module model_mod ! Storage for a random sequence for perturbing a single initial state type(random_seq_type) :: random_seq -! DART state vector contents are specified in the input.nml:&model_nml namelist. -integer, parameter :: max_state_variables = 10 -integer, parameter :: num_state_table_columns = 3 -character(len=vtablenamelength) :: variable_table( max_state_variables, num_state_table_columns ) -integer :: state_kinds_list( max_state_variables ) -logical :: update_var_list( max_state_variables ) - -! identifiers for variable_table -integer, parameter :: VAR_NAME_INDEX = 1 -integer, parameter :: VAR_QTY_INDEX = 2 -integer, parameter :: VAR_UPDATE_INDEX = 3 - ! things which can/should be in the model_nml integer :: assimilation_period_days = -1 integer :: assimilation_period_seconds = -1 real(r8) :: model_perturbation_amplitude = 0.2 logical :: update_dry_cell_walls = .false. -character(len=vtablenamelength) :: model_state_variables(max_state_variables * num_state_table_columns ) = ' ' +character(len=vtablenamelength) :: model_state_variables(MAX_STATE_VARIABLE_FIELDS) = ' ' integer :: debug = 0 ! turn up for more and more debug messages ! only valid values: native, big_endian, little_endian @@ -157,9 +147,6 @@ module model_mod ! is not being used. !------------------------------------------------------------------ -! Number of fields in the state vector -integer :: nfields - ! Grid parameters - the values will be read from a ! standard POP namelist and filled in here. @@ -190,8 +177,6 @@ module model_mod integer(i8) :: model_size ! the state vector length - - !------------------------------------------------ ! NOTE (dipole/tripole grids): since both of the dipole and tripole @@ -335,10 +320,6 @@ subroutine static_init_model() if (debug > 2) call write_grid_netcdf() ! DEBUG only if (debug > 2) call write_grid_interptest() ! DEBUG only -! verify that the model_state_variables namelist was filled in correctly. -! returns variable_table which has variable names, kinds and update strings. -call verify_state_variables(model_state_variables, nfields, variable_table, state_kinds_list, update_var_list) - ! in spite of the staggering, all grids are the same size ! and offset by half a grid cell. 4 are 3D and 1 is 2D. ! e.g. S,T,U,V = 256 x 225 x 70 @@ -355,10 +336,17 @@ subroutine static_init_model() ! Initialize the interpolation routines call init_interp() +if ( model_state_variables(1) == ' ' ) then ! no model_state_variables namelist provided + call use_default_state_variables( model_state_variables ) + string1 = 'model_nml:model_state_variables not specified - using default variables' + call error_handler(E_MSG,'static_init_model',string1,source,revision,revdate) +endif + !> @todo 'pop.r.nc' is hardcoded in dart_pop_mod.f90 -domain_id = add_domain('pop.r.nc', nfields, & - var_names = variable_table(1:nfields, VAR_NAME_INDEX), & - update_list = update_var_list(1:nfields)) +! parse_variables converts the character table that was read in from +! model_nml:model_state_variables and returns a state_var_type that +! can be passed to add_domain +domain_id = add_domain('pop.r.nc', parse_variables(model_state_variables)) model_size = get_domain_size(domain_id) if (do_output()) write(*,*) 'model_size = ', model_size @@ -841,7 +829,7 @@ subroutine model_interpolate(state_handle, ens_size, location, obs_type, expecte source, revision, revdate, text2=string2, text3=string3) endif - base_offset = get_index_start(domain_id, get_varid_from_kind(QTY_SEA_SURFACE_PRESSURE)) + base_offset = get_index_start(domain_id, get_varid_from_kind(domain_id, QTY_SEA_SURFACE_PRESSURE)) call lon_lat_interpolate(state_handle, ens_size, base_offset, llon, llat, & QTY_SEA_SURFACE_HEIGHT, 1, expected_obs, istatus, expected_mdt) @@ -862,10 +850,10 @@ subroutine model_interpolate(state_handle, ens_size, location, obs_type, expecte QTY_U_CURRENT_COMPONENT, & QTY_V_CURRENT_COMPONENT, & QTY_SEA_SURFACE_PRESSURE) - base_offset = get_index_start(domain_id, get_varid_from_kind(obs_type)) + base_offset = get_index_start(domain_id, get_varid_from_kind(domain_id, obs_type)) CASE (QTY_SEA_SURFACE_HEIGHT) - base_offset = get_index_start(domain_id, get_varid_from_kind(QTY_SEA_SURFACE_PRESSURE)) + base_offset = get_index_start(domain_id, get_varid_from_kind(domain_id, QTY_SEA_SURFACE_PRESSURE)) convert_to_ssh = .TRUE. ! simple linear transform of PSURF CASE DEFAULT @@ -1754,7 +1742,7 @@ subroutine get_state_meta_data(index_in, location, var_type) if ( .not. module_initialized ) call static_init_model call get_model_variable_indices(index_in, lon_index, lat_index, depth_index, var_id=var_id) -call get_state_kind(var_id, local_var) +local_var = get_kind_index(domain_id, var_id) if (is_on_ugrid(local_var)) then lon = ULON(lon_index, lat_index) @@ -1784,52 +1772,6 @@ subroutine get_state_meta_data(index_in, location, var_type) end subroutine get_state_meta_data -!-------------------------------------------------------------------- -!> given a DART kind, return the variable number (position in the list) - - -function get_varid_from_kind(dart_kind) - -integer, intent(in) :: dart_kind -integer :: get_varid_from_kind - -integer :: i - -do i = 1, get_num_variables(domain_id) - if (dart_kind == state_kinds_list(i)) then - get_varid_from_kind = i - return - endif -end do - -write(string1, *) 'Kind ', dart_kind, ' not found in state vector' -write(string2, *) 'AKA ', get_name_for_quantity(dart_kind), ' not found in state vector' -call error_handler(E_MSG,'get_varid_from_kind', string1, & - source, revision, revdate, text2=string2) - -get_varid_from_kind = -1 - -end function get_varid_from_kind - - -!------------------------------------------------------------------ -!> Given an integer index into the state vector structure, returns the kind, -!> and both the starting offset for this kind, as well as the offset into -!> the block of this kind. - - -subroutine get_state_kind(var_ind, var_type) - -integer, intent(in) :: var_ind -integer, intent(out) :: var_type - -if ( .not. module_initialized ) call static_init_model - -var_type = state_kinds_list(var_ind) - -end subroutine get_state_kind - - !------------------------------------------------------------------ !> Given an integer index into the state vector structure, returns the !> type, taking into account the ocean bottom and dry land. @@ -1845,7 +1787,7 @@ subroutine get_state_kind_inc_dry(index_in, var_type) if ( .not. module_initialized ) call static_init_model call get_model_variable_indices(index_in, lon_index, lat_index, depth_index, var_id=var_id) -call get_state_kind(var_id, var_type) +var_type = get_kind_index(domain_id, var_id) ! if on land or below ocean floor, replace type with dry land. if(is_dry_land(var_type, lon_index, lat_index, depth_index)) then @@ -2361,8 +2303,8 @@ subroutine compute_temperature(state_handle, ens_size, llon, llat, lheight, expe return endif -offset_salt = get_index_start(domain_id, get_varid_from_kind(QTY_SALINITY)) -offset_temp = get_index_start(domain_id, get_varid_from_kind(QTY_POTENTIAL_TEMPERATURE)) +offset_salt = get_index_start(domain_id, get_varid_from_kind(domain_id, QTY_SALINITY)) +offset_temp = get_index_start(domain_id, get_varid_from_kind(domain_id, QTY_POTENTIAL_TEMPERATURE)) ! salinity - in msu (kg/kg). converter will want psu (g/kg). call do_interp(state_handle, ens_size, offset_salt, hgt_bot, hgt_top, hgt_fract, llon, llat, & @@ -2746,103 +2688,6 @@ subroutine vert_convert(state_handle, location, obs_kind, istatus) end subroutine vert_convert -!------------------------------------------------------------------ -!> Verify that the namelist was filled in correctly, and check -!> that there are valid entries for the dart_kind. -!> Returns a table with columns: -!> -!> netcdf_variable_name ; dart_kind_string ; update_string - - -subroutine verify_state_variables( state_variables, ngood, table, kind_list, update_var ) - -character(len=*), intent(inout) :: state_variables(:) -integer, intent(out) :: ngood -character(len=*), intent(out) :: table(:,:) -integer, intent(out) :: kind_list(:) ! kind number -logical, optional, intent(out) :: update_var(:) ! logical update - -integer :: nrows, i -character(len=NF90_MAX_NAME) :: varname, dartstr, update - -if ( .not. module_initialized ) call static_init_model - -nrows = size(table,1) - -ngood = 0 - -if ( state_variables(1) == ' ' ) then ! no model_state_variables namelist provided - call use_default_state_variables( state_variables ) - string1 = 'model_nml:model_state_variables not specified using default variables' - call error_handler(E_MSG,'verify_state_variables',string1,source,revision,revdate) -endif - -MyLoop : do i = 1, nrows - - varname = trim(state_variables(3*i -2)) - dartstr = trim(state_variables(3*i -1)) - update = trim(state_variables(3*i )) - - call to_upper(update) - - table(i,1) = trim(varname) - table(i,2) = trim(dartstr) - table(i,3) = trim(update) - - if ( table(i,1) == ' ' .and. table(i,2) == ' ' .and. table(i,3) == ' ') exit MyLoop ! Found end of list. - - if ( table(i,1) == ' ' .or. table(i,2) == ' ' .or. table(i,3) == ' ' ) then - string1 = 'model_nml:model_state_variables not fully specified' - call error_handler(E_ERR,'verify_state_variables',string1,source,revision,revdate) - endif - - ! Make sure DART kind is valid - - kind_list(i) = get_index_for_quantity(dartstr) - if( kind_list(i) < 0 ) then - write(string1,'(''there is no obs_kind <'',a,''> in obs_kind_mod.f90'')') trim(dartstr) - call error_handler(E_ERR,'verify_state_variables',string1,source,revision,revdate) - endif - - ! Make sure the update variable has a valid name - - if ( present(update_var) )then - SELECT CASE (update) - CASE ('UPDATE') - update_var(i) = .true. - CASE ('NO_COPY_BACK') - update_var(i) = .false. - CASE DEFAULT - write(string1,'(A)') 'only UPDATE or NO_COPY_BACK supported in model_state_variable namelist' - write(string2,'(6A)') 'you provided : ', trim(varname), ', ', trim(dartstr), ', ', trim(update) - call error_handler(E_ERR,'verify_state_variables',string1,source,revision,revdate, text2=string2) - END SELECT - endif - - ! Record the contents of the DART state vector - - if (do_output()) then - write(string1,'(A,I2,6A)') 'variable ',i,' is ',trim(varname), ', ', trim(dartstr), ', ', trim(update) - call error_handler(E_MSG,'verify_state_variables',string1,source,revision,revdate) - endif - - ngood = ngood + 1 -enddo MyLoop - -! check to see if temp and salinity are both in the state otherwise you will not -! be able to interpolate in XXX subroutine -if ( any(kind_list == QTY_SALINITY) ) then - ! check to see that temperature is also in the variable list - if ( .not. any(kind_list == QTY_POTENTIAL_TEMPERATURE) ) then - write(string1,'(A)') 'in order to compute temperature you need to have both ' - write(string2,'(A)') 'QTY_SALINITY and QTY_POTENTIAL_TEMPERATURE in the model state' - call error_handler(E_ERR,'verify_state_variables',string1,source,revision,revdate, text2=string2) - endif -endif - -end subroutine verify_state_variables - - !------------------------------------------------------------------ !> Default state_variables from the original pop model_mod. Must !> keep in the same order to be consistent with previous versions. @@ -2853,7 +2698,7 @@ subroutine use_default_state_variables( state_variables ) character(len=*), intent(inout) :: state_variables(:) ! strings must all be the same length for the gnu compiler -state_variables( 1:5*num_state_table_columns ) = & +state_variables( 1:15 ) = & (/ 'SALT_CUR ', 'QTY_SALINITY ', 'UPDATE ', & 'TEMP_CUR ', 'QTY_POTENTIAL_TEMPERATURE ', 'UPDATE ', & 'UVEL_CUR ', 'QTY_U_CURRENT_COMPONENT ', 'UPDATE ', & @@ -2868,8 +2713,3 @@ end subroutine use_default_state_variables end module model_mod -! -! $URL$ -! $Id$ -! $Revision$ -! $Date$ diff --git a/models/aether_lat-lon/model_mod.f90 b/models/aether_lat-lon/model_mod.f90 index 87ac6c0e9c..39b76ed67e 100644 --- a/models/aether_lat-lon/model_mod.f90 +++ b/models/aether_lat-lon/model_mod.f90 @@ -61,7 +61,8 @@ module model_mod pert_model_copies, read_model_time, write_model_time, & init_time => fail_init_time, & init_conditions => fail_init_conditions, & - convert_vertical_obs, convert_vertical_state, adv_1step + convert_vertical_obs, convert_vertical_state, adv_1step, & + parse_variables_clamp, MAX_STATE_VARIABLE_FIELDS_CLAMP implicit none private @@ -97,17 +98,7 @@ module model_mod integer :: time_step_days = 0 integer :: time_step_seconds = 3600 -integer, parameter :: MAX_STATE_VARIABLES = 100 -integer, parameter :: NUM_STATE_TABLE_COLUMNS = 5 -character(len=vtablenamelength) :: variables(NUM_STATE_TABLE_COLUMNS,MAX_STATE_VARIABLES) = '' - -type :: var_type - integer :: count - character(len=64), allocatable :: names(:) - integer, allocatable :: qtys(:) - real(r8), allocatable :: clamp_values(:, :) - logical, allocatable :: updates(:) -end type var_type +character(len=vtablenamelength) :: variables(MAX_STATE_VARIABLE_FIELDS_CLAMP) = '' namelist /model_nml/ template_file, time_step_days, time_step_seconds, variables @@ -157,7 +148,6 @@ module model_mod subroutine static_init_model() integer :: iunit, io -type(var_type) :: var module_initialized = .true. @@ -179,8 +169,6 @@ subroutine static_init_model() lat_start = lats(1) lat_delta = lats(2) - lats(1) -var = assign_var(variables, MAX_STATE_VARIABLES) - ! This time is both the minimum time you can ask the model to advance ! (for models that can be advanced by filter) and it sets the assimilation ! window. All observations within +/- 1/2 this interval from the current @@ -188,10 +176,11 @@ subroutine static_init_model() ! feel free to hardcode it and remove from the namelist. assimilation_time_step = set_time(time_step_seconds, time_step_days) -! Define which variables are in the model state -! This is using add_domain_from_file (arg list matches) -dom_id = add_domain(template_file, var%count, var%names, var%qtys, & - var%clamp_values, var%updates) +! Define which variables are in the model state; +! parse_variables converts the character table that was read in from +! model_nml:model_state_variables and returns a state_var_type that +! can be passed to add_domain +dom_id = add_domain(template_file, parse_variables_clamp(variables)) call state_structure_info(dom_id) @@ -409,72 +398,6 @@ subroutine assign_dimensions() end subroutine assign_dimensions -!----------------------------------------------------------------------- -! Parse the table of variables characteristics into arrays for easier access. - -function assign_var(variables, MAX_STATE_VARIABLES) result(var) - -character(len=vtablenamelength), intent(in) :: variables(:, :) -integer, intent(in) :: MAX_STATE_VARIABLES - -type(var_type) :: var -integer :: ivar -character(len=vtablenamelength) :: table_entry - -!----------------------------------------------------------------------- -! Codes for interpreting the NUM_STATE_TABLE_COLUMNS of the variables table -integer, parameter :: NAME_INDEX = 1 ! ... variable name -integer, parameter :: QTY_INDEX = 2 ! ... DART qty -integer, parameter :: MIN_VAL_INDEX = 3 ! ... minimum value if any -integer, parameter :: MAX_VAL_INDEX = 4 ! ... maximum value if any -integer, parameter :: UPDATE_INDEX = 5 ! ... update (state) or not - -! Loop through the variables array to get the actual count of the number of variables -do ivar = 1, MAX_STATE_VARIABLES - ! If the element is an empty string, the loop has exceeded the extent of the variables - if (variables(1, ivar) == '') then - var%count = ivar-1 - exit - endif -enddo - -! Allocate the arrays in the var derived type -allocate(var%names(var%count), var%qtys(var%count), var%clamp_values(var%count, 2), var%updates(var%count)) - -do ivar = 1, var%count - - var%names(ivar) = trim(variables(NAME_INDEX, ivar)) - - table_entry = variables(QTY_INDEX, ivar) - call to_upper(table_entry) - - var%qtys(ivar) = get_index_for_quantity(table_entry) - - if (variables(MIN_VAL_INDEX, ivar) /= 'NA') then - read(variables(MIN_VAL_INDEX, ivar), '(d16.8)') var%clamp_values(ivar,1) - else - var%clamp_values(ivar,1) = MISSING_R8 - endif - - if (variables(MAX_VAL_INDEX, ivar) /= 'NA') then - read(variables(MAX_VAL_INDEX, ivar), '(d16.8)') var%clamp_values(ivar,2) - else - var%clamp_values(ivar,2) = MISSING_R8 - endif - - table_entry = variables(UPDATE_INDEX, ivar) - call to_upper(table_entry) - - if (table_entry == 'UPDATE') then - var%updates(ivar) = .true. - else - var%updates(ivar) = .false. - endif - -enddo - -end function assign_var - !----------------------------------------------------------------------- ! Extract state values needed by the interpolation from all ensemble members. diff --git a/models/cam-common-code/cam_common_code_mod.f90 b/models/cam-common-code/cam_common_code_mod.f90 index 4970e16474..dc635c268b 100644 --- a/models/cam-common-code/cam_common_code_mod.f90 +++ b/models/cam-common-code/cam_common_code_mod.f90 @@ -62,9 +62,9 @@ module cam_common_code_mod read_model_time, ref_model_top_pressure, ref_nlevels, scale_height, & set_vert_localization, vert_interp, vertical_localization_type, write_model_time -public :: nc_write_model_atts, grid_data, read_grid_info, set_cam_variable_info, & - MAX_STATE_VARIABLES, num_state_table_columns, common_initialized, & - MAX_PERT, shortest_time_between_assimilations, domain_id, & +public :: nc_write_model_atts, grid_data, read_grid_info, & + common_initialized, MAX_PERT, & + shortest_time_between_assimilations, domain_id, & ccustom_routine_to_generate_ensemble, & cfields_to_perturb, & cperturbation_amplitude, & @@ -107,8 +107,6 @@ module cam_common_code_mod ! info and is required for getting state variables. integer :: domain_id = -1 -integer, parameter :: MAX_STATE_VARIABLES = 100 -integer, parameter :: num_state_table_columns = 5 ! maximum number of fields you can list to be perturbed ! to generate an ensemble if starting from a single state. integer, parameter :: MAX_PERT = 100 @@ -172,94 +170,6 @@ module cam_common_code_mod contains - -!----------------------------------------------------------------------- -!> -!> Fill the array of requested variables, dart kinds, possible min/max -!> values and whether or not to update the field in the output file. -!> Then calls 'add_domain()' to tell the DART code which variables to -!> read into the state vector after this code returns. -!> -!>@param variable_array the list of variables and kinds from model_mod_nml -!>@param nfields the number of variable/Quantity pairs specified - -subroutine set_cam_variable_info(cam_template_filename, variable_array) - -character(len=*), intent(in) :: cam_template_filename -character(len=*), intent(in) :: variable_array(:) - -character(len=*), parameter :: routine = 'set_cam_variable_info:' - -integer :: i, nfields -integer, parameter :: MAX_STRING_LEN = 128 - -character(len=MAX_STRING_LEN) :: varname ! column 1, NetCDF variable name -character(len=MAX_STRING_LEN) :: dartstr ! column 2, DART Quantity -character(len=MAX_STRING_LEN) :: minvalstr ! column 3, Clamp min val -character(len=MAX_STRING_LEN) :: maxvalstr ! column 4, Clamp max val -character(len=MAX_STRING_LEN) :: updatestr ! column 5, Update output or not - -character(len=vtablenamelength) :: var_names(MAX_STATE_VARIABLES) = ' ' -logical :: update_list(MAX_STATE_VARIABLES) = .FALSE. -integer :: kind_list(MAX_STATE_VARIABLES) = MISSING_I -real(r8) :: clamp_vals(MAX_STATE_VARIABLES,2) = MISSING_R8 - -nfields = 0 -ParseVariables : do i = 1, MAX_STATE_VARIABLES - - varname = variable_array(num_state_table_columns*i-4) - dartstr = variable_array(num_state_table_columns*i-3) - minvalstr = variable_array(num_state_table_columns*i-2) - maxvalstr = variable_array(num_state_table_columns*i-1) - updatestr = variable_array(num_state_table_columns*i ) - - if ( varname == ' ' .and. dartstr == ' ' ) exit ParseVariables ! Found end of list. - - if ( varname == ' ' .or. dartstr == ' ' ) then - string1 = 'model_nml:model "state_variables" not fully specified' - call error_handler(E_ERR,routine,string1,source,revision,revdate) - endif - - ! Make sure DART kind is valid - - if( get_index_for_quantity(dartstr) < 0 ) then - write(string1,'(3A)') 'there is no obs_kind "', trim(dartstr), '" in obs_kind_mod.f90' - call error_handler(E_ERR,routine,string1,source,revision,revdate) - endif - - call to_upper(minvalstr) - call to_upper(maxvalstr) - call to_upper(updatestr) - - var_names( i) = varname - kind_list( i) = get_index_for_quantity(dartstr) - clamp_vals(i,1) = string_to_real(minvalstr) - clamp_vals(i,2) = string_to_real(maxvalstr) - update_list( i) = string_to_logical(updatestr, 'UPDATE') - - nfields = nfields + 1 - -enddo ParseVariables - -if (nfields == MAX_STATE_VARIABLES) then - write(string1,'(2A)') 'WARNING: There is a possibility you need to increase ', & - 'MAX_STATE_VARIABLES in the global variables in model_mod.f90' - - write(string2,'(A,i4,A)') 'WARNING: you have specified at least ', nfields, & - ' perhaps more' - - call error_handler(E_MSG,routine,string1,source,revision,revdate,text2=string2) -endif - -! CAM only has a single domain (only a single grid, no nests or multiple grids) - -domain_id = add_domain(cam_template_filename, nfields, var_names, kind_list, & - clamp_vals, update_list) - -end subroutine set_cam_variable_info - - - !----------------------------------------------------------------------- !> Read the data from the various cam grid arrays !> diff --git a/models/cam-fv/model_mod.f90 b/models/cam-fv/model_mod.f90 index edc90380e6..a2fd4adf8b 100644 --- a/models/cam-fv/model_mod.f90 +++ b/models/cam-fv/model_mod.f90 @@ -81,7 +81,9 @@ module model_mod QUAD_LOCATED_CELL_CENTERS use default_model_mod, only : adv_1step, nc_write_model_vars, & init_time => fail_init_time, & - init_conditions => fail_init_conditions + init_conditions => fail_init_conditions, & + parse_variables_clamp, & + MAX_STATE_VARIABLE_FIELDS_CLAMP use cam_common_code_mod, only : above_ramp_start, are_damping, build_cam_pressure_columns, build_heights, & cam_grid, cdebug_level, check_good_levels, cno_normalization_of_scale_heights, & @@ -94,8 +96,7 @@ module model_mod set_vert_localization, vert_interp, vertical_localization_type, write_model_time use cam_common_code_mod, only : nc_write_model_atts, grid_data, read_grid_info, & - set_cam_variable_info, MAX_STATE_VARIABLES, & - num_state_table_columns, MAX_PERT, & + MAX_PERT, & shortest_time_between_assimilations, domain_id, & cuse_log_vertical_scale, & cno_normalization_of_scale_heights, & @@ -180,8 +181,7 @@ module model_mod ! for no clamping, use the string 'NA' ! to have the assimilation change the variable use 'UPDATE', else 'NO_UPDATE' -character(len=vtablenamelength) :: state_variables(MAX_STATE_VARIABLES * & - num_state_table_columns ) = ' ' +character(len=vtablenamelength) :: state_variables(MAX_STATE_VARIABLE_FIELDS_CLAMP) = ' ' namelist /model_nml/ & cam_template_filename, & @@ -296,9 +296,10 @@ subroutine static_init_model() ! initialize global values that are used frequently call init_globals() -! read the namelist &model_nml :: state_variables -! to set up what will be read into the cam state vector -call set_cam_variable_info(cam_template_filename, state_variables) +! parse_variables converts the character table that was read in from +! model_nml:state_variables into a state_var_type that can be +! passed to add_domain +domain_id = add_domain(cam_template_filename, parse_variables_clamp(state_variables)) ! Verify that required variables are in the state vector. call verify_state_var_list diff --git a/models/cam-se/model_mod.f90 b/models/cam-se/model_mod.f90 index a51a8eaade..c2a228916b 100644 --- a/models/cam-se/model_mod.f90 +++ b/models/cam-se/model_mod.f90 @@ -82,7 +82,9 @@ module model_mod get_molar_mass, get_volume_mixing_ratio use default_model_mod, only : adv_1step, nc_write_model_vars, & init_time => fail_init_time, & - init_conditions => fail_init_conditions + init_conditions => fail_init_conditions, & + parse_variables_clamp, & + MAX_STATE_VARIABLE_FIELDS_CLAMP use cam_common_code_mod, only : above_ramp_start, are_damping, build_cam_pressure_columns, build_heights, & cam_grid, cdebug_level, check_good_levels, cno_normalization_of_scale_heights, & @@ -96,8 +98,7 @@ module model_mod use cam_common_code_mod, only : nc_write_model_atts, grid_data, read_grid_info, & - set_cam_variable_info, MAX_STATE_VARIABLES, & - num_state_table_columns, MAX_PERT, & + MAX_PERT, & shortest_time_between_assimilations, domain_id, & cuse_log_vertical_scale, & cno_normalization_of_scale_heights, & @@ -186,8 +187,7 @@ module model_mod ! for no clamping, use the string 'NA' ! to have the assimilation change the variable use 'UPDATE', else 'NO_UPDATE' -character(len=vtablenamelength) :: state_variables(MAX_STATE_VARIABLES * & - num_state_table_columns ) = ' ' +character(len=vtablenamelength) :: state_variables(MAX_STATE_VARIABLE_FIELDS_CLAMP) = ' ' namelist /model_nml/ & dry_mass_vertical_coordinate, & @@ -341,9 +341,11 @@ subroutine static_init_model() ! initialize global values that are used frequently call init_globals() -! read the namelist &model_nml :: state_variables -! to set up what will be read into the cam state vector -call set_cam_variable_info(cam_template_filename, state_variables) +! parse_variables converts the character table that was read in from +! model_nml:state_variables and returns a state_var_type +! that can be passed to add_domain +domain_id = add_domain(cam_template_filename, parse_variables_clamp(state_variables)) + call verify_state_var_list diff --git a/models/cice/model_mod.f90 b/models/cice/model_mod.f90 index 4346d584aa..08e838ff43 100644 --- a/models/cice/model_mod.f90 +++ b/models/cice/model_mod.f90 @@ -1,8 +1,6 @@ ! DART software - Copyright UCAR. This open source software is provided ! by UCAR, "as is", without charge, subject to all terms of use at ! http://www.image.ucar.edu/DAReS/DART/DART_download -! -! $Id$ module model_mod @@ -37,7 +35,8 @@ module model_mod nc_write_location use default_model_mod, only : init_time, init_conditions, adv_1step, & - nc_write_model_vars + nc_write_model_vars, parse_variables, & + MAX_STATE_VARIABLE_FIELDS use utilities_mod, only : register_module, error_handler, & E_ERR, E_MSG, nmlfileunit, get_unit, & @@ -115,7 +114,8 @@ module model_mod use state_structure_mod, only : add_domain, get_model_variable_indices, & get_num_variables, get_index_start, & - get_num_dims, get_domain_size, state_structure_info + get_num_dims, get_domain_size, state_structure_info, & + get_varid_from_kind, get_kind_index use typesizes use netcdf @@ -167,24 +167,12 @@ module model_mod ! Storage for a random sequence for perturbing a single initial state type(random_seq_type) :: random_seq -! DART state vector contents are specified in the input.nml:&model_nml namelist. -integer, parameter :: max_state_variables = 10 -integer, parameter :: num_state_table_columns = 3 -character(len=NF90_MAX_NAME) :: variable_table( max_state_variables, num_state_table_columns ) -integer :: state_kinds_list( max_state_variables ) -logical :: update_var_list( max_state_variables ) - -! identifiers for variable_table -integer, parameter :: VAR_NAME_INDEX = 1 -integer, parameter :: VAR_QTY_INDEX = 2 -integer, parameter :: VAR_UPDATE_INDEX = 3 - ! things which can/should be in the model_nml integer :: assimilation_period_days = 1 integer :: assimilation_period_seconds = 0 real(r8) :: model_perturbation_amplitude = 0.2 logical :: update_dry_cell_walls = .false. -character(len=metadatalength) :: model_state_variables(max_state_variables * num_state_table_columns ) = ' ' +character(len=metadatalength) :: model_state_variables(MAX_STATE_VARIABLE_FIELDS) = ' ' integer :: debug = 0 ! turn up for more and more debug messages ! valid values: native, big_endian, little_endian @@ -241,9 +229,6 @@ module model_mod ! is not being used. !------------------------------------------------------------------ -! Number of fields in the state vector -integer :: nfields - ! Grid parameters - the values will be read from a ! standard cice namelist and filled in here. @@ -410,11 +395,6 @@ subroutine static_init_model() if (debug > 2) call write_grid_netcdf() ! DEBUG only if (debug > 2) call write_grid_interptest() ! DEBUG only -! verify that the model_state_variables namelist was filled in correctly. -! returns variable_table which has variable names, kinds and update strings. -call verify_state_variables(model_state_variables, nfields, variable_table, & - state_kinds_list, update_var_list) - ! in spite of the staggering, all grids are the same size ! and offset by half a grid cell. 4 are 3D and 2 are 2D. ! e.g. aicen,vicen,vsnon = 256 x 225 x 5 @@ -428,13 +408,19 @@ subroutine static_init_model() ! Initialize the interpolation routines call init_interp() +if ( model_state_variables(1) == ' ' ) then ! no model_state_variables namelist provided + call use_default_state_variables( model_state_variables ) + string1 = 'model_nml:model_state_variables not specified - using default variables' + call error_handler(E_MSG,'static_init_model',string1,source,revision,revdate) +endif + ! Determine the shape of the variables from "cice.r.nc" ! The assimilate.csh, perfect_model.csh must ensure the cice restart file ! is linked to this filename. -domain_id = add_domain('cice.r.nc', nfields, & - var_names = variable_table(1:nfields, VAR_NAME_INDEX), & - kind_list = state_kinds_list(1:nfields), & - update_list = update_var_list(1:nfields)) +! parse_variables converts the character table that was read in from +! model_nml:model_state_variables and returns a state_var_type that can be +! passed to add_domain +domain_id = add_domain('cice.r.nc', parse_variables(model_state_variables)) if (debug > 2) call state_structure_info(domain_id) @@ -885,25 +871,25 @@ subroutine model_interpolate(state_handle, ens_size, location, obs_type, expecte SELECT CASE (obs_type) CASE (QTY_SEAICE_AGREG_THICKNESS ) ! these kinds require aggregating 3D vars to make a 2D var cat_signal = -1 ! for extra special procedure to aggregate - base_offset = get_index_start(domain_id, get_varid_from_kind(QTY_SEAICE_VOLUME)) + base_offset = get_index_start(domain_id, get_varid_from_kind(domain_id, QTY_SEAICE_VOLUME)) CASE (QTY_SEAICE_AGREG_SNOWDEPTH ) ! these kinds require aggregating 3D vars to make a 2D var cat_signal = -1 ! for extra special procedure to aggregate - base_offset = get_index_start(domain_id, get_varid_from_kind(QTY_SEAICE_SNOWVOLUME)) + base_offset = get_index_start(domain_id, get_varid_from_kind(domain_id, QTY_SEAICE_SNOWVOLUME)) CASE (QTY_SEAICE_AGREG_CONCENTR ) ! these kinds require aggregating a 3D var to make a 2D var cat_signal = 0 ! for aggregate variable, send signal to lon_lat_interp - base_offset = get_index_start(domain_id, get_varid_from_kind(QTY_SEAICE_CONCENTR)) + base_offset = get_index_start(domain_id, get_varid_from_kind(domain_id, QTY_SEAICE_CONCENTR)) CASE (QTY_SEAICE_AGREG_VOLUME ) ! these kinds require aggregating a 3D var to make a 2D var cat_signal = 0 ! for aggregate variable, send signal to lon_lat_interp - base_offset = get_index_start(domain_id, get_varid_from_kind(QTY_SEAICE_VOLUME)) + base_offset = get_index_start(domain_id, get_varid_from_kind(domain_id, QTY_SEAICE_VOLUME)) CASE (QTY_SEAICE_AGREG_SNOWVOLUME ) ! these kinds require aggregating a 3D var to make a 2D var cat_signal = 0 ! for aggregate variable, send signal to lon_lat_interp - base_offset = get_index_start(domain_id, get_varid_from_kind(QTY_SEAICE_SNOWVOLUME)) + base_offset = get_index_start(domain_id, get_varid_from_kind(domain_id, QTY_SEAICE_SNOWVOLUME)) CASE (QTY_SEAICE_AGREG_SURFACETEMP) ! FEI need aicen to average the temp, have not considered open water temp yet cat_signal = -3 - base_offset = get_index_start(domain_id, get_varid_from_kind(QTY_SEAICE_SURFACETEMP)) + base_offset = get_index_start(domain_id, get_varid_from_kind(domain_id, QTY_SEAICE_SURFACETEMP)) CASE (QTY_SOM_TEMPERATURE) ! these kinds are 1d variables cat_signal = 3 - base_offset = get_index_start(domain_id,get_varid_from_kind(QTY_SOM_TEMPERATURE)) + base_offset = get_index_start(domain_id,get_varid_from_kind(domain_id, QTY_SOM_TEMPERATURE)) CASE (QTY_SEAICE_CONCENTR , & ! these kinds have an additional dim for category QTY_SEAICE_FY , & QTY_SEAICE_VOLUME , & @@ -938,7 +924,7 @@ subroutine model_interpolate(state_handle, ens_size, location, obs_type, expecte QTY_SEAICE_SNOWENTHALPY003 ) ! move pointer to the particular category ! then treat as 2d field in lon_lat_interp - base_offset = get_index_start(domain_id, get_varid_from_kind(obs_type)) + base_offset = get_index_start(domain_id, get_varid_from_kind(domain_id, obs_type)) base_offset = base_offset + (cat_index-1) * Nx * Ny cat_signal = 1 ! now same as boring 2d field CASE ( QTY_U_SEAICE_COMPONENT , & ! these kinds are just 2D vars @@ -947,7 +933,7 @@ subroutine model_interpolate(state_handle, ens_size, location, obs_type, expecte QTY_SEAICE_ALBEDODIRNIR , & QTY_SEAICE_ALBEDOINDVIZ , & QTY_SEAICE_ALBEDOINDNIR ) - base_offset = get_index_start(domain_id, get_varid_from_kind(obs_type)) + base_offset = get_index_start(domain_id, get_varid_from_kind(domain_id, obs_type)) cat_signal = 2 ! also boring 2d field (treat same as cat_signal 1) CASE DEFAULT ! Not a legal type for interpolation, return istatus error @@ -965,12 +951,12 @@ subroutine model_interpolate(state_handle, ens_size, location, obs_type, expecte do icat = 1,Ncat !reads in aicen cat_signal_interm = 1 - base_offset = get_index_start(domain_id,get_varid_from_kind(QTY_SEAICE_CONCENTR)) + base_offset = get_index_start(domain_id,get_varid_from_kind(domain_id, QTY_SEAICE_CONCENTR)) base_offset = base_offset + (icat-1) * Nx * Ny call lon_lat_interpolate(state_handle, ens_size, base_offset, llon, llat, obs_type, cat_signal_interm, expected_conc, istatus) !reads in fyn cat_signal_interm = 1 - base_offset = get_index_start(domain_id,get_varid_from_kind(QTY_SEAICE_FY)) + base_offset = get_index_start(domain_id,get_varid_from_kind(domain_id, QTY_SEAICE_FY)) base_offset = base_offset + (icat-1) * Nx * Ny call lon_lat_interpolate(state_handle, ens_size, base_offset, llon, llat, obs_type, cat_signal_interm, expected_fy, istatus) temp = temp + expected_conc * expected_fy !sum(aicen*fyn) = FY % over ice @@ -993,12 +979,12 @@ subroutine model_interpolate(state_handle, ens_size, location, obs_type, expecte do icat = 1,Ncat !reads in aicen cat_signal_interm = 1 - base_offset = get_index_start(domain_id,get_varid_from_kind(QTY_SEAICE_CONCENTR)) + base_offset = get_index_start(domain_id,get_varid_from_kind(domain_id, QTY_SEAICE_CONCENTR)) base_offset = base_offset + (icat-1) * Nx * Ny call lon_lat_interpolate(state_handle, ens_size, base_offset, llon, llat, obs_type, cat_signal_interm, expected_conc, istatus) !reads in Tsfcn cat_signal_interm = 1 - base_offset = get_index_start(domain_id,get_varid_from_kind(QTY_SEAICE_SURFACETEMP)) + base_offset = get_index_start(domain_id,get_varid_from_kind(domain_id, QTY_SEAICE_SURFACETEMP)) base_offset = base_offset + (icat-1) * Nx * Ny call lon_lat_interpolate(state_handle, ens_size, base_offset, llon, llat, obs_type, cat_signal_interm, expected_tsfc, istatus) if (any(expected_conc<0.0) .or. any(expected_conc>1.0))then @@ -1036,7 +1022,7 @@ subroutine model_interpolate(state_handle, ens_size, location, obs_type, expecte if (cat_signal == -1) then ! we need to know the aggregate sea ice concentration for these special cases - base_offset = get_index_start(domain_id, get_varid_from_kind(QTY_SEAICE_CONCENTR)) + base_offset = get_index_start(domain_id, get_varid_from_kind(domain_id, QTY_SEAICE_CONCENTR)) call lon_lat_interpolate(state_handle, ens_size, base_offset, llon, llat, obs_type, cat_signal, expected_aggr_conc, istatus) expected_obs = expected_obs/max(expected_aggr_conc,1.0e-8) ! hope this is allowed so we never divide by zero @@ -1832,7 +1818,7 @@ subroutine get_state_meta_data(index_in, location, var_type) if ( .not. module_initialized ) call static_init_model call get_model_variable_indices(index_in, lon_index, lat_index, cat_index, var_id=var_id) -call get_state_kind(var_id, local_var) +local_var = get_kind_index(domain_id, var_id) if (is_on_ugrid(local_var)) then lon = ULON(lon_index, lat_index) @@ -1857,52 +1843,6 @@ end subroutine get_state_meta_data !-------------------------------------------------------------------- -function get_varid_from_kind(dart_kind) - -integer, intent(in) :: dart_kind -integer :: get_varid_from_kind - -! given a kind, return what variable number it is - -integer :: i - -do i = 1, get_num_variables(domain_id) - if (dart_kind == state_kinds_list(i)) then - get_varid_from_kind = i - return - endif -end do - -if (debug > 1) then - write(string1, *) 'Kind ', dart_kind, ' not found in state vector' - write(string2, *) 'AKA ', get_name_for_quantity(dart_kind), ' not found in state vector' - call error_handler(E_MSG,'get_varid_from_kind', string1, & - source, revision, revdate, text2=string2) -endif - -get_varid_from_kind = -1 - -end function get_varid_from_kind - - -!------------------------------------------------------------------ - -subroutine get_state_kind(var_ind, var_type) - integer, intent(in) :: var_ind - integer, intent(out) :: var_type - -! Given an integer index into the state vector structure, returns the kind, -! and both the starting offset for this kind, as well as the offset into -! the block of this kind. - -if ( .not. module_initialized ) call static_init_model - -var_type = state_kinds_list(var_ind) - -end subroutine get_state_kind - -!------------------------------------------------------------------ - subroutine get_state_kind_inc_dry(index_in, var_type) integer(i8), intent(in) :: index_in integer, intent(out) :: var_type @@ -1915,7 +1855,7 @@ subroutine get_state_kind_inc_dry(index_in, var_type) if ( .not. module_initialized ) call static_init_model call get_model_variable_indices(index_in, lon_index, lat_index, depth_index, var_id=var_id) -call get_state_kind(var_id, var_type) +var_type = get_kind_index(domain_id, var_id) ! if on land, replace type with dry land. if(is_dry_land(var_type, lon_index, lat_index)) then @@ -2587,94 +2527,6 @@ function read_model_time(filename) end function read_model_time -!------------------------------------------------------------------ -!> Verify that the namelist was filled in correctly, and check -!> that there are valid entries for the dart_kind. -!> Returns a table with columns: -!> -!> netcdf_variable_name ; dart_kind_string ; update_string -!> - -subroutine verify_state_variables( state_variables, ngood, table, kind_list, update_var ) - -character(len=*), intent(inout) :: state_variables(:) -integer, intent(out) :: ngood -character(len=*), intent(out) :: table(:,:) -integer, intent(out) :: kind_list(:) ! kind number -logical, optional, intent(out) :: update_var(:) ! logical update - -integer :: nrows, i -character(len=NF90_MAX_NAME) :: varname, dartstr, update - -if ( .not. module_initialized ) call static_init_model - -nrows = size(table,1) - -ngood = 0 - -!>@todo deprecate. Remove a hidden 'default' set of variables. -!>@ The default is provided in the input namelist. - -if ( state_variables(1) == ' ' ) then ! no model_state_variables namelist provided - call use_default_state_variables( state_variables ) - string1 = 'model_nml:model_state_variables not specified using default variables' - call error_handler(E_MSG,'verify_state_variables',string1,source,revision,revdate) -endif - -MyLoop : do i = 1, nrows - - varname = trim(state_variables(3*i -2)) - dartstr = trim(state_variables(3*i -1)) - update = trim(state_variables(3*i )) - - call to_upper(update) - - table(i,1) = trim(varname) - table(i,2) = trim(dartstr) - table(i,3) = trim(update) - - if ( table(i,1) == ' ' .and. table(i,2) == ' ' .and. table(i,3) == ' ') exit MyLoop - - if ( table(i,1) == ' ' .or. table(i,2) == ' ' .or. table(i,3) == ' ' ) then - string1 = 'model_nml:model_state_variables not fully specified' - call error_handler(E_ERR,'verify_state_variables',string1,source,revision,revdate) - endif - - ! Make sure DART kind is valid - - kind_list(i) = get_index_for_quantity(dartstr) - if( kind_list(i) < 0 ) then - write(string1,'(''there is no obs_kind <'',a,''> in obs_kind_mod.f90'')') trim(dartstr) - call error_handler(E_ERR,'verify_state_variables',string1,source,revision,revdate) - endif - - ! Make sure the update variable has a valid name - - if ( present(update_var) )then - SELECT CASE (update) - CASE ('UPDATE') - update_var(i) = .true. - CASE ('NO_COPY_BACK') - update_var(i) = .false. - CASE DEFAULT - write(string1,'(A)') 'only UPDATE or NO_COPY_BACK supported in model_state_variable namelist' - write(string2,'(6A)') 'you provided : ', trim(varname), ', ', trim(dartstr), ', ', trim(update) - call error_handler(E_ERR,'verify_state_variables',string1,source,revision,revdate, text2=string2) - END SELECT - endif - - ! Record the contents of the DART state vector - - if (do_output()) then - write(string1,'(A,I2,6A)') 'variable ',i,' is ',trim(varname), ', ', trim(dartstr), ', ', trim(update) - call error_handler(E_MSG,'verify_state_variables',string1,source,revision,revdate) - endif - - ngood = ngood + 1 -enddo MyLoop - -end subroutine verify_state_variables - !------------------------------------------------------------------ !> Default state_variables from model_mod. !>@todo DEPRECATE @@ -2684,7 +2536,7 @@ subroutine use_default_state_variables( state_variables ) character(len=*), intent(inout) :: state_variables(:) ! strings must all be the same length for the gnu compiler -state_variables( 1:5*num_state_table_columns ) = & +state_variables( 1:15 ) = & (/ 'CONCENTRATION ', 'QTY_SEAICE_CONCENTR ', 'UPDATE ', & 'ICEVOLUME ', 'QTY_SEAICE_VOLUME ', 'UPDATE ', & 'SNOWVOLUME ', 'QTY_SEAICE_SNOWVOLUME ', 'UPDATE ', & diff --git a/models/template/threed_model_mod.f90 b/models/template/threed_model_mod.f90 index 07976db5de..335e1d582a 100644 --- a/models/template/threed_model_mod.f90 +++ b/models/template/threed_model_mod.f90 @@ -9,7 +9,7 @@ module model_mod ! with the DART data assimilation infrastructure. Do not change the arguments ! for the public routines. -use types_mod, only : r8, i8, MISSING_R8 +use types_mod, only : r8, i8, MISSING_R8, vtablenamelength use time_manager_mod, only : time_type, set_time @@ -38,7 +38,8 @@ module model_mod use default_model_mod, only : pert_model_copies, read_model_time, write_model_time, & init_time => fail_init_time, & init_conditions => fail_init_conditions, & - convert_vertical_obs, convert_vertical_state, adv_1step + convert_vertical_obs, convert_vertical_state, adv_1step, & + parse_variables, MAX_STATE_VARIABLE_FIELDS implicit none private @@ -74,8 +75,9 @@ module model_mod character(len=256) :: template_file = 'model_restart.nc' integer :: time_step_days = 0 integer :: time_step_seconds = 3600 +character(len=vtablenamelength) :: state_variables(MAX_STATE_VARIABLE_FIELDS) = ' ' -namelist /model_nml/ template_file, time_step_days, time_step_seconds +namelist /model_nml/ template_file, time_step_days, time_step_seconds, state_variables contains @@ -108,9 +110,11 @@ subroutine static_init_model() assimilation_time_step = set_time(time_step_seconds, & time_step_days) - -! Define which variables are in the model state -dom_id = add_domain(template_file, num_vars=2, var_names=(/'Temp', 'Wind'/)) +! Define which variables are in the model state; +! parse_variables converts the character table that was read in from +! model_nml:model_state_variables to a state_var_type that can be passed +! to add_domain +dom_id = add_domain(template_file, parse_variables(state_variables)) end subroutine static_init_model diff --git a/models/template/work/threed_input.nml b/models/template/work/threed_input.nml index ca07cf2034..8c61116ac5 100644 --- a/models/template/work/threed_input.nml +++ b/models/template/work/threed_input.nml @@ -127,6 +127,9 @@ &model_nml time_step_days = 0, time_step_seconds = 3600 + state_variables = 'SALT_CUR ', 'QTY_SALINITY', 'UPDATE', + 'TEMP_CUR ', 'QTY_POTENTIAL_TEMPERATURE', 'UPDATE', + 'UVEL_CUR ', 'QTY_U_CURRENT_COMPONENT ', 'UPDATE', / &utilities_nml diff --git a/models/utilities/default_model_mod.f90 b/models/utilities/default_model_mod.f90 index d4046eb88e..49f9a0458c 100644 --- a/models/utilities/default_model_mod.f90 +++ b/models/utilities/default_model_mod.f90 @@ -19,14 +19,16 @@ module default_model_mod use utilities_mod, only : error_handler, E_ERR, E_MSG, nmlfileunit, & do_output, find_namelist_in_file, check_namelist_read, & - do_nml_file, do_nml_term + do_nml_file, do_nml_term, to_upper -use netcdf_utilities_mod, only : nc_check +use netcdf_utilities_mod, only : nc_check, NF90_MAX_NAME use ensemble_manager_mod, only : ensemble_type use dart_time_io_mod, only : read_model_time, write_model_time +use obs_kind_mod, only : get_index_for_quantity + implicit none private @@ -49,8 +51,25 @@ module default_model_mod convert_vertical_obs, & convert_vertical_state, & read_model_time, & ! from the dart_time_io module - write_model_time - + write_model_time, & + parse_variables, & + parse_variables_clamp, & + state_var_type, & + MAX_STATE_VARIABLE_FIELDS, & + MAX_STATE_VARIABLE_FIELDS_CLAMP + +type :: state_var_type + integer :: nvars = -1 + character(len=NF90_MAX_NAME), allocatable :: netcdf_var_names(:) + integer, allocatable :: qtys(:) + real(r8), allocatable :: clamp_values(:, :) + logical, allocatable :: updates(:) +end type state_var_type + +! note _fields is 3* and _fields_clamp is 5*MAX_STATE_VARIABLES +integer, parameter :: MAX_STATE_VARIABLES = 1000 +integer, parameter :: MAX_STATE_VARIABLE_FIELDS = 3000 +integer, parameter :: MAX_STATE_VARIABLE_FIELDS_CLAMP = 5000 character(len=*), parameter :: source = 'utilities/default_model_mod.f90' contains @@ -82,6 +101,7 @@ subroutine init_conditions(x) end subroutine init_conditions +!------------------------------------------------------------------ !------------------------------------------------------------------ subroutine fail_init_conditions(x) @@ -283,6 +303,167 @@ subroutine pert_model_copies(state_ens_handle, ens_size, pert_amp, interf_provid end subroutine pert_model_copies +!-------------------------------------------------------------------- + +!> Parses the character table that was read in from +!> model_nml:model_state_variables and returns a +!> state_var_type state_vars with nvars ; netcdf variable names ; +!> qtys (kinds) ; clamp values ; updates +!> that there are valid entries for the dart_kind. +! +!> Verifies that the namelist entry was filled in correctly, and checks +!> that there are valid entries for the dart_kind. + +function parse_variables_clamp(vars_table) result(state_vars) + +character(len=*), intent(in) :: vars_table(MAX_STATE_VARIABLE_FIELDS_CLAMP) +type(state_var_type) :: state_vars + +character(len=NF90_MAX_NAME) :: netcdf_var_name, dart_qty_str, update +character(len=256) :: string1, string2 +integer :: i, ivar + +! Loop through the variables array to get the actual count of the number of variables +do ivar = 1, MAX_STATE_VARIABLES + ! If the first element in the row is an empty string, the loop has exceeded the extent of the variables + if (vars_table(5*ivar-4) == '') then + state_vars%nvars = ivar-1 + exit + endif +enddo + +if (state_vars%nvars >= MAX_STATE_VARIABLES -1) then + write(string1,*) 'nvars ', state_vars%nvars, ' >= MAX_STATE_VARIABLES-1', MAX_STATE_VARIABLES-1, 'increase MAX_STATE_VARIABLES' + call error_handler(E_ERR, string1, source) +endif + +! Allocate the arrays in the var derived type +allocate(state_vars%netcdf_var_names(state_vars%nvars), state_vars%qtys(state_vars%nvars), state_vars%clamp_values(state_vars%nvars, 2), state_vars%updates(state_vars%nvars)) + +RowsLoop : do i = 1, state_vars%nvars + + netcdf_var_name = trim(vars_table(5*i-4)) + state_vars%netcdf_var_names(i) = trim(netcdf_var_name) + + dart_qty_str = trim(vars_table(5*i-3)) + call to_upper(dart_qty_str) + ! Make sure DART qty is valid + state_vars%qtys(i) = get_index_for_quantity(dart_qty_str) + if( state_vars%qtys(i) < 0 ) then + write(string1,*) 'The quantity specified in the &model_nml "', trim(dart_qty_str), '", is not present in obs_kind_mod.f90' + call error_handler(E_ERR,'get_state_variables_clamp',string1) + endif + + if (vars_table(5*i-2) /= 'NA') then + read(vars_table(5*i-2), '(d16.8)') state_vars%clamp_values(i,1) + else + state_vars%clamp_values(i,1) = MISSING_R8 + endif + + if (vars_table(5*i-1) /= 'NA') then + read(vars_table(5*i-1), '(d16.8)') state_vars%clamp_values(i,2) + else + state_vars%clamp_values(i,2) = MISSING_R8 + endif + + update = trim(vars_table(5*i)) + call to_upper(update) + select case (update) + case ('UPDATE') + state_vars%updates(i) = .true. + case ('NO_COPY_BACK') + state_vars%updates(i) = .false. + case default + write(string1,'(A)') 'Invalid update variable in &model_nml:model_state_variable - only UPDATE or NO_COPY_BACK are supported' + write(string2,'(6A)') 'Issue: ', trim(netcdf_var_name), ', ', trim(dart_qty_str), ', ', trim(update) + call error_handler(E_ERR,'get_state_variables_clamp',string1, text2=string2) + end select + + ! Checking that the rows in the nml entry are all complete + if ( dart_qty_str == '' .or. vars_table(5*i-2) == '' .or. vars_table(5*i-1) == '' .or. update == '' ) then + string1 = 'model_nml:model_state_variables not fully specified' + call error_handler(E_ERR, 'get_state_variables_clamp', string1) + endif + +enddo RowsLoop + +end function parse_variables_clamp + +!-------------------------------------------------------------------- + +!> Parses the character table that was read in from +!> model_nml:model_state_variables and returns a +!> state_var_type state_vars with nvars ; netcdf variable names ; +!> qtys (kinds) ; updates +! +!> Verifies that the namelist entry was filled in correctly, and checks +!> that there are valid entries for the dart_kind. + +function parse_variables(vars_table) result(state_vars) + +character(len=*), intent(in) :: vars_table(MAX_STATE_VARIABLE_FIELDS) +type(state_var_type) :: state_vars + +character(len=NF90_MAX_NAME) :: netcdf_var_name, dart_qty_str, update +character(len=256) :: string1, string2 +integer :: i, ivar + +! Loop through the variables array to get the actual count of the number of variables +do ivar = 1, MAX_STATE_VARIABLES + ! If the first element in the row is an empty string, the loop has exceeded the extent of the variables + if (vars_table(3*ivar-2) == '') then + state_vars%nvars = ivar-1 + exit + endif +enddo + +if (state_vars%nvars >= MAX_STATE_VARIABLES -1) then + write(string1,*) 'nvars ', state_vars%nvars, ' >= MAX_STATE_VARIABLES-1', MAX_STATE_VARIABLES-1, 'increase MAX_STATE_VARIABLES' + call error_handler(E_ERR, string1, source) +endif + +! Allocate the arrays in the var derived type +allocate(state_vars%netcdf_var_names(state_vars%nvars), state_vars%qtys(state_vars%nvars), state_vars%updates(state_vars%nvars)) + +RowsLoop : do i = 1, state_vars%nvars + + netcdf_var_name = trim(vars_table(3*i-2)) + state_vars%netcdf_var_names(i) = trim(netcdf_var_name) + + update = trim(vars_table(3*i)) + call to_upper(update) + + dart_qty_str = trim(vars_table(3*i-1)) + call to_upper(dart_qty_str) + + ! Checking that the rows in the nml entry are all complete + if ( dart_qty_str == '' .or. update == '' ) then + string1 = 'model_nml:model_state_variables not fully specified' + call error_handler(E_ERR, 'parse_variables', string1) + endif + + ! Make sure DART qty is valid + state_vars%qtys(i) = get_index_for_quantity(dart_qty_str) + if( state_vars%qtys(i) < 0 ) then + write(string1,'(3A)') 'The quantity specified in the &model_nml "', trim(dart_qty_str), '", is not present in obs_kind_mod.f90' + call error_handler(E_ERR,'parse_variables',string1) + endif + + select case (update) + case ('UPDATE') + state_vars%updates(i) = .true. + case ('NO_COPY_BACK') + state_vars%updates(i) = .false. + case default + write(string1,'(A)') 'Invalid update variable in &model_nml:model_state_variable - only UPDATE or NO_COPY_BACK are supported' + write(string2,'(6A)') 'Issue: ', trim(netcdf_var_name), ', ', trim(dart_qty_str), ', ', trim(update) + call error_handler(E_ERR,'parse_variables',string1, text2=string2) + end select + +enddo RowsLoop + +end function parse_variables + !=================================================================== ! End of model_mod !=================================================================== diff --git a/models/wrf_hydro/model_mod.f90 b/models/wrf_hydro/model_mod.f90 index 8f943989ea..8dc7b54ba3 100644 --- a/models/wrf_hydro/model_mod.f90 +++ b/models/wrf_hydro/model_mod.f90 @@ -44,7 +44,8 @@ module model_mod use dart_time_io_mod, only : write_model_time -use default_model_mod, only : adv_1step, nc_write_model_vars +use default_model_mod, only : adv_1step, nc_write_model_vars, parse_variables_clamp, & + MAX_STATE_VARIABLE_FIELDS_CLAMP use noah_hydro_mod, only : configure_lsm, configure_hydro, & n_link, linkLong, linkLat, linkAlt, get_link_tree, & @@ -115,15 +116,6 @@ module model_mod integer(i8) :: model_size type(time_type) :: time_step -! 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_STATEINDX = 5 ! ... update (state) or not -integer, parameter :: MAX_STATE_VARIABLES = 40 -integer, parameter :: NUM_STATE_TABLE_COLUMNS = 5 - integer :: domain_count integer :: idom, idom_hydro = -1, idom_parameters = -1, idom_lsm = -1 @@ -140,9 +132,9 @@ module model_mod character(len=256) :: perturb_distribution = 'lognormal' character(len=obstypelength) :: & - lsm_variables( NUM_STATE_TABLE_COLUMNS,MAX_STATE_VARIABLES) = '', & - hydro_variables(NUM_STATE_TABLE_COLUMNS,MAX_STATE_VARIABLES) = '', & - parameters( NUM_STATE_TABLE_COLUMNS,MAX_STATE_VARIABLES) = '' + lsm_variables(MAX_STATE_VARIABLE_FIELDS_CLAMP) = '', & + hydro_variables(MAX_STATE_VARIABLE_FIELDS_CLAMP) = '', & + parameters(MAX_STATE_VARIABLE_FIELDS_CLAMP) = '' namelist /model_nml/ assimilation_period_days, & assimilation_period_seconds, & @@ -177,16 +169,8 @@ subroutine static_init_model() character(len=*), parameter :: routine = 'static_init_model' integer :: iunit, io, domainID -integer :: n_lsm_fields -integer :: n_hydro_fields -integer :: n_parameters integer :: vsize -character(len=obstypelength) :: var_names(MAX_STATE_VARIABLES) -real(r8) :: var_ranges(MAX_STATE_VARIABLES,2) -logical :: var_update(MAX_STATE_VARIABLES) -integer :: var_qtys( MAX_STATE_VARIABLES) - character(len=256) :: domain_name if ( module_initialized ) return ! only need to do this once @@ -250,13 +234,9 @@ subroutine static_init_model() call configure_hydro() call read_hydro_global_atts(domain_shapefiles(domainID)) - call verify_variables(hydro_variables, domain_shapefiles(domainID), & - n_hydro_fields, var_names, var_qtys, var_ranges, var_update) + idom_hydro = add_domain(domain_shapefiles(domainID), & - n_hydro_fields, var_names, & - kind_list=var_qtys, & - clamp_vals=var_ranges(1:n_hydro_fields,:), & - update_list=var_update) + parse_variables_clamp(hydro_variables)) if (debug > 99) call state_structure_info(idom_hydro) @@ -273,13 +253,9 @@ subroutine static_init_model() elseif (index(domain_name,'PARAMETER') > 0) then - call verify_variables(parameters, domain_shapefiles(domainID), n_parameters, & - var_names, var_qtys, var_ranges, var_update) idom_parameters = add_domain(domain_shapefiles(domainID), & - n_parameters, var_names, & - kind_list=var_qtys, & - clamp_vals=var_ranges(1:n_parameters,:), & - update_list=var_update ) + parse_variables_clamp(parameters)) + if (debug > 99) call state_structure_info(idom_parameters) !>@todo check the size of the parameter variables against nlinks @@ -288,14 +264,10 @@ subroutine static_init_model() call configure_lsm(lsm_model_choice) call read_noah_global_atts(domain_shapefiles(domainID)) - call verify_variables(lsm_variables, domain_shapefiles(domainID), n_lsm_fields, & - var_names, var_qtys, var_ranges, var_update) - idom_lsm = add_domain(domain_shapefiles(domainID), & - n_lsm_fields, var_names, & - kind_list=var_qtys, & - clamp_vals=var_ranges(1:n_lsm_fields,:), & - update_list=var_update) - if (debug > 99) call state_structure_info(idom_lsm) + idom_hydro = add_domain(domain_shapefiles(domainID), & + parse_variables_clamp(lsm_variables)) + + if (debug > 99) call state_structure_info(idom_lsm) else @@ -1252,90 +1224,6 @@ end subroutine end_model ! End of the required interfaces !======================================================================= -!----------------------------------------------------------------------- -!> given the list of variables and a filename, check user input -!> return the handle to the open netCDF file and the number of variables -!> in this 'domain' - -subroutine verify_variables( variable_table, filename, ngood, & - var_names, var_qtys, var_ranges, var_update) - -character(len=*), intent(in) :: variable_table(:,:) -character(len=*), intent(in) :: filename -integer, intent(out) :: ngood -character(len=*), intent(out) :: var_names(:) -real(r8), intent(out) :: var_ranges(:,:) -logical, intent(out) :: var_update(:) -integer , intent(out) :: var_qtys(:) - -character(len=*), parameter :: routine = 'verify_variables' - -integer :: io, i, quantity -real(r8) :: minvalue, maxvalue - -character(len=NF90_MAX_NAME) :: varname -character(len=NF90_MAX_NAME) :: dartstr -character(len=NF90_MAX_NAME) :: minvalstring -character(len=NF90_MAX_NAME) :: maxvalstring -character(len=NF90_MAX_NAME) :: state_or_aux - -ngood = 0 -MyLoop : do i = 1, size(variable_table,2) - - varname = variable_table(VT_VARNAMEINDX,i) - dartstr = variable_table(VT_KINDINDX ,i) - minvalstring = variable_table(VT_MINVALINDX ,i) - maxvalstring = variable_table(VT_MAXVALINDX ,i) - state_or_aux = variable_table(VT_STATEINDX ,i) - - if ( varname == ' ' .and. dartstr == ' ' ) exit MyLoop ! Found end of list. - - if ( varname == ' ' .or. dartstr == ' ' ) then - string1 = 'model_nml: variable list not fully specified' - string2 = 'reading from "'//trim(filename)//'"' - call error_handler(E_ERR,routine, string1, & - source, revision, revdate, text2=string2) - endif - - ! The internal DART routines check if the variable name is valid. - - ! Make sure DART kind is valid - quantity = get_index_for_quantity(dartstr) - if( quantity < 0 ) then - write(string1,'(''there is no obs_kind "'',a,''" in obs_kind_mod.f90'')') & - trim(dartstr) - call error_handler(E_ERR,routine,string1,source,revision,revdate) - endif - - ! All good to here - fill the output variables - - ngood = ngood + 1 - var_names( ngood) = varname - var_qtys( ngood) = quantity - var_ranges(ngood,:) = (/ MISSING_R8, MISSING_R8 /) - var_update(ngood) = .false. ! at least initially - - ! convert the [min,max]valstrings to numeric values if possible - read(minvalstring,*,iostat=io) minvalue - if (io == 0) var_ranges(ngood,1) = minvalue - - read(maxvalstring,*,iostat=io) maxvalue - if (io == 0) var_ranges(ngood,2) = maxvalue - - call to_upper(state_or_aux) - if (state_or_aux == 'UPDATE') var_update(ngood) = .true. - -enddo MyLoop - -if (ngood == MAX_STATE_VARIABLES) then - string1 = 'WARNING: you may need to increase "MAX_STATE_VARIABLES"' - write(string2,'(''you have specified at least '',i4,'' perhaps more.'')') ngood - call error_handler(E_MSG,routine,string1,source,revision,revdate,text2=string2) -endif - -end subroutine verify_variables - - !----------------------------------------------------------------------- !> Sets the location information arrays for each domain !> Each location array is declared to be 3D to make it easy to use