diff --git a/models/mpas_atm/MPAS_model_modules/LICENSE b/models/mpas_atm/MPAS_model_modules/LICENSE new file mode 100644 index 0000000000..c8060c7f24 --- /dev/null +++ b/models/mpas_atm/MPAS_model_modules/LICENSE @@ -0,0 +1,39 @@ +Copyright (c) 2013-2019, Los Alamos National Security, LLC (LANS) (Ocean: LA-CC-13-047; +Land Ice: LA-CC-13-117) and the University Corporation for Atmospheric Research (UCAR). + +All rights reserved. + +LANS is the operator of the Los Alamos National Laboratory under Contract No. +DE-AC52-06NA25396 with the U.S. Department of Energy. UCAR manages the National +Center for Atmospheric Research under Cooperative Agreement ATM-0753581 with the +National Science Foundation. The U.S. Government has rights to use, reproduce, +and distribute this software. NO WARRANTY, EXPRESS OR IMPLIED IS OFFERED BY +LANS, UCAR OR THE GOVERNMENT AND NONE OF THEM ASSUME ANY LIABILITY FOR THE USE +OF THIS SOFTWARE. If software is modified to produce derivative works, such +modified software should be clearly marked, so as not to confuse it with the +version available from LANS and UCAR. + +Additionally, redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +1) Redistributions of source code must retain the above copyright notice, this +list of conditions and the following disclaimer. + +2) Redistributions in binary form must reproduce the above copyright notice, +this list of conditions and the following disclaimer in the documentation and/or +other materials provided with the distribution. + +3) None of the names of LANS, UCAR or the names of its contributors, if any, may +be used to endorse or promote products derived from this software without +specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR +ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/models/mpas_atm/get_coeff_mod.f90 b/models/mpas_atm/MPAS_model_modules/get_coeff_mod.f90 similarity index 98% rename from models/mpas_atm/get_coeff_mod.f90 rename to models/mpas_atm/MPAS_model_modules/get_coeff_mod.f90 index a1a17c0760..8b229d880b 100644 --- a/models/mpas_atm/get_coeff_mod.f90 +++ b/models/mpas_atm/MPAS_model_modules/get_coeff_mod.f90 @@ -1,8 +1,3 @@ -! This code may (or may not) be part of the MPAS distribution, -! So it is not protected by the DART copyright agreement. -! -! DART $Id$ - module get_coeff_mod use types_mod, only : r8 diff --git a/models/mpas_atm/get_geometry_mod.f90 b/models/mpas_atm/MPAS_model_modules/get_geometry_mod.f90 similarity index 93% rename from models/mpas_atm/get_geometry_mod.f90 rename to models/mpas_atm/MPAS_model_modules/get_geometry_mod.f90 index 8d0e44b566..abd8c387aa 100644 --- a/models/mpas_atm/get_geometry_mod.f90 +++ b/models/mpas_atm/MPAS_model_modules/get_geometry_mod.f90 @@ -1,15 +1,9 @@ -! This code may (or may not) be part of the MPAS distribution, -! So it is not protected by the DART copyright agreement. -! -! DART $Id$ - module get_geometry_mod use types_mod, only : r8 implicit none private - save !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Purpose: perform interpolation of scalar and vector functions in 2D @@ -105,9 +99,3 @@ subroutine mpas_cross_product_in_r3(p_1,p_2,p_out) end subroutine mpas_cross_product_in_r3 end module get_geometry_mod - -! -! $URL$ -! $Id$ -! $Revision$ -! $Date$ diff --git a/models/mpas_atm/get_reconstruct_mod.f90 b/models/mpas_atm/MPAS_model_modules/get_reconstruct_mod.f90 similarity index 95% rename from models/mpas_atm/get_reconstruct_mod.f90 rename to models/mpas_atm/MPAS_model_modules/get_reconstruct_mod.f90 index 0227ed76cd..b63110de3c 100644 --- a/models/mpas_atm/get_reconstruct_mod.f90 +++ b/models/mpas_atm/MPAS_model_modules/get_reconstruct_mod.f90 @@ -1,8 +1,3 @@ -! This code may (or may not) be part of the MPAS distribution, -! So it is not protected by the DART copyright agreement. -! -! DART $Id$ - module get_reconstruct_mod use types_mod, only : r8 @@ -10,6 +5,7 @@ module get_reconstruct_mod implicit none + private public :: get_reconstruct_init, get_reconstruct contains @@ -152,9 +148,3 @@ subroutine get_reconstruct(nData, latReconstruct, lonReconstruct, & end subroutine get_reconstruct end module get_reconstruct_mod - -! -! $URL$ -! $Id$ -! $Revision$ -! $Date$ diff --git a/models/mpas_atm/model_mod.f90 b/models/mpas_atm/model_mod.f90 index 02e7361c70..00cf397e57 100644 --- a/models/mpas_atm/model_mod.f90 +++ b/models/mpas_atm/model_mod.f90 @@ -5,12 +5,6 @@ module model_mod ! MPAS Atmosphere model interface to the DART data assimilation system. -! Code in this module is compiled with the DART executables. It isolates -! all information about the MPAS grids, model variables, and other details. -! There are a set of 16 subroutine interfaces that are required by DART; -! these cannot be changed. Additional public routines in this file can -! be used by converters and utilities and those interfaces can be anything -! that is useful to other pieces of code. ! This revision of the model_mod supports both a global MPAS grid and ! a regional grid. For the regional grid only observations which @@ -29,7 +23,8 @@ module model_mod ! by each to make sure the intended actions are what happens. use types_mod, only : r4, r8, i8, digits12, SECPERDAY, MISSING_R8, & - rad2deg, deg2rad, PI, MISSING_I, obstypelength + rad2deg, deg2rad, PI, MISSING_I, obstypelength, & + vtablenamelength use time_manager_mod, only : time_type, set_time, set_date, get_date, get_time,& print_time, print_date, set_calendar_type, & @@ -66,9 +61,9 @@ module model_mod xyz_get_close_type, xyz_get_close_init, xyz_get_close_destroy, & xyz_find_nearest, xyz_use_great_circle_dist -use utilities_mod, only : register_module, error_handler, get_unit, & +use utilities_mod, only : error_handler, get_unit, & E_ERR, E_WARN, E_MSG, E_ALLMSG, logfileunit, & - do_output, to_upper, nmlfileunit, & + to_upper, nmlfileunit, & find_namelist_in_file, check_namelist_read, & open_file, file_exist, find_textfile_dims, & file_to_text, close_file, do_nml_file, & @@ -110,28 +105,19 @@ module model_mod use ensemble_manager_mod, only : ensemble_type, get_my_num_vars, get_my_vars -use distributed_state_mod +use distributed_state_mod, only : get_state, get_state_array -! netcdf modules -use typesizes use netcdf -! RBF (radial basis function) modules, donated by LANL. currently deprecated -! in this version. they did the job but not as well as other techniques and -! at a much greater execution-time code. they were used to interpolate -! values at arbitrary locations, not just at cell centers. with too small -! a set of basis points the values were discontinuous at cell boundaries; -! with too many the values were too smoothed out. we went back to -! barycentric interpolation in triangles formed by the three cell centers -! that enclosed the given point. -use get_geometry_mod -use get_reconstruct_mod - use state_structure_mod, only : add_domain, get_model_variable_indices, & state_structure_info, get_index_start, get_index_end, & get_num_variables, get_domain_size, get_varid_from_varname, & get_variable_name, get_num_dims, get_dim_lengths, & - get_dart_vector_index, get_num_varids_from_kind + get_dart_vector_index, get_num_varids_from_kind, & + get_dim_name, get_varid_from_kind, get_num_domains + +use get_reconstruct_mod, only : get_reconstruct_init, get_reconstruct +use get_geometry_mod, only : get_geometry implicit none @@ -142,7 +128,6 @@ module model_mod ! routines in this list have code in this module public :: get_model_size, & - get_num_vars, & get_state_meta_data, & model_interpolate, & shortest_time_between_assimilations, & @@ -157,7 +142,7 @@ module model_mod read_model_time, & write_model_time -! code for these routines are in other modules +! pass through routines to the default model_mod public :: init_time, & init_conditions, & adv_1step, & @@ -165,44 +150,50 @@ module model_mod ! generally useful routines for various support purposes. ! the interfaces here can be changed as appropriate. - -public :: get_init_template_filename, & - analysis_file_to_statevector, & - statevector_to_analysis_file, & - statevector_to_boundary_file, & - get_analysis_time, & +public :: get_analysis_time, & get_grid_dims, & get_xland, & - get_surftype, & get_cell_center_coords, & get_bdy_mask, & - print_variable_ranges, & find_closest_cell_center, & - find_triangle, & - read_2d_from_nc_file, & - find_height_bounds, & cell_ok_to_interpolate, & - is_global_grid, & - uv_cell_to_edges - -! try adjusting what static_init_model does, before it is called. -! the main update_bc program, for example, can call these before -! calling static_init_model() and that will modify what it does. -public :: set_lbc_variables, & - force_u_into_state + uv_increments_cell_to_edges, & + uv_field_cell_to_edges, & + on_boundary_cell, & + on_boundary_edge, & + get_analysis_weight, & + cell_next_to_boundary_edge -! set_lbc_variables sets the lbc_variables string array -! force_u_into_state sets a logical add_u_to_state_list that forces u to be in state +public :: update_u_from_reconstruct, & + use_increments_for_u_update, & + anl_domid ! HK todo accessor functions? -! version controlled file description for error handling, do not edit character(len=*), parameter :: source = 'models/mpas_atm/model_mod.f90' -character(len=*), parameter :: revision = '' -character(len=*), parameter :: revdate = '' -! module global storage; maintains values between calls, accessible by -! any subroutine +! error codes: +integer, parameter :: GENERAL_ERROR = 1 ! general error +integer, parameter :: CRITICAL_ERROR = 99 ! general error in case something terrible goes wrong +integer, parameter :: VERTICAL_TOO_HIGH = 81 ! Vertical location too high +integer, parameter :: VERTICAL_TOO_LOW = 80 ! Vertical location too low +integer, parameter :: QTY_NOT_IN_STATE_VECTOR = 88 ! qty is not in the state vector +integer, parameter :: TK_COMPUTATION_ERROR = 89 ! tk cannot be computed +integer, parameter :: CELL_CENTER_NOT_FOUND = 11 ! Could not find the closest cell center that contains this lat/lon +integer, parameter :: SURFACE_OBS_TOO_FAR = 12 ! Surface obs too far away from model elevation +integer, parameter :: INTERPOLATION_MISSING_VALUE = 13 ! Missing value in interpolation +integer, parameter :: TRIANGLE_CELL_CENTER_NOT_FOUND = 14 ! Could not find the other two cell centers of the triangle that contains this lat/lon +integer, parameter :: TRIANGLE_CELL_CENTER_IN_BOUNDARY = 15 ! Cell centers of the triangle fall in the lateral boundary zone +integer, parameter :: VERTICAL_VELOCITY_NOT_AVAIL = 16 ! Do not know how to do vertical velocity for now +integer, parameter :: PRESSURE_COMPUTATION_ERROR = 17 ! Unable to compute pressure values +integer, parameter :: ILLEGAL_ALTITUDE = 18 +integer, parameter :: RBF_U_COMPUTATION_ERROR = 19 ! could not compute u using RBF code +integer, parameter :: INTERNAL_ERROR = 101 ! reached end of subroutine without finding an applicable case. +integer, parameter :: REJECT_OBS_USER_PRESSURE_LEVEL = 201 ! Reject observation from user specified pressure level +integer, parameter :: PRESSURE_NOT_MONOTONIC = 988 ! Pressure is not monotonically descreased with level +integer, parameter :: HEIGHT_TOO_LOW = 998 +integer, parameter :: HEIGHT_TOO_HIGH = 9998 + + character(len=512) :: string1, string2, string3 -character(len=256) :: locstring logical, save :: module_initialized = .false. ! length of an mpas (also wrf) time string: YYYY-MM-DD_hh:mm:ss @@ -226,7 +217,6 @@ module model_mod ! set in static_init_model() to 1e-5 or 1e-12 depending on compiled precision real(r8) :: roundoff -! Storage for a random sequence for perturbing a single initial state type(random_seq_type) :: random_seq ! Structure for computing distances to cell centers, and assorted arrays @@ -235,14 +225,10 @@ module model_mod type(xyz_location_type), allocatable :: cell_locs(:) logical :: search_initialized = .false. -! compile-time control over whether grid information is written to the -! diagnostic files or not. if it is, the files are self-contained (e.g. for -! ease of plotting), but are also much larger than they would be otherwise. -! change this private variable to control whether the grid information is -! written or not. -logical :: add_static_data_to_diags = .false. +integer :: anl_domid = -1 -! variables which are in the module namelist +!------------------------------------------------------------------------ +! variables which are in the model_nml namelist character(len=256) :: init_template_filename = 'mpas_init.nc' character(len=256) :: bdy_template_filename = '' integer :: vert_localization_coord = VERTISHEIGHT @@ -253,16 +239,7 @@ module model_mod character(len=32) :: calendar = 'Gregorian' real(r8) :: highest_obs_pressure_mb = 100.0_r8 ! do not assimilate obs higher than this level. real(r8) :: sfc_elev_max_diff = -1.0_r8 ! do not assimilate if |model - station| height is larger than this [m]. -integer :: xyzdebug = 0 -integer :: debug = 0 ! turn up for more and more debug messages -! this is not in the namelist or supported generally. -! (setting this to true avoids the surface elevation max diff -! test for elevation and surface pressure.) -logical :: always_assim_surf_altimeters = .false. - -integer :: anl_domid = -1 ! For state_structure_mod access -integer :: lbc_domid = -1 ! if .false. use U/V reconstructed winds tri interp at centers for wind forward ops ! if .true. use edge normal winds (u) with RBF functs for wind forward ops @@ -282,17 +259,10 @@ module model_mod ! assimilation, and use only the increments to update the edge winds logical :: use_increments_for_u_update = .true. -! if set by: call force_u_into_state() *BEFORE* calling static_init_model(), -! the code will add U (edge normal winds) to the mpas state vector even if it -! isn't in the namelist. -logical :: add_u_to_state_list = .false. - ! if > 0, amount of distance in fractional model levels that a vertical ! point can be above or below the top or bottom of the grid and still be -! evaluated without error. if extrapolate is true, extrapolate from the -! first or last model level. if extrapolate is false, use level 1 or N. +! evaluated without error. real(r8) :: outside_grid_level_tolerance = -1.0_r8 -logical :: extrapolate = .false. ! if the calling code updates an existing file it simply writes the state variable ! data. if it needs to create a file from scratch it calls nc_write_model_atts() @@ -306,25 +276,6 @@ module model_mod ! (backwards compatible with previous code), set this to .true. logical :: no_normalization_of_scale_heights = .true. -! for regional MPAS -real(r8) :: dxmax ! max distance between two adjacent cell centers in the mesh (in meters) - -! when updating boundary files for regional mpas, note whether the boundary -! file has the reconstructed winds (lbc_ur, lbc_vr) or not. (this is set by -! looking at the bdy template file and seeing if those variables are there.) -! if not, the other two options are ignored. -! if they are in the lbc file, then the other logicals control whether to use -! them instead of updating the U edge winds directly, and whether to use the -! reconstructed increments or not. -! the latter two options could be added to the namelist if someone wanted to -! explore the options for how the edge winds are updated in the boundary file. -! for now they're not - they're hardcoded true. -! note that these are for the boundary file update only - there are separate -! options for how to update the U winds in the main assimilation state. - -logical :: lbc_file_has_reconstructed_winds = .false. - - namelist /model_nml/ & init_template_filename, & vert_localization_coord, & @@ -333,69 +284,38 @@ module model_mod model_perturbation_amplitude, & log_p_vert_interp, & calendar, & - debug, & - xyzdebug, & use_u_for_wind, & use_rbf_option, & update_u_from_reconstruct, & use_increments_for_u_update, & highest_obs_pressure_mb, & outside_grid_level_tolerance, & - extrapolate, & sfc_elev_max_diff, & write_grid_to_diag_files, & no_normalization_of_scale_heights -! DART state vector contents are specified in the input.nml:&mpas_vars_nml namelist. -integer, parameter :: max_state_variables = 80 -integer, parameter :: num_state_table_columns = 2 -integer, parameter :: num_bounds_table_columns = 4 -character(len=NF90_MAX_NAME) :: mpas_state_variables(max_state_variables * num_state_table_columns ) = ' ' -character(len=NF90_MAX_NAME) :: mpas_state_bounds(num_bounds_table_columns, max_state_variables ) = ' ' -character(len=NF90_MAX_NAME) :: variable_table(max_state_variables, num_state_table_columns ) +!------------------------------------------------------------------------ -! this is special and not in the namelist. boundary files have a fixed -! set of variables with fixed names. -character(len=NF90_MAX_NAME) :: lbc_variables(max_state_variables) = '' +integer, parameter :: MAX_STATE_VARIABLES = 80 +integer, parameter :: NUM_STATE_TABLE_COLUMNS = 2 +integer, parameter :: NUM_BOUNDS_TABLE_COLUMNS = 4 !HK @todo get rid of clamp or fail +character(len=vtablenamelength) :: variable_table(MAX_STATE_VARIABLES, NUM_STATE_TABLE_COLUMNS ) + +!------------------------------------------------------------------------ +! mpas_vars_nml namelist +character(len=vtablenamelength) :: mpas_state_variables(MAX_STATE_VARIABLES * NUM_STATE_TABLE_COLUMNS ) = ' ' +character(len=vtablenamelength) :: mpas_state_bounds(NUM_BOUNDS_TABLE_COLUMNS, MAX_STATE_VARIABLES ) = ' ' namelist /mpas_vars_nml/ mpas_state_variables, mpas_state_bounds +!------------------------------------------------------------------------ -integer :: nfields +!------------------------------------------------------------------------ +! for regional MPAS +real(r8) :: dxmax ! max distance between two adjacent cell centers in the mesh (in meters) +!------------------------------------------------------------------------ -!>@todo FIXME - REMOVE AS MUCH OF THIS AS POSSIBLE. -!> some of this information is in the state structure now. -!> the duplicate progvar stuff should be removed and the -!> state routines used instead. this duplicates work and -!> makes us keep up code in 2 different places. - -! original code: -! Everything needed to describe a variable - -type progvartype - private - character(len=NF90_MAX_NAME) :: varname - character(len=NF90_MAX_NAME), dimension(NF90_MAX_VAR_DIMS) :: dimname - integer, dimension(NF90_MAX_VAR_DIMS) :: dimlens - integer :: xtype ! netCDF variable type (NF90_double, etc.) - integer :: numdims ! number of dims - excluding TIME - integer :: numvertical ! number of vertical levels in variable - integer :: numcells ! number of horizontal locations (cell centers) - integer :: numedges ! number of horizontal locations (edges for velocity components) - logical :: ZonHalf ! vertical coordinate has dimension nVertLevels - 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 - character(len=obstypelength) :: kind_string - logical :: clamping ! does variable need to be range-restricted before - real(r8) :: range(2) ! being stuffed back into MPAS analysis file. - logical :: out_of_range_fail ! is out of range fatal if range-checking? -end type progvartype - -type(progvartype), dimension(max_state_variables) :: progvar ! Grid parameters - the values will be read from an mpas analysis file. - integer :: nCells = -1 ! Total number of cells making up the grid integer :: nVertices = -1 ! Unique points in grid that are corners of cells integer :: nEdges = -1 ! Straight lines between vertices making up cells @@ -406,19 +326,6 @@ module model_mod integer :: nSoilLevels = -1 ! Number of soil layers ! scalar grid positions - -! FIXME: we read in a lot of metadata about the grids. if space becomes an -! issue we could consider reading in only the x,y,z arrays for all the items -! plus the radius, and then compute the lat/lon for locations needed by -! get_state_meta_data() on demand. most of the computations we need to do -! are actually easier in xyz coords (no issue with poles). - -! FIXME: it may be desirable to read in xCell(:), yCell(:), zCell(:) -! to keep from having to compute them on demand, especially since we -! have converted the radian lat/lon of the cell centers into degrees. -! we have to convert back, then take a few sin and cos to get xyz. -! time/space/accuracy tradeoff here. - real(r8), allocatable :: xVertex(:), yVertex(:), zVertex(:) real(r8), allocatable :: xEdge(:), yEdge(:), zEdge(:) real(r8), allocatable :: lonEdge(:) ! edge longitudes (degrees, original radians in file) @@ -432,8 +339,6 @@ module model_mod real(r8), allocatable :: zGridFace(:,:) ! geometric height at cell faces (nVertLevelsP1,nCells) real(r8), allocatable :: zGridCenter(:,:) ! geometric height at cell centers (nVertLevels, nCells) real(r8), allocatable :: zGridEdge(:,:) ! geometric height at edge centers (nVertLevels, nEdges) -!real(r8), allocatable :: zEdgeFace(:,:) ! geometric height at edges faces (nVertLevelsP1,nEdges) -!real(r8), allocatable :: zEdgeCenter(:,:) ! geometric height at edges faces (nVertLevels ,nEdges) integer, allocatable :: cellsOnVertex(:,:) ! list of cell centers defining a triangle integer, allocatable :: verticesOnCell(:,:) @@ -445,14 +350,13 @@ module model_mod ! Boundary information might be needed ... regional configuration? ! Read if available. - integer, allocatable :: bdyMaskCell(:) integer, allocatable :: bdyMaskEdge(:) integer, allocatable :: maxLevelCell(:) integer, allocatable :: surftype(:) ! ! surface type (land=0, water=1, seaice = 2) - for rttov -integer :: model_size ! the state vector length +integer(i8) :: model_size type(time_type) :: model_timestep ! smallest time to adv model ! useful flags in making decisions when searching for points, etc @@ -462,53 +366,13 @@ module model_mod logical :: has_uvreconstruct = .false. ! true = has reconstructed at centers ! Do we have any state vector items located on the cell edges? -! If not, avoid reading in or using the edge arrays to save space. -! FIXME: this should be set after looking at the fields listed in the -! namelist which are to be read into the state vector - if any of them -! are located on the edges then this flag should be changed to .true. -! however, the way the code is structured these arrays are allocated -! before the details of field list is examined. since right now the -! only possible field array that is on the edges is the 'u' edge normal -! winds, search specifically for that in the state field list and set -! this based on that. if any other data might be on edges, this routine -! will need to be updated: is_edgedata_in_state_vector() logical :: data_on_edges = .false. -! currently unused; for a regional model it is going to be necessary to know -! if the grid is continuous around longitudes (wraps in east-west) or not, -! and if it covers either of the poles. -character(len= 64) :: ew_boundary_type, ns_boundary_type - -! common names that call specific subroutines based on the arg types -INTERFACE vector_to_prog_var - MODULE PROCEDURE vector_to_1d_prog_var - MODULE PROCEDURE vector_to_2d_prog_var - MODULE PROCEDURE vector_to_3d_prog_var -END INTERFACE - -INTERFACE prog_var_to_vector - MODULE PROCEDURE prog_var_1d_to_vector - MODULE PROCEDURE prog_var_2d_to_vector - MODULE PROCEDURE prog_var_3d_to_vector -END INTERFACE - INTERFACE get_analysis_time MODULE PROCEDURE get_analysis_time_ncid MODULE PROCEDURE get_analysis_time_fname END INTERFACE -INTERFACE get_index_range - MODULE PROCEDURE get_index_range_int - MODULE PROCEDURE get_index_range_string -END INTERFACE - - -interface write_model_time - module procedure write_model_time_file - module procedure write_model_time_restart -end interface - - !------------------------------------------------ ! The regular grid used for triangle interpolation divides the sphere into @@ -526,61 +390,29 @@ module model_mod ! ??? is sufficient for ??? integer, parameter :: max_reg_list_num = 100 -! The triangle interpolation keeps a list of how many and which triangles -! overlap each regular lon-lat box. The number is stored in -! array triangle_num. The allocatable array -! triangle_list lists the uniquen index -! of each overlapping triangle. The entry in -! triangle_start for a given regular lon-lat box indicates -! where the list of triangles begins in the triangle_list. - -integer :: triangle_start(num_reg_x, num_reg_y) -integer :: triangle_num (num_reg_x, num_reg_y) = 0 -integer, allocatable :: triangle_list(:) - - contains -!================================================================== -! All the public REQUIRED interfaces come first - just by convention. -!================================================================== - - !------------------------------------------------------------------ - subroutine static_init_model() -!>@todo FIXME - can we add an optional arg here for update_bc use? - ! Called to do one time initialization of the model. -! -! All the grid information comes from the initialization of -! the dart_model_mod module. - -! Local variables - all the important ones have module scope - integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs character(len=NF90_MAX_NAME) :: varname,dimname character(len=obstypelength) :: kind_string integer :: ncid, VarID, numdims, varsize, dimlen -integer :: iunit, io, ivar, i, index1, indexN, iloc, kloc +integer :: iunit, io, ivar, i integer :: ss, dd, z1, m1 integer :: nDimensions, nVariables, nAttributes, unlimitedDimID, TimeDimID -integer :: cel1, cel2, lbc_nfields +integer :: nfields logical :: both -real(r8) :: variable_bounds(max_state_variables, 2) -integer :: variable_qtys(max_state_variables) +real(r8) :: variable_bounds(MAX_STATE_VARIABLES, 2) +integer :: variable_qtys(MAX_STATE_VARIABLES) character(len=*), parameter :: routine = 'static_init_model' -if ( module_initialized ) return ! only need to do this once. +if ( module_initialized ) return -! Since this routine calls other routines that could call this routine -! we'll say we've been initialized pretty dang early. module_initialized = .true. -! Print module information to log file and stdout. -call register_module(source, revision, revdate) - ! Read the DART namelist for this model call find_namelist_in_file('input.nml', 'model_nml', iunit) read(iunit, nml = model_nml, iostat = io) @@ -591,41 +423,120 @@ subroutine static_init_model() if (do_nml_term()) write( * , nml=model_nml) ! Read the MPAS variable list to populate DART state vector -! Intentionally do not try to dump them to the nml unit because -! they include large character arrays which output pages of white space. -! The routine that reads and parses this namelist will output what -! values it found into the log. call find_namelist_in_file('input.nml', 'mpas_vars_nml', iunit) read(iunit, nml = mpas_vars_nml, iostat = io) call check_namelist_read(iunit, io, 'mpas_vars_nml') -!--------------------------------------------------------------- -! Set the time step ... causes mpas namelists to be read. -! Ensures model_timestep is multiple of 'dynamics_timestep' +call set_calendar_type( calendar ) +model_timestep = set_model_time_step() +call get_time(model_timestep,ss,dd) +write(string1,*)'assimilation window is ',dd,' days ',ss,' seconds' +call error_handler(E_MSG,routine,string1,source) + -call set_calendar_type( calendar ) ! comes from model_mod_nml +ncid = nc_open_file_readonly(init_template_filename, routine) -model_timestep = set_model_time_step() +call verify_state_variables(ncid, init_template_filename, & + nfields, variable_qtys, variable_bounds) + +anl_domid = add_domain(init_template_filename, nfields, & + var_names = variable_table (1:nfields,1), & + kind_list = variable_qtys(1:nfields), & + clamp_vals = variable_bounds(1:nfields,:) ) -call get_time(model_timestep,ss,dd) ! set_time() assures the seconds [0,86400) +call read_grid() -write(string1,*)'assimilation window is ',dd,' days ',ss,' seconds' -call error_handler(E_MSG,routine,string1,source,revision,revdate) +call nc_close_file(ncid, routine) -!--------------------------------------------------------------- -! 1) get grid dimensions -! 2) allocate space for the grids -! 3) read them from the analysis file +model_size = get_domain_size(anl_domid) -ncid = nc_open_file_readonly(init_template_filename, routine) +! if you have at least one of these wind components in the state vector, +! you have to have them both. the subroutine will error out if only one +! is found and not both. +if (has_uvreconstruct) then + call winds_present(z1,m1,both) +endif + +! if you set the namelist to use the reconstructed cell center winds, +! they have to be in the state vector. +if (update_u_from_reconstruct .and. .not. has_uvreconstruct) then + write(string1,*) 'update_u_from_reconstruct cannot be True' + write(string2,*) 'because state vector does not contain U/V reconstructed winds' + call error_handler(E_ERR,'static_init_model',string1,source, & + text2=string2) +endif + +! if you set the namelist to update the edge normal winds based on the +! updated cell center winds, and you also have the edge normal winds in +! the state vector, warn that the edge normal winds will be ignored +! when going back to the mpas_analysis.nc file. not an error, but the +! updates to the edge normal vectors won't affect the results. +if (update_u_from_reconstruct .and. has_edge_u) then + write(string1,*) 'edge normal winds (u) in MPAS file will be updated based on U/V reconstructed winds' + write(string2,*) 'and not from the updated edge normal values in the state vector' + write(string3,*) 'because update_u_from_reconstruct is True' + call error_handler(E_MSG,'static_init_model',string1,source, & + text2=string2, text3=string3) +endif + +! there are two choices when updating the edge normal winds based on the +! winds at the cell centers. one is a direct interpolation of the values; +! the other is to read in the original cell centers, compute the increments +! changed by the interpolation, and then add or substract only the increments +! from the original edge normal wind values. +if (update_u_from_reconstruct) then + if (use_increments_for_u_update) then + write(string1,*) 'edge normal winds (u) in MPAS file will be updated based on the difference' + write(string2,*) 'between the original U/V reconstructed winds and the updated values' + write(string3,*) 'because use_increments_for_u_update is True' + else + write(string1,*) 'edge normal winds (u) in MPAS file will be updated by averaging the' + write(string2,*) 'updated U/V reconstructed winds at the corresponding cell centers' + write(string3,*) 'because use_increments_for_u_update is False' + endif + call error_handler(E_MSG,'static_init_model',string1,source, & + text2=string2, text3=string3) +endif + +! Required state vector fields +if ((get_num_varids_from_kind(anl_domid, QTY_POTENTIAL_TEMPERATURE) < 0) .or. & + (get_num_varids_from_kind(anl_domid, QTY_DENSITY) < 0) .or. & + (get_num_varids_from_kind(anl_domid, QTY_VAPOR_MIXING_RATIO) < 0)) then + write(string1, *) 'State vector is missing one or more of the following fields:' + write(string2, *) 'Potential Temperature (theta), Density (rho), Vapor Mixing Ratio (qv).' + write(string3, *) 'Cannot convert between height/pressure nor compute sensible temps.' + call error_handler(E_ERR,'static_init_model',string1,source, & + text2=string2, text3=string3) +endif + +call set_vertical_localization_coord(vert_localization_coord) + +if (r8 == digits12) then ! double precision + roundoff = 1.0e-12_r8 +else ! single precision + roundoff = 1.0e-5_r8 +endif + +end subroutine static_init_model + +!------------------------------------------------------------------ +subroutine read_grid() + +integer :: cel1, cel2 +integer :: iloc, kloc +integer :: ncid +character(len=*), parameter :: routine = 'read_grid' -! move this up - some of the routines below depend on having -! the variable table parsed already, and U added if needed. -call verify_state_variables( mpas_state_variables, ncid, init_template_filename, & - nfields, variable_table) +ncid = nc_open_file_readonly(init_template_filename, routine) -! get sizes -call read_grid_dims(ncid) +nCells = nc_get_dimension_size(ncid, 'nCells', routine) +nVertices = nc_get_dimension_size(ncid, 'nVertices', routine) +nEdges = nc_get_dimension_size(ncid, 'nEdges', routine) +maxEdges = nc_get_dimension_size(ncid, 'maxEdges', routine) +nVertLevels = nc_get_dimension_size(ncid, 'nVertLevels', routine) +nVertLevelsP1 = nc_get_dimension_size(ncid, 'nVertLevelsP1', routine) +vertexDegree = nc_get_dimension_size(ncid, 'vertexDegree', routine) +nSoilLevels = nc_get_dimension_size(ncid, 'nSoilLevels', routine) allocate(latCell(nCells), lonCell(nCells)) allocate(zGridFace(nVertLevelsP1, nCells)) @@ -644,9 +555,9 @@ subroutine static_init_model() allocate(edgeNormalVectors(3, nEdges)) allocate(xVertex(nVertices), yVertex(nVertices), zVertex(nVertices)) -! see if U is in the state vector list. if not, don't read in or +! see if any variables on Edges are in the state vector. If not, do not read in or ! use any of the Edge arrays to save space. -data_on_edges = is_edgedata_in_state_vector(variable_table, lbc_variables) +data_on_edges = is_edgedata_in_state_vector() if(data_on_edges) then allocate(zGridEdge(nVertLevels, nEdges)) @@ -684,813 +595,324 @@ subroutine static_init_model() zGridEdge(iloc,kloc) = zGridCenter(iloc,cel2) else !this is bad... write(string1,*)'Edge ',kloc,' at vertlevel ',iloc,' has no neighbouring cells!' - call error_handler(E_ERR,'static_init_model', string1, source, revision, revdate) + call error_handler(E_ERR,'static_init_model', string1, source) endif enddo enddo endif -!--------------------------------------------------------------- -! Compile the list of model variables to use in the creation -! of the DART state vector. -! -! THIS CODE SHOULD BE REMOVED - it is done by the add_domain code. -! - - -TimeDimID = FindTimeDimension( ncid ) - -if (TimeDimID < 0 ) then - write(string1,*)'unable to find a dimension named Time.' - call error_handler(E_MSG,'static_init_model', string1, source, revision, revdate) -endif - -call nc_check(nf90_Inquire(ncid,nDimensions,nVariables,nAttributes,unlimitedDimID), & - 'static_init_model', 'inquire '//trim(init_template_filename)) - -if ( (TimeDimID > 0) .and. (unlimitedDimID > 0) .and. (TimeDimID /= unlimitedDimID)) then - write(string1,*)'IF Time is not the unlimited dimension, I am lost.' - call error_handler(E_MSG,'static_init_model', string1, source, revision, revdate) -endif - -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_index_for_quantity( progvar(ivar)%kind_string ) - progvar(ivar)%numdims = 0 - progvar(ivar)%numvertical = 1 - progvar(ivar)%dimlens = MISSING_I - progvar(ivar)%numcells = MISSING_I - progvar(ivar)%numedges = MISSING_I +call nc_close_file(ncid, routine) - string2 = trim(init_template_filename)//' '//trim(varname) - call nc_check(nf90_inq_varid(ncid, trim(varname), VarID), & - 'static_init_model', 'inq_varid '//trim(string2)) +end subroutine read_grid - call nc_check(nf90_inquire_variable(ncid, VarID, xtype=progvar(ivar)%xtype, & - dimids=dimIDs, ndims=numdims), 'static_init_model', 'inquire '//trim(string2)) - ! Since we are not concerned with the TIME dimension, we need to skip it. - ! When the variables are read, only a single timestep is ingested into - ! the DART state vector. +!------------------------------------------------------------------ - varsize = 1 - dimlen = 1 - DimensionLoop : do i = 1,numdims +subroutine get_state_meta_data(index_in, location, qty) - if (dimIDs(i) == TimeDimID) cycle DimensionLoop +! Given an index into the state vector, return its location and +! optionally the qty - write(string1,'(''inquire dimension'',i2,A)') i,trim(string2) - call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), len=dimlen, name=dimname), & - 'static_init_model', string1) +integer(i8), intent(in) :: index_in +type(location_type), intent(out) :: location +integer, optional, intent(out) :: qty - progvar(ivar)%numdims = progvar(ivar)%numdims + 1 - progvar(ivar)%dimlens(i) = dimlen - progvar(ivar)%dimname(i) = trim(dimname) - varsize = varsize * dimlen +integer :: index, i, j, k, dom_id, var_id, qty_local, vert, cell, level +real(r8) :: lon, lat, height +type(location_type) :: new_location - select case ( dimname(1:6) ) - case ('nCells') - progvar(ivar)%numcells = dimlen - case ('nEdges') - progvar(ivar)%numedges = dimlen - case ('nVertL') ! nVertLevels, nVertLevelsP1, nVertLevelsP2 - progvar(ivar)%numvertical = dimlen - case ('nSoilL') ! nSoilLevels - progvar(ivar)%numvertical = dimlen - end select +if ( .not. module_initialized ) call static_init_model - enddo DimensionLoop +call get_model_variable_indices(index_in, i, j, k, dom_id=dom_id, var_id=var_id, kind_index=qty_local) - ! this call sets: clamping, bounds, and out_of_range_fail in the progvar entry - call get_variable_bounds(mpas_state_bounds, ivar) +! Need to find if the variable is 2d (1d netcdf) or 3d (2d netcdf) +! Need to find if variable is on cells or edges +! Need to find if variable is on nVertLevels or nVertLevelsP1 - if (progvar(ivar)%numvertical == nVertLevels) then - progvar(ivar)%ZonHalf = .TRUE. +if (get_num_dims(dom_id, var_id) == 1) then ! 2d (1d netcdf) variable: (cells) + cell = i + level = 1 + vert = VERTISSURFACE + if (on_edge(dom_id, var_id)) then + lon = lonEdge(cell) + lat = latEdge(cell) + height = zGridEdge(level, cell) else - progvar(ivar)%ZonHalf = .FALSE. + lon = lonCell(cell) + lat = latCell(cell) + height = zGridFace(level, cell) !HK 2D on face in original code endif - if (varname == 'u') has_edge_u = .true. - if (varname == 'uReconstructZonal' .or. & - varname == 'uReconstructMeridional') has_uvreconstruct = .true. +else ! 2d (1d netcdf) variable: (levels, cells) + cell = j + level = i + vert = VERTISHEIGHT + if (on_edge(dom_id, var_id)) then + lon = lonEdge(cell) + lat = latEdge(cell) + height = zGridEdge(level, cell) + else + lon = lonCell(cell) + lat = latCell(cell) + if (get_dim_name(dom_id, var_id, 1) == 'nVertLevels') then + height = zGridCenter(level, cell) + else + height = zGridFace(level, cell) + endif + endif +endif - progvar(ivar)%varsize = varsize - progvar(ivar)%index1 = index1 - progvar(ivar)%indexN = index1 + varsize - 1 - index1 = index1 + varsize ! sets up for next variable +location = set_location(lon, lat, height, vert) - if ( debug > 11 .and. do_output()) call dump_progvar(ivar) +if (present(qty)) then + qty = qty_local +endif -enddo +end subroutine get_state_meta_data -call nc_close_file(ncid, routine) +!------------------------------------------------------------------ +function on_edge(dom, var) +integer, intent(in) :: dom, var +logical :: on_edge -! FIXME: moved below to after add_domain() calls -!model_size = progvar(nfields)%indexN - -if ( debug > 0 .and. do_output()) then - write(logfileunit,*) - write( * ,*) - write(logfileunit,'(" static_init_model: nCells, nEdges, nVertices, nVertLevels =",4(1x,i9))') & - nCells, nEdges, nVertices, nVertLevels - write( * ,'(" static_init_model: nCells, nEdges, nVertices, nVertLevels =",4(1x,i9))') & - nCells, nEdges, nVertices, nVertLevels -! write(logfileunit, *)'static_init_model: model_size = ', model_size -! write( * , *)'static_init_model: model_size = ', model_size - if ( global_grid ) then - write(logfileunit, *)'static_init_model: grid is a global grid ' - write( * , *)'static_init_model: grid is a global grid ' - else - write(logfileunit, *)'static_init_model: grid is NOT a global grid. Lateral boundaries exist ' - write( * , *)'static_init_model: grid is NOT a global grid. Lateral boundaries exist ' - endif - if ( all_levels_exist_everywhere ) then - write(logfileunit, *)'static_init_model: all cells have same number of vertical levels ' - write( * , *)'static_init_model: all cells have same number of vertical levels ' - else - write(logfileunit, *)'static_init_model: cells have varying number of vertical levels ' - write( * , *)'static_init_model: cells have varying number of vertical levels ' - endif +if (get_num_dims(dom, var) == 1) then + on_edge = ( get_dim_name(dom, var, 1) == 'nEdges' ) +else + on_edge = ( get_dim_name(dom, var, 2) == 'nEdges' ) endif +end function on_edge -! set the domain(s) in the state structure here -variable_bounds(1:nfields, 1) = progvar(1:nfields)%range(1) -variable_bounds(1:nfields, 2) = progvar(1:nfields)%range(2) -variable_qtys(1:nfields) = progvar(1:nfields)%dart_kind +!------------------------------------------------------------------ +!> given a state ensemble_handle, a location, and a QTY_xxx, return the +!> interpolated value at that location, and an error code. 0 is success, +!> anything positive is an error. (negative reserved for system use) +subroutine model_interpolate(state_handle, ens_size, location, qty, expected_obs, istatus) -anl_domid = add_domain(init_template_filename, nfields, & - var_names = variable_table (1:nfields,1), & - kind_list = variable_qtys(1:nfields), & - clamp_vals = variable_bounds(1:nfields,:) ) +type(ensemble_type), intent(in) :: state_handle +integer, intent(in) :: ens_size +type(location_type), intent(in) :: location +integer, intent(in) :: qty +real(r8), intent(out) :: expected_obs(ens_size) +integer, intent(out) :: istatus(ens_size) -model_size = get_domain_size(anl_domid) -if ( debug > 4 .and. do_output()) print*,'model_size(anl_domid)=',model_size ! HA - -lbc_nfields = 0 - -! if we have a lateral boundary file, add it to the domain -! so we have access to the corresponding lbc_xxx fields. -!>@todo FIXME: if we want to do increments, we could also add a -! third domain which is the original forecast fields before -! the assimilation (so we can compute increments) -if (.not. global_grid .and. lbc_variables(1) /= '') then - ! regional: count number of lbc fields to read in - COUNTUP: do i=1, max_state_variables - if (lbc_variables(i) /= '') then - lbc_nfields = lbc_nfields + 1 - else - exit COUNTUP - endif - enddo COUNTUP - if( debug > 4 .and. do_output()) print*, 'Regional: number of lbc fields to read in = ', lbc_nfields - lbc_domid = add_domain(bdy_template_filename, lbc_nfields, & - var_names = lbc_variables) - ! FIXME clamp_vals = variable_bounds(1:nfields,:) ) - if( debug > 4 .and. do_output()) print*, 'model_size, lbc_domid =',model_size,lbc_domid - - model_size = model_size + get_domain_size(lbc_domid) -else - lbc_domid = -1 -endif -if ( debug > 4 .and. do_output()) then - call state_structure_info(anl_domid) - print*, 'model_size =',model_size - if(lbc_domid >= 0) print*, 'get_domain_size(lbc_domid):',get_domain_size(lbc_domid) -endif -!if ( debug > 4 .and. do_output() .and. lbc_domid >= 0) call state_structure_info(lbc_domid) +type(location_type) :: location_tmp(ens_size) +integer :: ivar, obs_kind +integer :: tvars(3) +integer :: cellid +logical :: goodkind, surface_obs +real(r8) :: lpres(ens_size), values(3, ens_size) +real(r8) :: llv(3) ! lon/lat/vert +integer :: e, verttype +real(r8) :: single_expected_obs +integer :: single_istatus -! do some sanity checking here: +if ( .not. module_initialized ) call static_init_model -! if you have at least one of these wind components in the state vector, -! you have to have them both. the subroutine will error out if only one -! is found and not both. -if (has_uvreconstruct) then - call winds_present(z1,m1,both) -endif +expected_obs = MISSING_R8 +istatus = INTERNAL_ERROR -! if you set the namelist to use the reconstructed cell center winds, -! they have to be in the state vector. -if (update_u_from_reconstruct .and. .not. has_uvreconstruct) then - write(string1,*) 'update_u_from_reconstruct cannot be True' - write(string2,*) 'because state vector does not contain U/V reconstructed winds' - call error_handler(E_ERR,'static_init_model',string1,source,revision,revdate, & - text2=string2) -endif +llv = get_location(location) +verttype = nint(query_location(location)) +surface_obs = (verttype == VERTISSURFACE) -! if you set the namelist to update the edge normal winds based on the -! updated cell center winds, and you also have the edge normal winds in -! the state vector, warn that the edge normal winds will be ignored -! when going back to the mpas_analysis.nc file. not an error, but the -! updates to the edge normal vectors won't affect the results. -if (update_u_from_reconstruct .and. has_edge_u) then - write(string1,*) 'edge normal winds (u) in MPAS file will be updated based on U/V reconstructed winds' - write(string2,*) 'and not from the updated edge normal values in the state vector' - write(string3,*) 'because update_u_from_reconstruct is True' - call error_handler(E_MSG,'static_init_model',string1,source,revision,revdate, & - text2=string2, text3=string3) +cellid = cell_ok_to_interpolate(location) +if (cellid < 1) then + istatus = CELL_CENTER_NOT_FOUND + return endif -! there are two choices when updating the edge normal winds based on the -! winds at the cell centers. one is a direct interpolation of the values; -! the other is to read in the original cell centers, compute the increments -! changed by the interpolation, and then add or substract only the increments -! from the original edge normal wind values. -if (update_u_from_reconstruct) then - if (use_increments_for_u_update) then - write(string1,*) 'edge normal winds (u) in MPAS file will be updated based on the difference' - write(string2,*) 'between the original U/V reconstructed winds and the updated values' - write(string3,*) 'because use_increment_for_u_update is True' - else - write(string1,*) 'edge normal winds (u) in MPAS file will be updated by averaging the' - write(string2,*) 'updated U/V reconstructed winds at the corresponding cell centers' - write(string3,*) 'because use_increment_for_u_update is False' +! Reject obs if the station height is far way from the model terrain. (but not for rttov!) +if (.not. required_for_rttov(qty)) then + if(surface_obs .and. sfc_elev_max_diff >= 0) then + if(abs(llv(3) - zGridFace(1,cellid)) > sfc_elev_max_diff) then + istatus = SURFACE_OBS_TOO_FAR + return + endif endif - call error_handler(E_MSG,'static_init_model',string1,source,revision,revdate, & - text2=string2, text3=string3) endif -! basically we cannot do much without having at least these -! three fields in the state vector. refuse to go further -! if these are not present: -!print *, get_num_varids_from_kind(anl_domid, QTY_POTENTIAL_TEMPERATURE) -!print *, get_num_varids_from_kind(anl_domid, QTY_DENSITY) -!print *, get_num_varids_from_kind(anl_domid, QTY_VAPOR_MIXING_RATIO) -if ((get_num_varids_from_kind(anl_domid, QTY_POTENTIAL_TEMPERATURE) < 0) .or. & - (get_num_varids_from_kind(anl_domid, QTY_DENSITY) < 0) .or. & - (get_num_varids_from_kind(anl_domid, QTY_VAPOR_MIXING_RATIO) < 0)) then - write(string1, *) 'State vector is missing one or more of the following fields:' - write(string2, *) 'Potential Temperature (theta), Density (rho), Vapor Mixing Ratio (qv).' - write(string3, *) 'Cannot convert between height/pressure nor compute sensible temps.' - call error_handler(E_ERR,'static_init_model',string1,source,revision,revdate, & - text2=string2, text3=string3) +if (.not. qty_ok_to_interpolate(qty)) then + istatus = QTY_NOT_IN_STATE_VECTOR + return endif - -if (extrapolate) then - call error_handler(E_MSG,'static_init_model',& - 'extrapolate not supported yet; will use level 1 or N values', & - source, revision, revdate) +! Reject obs above a user specified pressure level +if (highest_obs_pressure_mb > 0.0) then + if (.not. surface_obs) then + if (.not.(is_vertical(location, "UNDEFINED"))) then + call compute_pressure_at_loc(state_handle, ens_size, location, lpres, istatus) + if (fails(istatus)) then + istatus = PRESSURE_COMPUTATION_ERROR + return + endif + do e = 1, ens_size + if (lpres(e) < highest_obs_pressure_mb * 100.0_r8 ) then + istatus(e) = REJECT_OBS_USER_PRESSURE_LEVEL + endif + enddo + if (fails(istatus)) then + istatus = REJECT_OBS_USER_PRESSURE_LEVEL + return + endif + endif + endif endif -! tell the location module how we want to localize in the vertical -call set_vertical_localization_coord(vert_localization_coord) +! HK @todo 2M temp, 10M winds. #768 -! set an appropriate value for roundoff tests based -! on this code being compiled single or double precision. -! set to 1e-5 (for single) or 1e-12 (for double precision). -if (r8 == digits12) then - roundoff = 1.0e-12_r8 -else - roundoff = 1.0e-5_r8 -endif - -end subroutine static_init_model - - -!------------------------------------------------------------------ - -subroutine get_state_meta_data(index_in, location, var_type) - -! given an index into the state vector, return its location and -! if given, the var kind. despite the name, var_type is a generic -! kind, like those in obs_kind/obs_kind_mod.f90, starting with QTY_ - -! passed variables -integer(i8), intent(in) :: index_in -type(location_type), intent(out) :: location -integer, optional, intent(out) :: var_type - -! Local variables - -integer :: nzp, iloc, vloc, nf, ndim -real(r8) :: height -type(location_type) :: new_location - -if ( .not. module_initialized ) call static_init_model - -! get the local indicies and type from dart index. kloc is a dummy variable for this subroutine - -call find_mpas_indices(index_in, iloc, vloc, ndim, nf) - -nzp = progvar(nf)%numvertical - -! the zGrid array contains the location of the cell top and bottom faces, so it has one -! more value than the number of cells in each column. for locations of cell centers -! you have to take the midpoint of the top and bottom face of the cell. -if (progvar(nf)%numedges /= MISSING_I) then - if (.not. data_on_edges) then - call error_handler(E_ERR, 'get_state_meta_data', & - 'Internal error: numedges present but data_on_edges false', & - source, revision, revdate, text2='variable '//trim(progvar(nf)%varname)) - endif - if ( progvar(nf)%ZonHalf ) then - height = zGridEdge(vloc,iloc) - else - call error_handler(E_ERR, 'get_state_meta_data', 'no support for edges at face heights', & - source, revision, revdate) - endif - - if (nzp <= 1) then - location = set_location(lonEdge(iloc),latEdge(iloc), height, VERTISSURFACE) - else - location = set_location(lonEdge(iloc),latEdge(iloc), height, VERTISHEIGHT) - endif -else - if ( progvar(nf)%ZonHalf ) then - height = zGridCenter(vloc,iloc) - else if (nzp <= 1) then - height = zGridFace(1,iloc) - else - height = zGridFace(vloc,iloc) - endif - - if (nzp <= 1) then - location = set_location(lonCell(iloc),latCell(iloc), height, VERTISSURFACE) - else - location = set_location(lonCell(iloc),latCell(iloc), height, VERTISHEIGHT) - endif -endif - -! we are no longer passing in an ensemble handle to this routine, so we -! cannot do vertical conversion here. assim_tools will call vertical conversion -! on the obs and on the state. - -if (debug > 20 .and. do_output()) then - - write(*,'("INDEX_IN / IVAR : ",(i10,2x),(i5,2x))') index_in, nf - write(*,'(" ILOC, VLOC: ",2(i5,2x))') iloc, vloc - write(*,'(" LON/LAT/HGT: ",3(f12.3,2x))') lonCell(iloc), latCell(iloc), height - -endif - -if (present(var_type)) then - var_type = progvar(nf)%dart_kind -endif - -end subroutine get_state_meta_data - -!------------------------------------------------------------------ -!> given an index into the state vector, return with the cellid -!> and the vertical level if this is a 2d variable. also return -!> the dimensionality, and optionally the progvar index. - -subroutine find_mpas_indices(index_in, cellid, vert_level, ndim, nf) -integer(i8), intent(in) :: index_in -integer, intent(out) :: cellid -integer, intent(out) :: vert_level -integer, intent(out), optional :: ndim -integer, intent(out), optional :: nf - -integer :: i, j, k ! Indices into variable (note k is not used in MPAS) -integer :: nzp, iloc, vloc, nnf - -! get the local indicies and type from dart index. 'k' is a dummy variable for this subroutine - -call get_model_variable_indices(index_in, i, j, k, var_id=nnf) - -if (progvar(nnf)%numdims == 2) then ! variable(vcol, iloc) - vert_level = i - cellid = j -elseif (progvar(nnf)%numdims == 1) then ! variable(iloc) - cellid = i - vert_level = 1 -else - call error_handler(E_ERR, 'find_mpas_indices ', 'expecting 1D or 2D variable') -endif +select case (qty) -if (present(ndim)) ndim = progvar(nnf)%numdims -if (present(nf)) nf = nnf - -end subroutine find_mpas_indices - -!------------------------------------------------------------------ -!> given a domain and varid, return 3 dimensions - setting to 1 if -!> not present. - -subroutine find_mpas_dims(domid, ivar, ndims, dims) -integer, intent(in) :: domid -integer, intent(in) :: ivar -integer, intent(out) :: ndims -integer, intent(out) :: dims(3) - -! start with everything 1, then set values which are different -dims(:) = 1 - -! get the domain info -ndims = get_num_dims(domid, ivar) -dims(1:ndims) = get_dim_lengths(domid, ivar) - -end subroutine find_mpas_dims - -!------------------------------------------------------------------ -!> given a state vector, a location, and a QTY_xxx, return the -!> interpolated value at that location, and an error code. 0 is success, -!> anything positive is an error. (negative reserved for system use) -!> -!> ERROR codes: -!> -!> ISTATUS = 1: general error for rttov - at least one of the input variables goes wrong -!> ISTATUS = 99: general error in case something terrible goes wrong... -!> ISTATUS = 81: Vertical location too high -!> ISTATUS = 80: Vertical location too low -!> ISTATUS = 88: this kind is not in the state vector -!> ISTATUS = 89: tk cannot be computed. -!> ISTATUS = 11: Could not find the closest cell center that contains this lat/lon -!> ISTATUS = 12: Surface obs too far away from model elevation -!> ISTATUS = 13: Missing value in interpolation. -!> ISTATUS = 14: Could not find the other two cell centers of the triangle that contains this lat/lon -!> ISTATUS = 15: Cell centers of the triangle fall in the lateral boundary zone -!> ISTATUS = 16: Don't know how to do vertical velocity for now -!> ISTATUS = 17: Unable to compute pressure values -!> ISTATUS = 18: altitude illegal -!> ISTATUS = 19: could not compute u using RBF code -!> ISTATUS = 101: Internal error; reached end of subroutine without -!> finding an applicable case. -!> ISTATUS = 201: Reject observation from user specified pressure level -!> ISTATUS = 988: pressure is not monotonically descreased with level. -!> -!> Debugging options: -!> 0: No prints, whatsoever. -!> 1: Print only when each obs is rejected. -!> 2: Print the basic info on each obs. -!> 3: Print the base location of the localization. -!> 5: Print info on the localized obs and the model grids -!> 9: Print info on wind updates -!> 10: Print info on each obs in detail. -!> 11: Print info on localized states -!> 12: Print info on each ensemble member in detail. -!> 13: Print detailed info on triangle and level search - - -subroutine model_interpolate(state_handle, ens_size, location, obs_type, expected_obs, istatus) - -! passed variables - -type(ensemble_type), intent(in) :: state_handle -integer, intent(in) :: ens_size -type(location_type), intent(in) :: location -integer, intent(in) :: obs_type -real(r8), intent(out) :: expected_obs(ens_size) -integer, intent(out) :: istatus(ens_size) - -! local storage - -type(location_type) :: location_tmp(ens_size) -integer :: ivar, obs_kind -integer :: tvars(3) -integer :: cellid -logical :: goodkind, surface_obs -real(r8) :: lpres(ens_size), values(3, ens_size) -real(r8) :: llv(3) ! lon/lat/vert -integer :: e, verttype - -if ( .not. module_initialized ) call static_init_model - -expected_obs = MISSING_R8 -istatus = 0 ! must be positive (and integer) - -! rename for sanity - we can't change the argument names -! to this subroutine, but this really is a kind. -obs_kind = obs_type - -! write the location information into a string for error messages and debugging -call write_location(0,location,charstring=locstring) - -llv = get_location(location) -verttype = nint(query_location(location)) -surface_obs = (verttype == VERTISSURFACE) - -! this routine returns the cellid for a global mpas grid, same as -! find_closest_cell_center(). -! for a regional grid it only returns a good cellid if the closest cell center -! AND the other 2 triangle points surrounding this location are completely inside -! the grid and none of the vertices are in the boundary region. -cellid = cell_ok_to_interpolate(location) -if (cellid < 1) then - !print *, 'model_interpolate: lon/lat is outside the domain: ', llv(1), llv(2) - istatus = 11 - goto 100 -endif - -if (debug > 1 .and. do_output()) & - print *, 'model_interpolate for obs_kind:',& - trim(get_name_for_quantity(obs_kind)),' at ',trim(locstring),' cellid:',cellid - -! FIXME see issue #96 - remove all but surface elevation -! pass surface variables for rttov - it should be ok as these are diagnostic (except for surface elevation). -if((obs_kind == QTY_SURFACE_PRESSURE) .or. (obs_kind == QTY_SURFACE_ELEVATION) .or. & - (obs_kind == QTY_2M_TEMPERATURE) .or. (obs_kind == QTY_2M_SPECIFIC_HUMIDITY) .or. & - (obs_kind == QTY_SKIN_TEMPERATURE) .or. (obs_kind == QTY_SURFACE_TYPE) .or. & - (obs_kind == QTY_CLOUD_FRACTION) .or. & - (obs_kind == QTY_10M_U_WIND_COMPONENT) .or. (obs_kind == QTY_10M_V_WIND_COMPONENT)) then - istatus = 0 - if ( debug > 1 .and. do_output()) print *, & - 'model_interpolate: pass sfc_elev_max_diff check for ',trim(get_name_for_quantity(obs_kind)) -else -! Reject obs if the station height is far way from the model terrain. - if(is_vertical(location, "SURFACE").and. sfc_elev_max_diff >= 0) then - if(abs(llv(3) - zGridFace(1,cellid)) > sfc_elev_max_diff) then - istatus = 12 - if (debug > 0 .and. do_output()) then - print*, 'model_interpolate: Skip due to abs(dz) > sfc_elev_max_diff:', & - llv(3)-zGridFace(1,cellid), trim(get_name_for_quantity(obs_kind)),' at ',trim(locstring) - goto 100 + case (QTY_U_WIND_COMPONENT) + if (use_u_for_wind .and. has_edge_u) then + call compute_u_with_rbf(state_handle, ens_size, location, .TRUE., expected_obs, istatus) + else + tvars(1) = get_varid_from_kind(anl_domid, qty) + call compute_scalar_with_barycentric(state_handle, ens_size, location, 1, tvars(1), expected_obs, istatus) + endif + + case (QTY_V_WIND_COMPONENT) + if (use_u_for_wind .and. has_edge_u) then + call compute_u_with_rbf(state_handle, ens_size, location, .FALSE., expected_obs, istatus) + else + tvars(1) = get_varid_from_kind(anl_domid, qty) + call compute_scalar_with_barycentric(state_handle, ens_size, location, 1, tvars(1), expected_obs, istatus) endif - endif - endif -endif - -! see if observation quantity is in the state vector. this sets an -! error code and returns without a fatal error if answer is no. -! exceptions: the state vector has potential temp, but we can -! compute sensible temperature from pot_temp, rho, and qv. -! also there are options for the winds because mpas has both -! winds on the cell edges (normal only) and reconstructed winds -! at the cell centers (U,V). there are namelist options to control -! which to use if both are in the state vector. -! As another note, mpas defines the model variable qv as water vapor -! mixing ratio while it defines q2 as 2-meter specific humidity, -! not 2-m water vapor mixing ratio, so q2 should be specified as -! QTY_2M_SPECIFIC_HUMIDITY in mpas_state_variables in &mpas_vars_nml. - -! is this field in the state? -ivar = get_progvar_index_from_kind(obs_kind) -if (ivar > 0) then - goodkind = .true. ! yes - -else - goodkind = .false. ! but check for exceptions - - ! exceptions if the kind isn't directly - ! a field in the state vector: - select case (obs_kind) - case (QTY_TEMPERATURE, QTY_2M_TEMPERATURE) - goodkind = .true. - case (QTY_SURFACE_ELEVATION, QTY_GEOPOTENTIAL_HEIGHT) - goodkind = .true. - case (QTY_PRESSURE) ! surface pressure should be in the state - goodkind = .true. - case (QTY_SKIN_TEMPERATURE, QTY_SURFACE_TYPE, QTY_CLOUD_FRACTION) - goodkind = .true. - case (QTY_SPECIFIC_HUMIDITY, QTY_2M_SPECIFIC_HUMIDITY) - goodkind = .true. - case (QTY_U_WIND_COMPONENT, QTY_V_WIND_COMPONENT) - ! if the reconstructed winds at the cell centers aren't there, - ! we can use the edge normal winds, if the user allows it. - if (get_progvar_index_from_kind(QTY_EDGE_NORMAL_SPEED) > 0 & - .and. use_u_for_wind) goodkind = .true. - end select -endif - -if (debug > 10 .and. do_output()) & - print *, 'model_interpolate: ivar, goodkind? ', ivar, goodkind,' ',& - trim(get_name_for_quantity(obs_kind)),' at ',trim(locstring) - -! this kind is not in the state vector and it isn't one of the exceptions -! that we know how to handle. -if (.not. goodkind) then - istatus(:) = 88 - if (debug > 4 .and. do_output()) print *, 'model_interpolate: kind rejected', obs_kind - goto 100 -endif - -! Reject obs above a user specified pressure level. -! this is expensive - only do it if users want to reject observations -! at the top of the model. negative values mean ignore this test. -if (.not.surface_obs) then + case (QTY_TEMPERATURE) + ! need to get potential temp, pressure, qv here, but can + ! use same weights, so push all three types into the subroutine. + tvars(1) = get_varid_from_kind(anl_domid, QTY_POTENTIAL_TEMPERATURE) + tvars(2) = get_varid_from_kind(anl_domid, QTY_DENSITY) + tvars(3) = get_varid_from_kind(anl_domid, QTY_VAPOR_MIXING_RATIO) + + call compute_scalar_with_barycentric(state_handle, ens_size, location, 3, tvars, values, istatus) + if (fails(istatus)) then + istatus = TK_COMPUTATION_ERROR + return ! early return to not pass nonsense values to theta_to_tk + endif + ! convert pot_temp, density, vapor mixing ratio into sensible temperature + expected_obs(:) = theta_to_tk(ens_size, values(1, :), values(2, :), values(3, :), istatus) -if (highest_obs_pressure_mb > 0.0) then - if (.not.(is_vertical(location, "SURFACE")) .and. (.not.(is_vertical(location, "UNDEFINED")))) then + case (QTY_PRESSURE) + call compute_pressure_at_loc(state_handle, ens_size, location, expected_obs, istatus) + if (fails(istatus)) istatus = PRESSURE_COMPUTATION_ERROR - call compute_pressure_at_loc(state_handle, ens_size, location, lpres, istatus) - where (lpres < highest_obs_pressure_mb * 100.0_r8) - ! Exclude from assimilation the obs above a user specified level - istatus(:) = 201 - end where + case (QTY_GEOPOTENTIAL_HEIGHT) + location_tmp = location + call convert_vert(state_handle, ens_size, location_tmp, QTY_GEOPOTENTIAL_HEIGHT, VERTISHEIGHT, istatus) - if (debug > 10 .and. do_output()) then do e = 1, ens_size - if (istatus(e) == 201) print *, 'ens ', e, ' rejected, pressure < upper limit', lpres(e), highest_obs_pressure_mb + if(istatus(e) == 0) expected_obs(e) = query_location(location_tmp(e), 'VLOC') enddo - endif - if (all(istatus /= 0)) goto 100 ! if everyone has failed, we can quit - - endif -endif - -if (debug > 9 .and. do_output()) then - print *, 'high pressure test is passed, ready to interpolate kind ', obs_kind -endif - -endif !(.not.surface_obs) then - -! Not prepared to do W interpolation at this time -if(obs_kind == QTY_VERTICAL_VELOCITY) then - if (debug > 0 .and. do_output()) print *, 'model_interpolate: code does not handle vertical velocity yet' - istatus(:) = 16 - goto 100 -endif - -! winds -if ((obs_kind == QTY_U_WIND_COMPONENT .or. & - obs_kind == QTY_V_WIND_COMPONENT) .and. has_edge_u .and. use_u_for_wind) then - if (obs_kind == QTY_U_WIND_COMPONENT) then - ! return U - call compute_u_with_rbf(state_handle, ens_size, location, .TRUE., expected_obs, istatus) - else - ! return V - call compute_u_with_rbf(state_handle, ens_size, location, .FALSE., expected_obs, istatus) - endif - if (debug > 9 .and. do_output()) print *, 'model_interpolate: u_with_rbf ', obs_kind, istatus(1), expected_obs(1) - if (all(istatus /= 0)) goto 100 ! if any member has failed, we can exit -else if (obs_kind == QTY_TEMPERATURE) then - ! need to get potential temp, pressure, qv here, but can - ! use same weights, so push all three types into the subroutine. - tvars(1) = get_progvar_index_from_kind(QTY_POTENTIAL_TEMPERATURE) - tvars(2) = get_progvar_index_from_kind(QTY_DENSITY) - tvars(3) = get_progvar_index_from_kind(QTY_VAPOR_MIXING_RATIO) - - call compute_scalar_with_barycentric(state_handle, ens_size, location, 3, tvars, values, istatus) - if (all(istatus /= 0)) goto 100 - - ! convert pot_temp, density, vapor mixing ratio into sensible temperature - expected_obs(:) = theta_to_tk(ens_size, values(1, :), values(2, :), values(3, :), istatus) - - if (debug > 9 .and. do_output()) & - print *, 'model_interpolate: TEMPERATURE ', istatus(1), expected_obs(1), trim(locstring) - if ( all(istatus /= 0 ) ) goto 100 - -else if (obs_kind == QTY_PRESSURE) then - call compute_pressure_at_loc(state_handle, ens_size, location, expected_obs, istatus) - - if (debug > 9 .and. do_output()) & - print *, 'model_interpolate: PRESSURE ', istatus(1), expected_obs(1), trim(locstring) - if ( all(istatus /= 0 ) ) goto 100 - -else if (obs_kind == QTY_GEOPOTENTIAL_HEIGHT) then - location_tmp = location - call convert_vert_distrib(state_handle, ens_size, location_tmp, QTY_GEOPOTENTIAL_HEIGHT, VERTISHEIGHT, istatus) - - do e = 1, ens_size - if(istatus(e) == 0) expected_obs(e) = query_location(location_tmp(e), 'VLOC') - enddo - - if ( all(istatus /= 0 ) ) goto 100 - -else if (obs_kind == QTY_VAPOR_MIXING_RATIO .or. obs_kind == QTY_2M_SPECIFIC_HUMIDITY .or. & - obs_kind == QTY_CLOUDWATER_MIXING_RATIO .or. obs_kind == QTY_RAINWATER_MIXING_RATIO .or. & - obs_kind == QTY_ICE_MIXING_RATIO .or. obs_kind == QTY_SNOW_MIXING_RATIO .or. & - obs_kind == QTY_GRAUPEL_MIXING_RATIO .or. obs_kind == QTY_CLOUD_FRACTION ) then - tvars(1) = get_progvar_index_from_kind(obs_kind) - call compute_scalar_with_barycentric(state_handle, ens_size, location, 1, tvars(1), values(1,:), istatus) - expected_obs = values(1, :) - - ! Don't accept negative hydrometeors - do e = 1, ens_size - if(istatus(e) == 0) then - expected_obs(e) = max(values(1, e),0.0_r8) - if(obs_kind == QTY_CLOUD_FRACTION) expected_obs(e) = min(values(1, e),1.0_r8) + ! Things that cannot be negative + case (QTY_VAPOR_MIXING_RATIO, & + QTY_CLOUDWATER_MIXING_RATIO, & + QTY_ICE_MIXING_RATIO, & + QTY_GRAUPEL_MIXING_RATIO, & + QTY_2M_SPECIFIC_HUMIDITY, & + QTY_RAINWATER_MIXING_RATIO, & + QTY_SNOW_MIXING_RATIO, & + QTY_CLOUD_FRACTION ) + + tvars(1) = get_varid_from_kind(anl_domid, qty) + call compute_scalar_with_barycentric(state_handle, ens_size, location, 1, tvars(1), expected_obs, istatus) + if (fails(istatus)) then + istatus = GENERAL_ERROR + return endif - enddo - - if (debug > 9 .and. do_output()) then - print *, 'model_interpolate: obs_kind,name,istatus,expected_obs,location =', & - obs_kind,trim(get_name_for_quantity(obs_kind)),istatus(1),expected_obs(1),trim(locstring) - endif - if ( all(istatus /= 0 ) ) goto 100 - -else if (obs_kind == QTY_SPECIFIC_HUMIDITY) then - tvars(1) = get_progvar_index_from_kind(QTY_VAPOR_MIXING_RATIO) - call compute_scalar_with_barycentric(state_handle, ens_size, location, 1, tvars(1), values(1,:), istatus) - expected_obs = values(1, :) - - ! compute vapor pressure, then: sh = vp / (1.0 + vp) - do e = 1, ens_size - if (istatus(e) == 0) then - if (expected_obs(e) >= 0.0_r8) then - expected_obs(e) = expected_obs(e) / (1.0_r8 + expected_obs(e)) - else - expected_obs(e) = 1.0e-12_r8 - endif - endif - enddo - if ( all(istatus /= 0 ) ) goto 100 - - if (debug > 9 .and. do_output()) & - print *, 'model_interpolate: ', trim(get_name_for_quantity(obs_kind)), istatus(1), & - expected_obs(1), trim(locstring) - -! Only for the variables NOT included in the dart state vector. -! Anything included in the state_handle should not go into the case right down here. -! ex. QTY_SURFACE_PRESSURE will be taken care of in the else statement (one further down) -! as a generic interpolation case thru compute_scalar_with_barycentric -! because that obs_kind is included in mpas_state_variables in &mpas_vars_nml (e.g. state vector). -else if (obs_kind == QTY_SURFACE_ELEVATION .or. & - obs_kind == QTY_SKIN_TEMPERATURE .or. obs_kind == QTY_SURFACE_TYPE ) then - - if (obs_kind == QTY_SURFACE_ELEVATION) then - call compute_surface_data_with_barycentric(zGridFace(1,:), location, expected_obs(1), istatus(1)) - else if (obs_kind == QTY_SKIN_TEMPERATURE) then - call compute_surface_data_with_barycentric(skintemp(:), location, expected_obs(1), istatus(1)) - else if (obs_kind == QTY_SURFACE_TYPE) then - call get_surftype(nCells,surftype) - call compute_surface_data_with_barycentric(surftype(:)*1.0_r8, location, expected_obs(1), istatus(1)) - endif + expected_obs = max(expected_obs, 0.0_r8) + if (qty == QTY_CLOUD_FRACTION) expected_obs = min(expected_obs, 1.0_r8) - expected_obs(2:ens_size) = expected_obs(1) - istatus(2:ens_size) = istatus(1) + case(QTY_SPECIFIC_HUMIDITY) + tvars(1) = get_varid_from_kind(anl_domid, QTY_VAPOR_MIXING_RATIO) + call compute_scalar_with_barycentric(state_handle, ens_size, location, 1, tvars(1), expected_obs, istatus) - if (debug > 9 .and. do_output()) & - print *, 'model_interpolate: ', trim(get_name_for_quantity(obs_kind)), ' ', istatus(1), & - expected_obs(1), trim(locstring) - if ( all(istatus /= 0 ) ) goto 100 - -else - ! all other kinds come here. - ! direct interpolation: kind is in the state vector and no clamping or other conversions needed - - tvars(1) = ivar - call compute_scalar_with_barycentric(state_handle, ens_size, location, 1, tvars(1), values(1,:), istatus) - expected_obs = values(1, :) - - if (debug > 9 .and. do_output()) & - print*, 'generic interpolation in compute_scalar_with_barycentric: ', obs_kind, & - trim(get_name_for_quantity(obs_kind)) - if (debug > 11 .and. do_output()) then - do e = 1, ens_size - print*, 'member ',e, ' istatus(e)=',istatus(e), ', expected_obs(e)=', expected_obs(e) - enddo - endif - if ( all(istatus /= 0 ) ) goto 100 + ! compute vapor pressure, then: sh = vp / (1.0 + vp) + do e = 1, ens_size + if (istatus(e) == 0) then + if (expected_obs(e) >= 0.0_r8) then + expected_obs(e) = expected_obs(e) / (1.0_r8 + expected_obs(e)) + else + expected_obs(e) = 1.0e-12_r8 + endif + endif + enddo -endif + ! qtys not included in the dart state vector, but static across the ensemble + case (QTY_SURFACE_ELEVATION) + call compute_surface_data_with_barycentric(zGridFace(1,:), location, single_expected_obs, single_istatus) + expected_obs(:) = single_expected_obs + istatus(:) = single_istatus + + case (QTY_SKIN_TEMPERATURE) + call compute_surface_data_with_barycentric(skintemp(:), location, single_expected_obs, single_istatus) + expected_obs(:) = single_expected_obs + istatus(:) = single_istatus + + case (QTY_SURFACE_TYPE) + call get_surftype(nCells,surftype) + call compute_surface_data_with_barycentric(surftype(:)*1.0_r8, location, single_expected_obs, single_istatus) + expected_obs(:) = single_expected_obs + istatus(:) = single_istatus + + case default ! regular interpolation + tvars(1) = get_varid_from_kind(anl_domid, qty) + call compute_scalar_with_barycentric(state_handle, ens_size, location, 1, tvars(1), expected_obs, istatus) +end select -100 continue +end subroutine model_interpolate -! this is for debugging - when we're confident the code is returning -! consistent values and rc codes, both these tests can be removed for speed. -! also optionally check for generated NaNs for now. +!------------------------------------------------------------------ +function fails(e) result(fail) -do e = 1, ens_size +integer, intent(in) :: e(:) +logical :: fail - if ((istatus(e) < 0) .or. & - (istatus(e) /= 0 .and. expected_obs(e) /= MISSING_R8) .or. & - (istatus(e) == 0 .and. expected_obs(e) == MISSING_R8)) then +fail = (any(e /= 0)) - write(string2,*) 'member ',e,' obs_kind', obs_kind,' value = ', expected_obs(e), & - ' istatus = ', istatus(e), ' cellid:',cellid - write(string3,*) trim(locstring) +end function fails - if (istatus(e) < 0) then - write(string1,*) 'interp routine returned a negative status which is an illegal value' - else if (istatus(e) /= 0 .and. expected_obs(e) /= MISSING_R8) then - write(string1,*) 'interp routine returned a bad status but not a MISSING_R8 value' - expected_obs(e) = MISSING_R8 - else - write(string1,*) 'interp routine returned a good status but set value to MISSING_R8' - endif +!------------------------------------------------------------------ +function qty_ok_to_interpolate(qty) result(qty_ok) - call error_handler(E_ALLMSG,'model_interpolate', string1, source,revision,revdate, & - text2=string2, text3=string3) - endif +integer, intent(in) :: qty +logical :: qty_ok - ! the only portable, reliable test for NaNs we know - if some number is neither - ! less than nor equal to/greater than 0, it must be a NaN. all numerical comparisons - ! fail if one or more of the operands are NaN. - - if (.not. expected_obs(e) < 0 .and. .not. expected_obs(e) >= 0) then - write(string1,*) 'Skip member ', e, ' expected obs may be NaN: ', expected_obs(e), & - ' obs_kind', obs_kind,' istatus = ', istatus(e) - write(string2,*) 'cellid:',cellid,' location ', trim(locstring) - expected_obs(e) = MISSING_R8 - istatus(e) = 99 - call error_handler(E_ERR,'model_interpolate', string1, source,revision,revdate,text2=string2) - endif +integer :: varid - if (debug > 11 .and. do_output()) then - write(string2,*) 'Completed for member ',e,' obs_kind', obs_kind,' expected_obs = ', expected_obs(e) - write(string3,*) 'istatus = ', istatus(e), ' at ', trim(locstring) - call error_handler(E_ALLMSG, 'model_interpolate', string2, source, revision, revdate, text2=string3) - endif -enddo +qty_ok = .false. +varid = get_varid_from_kind(anl_domid, qty) -end subroutine model_interpolate +if (varid > 0) then ! in the state vector + qty_ok = .true. + return +endif +select case (qty) +case (QTY_TEMPERATURE, QTY_2M_TEMPERATURE) + qty_ok = .true. +case (QTY_SURFACE_ELEVATION, QTY_GEOPOTENTIAL_HEIGHT) + qty_ok = .true. +case (QTY_PRESSURE) ! surface pressure should be in the state !HK @todo check "should" + qty_ok = .true. +case (QTY_SKIN_TEMPERATURE, QTY_SURFACE_TYPE, QTY_CLOUD_FRACTION) + qty_ok = .true. +case (QTY_SPECIFIC_HUMIDITY, QTY_2M_SPECIFIC_HUMIDITY) + qty_ok = .true. +case (QTY_U_WIND_COMPONENT, QTY_V_WIND_COMPONENT) + if (get_varid_from_kind(anl_domid, QTY_EDGE_NORMAL_SPEED) > 0 .and. use_u_for_wind) qty_ok = .true. +case default + qty_ok = .false. +end select +end function qty_ok_to_interpolate !------------------------------------------------------------------ subroutine nc_write_model_atts(ncid, domain_id) @@ -1503,44 +925,19 @@ subroutine nc_write_model_atts(ncid, domain_id) real(r8), allocatable :: data1d(:) integer :: ncid2 - -!------------------------------------------------------------------------------- -! put file into define mode. -!------------------------------------------------------------------------------- - call nc_begin_define_mode(ncid) -!------------------------------------------------------------------------------- -! Write Global Attributes -!------------------------------------------------------------------------------- - call nc_add_global_creation_time(ncid) - call nc_add_global_attribute(ncid, "model_source", source) -call nc_add_global_attribute(ncid, "model_revision", revision) -call nc_add_global_attribute(ncid, "model_revdate", revdate) - call nc_add_global_attribute(ncid, "model", "MPAS_ATM") - -!---------------------------------------------------------------------------- -! if not adding grid info, return here -!---------------------------------------------------------------------------- - -if (.not. write_grid_to_diag_files) then +if (.not. write_grid_to_diag_files) then ! early return, no need to write grid info call nc_end_define_mode(ncid) call nc_synchronize_file(ncid) return endif -!---------------------------------------------------------------------------- -! Everything below here is static grid info -!---------------------------------------------------------------------------- - -!---------------------------------------------------------------------------- ! Dimensions -!---------------------------------------------------------------------------- - call nc_define_dimension(ncid, 'nCells', nCells, routine) call nc_define_dimension(ncid, 'nEdges', nEdges, routine) call nc_define_dimension(ncid, 'nVertLevels', nVertLevels, routine) @@ -1551,9 +948,8 @@ subroutine nc_write_model_atts(ncid, domain_id) call nc_define_dimension(ncid, 'nVertices', nVertices, routine) call nc_define_dimension(ncid, 'VertexDegree', VertexDegree, routine) -!---------------------------------------------------------------------------- ! Coordinate Variables and the Attributes -!---------------------------------------------------------------------------- + ! Cell Longitudes call nc_define_real_variable(ncid, 'lonCell', 'nCells', routine) call nc_add_attribute_to_variable(ncid, 'lonCell', 'long_name', 'cell center longitudes', routine) @@ -1611,16 +1007,9 @@ subroutine nc_write_model_atts(ncid, domain_id) call nc_define_real_variable(ncid, 'areaCell', 'nCells', routine) -!---------------------------------------------------------------------------- -! Finished with dimension/variable definitions, must end 'define' mode to fill. -!---------------------------------------------------------------------------- call nc_end_define_mode(ncid) -!---------------------------------------------------------------------------- -! Fill the coordinate variables -!---------------------------------------------------------------------------- - call nc_put_variable(ncid, 'lonCell', lonCell, routine) call nc_put_variable(ncid, 'latCell', latCell, routine) @@ -1635,11 +1024,12 @@ subroutine nc_write_model_atts(ncid, domain_id) call nc_put_variable(ncid, 'verticesOnCell', verticesOnCell, routine) call nc_put_variable(ncid, 'cellsOnVertex', cellsOnVertex, routine) +call nc_synchronize_file(ncid) ! Flush the buffer and leave netCDF file open + !---------------------------------------------------------------------------- ! DART has not read these in, so we have to read them from the input file ! and copy them to the DART output file. !---------------------------------------------------------------------------- - ncid2 = nc_open_file_readonly(init_template_filename, routine) allocate(data1d(nCells)) @@ -1670,10 +1060,6 @@ subroutine nc_write_model_atts(ncid, domain_id) call nc_close_file(ncid2) -!------------------------------------------------------------------------------- -! Flush the buffer and leave netCDF file open -!------------------------------------------------------------------------------- -call nc_synchronize_file(ncid) end subroutine nc_write_model_atts @@ -1681,36 +1067,15 @@ end subroutine nc_write_model_atts function get_model_size() -! Returns the size of the model as an integer. -! Required for all applications. - integer(i8) :: get_model_size if ( .not. module_initialized ) call static_init_model -!print *, 'model_size = ', model_size get_model_size = model_size end function get_model_size !------------------------------------------------------------------ -!> Returns the number of variables as an integer. - -function get_num_vars() - -integer :: get_num_vars - -if ( .not. module_initialized ) call static_init_model - -get_num_vars = nfields - -end function get_num_vars - -!------------------------------------------------------------------ -!> Returns the the time step of the model; the smallest increment -!> in time that the model is capable of advancing the state in a given -!> implementation. This interface is required for all applications. - function shortest_time_between_assimilations() type(time_type) :: shortest_time_between_assimilations @@ -1772,25 +1137,8 @@ subroutine pert_model_copies(ens_handle, ens_size, pert_amp, interf_provided) integer :: i, j integer(i8), allocatable :: var_list(:) - - -! Perturbs a model state copies for generating initial ensembles. -! A model may choose to provide a NULL INTERFACE by returning -! .false. for the interf_provided argument. This indicates to -! the filter that if it needs to generate perturbed states, -! it may do so by adding a perturbation to each model state -! variable independently. The interf_provided argument -! should be returned as .true. if the model wants to do its own -! perturbing of states. - interf_provided = .true. -!>@todo If MPAS ever supports more than a single domain then -!>look at the wrf model_mod code for how to change this. you -!>have to separate out the total number of variables across -!>all domains for the min/max part, and then loop over only -!>the number of variables in each domain in the second part. - num_variables = get_num_variables(anl_domid) ! Get min and max of each variable in each domain @@ -1863,93 +1211,62 @@ subroutine get_close_obs(gc, base_loc, base_type, locs, loc_qtys, loc_types, & integer :: ztypeout -integer :: t_ind, istatus1, istatus2, k, istat_arr(1) +integer :: t_ind, istatus(1), k integer :: base_which, local_obs_which, base_qty real(r8), dimension(3) :: base_llv, local_obs_llv ! lon/lat/vert type(location_type) :: local_obs_loc, location_arr(1) -! timing -real(digits12) :: t_base, t_base2, interval - -real(r8) :: hor_dist -hor_dist = 1.0e9_r8 - -! Initialize variables to missing status num_close = 0 -istatus1 = 0 -istatus2 = 0 ! Convert base_obs vertical coordinate to requested vertical coordinate if necessary - base_llv = get_location(base_loc) base_which = nint(query_location(base_loc)) ztypeout = vert_localization_coord if (vertical_localization_on()) then - if (base_llv(3) == MISSING_R8) then - istatus1 = 1 + if (base_llv(3) == MISSING_R8) then !HK @todo what about VERTISUNDEF? + return else if (base_which /= vert_localization_coord .and. base_which /= VERTISUNDEF) then base_qty = get_quantity_for_type_of_obs(base_type) location_arr(1) = base_loc - call convert_vert_distrib(state_handle, 1, location_arr, base_qty, vert_localization_coord, istat_arr) - istatus1 = istat_arr(1) + call convert_vert(state_handle, 1, location_arr, base_qty, vert_localization_coord, istatus) + if (istatus(1) /= 0) return base_loc = location_arr(1) - if(debug > 4 .and. do_output()) then - call write_location(0,base_loc,charstring=string1) - call error_handler(E_ALLMSG, 'get_close_obs: base_loc',string1,source, revision, revdate) - endif endif endif -if (istatus1 == 0) then - - ! Loop over potentially close subset of obs priors or state variables - ! This way, we are decreasing the number of distance computations that will follow. - ! This is a horizontal-distance operation and we don't need to have the relevant vertical - ! coordinate information yet (for locs). - call loc_get_close_obs(gc, base_loc, base_type, locs, loc_qtys, loc_types, & - num_close, close_ind) - - do k = 1, num_close - - t_ind = close_ind(k) - local_obs_loc = locs(t_ind) - local_obs_which = nint(query_location(local_obs_loc)) - - ! Convert local_obs vertical coordinate to requested vertical coordinate if necessary. - ! This should only be necessary for obs priors, as state location information already - ! contains the correct vertical coordinate (filter_assim's call to get_state_meta_data). - if ((base_which /= VERTISUNDEF) .and. (vertical_localization_on())) then - if (local_obs_which /= vert_localization_coord .and. local_obs_which /= VERTISUNDEF) then - location_arr(1) = local_obs_loc - call convert_vert_distrib(state_handle, 1, location_arr, loc_qtys(t_ind), vert_localization_coord, istat_arr) - istatus2 = istat_arr(1) - locs(t_ind) = location_arr(1) - else - istatus2 = 0 - endif - endif +! Loop over potentially close subset of obs priors or state variables +call loc_get_close_obs(gc, base_loc, base_type, locs, loc_qtys, loc_types, num_close, close_ind) - if (present(dist)) then - ! Compute distance - set distance to a very large value if vert coordinate is missing - ! or vert_interpolate returned error (istatus2=1) - local_obs_llv = get_location(local_obs_loc) - if ( (vertical_localization_on() .and. & - (local_obs_llv(3) == MISSING_R8)) .or. (istatus2 /= 0) ) then - dist(k) = 1.0e9_r8 - else - dist(k) = get_dist(base_loc, locs(t_ind), base_type, loc_qtys(t_ind)) - endif +do k = 1, num_close + + t_ind = close_ind(k) + local_obs_loc = locs(t_ind) + local_obs_which = nint(query_location(local_obs_loc)) + + ! Convert local_obs vertical coordinate to requested vertical coordinate if necessary. + if ((base_which /= VERTISUNDEF) .and. (vertical_localization_on())) then + if (local_obs_which /= vert_localization_coord .and. local_obs_which /= VERTISUNDEF) then + location_arr(1) = local_obs_loc + call convert_vert(state_handle, 1, location_arr, loc_qtys(t_ind), vert_localization_coord, istatus) + locs(t_ind) = location_arr(1) + else + istatus = 0 endif + else + istatus = 0 + endif - enddo -endif + if (present(dist)) then + if (istatus(1) == 0) then + dist(k) = get_dist(base_loc, locs(t_ind), base_type, loc_qtys(t_ind)) + else + dist(k) = 1.0e9_r8 + endif + endif -if ((debug > 2) .and. do_output()) then - call write_location(0,base_loc,charstring=string2) - print *, 'get_close_obs: nclose, base_loc ', num_close, trim(string2) -endif +enddo end subroutine get_close_obs @@ -1961,10 +1278,6 @@ end subroutine get_close_obs subroutine get_close_state(gc, base_loc, base_type, locs, loc_qtys, loc_indx, & num_close, close_ind, dist, state_handle) -!>@todo FIXME this is working on state vector items. if a vertical -!>conversion is needed, it doesn't need to interpolate. it can compute -!>the location using the logic that get_state_meta_data() uses. - type(get_close_type), intent(in) :: gc type(location_type), intent(inout) :: base_loc, locs(:) integer, intent(in) :: base_type, loc_qtys(:) @@ -1975,955 +1288,71 @@ subroutine get_close_state(gc, base_loc, base_type, locs, loc_qtys, loc_indx, & integer :: ztypeout -integer :: t_ind, istatus1, istatus2, k, istat_arr(1) +integer :: t_ind, istatus(1), k integer :: base_which, local_obs_which, base_qty real(r8), dimension(3) :: base_llv, local_obs_llv ! lon/lat/vert type(location_type) :: local_obs_loc, location_arr(1) -! timing -real(digits12) :: t_base, t_base2, interval - - -real(r8) :: hor_dist -hor_dist = 1.0e9_r8 - -! Initialize variables to missing status num_close = 0 -istatus1 = 0 -istatus2 = 0 ! Convert base_obs vertical coordinate to requested vertical coordinate if necessary - base_llv = get_location(base_loc) base_which = nint(query_location(base_loc)) ztypeout = vert_localization_coord if (vertical_localization_on()) then - if (base_llv(3) == MISSING_R8) then - istatus1 = 1 + if (base_llv(3) == MISSING_R8) then !HK @todo what about VERTISUNDEF? + return else if (base_which /= vert_localization_coord .and. base_which /= VERTISUNDEF) then base_qty = get_quantity_for_type_of_obs(base_type) location_arr(1) = base_loc - call convert_vert_distrib(state_handle, 1, location_arr, base_qty, vert_localization_coord, istat_arr) - istatus1 = istat_arr(1) + call convert_vert(state_handle, 1, location_arr, base_qty, vert_localization_coord, istatus) + if (istatus(1) /= 0) return base_loc = location_arr(1) - if(debug > 9 .and. do_output()) then - print*, 'get_close_state: istatus1 from convert_vert_distrib for base_loc, base_which = ',& - istatus1,base_which - call write_location(0,base_loc,charstring=string1) - call error_handler(E_ALLMSG, 'get_close_state: base_loc',string1,source, revision, revdate) - endif - endif -endif - -if (istatus1 == 0) then - - ! Loop over potentially close subset of obs priors or state variables - ! This way, we are decreasing the number of distance computations that will follow. - ! This is a horizontal-distance operation and we don't need to have the relevant vertical - ! coordinate information yet (for locs). - call loc_get_close_state(gc, base_loc, base_type, locs, loc_qtys, loc_indx, & - num_close, close_ind) - - do k = 1, num_close - - t_ind = close_ind(k) - local_obs_loc = locs(t_ind) - local_obs_which = nint(query_location(local_obs_loc)) - - ! Convert local_obs vertical coordinate to requested vertical coordinate if necessary. - ! This should only be necessary for obs priors, as state location information already - ! contains the correct vertical coordinate (filter_assim's call to get_state_meta_data). - if ((base_which /= VERTISUNDEF) .and. (vertical_localization_on())) then - if (local_obs_which /= vert_localization_coord .and. local_obs_which /= VERTISUNDEF) then - location_arr(1) = local_obs_loc - call convert_vert_distrib_state(state_handle, 1, location_arr, loc_qtys(t_ind), & - loc_indx(t_ind), vert_localization_coord, istat_arr) - istatus2 = istat_arr(1) - locs(t_ind) = location_arr(1) - if(istatus2 /= 0 .and. debug > 9 .and. do_output()) then - call write_location(0,local_obs_loc,charstring=string1) - call error_handler(E_ALLMSG, 'get_close_state: local_obs_loc',string1,source, revision, revdate) - endif - else - istatus2 = 0 - endif - endif - - if (present(dist)) then - ! Compute distance - set distance to a very large value if vert coordinate is missing - ! or vert_interpolate returned error (istatus2=1) - local_obs_llv = get_location(local_obs_loc) - if ( (vertical_localization_on() .and. & - (local_obs_llv(3) == MISSING_R8)) .or. (istatus2 /= 0) ) then - dist(k) = 1.0e9_r8 - else - dist(k) = get_dist(base_loc, locs(t_ind), base_type, loc_qtys(t_ind)) - endif - endif - - enddo -endif - -if ((debug > 2) .and. do_output()) then - call write_location(0,base_loc,charstring=string2) - print *, 'get_close_state: nclose, base_loc ', num_close, trim(string2) -endif - -end subroutine get_close_state - - -!================================================================== -! The (model-specific) additional public interfaces come next -! (these are not required by dart but are used by other programs) -!================================================================== - -subroutine get_init_template_filename( filename ) - -! return the name of the template filename that was set -! in the model_nml namelist (ex. init.nc) - -character(len=*), intent(OUT) :: filename - -if ( .not. module_initialized ) call static_init_model - -filename = trim(init_template_filename) - -end subroutine get_init_template_filename - - -!------------------------------------------------------------------- -! modify what static_init_model does. this *must* be called before -! calling static_init_model(). -! the boundary file variables are fixed by the model and so we -! don't allow the user to set them via namelist - -subroutine set_lbc_variables(template_filename) - -character(len=*), intent(in) :: template_filename - -integer :: ncid - -bdy_template_filename = template_filename - -! this initial list always exists. hardcode them for now, -! and query to see if the reconstructed winds are there or not. - -lbc_variables(1) = 'lbc_qc' -lbc_variables(2) = 'lbc_qr' -lbc_variables(3) = 'lbc_qv' -lbc_variables(4) = 'lbc_rho' -lbc_variables(5) = 'lbc_theta' -lbc_variables(6) = 'lbc_u' -lbc_variables(7) = 'lbc_w' - -ncid = nc_open_file_readonly(template_filename, 'set_lbc_variables') -if (nc_variable_exists(ncid, 'lbc_ur')) then - lbc_variables(8) = 'lbc_ur' - lbc_variables(9) = 'lbc_vr' - lbc_file_has_reconstructed_winds = .true. -endif -call nc_close_file(ncid) - -end subroutine set_lbc_variables - - -!------------------------------------------------------------------- -! modify what static_init_model does. this *must* be called before -! calling static_init_model(). -! set a logical add_u_to_state_list that forces u to be in state - -subroutine force_u_into_state() - -add_u_to_state_list = .true. - -end subroutine force_u_into_state - - -!------------------------------------------------------------------- - -!>@todo FIXME: -!> i believe no one is calling this routine anymore. we should remove -!> it from the public list and see what breaks. if nothing, remove it. -!> -!> these need to be replaced by calls to the state structure. -!> at the moment they are only used by the postprocessing program that -!> moves the wind increments from cell centers to edges - it's not needed -!> by anything in filter. it's a holdover from code before the state structure -!> was a general facility. - -subroutine analysis_file_to_statevector(filename, state_vector, model_time) - -! Reads the current time and state variables from a mpas analysis -! file and packs them into a dart state vector. - -character(len=*), intent(in) :: filename -real(r8), intent(inout) :: state_vector(:) -type(time_type), intent(out) :: model_time - -! temp space to hold data while we are reading it -integer :: i, ivar -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, mystart, mycount -character(len=NF90_MAX_NAME) :: varname -integer :: VarID, ncNdims, dimlen -integer :: ncid, TimeDimID, TimeDimLength -character(len=256) :: myerrorstring -character(len=*), parameter :: routine = 'analysis_file_to_statevector' - -if ( .not. module_initialized ) call static_init_model - -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,'analysis_file_to_statevector',string1,source,revision,revdate) -endif - -ncid = nc_open_file_readonly(filename, routine) - -model_time = get_analysis_time(ncid, filename) - -! Start counting and filling the state vector one item at a time, -! repacking the Nd arrays into a single 1d list of numbers. - -! The DART prognostic variables are only defined for a single time. -! We already checked the assumption that variables are xy2d or xyz3d ... -! IF the netCDF variable has a TIME dimension, it must be the last dimension, -! and we need to read the LAST timestep and effectively squeeze out the -! singleton dimension when we stuff it into the DART state vector. - -TimeDimID = FindTimeDimension( ncid ) - -if ( TimeDimID > 0 ) then - call nc_check(nf90_inquire_dimension(ncid, TimeDimID, len=TimeDimLength), & - 'analysis_file_to_statevector', 'inquire timedimlength '//trim(filename)) -else - TimeDimLength = 0 -endif - -do ivar=1, nfields - - varname = progvar(ivar)%varname - myerrorstring = trim(filename)//' '//trim(varname) - - ! determine the shape of the netCDF variable - - call nc_check(nf90_inq_varid(ncid, varname, VarID), & - 'analysis_file_to_statevector', 'inq_varid '//trim(myerrorstring)) - - call nc_check(nf90_inquire_variable(ncid,VarID,dimids=dimIDs,ndims=ncNdims), & - 'analysis_file_to_statevector', 'inquire '//trim(myerrorstring)) - - mystart = 1 ! These are arrays, actually. - mycount = 1 - - ! Only checking the shape of the variable - sans TIME - DimCheck : do i = 1,progvar(ivar)%numdims - - write(string1,'(''inquire dimension'',i2,A)') i,trim(myerrorstring) - call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), len=dimlen), & - 'analysis_file_to_statevector', string1) - - mycount(i) = dimlen - - enddo DimCheck - - where(dimIDs == TimeDimID) mystart = TimeDimLength ! pick the latest time - where(dimIDs == TimeDimID) mycount = 1 ! only use one time - - if ((debug > 0) .and. do_output()) then - write(*,*)'analysis_file_to_statevector '//trim(varname)//' start = ',mystart(1:ncNdims) - write(*,*)'analysis_file_to_statevector '//trim(varname)//' count = ',mycount(1:ncNdims) - endif - - if (ncNdims == 1) then - - ! If the single dimension is TIME, we only need a scalar. - ! Pretty sure this cannot happen ... - allocate(data_1d_array(mycount(1))) - call nc_check(nf90_get_var(ncid, VarID, data_1d_array, & - start=mystart(1:ncNdims), count=mycount(1:ncNdims)), & - 'analysis_file_to_statevector', 'get_var '//trim(varname)) - - call prog_var_to_vector(data_1d_array, state_vector, ivar) - deallocate(data_1d_array) - - elseif (ncNdims == 2) then - - allocate(data_2d_array(mycount(1), mycount(2))) - call nc_check(nf90_get_var(ncid, VarID, data_2d_array, & - start=mystart(1:ncNdims), count=mycount(1:ncNdims)), & - 'analysis_file_to_statevector', 'get_var '//trim(varname)) - - call prog_var_to_vector(data_2d_array, state_vector, ivar) - deallocate(data_2d_array) - - elseif (ncNdims == 3) then - - allocate(data_3d_array(mycount(1), mycount(2), mycount(3))) - call nc_check(nf90_get_var(ncid, VarID, data_3d_array, & - start=mystart(1:ncNdims), count=mycount(1:ncNdims)), & - 'analysis_file_to_statevector', 'get_var '//trim(varname)) - - call prog_var_to_vector(data_3d_array, state_vector, ivar) - deallocate(data_3d_array) - - else - write(string1, *) 'no support for data array of dimension ', ncNdims - call error_handler(E_ERR,'analysis_file_to_statevector', string1, & - source,revision,revdate) - endif - -enddo - -call nc_close_file(ncid, routine) - -end subroutine analysis_file_to_statevector - - -!------------------------------------------------------------------- - -!> @todo FIXME this routine is a holdover from the days before -!> we had a state structure module. it should be replaced by -!> something like write_state and other calls to the state structure -!> module. -!> -!> but there is one twist here. this code is only being called -!> by update_mpas_states, which does a simple ncks copy from the -!> output of filter (netcdf with only the state vector fields in it) -!> to a full mpas restart file with grid info, other fields, etc. -!> -!> however it also averages the cell center winds and projects them -!> onto the cell edges and updates 'u' - which is usually not in the -!> model state, so it isn't in the input file, only the output file. -!> -!> the normal i/o routines wouldn't write it. they would write -!> everything else, so read_state()/write_state() plus one additional -!> 'fix_u' routine would seem to suffice. TBD. -!> - -subroutine statevector_to_analysis_file(state_vector, ncid, filename) - -! Writes the current time and state variables from a dart state -! vector (1d array) into a mpas netcdf analysis file. - -real(r8), intent(in) :: state_vector(:) -integer, intent(in) :: ncid -character(len=*), intent(in) :: filename - -! temp space to hold data while we are writing it -integer :: i, ivar -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, mystart, mycount -character(len=NF90_MAX_NAME) :: varname -integer :: VarID, ncNdims, dimlen -integer :: TimeDimID, TimeDimLength -logical :: done_winds -type(time_type) :: model_time - -if ( .not. module_initialized ) call static_init_model - -TimeDimID = FindTimeDimension( ncid ) - -if ( TimeDimID > 0 ) then - call nc_check(nf90_inquire_dimension(ncid, TimeDimID, len=TimeDimLength), & - 'statevector_to_analysis_file', 'inquire timedimlength '//trim(filename)) -else - TimeDimLength = 0 -endif - -done_winds = .false. -PROGVARLOOP : do ivar=1, nfields - - varname = progvar(ivar)%varname - string2 = trim(filename)//' '//trim(varname) - - if (( varname == 'uReconstructZonal' .or. & - varname == 'uReconstructMeridional' ) .and. update_u_from_reconstruct ) then - if (done_winds) cycle PROGVARLOOP - - ! this routine updates the edge winds from both the zonal and meridional - ! fields, so only call it once. - call update_wind_components(ncid, filename, state_vector, use_increments_for_u_update) - done_winds = .true. - cycle PROGVARLOOP - endif - if ( varname == 'u' .and. update_u_from_reconstruct ) then - write(string1, *) 'skipping update of edge normal winds (u) because' - write(string2, *) 'update_u_from_reconstruct is True' - call error_handler(E_MSG,'statevector_to_analysis_file',string1,& - source,revision,revdate, text2=string2) - cycle PROGVARLOOP - endif - - ! Ensure netCDF variable is conformable with progvar quantity. - ! The TIME and Copy dimensions are intentionally not queried - ! by looping over the dimensions stored in the progvar type. - - call nc_check(nf90_inq_varid(ncid, varname, VarID), & - 'statevector_to_analysis_file', 'inq_varid '//trim(string2)) - - call nc_check(nf90_inquire_variable(ncid,VarID,dimids=dimIDs,ndims=ncNdims), & - 'statevector_to_analysis_file', 'inquire '//trim(string2)) - - mystart = 1 ! These are arrays, actually. - mycount = 1 - DimCheck : do i = 1,progvar(ivar)%numdims - - write(string1,'(''inquire dimension'',i2,A)') i,trim(string2) - call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), len=dimlen), & - 'statevector_to_analysis_file', string1) - - mycount(i) = dimlen - - enddo DimCheck - - - where(dimIDs == TimeDimID) mystart = TimeDimLength - where(dimIDs == TimeDimID) mycount = 1 ! only the latest one - - if ((debug > 0) .and. do_output()) then - write(*,*)'statevector_to_analysis_file '//trim(varname)//' start is ',mystart(1:ncNdims) - write(*,*)'statevector_to_analysis_file '//trim(varname)//' count is ',mycount(1:ncNdims) - endif - -!> @todo FIXME the clamping can be done on the 1d array by getting -!> the start/end index from the state vector. you can also call nf90_put_var() -!> with a 1d conformable array without allocating, copying, writing, and freeing -!> extra space. or call write_state() which does all this for us. - - if (progvar(ivar)%numdims == 1) then - allocate(data_1d_array(mycount(1))) - call vector_to_prog_var(state_vector, ivar, data_1d_array) - - ! did the user specify lower and/or upper bounds for this variable? - ! if so, follow the instructions to either fail on out-of-range values, - ! or set out-of-range values to the given min or max vals - if ( progvar(ivar)%clamping ) then - call do_clamping(progvar(ivar)%out_of_range_fail, progvar(ivar)%range, & - progvar(ivar)%numdims, varname, array_1d = data_1d_array) - endif - - call nc_check(nf90_put_var(ncid, VarID, data_1d_array, & - start=mystart(1:ncNdims), count=mycount(1:ncNdims)), & - 'statevector_to_analysis_file', 'put_var '//trim(varname)) - deallocate(data_1d_array) - - elseif (progvar(ivar)%numdims == 2) then - - allocate(data_2d_array(mycount(1), mycount(2))) - call vector_to_prog_var(state_vector, ivar, data_2d_array) - - ! did the user specify lower and/or upper bounds for this variable? - ! if so, follow the instructions to either fail on out-of-range values, - ! or set out-of-range values to the given min or max vals - if ( progvar(ivar)%clamping ) then - call do_clamping(progvar(ivar)%out_of_range_fail, progvar(ivar)%range, & - progvar(ivar)%numdims, varname, array_2d = data_2d_array) - endif - - call nc_check(nf90_put_var(ncid, VarID, data_2d_array, & - start=mystart(1:ncNdims), count=mycount(1:ncNdims)), & - 'statevector_to_analysis_file', 'put_var '//trim(varname)) - deallocate(data_2d_array) - - elseif (progvar(ivar)%numdims == 3) then - - allocate(data_3d_array(mycount(1), mycount(2), mycount(3))) - call vector_to_prog_var(state_vector, ivar, data_3d_array) - - ! did the user specify lower and/or upper bounds for this variable? - ! if so, follow the instructions to either fail on out-of-range values, - ! or set out-of-range values to the given min or max vals - if ( progvar(ivar)%clamping ) then - call do_clamping(progvar(ivar)%out_of_range_fail, progvar(ivar)%range, & - progvar(ivar)%numdims, varname, array_3d = data_3d_array) - endif - - call nc_check(nf90_put_var(ncid, VarID, data_3d_array, & - start=mystart(1:ncNdims), count=mycount(1:ncNdims)), & - 'statevector_to_analysis_file', 'put_var '//trim(varname)) - deallocate(data_3d_array) - - else - write(string1, *) 'no support for data array of dimension ', ncNdims - call error_handler(E_ERR,'statevector_to_analysis_file', string1, & - source,revision,revdate) endif - -enddo PROGVARLOOP - - -end subroutine statevector_to_analysis_file - -!------------------------------------------------------------------ - -!> regional mpas only: -!> update the boundary fields based on the analysis file values -!> in the boundary region, which were updated by blending the analysis -!> and the original (e.g., prior) boundary values. -!> There are options to update edge winds directly (by blending 'u'), -!> or to update reconstructed winds at cell centers first then project them onto edges. -!> When the edge wind is updated, it can be replaced by the blended value, -!> or modified with the increments (from cell-center winds). -!> state_vector - read from both ncid_a and ncid_b. -!> ncid_b = lbc to be overwritten with blended fields in the boundary zone. -!> ncid_a = analysis - either init or restart. - -subroutine statevector_to_boundary_file(state_vector, ncid_b, ncid_a, & - lbc_update_from_reconstructed_winds, lbc_update_winds_from_increments, idebug) - -real(r8), intent(inout) :: state_vector(:) -integer, intent(in) :: ncid_b, ncid_a -logical, intent(in) :: lbc_update_from_reconstructed_winds -logical, intent(in) :: lbc_update_winds_from_increments -integer, intent(in) :: idebug - -integer :: i, a_ivar, b_ivar, ivar, ivar_u, avar_u -integer(i8) :: a_index, b_index, l, sb_index, eb_index -integer :: cellid, vert_level, ndims, nvars, dims(3), col -integer :: adims, dima(3), iup ! HA -integer :: edgeid -real(r8) :: weight -real(r8), allocatable :: lbc_u(:,:), lbc_ucell(:,:), lbc_vcell(:,:) -real(r8), allocatable :: delta_u(:,:), old_lbc_ucell(:,:), old_lbc_vcell(:,:) -real(r8), allocatable :: inc_lbc_ucell(:,:), inc_lbc_vcell(:,:) - -character(len=NF90_MAX_NAME) :: avarname, bvarname -character(len=*), parameter :: routine = 'statevector_to_boundary_file' - -nvars = get_num_variables(lbc_domid) - -write(string1, *) 'lbc_update_from_reconstructed_winds = ', lbc_update_from_reconstructed_winds -write(string2, *) 'lbc_update_winds_from_increments = ',lbc_update_winds_from_increments -call error_handler(E_MSG,'statevector_to_boundary_file',string1,& - source,revision,revdate, text2=string2) - -! save a copy of the reconstructed cell winds in separate arrays -! if we are doing an incremental update of the edge normal winds. -if (lbc_update_winds_from_increments) then - if (.not. lbc_file_has_reconstructed_winds) then - write(string1, *) 'Cannot update edge winds from increments because the boundary file does not contain the reconstructed winds (lbc_ur, lbc_vr)' - write(string2, *) 'lbc_update_winds_from_increments should be .false.' - call error_handler(E_MSG,'statevector_to_boundary_file',string1,& - source,revision,revdate, text2=string2) - endif - - allocate(old_lbc_ucell(nVertLevels, nCells)) - allocate(old_lbc_vcell(nVertLevels, nCells)) - allocate( delta_u(nVertLevels, nEdges)) - - ivar = get_varid_from_varname(lbc_domid, 'lbc_ur') - call bdy_vector_to_prog_var(state_vector, ivar, old_lbc_ucell) - - ivar = get_varid_from_varname(lbc_domid, 'lbc_vr') - call bdy_vector_to_prog_var(state_vector, ivar, old_lbc_vcell) endif -! for each cell in the grid, find the analysis in the -! boundary region and blend them with prior lbc values. -CELLS: do cellid = 1, nCells - - ! Soyoung: We blend the analysis in the boundary zone only. - if (.not. on_boundary_cell(cellid)) cycle CELLS - - ! 1.0 is interior, 0.0 is exterior boundary - weight = get_analysis_weight(cellid) - - ! do all variables associated with this cellid. - - VARLOOP: do b_ivar = 1, nvars - - bvarname = get_variable_name(lbc_domid, b_ivar) - if (bvarname(1:4) /= 'lbc_') then - write(string1, *) 'skipping update of boundary variable ', trim(bvarname) - write(string2, *) 'because the name does not start with "lbc"' - call error_handler(E_MSG,'statevector_to_boundary_file',string1,& - source,revision,revdate, text2=string2) - cycle VARLOOP - endif - - ! skip edge normal 'U' winds here - they will - ! be handled in a separate code section below. - if (bvarname == 'lbc_u') cycle VARLOOP - - ! get corresponding field in analysis domain - avarname = trim(bvarname(5:)) - - ! reconstructed cell-center winds have different names in the lbc file. - if (bvarname == 'lbc_ur') avarname = 'uReconstructZonal' - if (bvarname == 'lbc_vr') avarname = 'uReconstructMeridional' - - a_ivar = get_varid_from_varname(anl_domid, avarname) - - !if(do_output().and.idebug > 4) print*, 'statevector_to_boundary_file: ', & - ! trim(bvarname), b_ivar, trim(avarname), a_ivar - - call find_mpas_dims(lbc_domid, b_ivar, ndims, dims) - - ! HA: double-check if dimensions are the same between lbc_domid and anl_domid. - call find_mpas_dims(anl_domid, a_ivar, adims, dima) - if(dims(1) /= dima(1) .or. dims(2) /= dima(2)) then - write(string1, *) 'Dimension mismatches:',dims,' vs.',dima - call error_handler(E_ERR,'statevector_to_boundary_file',string1,& - source,revision,revdate) - exit - endif - - ! loop over vert_levels. - THISCOL: do col=1, dims(1) - a_index = get_dart_vector_index(col, cellid, 1, anl_domid, a_ivar) - b_index = get_dart_vector_index(col, cellid, 1, lbc_domid, b_ivar) - - ! compute (1-w)*x_lbc + w*x_anl - state_vector(a_index) = (1.0_r8 - weight) * state_vector(b_index) + & - weight * state_vector(a_index) - - enddo THISCOL - - enddo VARLOOP - -enddo CELLS - -!> here is where we fix up the U edge normal winds -!> two options - do them directly, or compute increments from the -!> reconstructed winds and update from them. - -if (.not. lbc_update_from_reconstructed_winds) then - - ! this is the prior u, not updated yet - a_ivar = get_varid_from_varname(anl_domid, 'u') ! analysis edge winds - b_ivar = get_varid_from_varname(lbc_domid, 'lbc_u') ! prior edge winds in the lbc file - - call find_mpas_dims(lbc_domid, b_ivar, ndims, dims) - - ! for each edge in the grid, find the ones which are in the - ! boundary region and blend their values. - EDGES: do edgeid = 1, nEdges - - if (.not. on_boundary_edge(edgeid)) cycle EDGES - - ! 1.0 is interior, 0.0 is exterior boundary - weight = get_analysis_weight(edgeid,.false.) - - ! loop over vert_levels. - THATCOL: do col=1, dims(1) - a_index = get_dart_vector_index(col, edgeid, 1, anl_domid, a_ivar) - b_index = get_dart_vector_index(col, edgeid, 1, lbc_domid, b_ivar) - - ! compute (1-w)*x_lbc + w*x_anl - state_vector(a_index) = (1.0_r8 - weight) * state_vector(b_index) + & - weight * state_vector(a_index) - enddo THATCOL - - enddo EDGES - -else ! do the increment process +! Loop over potentially close subset of obs priors or state variables +call loc_get_close_state(gc, base_loc, base_type, locs, loc_qtys, loc_indx, num_close, close_ind) - ! We only blended cell-center fields (e.g. looping over all the variables in the CELLS loop above, but not over 'u' in nEdges). - ! Now we compute diffs (or increments) between the blended ur (vr) and the prior lbc_ur (vr). - - allocate( lbc_u(nVertLevels, nEdges)) - allocate( lbc_ucell(nVertLevels, nCells)) - allocate( lbc_vcell(nVertLevels, nCells)) - - ! these analyses have been blended in the boundary zone already. - ivar = get_varid_from_varname(anl_domid, 'uReconstructZonal') - call vector_to_prog_var(state_vector, ivar, lbc_ucell) - - ivar = get_varid_from_varname(anl_domid, 'uReconstructMeridional') - call vector_to_prog_var(state_vector, ivar, lbc_vcell) - - ! this is the analysis u, not blended in the boundary zone yet. - avar_u = get_varid_from_varname(anl_domid, 'u') - call vector_to_prog_var(state_vector, avar_u, lbc_u) - - if (idebug > 9 .and. do_output()) print *, 'MIN/MAX lbc_u before update:',MINVAL(lbc_u),MAXVAL(lbc_u) - - if (lbc_update_winds_from_increments) then - - ! project analysis increments at cell centers onto the edges. - - allocate(inc_lbc_ucell(nVertLevels, nCells)) - allocate(inc_lbc_vcell(nVertLevels, nCells)) - - inc_lbc_ucell = lbc_ucell - old_lbc_ucell - inc_lbc_vcell = lbc_vcell - old_lbc_vcell - - call uv_cell_to_edges(inc_lbc_ucell, inc_lbc_vcell, delta_u) - - ! Soyoung: Add the blended u increments back to lbc_u in the boundary zone. - ! We should not change the analysis u in the interior domain, but - ! We should also check bdyMaskCell for the two adjacent cells as - ! bdyMaskEdge is assigned with the lower mask value between the two - ! cells. - ! Ex) An edge between cell1 (w/ bdyMaskCell = 0) and cell2 (w/ bdyMaskCell = 1) - ! has bdyMaskEdge = 0. In this case, even if bdyMaskEdge of the edge is zero, - ! cell2 has been updated in the CELLS loop above, so the edge has to be updated. - - iup = 0 - IEDGE: do edgeid = 1, nEdges - - if (.not. on_boundary_edge(edgeid) .and. & - .not. on_boundary_cell(cellsOnEdge(1,edgeid)) .and. & - .not. on_boundary_cell(cellsOnEdge(2,edgeid)) ) cycle IEDGE - - lbc_u(:,edgeid) = lbc_u(:,edgeid) + delta_u(:,edgeid) - - if (idebug > 9 .and. .not.on_boundary_edge(edgeid)) & - print*, iup, edgeid, delta_u(1,edgeid) - - enddo IEDGE - - if (idebug > 9 .and. do_output()) print*, 'MIN/MAX delta_u: ',MINVAL(delta_u),MAXVAL(delta_u) - - deallocate(old_lbc_ucell, old_lbc_vcell, delta_u) - - else - - ! just replace, no increments - ! project the diffs (or increments) onto the edges by calling uv_cell_to_edges. - call uv_cell_to_edges(lbc_ucell, lbc_vcell, lbc_u, .true.) - - endif - if (idebug > 9 .and. do_output()) print *, 'MIN/MAX lbc_u after update:',MINVAL(lbc_u),MAXVAL(lbc_u) - - ! put lbc_u array data back into the state_vector - - sb_index = get_index_start(anl_domid, avar_u) - eb_index = get_index_end (anl_domid, avar_u) - state_vector(sb_index:eb_index) = reshape(lbc_u, (/eb_index-sb_index+1/) ) - - deallocate(lbc_u, lbc_ucell, lbc_vcell) - -endif ! U updates - - -! for each boundary variable, write it to the output file. -VARLOOP2: do b_ivar = 1, nvars - - bvarname = get_variable_name(lbc_domid, b_ivar) - if (bvarname(1:4) /= 'lbc_') then - write(string1, *) 'skipping update of boundary variable ', trim(bvarname) - write(string2, *) 'because the name does not start with "lbc"' - call error_handler(E_MSG,'statevector_to_boundary_file',string1,& - source,revision,revdate, text2=string2) - cycle VARLOOP2 - endif - - ! by default strip off 'lbc_' from boundary file name - ! and update that field name for the analysis file, - ! unless it doesn't follow the pattern - avarname = trim(bvarname(5:)) - - if (bvarname == 'lbc_ur') avarname = 'uReconstructZonal' - if (bvarname == 'lbc_vr') avarname = 'uReconstructMeridional' - - ! Soyoung - we blended the analysis vector in the boundary zone. - a_ivar = get_varid_from_varname(anl_domid, avarname) +do k = 1, num_close - ! nsc - it's possible we could remove this line and the - ! reshape()s from put_variable lines below. - call find_mpas_dims(lbc_domid, b_ivar, ndims, dims) + t_ind = close_ind(k) + local_obs_loc = locs(t_ind) + local_obs_which = nint(query_location(local_obs_loc)) - sb_index = get_index_start(anl_domid, a_ivar) - eb_index = get_index_end (anl_domid, a_ivar) - - if (idebug > 9 .and. do_output()) print *, 'updating ', trim(bvarname), ' min, max:',& - minval(state_vector(sb_index:eb_index)), maxval(state_vector(sb_index:eb_index)) - !, ' with length ', eb_index - sb_index + 1 - - ! Soyoung - Now, the lbc file has the analysis in the interior and the blended values - ! in the boundary zone. In other words, the lbc file is updated not only in - ! the boundary, but over the entire domain (although it is used only for the - ! boundary zone). - call nc_put_variable(ncid_b, bvarname, reshape(state_vector(sb_index:eb_index), dims), routine) -! call nc_put_variable(ncid_a, avarname, reshape(state_vector(sb_index:eb_index), dims), routine) - -enddo VARLOOP2 - -end subroutine statevector_to_boundary_file - -!------------------------------------------------------------------ - -subroutine do_clamping(out_of_range_fail, range, dimsize, varname, array_1d, array_2d, array_3d) - logical, intent(in) :: out_of_range_fail - real(r8), intent(in) :: range(2) - integer, intent(in) :: dimsize - character(len=*), intent(in) :: varname - real(r8),optional,intent(inout) :: array_1d(:), array_2d(:,:), array_3d(:,:,:) - -! FIXME: This should be replaced by the clamp_variable in direct_netcdf_mpi. -! clamp_variable currently works on one dimensional arrays for variables. This -! only becomes an issue in statevector_to_analysis_file which requires the dimension -! of the variable to output to a netcdf file, however, this information is -! already stored in the state_structure. - -! for a given directive and range, do the data clamping for the given -! input array. only one of the optional array args should be specified - the -! one which matches the given dimsize. this still has replicated sections for -! each possible dimensionality (which so far is only 1 to 3 - add 4-7 only -! if needed) but at least it is isolated to this subroutine. - -! these sections should all be identical except for the array_XX specified. -! if anyone can figure out a way to defeat fortran's strong typing for arrays -! so we don't have to replicate each of these sections, i'll buy you a cookie. -! (sorry, you can't suggest using the preprocessor, which is the obvious -! solution. up to now we have avoided any preprocessed code in the entire -! system. if we cave at some future point this routine is a prime candidate -! to autogenerate.) - -if (dimsize == 1) then - if (.not. present(array_1d)) then - call error_handler(E_ERR, 'do_clamping', 'Internal error. Should not happen', & - source,revision,revdate, text2='array_1d not present for 1d case') - endif - - ! is lower bound set - if ( range(1) /= missing_r8 ) then - - if ( out_of_range_fail ) then - if ( minval(array_1d) < range(1) ) then - write(string1, *) 'min data val = ', minval(array_1d), & - 'min data bounds = ', range(1) - call error_handler(E_ERR, 'statevector_to_analysis_file', & - 'Variable '//trim(varname)//' failed lower bounds check.', & - source,revision,revdate) - endif + ! Convert local_obs vertical coordinate to requested vertical coordinate if necessary. + if ((base_which /= VERTISUNDEF) .and. (vertical_localization_on())) then + if (local_obs_which /= vert_localization_coord .and. local_obs_which /= VERTISUNDEF) then + location_arr(1) = local_obs_loc + call convert_vert_state(state_handle, 1, location_arr, loc_qtys(t_ind), & + loc_indx(t_ind), vert_localization_coord, istatus) + locs(t_ind) = location_arr(1) else - where ( array_1d < range(1) ) array_1d = range(1) + istatus = 0 endif - - endif ! min range set - - ! is upper bound set - if ( range(2) /= missing_r8 ) then - - if ( out_of_range_fail ) then - if ( maxval(array_1d) > range(2) ) then - write(string1, *) 'max data val = ', maxval(array_1d), & - 'max data bounds = ', range(2) - call error_handler(E_ERR, 'statevector_to_analysis_file', & - 'Variable '//trim(varname)//' failed upper bounds check.', & - source,revision,revdate, text2=string1) - endif - else - where ( array_1d > range(2) ) array_1d = range(2) - endif - - endif ! max range set - - write(string1, '(A,A32,2F16.7)') 'BOUND min/max ', trim(varname), & - minval(array_1d), maxval(array_1d) - call error_handler(E_MSG, '', string1, source,revision,revdate) - -else if (dimsize == 2) then - if (.not. present(array_2d)) then - call error_handler(E_ERR, 'do_clamping', 'Internal error. Should not happen', & - source,revision,revdate, text2='array_2d not present for 2d case') + else + istatus = 0 endif - ! is lower bound set - if ( range(1) /= missing_r8 ) then - - if ( out_of_range_fail ) then - if ( minval(array_2d) < range(1) ) then - write(string1, *) 'min data val = ', minval(array_2d), & - 'min data bounds = ', range(1) - call error_handler(E_ERR, 'statevector_to_analysis_file', & - 'Variable '//trim(varname)//' failed lower bounds check.', & - source,revision,revdate) - endif + if (present(dist)) then + if (istatus(1) == 0) then + dist(k) = get_dist(base_loc, locs(t_ind), base_type, loc_qtys(t_ind)) else - where ( array_2d < range(1) ) array_2d = range(1) + dist(k) = 1.0e9_r8 endif - - endif ! min range set - - ! is upper bound set - if ( range(2) /= missing_r8 ) then - - if ( out_of_range_fail ) then - if ( maxval(array_2d) > range(2) ) then - write(string1, *) 'max data val = ', maxval(array_2d), & - 'max data bounds = ', range(2) - call error_handler(E_ERR, 'statevector_to_analysis_file', & - 'Variable '//trim(varname)//' failed upper bounds check.', & - source,revision,revdate, text2=string1) - endif - else - where ( array_2d > range(2) ) array_2d = range(2) - endif - - endif ! max range set - - write(string1, '(A,A32,2F16.7)') 'BOUND min/max ', trim(varname), & - minval(array_2d), maxval(array_2d) - call error_handler(E_MSG, '', string1, source,revision,revdate) - -else if (dimsize == 3) then - if (.not. present(array_3d)) then - call error_handler(E_ERR, 'do_clamping', 'Internal error. Should not happen', & - source,revision,revdate, text2='array_3d not present for 3d case') endif - ! is lower bound set - if ( range(1) /= missing_r8 ) then - - if ( out_of_range_fail ) then - if ( minval(array_3d) < range(1) ) then - write(string1, *) 'min data val = ', minval(array_3d), & - 'min data bounds = ', range(1) - call error_handler(E_ERR, 'statevector_to_analysis_file', & - 'Variable '//trim(varname)//' failed lower bounds check.', & - source,revision,revdate) - endif - else - where ( array_3d < range(1) ) array_3d = range(1) - endif - - endif ! min range set - - ! is upper bound set - if ( range(2) /= missing_r8 ) then - - if ( out_of_range_fail ) then - if ( maxval(array_3d) > range(2) ) then - write(string1, *) 'max data val = ', maxval(array_3d), & - 'max data bounds = ', range(2) - call error_handler(E_ERR, 'statevector_to_analysis_file', & - 'Variable '//trim(varname)//' failed upper bounds check.', & - source,revision,revdate, text2=string1) - endif - else - where ( array_3d > range(2) ) array_3d = range(2) - endif - - endif ! max range set - - write(string1, '(A,A32,2F16.7)') 'BOUND min/max ', trim(varname), & - minval(array_3d), maxval(array_3d) - call error_handler(E_MSG, '', string1, source,revision,revdate) + enddo -else - write(string1, *) 'dimsize of ', dimsize, ' found where only 1-3 expected' - call error_handler(E_MSG, 'do_clamping', 'Internal error, should not happen', & - source,revision,revdate, text2=string1) -endif ! dimsize +end subroutine get_close_state -end subroutine do_clamping -!------------------------------------------------------------------ +!================================================================== +! The (model-specific) additional public interfaces come next +! (these are not required by dart but are used by other programs) +!================================================================== function get_analysis_time_ncid( ncid, filename ) @@ -2951,7 +1380,7 @@ function get_analysis_time_ncid( ncid, filename ) if (numdims /= 2) then write(string1,*) 'xtime variable has unknown shape in ', trim(filename) - call error_handler(E_ERR,'get_analysis_time',string1,source,revision,revdate) + call error_handler(E_ERR,'get_analysis_time',string1,source) endif call nc_check( nf90_inquire_dimension(ncid, dimIDs(1), len=idims(1)), & @@ -2962,7 +1391,7 @@ function get_analysis_time_ncid( ncid, filename ) if (idims(2) /= 1) then write(string1,*) 'multiple timesteps (',idims(2),') in file ', trim(filename) write(string2,*) 'We are using the LAST one, presumably, the LATEST timestep.' - call error_handler(E_MSG,'get_analysis_time',string1,source,revision,revdate,text2=string2) + call error_handler(E_MSG,'get_analysis_time',string1,source,text2=string2) endif ! Get the highest ranking time ... the last one, basically. @@ -2972,11 +1401,6 @@ function get_analysis_time_ncid( ncid, filename ) get_analysis_time_ncid = string_to_time(timestring) -if ((debug > 0) .and. do_output()) then - call print_date(get_analysis_time_ncid, 'get_analysis_time:model date') - call print_time(get_analysis_time_ncid, 'get_analysis_time:model time') -endif - end function get_analysis_time_ncid @@ -3000,14 +1424,14 @@ function get_analysis_time_fname(filename) ! string from the filename itself. if ( .not. file_exist(filename) ) then write(string1,*) 'file ', trim(filename),' does not exist.' - call error_handler(E_ERR,'get_analysis_time',string1,source,revision,revdate) + call error_handler(E_ERR,'get_analysis_time',string1,source) endif ! find the first digit and use that as the start of the string conversion i = scan(filename, "0123456789") if (i <= 0) then write(string1,*) 'cannot find time string in name ', trim(filename) - call error_handler(E_ERR,'get_analysis_time',string1,source,revision,revdate) + call error_handler(E_ERR,'get_analysis_time',string1,source) endif get_analysis_time_fname = string_to_time(filename(i:i+TIMELEN-1)) @@ -3015,44 +1439,15 @@ function get_analysis_time_fname(filename) end function get_analysis_time_fname -!------------------------------------------------------------------ - -subroutine write_model_time_file(time_filename, model_time, adv_to_time) - character(len=*), intent(in) :: time_filename - type(time_type), intent(in) :: model_time - type(time_type), intent(in), optional :: adv_to_time - -integer :: iunit -character(len=TIMELEN) :: timestring -type(time_type) :: deltatime - -iunit = open_file(time_filename, action='write') - -timestring = time_to_string(model_time) -write(iunit, '(A)') timestring - -if (present(adv_to_time)) then - timestring = time_to_string(adv_to_time) - write(iunit, '(A)') timestring - - deltatime = adv_to_time - model_time - timestring = time_to_string(deltatime, interval=.true.) - write(iunit, '(A)') timestring -endif - -call close_file(iunit) - -end subroutine write_model_time_file - !----------------------------------------------------------------------- -subroutine write_model_time_restart(ncid, dart_time) +subroutine write_model_time(ncid, dart_time) integer, intent(in) :: ncid !< netcdf file handle type(time_type), intent(in) :: dart_time integer :: year, month, day, hour, minute, second character(len=64) :: timestring -character(len=*), parameter :: routine = 'write_model_time_restart' +character(len=*), parameter :: routine = 'write_model_time' timestring = '' call get_date(dart_time, year, month, day, hour, minute, second) @@ -3081,15 +1476,12 @@ subroutine write_model_time_restart(ncid, dart_time) call nc_put_variable(ncid, 'xtime', timestring, routine) -end subroutine write_model_time_restart +end subroutine write_model_time !------------------------------------------------------------------ subroutine get_grid_dims(Cells, Vertices, Edges, VertLevels, VertexDeg, SoilLevels) -! public routine for returning the counts of various things in the grid -! - integer, intent(out) :: Cells ! Total number of cells making up the grid integer, intent(out) :: Vertices ! Unique points in grid which are corners of cells integer, intent(out) :: Edges ! Straight lines between vertices making up cells @@ -3190,49 +1582,6 @@ end subroutine get_bdy_mask ! The (model-specific) private interfaces come last !================================================================== - -!------------------------------------------------------------------ -!> convert time type into a character string with the -!> format of YYYY-MM-DD_hh:mm:ss - -function time_to_string(t, interval) - - type(time_type), intent(in) :: t - logical, intent(in), optional :: interval -character(len=TIMELEN) :: time_to_string - -integer :: iyear, imonth, iday, ihour, imin, isec -integer :: ndays, nsecs -logical :: dointerval - - dointerval = .false. -if (present(interval)) dointerval = interval - -! for interval output, output the number of days, then hours, mins, secs -! for date output, use the calendar routine to get the year/month/day hour:min:sec -if (dointerval) then - call get_time(t, nsecs, ndays) - if (ndays > 99) then - write(string1, *) 'interval number of days is ', ndays - call error_handler(E_ERR,'time_to_string', 'interval days cannot be > 99', & - source, revision, revdate, text2=string1) - endif - ihour = nsecs / 3600 - nsecs = nsecs - (ihour * 3600) - imin = nsecs / 60 - nsecs = nsecs - (imin * 60) - isec = nsecs - write(time_to_string, '(I2.2,3(A1,I2.2))') & - ndays, '_', ihour, ':', imin, ':', isec -else - call get_date(t, iyear, imonth, iday, ihour, imin, isec) - write(time_to_string, '(I4.4,5(A1,I2.2))') & - iyear, '-', imonth, '-', iday, '_', ihour, ':', imin, ':', isec -endif - -end function time_to_string - - !------------------------------------------------------------------ function string_to_time(s) @@ -3256,51 +1605,14 @@ end function string_to_time function set_model_time_step() -! the static_init_model ensures that the model namelists are read. - type(time_type) :: set_model_time_step if ( .not. module_initialized ) call static_init_model -! these are from the namelist set_model_time_step = set_time(assimilation_period_seconds, assimilation_period_days) end function set_model_time_step - -!------------------------------------------------------------------ -!> Read the grid dimensions from the MPAS netcdf file. -!> - -subroutine read_grid_dims(ncid) -integer, intent(in) :: ncid - -character(len=*), parameter :: routine = 'read_grid_dims' - -nCells = nc_get_dimension_size(ncid, 'nCells', routine) -nVertices = nc_get_dimension_size(ncid, 'nVertices', routine) -nEdges = nc_get_dimension_size(ncid, 'nEdges', routine) -maxEdges = nc_get_dimension_size(ncid, 'maxEdges', routine) -nVertLevels = nc_get_dimension_size(ncid, 'nVertLevels', routine) -nVertLevelsP1 = nc_get_dimension_size(ncid, 'nVertLevelsP1', routine) -vertexDegree = nc_get_dimension_size(ncid, 'vertexDegree', routine) -nSoilLevels = nc_get_dimension_size(ncid, 'nSoilLevels', routine) - -if (debug > 4 .and. do_output()) then - write(*,*) - write(*,*)'read_grid_dims: nCells is ', nCells - write(*,*)'read_grid_dims: nVertices is ', nVertices - write(*,*)'read_grid_dims: nEdges is ', nEdges - write(*,*)'read_grid_dims: maxEdges is ', maxEdges - write(*,*)'read_grid_dims: nVertLevels is ', nVertLevels - write(*,*)'read_grid_dims: nVertLevelsP1 is ', nVertLevelsP1 - write(*,*)'read_grid_dims: vertexDegree is ', vertexDegree - write(*,*)'read_grid_dims: nSoilLevels is ', nSoilLevels -endif - -end subroutine read_grid_dims - - !------------------------------------------------------------------ !> Read the grid values in from the MPAS netcdf file. !> @@ -3317,526 +1629,95 @@ subroutine get_grid(ncid) ! MPAS locations are in radians - at this point DART needs degrees. ! watch out for tiny rounding errors and clamp to exactly +/- 90 -latCell = latCell * rad2deg -lonCell = lonCell * rad2deg -where (latCell > 90.0_r8) latCell = 90.0_r8 -where (latCell < -90.0_r8) latCell = -90.0_r8 - -call nc_get_variable(ncid, 'dcEdge', dcEdge, routine) -call nc_get_variable(ncid, 'zgrid', zGridFace, routine) -call nc_get_variable(ncid, 'cellsOnVertex', cellsOnVertex, routine) -call nc_get_variable(ncid, 'xland', xland, routine) -call nc_get_variable(ncid, 'seaice', seaice, routine) -call nc_get_variable(ncid, 'skintemp', skintemp, routine) - -dxmax = maxval(dcEdge) ! max grid resolution in meters - -call nc_get_variable(ncid, 'edgeNormalVectors', edgeNormalVectors, routine) -call nc_get_variable(ncid, 'nEdgesOnCell', nEdgesOnCell, routine) -call nc_get_variable(ncid, 'edgesOnCell', edgesOnCell, routine) -call nc_get_variable(ncid, 'cellsOnEdge', cellsOnEdge, routine) -call nc_get_variable(ncid, 'verticesOnCell', verticesOnCell, routine) - -if(data_on_edges) then - call nc_get_variable(ncid, 'lonEdge', lonEdge, routine) - call nc_get_variable(ncid, 'latEdge', latEdge, routine) - - latEdge = latEdge * rad2deg - lonEdge = lonEdge * rad2deg - where (latEdge > 90.0_r8) latEdge = 90.0_r8 - where (latEdge < -90.0_r8) latEdge = -90.0_r8 - - call nc_get_variable(ncid, 'xEdge', xEdge, routine) - call nc_get_variable(ncid, 'yEdge', yEdge, routine) - call nc_get_variable(ncid, 'zEdge', zEdge, routine) -endif - -call nc_get_variable(ncid, 'xVertex', xVertex, routine) -call nc_get_variable(ncid, 'yVertex', yVertex, routine) -call nc_get_variable(ncid, 'zVertex', zVertex, routine) - -if (nc_variable_exists(ncid, 'maxLevelCell')) then - allocate(maxLevelCell(nCells)) - call nc_get_variable(ncid, 'maxLevelCell', maxlevelCell, routine) - all_levels_exist_everywhere = .false. -endif - - -! A little sanity check - -if ( debug > 9 .and. do_output() ) then - - write(*,*) - write(*,*)'latCell range ',minval(latCell), maxval(latCell) - write(*,*)'lonCell range ',minval(lonCell), maxval(lonCell) - write(*,*)'zgrid range ',minval(zGridFace), maxval(zGridFace) - write(*,*)'cellsOnVertex range ',minval(cellsOnVertex), maxval(cellsOnVertex) - write(*,*)'edgeNormalVectors range ',minval(edgeNormalVectors), maxval(edgeNormalVectors) - write(*,*)'nEdgesOnCell range ',minval(nEdgesOnCell), maxval(nEdgesOnCell) - write(*,*)'EdgesOnCell range ',minval(EdgesOnCell), maxval(EdgesOnCell) - write(*,*)'cellsOnEdge range ',minval(cellsOnEdge), maxval(cellsOnEdge) - write(*,*)'xland range ',minval(xland), maxval(xland) - write(*,*)'seaice range ',minval(seaice), maxval(seaice) ! for rttov - write(*,*)'skintemp range ',minval(skintemp), maxval(skintemp) ! for rttov - if(data_on_edges) then - write(*,*)'latEdge range ',minval(latEdge), maxval(latEdge) - write(*,*)'lonEdge range ',minval(lonEdge), maxval(lonEdge) - write(*,*)'xEdge range ',minval(xEdge), maxval(xEdge) - write(*,*)'yEdge range ',minval(yEdge), maxval(yEdge) - write(*,*)'zEdge range ',minval(zEdge), maxval(zEdge) - endif - write(*,*)'xVertex range ',minval(xVertex), maxval(xVertex) - write(*,*)'yVertex range ',minval(yVertex), maxval(yVertex) - write(*,*)'zVertex range ',minval(zVertex), maxval(zVertex) - write(*,*)'verticesOnCell range ',minval(verticesOnCell), maxval(verticesOnCell) - if (allocated(bdyMaskCell)) & - write(*,*)'bdyMaskCell range ',minval(bdyMaskCell), maxval(bdyMaskCell) - if (allocated(bdyMaskEdge)) & - write(*,*)'bdyMaskEdge range ',minval(bdyMaskEdge), maxval(bdyMaskEdge) - if (allocated(maxLevelCell)) & - write(*,*)'maxLevelCell range ',minval(maxLevelCell), maxval(maxLevelCell) - -endif - -end subroutine get_grid - - -!------------------------------------------------------------------ - -subroutine update_wind_components(ncid, filename, state_vector, use_increments_for_u_update) - - integer, intent(in) :: ncid - character(len=*), intent(in) :: filename - real(r8), intent(in) :: state_vector(:) - logical, intent(in) :: use_increments_for_u_update - -! the winds pose a special problem because the model uses the edge-normal component -! of the winds at the center of the cell edges at half levels in the vertical ('u'). -! the output files from the model can include interpolated Meridional and Zonal -! winds as prognostic fields, which are easier for us to use when computing -! forward operator values. but in the end we need to update 'u' in the output -! model file. - -! this routine is only called when 'u' is not being directly updated by the -! assimilation, and the updated cell center values need to be converted back -! to update 'u'. there are several choices for how to do this and most are -! controlled by namelist settings. - -! If 'use_increments_for_u_update' is .true.: -! Read in the previous reconstructed winds from the original mpas netcdf file -! and compute what increments (changes in values) were added by the assimilation. -! Read in the original edge normal wind 'u' field from that same mpas netcdf -! file and add the interpolated increments to compute the updated 'u' values. -! (note that we can't use the DART Prior_Diag.nc file to get the previous -! values if we're using Prior inflation, because the diagnostic values are -! written out after inflation is applied.) - -! If 'use_increments_for_u_update' is .false.: -! use the Zonal/Meridional cell center values directly. The edge normal winds -! are each directly between 2 cell centers, so average the components normal -! to the edge direction. don't read in the previous values at the cell centers -! or the edge normal winds. - -! there are several changes here from previous versions: -! 1. it requires both zonal and meridional fields to be there. it doesn't -! make sense to assimilate with only one component of the winds. -! 2. i removed the 'return if this has been called already' flag. -! this is a generic utility routine. if someone wrote a main program -! that wanted to cycle over multiple files in a loop, this would have -! only updated the first file and silently returned for all the rest. -! 3. i moved the netcdf code for 3 identical operations into a subroutine. -! the code is easier to read and it's less likely to make a mistake if -! the code needs to be changed (one place vs three). - -! space to hold existing data from the analysis file -real(r8), allocatable :: u(:,:) ! u(nVertLevels, nEdges) -real(r8), allocatable :: ucell(:,:) ! uReconstructZonal(nVertLevels, nCells) -real(r8), allocatable :: vcell(:,:) ! uReconstructMeridional(nVertLevels, nCells) -real(r8), allocatable :: data_2d_array(:,:) ! temporary -logical :: both -integer :: zonal, meridional - -if ( .not. module_initialized ) call static_init_model - -! get the ivar values for the zonal and meridional wind fields -call winds_present(zonal,meridional,both) -if (.not. both) call error_handler(E_ERR, 'update_wind_components', & - 'internal error: wind fields not found', source, revision, revdate) - -allocate( u(nVertLevels, nEdges)) -allocate(ucell(nVertLevels, nCells)) -allocate(vcell(nVertLevels, nCells)) - -! if doing increments, read in 'u' (edge normal winds), plus the uReconstructZonal -! and uReconstructMeridional fields from the mpas analysis netcdf file. - -if (use_increments_for_u_update) then - call read_2d_from_nc_file(ncid, filename, 'u', u) - call read_2d_from_nc_file(ncid, filename, 'uReconstructZonal', ucell) - call read_2d_from_nc_file(ncid, filename, 'uReconstructMeridional', vcell) - - if ((debug > 8) .and. do_output()) then - write(*,*) - write(*,*)'update_winds: org u range ',minval(u), maxval(u) - write(*,*)'update_winds: org zonal range ',minval(ucell), maxval(ucell) - write(*,*)'update_winds: org meridional range ',minval(vcell), maxval(vcell) - endif - - ! compute the increments compared to the updated values in the state vector - - allocate(data_2d_array(nVertLevels, nCells)) - - ! write updated reconstructed winds to restart file for later use. - call vector_to_prog_var(state_vector, zonal, data_2d_array) - call put_u(ncid, filename, data_2d_array, 'uReconstructZonal') - ucell = data_2d_array - ucell - - call vector_to_prog_var(state_vector, meridional, data_2d_array) - call put_u(ncid, filename, data_2d_array, 'uReconstructMeridional') - vcell = data_2d_array - vcell - - deallocate(data_2d_array) - - ! this is by nedges, not ncells as above - allocate(data_2d_array(nVertLevels, nEdges)) - call uv_cell_to_edges(ucell, vcell, data_2d_array) - u(:,:) = u(:,:) + data_2d_array(:,:) - deallocate(data_2d_array) - - if ((debug > 8) .and. do_output()) then - write(*,*) - write(*,*)'update_winds: u increment range ',minval(ucell), maxval(ucell) - write(*,*)'update_winds: v increment range ',minval(vcell), maxval(vcell) - endif - -else - - ! The state vector has updated zonal and meridional wind components. - ! put them directly into the arrays. these are the full values, not - ! just increments. - call vector_to_prog_var(state_vector, zonal, ucell) - call vector_to_prog_var(state_vector, meridional, vcell) - - call uv_cell_to_edges(ucell, vcell, u, .true.) - - if ((debug > 8) .and. do_output()) then - write(*,*) - write(*,*)'update_winds: u values range ',minval(ucell), maxval(ucell) - write(*,*)'update_winds: v values range ',minval(vcell), maxval(vcell) - endif - - call put_u(ncid, filename, ucell, 'uReconstructZonal') - call put_u(ncid, filename, vcell, 'uReconstructMeridional') -endif - -if ((debug > 8) .and. do_output()) then - write(*,*) - write(*,*)'update_winds: u after update:',minval(u), maxval(u) -endif - -! Write back to the mpas analysis file. - -call put_u(ncid, filename, u, 'u') - -deallocate(ucell, vcell, u) - -end subroutine update_wind_components - - -!------------------------------------------------------------------ - -subroutine read_2d_from_nc_file(ncid, filename, varname, data) - integer, intent(in) :: ncid - character(len=*), intent(in) :: filename - character(len=*), intent(in) :: varname - real(r8), intent(out) :: data(:,:) - -! -! Read the values for all dimensions but the time dimension. -! Only read the last time (if more than 1 present) -! - -integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs, mystart, mycount -character(len=NF90_MAX_NAME) :: dimname -integer :: VarID, numdims, dimlen, i - -call nc_check(nf90_inq_varid(ncid, varname, VarID), & - 'read_from_nc_file', & - 'inq_varid '//trim(varname)//' '//trim(filename)) - -call nc_check(nf90_inquire_variable(ncid, VarID, dimids=dimIDs, ndims=numdims), & - 'read_from_nc_file', & - 'inquire '//trim(varname)//' '//trim(filename)) - -do i=1, numdims - write(string1,*)'inquire length for dimension ',i - call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), len=dimlen, name=dimname), & - 'read_2d_from_nc_file', & - trim(string1)//' '//trim(filename)) - if (trim(dimname) == 'Time') then - mystart(i) = dimlen - mycount(numdims) = 1 - else - mystart(i) = 1 - mycount(i) = dimlen - endif -enddo - -call nc_check( nf90_get_var(ncid, VarID, data, & - start=mystart(1:numdims), count=mycount(1:numdims)), & - 'read_2d_from_nc_file', & - 'get_var u '//trim(filename)) - -end subroutine read_2d_from_nc_file - - -!------------------------------------------------------------------ - -subroutine put_u(ncid, filename, u, vchar) - character(len=*) :: vchar - integer, intent(in) :: ncid - character(len=*), intent(in) :: filename - real(r8), intent(in) :: u(:,:) ! u(nVertLevels, nEdges) - -! Put the newly updated 'u' field back into the netcdf file. - -integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs, mystart, mycount, numu -integer :: VarID, numdims, nDimensions, nVariables, nAttributes, unlimitedDimID -integer :: ntimes, i - -if ( .not. module_initialized ) call static_init_model - -string2 = trim(filename)//' '//trim(vchar) - -call nc_check(nf90_Inquire(ncid,nDimensions,nVariables,nAttributes,unlimitedDimID), & - 'put_u', 'inquire '//trim(filename)) - -call nc_check(nf90_inquire_dimension(ncid, unlimitedDimID, len=ntimes), & - 'put_u', 'inquire time dimension length '//trim(filename)) - -call nc_check(nf90_inq_varid(ncid, trim(vchar), VarID), & - 'put_u', 'inq_varid '//trim(string2)) - -call nc_check(nf90_inquire_variable(ncid, VarID, dimids=dimIDs, ndims=numdims), & - 'put_u', 'inquire '//trim(string2)) - -do i=1, numdims - write(string1,'(''inquire dimension'',i2,A)') i,trim(string2) - call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), len=numu(i)), & - 'put_u', string1) -enddo - -! for all but the time dimension, read all the values. -! for time read only the last one (if more than 1 present) -mystart = 1 -mystart(numdims) = ntimes -mycount = numu -mycount(numdims) = 1 - -call nc_check(nf90_put_var(ncid, VarID, u, start=mystart, count=mycount), & - 'put_u', 'put_var '//trim(string2)) - - -! A little sanity check - -if ((debug > 8) .and. do_output()) then - - write(*,*) - write(*,*) trim(vchar),' range ',minval(u), maxval(u) - -endif - -end subroutine put_u - - -!------------------------------------------------------------------ - -subroutine vector_to_1d_prog_var(x, ivar, data_1d_array) - -! convert the values from a 1d array, starting at an offset, -! into a 1d array. - -real(r8), dimension(:), intent(in) :: x -integer, intent(in) :: ivar -real(r8), dimension(:), intent(out) :: data_1d_array - -integer :: start_offset, end_offset - -start_offset = progvar(ivar)%index1 -end_offset = start_offset + size(data_1d_array) - 1 - -data_1d_array = x(start_offset:end_offset) - -end subroutine vector_to_1d_prog_var - - -!------------------------------------------------------------------ - -subroutine vector_to_2d_prog_var(x, ivar, data_2d_array) - -! convert the values from a 1d array, starting at an offset, -! into a 2d array. - -real(r8), dimension(:), intent(in) :: x -integer, intent(in) :: ivar -real(r8), dimension(:,:), intent(out) :: data_2d_array - -integer :: start_offset, end_offset - -start_offset = get_index_start(anl_domid, ivar) -end_offset = get_index_end(anl_domid, ivar) - -data_2d_array = reshape(x(start_offset:end_offset), shape(data_2d_array)) - -end subroutine vector_to_2d_prog_var - - -!------------------------------------------------------------------ - -subroutine bdy_vector_to_prog_var(x, ivar, data_2d_array) - -! convert the values from a 1d array, starting at an offset, -! into a 2d array. - -real(r8), dimension(:), intent(in) :: x -integer, intent(in) :: ivar -real(r8), dimension(:,:), intent(out) :: data_2d_array - -integer :: start_offset, end_offset - -start_offset = get_index_start(lbc_domid, ivar) -end_offset = get_index_end(lbc_domid, ivar) - -data_2d_array = reshape(x(start_offset:end_offset), shape(data_2d_array)) - -end subroutine bdy_vector_to_prog_var - - -!------------------------------------------------------------------ - -subroutine vector_to_3d_prog_var(x, ivar, data_3d_array) - -! convert the values from a 1d array, starting at an offset, -! into a 3d array. - -real(r8), dimension(:), intent(in) :: x -integer, intent(in) :: ivar -real(r8), dimension(:,:,:), intent(out) :: data_3d_array - -integer :: start_offset, end_offset - -start_offset = progvar(ivar)%index1 -end_offset = start_offset + size(data_3d_array) - 1 - -data_3d_array = reshape(x(start_offset:end_offset), shape(data_3d_array)) - -end subroutine vector_to_3d_prog_var - - -!------------------------------------------------------------------ - -subroutine prog_var_1d_to_vector(data_1d_array, x, ivar) - -! convert the values from a 1d array into a 1d array -! starting at an offset. - -real(r8), dimension(:), intent(in) :: data_1d_array -real(r8), dimension(:), intent(inout) :: x -integer, intent(in) :: ivar - -integer :: start_offset, end_offset - -start_offset = progvar(ivar)%index1 -end_offset = start_offset + size(data_1d_array) - 1 - -x(start_offset:end_offset) = data_1d_array - -end subroutine prog_var_1d_to_vector - - -!------------------------------------------------------------------ - -subroutine prog_var_2d_to_vector(data_2d_array, x, ivar) - -! convert the values from a 2d array into a 1d array -! starting at an offset. - -real(r8), dimension(:,:), intent(in) :: data_2d_array -real(r8), dimension(:), intent(inout) :: x -integer, intent(in) :: ivar - -integer :: start_offset, end_offset - -start_offset = progvar(ivar)%index1 -end_offset = start_offset + size(data_2d_array) - 1 - -x(start_offset:end_offset) = reshape(data_2d_array, (/ size(data_2d_array) /) ) - -end subroutine prog_var_2d_to_vector - +latCell = latCell * rad2deg +lonCell = lonCell * rad2deg +where (latCell > 90.0_r8) latCell = 90.0_r8 +where (latCell < -90.0_r8) latCell = -90.0_r8 -!------------------------------------------------------------------ +call nc_get_variable(ncid, 'dcEdge', dcEdge, routine) +call nc_get_variable(ncid, 'zgrid', zGridFace, routine) +call nc_get_variable(ncid, 'cellsOnVertex', cellsOnVertex, routine) +call nc_get_variable(ncid, 'xland', xland, routine) +call nc_get_variable(ncid, 'seaice', seaice, routine) +call nc_get_variable(ncid, 'skintemp', skintemp, routine) -subroutine prog_var_3d_to_vector(data_3d_array, x, ivar) +dxmax = maxval(dcEdge) ! max grid resolution in meters -! convert the values from a 2d array into a 1d array -! starting at an offset. +call nc_get_variable(ncid, 'edgeNormalVectors', edgeNormalVectors, routine) +call nc_get_variable(ncid, 'nEdgesOnCell', nEdgesOnCell, routine) +call nc_get_variable(ncid, 'edgesOnCell', edgesOnCell, routine) +call nc_get_variable(ncid, 'cellsOnEdge', cellsOnEdge, routine) +call nc_get_variable(ncid, 'verticesOnCell', verticesOnCell, routine) -real(r8), dimension(:,:,:), intent(in) :: data_3d_array -real(r8), dimension(:), intent(inout) :: x -integer, intent(in) :: ivar +if(data_on_edges) then + call nc_get_variable(ncid, 'lonEdge', lonEdge, routine) + call nc_get_variable(ncid, 'latEdge', latEdge, routine) -integer :: start_offset, end_offset + latEdge = latEdge * rad2deg + lonEdge = lonEdge * rad2deg + where (latEdge > 90.0_r8) latEdge = 90.0_r8 + where (latEdge < -90.0_r8) latEdge = -90.0_r8 -start_offset = progvar(ivar)%index1 -end_offset = start_offset + size(data_3d_array) - 1 + call nc_get_variable(ncid, 'xEdge', xEdge, routine) + call nc_get_variable(ncid, 'yEdge', yEdge, routine) + call nc_get_variable(ncid, 'zEdge', zEdge, routine) +endif -x(start_offset:end_offset) = reshape(data_3d_array, (/ size(data_3d_array) /)) +call nc_get_variable(ncid, 'xVertex', xVertex, routine) +call nc_get_variable(ncid, 'yVertex', yVertex, routine) +call nc_get_variable(ncid, 'zVertex', zVertex, routine) -end subroutine prog_var_3d_to_vector +if (nc_variable_exists(ncid, 'maxLevelCell')) then + allocate(maxLevelCell(nCells)) + call nc_get_variable(ncid, 'maxLevelCell', maxlevelCell, routine) + all_levels_exist_everywhere = .false. +endif +end subroutine get_grid !------------------------------------------------------------------ +subroutine verify_state_variables(ncid, filename, ngood, qty_list, variable_bounds) -!>@todo fill in the inputs we need for the add_domain() routine - -subroutine verify_state_variables( state_variables, ncid, filename, ngood, table ) - -character(len=*), dimension(:), intent(in) :: state_variables integer, intent(in) :: ncid character(len=*), intent(in) :: filename integer, intent(out) :: ngood -character(len=*), dimension(:,:), intent(out) :: table +integer, intent(out) :: qty_list(:) +real(r8), intent(out) :: variable_bounds(:,:) -integer :: nrows, ncols, i, j, VarID +integer :: i, j, VarID integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs character(len=NF90_MAX_NAME) :: varname, dimname character(len=NF90_MAX_NAME) :: dartstr +real(r8) :: lower_bound, upper_bound +character(len=10) :: bound integer :: dimlen, numdims logical :: u_already_in_list logical :: failure if ( .not. module_initialized ) call static_init_model -failure = .FALSE. ! perhaps all with go well +failure = .FALSE. ! this is for unsupported dimensions u_already_in_list = .FALSE. -nrows = size(table,1) -ncols = size(table,2) +variable_bounds(:,:) = MISSING_R8 ngood = 0 -MyLoop : do i = 1, nrows +MyLoop : do i = 1, MAX_STATE_VARIABLES + + varname = trim(mpas_state_variables(2*i -1)) + dartstr = trim(mpas_state_variables(2*i )) + variable_table(i,1) = trim(varname) + variable_table(i,2) = trim(dartstr) !HK @todo to_upper - varname = trim(state_variables(2*i -1)) - dartstr = trim(state_variables(2*i )) - table(i,1) = trim(varname) - table(i,2) = trim(dartstr) + if (varname == 'u') has_edge_u = .true. + if (varname == 'uReconstructZonal' .or. & + varname == 'uReconstructMeridional') has_uvreconstruct = .true. - if ( table(i,1) == ' ' .and. table(i,2) == ' ' ) exit MyLoop ! Found end of list. + if ( variable_table(i,1) == ' ' .and. variable_table(i,2) == ' ' ) exit MyLoop ! Found end of list. - if ( table(i,1) == ' ' .or. table(i,2) == ' ' ) then + if ( variable_table(i,1) == ' ' .or. variable_table(i,2) == ' ' ) then string1 = 'mpas_vars_nml:model state_variables not fully specified' - call error_handler(E_ERR,'verify_state_variables',string1,source,revision,revdate) + call error_handler(E_ERR,'verify_state_variables',string1,source) endif ! Make sure variable exists in model analysis variable list @@ -3847,7 +1728,7 @@ subroutine verify_state_variables( state_variables, ncid, filename, ngood, table 'verify_state_variables', trim(string2)) ! Make sure variable is defined by (Time,nCells) or (Time,nCells,vertical) - ! unable to support Edges or Vertices at this time. + ! unable to support Edges or Vertices at this time.m !HK @todo this is not true. call nc_check(nf90_inquire_variable(ncid, VarID, dimids=dimIDs, ndims=numdims), & 'verify_state_variables', 'inquire '//trim(string1)) @@ -3872,7 +1753,7 @@ subroutine verify_state_variables( state_variables, ncid, filename, ngood, table ! supported - do nothing case default write(string2,'(''unsupported dimension '',a,'' in '',a)') trim(dimname),trim(string1) - call error_handler(E_MSG,'verify_state_variables',string2,source,revision,revdate) + call error_handler(E_MSG,'verify_state_variables',string2,source) failure = .TRUE. end select @@ -3880,226 +1761,74 @@ subroutine verify_state_variables( state_variables, ncid, filename, ngood, table if (failure) then string2 = 'unsupported dimension(s) are fatal' - call error_handler(E_ERR,'verify_state_variables',string2,source,revision,revdate) + call error_handler(E_ERR,'verify_state_variables',string2,source) endif - ! Make sure DART kind is valid - - if( get_index_for_quantity(dartstr) < 0 ) then + qty_list(i) = get_index_for_quantity(dartstr) + if( qty_list(i) < 0 ) then write(string1,'(''there is no qty <'',a,''> in obs_kind_mod.f90'')') trim(dartstr) - call error_handler(E_ERR,'verify_state_variables',string1,source,revision,revdate) - endif - - ! Record the contents of the DART state vector - - if ((debug > 0) .and. do_output()) then - write(logfileunit,*)'variable ',i,' is ',trim(table(i,1)), ' ', trim(table(i,2)) - write( * ,*)'variable ',i,' is ',trim(table(i,1)), ' ', trim(table(i,2)) + call error_handler(E_ERR,'verify_state_variables',string1,source) endif ! Keep track of whether U (edge normal winds) is part of the user-specified field list - if ((table(i, 1) == 'u') .or. (table(i, 1) == 'U')) u_already_in_list = .true. + if ((variable_table(i, 1) == 'u') .or. (variable_table(i, 1) == 'U')) u_already_in_list = .true. + + ! check variable bounds + do j = 1, MAX_STATE_VARIABLES + if (trim(mpas_state_bounds(1, j)) == ' ') exit + if (trim(mpas_state_bounds(1, j)) == trim(varname)) then + + bound = trim(mpas_state_bounds(2, j)) + if ( bound /= 'NULL' .and. bound /= '' ) then + read(bound,'(d16.8)') lower_bound + else + lower_bound = MISSING_R8 + endif + + bound = trim(mpas_state_bounds(3, j)) + if ( bound /= 'NULL' .and. bound /= '' ) then + read(bound,'(d16.8)') upper_bound + else + upper_bound = MISSING_R8 + endif + variable_bounds(i,1) = lower_bound + variable_bounds(i,2) = upper_bound + endif + enddo ngood = ngood + 1 enddo MyLoop -if (ngood == nrows) then - string1 = 'WARNING: There is a possibility you need to increase ''max_state_variables''' +if (ngood == MAX_STATE_VARIABLES) then + string1 = 'WARNING: There is a possibility you need to increase ''MAX_STATE_VARIABLES''' write(string2,'(''WARNING: you have specified at least '',i4,'' perhaps more.'')')ngood - call error_handler(E_MSG,'verify_state_variables',string1,source,revision,revdate,text2=string2) -endif - -! if this flag is true and the user hasn't said U should be in the state, -! add it to the list. -if (add_u_to_state_list .and. .not. u_already_in_list) then - ngood = ngood + 1 - table(ngood,1) = "u" - table(ngood,2) = "QTY_EDGE_NORMAL_SPEED" + call error_handler(E_MSG,'verify_state_variables',string1,source,text2=string2) endif end subroutine verify_state_variables !------------------------------------------------------------------ +function is_edgedata_in_state_vector() -!> @todo FIXME if you can call this *after* add_domain() has been -!> called, then we could use state structure calls for this. -!> -!> but right now, it's called first. so pass in the namelist -!> and boundary lists and base the decision on those. -!> -!> this routine can only be called after the namelist is read. -!> also, it's called BY static_init_model() so it can't call -!> back into it. -!> - -function is_edgedata_in_state_vector(state_list, bdy_list) - character(len=*), intent(in) :: state_list(max_state_variables, num_state_table_columns) - character(len=*), intent(in) :: bdy_list(max_state_variables) - logical :: is_edgedata_in_state_vector - -integer :: i - -StateLoop : do i = 1, max_state_variables - - if (state_list(i, 1) == ' ') exit StateLoop ! Found end of list. - - if (state_list(i, 1) == 'u') then - is_edgedata_in_state_vector = .true. - return - endif - -enddo StateLoop +logical :: is_edgedata_in_state_vector -! if U is not in the state, does it matter if U is -! in the boundary file? yes, return true if so. -BdyLoop : do i = 1, max_state_variables - - if (bdy_list(i) == ' ') exit BdyLoop ! Found end of list. - - if (bdy_list(i) == 'lbc_u') then - is_edgedata_in_state_vector = .true. - return - endif - -enddo BdyLoop +integer :: dom, var, dim is_edgedata_in_state_vector = .false. -end function is_edgedata_in_state_vector - - -!------------------------------------------------------------------ - -subroutine dump_progvar(ivar) - -! dump the contents of the metadata for an individual entry. -! expected to be called in a loop or called for entries of interest. - -integer, intent(in) :: ivar - -!%! type progvartype -!%! private -!%! character(len=NF90_MAX_NAME) :: varname -!%! integer, dimension(NF90_MAX_VAR_DIMS) :: dimlens -!%! integer :: xtype ! netCDF variable type (NF90_double, etc.) -!%! integer :: numdims ! number of dims - excluding TIME -!%! integer :: numvertical ! number of vertical levels in variable -!%! integer :: numcells ! number of horizontal locations (typically cell centers) -!%! integer :: numedges -!%! logical :: ZonHalf ! vertical coordinate has dimension nVertLevels -!%! 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 -!%! character(len=obstypelength) :: kind_string -!%! logical :: clamping ! does variable need to be range-restricted before -!%! real(r8) :: range(2) ! being stuffed back into MPAS analysis file. -!%! end type progvartype - -integer :: i - -! take care of parallel runs where we only want a single copy of -! the output. -if (.not. do_output()) return - -write(logfileunit,*) -write( * ,*) -write(logfileunit,*) 'variable number ',ivar,' is ',trim(progvar(ivar)%varname) -write( * ,*) 'variable number ',ivar,' is ',trim(progvar(ivar)%varname) -write(logfileunit,*) ' xtype ',progvar(ivar)%xtype -write( * ,*) ' xtype ',progvar(ivar)%xtype -write(logfileunit,*) ' dimlens ',progvar(ivar)%dimlens(1:progvar(ivar)%numdims) -write( * ,*) ' dimlens ',progvar(ivar)%dimlens(1:progvar(ivar)%numdims) -write(logfileunit,*) ' numdims ',progvar(ivar)%numdims -write( * ,*) ' numdims ',progvar(ivar)%numdims -write(logfileunit,*) ' numvertical ',progvar(ivar)%numvertical -write( * ,*) ' numvertical ',progvar(ivar)%numvertical -write(logfileunit,*) ' numcells ',progvar(ivar)%numcells -write( * ,*) ' numcells ',progvar(ivar)%numcells -write(logfileunit,*) ' numedges ',progvar(ivar)%numedges -write( * ,*) ' numedges ',progvar(ivar)%numedges -write(logfileunit,*) ' ZonHalf ',progvar(ivar)%ZonHalf -write( * ,*) ' ZonHalf ',progvar(ivar)%ZonHalf -write(logfileunit,*) ' varsize ',progvar(ivar)%varsize -write( * ,*) ' varsize ',progvar(ivar)%varsize -write(logfileunit,*) ' index1 ',progvar(ivar)%index1 -write( * ,*) ' index1 ',progvar(ivar)%index1 -write(logfileunit,*) ' indexN ',progvar(ivar)%indexN -write( * ,*) ' indexN ',progvar(ivar)%indexN -write(logfileunit,*) ' dart_kind ',progvar(ivar)%dart_kind -write( * ,*) ' dart_kind ',progvar(ivar)%dart_kind -write(logfileunit,*) ' kind_string ',progvar(ivar)%kind_string -write( * ,*) ' kind_string ',progvar(ivar)%kind_string -write(logfileunit,*) ' clamping ',progvar(ivar)%clamping -write( * ,*) ' clamping ',progvar(ivar)%clamping -write(logfileunit,*) ' clamp range ',progvar(ivar)%range -write( * ,*) ' clamp range ',progvar(ivar)%range -write(logfileunit,*) ' clamp fail ',progvar(ivar)%out_of_range_fail -write( * ,*) ' clamp fail ',progvar(ivar)%out_of_range_fail -do i = 1,progvar(ivar)%numdims - write(logfileunit,*) ' dimension/length/name ',i,progvar(ivar)%dimlens(i),trim(progvar(ivar)%dimname(i)) - write( * ,*) ' dimension/length/name ',i,progvar(ivar)%dimlens(i),trim(progvar(ivar)%dimname(i)) -enddo - -end subroutine dump_progvar - -!------------------------------------------------------------------ - -subroutine print_variable_ranges(x,title) - -! given a state vector, print out the min and max -! data values for the variables in the vector. - -real(r8), intent(in) :: x(:) -character(len=*), optional, intent(in) :: title - -integer :: ivar - -if (present(title)) write(*,*) trim(title) - -do ivar = 1, nfields - call print_minmax(ivar, x) +do dom = 1, get_num_domains() + do var = 1, get_num_variables(dom) + do dim = 1, get_num_dims(dom, var) + if (get_dim_name(dom, var, dim) == 'nEdges') then + is_edgedata_in_state_vector = .true. + return + endif + enddo + enddo enddo -end subroutine print_variable_ranges - -!------------------------------------------------------------------ - -subroutine print_minmax(ivar, x) - -! given an index and a state vector, print out the min and max -! data values for the items corresponding to that progvar index. - -integer, intent(in) :: ivar -real(r8), intent(in) :: x(:) - -write(string1, '(A,A32,2F16.7)') 'data min/max ', trim(progvar(ivar)%varname), & - minval(x(progvar(ivar)%index1:progvar(ivar)%indexN)), & - maxval(x(progvar(ivar)%index1:progvar(ivar)%indexN)) - -call error_handler(E_MSG, '', string1, source,revision,revdate) - -end subroutine print_minmax - - -!------------------------------------------------------------------ - -function FindTimeDimension(ncid) result(timedimid) - -! Find the Time Dimension ID in a netCDF file. -! If there is none - (spelled the obvious way) - the routine -! returns a negative number. You don't HAVE to have a TIME dimension. - -integer :: timedimid -integer, intent(in) :: ncid - -integer :: nc_rc - -TimeDimID = -1 ! same as the netCDF library routines. -nc_rc = nf90_inq_dimid(ncid,'Time',dimid=TimeDimID) - -end function FindTimeDimension +end function is_edgedata_in_state_vector !------------------------------------------------------------------ @@ -4128,211 +1857,10 @@ subroutine winds_present(zonal,meridional,both) ! only one present - error. write(string1,*) 'both components for U/V reconstructed winds must be in state vector' write(string2,*) 'cannot have only one of uReconstructMeridional and uReconstructZonal' -call error_handler(E_ERR,'winds_present',string1,source,revision,revdate,text2=string2) +call error_handler(E_ERR,'winds_present',string1,source,text2=string2) end subroutine winds_present - -!------------------------------------------------------------ -subroutine get_variable_bounds(bounds_table, ivar) - -! matches MPAS variable name in bounds table to assign -! the bounds if they exist. otherwise sets the bounds -! to missing_r8 -! -! SYHA (May-30-2013) -! Adopted from wrf/model_mod.f90 after adding mpas_state_bounds in mpas_vars_nml. - -! bounds_table is the global mpas_state_bounds -character(len=*), intent(in) :: bounds_table(num_bounds_table_columns, max_state_variables) -integer, intent(in) :: ivar - -! local variables -character(len=50) :: bounds_varname, bound -character(len=10) :: clamp_or_fail -real(r8) :: lower_bound, upper_bound -integer :: n - -n = 1 -do while ( trim(bounds_table(1,n)) /= 'NULL' .and. trim(bounds_table(1,n)) /= '' ) - - bounds_varname = trim(bounds_table(1,n)) - - if ( bounds_varname == trim(progvar(ivar)%varname) ) then - - bound = trim(bounds_table(2,n)) - if ( bound /= 'NULL' .and. bound /= '' ) then - read(bound,'(d16.8)') lower_bound - else - lower_bound = missing_r8 - endif - - bound = trim(bounds_table(3,n)) - if ( bound /= 'NULL' .and. bound /= '' ) then - read(bound,'(d16.8)') upper_bound - else - upper_bound = missing_r8 - endif - - ! How do we want to handle out of range values? - ! Set them to predefined limits (clamp) or simply fail (fail). - clamp_or_fail = trim(bounds_table(4,n)) - if ( clamp_or_fail == 'NULL' .or. clamp_or_fail == '') then - write(string1, *) 'instructions for CLAMP_or_FAIL on ', & - trim(bounds_varname), ' are required' - call error_handler(E_ERR,'get_variable_bounds',string1, & - source,revision,revdate) - else if ( clamp_or_fail == 'CLAMP' ) then - progvar(ivar)%out_of_range_fail = .FALSE. - else if ( clamp_or_fail == 'FAIL' ) then - progvar(ivar)%out_of_range_fail = .TRUE. - else - write(string1, *) 'last column must be "CLAMP" or "FAIL" for ', & - trim(bounds_varname) - call error_handler(E_ERR,'get_variable_bounds',string1, & - source,revision,revdate, text2='found '//trim(clamp_or_fail)) - endif - - ! Assign the clamping information into the variable - progvar(ivar)%clamping = .true. - progvar(ivar)%range = (/ lower_bound, upper_bound /) - - if ((debug > 0) .and. do_output()) then - write(*,*) 'In get_variable_bounds assigned ', trim(progvar(ivar)%varname) - write(*,*) ' clamping range ',progvar(ivar)%range - endif - - ! we found the progvar entry and set the values. return here. - return - endif - - n = n + 1 - -enddo !n - -! we got through all the entries in the bounds table and did not -! find any instructions for this variable. set the values to indicate -! we are not doing any processing when we write the updated state values -! back to the model restart file. - -progvar(ivar)%clamping = .false. -progvar(ivar)%range = missing_r8 -progvar(ivar)%out_of_range_fail = .false. ! should be unused so setting shouldn't matter - -return - -end subroutine get_variable_bounds - -!------------------------------------------------------------ - -subroutine get_index_range_string(string,index1,indexN) - -! Determine where a particular DART kind (string) exists in the -! DART state vector. - -character(len=*), intent(in) :: string -integer, intent(out) :: index1 -integer, optional, intent(out) :: indexN - -integer :: i - -index1 = 0 -if (present(indexN)) indexN = 0 - -FieldLoop : do i=1,nfields - if (progvar(i)%kind_string /= trim(string)) cycle FieldLoop - index1 = progvar(i)%index1 - if (present(indexN)) indexN = progvar(i)%indexN - exit FieldLoop -enddo FieldLoop - -if (index1 == 0) then - write(string1,*) 'Problem, cannot find indices for '//trim(string) - call error_handler(E_ERR,'get_index_range_string',string1,source,revision,revdate) -endif -end subroutine get_index_range_string - - -!------------------------------------------------------------------ - -subroutine get_index_range_int(dartkind,index1,indexN) - -! Determine where a particular DART kind (integer) exists in the -! DART state vector. - -integer , intent(in) :: dartkind -integer(i8), intent(out) :: index1 -integer(i8), optional, intent(out) :: indexN - -integer :: i -character(len=obstypelength) :: string - -index1 = 0 -if (present(indexN)) indexN = 0 - -FieldLoop : do i=1,nfields - if (progvar(i)%dart_kind /= dartkind) cycle FieldLoop - index1 = progvar(i)%index1 - if (present(indexN)) indexN = progvar(i)%indexN - exit FieldLoop -enddo FieldLoop - -string = get_name_for_quantity(dartkind) - -if (index1 == 0) then - write(string1,*) 'Problem, cannot find indices for kind ',dartkind,trim(string) - call error_handler(E_ERR,'get_index_range_int',string1,source,revision,revdate) -endif - -end subroutine get_index_range_int - - -!------------------------------------------------------------------ - -function get_progvar_index_from_kind(dartkind) - -! Determine what index a particular DART kind (integer) is in the -! progvar array. -integer :: get_progvar_index_from_kind -integer, intent(in) :: dartkind - -integer :: i - -FieldLoop : do i=1,nfields - if (progvar(i)%dart_kind /= dartkind) cycle FieldLoop - get_progvar_index_from_kind = i - return -enddo FieldLoop - -get_progvar_index_from_kind = -1 - -end function get_progvar_index_from_kind - - -!------------------------------------------------------------------ - -function get_index_from_varname(varname) - -! Determine what index corresponds to the given varname -! if name not in state vector, return -1 -- not an error. - -integer :: get_index_from_varname -character(len=*), intent(in) :: varname - -integer :: i - -FieldLoop : do i=1,nfields - if (trim(progvar(i)%varname) == trim(varname)) then - get_index_from_varname = i - return - endif -enddo FieldLoop - -get_index_from_varname = -1 -return - -end function get_index_from_varname - !------------------------------------------------------------------ subroutine compute_pressure_at_loc(state_handle, ens_size, location, ploc, istatus) @@ -4357,12 +1885,12 @@ subroutine compute_pressure_at_loc(state_handle, ens_size, location, ploc, istat real(r8), dimension(3, ens_size) :: values real(r8), dimension(ens_size) :: tk type(location_type), dimension(ens_size) :: new_location -integer :: ivars(3) +integer :: ivars(1) integer :: e ! loop index ! default is failure ploc = MISSING_R8 -istatus = 99 +istatus = CRITICAL_ERROR ! base location llv = get_location(location) @@ -4375,7 +1903,7 @@ subroutine compute_pressure_at_loc(state_handle, ens_size, location, ploc, istat new_location = location ! the quantity doesn't matter for this call but we have to pass in something. - call convert_vert_distrib(state_handle, ens_size, new_location, QTY_TEMPERATURE, VERTISPRESSURE, istatus) + call convert_vert(state_handle, ens_size, new_location, QTY_TEMPERATURE, VERTISPRESSURE, istatus) do e = 1, ens_size if(istatus(e) == 0) then @@ -4386,35 +1914,14 @@ subroutine compute_pressure_at_loc(state_handle, ens_size, location, ploc, istat else if(is_vertical(location, "SURFACE")) then - ivars(1) = get_progvar_index_from_kind(QTY_SURFACE_PRESSURE) + ivars(1) = get_varid_from_kind(anl_domid, QTY_SURFACE_PRESSURE) if ( ivars(1) >= 0 ) then call compute_scalar_with_barycentric(state_handle, ens_size, location, 1, ivars(1), values(1,:), istatus) where (istatus == 0) ploc(:) = values(1,:) else - istatus = 88 ! required quantity not in state vector - return - -!%! !>original code: -!%! !>@todo FIXME: do we really want to do this if the vert is surface and -!%! !> the surface pressure field is not in the state? this is going to return -!%! !> the pressure at the midpoint of the first level, is it not? -!%! -!%! new_location(1) = set_location(llv(1), llv(2), 1.0_r8, VERTISLEVEL) -!%! -!%! ! Need to get base offsets for the potential temperature, density, and water -!%! ! vapor mixing fields in the state vector -!%! ivars(1) = get_progvar_index_from_kind(QTY_POTENTIAL_TEMPERATURE) -!%! ivars(2) = get_progvar_index_from_kind(QTY_DENSITY) -!%! ivars(3) = get_progvar_index_from_kind(QTY_VAPOR_MIXING_RATIO) -!%! -!%! call compute_scalar_with_barycentric (state_handle, ens_size, new_location(1), 3, ivars, values, istatus) -!%! if ( all(istatus /= 0) ) return -!%! -!%! ! Convert surface theta, rho, qv into pressure -!%! call compute_full_pressure(ens_size, values(1, :), values(2, :), values(3, :), ploc(:), tk(:), istatus(:) ) - + istatus = QTY_NOT_IN_STATE_VECTOR endif else if(is_vertical(location, "UNDEFINED")) then ! not error, but no exact vert loc either @@ -4423,12 +1930,7 @@ subroutine compute_pressure_at_loc(state_handle, ens_size, location, ploc, istat else call error_handler(E_ERR, 'compute_pressure:', 'internal error: unknown type of vertical', & - source, revision, revdate) -endif - -if(debug > 9 .and. do_output()) then - print *, 'compute_pressure_at_loc: base location ',llv(1:3) - print *, 'compute_pressure_at_loc: istatus, pressure ', istatus(1), ploc(1) + source) endif end subroutine compute_pressure_at_loc @@ -4449,7 +1951,7 @@ subroutine convert_vertical_obs(state_handle, num, locs, loc_qtys, loc_types, & integer :: i do i=1, num - call convert_vert_distrib(state_handle, 1, locs(i:i), loc_qtys(i), which_vert, status(i:i)) + call convert_vert(state_handle, 1, locs(i:i), loc_qtys(i), which_vert, status(i:i)) enddo end subroutine convert_vertical_obs @@ -4472,13 +1974,13 @@ subroutine convert_vertical_state(state_handle, num, locs, loc_qtys, loc_indx, & !> we are using code from get_state_meta_data to get !> the indices for cell id and vert level. this calls a -!> modified convert_vert_distrib() routine that doesn't +!> modified convert_vert() routine that doesn't !> have to search for the cell centers. istatus = 0 do i=1, num - call convert_vert_distrib_state(state_handle, 1, locs(i:i), loc_qtys(i), & + call convert_vert_state(state_handle, 1, locs(i:i), loc_qtys(i), & loc_indx(i), which_vert, status) ! save the first error we see - but continue to convert the rest @@ -4491,30 +1993,7 @@ end subroutine convert_vertical_state !------------------------------------------------------------------ !> code to convert an observation location's vertical coordinate type. - -subroutine convert_vert_distrib(state_handle, ens_size, location, obs_kind, ztypeout, istatus) - -! This subroutine converts a given obs vertical coordinate to -! the vertical localization coordinate type requested through the -! model_mod namelist. -! -! Notes: (1) obs_kind is only necessary to check whether the ob -! is an identity ob. -! (2) This subroutine can convert both obs' and state points' -! vertical coordinates. Remember that state points get -! their DART location information from get_state_meta_data -! which is called by filter_assim during the assimilation -! process. -! (3) state_handle is the relevant DART state vector for carrying out -! computations necessary for the vertical coordinate -! transformations. As the vertical coordinate is only used -! in distance computations, this is actually the "expected" -! vertical coordinate, so that computed distance is the -! "expected" distance. Thus, under normal circumstances, -! the state that is supplied to convert_vert_distrib should be the -! ensemble mean. Nevertheless, the subroutine has the -! functionality to operate on any DART state vector that -! is supplied to it. +subroutine convert_vert(state_handle, ens_size, location, obs_kind, ztypeout, istatus) type(ensemble_type), intent(in) :: state_handle integer, intent(in) :: ens_size @@ -4523,9 +2002,6 @@ subroutine convert_vert_distrib(state_handle, ens_size, location, obs_kind, ztyp integer, intent(in) :: ztypeout integer, intent(out) :: istatus(ens_size) -! zin and zout are the vert values coming in and going out. -! ztype{in,out} are the vert types as defined by the 3d sphere -! locations mod (location/threed_sphere/location_mod.f90) real(r8), dimension(3, ens_size) :: llv_loc real(r8), dimension(3,ens_size) :: zk_mid, values, fract, fdata integer, dimension(3,ens_size) :: k_low, k_up @@ -4538,8 +2014,7 @@ subroutine convert_vert_distrib(state_handle, ens_size, location, obs_kind, ztyp integer :: ztypein, i, e, n integer :: c(3), ivars(3) -! assume failure. -istatus = 1 +istatus = GENERAL_ERROR ! initialization k_low = 0.0_r8 @@ -4553,18 +2028,10 @@ subroutine convert_vert_distrib(state_handle, ens_size, location, obs_kind, ztyp location(:) = location(1) endif -! if the existing coord is already in the requested vertical units -! or if the vert is 'undef' which means no specifically defined -! vertical coordinate, return now. ztypein = nint(query_location(location(1), 'which_vert')) if ((ztypein == ztypeout) .or. (ztypein == VERTISUNDEF)) then istatus = 0 return -else - if (debug > 9 .and. do_output()) then - write(string1,'(A,3X,2I3)') 'ztypein, ztypeout:',ztypein,ztypeout - call error_handler(E_MSG, 'convert_vert_distrib',string1,source, revision, revdate) - endif endif ! we do need to convert the vertical. start by @@ -4573,9 +2040,6 @@ subroutine convert_vert_distrib(state_handle, ens_size, location, obs_kind, ztyp llv_loc(:, e) = get_location(location(e)) enddo -! the routines below will use zin as the incoming vertical value -! and zout as the new outgoing one. start out assuming failure -! (zout = missing) and wait to be pleasantly surprised when it works. zin(:) = llv_loc(3, :) zout(:) = missing_r8 @@ -4583,7 +2047,7 @@ subroutine convert_vert_distrib(state_handle, ens_size, location, obs_kind, ztyp ! with the requested type as out. do e = 1, ens_size if (zin(e) == missing_r8) then - location(e) = set_location(llv_loc(1, e),llv_loc(2, e),missing_r8,ztypeout) + location(e) = set_location(llv_loc(1, e), llv_loc(2, e), MISSING_R8, ztypeout) endif enddo ! if the entire ensemble has missing vertical values we can return now. @@ -4593,123 +2057,86 @@ subroutine convert_vert_distrib(state_handle, ens_size, location, obs_kind, ztyp return endif -! Convert the incoming vertical type (ztypein) into the vertical -! localization coordinate given in the namelist (ztypeout). -! Various incoming vertical types (ztypein) are taken care of -! inside find_vert_level. So we only check ztypeout here. - -! convert into: select case (ztypeout) - ! ------------------------------------------------------------ - ! outgoing vertical coordinate should be 'model level number' - ! ------------------------------------------------------------ case (VERTISLEVEL) - ! Identify the three cell ids (c) in the triangle enclosing the obs and - ! the vertical indices for the triangle at two adjacent levels (k_low and k_up) - ! and the fraction (fract) for vertical interpolation. - - call find_triangle (location(1), n, c, weights, istatus(1)) - if(istatus(1) /= 0) then - istatus(:) = istatus(1) - location(:) = set_location(llv_loc(1, 1),llv_loc(2, 1),missing_r8,ztypeout) - return - endif - call find_vert_indices (state_handle, ens_size, location(1), n, c, k_low, k_up, fract, istatus) - if( all(istatus /= 0) ) then - location(:) = set_location(llv_loc(1, 1),llv_loc(2, 1),missing_r8,ztypeout) - return - endif + ! Identify the three cell ids (c) in the triangle enclosing the obs and + ! the vertical indices for the triangle at two adjacent levels (k_low and k_up) + ! and the fraction (fract) for vertical interpolation. - zk_mid = k_low + fract - do e = 1, ens_size - if(istatus(e) == 0) then - zout(e) = sum(weights(:) * zk_mid(:, e)) + call find_triangle (location(1), n, c, weights, istatus(1)) + if(istatus(1) /= 0) then + istatus(:) = istatus(1) + location(:) = set_location(llv_loc(1, 1), llv_loc(2, 1), MISSING_R8, ztypeout) + return + endif + call find_vert_level (state_handle, ens_size, location(1), n, c, -1, k_low, k_up, fract, istatus) !HK @todo needs qty to indicate center + if( all(istatus /= 0) ) then + location(:) = set_location(llv_loc(1, 1), llv_loc(2, 1), MISSING_R8, ztypeout) + return endif - enddo - if (debug > 9 .and. do_output()) then - write(string2,'("Zk for member 1:", 3F8.2," => ",F8.2)') zk_mid(:,1),zout(1) - call error_handler(E_MSG, 'convert_vert_distrib',string2,source, revision, revdate) - endif + zk_mid = k_low + fract + do e = 1, ens_size + if(istatus(e) == 0) then + zout(e) = sum(weights(:) * zk_mid(:, e)) + endif + enddo - ! ------------------------------------------------------------ - ! outgoing vertical coordinate should be 'pressure' in Pa - ! ------------------------------------------------------------ case (VERTISPRESSURE) - ! Need to get base offsets for the potential temperature, density, and water - ! vapor mixing fields in the state vector - ivars(1) = get_progvar_index_from_kind(QTY_POTENTIAL_TEMPERATURE) - ivars(2) = get_progvar_index_from_kind(QTY_DENSITY) - ivars(3) = get_progvar_index_from_kind(QTY_VAPOR_MIXING_RATIO) - - if (any(ivars(1:3) < 0)) then - write(string1,*) 'Internal error, cannot find one or more of: theta, rho, qv' - call error_handler(E_ERR, 'convert_vert_distrib',string1,source, revision, revdate) - endif - - ! Get theta, rho, qv at the interpolated location - call compute_scalar_with_barycentric (state_handle, ens_size, location(1), 3, ivars, values, istatus) - if( all(istatus /= 0) ) then - location(:) = set_location(llv_loc(1, 1),llv_loc(2, 1),missing_r8,ztypeout) - return - endif + ivars(1) = get_varid_from_kind(anl_domid, QTY_POTENTIAL_TEMPERATURE) + ivars(2) = get_varid_from_kind(anl_domid, QTY_DENSITY) + ivars(3) = get_varid_from_kind(anl_domid, QTY_VAPOR_MIXING_RATIO) - ! Convert theta, rho, qv into pressure - call compute_full_pressure(ens_size, values(1,:), values(2,:), values(3,:), zout, tk, istatus) - if( all(istatus /= 0) ) then - if (debug > 1 .and. do_output()) then - write(string2,'("Failed in compute_full_pressure in VERTISPRESSURE")') - call error_handler(E_ALLMSG, 'convert_vert_distrib',string2,source, revision, revdate) + if (any(ivars(1:3) < 0)) then + write(string1,*) 'Internal error, cannot find one or more of: theta, rho, qv' + call error_handler(E_ERR, 'convert_vert',string1,source) endif - endif - ! ------------------------------------------------------------ - ! outgoing vertical coordinate should be 'height' in meters - ! ------------------------------------------------------------ - case (VERTISHEIGHT) + ! Get theta, rho, qv at the interpolated location + call compute_scalar_with_barycentric (state_handle, ens_size, location(1), 3, ivars, values, istatus) + if( all(istatus /= 0) ) then + location(:) = set_location(llv_loc(1, 1), llv_loc(2, 1), MISSING_R8, ztypeout) + return + endif - call find_triangle (location(1), n, c, weights, istatus(1)) - if(istatus(1) /= 0) then - istatus(:) = istatus(1) - location(:) = set_location(llv_loc(1, 1),llv_loc(2, 1),missing_r8,ztypeout) - return - endif - call find_vert_indices (state_handle, ens_size, location(1), n, c, k_low, k_up, fract, istatus) - if( all(istatus /= 0) ) then - location(:) = set_location(llv_loc(1, 1),llv_loc(2, 1),missing_r8,ztypeout) - return - endif + ! Convert theta, rho, qv into pressure + call compute_full_pressure(ens_size, values(1,:), values(2,:), values(3,:), zout, tk, istatus) - fdata = 0.0_r8 - do i = 1, n - where (istatus == 0) - fdata(i, :) = zGridCenter(k_low(i, :),c(i))*(1.0_r8 - fract(i, :)) + & - zGridCenter(k_up (i, :),c(i))*fract(i, :) - end where - enddo + case (VERTISHEIGHT) - ! now have vertically interpolated values at cell centers. - ! use horizontal weights to compute value at interp point. - do e = 1, ens_size - if(istatus(e) == 0) then - zout(e) = sum(weights(:) * fdata(:,e)) - else - zout(e) = missing_r8 + call find_triangle (location(1), n, c, weights, istatus(1)) + if(istatus(1) /= 0) then + istatus(:) = istatus(1) + location(:) = set_location(llv_loc(1, 1),llv_loc(2, 1),missing_r8,ztypeout) + return + endif + call find_vert_level (state_handle, ens_size, location(1), n, c, -1, k_low, k_up, fract, istatus) + if( all(istatus /= 0) ) then + location(:) = set_location(llv_loc(1, 1),llv_loc(2, 1),missing_r8,ztypeout) + return endif - enddo - if (debug > 9 .and. do_output()) then - write(string2,'("zout_in_height [m] for member 1:",F10.2)') zout(1) - call error_handler(E_MSG, 'convert_vert_distrib',string2,source, revision, revdate) - endif + fdata = 0.0_r8 + do i = 1, n + where (istatus == 0) + fdata(i, :) = zGridCenter(k_low(i, :),c(i))*(1.0_r8 - fract(i, :)) + & + zGridCenter(k_up (i, :),c(i))*fract(i, :) + end where + enddo + ! now have vertically interpolated values at cell centers. + ! use horizontal weights to compute value at interp point. + do e = 1, ens_size + if(istatus(e) == 0) then + zout(e) = sum(weights(:) * fdata(:,e)) + else + zout(e) = missing_r8 + endif + enddo - ! ------------------------------------------------------------ - ! outgoing vertical coordinate should be 'scale height' (a ratio) - ! ------------------------------------------------------------ case (VERTISSCALEHEIGHT) ! Scale Height is defined as: log(pressure) @@ -4732,14 +2159,12 @@ subroutine convert_vert_distrib(state_handle, ens_size, location, obs_kind, ztyp if (at_surf .and. do_norm) then zout = 0.0_r8 istatus(:) = 0 - goto 101 + return endif - ! Base offsets for the potential temperature, density, and water - ! vapor mixing fields in the state vector. - ivars(1) = get_progvar_index_from_kind(QTY_POTENTIAL_TEMPERATURE) - ivars(2) = get_progvar_index_from_kind(QTY_DENSITY) - ivars(3) = get_progvar_index_from_kind(QTY_VAPOR_MIXING_RATIO) + ivars(1) = get_varid_from_kind(anl_domid, QTY_POTENTIAL_TEMPERATURE) + ivars(2) = get_varid_from_kind(anl_domid, QTY_DENSITY) + ivars(3) = get_varid_from_kind(anl_domid, QTY_VAPOR_MIXING_RATIO) if (at_surf .or. do_norm) then ! we will need surface pressure @@ -4747,18 +2172,12 @@ subroutine convert_vert_distrib(state_handle, ens_size, location, obs_kind, ztyp surfloc(1) = set_location(llv_loc(1, 1), llv_loc(2, 1), 1.0_r8, VERTISLEVEL) call compute_scalar_with_barycentric (state_handle, ens_size, surfloc(1), 3, ivars, values, istatus) if( all(istatus /= 0) ) then - location(:) = set_location(llv_loc(1, 1),llv_loc(2, 1),missing_r8,ztypeout) + location(:) = set_location(llv_loc(1, 1), llv_loc(2, 1), MISSING_R8, ztypeout) return endif ! Convert surface theta, rho, qv into pressure call compute_full_pressure(ens_size, values(1,:), values(2,:), values(3,:), surfp, tk, istatus) - if( all(istatus /= 0) ) then - if (debug > 1 .and. do_output()) then - write(string2,'("Failed in compute_full_pressure in surfp")') - call error_handler(E_MSG, 'convert_vert_distrib',string2,source, revision, revdate) - endif - endif endif @@ -4769,12 +2188,6 @@ subroutine convert_vert_distrib(state_handle, ens_size, location, obs_kind, ztyp ! Convert theta, rho, qv into pressure call compute_full_pressure(ens_size, values(1,:), values(2,:), values(3,:), fullp, tk, istatus) - if( all(istatus /= 0) ) then - if (debug > 1 .and. do_output()) then - write(string2,'("Failed in compute_full_pressure in fullp")') - call error_handler(E_MSG, 'convert_vert_distrib',string2,source, revision, revdate) - endif - endif endif @@ -4803,87 +2216,48 @@ subroutine convert_vert_distrib(state_handle, ens_size, location, obs_kind, ztyp end where endif -101 continue - - if (debug > 9 .and. do_output()) then - write(string2,'("zout_in_scaleheight for member 1:",F10.2)') zout(1) - call error_handler(E_MSG, 'convert_vert_distrib',string2,source, revision, revdate) - endif - - ! ------------------------------------------------------- - ! outgoing vertical coordinate is unrecognized - ! ------------------------------------------------------- case default write(string1,*) 'Requested vertical coordinate not recognized: ', ztypeout - call error_handler(E_ERR,'convert_vert_distrib', string1, & - source, revision, revdate) + call error_handler(E_ERR,'convert_vert', string1, & + source) end select ! outgoing vert type ! Returned location do e = 1, ens_size if(istatus(e) == 0) then - location(e) = set_location(llv_loc(1, e),llv_loc(2, e),zout(e),ztypeout) + location(e) = set_location(llv_loc(1, e), llv_loc(2, e), zout(e), ztypeout) else - location(e) = set_location(llv_loc(1, e),llv_loc(2, e),missing_r8,ztypeout) + location(e) = set_location(llv_loc(1, e), llv_loc(2, e), MISSING_R8, ztypeout) endif enddo - -end subroutine convert_vert_distrib +end subroutine convert_vert !------------------------------------------------------------------ -!> code to convert an state location's vertical coordinate type. - -subroutine convert_vert_distrib_state(state_handle, ens_size, location, quantity, state_indx, ztypeout, istatus) - -! This subroutine converts a given state vertical coordinate to -! the vertical localization coordinate type requested through the -! model_mod namelist. -! -! (3) state_handle is the relevant DART state vector for carrying out -! computations necessary for the vertical coordinate -! transformations. As the vertical coordinate is only used -! in distance computations, this is actually the "expected" -! vertical coordinate, so that computed distance is the -! "expected" distance. Thus, under normal circumstances, -! state_handle that is supplied to convert_vert_distrib should be the -! ensemble mean. Nevertheless, the subroutine has the -! functionality to operate on any DART state vector that -! is supplied to it. +! given an index into state, convert the location into requested vertical type +subroutine convert_vert_state(state_handle, ens_size, location, quantity, state_indx, ztypeout, istatus) type(ensemble_type), intent(in) :: state_handle integer, intent(in) :: ens_size type(location_type), intent(inout) :: location(ens_size) ! because the verticals may differ integer, intent(in) :: quantity integer(i8), intent(in) :: state_indx -integer, intent(in) :: ztypeout +integer, intent(in) :: ztypeout ! requested vertical type integer, intent(out) :: istatus(ens_size) -! zin and zout are the vert values coming in and going out. -! ztype{in,out} are the vert types as defined by the 3d sphere -! locations mod (location/threed_sphere/location_mod.f90) real(r8), dimension(3, ens_size) :: llv_loc -real(r8), dimension(3,ens_size) :: zk_mid, values, fract, fdata -integer, dimension(3,ens_size) :: k_low, k_up +real(r8), dimension(3,ens_size) :: values real(r8), dimension(ens_size) :: zin, zout real(r8), dimension(ens_size) :: tk, fullp, surfp -logical :: at_surf, do_norm, on_bound +logical :: at_surf, do_norm type(location_type), dimension(ens_size) :: surfloc -real(r8) :: weights(3) integer :: ztypein, i, e, n, ndim integer :: c(3), ivars(3), vert_level -integer :: cellid !SYHA - -! assume failure. -istatus = 1 - -! initialization -k_low = 0.0_r8 -k_up = 0.0_r8 -weights = 0.0_r8 +integer :: cellid +istatus = GENERAL_ERROR ! if the existing coord is already in the requested vertical units ! or if the vert is 'undef' which means no specifically defined @@ -4892,148 +2266,78 @@ subroutine convert_vert_distrib_state(state_handle, ens_size, location, quantity if ((ztypein == ztypeout) .or. (ztypein == VERTISUNDEF)) then istatus = 0 return -else - if (debug > 9 .and. do_output()) then - write(string1,'(A,3X,2I3)') 'ztypein, ztypeout:',ztypein,ztypeout - call error_handler(E_MSG, 'convert_vert_distrib_state',string1,source, revision, revdate) - endif endif -!> assume that all locations have the same incoming lat/lon and level. -!> depending on the output vert type each member might have a different -!> vertical value. - ! unpack the incoming location(s) do e = 1, ens_size llv_loc(:, e) = get_location(location(e)) enddo -!> state_indx is i8, intent(in). -!> cellid and vert_level are integer, intent(out) call find_mpas_indices(state_indx, cellid, vert_level, ndim) -if (debug > 9 .and. do_output()) print*,'convert_vert_distrib_state: ndim=', ndim - -! the routines below will use zin as the incoming vertical value -! and zout as the new outgoing one. start out assuming failure -! (zout = missing) and wait to be pleasantly surprised when it works. zin(:) = vert_level ztypein = VERTISLEVEL + zout(:) = missing_r8 ! if the vertical is missing to start with, return it the same way ! with the requested type as out. do e = 1, ens_size if (zin(e) == missing_r8) then - location(e) = set_location(llv_loc(1, e),llv_loc(2, e),missing_r8,ztypeout) + location(e) = set_location(llv_loc(1, e), llv_loc(2, e), MISSING_R8, ztypeout) endif enddo -! if the entire ensemble has missing vertical values we can return now. -! otherwise we need to continue to convert the members with good vertical values. -! boundary cells will be updated by the assimilation. -! if all the vertical localization coord values are missing, -! we don't call this routine again, and return. -if (all(zin == missing_r8)) then ! .or. on_bound) then + +if (all(zin == MISSING_R8)) then istatus(:) = 0 return endif -! Convert the incoming vertical type (ztypein) into the vertical -! localization coordinate given in the namelist (ztypeout). -! Because this is only for state vector locations, we have already -! computed the vertical level, so all these conversions are from -! model level to something. - -! convert into: select case (ztypeout) - ! ------------------------------------------------------------ - ! outgoing vertical coordinate should be 'model level number' - ! ------------------------------------------------------------ case (VERTISLEVEL) - ! we have the vert_level and cellid - no need to call find_triangle or find_vert_indices - - zout(:) = vert_level - istatus(:) = 0 - if (debug > 9 .and. do_output()) then - write(string2,'("zout_in_level for member 1:",F10.2)') zout(1) - call error_handler(E_MSG, 'convert_vert_distrib_state',string2,source, revision, revdate) - endif + zout(:) = vert_level + istatus(:) = 0 - ! ------------------------------------------------------------ - ! outgoing vertical coordinate should be 'pressure' in Pa - ! ------------------------------------------------------------ case (VERTISPRESSURE) - !>@todo FIXME - this is the original code from the - !> observation version which does horizontal interpolation. - !> in this code we know we are on a state vector location - !> so no interp is needed. we should be able to make this - !> more computationally efficient. - - ! Need to get base offsets for the potential temperature, density, and water - ! vapor mixing fields in the state vector - ivars(1) = get_progvar_index_from_kind(QTY_POTENTIAL_TEMPERATURE) - ivars(2) = get_progvar_index_from_kind(QTY_DENSITY) - ivars(3) = get_progvar_index_from_kind(QTY_VAPOR_MIXING_RATIO) - - if (any(ivars(1:3) < 0)) then - write(string1,*) 'Internal error, cannot find one or more of: theta, rho, qv' - call error_handler(E_ERR, 'convert_vert_distrib_state',string1,source, revision, revdate) - endif - - ! Get theta, rho, qv at the interpolated location - pass in cellid we have already located - ! to save the search time. - call compute_scalar_with_barycentric (state_handle, ens_size, location(1), 3, ivars, values, istatus, cellid) - if( all(istatus /= 0) ) then - location(:) = set_location(llv_loc(1, 1),llv_loc(2, 1),missing_r8,ztypeout) - return - endif + !>@todo FIXME - this is the original code from the + !> observation version which does horizontal interpolation. + !> in this code we know we are on a state vector location + !> so no interp is needed. we should be able to make this + !> more computationally efficient. + + ivars(1) = get_varid_from_kind(anl_domid, QTY_POTENTIAL_TEMPERATURE) + ivars(2) = get_varid_from_kind(anl_domid, QTY_DENSITY) + ivars(3) = get_varid_from_kind(anl_domid, QTY_VAPOR_MIXING_RATIO) + + ! Get theta, rho, qv at the interpolated location - pass in cellid we have already located + ! to save the search time. + call compute_scalar_with_barycentric (state_handle, ens_size, location(1), 3, ivars, values, istatus, cellid) + if( all(istatus /= 0) ) then + location(:) = set_location(llv_loc(1, 1),llv_loc(2, 1),missing_r8,ztypeout) + return + endif - ! Convert theta, rho, qv into pressure - call compute_full_pressure(ens_size, values(1,:), values(2,:), values(3,:), zout, tk, istatus) - if( all(istatus /= 0) ) then - if (debug > 1 .and. do_output()) then - write(string2,'("Failed in compute_full_pressure")') - call error_handler(E_MSG, 'convert_vert_distrib_state',string2,source, revision, revdate) - endif - else - if (debug > 9 .and. do_output()) then - write(string2,'("zout[Pa] for member 1,theta,rho,qv,ier:",3F10.2,F12.8,I5)') zout(1),values(1:3,1),istatus(1) - call error_handler(E_MSG, 'convert_vert_distrib_state',string2,source, revision, revdate) - endif - endif ! istatus + ! Convert theta, rho, qv into pressure + call compute_full_pressure(ens_size, values(1,:), values(2,:), values(3,:), zout, tk, istatus) - ! ------------------------------------------------------------ - ! outgoing vertical coordinate should be 'height' in meters - ! ------------------------------------------------------------ case (VERTISHEIGHT) - ! surface obs should use the lower face of the first level. the rest - ! of the quantities should use the level centers. - if ( ndim == 1 ) then - zout(:) = zGridFace(1, cellid) - istatus(:) = 0 - else - zout(:) = zGridCenter(vert_level, cellid) - if ( quantity == QTY_VERTICAL_VELOCITY ) zout(:) = zGridFace(vert_level, cellid) - if ( quantity == QTY_EDGE_NORMAL_SPEED ) zout(:) = zGridEdge(vert_level, cellid) - istatus(:) = 0 - endif - - if (debug > 9 .and. do_output()) then - write(string2,'("zout[m] for member 1:",F10.2)') zout(1) - call error_handler(E_MSG, 'convert_vert_distrib_state',string2,source, revision, revdate) - endif - + if ( ndim == 1 ) then + zout(:) = zGridFace(1, cellid) + istatus(:) = 0 + else + zout(:) = zGridCenter(vert_level, cellid) + if ( quantity == QTY_VERTICAL_VELOCITY ) zout(:) = zGridFace(vert_level, cellid) + if ( quantity == QTY_EDGE_NORMAL_SPEED ) zout(:) = zGridEdge(vert_level, cellid) + istatus(:) = 0 + endif - ! ------------------------------------------------------------ - ! outgoing vertical coordinate should be 'scale height' (a ratio) - ! ------------------------------------------------------------ case (VERTISSCALEHEIGHT) - !>@todo FIXME - !> whatever we do for pressure, something similar here + !>@todo FIXME + !> whatever we do for pressure, something similar here ! Scale Height is defined as: log(pressure) ! if namelist item: no_normalization_of_scale_heights = .true. @@ -5055,14 +2359,14 @@ subroutine convert_vert_distrib_state(state_handle, ens_size, location, quantity if (at_surf .and. do_norm) then zout = 0.0_r8 istatus(:) = 0 - goto 101 + return endif ! Base offsets for the potential temperature, density, and water ! vapor mixing fields in the state vector. - ivars(1) = get_progvar_index_from_kind(QTY_POTENTIAL_TEMPERATURE) - ivars(2) = get_progvar_index_from_kind(QTY_DENSITY) - ivars(3) = get_progvar_index_from_kind(QTY_VAPOR_MIXING_RATIO) + ivars(1) = get_varid_from_kind(anl_domid, QTY_POTENTIAL_TEMPERATURE) + ivars(2) = get_varid_from_kind(anl_domid, QTY_DENSITY) + ivars(3) = get_varid_from_kind(anl_domid, QTY_VAPOR_MIXING_RATIO) if (at_surf .or. do_norm) then ! we will need surface pressure @@ -5076,17 +2380,6 @@ subroutine convert_vert_distrib_state(state_handle, ens_size, location, quantity ! Convert surface theta, rho, qv into pressure call compute_full_pressure(ens_size, values(1,:), values(2,:), values(3,:), surfp, tk, istatus) - if( all(istatus /= 0) ) then - if (debug > 1 .and. do_output()) then - write(string2,'("Failed in compute_full_pressure")') - call error_handler(E_MSG, 'convert_vert_distrib_state',string2,source, revision, revdate) - endif - else - if (debug > 9 .and. do_output()) then - write(string2,'("zout_psfc,theta,rho,qv:",3F10.2,F12.8)') surfp(1), values(1:3,1) - call error_handler(E_MSG, 'convert_vert_distrib_state',string2,source, revision, revdate) - endif - endif ! istatus endif @@ -5097,23 +2390,9 @@ subroutine convert_vert_distrib_state(state_handle, ens_size, location, quantity ! Convert theta, rho, qv into pressure call compute_full_pressure(ens_size, values(1,:), values(2,:), values(3,:), fullp, tk, istatus) - if( all(istatus /= 0) ) then - if (debug > 1 .and. do_output()) then - write(string2,'("Failed in compute_full_pressure")') - call error_handler(E_MSG, 'convert_vert_distrib_state',string2,source, revision, revdate) - endif - else - if (debug > 9 .and. do_output()) then - write(string2,'("zout [Pa] for member 1,theta,rho,qv:",3F10.2,F12.8)') fullp(1), values(1:3,1) - call error_handler(E_MSG, 'convert_vert_distrib_state',string2,source, revision, revdate) - endif - endif ! istatus endif - ! we have what we need now. figure out the case and set zout to the right values. - ! we've already taken care of the (at_surf .and. do_norm) case, so that simplifies - ! the tests here. if (at_surf) then where (surfp /= MISSING_R8) zout = log(surfp) @@ -5136,42 +2415,77 @@ subroutine convert_vert_distrib_state(state_handle, ens_size, location, quantity end where endif -101 continue - - if (debug > 9 .and. do_output()) then - write(string2,'("zout_in_scaleheight for member 1:",F10.2)') zout(1) - call error_handler(E_MSG, 'convert_vert_distrib_state',string2,source, revision, revdate) - endif - - ! ------------------------------------------------------- - ! outgoing vertical coordinate is unrecognized - ! ------------------------------------------------------- case default write(string1,*) 'Requested vertical coordinate not recognized: ', ztypeout - call error_handler(E_ERR,'convert_vert_distrib_state', string1, & - source, revision, revdate) + call error_handler(E_ERR,'convert_vert_state', string1, & + source) end select ! outgoing vert type -! Returned location do e = 1, ens_size if(istatus(e) == 0) then - location(e) = set_location(llv_loc(1, e),llv_loc(2, e),zout(e),ztypeout) + location(e) = set_location(llv_loc(1, e), llv_loc(2, e), zout(e), ztypeout) else - location(e) = set_location(llv_loc(1, e),llv_loc(2, e),missing_r8,ztypeout) + location(e) = set_location(llv_loc(1, e), llv_loc(2, e), MISSING_R8, ztypeout) endif enddo -end subroutine convert_vert_distrib_state +end subroutine convert_vert_state !------------------------------------------------------------------- +function required_for_rttov(qty) +integer, intent(in) :: qty +logical :: required_for_rttov -!================================================================== -! The following (private) interfaces are used for triangle interpolation -!================================================================== +select case (qty) + +!HK @todo what is a surface qty that is not required for RTTOV? + +case (QTY_SURFACE_PRESSURE, & + QTY_SURFACE_ELEVATION, & + QTY_2M_TEMPERATURE, & + QTY_2M_SPECIFIC_HUMIDITY, & + QTY_SKIN_TEMPERATURE, & + QTY_SURFACE_TYPE, & + QTY_CLOUD_FRACTION, & + QTY_10M_U_WIND_COMPONENT, & + QTY_10M_V_WIND_COMPONENT) + + required_for_rttov = .true. + +case default + + required_for_rttov = .false. + +end select + +end function required_for_rttov + +!------------------------------------------------------------------ +subroutine find_mpas_indices(index_in, cellid, vert_level, ndim) + +integer(i8), intent(in) :: index_in +integer, intent(out) :: cellid +integer, intent(out) :: vert_level +integer, intent(out) :: ndim + +integer :: i, j, k ! Indices into variable (note k is not used in MPAS) +integer :: nzp, iloc, vloc, nnf + +call get_model_variable_indices(index_in, i, j, k, var_id=nnf) +ndim = get_num_dims(anl_domid, nnf) + +if ( ndim == 2) then ! 3d (2d netcdf ) variable(vcol, iloc) + vert_level = i + cellid = j +else ! 2d (1d netcdf) variable: (cells) + cellid = i + vert_level = 1 +endif +end subroutine find_mpas_indices !------------------------------------------------------------------ !> Finds position of a given height in an array of height grid points and returns @@ -5217,9 +2531,7 @@ subroutine find_height_bounds(height, nbounds, bounds, lower, upper, fract, ier) ! too low? if(height < bounds(1)) then if(outside_grid_level_tolerance <= 0.0_r8) then - ier = 998 - if ((debug > 1) .and. do_output()) print*,'find_height_bounds: height < bounds(1) ', & - height, bounds(1) + ier = HEIGHT_TOO_LOW else ! let points outside the grid but "close enough" use the ! bottom level as the height. hlevel should come back < 1 @@ -5235,7 +2547,7 @@ subroutine find_height_bounds(height, nbounds, bounds, lower, upper, fract, ier) ! too high? if(height > bounds(nbounds)) then if(outside_grid_level_tolerance <= 0.0_r8) then - ier = 9998 + ier = HEIGHT_TOO_HIGH else ! let points outside the grid but "close enough" use the ! top level as the height. hlevel should come back > nbounds @@ -5277,7 +2589,7 @@ end subroutine find_height_bounds !> has one level, both level numbers are identical ... 1 and the !> fractions are identical ... 0.0 -subroutine find_vert_level(state_handle, ens_size, loc, nc, ids, oncenters, lower, upper, fract, ier) +subroutine find_vert_level(state_handle, ens_size, loc, nc, ids, var_id, lower, upper, fract, ier) ! note that this code handles data at cell centers, at edges, but not ! data on faces. so far we don't have any on faces. @@ -5286,18 +2598,20 @@ subroutine find_vert_level(state_handle, ens_size, loc, nc, ids, oncenters, lowe type(ensemble_type), intent(in) :: state_handle type(location_type), intent(in) :: loc integer, intent(in) :: ens_size -integer, intent(in) :: nc, ids(:) -logical, intent(in) :: oncenters -integer, intent(out) :: lower(:, :), upper(:, :) ! ens_size -real(r8), intent(out) :: fract(:, :) -integer, intent(out) :: ier(:) +integer, intent(in) :: nc, ids(nc) +integer, intent(in) :: var_id ! state variable id +integer, intent(out) :: lower(nc, ens_size), upper(nc, ens_size) +real(r8), intent(out) :: fract(nc, ens_size) +integer, intent(out) :: ier(ens_size) real(r8) :: lat, lon, vert, llv(3) real(r8) :: vert_array(ens_size) -integer :: track_ier(ens_size) -integer(i8) :: pt_base_offset, density_base_offset, qv_base_offset integer :: verttype, i integer :: e +integer :: num_levs +logical :: edgy +real(r8) :: vert_levels(nVertLevelsP1) + ! Initialization ier = 0 @@ -5318,140 +2632,93 @@ subroutine find_vert_level(state_handle, ens_size, loc, nc, ids, oncenters, lowe !is_vertical: VERTISPRESSURE !is_vertical: VERTISHEIGHT -! unpack the location into local vars -! I think you can do this with the first ensemble member? -! Because they are the same horizontally? llv = get_location(loc) lon = llv(1) lat = llv(2) vert = llv(3) verttype = nint(query_location(loc)) -! these first 3 types need no cell/edge location information. -if ((debug > 9) .and. do_output()) then - write(string2,'("vert, which_vert:",3F20.12,I5)') lon,lat,vert,verttype - call error_handler(E_MSG, 'find_vert_level',string2,source, revision, revdate) -endif +! Need - number of levels +! - edges +! - centers -! VERTISSURFACE describes variables on A surface (not necessarily THE surface) -! VERTISUNDEF describes no defined vertical location (e.g. vertically integrated vals) -if(verttype == VERTISSURFACE .or. verttype == VERTISUNDEF) then ! same across the ensemble - lower(1:nc, :) = 1 - upper(1:nc, :) = 1 - fract(1:nc, :) = 0.0_r8 - ier = 0 - return +if (get_dim_name(anl_domid, var_id, 1) == 'nVertLevels') then + num_levs = nVertLevels +else + num_levs = nVertLevelsP1 endif -! model level numbers (supports fractional levels) -if(verttype == VERTISLEVEL) then - ! FIXME: if this is W, the top is nVertLevels+1 - if (vert > nVertLevels) then - ier(:) = 81 - return - endif - - ! at the top we have to make the fraction 1.0; all other - ! integral levels can be 0.0 from the lower level. - if (vert == nVertLevels) then - lower(1:nc, :) = nint(vert) - 1 ! round down - upper(1:nc, :) = nint(vert) - fract(1:nc, :) = 1.0_r8 +! VERTISSURFACE describes variables on A surface (not necessarily THE surface) !HK @todo what does this mean? +! VERTISUNDEF describes no defined vertical location (e.g. vertically integrated vals) +select case (verttype) + case (VERTISSURFACE, VERTISUNDEF) ! same across the ensemble + lower(1:nc, :) = 1 + upper(1:nc, :) = 1 + fract(1:nc, :) = 0.0_r8 ier = 0 - return - endif - lower(1:nc, :) = aint(vert) ! round down - upper(1:nc, :) = lower+1 - fract(1:nc, :) = vert - lower - ier = 0 - return - -endif - -! ok, now we need to know where we are in the grid for heights or pressures -! as the vertical coordinate. - -! Vertical interpolation for pressure coordinates -if(verttype == VERTISPRESSURE ) then - - track_ier = 0 - vert_array = vert - - ! Need to get base offsets for the potential temperature, density, and water - ! vapor mixing fields in the state vector - call get_index_range(QTY_POTENTIAL_TEMPERATURE, pt_base_offset) - call get_index_range(QTY_DENSITY, density_base_offset) - call get_index_range(QTY_VAPOR_MIXING_RATIO, qv_base_offset) - - do i=1, nc - call find_pressure_bounds(state_handle, ens_size, vert_array, ids(i), nVertLevels, & - pt_base_offset, density_base_offset, qv_base_offset, & - lower(i, :), upper(i, :), fract(i, :), ier) + case (VERTISLEVEL) ! model level numbers (supports fractional levels) + if (vert > num_levs) then + ier(:) = VERTICAL_TOO_HIGH + return + endif - !if(debug > 9) print '(A,5I5,F10.4)', & - ! ' after find_pressure_bounds: ier, i, cellid, lower, upper, fract = ', & - ! ier, i, ids(i), lower(i, :), upper(i, :), fract(i, :) - !if(debug > 5 .and. ier(e) /= 0) print '(A,4I5,F10.4,I3,F10.4)', & - !'fail in find_pressure_bounds: ier, nc, i, id, vert, lower, fract: ', & - ! ier, nc, i, ids(i), vert_array, lower(i, :), fract(i, :) + ! at the top we have to make the fraction 1.0; all other + ! integral levels can be 0.0 from the lower level. + if (vert == num_levs) then + lower(1:nc, :) = nint(vert) - 1 ! round down + upper(1:nc, :) = nint(vert) + fract(1:nc, :) = 1.0_r8 + ier = 0 + return + endif - ! we are inside a loop over each corner. consolidate error codes - ! so that we return an error for that ensemble member if any - ! of the corners fails the pressure bounds test. - where (ier /= 0 .and. track_ier == 0) track_ier = ier + lower(1:nc, :) = aint(vert) ! round down + upper(1:nc, :) = lower+1 + fract(1:nc, :) = vert - lower + ier = 0 - enddo + case (VERTISPRESSURE ) - ier = track_ier - return -endif + vert_array(:) = vert -!HK I believe height is the same across the ensemble, so you can use -! vert and not vert_array: + do i=1, nc + call find_pressure_bounds(state_handle, ens_size, vert_array, ids(i), num_levs, & + lower(i, :), upper(i, :), fract(i, :), ier) -! grid is in height, so this needs to know which cell to index into -! for the column of heights and call the bounds routine. -if(verttype == VERTISHEIGHT) then - ! For height, can do simple vertical search for interpolation for now - ! Get the lower and upper bounds and fraction for each column - track_ier = 0 + if(fails(ier)) return - do i=1, nc - if (oncenters) then - call find_height_bounds(vert, nVertLevels, zGridCenter(:, ids(i)), & - lower(i, :), upper(i, :), fract(i, :), ier) - else + enddo + + case (VERTISHEIGHT) + ! For height, can do simple vertical search for interpolation for now + ! Get the lower and upper bounds and fraction for each column + + edgy = on_edge(anl_domid, var_id) + do i=1, nc - if (.not. data_on_edges) then - call error_handler(E_ERR, 'find_vert_level', & - 'Internal error: oncenters false but data_on_edges false', & - source, revision, revdate) + if (edgy) then + vert_levels(1:num_levs) = zGridEdge(:, ids(i)) + elseif ( num_levs == nVertLevels ) then + vert_levels(1:num_levs) = zGridCenter(:, ids(i)) + else + vert_levels(1:num_levs) = zGridFace(:, ids(i)) endif - call find_height_bounds(vert, nVertLevels, zGridEdge(:, ids(i)), & + call find_height_bounds(vert, num_levs, vert_levels(1:num_levs), & lower(i, :), upper(i, :), fract(i, :), ier) - endif - - ! we are inside a loop over each corner. consolidate error codes - ! so that we return an error for that ensemble member if any - ! of the corners fails the pressure bounds test. - where (ier /= 0 .and. track_ier == 0) track_ier = ier - - enddo - if ((debug > 9) .and. do_output()) print*,'find_vert_level is done' + if(fails(ier)) return - ier = track_ier - return + enddo -endif + end select end subroutine find_vert_level !------------------------------------------------------------------ + subroutine find_pressure_bounds(state_handle, ens_size, p, cellid, nbounds, & - pt_base_offset, density_base_offset, qv_base_offset, & lower, upper, fract, ier) ! Finds vertical interpolation indices and fraction for a quantity with @@ -5464,11 +2731,10 @@ subroutine find_pressure_bounds(state_handle, ens_size, p, cellid, nbounds, & integer, intent(in) :: ens_size real(r8), intent(in) :: p(ens_size) integer, intent(in) :: cellid -integer, intent(in) :: nbounds ! number of vertical levels? -integer(i8), intent(in) :: pt_base_offset, density_base_offset, qv_base_offset -integer, intent(out) :: lower(:), upper(:) ! ens_size -real(r8), intent(out) :: fract(:) ! ens_size -integer, intent(out) :: ier(:) ! ens_size +integer, intent(in) :: nbounds ! number of vertical levels !HK @todo why not call it num_levels? +integer, intent(out) :: lower(ens_size), upper(ens_size) +real(r8), intent(out) :: fract(ens_size) +integer, intent(out) :: ier(ens_size) integer :: i, ier2 real(r8) :: pr @@ -5488,18 +2754,12 @@ subroutine find_pressure_bounds(state_handle, ens_size, p, cellid, nbounds, & ier = 0 ! Find the lowest pressure -call get_interp_pressure(state_handle, ens_size, pt_base_offset, density_base_offset, & - qv_base_offset, cellid, 1, nbounds, pressure(1, :), temp_ier) -if(debug > 10 .and. do_output()) & - print *, 'find_pressure_bounds k=1: p, temp_ier, ier', 1, pressure(1,:), temp_ier, ier +call get_interp_pressure(state_handle, ens_size, cellid, 1, pressure(1, :), temp_ier) where(ier(:) == 0) ier(:) = temp_ier(:) ! Get the highest pressure level -call get_interp_pressure(state_handle, ens_size, pt_base_offset, density_base_offset, & - qv_base_offset, cellid, nbounds, nbounds, pressure(nbounds, :), temp_ier) -if(debug > 10 .and. do_output()) & - print *, 'find_pressure_bounds k=N: p, temp_ier, ier', nbounds, pressure(nbounds,:), temp_ier, ier +call get_interp_pressure(state_handle, ens_size, cellid, nbounds, pressure(nbounds, :), temp_ier) where(ier(:) == 0) ier(:) = temp_ier(:) @@ -5526,24 +2786,13 @@ subroutine find_pressure_bounds(state_handle, ens_size, p, cellid, nbounds, & do i = 2, nbounds ! we've already done this call for level == nbounds if (i /= nbounds) then - call get_interp_pressure(state_handle, ens_size, pt_base_offset, density_base_offset, & - qv_base_offset, cellid, i, nbounds, pressure(i, :), temp_ier) - if(debug > 11 .and. do_output()) & - print *, 'find_pressure_bounds k=?: p, temp_ier, ier ', i, pressure(i,:), temp_ier, ier - + call get_interp_pressure(state_handle, ens_size, cellid, i, pressure(i, :), temp_ier) where (ier(:) == 0) ier(:) = temp_ier(:) endif ! Check if pressure is not monotonically descreased with level. if(any(pressure(i, :) > pressure(i-1, :))) then - if (debug > 0 .and. do_output()) then - write(*, *) 'pressure at the level is larger than the level below at cellid', cellid - do e=1, ens_size - write(*, '(A,3I4,F9.2,A,F9.2)') & - 'ens#, level, level-1, p(level), p(level-1) ', e,i,i-1,pressure(i,e),' ?>? ',pressure(i-1,e) - enddo - endif - where(pressure(i, :) > pressure(i-1, :)) ier(:) = 988 + where(pressure(i, :) > pressure(i-1, :)) ier(:) = PRESSURE_NOT_MONOTONIC endif ! each ensemble member could have a vertical between different levels, @@ -5570,11 +2819,6 @@ subroutine find_pressure_bounds(state_handle, ens_size, p, cellid, nbounds, & else fract(e) = (p(e) - pressure(i-1,e)) / (pressure(i,e) - pressure(i-1,e)) endif - - if ((debug > 9) .and. do_output()) print '(A,3F10.2,2I4,F10.2)', & - "find_pressure_bounds: p_in, pr(i-1), pr(i), lower, upper, fract = ", & - p(e), pressure(i-1,e), pressure(i,e), lower(e), upper(e), fract(e) - endif enddo if(all(found_level)) return @@ -5584,35 +2828,36 @@ end subroutine find_pressure_bounds !------------------------------------------------------------------ -subroutine get_interp_pressure(state_handle, ens_size, pt_offset, density_offset, qv_offset, & - cellid, lev, nlevs, pressure, ier) +subroutine get_interp_pressure(state_handle, ens_size, cellid, lev, pressure, ier) ! Finds the value of pressure at a given point at model level lev type(ensemble_type), intent(in) :: state_handle integer, intent(in) :: ens_size -integer(i8), intent(in) :: pt_offset, density_offset, qv_offset integer, intent(in) :: cellid -integer, intent(in) :: lev, nlevs -real(r8), intent(out) :: pressure(:) -integer, intent(out) :: ier(:) +integer, intent(in) :: lev +real(r8), intent(out) :: pressure(ens_size) +integer, intent(out) :: ier(ens_size) -integer :: offset +integer(i8) :: pt_idx, density_idx, qv_idx +integer :: dummyk real(r8) :: pt(ens_size), density(ens_size), qv(ens_size), tk(ens_size) integer :: e -if(debug > 9 .and. do_output()) print*, 'get_interp_pressure for k = ',lev - ! Get the values of potential temperature, density, and vapor -offset = (cellid - 1) * nlevs + lev - 1 -pt = get_state(pt_offset + offset, state_handle) -density = get_state(density_offset + offset, state_handle) -qv = get_state(qv_offset + offset, state_handle) +pt_idx = get_dart_vector_index(lev, cellid, dummyk, anl_domid, get_varid_from_kind(anl_domid, QTY_POTENTIAL_TEMPERATURE)) +density_idx = get_dart_vector_index(lev, cellid, dummyk, anl_domid, get_varid_from_kind(anl_domid, QTY_DENSITY)) +qv_idx = get_dart_vector_index(lev, cellid, dummyk, anl_domid, get_varid_from_kind(anl_domid, QTY_VAPOR_MIXING_RATIO)) + +pt = get_state(pt_idx, state_handle) +density = get_state(density_idx, state_handle) +qv = get_state(qv_idx, state_handle) ! Initialization ier = 0 ! Error if any of the values are missing; probably will be all or nothing +!HK @todo why would the state be missing? do e = 1, ens_size if(pt(e) == MISSING_R8 .or. density(e) == MISSING_R8 .or. qv(e) == MISSING_R8) then ier(e) = 2 @@ -5622,17 +2867,6 @@ subroutine get_interp_pressure(state_handle, ens_size, pt_offset, density_offset ! Convert theta, rho, qv into pressure call compute_full_pressure(ens_size, pt(:), density(:), qv(:), pressure(:), tk(:), ier) !HK where is tk used? -if( all(ier/= 0) ) then - if (debug > 1 .and. do_output()) then - write(string2,'("Failed in compute_full_pressure")') - call error_handler(E_ALLMSG, 'get_interp_pressure',string2,source, revision, revdate) - endif -else - if((debug > 11) .and. do_output()) then - write(*,*) 'get_interp_pressure: cellid, lev', cellid, lev - write(*,*) 'get_interp_pressure: theta,rho,qv,p,tk', pt(1), density(1), qv(1), pressure(1), tk(1) - endif -endif !( all(istatus /= 0) ) end subroutine get_interp_pressure @@ -5732,7 +2966,7 @@ end subroutine get_barycentric_weights !------------------------------------------------------------ !> this routine computes 1 or more values at a single location, !> for each ensemble member. "n" is the number of different -!> quantities to compute, ival is an array of 'n' progval() indices +!> quantities to compute, ival is an array of 'n' state indices !> to indicate which quantities to compute. !> dval(n, ens_size) are the output values, and ier(ens_size) are !> the success/error returns for each ensemble member. @@ -5741,7 +2975,7 @@ subroutine compute_scalar_with_barycentric(state_handle, ens_size, loc, n, ival, type(ensemble_type), intent(in) :: state_handle integer, intent(in) :: ens_size -type(location_type), intent(in) :: loc !(ens_size) +type(location_type), intent(in) :: loc integer, intent(in) :: n integer, intent(in) :: ival(n) real(r8), intent(out) :: dval(n, ens_size) @@ -5752,92 +2986,69 @@ subroutine compute_scalar_with_barycentric(state_handle, ens_size, loc, n, ival, real(r8), dimension(ens_size) :: lowval, uppval real(r8) :: weights(3) integer :: c(3), nvert, k, i, nc -integer(i8) :: index1, low_offset, upp_offset, low_state_indx, upp_state_indx +integer(i8) :: low_state_indx(ens_size), upp_state_indx(ens_size) +integer :: jdummy, kdummy integer :: lower(3, ens_size), upper(3,ens_size) integer :: e, e2, thislower, thisupper logical :: did_member(ens_size) -! assume failure dval = MISSING_R8 -ier = 88 ! field not in state vector +ier = QTY_NOT_IN_STATE_VECTOR ! make sure we have all good field indices first if (any(ival < 0)) return -call find_triangle (loc, nc, c, weights, ier(1), this_cellid) +call find_triangle (loc, nc, c, weights, ier(1), this_cellid) !HK c is cell centers if(ier(1) /= 0) then ier(:) = ier(1) return endif ! If the field is on a single level, lower and upper are both 1 -call find_vert_indices (state_handle, ens_size, loc, nc, c, lower, upper, fract, ier) -if(all(ier /= 0)) return +! HK @todo this is assuming the array qtys are all on the same zgrid. +! ival(1) to find levels for array of qtys +call find_vert_level(state_handle, ens_size, loc, nc, c, ival(1), lower, upper, fract, ier) -! for each field to compute at this location: -do k=1, n - ! get the starting index in the state vector - index1 = progvar(ival(k))%index1 - nvert = progvar(ival(k))%numvertical +if(fails(ier)) return - ! for each corner: could be 1 if location is directly at a vertex, but - ! normally is 3 for the enclosing triangle made up of cell centers. - do i = 1, nc - ! go around triangle and interpolate in the vertical - ! c(3) are the cell ids +! for each field to compute at this location: !HK @todo why loop in here? why not outside? +do k = 1, n - low_offset = (c(i)-1) * nvert !FIXME low_offset and upp_offset are the same? - upp_offset = (c(i)-1) * nvert - - if( nvert == 1 ) then ! fields on a surface (1-D, one level, ...) - - ! Because 'lower(i,1)-1' evaluates to zero for 1-D, no need to add it - fdata(i,:) = get_state(index1 + low_offset, state_handle) + if( get_num_dims(anl_domid, ival(k)) == 1 ) then ! 1-D field + do i = 1, nc + ! go around triangle and interpolate in the vertical + ! c(3) are the cell ids + low_state_indx(1) = get_dart_vector_index(c(i), jdummy, kdummy, anl_domid, ival(k)) + fdata(i,:) = get_state(low_state_indx(1), state_handle) + enddo - else ! 2-D fields + else ! 2-D field - did_member(:) = .false. + do i = 1, nc members: do e = 1, ens_size - if (did_member(e)) cycle members - - !> minimize the number of times we call get_state() by - !> doing all the ensemble members which are between the same - !> two vertical levels. this is true most of the time. - !> in some cases it could be 2 or 3 different pairs of levels because - !> of differences in vertical conversion that depends on per-member fields. - low_state_indx = index1 + low_offset + lower(i,e)-1 - lowval(:) = get_state(low_state_indx, state_handle) - upp_state_indx = index1 + upp_offset + upper(i,e)-1 - uppval(:) = get_state(upp_state_indx, state_handle) - - thislower = lower(i, e) - thisupper = upper(i, e) - - ! for all remaining ensemble members, use these values if the lower and - ! upper level numbers are the same. fract() will vary with member. - do e2=e, ens_size - if (thislower == lower(i, e2) .and. thisupper == upper(i, e2)) then - fdata(i, e2) = lowval(e2)*(1.0_r8 - fract(i, e2)) + uppval(e2)*fract(i, e2) - did_member(e2) = .true. - endif - enddo + low_state_indx(e) = get_dart_vector_index(lower(i,e), c(i), kdummy, anl_domid, ival(k)) + upp_state_indx(e) = get_dart_vector_index(upper(i,e), c(i), kdummy, anl_domid, ival(k)) + enddo members - endif ! 2-D fields + call get_state_array(lowval(:), low_state_indx, state_handle) + call get_state_array(uppval(:), upp_state_indx, state_handle) - enddo ! corners + fdata(i, :) = lowval(:)*(1.0_r8 - fract(i, :)) + uppval(:)*fract(i,:) + enddo ! corners + + endif ! now have vertically interpolated values at cell centers. ! use weights to compute value at interp point. do e = 1, ens_size - if (ier(e) /= 0) cycle - dval(k, e) = sum(weights(1:nc) * fdata(1:nc, e)) + dval(k, e) = sum(weights(1:nc) * fdata(1:nc, e)) enddo - enddo + end subroutine compute_scalar_with_barycentric !------------------------------------------------------------ @@ -5847,7 +3058,7 @@ end subroutine compute_scalar_with_barycentric subroutine compute_surface_data_with_barycentric(var1d, loc, dval, ier, this_cellid) type(location_type), intent(in) :: loc -real(r8), intent(in) :: var1d(:) +real(r8), intent(in) :: var1d(:) ! whole static variable. HK why pass this in? real(r8), intent(out) :: dval integer, intent(out) :: ier integer, optional, intent(in) :: this_cellid @@ -5868,10 +3079,6 @@ subroutine compute_surface_data_with_barycentric(var1d, loc, dval, ier, this_cel ! use weights to compute value at interp point. dval = sum(weights(1:nc) * fdata(1:nc)) -if(debug > 9 .and. do_output()) & - print '(A,7f12.5)','compute_surface_data_with_barycentric: corner vals, weights, result: ',& - fdata(:),weights(:),dval - end subroutine compute_surface_data_with_barycentric !------------------------------------------------------------ @@ -5920,11 +3127,7 @@ subroutine find_triangle(loc, nc, c, weights, ier, this_cellid) cellid = find_closest_cell_center(lat, lon) endif -if ((xyzdebug > 1) .and. do_output()) & - print *, 'ft: closest cell center for lon/lat: ', lon, lat, cellid - if (cellid < 1) then - if(xyzdebug > 0) print *, 'ft: closest cell center for lon/lat: ', lon, lat, cellid ier = 11 return endif @@ -5933,8 +3136,6 @@ subroutine find_triangle(loc, nc, c, weights, ier, this_cellid) ! closest vertex to given point. closest_vert = closest_vertex_ll(cellid, lat, lon) -if ((xyzdebug > 5) .and. do_output()) & - print *, 'ft: closest vertex for lon/lat: ', lon, lat, closest_vert ! collect the neighboring cell ids and vertex numbers ! this 2-step process avoids us having to read in the @@ -5946,12 +3147,10 @@ subroutine find_triangle(loc, nc, c, weights, ier, this_cellid) vindex = 1 nedges = nEdgesOnCell(cellid) do i=1, nedges -!print *, 'ft: i: ', i edgeid = edgesOnCell(i, cellid) -!print *, 'ft: edgeid: ', edgeid if (.not. global_grid .and. & (cellsOnEdge(1, edgeid) <= 0 .or. cellsOnEdge(2, edgeid) <= 0)) then - ier = 14 + ier = TRIANGLE_CELL_CENTER_NOT_FOUND return endif if (cellsOnEdge(1, edgeid) /= cellid) then @@ -5960,7 +3159,6 @@ subroutine find_triangle(loc, nc, c, weights, ier, this_cellid) neighborcells(i) = cellsOnEdge(2, edgeid) endif verts(i) = verticesOnCell(i, cellid) -!print *, 'ft: verts: ', verts(i), closest_vert if (verts(i) == closest_vert) vindex = i call latlon_to_xyz(latCell(neighborcells(i)), lonCell(neighborcells(i)), & xdata(i), ydata(i), zdata(i)) @@ -5972,7 +3170,6 @@ subroutine find_triangle(loc, nc, c, weights, ier, this_cellid) ! and the observation point call latlon_to_xyz(lat, lon, r(1), r(2), r(3)) -!print *, 'ft: lat/lon: ', lat, lon if (all(abs(t1-r) < roundoff)) then ! Located at a grid point (counting roundoff errors) @@ -6010,65 +3207,14 @@ subroutine find_triangle(loc, nc, c, weights, ier, this_cellid) endif enddo findtri if (.not. foundit) then - ier = 14 ! 11? + ier = TRIANGLE_CELL_CENTER_NOT_FOUND return endif endif ! horizontal index search is done now. -if (ier /= 0) return - -if (debug > 12 .and. do_output()) then - write(string3,*) 'ier = ',ier, ' triangle = ',c(1:nc), ' weights = ',weights(1:nc) - call error_handler(E_MSG, 'find_triangle', string3, source, revision, revdate) -endif end subroutine find_triangle -!------------------------------------------------------------ -!> This routine adds some error checking because find_vert_level -!> has some early return statements that prevent having the debugging -!> information in it. - -subroutine find_vert_indices (state_handle, ens_size, loc, nc, c, lower, upper, fract, ier) - -type(ensemble_type), intent(in) :: state_handle -type(location_type), intent(in) :: loc -integer, intent(in) :: ens_size -integer, intent(in) :: nc -integer, intent(in) :: c(:) -integer, intent(out) :: lower(:, :), upper(:, :) ! ens_size -real(r8), intent(out) :: fract(:, :) ! ens_size -integer, intent(out) :: ier(:) ! ens_size - -integer :: e - -! initialization -lower = MISSING_I -upper = MISSING_I -fract = 0.0_r8 - -! need vert index for the vertical level -call find_vert_level(state_handle, ens_size, loc, nc, c, .true., lower, upper, fract, ier) - -if (debug > 9 .and. do_output()) then - write(string3,*) 'ier = ',ier(1), ' triangle = ',c(1:nc), ' vert_index = ',lower(1:nc, 1)+fract(1:nc, 1), ' nc = ', nc - call error_handler(E_MSG, 'find_vert_indices', string3, source, revision, revdate) -endif - -if (debug > 12 .and. do_output()) then - do e = 1, ens_size - if(ier(e) /= 0) then - print *, 'find_vert_indices: e = ', e, ' nc = ', nc, ' ier = ', ier(e) - print *, 'find_vert_indices: c = ', c - print *, 'find_vert_indices: lower = ', lower(1:nc, e) - print *, 'find_vert_indices: upper = ', upper(1:nc, e) - print *, 'find_vert_indices: fract = ', fract(1:nc, e) - endif - enddo -endif - -end subroutine find_vert_indices - !------------------------------------------------------------ subroutine compute_u_with_rbf(state_handle, ens_size, loc, zonal, uval, ier) @@ -6077,8 +3223,8 @@ subroutine compute_u_with_rbf(state_handle, ens_size, loc, zonal, uval, ier) integer, intent(in) :: ens_size type(location_type), intent(in) :: loc logical, intent(in) :: zonal -real(r8), intent(out) :: uval(:) ! ens_size -integer, intent(out) :: ier(:) ! ens_size +real(r8), intent(out) :: uval(ens_size) +integer, intent(out) :: ier(ens_size) ! max edges we currently use is ~50, but overestimate for now integer, parameter :: listsize = 200 @@ -6092,34 +3238,25 @@ subroutine compute_u_with_rbf(state_handle, ens_size, loc, zonal, uval, ier) real(r8) :: ureconstructzonal, ureconstructmeridional real(r8) :: datatangentplane(3,2) real(r8) :: coeffs_reconstruct(3,listsize) -integer(i8) :: upindx, lowindx -integer :: index1, progindex, cellid, vertexid +integer(i8) :: upindx(ens_size), lowindx(ens_size) +integer :: cellid, vertexid real(r8) :: lat, lon, vert, llv(3), fract(listsize, ens_size), lowval(ens_size), uppval(ens_size) integer :: verttype, lower(listsize, ens_size), upper(listsize, ens_size), ncells, celllist(listsize) +integer :: var_id, dummy integer :: e ! loop index -! Initialization ier = 0 uval = MISSING_R8 -! FIXME: it would be great to make this cache the last value and -! if the location is the same as before and it's asking for V now -! instead of U, skip the expensive computation. however, given -! how we currently distribute observations the V wind obs will -! almost certainly be given to a different task. if that changes -! in some future version of dart, revisit this code as well. - -progindex = get_index_from_varname('u') -if (progindex < 0 .or. .not. data_on_edges) then - ! cannot compute u if it isn't in the state vector, or if we - ! haven't read in the edge data (which shouldn't happen if +var_id = get_varid_from_varname(anl_domid, 'u') +if (var_id < 0 .or. .not. data_on_edges) then + ! cannot compute u if it is not in the state vector, or if we + ! have not read in the edge data (which shouldn't happen if ! u is in the state vector. - ier = 18 + ier = RBF_U_COMPUTATION_ERROR return endif -index1 = progvar(progindex)%index1 -nvert = progvar(progindex)%numvertical ! unpack the location into local vars llv = get_location(loc) ! I believe this is the same across the ensemble @@ -6129,33 +3266,32 @@ subroutine compute_u_with_rbf(state_handle, ens_size, loc, zonal, uval, ier) verttype = nint(query_location(loc)) call find_surrounding_edges(lat, lon, nedges, edgelist, cellid, vertexid) -if (nedges <= 0) then - ! we are on a boundary, no interpolation - ier = 18 - return -endif +if (nedges <= 0) then ! we are on a boundary, no interpolation + ier = 18 + return + endif if (verttype == VERTISPRESSURE) then ! get all cells which share any edges with the 3 cells which share ! the closest vertex. call make_cell_list(vertexid, 3, ncells, celllist) - call find_vert_level(state_handle, ens_size, loc, ncells, celllist, .true., & - lower, upper, fract, ier) + call find_vert_level(state_handle, ens_size, loc, ncells, celllist, -1, & !HK @todo this varid needs to signify pressure + lower, upper, fract, ier) - if (all(ier /= 0)) return + if (fails(ier)) return ! now have pressure at all cell centers - need to interp to get pressure ! at edge centers. call move_pressure_to_edges(ncells, celllist, lower, upper, fract, & - nedges, edgelist, ier) - if (all(ier /= 0)) return + nedges, edgelist, ier) + if (fails(ier)) return -else +else !HK @todo is the only other option VERTISHEIGHT? ! need vert index for the vertical level - call find_vert_level(state_handle, ens_size, loc, nedges, edgelist, .false., & - lower, upper, fract, ier) - if (all(ier /= 0)) return + call find_vert_level(state_handle, ens_size, loc, nedges, edgelist, var_id, & + lower, upper, fract, ier) + if (fails(ier)) return endif ! the rbf code needs (their names == our names): @@ -6169,11 +3305,12 @@ subroutine compute_u_with_rbf(state_handle, ens_size, loc, zonal, uval, ier) ! ryan has size(xEdge, 1) but this is a 1d array, so a dimension shouldn't be needed ! was possibly required for ifort v17 on cheyenne if (edgelist(i) > size(xEdge,1)) then - write(string1, *) 'edgelist has index larger than edge count', & - i, edgelist(i), size(xEdge,1) + write(string1, *) 'edgelist has index larger than edge count', & + i, edgelist(i), size(xEdge,1) call error_handler(E_ERR, 'compute_u_with_rbf', 'internal error', & - source, revision, revdate, text2=string1) + source, text2=string1) endif + xdata(i) = xEdge(edgelist(i)) ydata(i) = yEdge(edgelist(i)) zdata(i) = zEdge(edgelist(i)) @@ -6182,21 +3319,19 @@ subroutine compute_u_with_rbf(state_handle, ens_size, loc, zonal, uval, ier) edgenormals(j, i) = edgeNormalVectors(j, edgelist(i)) enddo - !>@todo lower (upper) could be different levels in pressure - !lowval = x(index1 + (edgelist(i)-1) * nvert + lower(i)-1) - lowindx = int(index1,i8) + int((edgelist(i)-1) * nvert,i8) + int(lower(i,1)-1,i8) - lowval = get_state(lowindx, state_handle) - - !uppval = x(index1 + (edgelist(i)-1) * nvert + upper(i)-1) - upindx = int(index1,i8) + int((edgelist(i)-1) * nvert,i8) + int(upper(i,1)-1,i8) - uppval = get_state(upindx, state_handle) + ! Lower/upper could be different levels in pressure + do e = 1, ens_size + lowindx = get_dart_vector_index(lower(i, e), edgelist(i), dummy,anl_domid, var_id) !HK @todo check x,y,z + upindx = get_dart_vector_index(upper(i, e), edgelist(i), dummy,anl_domid, var_id) !HK @todo check x,y,z + enddo + call get_state_array(lowval, lowindx, state_handle) + call get_state_array(uppval, upindx, state_handle) veldata(i, :) = lowval*(1.0_r8 - fract(i, :)) + uppval*fract(i, :) enddo - ! this intersection routine assumes the points are on the surface ! of a sphere. they do NOT try to make a plane between the three ! points. the difference is small, granted. @@ -6204,8 +3339,8 @@ subroutine compute_u_with_rbf(state_handle, ens_size, loc, zonal, uval, ier) ! call a simple subroutine to define vectors in the tangent plane call get_geometry(nedges, xdata, ydata, zdata, & - xreconstruct, yreconstruct, zreconstruct, edgenormals, & - on_a_sphere, datatangentplane) + xreconstruct, yreconstruct, zreconstruct, edgenormals, & + on_a_sphere, datatangentplane) ! calculate coeffs_reconstruct call get_reconstruct_init(nedges, xdata, ydata, zdata, & @@ -6215,14 +3350,14 @@ subroutine compute_u_with_rbf(state_handle, ens_size, loc, zonal, uval, ier) ! do the reconstruction do e = 1, ens_size call get_reconstruct(nedges, lat*deg2rad, lon*deg2rad, & - coeffs_reconstruct, on_a_sphere, veldata(:, e), & - ureconstructx, ureconstructy, ureconstructz, & - ureconstructzonal, ureconstructmeridional) + coeffs_reconstruct, on_a_sphere, veldata(:, e), & + ureconstructx, ureconstructy, ureconstructz, & + ureconstructzonal, ureconstructmeridional) if (zonal) then - uval = ureconstructzonal + uval(e) = ureconstructzonal else - uval = ureconstructmeridional + uval(e) = ureconstructmeridional endif enddo @@ -6331,7 +3466,7 @@ subroutine find_surrounding_edges(lat, lon, nedges, edge_list, cellid, vertexid) case default call error_handler(E_ERR, 'find_surrounding_edges', 'bad use_rbf_option value', & - source, revision, revdate) + source) end select ! Check if any of edges are located in the boundary zone. @@ -6339,9 +3474,6 @@ subroutine find_surrounding_edges(lat, lon, nedges, edge_list, cellid, vertexid) if (on_boundary_edgelist(edge_list)) then nedges = -1 edge_list(:) = -1 - if(debug > 0 .and. do_output()) call error_handler(E_MSG, 'find_surrounding_edges', 'edges in the boundary', & - source, revision, revdate) - return endif end subroutine find_surrounding_edges @@ -6401,15 +3533,10 @@ function find_closest_cell_center(lat, lon) !> happen in the regional case. if (rc /= 0 .or. closest_cell < 0) then - if (debug > 0 .or. global_grid) then - print *, 'cannot find nearest cell to lon, lat: ', lon, lat - endif find_closest_cell_center = -1 return endif -if (debug > 9 .and. do_output()) print *, 'find_closest_cell_center: lat/lon closest to cellid, with lat/lonCell: ', & - lat, lon, closest_cell, latCell(closest_cell), lonCell(closest_cell) ! do allow boundary cells to be returned. find_closest_cell_center = closest_cell @@ -6468,20 +3595,10 @@ subroutine set_global_grid(ncid) else string1 = 'MPAS running in regional (limited-area) mode' endif -call error_handler(E_MSG,'set_global_grid',string1,source,revision,revdate) +call error_handler(E_MSG,'set_global_grid',string1,source) end subroutine set_global_grid -!------------------------------------------------------------ -!> accessor function for global_grid module variable - -function is_global_grid() -logical :: is_global_grid - -is_global_grid = global_grid - -end function is_global_grid - !------------------------------------------------------------ !> find the closest cell center. if global grid, return it. !> if regional grid, continue to be sure that all three corners @@ -6576,6 +3693,23 @@ function on_boundary_edge(edgeid) end function on_boundary_edge +!------------------------------------------------------------ +!> Determine if this edge is on the boundary + +function cell_next_to_boundary_edge(edgeid) + +integer, intent(in) :: edgeid +logical :: cell_next_to_boundary_edge + +if (global_grid .or. .not. allocated(bdyMaskEdge)) return + +cell_next_to_boundary_edge = .false. +cell_next_to_boundary_edge = on_boundary_cell(cellsOnEdge(1,edgeid)) +cell_next_to_boundary_edge = on_boundary_cell(cellsOnEdge(2,edgeid)) + +end function cell_next_to_boundary_edge + + !------------------------------------------------------------ !> Determine if this edge is the outermost edge (bdyMaskEdge = 7) @@ -6685,9 +3819,6 @@ function closest_vertex_ll(cellid, lat, lon) call latlon_to_xyz(lat, lon, px, py, pz) closest_vertex_ll = closest_vertex_xyz(cellid, px, py, pz) -if ((closest_vertex_ll < 0) .and. & - (debug > 0) .and. do_output()) & - print *, 'cannot find nearest vertex to lon, lat: ', lon, lat end function closest_vertex_ll @@ -6893,7 +4024,7 @@ subroutine make_cell_list(vertexid, degree, ncells, cell_list) if (degree < 1 .or. degree > 3) then call error_handler(E_ERR, 'make_cell_list', 'internal error, must have 1 <= degree <= 3', & - source, revision, revdate) + source) endif ! use the cellsOnVertex() array to find the three cells @@ -7006,7 +4137,7 @@ subroutine move_pressure_to_edges(ncells, celllist, lower, upper, fract, nedges, if (x1 < 0.0_r8 .or. x2 < 0.0_r8) then write(string1, *) 'neither cell found in celllist', c1, c2, ' edge is ', edgelist(i) call error_handler(E_ERR, 'move_pressure_to_edges', 'cell id not found in list', & - source, revision, revdate, text2=string1) + source, text2=string1) endif ! FIXME: should this be a log interpolation since we know this is going ! to be used for pressure only? it's in the horizontal so the values @@ -7055,10 +4186,9 @@ function extrapolate_level(h, lb, ub) write(string1, *) h, ' not outside range of ', lb, ' and ', ub call error_handler(E_ERR, 'extrapolate_level', & 'extrapolate code called for height not out of range', & - source, revision, revdate, text2=string1) + source, text2=string1) endif -if (debug > 1 .and. do_output()) print *, 'extrapolate: ', h, lb, ub, extrapolate_level end function extrapolate_level !------------------------------------------------------------ @@ -7088,32 +4218,6 @@ end subroutine latlon_to_xyz !------------------------------------------------------------ -subroutine xyz_to_latlon(x, y, z, lat, lon) - -! Given a cartesian x, y, z coordinate relative to the origin -! at the center of the earth, using a fixed radius specified -! by MPAS (in the grid generation step), return the corresponding -! lat, lon location in degrees. - -real(r8), intent(in) :: x, y, z -real(r8), intent(out) :: lat, lon - -real(r8) :: rlat, rlon - -! right now this is only needed for debugging messages. -! the arc versions of routines are expensive. - -rlat = PI/2.0_r8 - acos(z/radius) -rlon = atan2(y,x) -if (rlon < 0) rlon = rlon + PI*2 - -lat = rlat * rad2deg -lon = rlon * rad2deg - -end subroutine xyz_to_latlon - -!------------------------------------------------------------ - subroutine inside_triangle(t1, t2, t3, r, lat, lon, inside, weights) ! given 3 corners of a triangle and an xyz point, compute whether @@ -7166,134 +4270,68 @@ subroutine inside_triangle(t1, t2, t3, r, lat, lon, inside, weights) end subroutine inside_triangle !------------------------------------------------------------ +function uv_increments_cell_to_edges(zonal_wind, meridional_wind) result(u_inc) -function vector_magnitude(a) - -! Given a cartesian vector, compute the magnitude - -real(r8), intent(in) :: a(3) -real(r8) :: vector_magnitude - -vector_magnitude = sqrt(a(1)*a(1) + a(2)*a(2) + a(3)*a(3)) - -end function vector_magnitude - -!------------------------------------------------------------ - -subroutine vector_cross_product(a, b, r) - -! Given 2 cartesian vectors, compute the cross product of a x b - -real(r8), intent(in) :: a(3), b(3) -real(r8), intent(out) :: r(3) - -r(1) = a(2)*b(3) - a(3)*b(2) -r(2) = a(3)*b(1) - a(1)*b(3) -r(3) = a(1)*b(2) - a(2)*b(1) - -end subroutine vector_cross_product - -!------------------------------------------------------------ - -function vector_dot_product(a, b) - -! Given 2 cartesian vectors, compute the dot product of a . b - -real(r8), intent(in) :: a(3), b(3) -real(r8) :: vector_dot_product - -vector_dot_product = a(1)*b(1) + a(2)*b(2) + a(3)*b(3) - -end function vector_dot_product - -!------------------------------------------------------------ - -subroutine vector_projection(a, b, r) - -! Given 2 cartesian vectors, project a onto b - -real(r8), intent(in) :: a(3), b(3) -real(r8), intent(out) :: r(3) - -real(r8) :: ab_over_bb - -ab_over_bb = vector_dot_product(a, b) / vector_dot_product(b, b) -r = (ab_over_bb) * b - -end subroutine vector_projection - -!------------------------------------------------------------ - -subroutine determinant3(a, r) - -! Given a 3x3 matrix, compute the determinant - -real(r8), intent(in) :: a(3,3) -real(r8), intent(out) :: r - -r = a(1,1)*(a(2,2)*a(3,3) - (a(3,2)*a(2,3))) + & - a(2,1)*(a(3,2)*a(1,3) - (a(3,3)*a(1,2))) + & - a(3,1)*(a(1,2)*a(2,3) - (a(2,2)*a(1,3))) +! Project u, v wind increments at cell centers onto the edges. +! Increments at the outermost edge are set to 0.0 because we do not update/compute +! uedge in the outermost edge in the regional MPAS. -end subroutine determinant3 +real(r8), intent(in) :: zonal_wind(:,:) ! u wind increments from filter +real(r8), intent(in) :: meridional_wind(:,:) ! v wind increments from filter +real(r8) :: u_inc(size(zonal_wind, 1), size(zonal_wind, 2)) ! normal velocity increment on the edges -!------------------------------------------------------------ +real(r8) :: uedge(size(zonal_wind, 1), size(zonal_wind, 2)) -subroutine invert3(a, r) +if ( .not. module_initialized ) call static_init_model -! Given a 3x3 matrix, compute the inverse +uedge(:,:) = 0.0_r8 ! initialize increments to 0 for outermost edge +call project_uv_cell_to_edges(zonal_wind, meridional_wind, uedge) +u_inc = uedge -real(r8), intent(in) :: a(3,3) -real(r8), intent(out) :: r(3,3) +end function uv_increments_cell_to_edges -real(r8) :: det, b(3,3) +!------------------------------------------------------------------ -call determinant3(a, det) -if (det == 0.0_r8) then - print *, 'matrix cannot be inverted' - r = 0.0_r8 - return -endif +subroutine uv_field_cell_to_edges(zonal_wind, meridional_wind, uedge) -b(1,1) = a(2,2)*a(3,3) - a(3,2)*a(2,3) -b(2,1) = a(3,1)*a(2,3) - a(2,1)*a(3,3) -b(3,1) = a(2,1)*a(3,2) - a(3,1)*a(2,2) +! Replaces a edge normal velocity field, by projecting u, v wind fields at cell +! centers onto the edges. -b(1,2) = a(3,2)*a(1,3) - a(1,2)*a(3,3) -b(2,2) = a(1,1)*a(3,3) - a(3,1)*a(1,3) -b(3,2) = a(3,1)*a(1,2) - a(1,1)*a(3,2) +real(r8), intent(in) :: zonal_wind(:,:) ! u wind from filter +real(r8), intent(in) :: meridional_wind(:,:) ! v wind from filter +real(r8), intent(inout):: uedge(:,:) ! normal velocity on the edges -b(1,3) = a(1,2)*a(2,3) - a(2,2)*a(1,3) -b(2,3) = a(1,3)*a(2,1) - a(1,1)*a(2,3) -b(3,3) = a(1,1)*a(2,2) - a(2,1)*a(1,2) +integer, parameter :: R3 = 3 +real(r8) :: east(R3,nCells), north(R3,nCells) +real(r8) :: lonCell_rad(nCells), latCell_rad(nCells) +integer :: iCell, iEdge, cell1, cell2 -r = b / det +if ( .not. module_initialized ) call static_init_model -end subroutine invert3 +call project_uv_cell_to_edges(zonal_wind, meridional_wind, uedge) -!------------------------------------------------------------ +end subroutine uv_field_cell_to_edges + +!------------------------------------------------------------------ -subroutine uv_cell_to_edges(zonal_wind, meridional_wind, uedge, full_u) +subroutine project_uv_cell_to_edges(zonal_wind, meridional_wind, uedge) -! Project u, v wind (increments) at cell centers onto the edges. +! Replaces a edge normal velocity , by projecting u, v wind at cell +! centers onto the edges. +! The original field on the outermost edge is left unchanged, because +! we do not update/compute uedge in the outermost edge in the regional MPAS. +! ! FIXME: ! we can hard-code R3 here since it comes from the (3d) x/y/z cartesian coordinate. ! We read cellsOnEdge and edgeNormalVectors in get_grid. -! Here "uedge" is an edge normal wind which can be either a full field or an increment -! depending on the optional input (full_u). -! if full_u = .true., then the edge wind is replaced by averaging zonal and meridional -! winds at cell centers. -! if full_u does not exist, uedge is the analysis increment from the updated cell-center winds. -! We do not update/compute uedge in the outermost edge in the regional MPAS. +! ! This routine followed the updating part in tend_toEdges in ! MPAS/src/core_atmosphere/physics/mpas_atmphys_todynamics.F. -real(r8), intent(in) :: zonal_wind(:,:) ! u wind updated from filter -real(r8), intent(in) :: meridional_wind(:,:) ! v wind updated from filter -real(r8), intent(inout):: uedge(:,:) ! normal velocity (increment) on the edges -logical, intent(in), optional :: full_u ! compute a full field, not an increment +real(r8), intent(in) :: zonal_wind(:,:) ! u wind from filter +real(r8), intent(in) :: meridional_wind(:,:) ! v wind from filter +real(r8), intent(inout):: uedge(:,:) ! normal velocity on the edges -! Local variables integer, parameter :: R3 = 3 real(r8) :: east(R3,nCells), north(R3,nCells) real(r8) :: lonCell_rad(nCells), latCell_rad(nCells) @@ -7301,11 +4339,6 @@ subroutine uv_cell_to_edges(zonal_wind, meridional_wind, uedge, full_u) if ( .not. module_initialized ) call static_init_model -! Initialization for increments -if ( .not. present(full_u)) then - uedge(:,:) = 0.0_r8 -endif - ! Back to radians (locally) lonCell_rad = lonCell*deg2rad latCell_rad = latCell*deg2rad @@ -7313,15 +4346,15 @@ subroutine uv_cell_to_edges(zonal_wind, meridional_wind, uedge, full_u) ! Compute unit vectors in east and north directions for each cell: do iCell = 1, nCells - east(1,iCell) = -sin(lonCell_rad(iCell)) - east(2,iCell) = cos(lonCell_rad(iCell)) - east(3,iCell) = 0.0_r8 - call r3_normalize(east(1,iCell), east(2,iCell), east(3,iCell)) + east(1,iCell) = -sin(lonCell_rad(iCell)) + east(2,iCell) = cos(lonCell_rad(iCell)) + east(3,iCell) = 0.0_r8 + call r3_normalize(east(1,iCell), east(2,iCell), east(3,iCell)) - north(1,iCell) = -cos(lonCell_rad(iCell))*sin(latCell_rad(iCell)) - north(2,iCell) = -sin(lonCell_rad(iCell))*sin(latCell_rad(iCell)) - north(3,iCell) = cos(latCell_rad(iCell)) - call r3_normalize(north(1,iCell), north(2,iCell), north(3,iCell)) + north(1,iCell) = -cos(lonCell_rad(iCell))*sin(latCell_rad(iCell)) + north(2,iCell) = -sin(lonCell_rad(iCell))*sin(latCell_rad(iCell)) + north(3,iCell) = cos(latCell_rad(iCell)) + call r3_normalize(north(1,iCell), north(2,iCell), north(3,iCell)) enddo @@ -7347,11 +4380,11 @@ subroutine uv_cell_to_edges(zonal_wind, meridional_wind, uedge, full_u) endif ! if(.not.on_outermost_edge(iEdge)) then enddo ! do iEdge = 1, nEdges -end subroutine uv_cell_to_edges - - +end subroutine project_uv_cell_to_edges + !------------------------------------------------------------------ + !================================================================== ! The following (private) routines were borrowed from the MPAS code !================================================================== @@ -7396,18 +4429,6 @@ function theta_to_tk (ens_size, theta, rho, qv, istatus) qv_nonzero = max(qv,0.0_r8) theta_to_tk = missing_r8 - -if ( debug > 0 .and. do_output()) then - if (any(istatus /= 0)) then - print *, 'theta_to_tk - nonzero istatus coming in' - do e = 1, ens_size - if (istatus(e) /= 0) then - write(string2, *) 'member ', e, ' incoming istatus = ', istatus(e) - call error_handler(E_ALLMSG, 'theta_to_tk',string2,source, revision, revdate) - endif !(istatus(e) /= 0) then - enddo ! ens_size - endif -endif where (istatus == 0) theta_m = (1.0_r8 + rvord * qv_nonzero)*theta @@ -7421,7 +4442,7 @@ function theta_to_tk (ens_size, theta, rho, qv, istatus) elsewhere - istatus = 89 + istatus = TK_COMPUTATION_ERROR endwhere @@ -7462,36 +4483,10 @@ subroutine compute_full_pressure(ens_size, theta, rho, qv, pressure, tk, istatus pressure = rho * rgas * tk * (1.0_r8 + rvord * qv_nonzero) end where -if ( debug > 1 ) then - if( any(istatus /= 0) ) then - do e = 1, ens_size - if (istatus(e) /= 0) then - write(string2,'("Failed in member,istatus,P[Pa],tk,theta,rho,qv:",2I4,4F12.2,F15.6)') & - e, istatus(e), pressure(e), tk(e), theta(e), rho(e), qv(e) - call error_handler(E_ALLMSG, 'compute_full_pressure',string2,source, revision, revdate) - endif !(istatus(e) /= 0) then - enddo ! ens_size - endif ! debug -else - if ((debug > 12) .and. do_output()) print *, 'compute_full_pressure: theta,r,q,p,tk,istatus =', & - theta, rho, qv, pressure, tk, istatus -endif - end subroutine compute_full_pressure -!-------------------------------------------------------------------- -!> pass the vertical localization coordinate to assim_tools_mod -function query_vert_localization_coord() - -integer :: query_vert_localization_coord - -query_vert_localization_coord = vert_localization_coord - -end function query_vert_localization_coord - !-------------------------------------------------------------------- !> read the time from the input file -!> stolen get_analysis_time_fname function read_model_time(filename) character(len=256), intent(in) :: filename @@ -7548,7 +4543,8 @@ subroutine set_wrf_date (tstring, year, month, day, hour, minute, second) end subroutine set_wrf_date -!=================================================================== -! End of model_mod -!=================================================================== +!----------------------------------------------------------------------! + + end module model_mod + diff --git a/models/mpas_atm/mpas_dart_obs_preprocess.f90 b/models/mpas_atm/mpas_dart_obs_preprocess.f90 index 68315475a0..7123db1d11 100644 --- a/models/mpas_atm/mpas_dart_obs_preprocess.f90 +++ b/models/mpas_atm/mpas_dart_obs_preprocess.f90 @@ -71,7 +71,7 @@ program mpas_dart_obs_preprocess RADIOSONDE_V_WIND_COMPONENT, SAT_U_WIND_COMPONENT, SAT_V_WIND_COMPONENT use model_mod, only : static_init_model, get_grid_dims, get_xland, & model_interpolate, find_closest_cell_center, & - cell_ok_to_interpolate, is_global_grid, & + cell_ok_to_interpolate, & get_bdy_mask, get_cell_center_coords use ensemble_manager_mod, only : ensemble_type, init_ensemble_manager, end_ensemble_manager use netcdf @@ -80,8 +80,6 @@ program mpas_dart_obs_preprocess ! version controlled file description for error handling, do not edit character(len=*), parameter :: source = 'models/mpas_atm/mpas_dart_obs_preprocess.f90' -character(len=*), parameter :: revision = '' -character(len=*), parameter :: revdate = '' ! ---------------------------------------------------------------------- ! Declare namelist parameters @@ -988,7 +986,7 @@ subroutine read_and_parse_input_seq(filename, landmask, obs_bdy_dist, siglevel,& character(len=129) :: qcmeta integer :: fid, var_id, okind, cellid, dsec, nobs, nth_obs integer :: bsec, bday, esec, eday, num_excluded_bytime -logical :: file_exist, last_obs, input_ncep_qc, global +logical :: file_exist, last_obs, input_ncep_qc real(r8), allocatable :: qc(:) real(r8) :: llv_loc(3) @@ -1010,8 +1008,6 @@ subroutine read_and_parse_input_seq(filename, landmask, obs_bdy_dist, siglevel,& call init_obs(prev_obs, get_num_copies(seq), get_num_qc(seq)) allocate(qc(get_num_qc(seq))) -global = is_global_grid() - input_ncep_qc = .false. qcmeta = get_qc_meta_data(seq, 1) if ( trim(adjustl(qcmeta)) == 'NCEP QC index' ) input_ncep_qc = .true. @@ -1547,7 +1543,7 @@ subroutine superob_aircraft_data(seq, ncell, atime, vdist, iqc_thres, ptop, obsk icell = find_closest_cell_center(airobs(k)%lat, airobs(k)%lon) if(icell < 1) then write(string,*) 'Cannot find any cell for this obs at ', airobs(k)%lat, airobs(k)%lon - call error_handler(E_MSG, routine, source, revision, revdate, text2=string) + call error_handler(E_MSG, routine, source, text2=string) endif if(airobs(k)%pressure > ps) then ik = 1 @@ -1919,7 +1915,7 @@ subroutine superob_sat_wind_data(seq, ncell, atime, vdist, iqc_thres, ptop) if(icell < 1) then write(string,*) 'Cannot find any cell for this obs at ',& satobs(k)%lat,satobs(k)%lon - call error_handler(E_MSG, routine, source, revision, revdate, text2=string) + call error_handler(E_MSG, routine, source, text2=string) endif if(satobs(k)%pressure > ps) then ik = 1 diff --git a/models/mpas_atm/update_bc.f90 b/models/mpas_atm/update_bc.f90 index 1fad8858e9..1a06bbd231 100644 --- a/models/mpas_atm/update_bc.f90 +++ b/models/mpas_atm/update_bc.f90 @@ -24,175 +24,286 @@ program update_bc use utilities_mod, only : initialize_utilities, finalize_utilities, & find_namelist_in_file, check_namelist_read, & logfileunit, open_file, close_file, & - get_next_filename, E_ERR, error_handler + get_next_filename, E_ERR, E_MSG, error_handler use time_manager_mod, only : time_type, print_time, print_date, operator(-), & get_time, get_date, operator(/=) -use direct_netcdf_mod,only : read_transpose, read_variables -use model_mod, only : static_init_model, statevector_to_analysis_file, & - get_model_size, get_init_template_filename, & - get_analysis_time, statevector_to_boundary_file, & - print_variable_ranges, & ! , get_num_vars - set_lbc_variables, force_u_into_state -use state_structure_mod, only : get_num_variables, get_domain_size +use model_mod, only : static_init_model, & + get_model_size, & + get_analysis_time, & + uv_increments_cell_to_edges, & + uv_field_cell_to_edges, & + get_analysis_weight, & + on_boundary_cell, & + on_boundary_edge, & + cell_next_to_boundary_edge, & + get_grid_dims + + +use state_structure_mod, only : get_num_variables, get_domain_size, & + get_variable_name + use netcdf_utilities_mod, only : nc_open_file_readonly, & nc_open_file_readwrite, & nc_get_dimension_size, & - nc_close_file - -use netcdf + nc_close_file, & + NF90_MAX_NAME, & + nc_get_variable, nc_put_variable implicit none -! version controlled file description for error handling, do not edit character(len=*), parameter :: source = 'models/mpas_atm/update_bc.f90' -character(len=*), parameter :: revision = '' -character(len=*), parameter :: revdate = '' !------------------------------------------------------------------ ! The namelist variables !------------------------------------------------------------------ - character(len=256) :: update_analysis_file_list = 'filter_in.txt' character(len=256) :: update_boundary_file_list = 'boundary_inout.txt' -integer :: debug = 0 logical :: lbc_update_from_reconstructed_winds = .true. logical :: lbc_update_winds_from_increments = .true. -namelist /update_bc_nml/ update_analysis_file_list, update_boundary_file_list, debug, & +namelist /update_bc_nml/ update_analysis_file_list, update_boundary_file_list, & lbc_update_from_reconstructed_winds, lbc_update_winds_from_increments !---------------------------------------------------------------------- character (len=256) :: next_infile, next_outfile character (len=256) :: bdy_template_filename -character (len=256) :: static_filename -character (len=512) :: string1 -integer :: iunit, io, x_size, nanlvars, nbdyvars -integer :: d1size, d2size -integer :: ncAnlID, ncBdyID, istatus +character (len=512) :: string1, string2 +integer :: iunit, io +integer :: ncAnlID, ncBdyID integer :: filenum -real(r8), allocatable :: statevector(:) type(time_type) :: model_time type(time_type) :: state_time integer :: nCellsA = -1 ! Total number of cells in ncAnlID integer :: nCellsB = -1 ! Total number of cells in ncBdyID integer :: nVertLevelsA = -1 ! Total number of levels in ncAnlID integer :: nVertLevelsB = -1 ! Total number of levels in ncBdyID -!---------------------------------------------------------------------- - -call initialize_utilities(progname='update_bc') -! Read the namelist to get the input filename. +integer :: nCells = -1 ! Total number of cells making up the grid +integer :: nVertices = -1 ! Unique points in grid that are corners of cells +integer :: nEdges = -1 ! Straight lines between vertices making up cells +integer :: nVertLevels = -1 ! Vertical levels; count of vert cell centers +integer :: vertexDegree = -1 ! Max number of cells/edges that touch any vertex +integer :: nSoilLevels = -1 ! Number of soil layers -call find_namelist_in_file("input.nml", "update_bc_nml", iunit) -read(iunit, nml = update_bc_nml, iostat = io) -call check_namelist_read(iunit, io, "update_bc_nml") +integer :: cellid, edgeid, i -!---------------------------------------------------------------------- -! Call model_mod:static_init_model() which reads the model namelists -! to set grid sizes, etc. -!---------------------------------------------------------------------- +logical :: lbc_file_has_reconstructed_winds ! HK @todo may always be false +character(len=NF90_MAX_NAME), dimension(9) :: lbc_variables +real(r8), allocatable, dimension(:,:) :: lbc_u, lbc_ucell, lbc_vcell, inc_lbc_ucell, inc_lbc_vcell +real(r8), allocatable, dimension(:,:) :: old_lbc_ucell, old_lbc_vcell, delta_u +real(r8), allocatable, dimension(:,:) :: a_var_data, b_var_data, var_data +real(r8) :: weight -! get the first member file to use as a template -bdy_template_filename = get_next_filename(update_boundary_file_list, 1) -! Note that force_u_into_state should be called before static_init_model, which is unusual. -call force_u_into_state() -call set_lbc_variables(bdy_template_filename) +character(len=NF90_MAX_NAME) :: bvarname, avarname -call static_init_model() -call get_init_template_filename(static_filename) +!---------------------------------------------------------------------- -x_size = get_model_size() -allocate(statevector(x_size)) +call initialize_utilities(progname=source) -! use get_num_variables() after setting the domains -! separate number for analysis file, boundary file -nanlvars = get_num_variables(1) -nbdyvars = get_num_variables(2) +! Read the namelist to get the input filename. +call find_namelist_in_file("input.nml", "update_bc_nml", iunit) +read(iunit, nml = update_bc_nml, iostat = io) +call check_namelist_read(iunit, io, "update_bc_nml") -write(*,*) -write(*,*) 'update_bc: Updating ',nbdyvars,' variables' +call static_init_model() +call get_grid_dims(nCells, nVertices, nEdges, nVertLevels, vertexDegree, nSoilLevels) + +lbc_file_has_reconstructed_winds = .false. +lbc_variables(1) = 'lbc_qc' +lbc_variables(2) = 'lbc_qr' +lbc_variables(3) = 'lbc_qv' +lbc_variables(4) = 'lbc_rho' +lbc_variables(5) = 'lbc_theta' +lbc_variables(6) = 'lbc_u' +lbc_variables(7) = 'lbc_w' +lbc_variables(8) = 'lbc_ur' ! may not be in the file +lbc_variables(9) = 'lbc_vr' ! may not be in the file !---------------------------------------------------------------------- ! Reads lists of input mpas (prior) and filter (analysis) files +!HK @todo loop around files, why not run this code in parallel? !---------------------------------------------------------------------- filenum = 1 -fileloop: do ! until out of files - - ! get a file name from the list, one at a time. - next_infile = get_next_filename(update_analysis_file_list, filenum) - next_outfile = get_next_filename(update_boundary_file_list, filenum) - if (next_infile == '' .or. next_outfile == '') exit fileloop - - !---------------------------------------------------------------------- - ! Reads input lbc (prior) and filter (analysis) files - !---------------------------------------------------------------------- - - ncAnlID = nc_open_file_readonly(next_infile, 'update_bc - open readonly') - ncBdyID = nc_open_file_readwrite(next_outfile, 'update_bc - open readwrite') - - !---------------------------------------------------------------------- - ! Read the model time - !---------------------------------------------------------------------- - model_time = get_analysis_time(ncBdyID, next_outfile) - state_time = get_analysis_time(ncAnlID, next_infile) - call print_time(state_time,'DART current time') - call print_time(model_time,'mpas current time') - - if ( model_time /= state_time ) then - call print_time(state_time,'DART current time',logfileunit) - call print_time(model_time,'mpas current time',logfileunit) - write(string1,*) trim(next_infile),' current time must equal model time' - call error_handler(E_ERR,'update_bc',string1,source,revision,revdate) - endif - - !---------------------------------------------------------------------- - ! Check dimension size in both files - !---------------------------------------------------------------------- - nCellsA = nc_get_dimension_size(ncAnlID, 'nCells', 'update_bc') ! Ha - nCellsB = nc_get_dimension_size(ncBdyID, 'nCells', 'update_bc') ! Ha - nVertLevelsA = nc_get_dimension_size(ncAnlID, 'nVertLevels', 'update_bc') ! Ha - nVertLevelsB = nc_get_dimension_size(ncBdyID, 'nVertLevels', 'update_bc') ! Ha - print*,'nCells, nVertLevels:', nCellsA, nVertLevelsA,' in ',trim(next_infile) ! Ha - - if((nCellsA /= nCellsB) .or. (nVertLevelsA /= nVertLevelsB)) then ! Ha - print*,'nCells, nVertLevels:', nCellsB, nVertLevelsB,' in ',trim(next_outfile) - write(string1,*) 'Domain size mismatches' - call error_handler(E_ERR,'update_bc',string1,source,revision,revdate) - endif - - !---------------------------------------------------------------------- - ! Read analysis state vector (assuming to be available at the model time) - !---------------------------------------------------------------------- - d1size = get_domain_size(1) - d2size = get_domain_size(2) - - call read_variables(ncAnlID, statevector(1:d1size), 1, nanlvars, domain=1) - call read_variables(ncBdyID, statevector(d1size+1:d1size+d2size), 1, nbdyvars, domain=2) - - !---------------------------------------------------------------------- - ! update the current model state vector - !---------------------------------------------------------------------- - write(*,*) 'Updating boundary variables in ',trim(next_outfile) - call statevector_to_boundary_file(statevector, ncBdyID, ncAnlID, & - lbc_update_from_reconstructed_winds, lbc_update_winds_from_increments, debug) - - !---------------------------------------------------------------------- - ! Log what we think we're doing, and exit. - !---------------------------------------------------------------------- - +fileloop: do ! until out of files + + ! get a file name from the list, one at a time. + next_infile = get_next_filename(update_analysis_file_list, filenum) + next_outfile = get_next_filename(update_boundary_file_list, filenum) + if (next_infile == '' .or. next_outfile == '') exit fileloop + + !---------------------------------------------------------------------- + ! Reads input lbc (prior) and filter (analysis) files + !---------------------------------------------------------------------- + ncAnlID = nc_open_file_readonly(next_infile, 'update_bc - open readonly') ! analysis file from DART (ouput from filter for anl_dom) + ncBdyID = nc_open_file_readwrite(next_outfile, 'update_bc - open readwrite') ! prior boundary, original mpas file + + !---------------------------------------------------------------------- + ! Read the model time + !---------------------------------------------------------------------- + state_time = get_analysis_time(ncAnlID, next_infile) + model_time = get_analysis_time(ncBdyID, next_outfile) + call print_time(state_time,'DART current time') + call print_time(model_time,'mpas current time') + + if ( model_time /= state_time ) then + call print_time(state_time,'DART current time',logfileunit) + call print_time(model_time,'mpas current time',logfileunit) + write(string1,*) trim(next_infile),' current time must equal model time' + call error_handler(E_ERR,source,string1,source) + endif + + !---------------------------------------------------------------------- + ! Check dimension size in both files + !---------------------------------------------------------------------- + nCellsA = nc_get_dimension_size(ncAnlID, 'nCells', source) + nCellsB = nc_get_dimension_size(ncBdyID, 'nCells', source) + nVertLevelsA = nc_get_dimension_size(ncAnlID, 'nVertLevels', source) + nVertLevelsB = nc_get_dimension_size(ncBdyID, 'nVertLevels', source) + + if((nCellsA /= nCellsB) .or. (nVertLevelsA /= nVertLevelsB)) then + ! HK @todo also check against static_init_model values + write(string1,*) 'Domain size mismatches' + call error_handler(E_ERR,'update_bc',string1,source) + endif + + if (lbc_update_from_reconstructed_winds) then ! save a copy of the reconstruced cell winds + + if (.not. lbc_file_has_reconstructed_winds) then + write(string1, *) 'Cannot update edge winds from increments because the boundary file does not contain the reconstructed winds (lbc_ur, lbc_vr)' + write(string2, *) 'lbc_update_winds_from_increments should be .false.' + call error_handler(E_MSG,'statevector_to_boundary_file',string1,& + source, text2=string2) + endif + + allocate(old_lbc_ucell(nVertLevels, nCells)) + allocate(old_lbc_vcell(nVertLevels, nCells)) + allocate( delta_u(nVertLevels, nEdges)) + + call nc_get_variable(ncBdyID, 'uReconstructMeridional', old_lbc_ucell) + call nc_get_variable(ncBdyID, 'uReconstructZonal', old_lbc_vcell) + + endif + + + ! Update variables except 'u' from the analysis + allocate(a_var_data(nVertLevels, nCells)) + allocate(b_var_data(nVertLevels, nCells)) + allocate(var_data(nVertLevels, nCells)) + + VARLOOP: do i = 1, size(lbc_variables) + bvarname = lbc_variables(i) + avarname = trim(bvarname(5:)) !corresponding field in analysis domain + + ! skip edge normal winds + if (bvarname == 'u') cycle VARLOOP + + ! reconstructed cell-center winds have different names in the lbc file. + if (bvarname == 'lbc_ur') avarname = 'uReconstructZonal' + if (bvarname == 'lbc_vr') avarname = 'uReconstructMeridional' + + call nc_get_variable(ncAnlID, avarname, a_var_data) + call nc_get_variable(ncBdyID, bvarname, b_var_data) + + ! for each cell in the grid, find the analysis in the + ! boundary region and blend them with prior lbc values. + CELLS: do cellid = 1, nCells + + if (.not. on_boundary_cell(cellid)) cycle CELLS + + weight = get_analysis_weight(cellid) ! 1.0 is interior, 0.0 is exterior boundary + var_data(:, cellid) = (1.0_r8 - weight) * b_var_data(:, cellid) + weight * a_var_data(:, cellid) + + enddo CELLS + + call nc_put_variable(ncBdyID, bvarname, var_data) + + enddo VARLOOP + + deallocate(a_var_data, b_var_data, var_data) + + ! deal with u edge winds + if (.not. lbc_update_from_reconstructed_winds) then ! update u edge winds directly + + allocate(a_var_data(nVertLevels, nEdges)) + allocate(b_var_data(nVertLevels, nEdges)) + allocate(var_data(nVertLevels, nEdges)) + + call nc_get_variable(ncAnlID, 'u', a_var_data) ! analysis edge winds + call nc_get_variable(ncBdyID, 'lbc_u', b_var_data) ! prior edge winds in the lbc file + + ! for each edge in the grid, find the ones which are in the + ! boundary region and blend their values. + EDGES: do edgeid = 1, nEdges + + if (.not. on_boundary_edge(edgeid)) cycle EDGES + + weight = get_analysis_weight(edgeid, .false.) ! 1.0 is interior, 0.0 is exterior boundary + var_data(:, edgeid) = (1.0_r8 - weight) * b_var_data(:,edgeid) + weight * a_var_data(:,edgeid) + + enddo EDGES + + call nc_put_variable(ncBdyID, 'lbc_u', var_data) + + else ! update u edge winds from reconstructed winds + + allocate( lbc_u(nVertLevels, nEdges)) + allocate( lbc_ucell(nVertLevels, nCells)) + allocate( lbc_vcell(nVertLevels, nCells)) + + call nc_get_variable(ncBdyID, 'lbc_ur', lbc_vcell) ! already blended in VARLOOP + call nc_get_variable(ncBdyID, 'lbc_vr', lbc_ucell) ! already blended in VALOOP + + if (lbc_update_winds_from_increments) then + + call nc_get_variable(ncBdyID, 'lbc_u', lbc_u) ! not blended + allocate(inc_lbc_ucell(nVertLevels, nCells)) + allocate(inc_lbc_vcell(nVertLevels, nCells)) + allocate(delta_u(nVertLevels, nEdges)) + + inc_lbc_ucell = lbc_ucell - old_lbc_ucell + inc_lbc_vcell = lbc_vcell - old_lbc_vcell + + delta_u = uv_increments_cell_to_edges(inc_lbc_ucell, inc_lbc_vcell) + + IEDGE: do edgeid = 1, nEdges + + ! Soyoung: Add the blended u increments back to lbc_u in the boundary zone. + ! We should not change the analysis u in the interior domain, but + ! We should also check bdyMaskCell for the two adjacent cells as + ! bdyMaskEdge is assigned with the lower mask value between the two + ! cells. + ! Ex) An edge between cell1 (w/ bdyMaskCell = 0) and cell2 (w/ bdyMaskCell = 1) + ! has bdyMaskEdge = 0. In this case, even if bdyMaskEdge of the edge is zero, + ! cell2 has been updated in the CELLS loop above, so the edge has to be updated. + + if (.not. on_boundary_edge(edgeid) .and. .not. cell_next_to_boundary_edge(edgeid)) cycle IEDGE + + lbc_u(:,edgeid) = lbc_u(:,edgeid) + delta_u(:,edgeid) + + enddo IEDGE + + call nc_put_variable(ncBdyID, 'lbc_u', lbc_u) + deallocate(old_lbc_ucell, old_lbc_vcell, inc_lbc_ucell,inc_lbc_vcell, delta_u) + + else ! just replace, no increments + + call uv_field_cell_to_edges(lbc_ucell, lbc_vcell, lbc_u) + call nc_put_variable(ncBdyID, 'lbc_u', var_data) + + endif + + endif + call print_date( model_time,'update_bc:model date') call print_time( model_time,'update_bc:model time') call print_date( model_time,'update_bc:model date',logfileunit) call print_time( model_time,'update_bc:model time',logfileunit) - ! Because the files were open with the nc_open_file...() routines, - ! it is not necessary to supply the filenames for error msg purposes. - - call nc_close_file(ncAnlID,'update_bc') - call nc_close_file(ncBdyID,'update_bc') + call nc_close_file(ncAnlID, source) + call nc_close_file(ncBdyID, source) filenum = filenum + 1 diff --git a/models/mpas_atm/update_mpas_states.f90 b/models/mpas_atm/update_mpas_states.f90 index 0c0141f88c..83f79a5502 100644 --- a/models/mpas_atm/update_mpas_states.f90 +++ b/models/mpas_atm/update_mpas_states.f90 @@ -14,7 +14,10 @@ program update_mpas_states ! The update_mpas_states_nml namelist defines the input and output file ! name lists for all ensemble members. ! The input list should be matched with output_state_file_list in &filter_nml. -! +! +! variables that are not wind, copied from one file to the another +! variables that are wind, reconstructed +! ! author: Soyoung Ha 23 Aug 16 ! Updated in 4 May 2017 for the Manhattan release !---------------------------------------------------------------------- @@ -22,75 +25,62 @@ program update_mpas_states use types_mod, only : r8 use utilities_mod, only : initialize_utilities, finalize_utilities, & find_namelist_in_file, check_namelist_read, & - logfileunit, get_next_filename, E_ERR, error_handler + logfileunit, get_next_filename, E_ERR, E_MSG, error_handler use time_manager_mod, only : time_type, print_time, print_date, operator(-), & get_time, get_date, operator(/=) -use direct_netcdf_mod,only : read_transpose, read_variables -use model_mod, only : static_init_model, statevector_to_analysis_file, & +use model_mod, only : static_init_model, & get_model_size, & - get_num_vars, get_analysis_time, & - print_variable_ranges + get_analysis_time, update_u_from_reconstruct, & + use_increments_for_u_update, uv_increments_cell_to_edges, & + uv_field_cell_to_edges, dom_id => anl_domid + +use state_structure_mod, only : get_num_variables, get_variable_name, & + get_variable_size, get_varid_from_varname, & + get_dim_lengths use netcdf_utilities_mod, only : nc_open_file_readonly, & nc_open_file_readwrite, & + nc_get_variable, nc_put_variable, & nc_close_file - -use netcdf - + implicit none -! version controlled file description for error handling, do not edit character(len=*), parameter :: source = 'models/mpas_atm/update_mpas_states.f90' -character(len=*), parameter :: revision = '' -character(len=*), parameter :: revdate = '' !------------------------------------------------------------------ ! The namelist variables !------------------------------------------------------------------ - character(len=256) :: update_input_file_list = 'filter_out.txt' character(len=256) :: update_output_file_list = 'filter_in.txt' -logical :: print_data_ranges = .true. -integer :: dom_id = 1 -namelist /update_mpas_states_nml/ update_input_file_list, update_output_file_list, dom_id, & - print_data_ranges +namelist /update_mpas_states_nml/ update_input_file_list, update_output_file_list + !---------------------------------------------------------------------- character (len=256) :: next_infile, next_outfile character (len=512) :: string1 -integer :: iunit, io, x_size, nvars +integer :: iunit, io integer :: ncAnlID, ncBckID, istatus -integer :: filenum -real(r8), allocatable :: statevector(:) +integer :: filenum, i +real(r8), allocatable :: variable(:) type(time_type) :: model_time type(time_type) :: state_time +integer :: dims(2) !(nVertLevels, nEdges | nCells) +real(r8), allocatable :: u(:,:), ucell(:,:), vcell(:,:) +real(r8), allocatable :: ucell_dart(:,:), vcell_dart(:,:), increments(:,:) !---------------------------------------------------------------------- call initialize_utilities(progname=source) ! Read the namelist to get the input filename. - call find_namelist_in_file("input.nml", "update_mpas_states_nml", iunit) read(iunit, nml = update_mpas_states_nml, iostat = io) call check_namelist_read(iunit, io, "update_mpas_states_nml") -!---------------------------------------------------------------------- -! Call model_mod:static_init_model() which reads the model namelists -! to set grid sizes, etc. -!---------------------------------------------------------------------- - call static_init_model() -nvars = get_num_vars() -x_size = get_model_size() -allocate(statevector(x_size)) - -write(*,*) -write(*,*) 'update_mpas_states: Updating ',nvars,' variables' -write(*,*) - !---------------------------------------------------------------------- ! Reads lists of input mpas (prior) and filter (analysis) files +! HK @todo this is bad to have a serial loop around files. !---------------------------------------------------------------------- filenum = 1 fileloop: do ! until out of files @@ -100,12 +90,9 @@ program update_mpas_states next_outfile = get_next_filename(update_output_file_list, filenum) if (next_infile == '' .or. next_outfile == '') exit fileloop - ncAnlID = nc_open_file_readonly(next_infile, 'update_mpas_states - open readonly') - ncBckID = nc_open_file_readwrite(next_outfile, 'update_mpas_states - open readwrite') + ncAnlID = nc_open_file_readonly(next_infile, 'update_mpas_states - open readonly') ! analysis from DART + ncBckID = nc_open_file_readwrite(next_outfile, 'update_mpas_states - open readwrite') ! background, original mpas file - !---------------------------------------------------------------------- - ! Read the model time - !---------------------------------------------------------------------- model_time = get_analysis_time(ncBckID, next_outfile) state_time = get_analysis_time(ncAnlID, next_infile) call print_time(state_time,'DART current time') @@ -115,48 +102,105 @@ program update_mpas_states call print_time(state_time,'DART current time',logfileunit) call print_time(model_time,'mpas current time',logfileunit) write(string1,*) trim(next_infile),' current time must equal model time' - call error_handler(E_ERR,'update_mpas_states',string1,source,revision,revdate) + call error_handler(E_ERR,'update_mpas_states',string1,source) endif - !---------------------------------------------------------------------- - ! Read analysis state vector (assuming to be available at the model time) - !---------------------------------------------------------------------- - call read_variables(ncAnlID, statevector, 1, nvars, dom_id) - - !---------------------------------------------------------------------- - ! if requested, print out the data ranges variable by variable - ! (note if we are clamping data values, that happens in the conversion - ! routine below and these values are before the clamping happens.) - !---------------------------------------------------------------------- - if (print_data_ranges) then - write(*,*) - write(*,*) ' Input: ', trim(next_infile) - write(*,*) 'Output: ', trim(next_outfile) - call print_variable_ranges(statevector, 'Analysis states') - write(*,*) - endif + ! copy variables that are not wind from analysis to background + varloop: do i = 1, get_num_variables(dom_id) + + allocate(variable(get_variable_size(dom_id, i))) + + if (get_variable_name(dom_id,i) == 'uReconstructZonal' .or. & + get_variable_name(dom_id,i) == 'uReconstructMeridional'.or. & + get_variable_name(dom_id,i) == 'u') cycle varloop + + call nc_get_variable(ncAnlID, get_variable_name(dom_id,i), variable) + call nc_put_variable(ncBckID, get_variable_name(dom_id,i), variable) + + deallocate(variable) + + enddo varloop - !---------------------------------------------------------------------- - ! update the current model state vector - !---------------------------------------------------------------------- - write(*,*) 'Overwriting states in ',trim(next_outfile) - call statevector_to_analysis_file(statevector, ncBckID, next_outfile) + ! deal with wind + if (update_u_from_reconstruct) then - !---------------------------------------------------------------------- - ! Log what we think we're doing, and exit. - !---------------------------------------------------------------------- + if (use_increments_for_u_update) then + ! Read in the previous reconstructed winds from the original mpas netcdf file + ! and compute what increments (changes in values) were added by the assimilation. + ! Read in the original edge normal wind 'u' field from that same mpas netcdf + ! file and add the interpolated increments to compute the updated 'u' values. + + ! read in u, uReconstrtuctZonal, uReconstructMeridional from background + ! read in uReconstrtuctZonal, uReconstructMeridional from analysis + + dims = get_dim_lengths(dom_id, get_varid_from_varname(dom_id, 'u')) + allocate(u(dims(1), dims(2)), increments(dims(1), dims(2))) + call nc_get_variable(ncBckID, 'u', u) + + dims = get_dim_lengths(dom_id, get_varid_from_varname(dom_id, 'uReconstructZonal')) + allocate(ucell(dims(1), dims(2)), ucell_dart(dims(1), dims(2))) + call nc_get_variable(ncBckID, 'uReconstructZonal', ucell) + call nc_get_variable(ncAnlID, 'uReconstructZonal', ucell_dart) + ucell = ucell_dart - ucell ! u increments + + dims = get_dim_lengths(dom_id, get_varid_from_varname(dom_id, 'uReconstructMeridional')) + allocate(vcell(dims(1), dims(2)), vcell_dart(dims(1), dims(2))) + call nc_get_variable(ncBckID, 'uReconstructMeridional', vcell) + call nc_get_variable(ncAnlID, 'uReconstructMeridional', vcell_dart) + vcell = vcell_dart - vcell ! v increments + + u = u + uv_increments_cell_to_edges(ucell, vcell) + + call nc_put_variable(ncBckID, 'u', u) + call nc_put_variable(ncBckID, 'uReconstructZonal', ucell_dart) + call nc_put_variable(ncBckID, 'uReconstructMeridional', vcell_dart) + + deallocate(u, increments, ucell, vcell, ucell_dart, vcell_dart) + + else + ! The state vector has updated zonal and meridional wind components. + ! put them directly into the arrays. These are the full values, not + ! just increments. + dims = get_dim_lengths(dom_id, get_varid_from_varname(dom_id, 'u')) + allocate(u(dims(1), dims(2))) + dims = get_dim_lengths(dom_id, get_varid_from_varname(dom_id, 'uReconstructZonal')) + allocate(ucell_dart(dims(1), dims(2))) + call nc_get_variable(ncAnlID, 'uReconstructZonal', ucell_dart) + dims = get_dim_lengths(dom_id, get_varid_from_varname(dom_id, 'uReconstructMeridional')) + allocate(vcell_dart(dims(1), dims(2))) + call nc_get_variable(ncAnlID, 'uReconstructMeridional', vcell_dart) + + call uv_field_cell_to_edges(ucell_dart, vcell_dart, u) + + call nc_put_variable(ncBckID, 'u', variable) + call nc_put_variable(ncBckID, 'uReconstructZonal', variable) + call nc_put_variable(ncBckID, 'uReconstructMeridional', variable) + + deallocate(u, ucell_dart, vcell_dart) + + endif + + else + ! copy u from analysis to background + allocate(variable(get_variable_size(1, get_varid_from_varname(dom_id, 'u')))) + + call nc_get_variable(ncAnlID, 'u', variable) + call nc_put_variable(ncBckID, 'u', variable) + + deallocate(variable) + + endif + + call error_handler(E_MSG, 'Overwritten states in ',trim(next_outfile), source) call print_date( model_time,'update_mpas_states:model date') call print_time( model_time,'update_mpas_states:model time') call print_date( model_time,'update_mpas_states:model date',logfileunit) call print_time( model_time,'update_mpas_states:model time',logfileunit) - ! Because the files were open with the nc_open_file...() routines, - ! it is not necessary to supply the filenames for error msg purposes. - call nc_close_file(ncAnlID,'update_mpas_states') call nc_close_file(ncBckID,'update_mpas_states') - + filenum = filenum + 1 end do fileloop diff --git a/models/mpas_atm/work/input.nml b/models/mpas_atm/work/input.nml index 10a2171162..6c4fb8f0e0 100644 --- a/models/mpas_atm/work/input.nml +++ b/models/mpas_atm/work/input.nml @@ -169,7 +169,7 @@ &model_nml - init_template_filename = 'init.nc' + init_template_filename = 'restart.2014-09-15_00.00.00.nc' assimilation_period_days = 0 assimilation_period_seconds = 21600 model_perturbation_amplitude = 0.0001 @@ -178,7 +178,6 @@ highest_obs_pressure_mb = 1.0 sfc_elev_max_diff = 100. log_p_vert_interp = .false. - debug = 0 use_u_for_wind = .false. use_rbf_option = 2 update_u_from_reconstruct = .true. @@ -211,6 +210,7 @@ 'w', 'QTY_VERTICAL_VELOCITY', 'qv', 'QTY_VAPOR_MIXING_RATIO', 'surface_pressure', 'QTY_SURFACE_PRESSURE', + 't2m' 'QTY_2M_TEMPERATURE', mpas_state_bounds = 'qv','0.0','NULL','CLAMP', / @@ -366,13 +366,13 @@ quantity_of_interest = 'QTY_VAPOR_MIXING_RATIO' &model_mod_check_nml - input_state_files = 'init.nc' + input_state_files = 'restart.2014-09-15_00.00.00.nc' output_state_files = 'check_me.nc' - test1thru = 0 - run_tests = 1,2,3,4,5,7 + test1thru = -1 + run_tests = 1,2,3,4 x_ind = 300 loc_of_interest = 240.0, 0.0, 10000.0 - quantity_of_interest = 'QTY_U_WIND_COMPONENT' + quantity_of_interest = 'QTY_V_WIND_COMPONENT' interp_test_lonrange = 0.0, 359.0 interp_test_dlon = 1.0 interp_test_latrange = -89.0, 89.0