[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