[Dart-dev] [5871] DART/branches/development/models/mpas_atm/model_mod.f90: Restructure code, remove unused routines, remove some debugging
nancy at ucar.edu
nancy at ucar.edu
Tue Sep 18 15:32:24 MDT 2012
Revision: 5871
Author: nancy
Date: 2012-09-18 15:32:23 -0600 (Tue, 18 Sep 2012)
Log Message:
-----------
Restructure code, remove unused routines, remove some debugging
and development comments. Should give same results as previous
corrected versions (which were not committed). DOES NOT have the
latest speedup code - that will be committed separately.
Modified Paths:
--------------
DART/branches/development/models/mpas_atm/model_mod.f90
-------------- next part --------------
Modified: DART/branches/development/models/mpas_atm/model_mod.f90
===================================================================
--- DART/branches/development/models/mpas_atm/model_mod.f90 2012-09-18 21:05:24 UTC (rev 5870)
+++ DART/branches/development/models/mpas_atm/model_mod.f90 2012-09-18 21:32:23 UTC (rev 5871)
@@ -10,9 +10,17 @@
! $Revision$
! $Date$
-! This is the interface between the model model and DART.
+! MPAS Atmosphere model interfact 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.
-! Modules that are absolutely required for use are listed
+
+! Routines in other modules that are used here.
+
use types_mod, only : r4, r8, digits12, SECPERDAY, MISSING_R8, &
rad2deg, deg2rad, PI, MISSING_I
use time_manager_mod, only : time_type, set_time, set_date, get_date, get_time,&
@@ -23,15 +31,19 @@
use location_mod, only : location_type, get_dist, query_location, &
get_close_maxdist_init, get_close_type, &
- set_location, get_location, horiz_dist_only, &
+ set_location, get_location, horiz_dist_only, &
+ write_location, &
vert_is_undef, VERTISUNDEF, &
vert_is_surface, VERTISSURFACE, &
vert_is_level, VERTISLEVEL, &
vert_is_pressure, VERTISPRESSURE, &
vert_is_height, VERTISHEIGHT, &
+ vert_is_scale_height, VERTISSCALEHEIGHT, &
get_close_obs_init, get_close_obs_destroy, &
- loc_get_close_obs => get_close_obs
+ loc_get_close_obs
+ ! loc_get_close_obs => get_close_obs
+
use utilities_mod, only : register_module, error_handler, &
E_ERR, E_WARN, E_MSG, logfileunit, get_unit, &
nc_check, do_output, to_upper, nmlfileunit, &
@@ -57,8 +69,18 @@
use random_seq_mod, only: random_seq_type, init_random_seq, random_gaussian
+! 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
@@ -103,6 +125,8 @@
revision = '$Revision$', &
revdate = '$Date$'
+! module global storage; maintains values between calls, accessible by
+! any subroutine
character(len=256) :: string1, string2, string3
logical, save :: module_initialized = .false.
@@ -114,12 +138,15 @@
real(r8), parameter :: p0 = 100000.0_r8
real(r8), parameter :: rcv = rgas/(cp-rgas)
+! earth radius; needed to convert lat/lon to x,y,z cartesian coords.
! FIXME: one of the example ocean files had a global attr with 6371220.0
! instead of 1229. ??
real(r8), parameter :: radius = 6371229.0 ! meters
-! Storage for a random sequence for perturbing a single initial state
+! roundoff error
+real(r8), parameter :: roundoff = 1.0e-12_r8
+! 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
@@ -128,27 +155,31 @@
type(location_type), allocatable :: cell_locs(:)
integer, allocatable :: dummy(:), close_ind(:)
-! things which can/should be in the model_nml
-
+! variables which are in the module namelist
+integer :: vert_localization_coord = VERTISHEIGHT
integer :: assimilation_period_days = 0
integer :: assimilation_period_seconds = 60
real(r8) :: model_perturbation_amplitude = 0.0001 ! tiny amounts
-logical :: output_state_vector = .true.
+real(r8) :: highest_obs_pressure_mb = 100.0_r8 ! do not assimilate obs higher than this level.
+real(r8) :: cell_size_radians = 0.10_r8 ! size of largest cell in r
+logical :: output_state_vector = .false. ! output prognostic variables (if .false.)
+logical :: logp = .false. ! if true, interpolate pres in log space
integer :: debug = 0 ! turn up for more and more debug messages
character(len=32) :: calendar = 'Gregorian'
character(len=256) :: model_analysis_filename = 'mpas_analysis.nc'
character(len=256) :: grid_definition_filename = 'mpas_analysis.nc'
-! if .false. use the original triangle search code, no rbf available.
-! if .true. use search for closest vertex and has the rbf options
-logical :: use_new_code = .true.
+! FIXME: these are here to allow regression testing and should be
+! removed when we settle on the final version. for sure the 'old code'
+! for locating a point in a group of triangles should be removed
+! because there are known errors at high latitudes near the poles.
! 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
logical :: use_u_for_wind = .false.
! if using rbf, options 1,2,3 control how many points go into the rbf.
-! larger numbers are more points from further away
+! larger numbers use more points from further away
integer :: use_rbf_option = 2
! if .false. edge normal winds (u) should be in state vector and are written out directly
@@ -159,20 +190,18 @@
model_analysis_filename, &
grid_definition_filename, &
output_state_vector, &
+ vert_localization_coord, &
assimilation_period_days, &
assimilation_period_seconds, &
model_perturbation_amplitude, &
calendar, &
debug, &
- use_new_code, &
use_u_for_wind, &
use_rbf_option, &
- update_u_from_reconstruct
+ update_u_from_reconstruct, &
+ highest_obs_pressure_mb
-!------------------------------------------------------------------
-! DART state vector are specified in the input.nml:mpas_vars_nml namelist.
-!------------------------------------------------------------------
-
+! 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
character(len=NF90_MAX_NAME) :: mpas_state_variables(max_state_variables * num_state_table_columns ) = ' '
@@ -261,22 +290,23 @@
integer, allocatable :: boundaryVertex(:,:)
integer, allocatable :: maxLevelCell(:)
-real(r8), allocatable :: ens_mean(:) ! may be needed for forward ops
+real(r8), allocatable :: ens_mean(:) ! needed to convert vertical distances consistently
integer :: model_size ! the state vector length
type(time_type) :: model_timestep ! smallest time to adv model
+! useful flags in making decisions when searching for points, etc
logical :: global_grid = .true. ! true = the grid covers the sphere with no holes
logical :: all_levels_exist_everywhere = .true. ! true = cells defined at all levels
logical :: has_real_u = .false. ! true = has original u on edges
logical :: has_uvreconstruct = .false. ! true = has reconstructed at centers
-!------------------------------------------------------------------
-! The model analysis manager namelist variables
-!------------------------------------------------------------------
-
+! 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
@@ -611,27 +641,36 @@
endif
! do some sanity checking here?
-if (use_new_code) then
- 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
- if (update_u_from_reconstruct .and. has_real_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 real 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)
- endif
+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
+if (update_u_from_reconstruct .and. has_real_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 real 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)
+endif
-allocate( ens_mean(model_size) )
+! 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:
+if ((get_progvar_index_from_kind(KIND_POTENTIAL_TEMPERATURE) < 0) .or. &
+ (get_progvar_index_from_kind(KIND_DENSITY) < 0) .or. &
+ (get_progvar_index_from_kind(KIND_VAPOR_MIXING_RATIO) < 0)) then
+ write(string1, *) 'State vector is missing one or more of the following fields:'
+ write(string2, *) 'Potential Temperature, Density, Vapor Mixing Ratio.'
+ 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)
+endif
-! Initialize the interpolation data structures
-call init_interp()
+allocate( ens_mean(model_size) )
+
end subroutine static_init_model
@@ -653,7 +692,9 @@
integer :: nxp, nzp, iloc, vloc, nf, n
integer :: myindx
+integer :: ztypeout, istatus
real(r8) :: height
+type(location_type) :: new_location
if ( .not. module_initialized ) call static_init_model
@@ -709,8 +750,6 @@
endif
endif
-
-
if (progvar(nf)%numedges /= MISSING_I) then
if (nzp <= 1) then
location = set_location(lonEdge(iloc),latEdge(iloc), height, VERTISSURFACE)
@@ -725,8 +764,21 @@
endif
endif
-if (debug > 9) then
+! Let us return the vert location with the requested vertical localization coordinate
+! hoping that the code can run faster when same points are close to many obs
+! since vert_convert does not need to be called repeatedly in get_close_obs any more.
+! FIXME: we should test this to see if it's a win (in computation time). for obs
+! which are not dense relative to the grid, this might be slower than doing the
+! conversions on demand in the localization code (in get_close_obs()).
+if ( .not. horiz_dist_only .and. vert_localization_coord /= VERTISHEIGHT ) then
+ new_location = location
+ call vert_convert(ens_mean, new_location, progvar(nf)%dart_kind, vert_localization_coord, istatus)
+ if(istatus == 0) location = new_location
+endif
+
+if (debug > 12) then
+
write(*,'("INDEX_IN / myindx / IVAR / NX, NZ: ",2(i10,2x),3(i5,2x))') index_in, myindx, nf, nxp, nzp
write(*,'(" ILOC, KLOC: ",2(i5,2x))') iloc, vloc
write(*,'(" LON/LAT/HGT: ",3(f12.3,2x))') lonCell(iloc), latCell(iloc), height
@@ -748,122 +800,6 @@
! interpolated value at that location, and an error code. 0 is success,
! anything positive is an error. (negative reserved for system use)
!
-! for specific error codes, see the local_interpolate() comments
-
-! passed variables
-
-real(r8), intent(in) :: x(:)
-type(location_type), intent(in) :: location
-integer, intent(in) :: obs_type
-real(r8), intent(out) :: interp_val
-integer, intent(out) :: istatus
-
-! local storage
-
-integer :: obs_kinds(3), ivar
-real(r8) :: values(3)
-
-! call the normal interpolate code. if it fails because
-! the kind doesn't exist directly in the state vector, try a few
-! kinds we know how to convert.
-
-interp_val = MISSING_R8
-istatus = 888888 ! must be positive (and integer)
-
-obs_kinds(1) = obs_type
-
-! debug only
-if (.false.) then
- ivar = get_progvar_index_from_kind(obs_kinds(1))
- if ( ivar > 0 .and. debug > 7 ) call dump_progvar(ivar, x)
- ivar = get_progvar_index_from_kind(KIND_POTENTIAL_TEMPERATURE)
- if ( ivar > 0 .and. debug > 7 ) call dump_progvar(ivar, x)
- ivar = get_progvar_index_from_kind(KIND_DENSITY)
- if ( ivar > 0 .and. debug > 7 ) call dump_progvar(ivar, x)
- ivar = get_progvar_index_from_kind(KIND_VAPOR_MIXING_RATIO)
- if ( ivar > 0 .and. debug > 7 ) call dump_progvar(ivar, x)
- ivar = get_progvar_index_from_kind(KIND_EDGE_NORMAL_SPEED)
- if ( ivar > 0 .and. debug > 7 ) call dump_progvar(ivar, x)
-endif
-
-
-call local_interpolate(x, location, 1, obs_kinds, values, istatus)
-!print *, '0 local interpolate returns istatus = ', istatus
-if (istatus /= 88) then
- ! 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. FIXME.
- if (istatus /= 0 .and. values(1) /= MISSING_R8) then
- write(string1,*) 'interp routine returned a bad status but not a MISSING_R8 value'
- write(string2,*) 'value = ', values(1), ' istatus = ', istatus
- call error_handler(E_ERR,'model_interpolate',string1,source,revision,revdate, &
- text2=string2)
- endif
- if (istatus == 0 .and. values(1) == MISSING_R8) then
- write(string1,*) 'interp routine returned a good status but set value to MISSING_R8'
- call error_handler(E_ERR,'model_interpolate',string1,source,revision,revdate)
- endif
- interp_val = values(1)
-if (debug > 5) print *, 'returning value ', interp_val
- return
-endif
-
-if (obs_type == KIND_TEMPERATURE) then
- obs_kinds(1) = KIND_POTENTIAL_TEMPERATURE
- obs_kinds(2) = KIND_DENSITY
- obs_kinds(3) = KIND_VAPOR_MIXING_RATIO
-
- call local_interpolate(x, location, 3, obs_kinds, values, istatus)
-!print *, '1 local interpolate returns istatus = ', istatus
- if (istatus /= 0) then
- ! 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. FIXME.
- if (istatus /= 0 .and. .not. any(values == MISSING_R8)) then
- write(string1,*) 'interp routine returned a bad status but good values'
- call error_handler(E_ERR,'model_interpolate',string1,source,revision,revdate)
- endif
- if (istatus == 0 .and. any(values == MISSING_R8)) then
- write(string1,*) 'interp routine returned a good status but bad values'
- call error_handler(E_ERR,'model_interpolate',string1,source,revision,revdate)
- endif
- interp_val = MISSING_R8
- return
- endif
-
- ! compute sensible temperature as return value
- interp_val = theta_to_tk(values(1), values(2), values(3))
- return
-endif
-
-!print *, '2 local interpolate returns istatus = ', istatus
-
-! add other cases here for kinds we want to handle
-! if (istatus == 88) then
-! could try different interpolating different kinds here.
-! if (obs_type == KIND_xxx) then
-! code goes here
-! endif
-! endif
-
-! shouldn't get here.
-interp_val = MISSING_R8
-istatus = 100
-
-end subroutine model_interpolate
-
-
-!------------------------------------------------------------------
-
-subroutine local_interpolate(x, location, num_kinds, obs_kinds, interp_vals, istatus)
-
-! THIS VERSION IS ONLY CALLED INTERNALLY, so we can mess with the arguments.
-! For a given lat, lon, and height, interpolate the correct state values
-! to that location for the filter from the model state vectors. this version
-! can take multiple kinds at a single location. there is only a single return
-! istatus code - 0 is all interpolations worked, positive means one or more
-! failed.
-!
! ERROR codes:
!
! ISTATUS = 99: general error in case something terrible goes wrong...
@@ -879,299 +815,156 @@
! finding an applicable case.
!
-! Passed variables
+! passed variables
real(r8), intent(in) :: x(:)
type(location_type), intent(in) :: location
-integer, intent(in) :: num_kinds
-integer, intent(in) :: obs_kinds(:)
-real(r8), intent(out) :: interp_vals(:)
+integer, intent(in) :: obs_type
+real(r8), intent(out) :: interp_val
integer, intent(out) :: istatus
-! Local storage
+! local storage
-real(r8) :: loc_array(3), llon, llat, lheight, fract, v_interp
-integer :: i, j, ivar(1), ier, lower, upper, this_kind
-integer :: pt_base_offset, density_base_offset, qv_base_offset
+integer :: ivar, ier, obs_kind
+integer :: tvars(3)
+logical :: badkind
+real(r8) :: values(3), loc_array(3), lpres
+type(location_type) :: new_location
-real(r8) :: weights(3), lower_interp, upper_interp, ltemp, utemp, values(3)
-integer :: tri_indices(3), base_offset(num_kinds), tvars(3)
-
if ( .not. module_initialized ) call static_init_model
-! Assume failure. Set return val to missing, then the code can
-! just set istatus to something indicating why it failed, and return.
-! If the interpolation is good, the interp_val will be set to the
-! good value. Make any error codes set here be in the 10s
+interp_val = MISSING_R8
+istatus = 99 ! must be positive (and integer)
-interp_vals(:) = MISSING_R8 ! the DART bad value flag
-istatus = 99 ! unknown error
+! rename for sanity - we can't change the argument names
+! to this subroutine, but this really is a kind.
+obs_kind = obs_type
-! see if all observation kinds are in the state vector. this sets an
-! error code and returns without a fatal error if any answer is no.
-! exceptions: if the vertical location is specified in pressure
-! we end up computing the sensible temperature as part of converting
-! to height (meters), so if the kind to be computed is sensible temp
-! (which is not in state vector; potential temp is), go ahead and
-! say yes, we know how to compute it. if the vert is any other units
-! then fail and let the calling code call in for the components
-! it needs to compute sensible temp itself. also, if the obs is a
-! wind (or ocean current) component, and the edge normal speeds are
-! in the state vector, we can do it.
+! 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 (highest_obs_pressure_mb > 0.0) then
+ ! Reject obs above a user specified pressure level
+ if(vert_is_pressure(location)) then
+ loc_array = get_location(location)
+ lpres = loc_array(3)
-oktointerp: do i=1, num_kinds
- ivar(1) = get_progvar_index_from_kind(obs_kinds(i))
- if (ivar(1) <= 0) then
- ! exceptions 1 and 2:
- ! FIXME: the new code cannot yet compute sensible temperature in one go. fail for that.
- ! remove the new code test once we add it.
- !if ((obs_kinds(i) == KIND_TEMPERATURE) .and. vert_is_pressure(location) .and. &
- ! .not. use_new_code) cycle oktointerp
- if ((obs_kinds(i) == KIND_TEMPERATURE) .and. vert_is_pressure(location)) cycle oktointerp
- if ((obs_kinds(i) == KIND_U_WIND_COMPONENT) .or. obs_kinds(i) == KIND_V_WIND_COMPONENT) then
- ivar(1) = get_progvar_index_from_kind(KIND_EDGE_NORMAL_SPEED)
- if (ivar(1) > 0 .and. use_u_for_wind) cycle oktointerp
+ else if(vert_is_height(location) .or. vert_is_level(location)) then
+ new_location = location
+ call vert_convert(ens_mean, new_location, obs_kind, VERTISPRESSURE, istatus)
+ if(istatus == 0) then
+ loc_array = get_location(new_location)
+ lpres = loc_array(3)
+ else
+ interp_val = MISSING_R8
+ istatus = 200
+ return
endif
- istatus = 88 ! this kind not in state vector
- return
+ else ! VERTISUNDEF or VERTISSURFACE
+ lpres = 200100. ! VERTISUNDEF should impact entire column
endif
-enddo oktointerp
-
-! Not prepared to do w interpolation at this time
-do i=1, num_kinds
- if(obs_kinds(i) == KIND_VERTICAL_VELOCITY) then
- istatus = 16
- return
+
+ if (lpres < highest_obs_pressure_mb * 100.0_r8) then
+ ! Exclude from assimilation the obs above a user specified level
+ interp_val = MISSING_R8
+ istatus = 201
+ return
endif
-enddo
+endif
-! Get the individual locations values
-loc_array = get_location(location)
-llon = loc_array(1)
-llat = loc_array(2)
-lheight = loc_array(3)
-if (debug > 5) print *, 'requesting interpolation at ', llon, llat, lheight
+! see if observation kind 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.
+badkind = .false.
+ivar = get_progvar_index_from_kind(obs_kind)
-! if (debug > 2) print *, 'base offset now ', base_offset(1)
+if (ivar <= 0) then
+ badkind = .true.
-! FIXME: this needs to play well with multiple types in the kinds list.
-! for now do only the first one. if 'u' is part of the state vector
-! and if they are asking for U/V components, use the new RBF code.
-! if u isn't in the state vector, default to trying to interpolate
-! in the reconstructed U/V components at the cell centers.
+ ! exceptions 1 and 2:
+ if (obs_kind == KIND_TEMPERATURE) badkind = .false.
-if (use_new_code) then
- kindloop: do i=1, num_kinds
- if ((obs_kinds(i) == KIND_U_WIND_COMPONENT .or. &
- obs_kinds(i) == KIND_V_WIND_COMPONENT) .and. has_real_u .and. use_u_for_wind) then
- if (obs_kinds(i) == KIND_U_WIND_COMPONENT) then
- call compute_u_with_rbf(x, location, .TRUE., interp_vals(i), istatus)
- else
- call compute_u_with_rbf(x, location, .FALSE., interp_vals(i), istatus)
- endif
+ if ((obs_kind == KIND_U_WIND_COMPONENT) .or. obs_kind == KIND_V_WIND_COMPONENT) then
+ ivar = get_progvar_index_from_kind(KIND_EDGE_NORMAL_SPEED)
+ if (ivar > 0 .and. use_u_for_wind) badkind = .false.
+ endif
+endif
- else if (obs_kinds(i) /= KIND_TEMPERATURE) then
- ivar(1) = get_progvar_index_from_kind(obs_kinds(i))
- call compute_scalar_with_barycentric(x, location, 1, ivar, interp_vals(i:i), istatus)
- if (istatus /= 0) return
-
- else if (obs_kinds(i) == KIND_TEMPERATURE) then
- ! need to get potential temp, pressure, qv here, but can
- ! use same weights. push this down into the scalar code for now.
- tvars(1) = get_progvar_index_from_kind(KIND_POTENTIAL_TEMPERATURE)
- tvars(2) = get_progvar_index_from_kind(KIND_DENSITY)
- tvars(3) = get_progvar_index_from_kind(KIND_VAPOR_MIXING_RATIO)
-
- call compute_scalar_with_barycentric(x, location, 3, tvars, values, istatus)
- if (istatus /= 0) return
-
- interp_vals(i) = theta_to_tk(values(1), values(2), values(3))
- endif
- enddo kindloop
+if (badkind) then
+ istatus = 88 ! this kind not in state vector
return
endif
-!! start of original code - should be unneeded if we decide to
-!! use the new code.
-if (debug > 4) print *, 'using original interpolation code'
-
-! Find the start and end offsets for these fields in the state vector x(:)
-do i=1, num_kinds
- if (obs_kinds(i) == KIND_TEMPERATURE) then
- base_offset(i) = 0 ! this won't be used in this case
- else
- call get_index_range(obs_kinds(i), base_offset(i))
- endif
-enddo
-
-! Find the indices of the three cell centers that surround this point in
-! the horizontal along with the barycentric weights.
-call get_cell_indices(llon, llat, tri_indices, weights, ier)
-if (debug > 5) write(*, *) 'tri_inds ', tri_indices
-if (debug > 5) write(*, *) 'weights ', weights
-
-! If istatus is not zero couldn't find a triangle, fail
-if(ier /= 0) then
- istatus = 11
+! Not prepared to do w interpolation at this time
+if(obs_kind == KIND_VERTICAL_VELOCITY) then
+ istatus = 16
return
endif
-! FIXME: cannot do both surface pressure and any other 3d field with
-! this code (yet). if num kinds > 1 and any are not surface pressure, this
-! should fail.
-! Surface pressure is a 2D field
-if(obs_kinds(1) == KIND_SURFACE_PRESSURE) then
- lheight = 1.0_r8
- ! the 1 below is max number of vert levels to examine
- call triangle_interp(x, base_offset(1), tri_indices, weights, &
- nint(lheight), 1, interp_vals(1), ier)
- if(ier /= 0) then
- istatus = 13
+if ((obs_kind == KIND_U_WIND_COMPONENT .or. &
+ obs_kind == KIND_V_WIND_COMPONENT) .and. has_real_u .and. use_u_for_wind) then
+ if (obs_kind == KIND_U_WIND_COMPONENT) then
+ ! return U
+ call compute_u_with_rbf(x, location, .TRUE., interp_val, istatus)
else
- istatus = 0
+ ! return V
+ call compute_u_with_rbf(x, location, .FALSE., interp_val, istatus)
endif
- return
-endif
+!print *, 'called u_with_rbf, kind, istatus: ', obs_kind, istatus
+ if (istatus /= 0) return
-! the code below has vert undef returning an error, and vert surface
-! returning level 1 (which i believe is correct - the vertical grid
-! is terrain following). FIXME:
-! vert undef could return anything in the column in theory, right?
+else if (obs_kind == KIND_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(KIND_POTENTIAL_TEMPERATURE)
+ tvars(2) = get_progvar_index_from_kind(KIND_DENSITY)
+ tvars(3) = get_progvar_index_from_kind(KIND_VAPOR_MIXING_RATIO)
-if (vert_is_undef(location)) then
- istatus = 12
- return
-endif
+ call compute_scalar_with_barycentric(x, location, 3, tvars, values, istatus)
+ if (istatus /= 0) return
+
+ ! convert pot_temp, density, vapor mixing ratio into sensible temperature
+ interp_val = theta_to_tk(values(1), values(2), values(3))
-if(vert_is_surface(location)) then
- do j=1, num_kinds
- lheight = 1.0_r8 ! first grid level is surface
- call triangle_interp(x, base_offset(j), tri_indices, weights, &
- nint(lheight), nVertLevels, interp_vals(j), ier)
- if(ier /= 0) then
- istatus = 13
- return
- endif
- enddo
- istatus = 0
- return
-endif
+!print *, 'called compute temp, kind, istatus: ', obs_kind, istatus
+else
-! If vertical is on a model level don't need vertical interpolation either
-! (FIXME: since the vertical value is a real, in the wrf model_mod they
-! support non-integer levels and interpolate in the vertical just like
-! any other vertical type. we could easily add that here.)
-if(vert_is_level(location)) then
- ! FIXME: if this is W, the top is nVertLevels+1
- if (lheight > nVertLevels) then
- istatus = 12
- return
- endif
- do j=1, num_kinds
- call triangle_interp(x, base_offset(j), tri_indices, weights, &
- nint(lheight), nVertLevels, interp_vals(j), ier)
- if(ier /= 0) then
- istatus = 13
- return
- endif
- enddo
- istatus = 0
- return
-endif
+ ! direct interpolation, kind is in the state vector
+ tvars(1) = ivar
+ call compute_scalar_with_barycentric(x, location, 1, tvars, values, istatus)
+ interp_val = values(1)
+ if (istatus /= 0) return
+!print *, 'called compute_w_bary, kind, istatus: ', obs_kind, istatus
-! Vertical interpolation for pressure coordinates
- if(vert_is_pressure(location) ) then
- ! Need to get base offsets for the potential temperature, density, and water
- ! vapor mixing fields in the state vector
- call get_index_range(KIND_POTENTIAL_TEMPERATURE, pt_base_offset)
- call get_index_range(KIND_DENSITY, density_base_offset)
- call get_index_range(KIND_VAPOR_MIXING_RATIO, qv_base_offset)
-!print *, 'bases: t/rho/v = ', pt_base_offset, density_base_offset, qv_base_offset
- call find_pressure_bounds(x, lheight, tri_indices, weights, nVertLevels, &
- pt_base_offset, density_base_offset, qv_base_offset, lower, upper, fract, &
- ltemp, utemp, ier)
- if(ier /= 0) then
- interp_vals = MISSING_R8
- istatus = 17
- return
- endif
- ! if any of the input kinds are sensible temperature, we had to compute
- ! that value already to convert the pressure to the model levels, so just
- ! interpolate in the vertical and we're done with that one.
- do j=1, num_kinds
- if (obs_kinds(j) == KIND_TEMPERATURE) then
- ! Already have both values, interpolate in the vertical here.
- if (ltemp == MISSING_R8 .or. utemp == MISSING_R8) then
- istatus = 12
- return
- endif
- interp_vals(j) = (1.0_r8 - fract) * ltemp + fract * utemp
- else
- ! Next interpolate the observed quantity to the horizontal point at both levels
- call triangle_interp(x, base_offset(j), tri_indices, weights, lower, nVertLevels, &
- lower_interp, ier)
- if(ier /= 0) then
- istatus = 13
- return
- endif
- call triangle_interp(x, base_offset(j), tri_indices, weights, upper, nVertLevels, &
- upper_interp, ier)
- if(ier /= 0) then
- istatus = 13
- return
- endif
-
- ! Got both values, interpolate and return
- if (lower_interp == MISSING_R8 .or. upper_interp == MISSING_R8) then
- istatus = 12
- return
- endif
- interp_vals(j) = (1.0_r8 - fract) * lower_interp + fract * upper_interp
- endif
- enddo
- istatus = 0
- return
endif
-! in this section, unlike the others, we loop 3 times adding successive
-! contributions to the return value. if there's an error the
-! result value has been set to 0 so we have to reset it to the
-! missing value instead of only setting the istatus code.
-if(vert_is_height(location)) then
- ! For height, can do simple vertical search for interpolation for now
- ! Get the lower and upper bounds and fraction for each column
- interp_vals(:) = 0.0_r8
- do i = 1, 3
- call find_height_bounds(lheight, nVertLevels, zGridCenter(:, tri_indices(i)), &
- lower, upper, fract, ier)
- if(ier /= 0) then
- istatus = 12
- interp_vals(:) = MISSING_R8
- return
- endif
- do j=1, num_kinds
- call vert_interp(x, base_offset(j), tri_indices(i), nVertLevels, lower, fract, v_interp, ier)
- if(ier /= 0) then
- istatus = 13
- interp_vals(:) = MISSING_R8
- return
- endif
- interp_vals(j) = interp_vals(j) + weights(i) * v_interp
- enddo
- enddo
- istatus = 0
- return
+! 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. FIXME.
+if (istatus /= 0 .and. interp_val /= MISSING_R8) then
+ write(string1,*) 'interp routine returned a bad status but not a MISSING_R8 value'
+ write(string2,*) 'value = ', interp_val, ' istatus = ', istatus
+ call error_handler(E_ERR,'model_interpolate',string1,source,revision,revdate, &
+ text2=string2)
endif
+if (istatus == 0 .and. interp_val == MISSING_R8) then
+ write(string1,*) 'interp routine returned a good status but set value to MISSING_R8'
+ call error_handler(E_ERR,'model_interpolate',string1,source,revision,revdate)
+endif
-! shouldn't get here.
-interp_vals(:) = MISSING_R8
-istatus = 101
+call write_location(0,location,charstring=string1)
+if (debug > 5) print *, 'x(1), kind, loc ', x(1), obs_kind, trim(string1)
+if (debug > 5) print *, 'returning value ', interp_val
-end subroutine local_interpolate
+end subroutine model_interpolate
!------------------------------------------------------------------
@@ -1756,7 +1549,7 @@
where(dimIDs == TimeDimID) mystart = timeindex
where(dimIDs == TimeDimID) mycount = 1
- if ( debug > 1 ) then
+ if ( debug > 9 ) then
write(*,*)'nc_write_model_vars '//trim(varname)//' start is ',mystart(1:ncNdims)
write(*,*)'nc_write_model_vars '//trim(varname)//' count is ',mycount(1:ncNdims)
endif
@@ -1882,6 +1675,7 @@
if ( .not. module_initialized ) call static_init_model
ens_mean = filter_ens_mean
+!print *, 'setting ens_mean, x(100) = ', ens_mean(100)
end subroutine ens_mean_for_model
@@ -1995,16 +1789,24 @@
subroutine get_close_obs(gc, base_obs_loc, base_obs_kind, &
obs_loc, obs_kind, num_close, close_ind, dist)
!
-! FIXME ... not tested.
-!
! Given a DART location (referred to as "base") and a set of candidate
! locations & kinds (obs, obs_kind), returns the subset close to the
! "base", their indices, and their distances to the "base" ...
! For vertical distance computations, general philosophy is to convert all
-! vertical coordinates to a common coordinate. This coordinate type is defined
-! in the namelist with the variable "vert_localization_coord".
+! vertical coordinates to a common coordinate. This coordinate type can be
+! defined in the namelist with the variable "vert_localization_coord".
+! But we first try a single coordinate type as the model level here.
+! FIXME: We need to add more options later.
+! Vertical conversion is carried out by the subroutine vert_convert.
+
+! Note that both base_obs_loc and obs_loc are intent(inout), meaning that these
+! locations are possibly modified here and returned as such to the calling routine.
+! The calling routine is always filter_assim and these arrays are local arrays
+! within filter_assim. In other words, these modifications will only matter within
+! filter_assim, but will not propagate backwards to filter.
+
type(get_close_type), intent(in) :: gc
type(location_type), intent(inout) :: base_obs_loc
integer, intent(in) :: base_obs_kind
@@ -2014,31 +1816,40 @@
integer, dimension(:), intent(out) :: close_ind
real(r8), dimension(:), intent(out) :: dist
-integer :: t_ind, istatus1, istatus2, k
+integer :: ztypeout
+integer :: t_ind, istatus1, istatus2, k, iv
integer :: base_which, local_obs_which
-real(r8), dimension(3) :: base_array, local_obs_array
+real(r8), dimension(3) :: base_xyz, local_obs_xyz
+real(r8), dimension(3) :: test_xyz
type(location_type) :: local_obs_loc
+real(r8) :: hor_dist
+hor_dist = 1.0e9_r8
+
! Initialize variables to missing status
num_close = 0
close_ind = -99
-dist = 1.0e9 !something big and positive (far away) in radians
+dist = 1.0e9_r8 !something big and positive (far away) in radians
istatus1 = 0
istatus2 = 0
! Convert base_obs vertical coordinate to requested vertical coordinate if necessary
-base_array = get_location(base_obs_loc)
+base_xyz = get_location(base_obs_loc)
base_which = nint(query_location(base_obs_loc))
-! fixme ...
+ztypeout = vert_localization_coord
+
if (.not. horiz_dist_only) then
-! if (base_which /= wrf%dom(1)%vert_coord) then
-! call vert_interpolate(ens_mean, base_obs_loc, base_obs_kind, istatus1)
-! elseif (base_array(3) == MISSING_R8) then
-! istatus1 = 1
-! endif
+ if (base_xyz(3) == MISSING_R8) then
+ istatus1 = 1
+ else if (base_which /= vert_localization_coord) then
+ call vert_convert(ens_mean, base_obs_loc, base_obs_kind, ztypeout, istatus1)
+ call write_location(0,base_obs_loc,charstring=string1)
+ if(debug > 5) &
+ call error_handler(E_MSG, 'get_close_obs: base_obs_loc',string1,source, revision, revdate)
+ endif
endif
if (istatus1 == 0) then
@@ -2060,22 +1871,31 @@
! 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 (.not. horiz_dist_only) then
- !fixme if (local_obs_which /= wrf%dom(1)%vert_coord) then
- !fixme call vert_interpolate(ens_mean, local_obs_loc, obs_kind(t_ind), istatus2)
- ! Store the "new" location into the original full local array
- obs_loc(t_ind) = local_obs_loc
- !fixme endif
+ if (local_obs_which /= vert_localization_coord) then
+ call vert_convert(ens_mean, local_obs_loc, obs_kind(t_ind), ztypeout, istatus2)
+ else
+ istatus2 = 0
+ endif
endif
! Compute distance - set distance to a very large value if vert coordinate is missing
! or vert_interpolate returned error (istatus2=1)
- local_obs_array = get_location(local_obs_loc)
- if (( (.not. horiz_dist_only) .and. &
- (local_obs_array(3) == MISSING_R8)) .or. &
- (istatus2 == 1) ) then
- dist(k) = 1.0e9
+ local_obs_xyz = get_location(local_obs_loc)
+ if (( (.not. horiz_dist_only) .and. &
+ (local_obs_xyz(3) == MISSING_R8)) .or. &
+ (istatus2 /= 0) ) then
+ dist(k) = 1.0e9_r8
else
- dist(k) = get_dist(base_obs_loc, local_obs_loc, base_obs_kind, obs_kind(t_ind))
+ dist(k) = get_dist(base_obs_loc, local_obs_loc, base_obs_kind, obs_kind(t_ind))
+! if (debug > 4 .and. k < 100 .and. my_task_id() == 0) then
+! print *, 'calling get_dist'
+! call write_location(0,base_obs_loc,charstring=string2)
+! call error_handler(E_MSG, 'get_close_obs: base_obs_loc',string2,source, revision, revdate)
+! call write_location(0,local_obs_loc,charstring=string2)
+! call error_handler(E_MSG, 'get_close_obs: local_obs_loc',string2,source, revision, revdate)
+! hor_dist = get_dist(base_obs_loc, local_obs_loc, base_obs_kind, obs_kind(t_ind), no_vert=.true.)
+! print *, 'hor/3d_dist for k =', k, ' is ', hor_dist,dist(k)
+! endif
endif
enddo
@@ -2295,7 +2115,7 @@
where(dimIDs == TimeDimID) mystart = TimeDimLength ! pick the latest time
where(dimIDs == TimeDimID) mycount = 1 ! only use one time
- if ( debug > 1 ) then
+ if ( debug > 9 ) then
write(*,*)'analysis_file_to_statevector '//trim(varname)//' start = ',mystart(1:ncNdims)
write(*,*)'analysis_file_to_statevector '//trim(varname)//' count = ',mycount(1:ncNdims)
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list