[Dart-dev] [3657] DART/trunk: Major update to the radar observation
module, along with matching documentation
nancy at ucar.edu
nancy at ucar.edu
Tue Nov 18 14:35:35 MST 2008
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20081118/43bb6f33/attachment-0001.html
-------------- next part --------------
Modified: DART/trunk/models/wrf/model_mod.f90
===================================================================
--- DART/trunk/models/wrf/model_mod.f90 2008-11-14 23:27:05 UTC (rev 3656)
+++ DART/trunk/models/wrf/model_mod.f90 2008-11-18 21:35:34 UTC (rev 3657)
@@ -1174,13 +1174,18 @@
! x: Full DART state vector relevant to what's being updated
! in the filter (mean or individual members).
-! istatus: Returned 0 if everything is OK, 1 if error occured.
-! -1 if the station height is lower than the lowest model level
-! while the station is located inside the horizontal model domain.
+! istatus: Returned 0 if everything is OK, >0 if error occured.
+! 1 = missing domain id
+! 2 = bad vertical level
+! 3 = unsupported obs kind
+! 10 = polar observation while restrict_polar namelist true
+! 99 = unknown reason (reached end of routine with missing_r8
+! as obs_val)
! modified 26 June 2006 to accomodate vortex attributes
! modified 13 December 2006 to accomodate changes for the mpi version
! modified 22 October 2007 to accomodate global WRF (3.0)
+! modified 18 November 2008 to allow unknown kinds to return without stopping
! arguments
real(r8), intent(in) :: x(:)
@@ -1400,7 +1405,7 @@
! Deal with missing vertical coordinates -- return with istatus .ne. 0
if(zloc == missing_r8) then
obs_val = missing_r8
- istatus = 1
+ istatus = 2
deallocate(v_h, v_p)
return
endif
@@ -2505,12 +2510,14 @@
!-----------------------------------------------------
! If obs_kind is not negative (for identity obs), or if it is not one of the above 15
- ! explicitly checked-for kinds, then return error message.
+ ! explicitly checked-for kinds, then set error istatus and missing_r8.
else
- write(errstring,*)'Obs kind not recognized for following kind: ',obs_kind
- call error_handler(E_ERR,'model_interpolate', errstring, &
- source, revision, revdate)
+ obs_val = missing_r8
+ istatus = 3
+ if (debug) print*, 'unrecognized obs KIND, value = ', obs_kind
+ deallocate(v_h, v_p)
+ return
end if
@@ -2581,8 +2588,11 @@
end if ! end of "if ( obs_kind < 0 )"
-! Now that we are done, check to see if a missing value somehow made it through
-if ( obs_val == missing_r8 ) istatus = 1
+! Now that we are done, check to see if a missing value somehow
+! made it through without already having set an error return code.
+if ( obs_val == missing_r8 .and. istatus == 0 ) then
+ istatus = 99
+endif
! Pring the observed value if in debug mode
if(debug) print*,'model_interpolate() return value= ',obs_val
Modified: DART/trunk/obs_def/obs_def_radar_mod.f90
===================================================================
--- DART/trunk/obs_def/obs_def_radar_mod.f90 2008-11-14 23:27:05 UTC (rev 3656)
+++ DART/trunk/obs_def/obs_def_radar_mod.f90 2008-11-18 21:35:34 UTC (rev 3657)
@@ -3,74 +3,107 @@
! University Corporation for Atmospheric Research
! Licensed under the GPL -- www.gpl.org/licenses/gpl.html
-! <next few lines under version control, do not edit>
-! $URL$
-! $Id$
-! $Revision$
-! $Date$
+!-----------------------------------------------------------------------------
+! DART radar observation module, including the observation operators for the
+! two primary radar-observation types -- Doppler velocity and reflectivity --
+! plus other utility subroutines and functions. A number of simplifications
+! are employed for the observation operators. Most notably, the model state
+! is mapped to a "point" observation, whereas a real radar observation is a
+! volumetric sample. The implications of this approximation have not been
+! investigated fully, so in the future it might be worth developing and testing
+! more sophisticated observation operators that produce volumetric power-
+! weighted samples.
+!
+! This module is able to compute reflectivity and precipitation fall speed
+! (needed for computing Doppler radial velocity) from the prognostic model
+! fields only for simple single-moment microphysics schemes such as the Kessler
+! and Lin schemes. If a more complicated microphysics scheme is used, then
+! reflectivity and fall speed must be accessible instead as diagnostic fields
+! in the model state.
+!
+! Author and Contact information:
+! Radar Science: David Dowell, ddowell at ucar.edu
+! DART Code: Nancy Collins, nancy at ucar.edu
+! Original DART/Radar work: Alain Caya
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
! BEGIN DART PREPROCESS KIND LIST
! DOPPLER_RADIAL_VELOCITY, KIND_VELOCITY
! RADAR_REFLECTIVITY, KIND_RADAR_REFLECTIVITY
+! PRECIPITATION_FALL_SPEED, KIND_POWER_WEIGHTED_FALL_SPEED
! END DART PREPROCESS KIND LIST
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
! BEGIN DART PREPROCESS USE OF SPECIAL OBS_DEF MODULE
-! use obs_def_radar_mod, only : write_rad_vel, read_rad_vel, &
-! interactive_rad_vel, get_expected_rad_vel, &
-! write_rad_ref, read_rad_ref, &
-! interactive_rad_ref, get_expected_rad_ref
+! use obs_def_radar_mod, only : write_radial_vel, read_radial_vel, &
+! interactive_radial_vel, get_expected_radial_vel, &
+! read_radar_ref, get_expected_radar_ref
! END DART PREPROCESS USE OF SPECIAL OBS_DEF MODULE
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
! BEGIN DART PREPROCESS GET_EXPECTED_OBS_FROM_DEF
-! case(DOPPLER_RADIAL_VELOCITY)
-! call get_expected_rad_vel(state, location, obs_def%key, obs_val, istatus)
-! case(RADAR_REFLECTIVITY)
-! call get_expected_rad_ref(state, location, obs_val, istatus)
+! case(DOPPLER_RADIAL_VELOCITY)
+! call get_expected_radial_vel(state, location, obs_def%key, obs_val, istatus)
+! case(RADAR_REFLECTIVITY)
+! call get_expected_radar_ref(state, location, obs_val, istatus)
! END DART PREPROCESS GET_EXPECTED_OBS_FROM_DEF
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
! BEGIN DART PREPROCESS READ_OBS_DEF
-! case(DOPPLER_RADIAL_VELOCITY)
-! call read_rad_vel(obs_def%key, ifile, fileformat)
-! case(RADAR_REFLECTIVITY)
-! call read_rad_ref(obs_val, obs_def%key, ifile, fileformat)
+! case(DOPPLER_RADIAL_VELOCITY)
+! call read_radial_vel(obs_def%key, ifile, fileformat)
+! case(RADAR_REFLECTIVITY)
+! call read_radar_ref(obs_val, obs_def%key)
! END DART PREPROCESS READ_OBS_DEF
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
! BEGIN DART PREPROCESS WRITE_OBS_DEF
-! case(DOPPLER_RADIAL_VELOCITY)
-! call write_rad_vel(obs_def%key, ifile, fileformat)
-! case(RADAR_REFLECTIVITY)
-! call write_rad_ref(obs_def%key, ifile, fileformat)
+! case(DOPPLER_RADIAL_VELOCITY)
+! call write_radial_vel(obs_def%key, ifile, fileformat)
+! case(RADAR_REFLECTIVITY)
+! continue
! END DART PREPROCESS WRITE_OBS_DEF
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
! BEGIN DART PREPROCESS INTERACTIVE_OBS_DEF
-! case(DOPPLER_RADIAL_VELOCITY)
-! call interactive_rad_vel(obs_def%key)
-! case(RADAR_REFLECTIVITY)
-! call interactive_rad_ref(obs_def%key)
+! case(DOPPLER_RADIAL_VELOCITY)
+! call interactive_radial_vel(obs_def%key)
+! case(RADAR_REFLECTIVITY)
+! continue
! END DART PREPROCESS INTERACTIVE_OBS_DEF
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
! BEGIN DART PREPROCESS MODULE CODE
module obs_def_radar_mod
-use types_mod, only : r8, missing_r8, ps0, PI, gravity, DEG2RAD
+use types_mod, only : r8, missing_r8, PI, deg2rad
use utilities_mod, only : register_module, error_handler, E_ERR, E_MSG, &
- check_namelist_read, find_namelist_in_file, &
+ check_namelist_read, find_namelist_in_file, &
nmlfileunit, do_output
use location_mod, only : location_type, write_location, read_location, &
- interactive_location
+ interactive_location, get_location
use assim_model_mod, only : interpolate
use obs_kind_mod, only : KIND_U_WIND_COMPONENT, KIND_V_WIND_COMPONENT, &
- KIND_TEMPERATURE, KIND_VERTICAL_VELOCITY, &
- KIND_RAINWATER_MIXING_RATIO, KIND_DENSITY, &
- KIND_GRAUPEL_MIXING_RATIO, KIND_SNOW_MIXING_RATIO
+ KIND_TEMPERATURE, KIND_VERTICAL_VELOCITY, &
+ KIND_RAINWATER_MIXING_RATIO, KIND_DENSITY, &
+ KIND_GRAUPEL_MIXING_RATIO, &
+ KIND_SNOW_MIXING_RATIO, &
+ KIND_POWER_WEIGHTED_FALL_SPEED, &
+ KIND_RADAR_REFLECTIVITY
implicit none
private
-public :: write_rad_vel, read_rad_vel, set_rad_vel, interactive_rad_vel, &
- write_rad_ref, read_rad_ref, interactive_rad_ref, &
- get_expected_rad_vel, get_expected_rad_ref, &
- get_obs_def_rad_ref, get_obs_def_rad_vel
+public :: read_radar_ref, get_expected_radar_ref, &
+ read_radial_vel, write_radial_vel, interactive_radial_vel, &
+ get_expected_radial_vel, get_obs_def_radial_vel, set_radial_vel
! version controlled file description for error handling, do not edit
character(len=128), parameter :: &
@@ -78,94 +111,235 @@
revision = "$Revision$", &
revdate = "$Date$"
-logical, save :: module_initialized = .false.
+logical :: module_initialized = .false.
-! Storage for the special information required for velocity observations of this type
-integer, parameter :: max_rad_vel_obs = 1000000
-type(location_type) :: rad_loc(max_rad_vel_obs)
-real(r8) :: direction(3,max_rad_vel_obs)
-real(r8) :: nyquistvel(max_rad_vel_obs)
+! Derived type for radial velocity. Contains auxiliary information stored
+! with each obs of this type; used to compute the forward operator.
+! See more extensive comments in the interactive_radial_vel() routine for
+! expected units, etc. Technically, the radar location is unused in the
+! forward operators currently in this code, but it may be useful for post
+! processing or diagnostics, especially if multiple radar locations are
+! in the same file.
+type radial_vel_type
+ private
+ type(location_type) :: radar_location ! location of radar
+ real(r8) :: beam_direction(3) ! direction of beam
+ real(r8) :: nyquist_velocity ! nyquist velocity
+end type radial_vel_type
-! Storage for the special information required for reflectivity observations of this type
-integer, parameter :: max_rad_ref_obs = 1000000
-!define here :: globally_defined_reflectivity_metadata_variable(max_rad_ref_obs)
+! Cumulative index into radial velocity metadata array
+integer :: velkeycount = 0
-integer :: velkeycount = 0 ! cumulative index into rad. velocity metadata
-integer :: refkeycount = 0 ! cumulative index into rad. reflectivity metadata
+! For error message content
+character(len=129) :: msgstring
-real(r8), parameter :: dief = 0.224_r8
+! Values which are initialized at run time so some can be changed by
+! namelist. After initialization, treated as parameters (values not changed).
-real(r8), parameter :: n0r = 8.0e6_r8, n0g = 4.0e4_r8, n0s = 3.0e6_r8
-real(r8), parameter :: rho_r = 1000.0_r8, rho_g = 917.0_r8, rho_s = 100.0_r8
+! Note that the value of gravity is hardcoded here. The value for gravity
+! used in the model should match this value. See additional comments
+! below in the initialize_constants() subroutine.
+real(r8) :: param_gravity ! gravitational acceleration (m s^-2)
+
+ ! empirical constants:
+real(r8) :: param_a ! 10^^18 * a (const in rain fall speed eqn)
+real(r8) :: param_b ! b (const in rain fall speed eqn)
+real(r8) :: param_c ! 10^^18 * c (const in snow fall speed eqn)
+real(r8) :: param_d ! d (const in snow fall speed eqn)
+
+real(r8) :: param_CD ! drag coefficient for graupel/hail
+real(r8) :: param_rhos0 ! reference density (kg m^-3) in emperical
+ ! dropsize-fall speed eqn
+real(r8) :: param_e ! parameter in graupel/hail fall speed eqn
+
+ ! results of gamma function applied to:
+real(r8) :: param_gam7b ! (7+b)
+real(r8) :: param_gam7d ! (7+d)
+real(r8) :: param_gam7f ! (7+0.5)
+
+ ! exponents in equations for:
+real(r8) :: param_powr ! rain fall speed
+real(r8) :: param_pows ! snow fall speed
+real(r8) :: param_powg_dry ! dry graupel/hail fall speed
+real(r8) :: param_powg_wet ! wet graupel/hail fall speed
+
+ ! parameters in the equations for:
+real(r8) :: param_fs_r ! rain fall speed
+real(r8) :: param_fs_wet_s ! wet snow fall speed
+real(r8) :: param_fs_dry_s ! dry snow fall speed
+real(r8) :: param_fs_wet_g ! wet graupel/hail fall speed
+real(r8) :: param_fs_dry_g ! dry graupel/hail fall speed
+
+ ! parameters in the equations for:
+real(r8) :: param_refl_r ! reflectivity from rain
+real(r8) :: param_refl_wet_s ! reflectivity from wet snow
+real(r8) :: param_refl_dry_s ! reflectivity from dry snow
+real(r8) :: param_refl_wet_g ! reflectivity from wet graupel/hail
+real(r8) :: param_refl_dry_g ! reflectivity from dry graupel/hail
+
!-------------------------------------------------------------
! Namelist with default values
!
-! convert_to_dbz: Convert final reflectivity to dB?
-! If .TRUE., then make sure dbz_threshold is specified correctly.
-! (Note that for now, this is reset to .TRUE. in initialize_module
-! as it is not clear how to assimilate Z.)
+! Obsolete: convert_to_dbz and dbz_threshold
+! convert_to_dbz and dbz_threshold have both been removed from the namelist.
+! Values will always be converted to dBZ, and threshold was only used to
+! ensure the log() call never saw a real 0.0_r8. Please remove these
+! values from your namelist to avoid a run-time error.
!
-! dbz_threshold: Only applied if convert_to_dbz is .TRUE.
-! Reflectivity below dbz_threshold is set to dbz_threshold, which is chosen
-! such to convert to the smallest possible dBZ.
-! The purpose is mainly computational: We want to avoid the occurence of log(0).
-! Some useful values: Z = 3.163 -> 5.0 dBZ
-! Z = 2.512 -> 4.0 dBZ
-! Z = 1.259 -> 1.0 dBZ
-! Z = 1.0233 -> 0.1 dBZ
-! Z = 0.1 -> -10.0 dBZ
-! Z = 0.01 -> -20.0 dBZ
-! Z = 0.001 -> -30.0 dBZ
+! There are two replicated sets of 3 namelist values below.
+! In each case, there is 1 logical value and 2 numeric values.
+! If the logical is false, the numeric values are ignored.
+! If true, then the 2 numeric values are:
+! 1) a threshold value. If the observation or forward operator value is
+! less than this threshold, it will be set to a different value.
+! 2) the value it should be set to.
+! The value is separate from the threshold to allow, for example, the option
+! of setting all values below -20 dBZ to -40 dBZ.
!
-! apply_ref_limit_to_obs: If .TRUE., replace reflectivities below "reflectivity_limit"
-! with "lowest_reflectivity". If FALSE, ignore "reflectivity_limit"
-! and "lowest_reflectivity".
+! If the observation value or the forward operator already has a value of
+! missing_r8 it is assumed either the istatus is marked as failed (for the
+! forward operator) or that the QC (quality control) flag is set to not
+! assimilate this observation and that value is left unchanged regardless of
+! the setting on the apply_ref_limit flag. Note however that it is not a
+! good idea to reset a good but small observation value to missing_r8 -- do
+! not use it as the lowest_reflectivity setting.
!
-! reflectivity_limit_obs: Reflectivity below this value is set to "lowest_reflectivity".
+! The next 3 namelist items apply to the incoming observation values.
+! They are in the namelist so they can be changed at runtime, instead
+! of set only when the observation file is originally generated.
!
-! lowest_reflectivity_obs: If reflectivity < reflectivity_limit, reflectivity is set to this value.
-! The default value of 'missing' is useful for perfect model experiments.
-! Suggested options: lowest_reflectivity_obs = missing_r8 or
-! lowest_reflectivity_obs = any real value.
+! apply_ref_limit_to_obs:
+! Logical. If .TRUE. replace all reflectivity values less than
+! "reflectivity_limit_obs" with "lowest_reflectivity_obs" value.
+! If .FALSE. leave all values as-is. Defaults to .FALSE.
!
-! apply_ref_limit_to_state: Similar to "apply_ref_limit_to_obs" but applied to the state
-! (forward operator).
+! reflectivity_limit_obs:
+! The threshold value. Observed reflectivity values less than this
+! threshold will be set to the "lowest_reflectivity_obs" value.
+! Units are dBZ. Defaults to -10.0.
!
-! reflectivity_limit_state: Similar to "reflectivity_limit_obs" for state.
+! lowest_reflectivity_obs:
+! The 'set-to' value. Observed reflectivity values less than the
+! threshold will be set to this value. Units are dBZ. Defaults to -10.0.
!
-! lowest_reflectivity_state: Similar to "lowest_reflectivity_obs" for state.
+! The next 3 namelist items apply to the forward operator values (the returned
+! value from each ensemble member predicting what the observation value should
+! be given the current state in this particular member).
+!
+! apply_ref_limit_to_fwd_op:
+! Same as "apply_ref_limit_to_obs", but for the forward operator.
+!
+! reflectivity_limit_fwd_op:
+! Same as "reflectivity_limit_obs", but for the forward operator values.
+!
+! lowest_reflectivity_fwd_op:
+! Same as "lowest_reflectivity_obs", but for the forward operator values.
+!
+! max_radial_vel_obs:
+! Integer value. Maximum number of observations of this type to support at
+! run time. This is combined total of all obs_seq files, for example the
+! observation diagnostic program potentially opens multiple obs_seq.final
+! files, or the obs merge program can also open multiple obs files.
+! Default is 1,000,000 observations.
+!
+! dielectric_factor:
+! According to Smith (1984), there are two choices for the dielectric
+! factor, depending on how the snow particle sizes are specified.
+! If melted raindrop diameters are used, then the factor is 0.224. If
+! equivalent ice sphere diameters are used, then the factor is 0.189.
+! Since the common convention is to use melted raindrop diameters, the
+! default here is 0.224.
+!
+! n0_rain, n0_graupel, n0_snow:
+! Intercept parameters (m^-4) for size distributions of each hydrometeor.
+! The defaults below are for the Lin et al. microphysics scheme
+! with the Hobbs settings for graupel/hail. (The Hobbs graupel settings
+! are also the default for the Lin scheme in WRF 2.2 and 3.0.)
+!
+! rho_rain, rho_graupel, rho_snow:
+! Density (kg m^-3) of each hydrometeor type. The defaults below are for the
+! Lin et al. microphysics scheme with the Hobbs setting for graupel/hail.
+!
+! allow_wet_graupel:
+! Logical. It is difficult to predict/diagnose whether graupel/hail has a
+! wet or dry surface. Even when the temperature is above freezing,
+! evaporation and/or absorption can still result in a dry surface. This
+! issue is important because the reflectivity from graupel with a wet surface
+! is significantly greater than that from graupel with a dry surface.
+! Currently, the user has two options for how to compute graupel
+! reflectivity. If allow_wet_graupel is .false. (the default), then graupel
+! is always assumed to be dry. If allow_wet_graupel is .true., then graupel
+! is assumed to be wet (dry) when the temperature is above (below) freezing.
+! A consequence is that a sharp gradient in reflectivity will be produced at
+! the freezing level. In the future, it might be better to provide the
+! option of having a transition layer.
+logical :: apply_ref_limit_to_obs = .false.
+real(r8) :: reflectivity_limit_obs = -10.0_r8
+real(r8) :: lowest_reflectivity_obs = -10.0_r8
+logical :: apply_ref_limit_to_fwd_op = .false.
+real(r8) :: reflectivity_limit_fwd_op = -10.0_r8
+real(r8) :: lowest_reflectivity_fwd_op = -10.0_r8
+integer :: max_radial_vel_obs = 1000000
+logical :: allow_wet_graupel = .false.
-logical :: convert_to_dbz = .true.
-real(r8) :: dbz_threshold = 0.001_r8
-logical :: apply_ref_limit_to_obs = .false.
-real(r8) :: reflectivity_limit_obs = 0.0_r8
-real(r8) :: lowest_reflectivity_obs = missing_r8
-logical :: apply_ref_limit_to_state = .false.
-real(r8) :: reflectivity_limit_state = 0.0_r8
-real(r8) :: lowest_reflectivity_state = missing_r8
-namelist /obs_def_radar_mod_nml/ convert_to_dbz, dbz_threshold, &
- apply_ref_limit_to_obs, reflectivity_limit_obs, &
- lowest_reflectivity_obs, &
- apply_ref_limit_to_state, reflectivity_limit_state, &
- lowest_reflectivity_state
+! Constants which depend on the microphysics scheme used. Should be set
+! in the namelist to match the case being run.
+real(r8) :: dielectric_factor = 0.224_r8
+real(r8) :: n0_rain = 8.0e6_r8
+real(r8) :: n0_graupel = 4.0e6_r8
+real(r8) :: n0_snow = 3.0e6_r8
+real(r8) :: rho_rain = 1000.0_r8
+real(r8) :: rho_graupel = 400.0_r8
+real(r8) :: rho_snow = 100.0_r8
+
+namelist /obs_def_radar_mod_nml/ apply_ref_limit_to_obs, &
+ reflectivity_limit_obs, &
+ lowest_reflectivity_obs, &
+ apply_ref_limit_to_fwd_op, &
+ reflectivity_limit_fwd_op, &
+ lowest_reflectivity_fwd_op, &
+ max_radial_vel_obs, &
+ allow_wet_graupel, &
+ dielectric_factor, &
+ n0_rain, &
+ n0_graupel, &
+ n0_snow, &
+ rho_rain, &
+ rho_graupel, &
+ rho_snow
+
+
+! Module global storage for auxiliary obs data, allocated in init routine
+type(radial_vel_type), allocatable :: radial_vel_data(:)
+
contains
!----------------------------------------------------------------------
+!----------------------------------------------------------------------
+! Start of executable routines
+!----------------------------------------------------------------------
+!----------------------------------------------------------------------
subroutine initialize_module
-implicit none
+! Called once to set values and allocate space
-integer :: iunit, io
+integer :: iunit, io, rc
+! Prevent multiple calls from executing this code more than once.
+if (module_initialized) return
+
+module_initialized = .true.
+
+! Log the version of this source file.
call register_module(source, revision, revdate)
-! Read the namelist entry
+! Read the namelist entry.
call find_namelist_in_file("input.nml", "obs_def_radar_mod_nml", iunit)
read(iunit, nml = obs_def_radar_mod_nml, iostat = io)
call check_namelist_read(iunit, io, "obs_def_radar_mod_nml")
@@ -174,662 +348,714 @@
if (do_output()) write(nmlfileunit, nml=obs_def_radar_mod_nml)
if (do_output()) write( * , nml=obs_def_radar_mod_nml)
-! This is temporary as it is not clear how to assimilate Z
-convert_to_dbz = .true.
+! Consistency warning; print a message if the thresholds and lower values
+! are going to be used and are different.
+call check_namelist_limits(apply_ref_limit_to_obs, &
+ reflectivity_limit_obs, &
+ lowest_reflectivity_obs, &
+ apply_ref_limit_to_fwd_op, &
+ reflectivity_limit_fwd_op, &
+ lowest_reflectivity_fwd_op)
-module_initialized = .true.
+! Allocate space for the auxiliary information associated with each obs
+! This code must be placed after reading the namelist, so the user can
+! increase or decrease the number of obs supported and use more or less
+! memory at run time.
+allocate(radial_vel_data(max_radial_vel_obs), stat = rc)
+if (rc /= 0) then
+ write(msgstring, *) 'initial allocation failed for radial vel obs data,', &
+ 'itemcount = ', max_radial_vel_obs
+ call error_handler(E_ERR,'initialize_module', msgstring, &
+ source, revision, revdate)
+endif
-end subroutine initialize_module
+! Set the module global values that do not change during the run.
+! This code uses some values which are set in the namelist, so this call
+! *must* happen after the namelist read above.
+call initialize_constants()
-!----------------------------------------------------------------------
+! Log the values used for the constants.
+call print_constants()
-subroutine write_rad_vel(velkey, ifile, fform)
+end subroutine initialize_module
-integer, intent(in) :: velkey, ifile
-character(len=*), intent(in), optional :: fform
-
-character(len=32) :: fileformat
-
-if ( .not. module_initialized ) call initialize_module
-
-fileformat = "ascii" ! supply default
-if(present(fform)) fileformat = trim(adjustl(fform))
-
-SELECT CASE (fileformat)
- CASE ("unf", "UNF", "unformatted", "UNFORMATTED")
- call write_location(ifile, rad_loc(velkey), fileformat)
- call write_orientation(ifile, direction(:,velkey), fileformat)
- call write_nyquistvel(ifile, nyquistvel(velkey), fileformat)
- ! Write out the obs_def velkey for this observation
- write(ifile) velkey
- CASE DEFAULT
- ! Write the 5 character identifier for verbose formatted output
- write(ifile, 11)
-11 format('platform')
- call write_location(ifile, rad_loc(velkey), fileformat)
- call write_orientation(ifile, direction(:,velkey), fileformat)
- call write_nyquistvel(ifile, nyquistvel(velkey), fileformat)
- ! Write out the obs_def velkey for this observation
- write(ifile, *) velkey
-END SELECT
-
-end subroutine write_rad_vel
-
!----------------------------------------------------------------------
-
-subroutine get_obs_def_rad_vel(velkey, radarloc, beamdir, velnyquist)
-
-integer, intent(in) :: velkey
-type(location_type), intent(out) :: radarloc
-real(r8), intent(out) :: beamdir(3)
-real(r8), intent(out) :: velnyquist
-
-radarloc = rad_loc(velkey)
-beamdir = direction(:,velkey)
-velnyquist = nyquistvel(velkey)
-
-end subroutine get_obs_def_rad_vel
-
!----------------------------------------------------------------------
+! Radial velocity section
+!----------------------------------------------------------------------
+!----------------------------------------------------------------------
-subroutine read_rad_vel(velkey, ifile, fform)
+subroutine read_radial_vel(velkey, ifile, fform)
-! Main read subroutine for the obs kind radial velocity metadata.
+! Main read subroutine for the radial velocity observation auxiliary data.
-! location: Refers to the lat/lon/height/vertical coordinate option
-! for the radar. Its type is defined through location_type in
-! location_mod. Uses the read_location function in
-! location_mod.
-!
-! orientation: This is a 3-element array specific to the radar module:
-! orientation(1) = sin(azimuth)*cos(elevation)
-! orientation(2) = cos(azimuth)*cos(elevation)
-! orientation(3) = sin(elevation)
-
integer, intent(out) :: velkey
integer, intent(in) :: ifile
character(len=*), intent(in), optional :: fform
character(len=8) :: header
-character(len=32) :: fileformat
-real(r8) :: orientation(3)
-type(location_type) :: location
-real(r8) :: nyquistv
+logical :: is_asciifile
+type(location_type) :: radar_location
+real(r8) :: beam_direction(3)
+real(r8) :: nyquist_velocity
+integer :: oldkey
+! radar_location: Refers to the lat/lon/height/vertical coordinate option
+! for the radar. Its type is defined through location_type in location_mod.
+! Uses the read_location function in location_mod.
+!
+! beam_direction: This is a 3-element array specific to the radar module:
+! beam_direction(1) = sin(azimuth)*cos(elevation)
+! beam_direction(2) = cos(azimuth)*cos(elevation)
+! beam_direction(3) = sin(elevation)
+
if ( .not. module_initialized ) call initialize_module
-fileformat = "ascii" ! supply default
-if(present(fform)) fileformat = trim(adjustl(fform))
+is_asciifile = ascii_file_format(fform)
-SELECT CASE (fileformat)
- CASE ("unf", "UNF", "unformatted", "UNFORMATTED")
- location = read_location(ifile, fileformat)
- orientation = read_orientation(ifile, fileformat)
- nyquistv = read_nyquistvel(ifile, fileformat)
- ! Read in the velkey for this particular observation
- read(ifile) velkey
- CASE DEFAULT
+if (is_asciifile) then
! Read the character identifier for verbose formatted output
- read(ifile, FMT='(a8)') header
+ read(ifile, FMT="(a8)") header
if(header /= 'platform') then
- call error_handler(E_ERR,'obs_def_radar_mod:read_rad_vel', &
- 'Expected location header "platform" in input file', &
+ call error_handler(E_ERR,'read_radial_vel', &
+ "Expected location header 'platform' in input file", &
source, revision, revdate)
endif
- location = read_location(ifile, fileformat)
- orientation = read_orientation(ifile, fileformat)
- nyquistv = read_nyquistvel(ifile, fileformat)
- ! Read in the velkey for this particular observation
- read(ifile, *) velkey
-END SELECT
+endif
-velkeycount = velkeycount + 1 ! the total velocity metadata key count from all sequences
-velkey = velkeycount ! copied to the output variable
+! read_location is a DART library routine that expects an optional string
+! arg, but the other two read routines are local to this module and we can
+! tell them exactly what format to be reading because we already know it.
+radar_location = read_location (ifile, fform)
+beam_direction = read_beam_direction (ifile, is_asciifile)
+nyquist_velocity = read_nyquist_velocity(ifile, is_asciifile)
-call set_rad_vel(velkey, location, orientation, nyquistv)
+! Read in the velkey for this particular observation, however, it will
+! be discarded and a new, unique key will be generated in the set routine.
+if (is_asciifile) then
+ read(ifile, *) oldkey
+else
+ read(ifile) oldkey
+endif
-end subroutine read_rad_vel
+! Generate new unique radial velocity observation key, and set the contents
+! of the private defined type.
+call set_radial_vel(velkey, radar_location, beam_direction, nyquist_velocity)
+end subroutine read_radial_vel
+
!----------------------------------------------------------------------
-subroutine set_rad_vel(velkey, rad_location, rad_orientation, rad_nyquistv)
+subroutine write_radial_vel(velkey, ifile, fform)
-integer, intent(in) :: velkey
-real(r8), intent(in) :: rad_orientation(3)
-type(location_type), intent(in) :: rad_location
-real(r8), intent(in) :: rad_nyquistv
+! Write radial velocity auxiliary information to the obs_seq file.
+integer, intent(in) :: velkey, ifile
+character(len=*), intent(in), optional :: fform
+
+logical :: is_asciifile
+type(location_type) :: radar_location
+real(r8) :: beam_direction(3)
+real(r8) :: nyquist_velocity
+
if ( .not. module_initialized ) call initialize_module
-!Make sure enough space is allocated
-if(velkey > max_rad_vel_obs) then
- write(*, *) 'velkey (',velkey,') exceeds max_rad_vel_obs (',max_rad_vel_obs,')'
- call error_handler(E_ERR,'read_rad_vel:set_rad_vel', &
- 'Increase max_rad_vel_obs.', source, revision, revdate)
+! Simple error check on key number before accessing the array
+call velkey_out_of_range(velkey)
+
+is_asciifile = ascii_file_format(fform)
+
+if (is_asciifile) then
+ ! Write the 5 character identifier for verbose formatted output
+ write(ifile, "('platform')")
endif
-rad_loc(velkey) = rad_location
-direction(:, velkey) = rad_orientation
-nyquistvel(velkey) = rad_nyquistv
+! Extract the values for this key and call the appropriate write routines.
+radar_location = radial_vel_data(velkey)%radar_location
+beam_direction(:) = radial_vel_data(velkey)%beam_direction(:)
+nyquist_velocity = radial_vel_data(velkey)%nyquist_velocity
-end subroutine set_rad_vel
+! write_location routine is part of the DART library and wants the optional
+! format string argument. The other two routines are local to this module,
+! and we have already figured out if it is a unformatted/binary file or
+! formatted/ascii, so go ahead and pass that info directly down to the routines.
+call write_location(ifile, radar_location, fform)
+call write_beam_direction(ifile, beam_direction(:), is_asciifile)
+call write_nyquist_velocity(ifile, nyquist_velocity, is_asciifile)
+! Write out the velkey used for this run, however this will be discarded
+! when this observation is read in and a new key will be generated.
+! (It may be useful for correlating error messages or identifying particular
+! observations so we are leaving it as part of the aux data.)
+if (is_asciifile) then
+ write(ifile, *) velkey
+else
+ write(ifile) velkey
+endif
+
+end subroutine write_radial_vel
+
!----------------------------------------------------------------------
-subroutine write_rad_ref(refkey, ifile, fform)
+function read_beam_direction(ifile, is_asciiformat)
-integer, intent(in) :: refkey, ifile
-character(len=*), intent(in), optional :: fform
+! Reads beam_direction from file that was written by write_beam_direction.
+! See read_radial_vel for additional discussion.
-character(len=32) :: fileformat
+integer, intent(in) :: ifile
+logical, intent(in) :: is_asciiformat
+real(r8) :: read_beam_direction(3)
+character(len=5) :: header
+real(r8) :: beam_direction(3)
+
if ( .not. module_initialized ) call initialize_module
-fileformat = "ascii" ! supply default
-if(present(fform)) fileformat = trim(adjustl(fform))
-! At this point, this is empty as there is no metadata for radar reflectivity
-! in the current obs sequence format. In the future, if more metadata pieces
-! are added, this part should be modified.
+if (is_asciiformat) then
+ read(ifile, "(a5)" ) header
-end subroutine write_rad_ref
+ if(header /= 'dir3d') then
+ write(msgstring,*)"Expected beam_direction header 'dir3d' in input file, got ", header
+ call error_handler(E_ERR, 'read_beam_direction', msgstring, source, revision, revdate)
+ endif
+ ! Now read the beam_direction data value into temporaries
+ read(ifile, *) beam_direction(1), beam_direction(2), beam_direction(3)
+else
+ ! No header label, just the binary direction values.
+ read(ifile) beam_direction(1), beam_direction(2), beam_direction(3)
+endif
-!----------------------------------------------------------------------
+! Check for illegal values
+if (minval(beam_direction) < -1.0_r8 .or. maxval(beam_direction) > 1.0_r8) then
+ write(msgstring,*) "beam_direction value must be between -1 and 1, got: ", &
+ beam_direction(1), beam_direction(2), beam_direction(3)
+ call error_handler(E_ERR, 'read_beam_direction', msgstring, &
+ source, revision, revdate)
+endif
-subroutine get_obs_def_rad_ref(refkey, is_this_dbz)
+! set function return value
+read_beam_direction(:) = beam_direction(:)
-integer, intent(in) :: refkey
-logical, intent(out) :: is_this_dbz
-
-is_this_dbz = convert_to_dbz
-
-end subroutine get_obs_def_rad_ref
+end function read_beam_direction
!----------------------------------------------------------------------
-subroutine read_rad_ref(obsvalue, refkey, ifile, fform)
+subroutine write_beam_direction(ifile, beam_direction, is_asciiformat)
-! Main read subroutine for the obs kind radar reflectivity metadata.
+! Writes beam_direction to obs file.
-! reftype: Denotes whether reflectivity obs units are in Z (reftype = .false.)
-! or are in dBZ (reftype = .true.)
+integer, intent(in) :: ifile
+real(r8), intent(in) :: beam_direction(3)
+logical, intent(in) :: is_asciiformat
-real(r8), intent(inout) :: obsvalue
-integer, intent(out) :: refkey
-integer, intent(in) :: ifile
-character(len=*), intent(in), optional :: fform
+if ( .not. module_initialized ) call initialize_module
-character(len=32) :: fileformat
-!logical :: reftype
+if (is_asciiformat) then
+ write(ifile, "('dir3d')" )
+ write(ifile, *) beam_direction(1), beam_direction(2), beam_direction(3)
+else
+ write(ifile) beam_direction(1), beam_direction(2), beam_direction(3)
+endif
-if ( .not. module_initialized ) call initialize_module
+end subroutine write_beam_direction
-fileformat = "ascii" ! supply default
-if(present(fform)) fileformat = trim(adjustl(fform))
+!----------------------------------------------------------------------
-! At this point, no metadata is available for radar reflectivity. For this
-! reason, internal counter "refkey" has also no consequence. In the
-! future, if more metadata pieces are added, this part should be modified.
-refkeycount = refkeycount + 1 ! the total reflectivity metadata key count from all sequences
-refkey = refkeycount ! copied to the output variable
+function read_nyquist_velocity(ifile, is_asciiformat)
-! Convert to Z and appy reflectivity limit if requested through the namelist
-! (Note that conversion to Z is suppressed for now by globally resetting
-! convert_to_dbz to .TRUE. after reading namelist)
-if (obsvalue /= missing_r8) then
- if (.not. convert_to_dbz) then
- obsvalue = 10.0_r8 ** (obsvalue/10.0_r8)
- endif
+! Reads Nyquist velocity from file that was written by write_nyquist_velocity.
+! See read_radial_vel for additional discussion.
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list