[Dart-dev] [4139] DART/trunk/obs_def: Updates from Glen to handle different types of radar
nancy at ucar.edu
nancy at ucar.edu
Fri Nov 6 14:02:21 MST 2009
Revision: 4139
Author: nancy
Date: 2009-11-06 14:02:21 -0700 (Fri, 06 Nov 2009)
Log Message:
-----------
Updates from Glen to handle different types of radar
physics schemes, and includes new support for obs of
fall velocity. Includes a couple new namelist items
and documentation. There may be additional changes
needed in this code to support more complex schemes,
but this should be a working baseline before it gets
more complicated.
Modified Paths:
--------------
DART/trunk/obs_def/obs_def_radar_mod.f90
DART/trunk/obs_def/obs_def_radar_mod.html
DART/trunk/obs_def/obs_def_radar_mod.nml
-------------- next part --------------
Modified: DART/trunk/obs_def/obs_def_radar_mod.f90
===================================================================
--- DART/trunk/obs_def/obs_def_radar_mod.f90 2009-11-06 20:13:09 UTC (rev 4138)
+++ DART/trunk/obs_def/obs_def_radar_mod.f90 2009-11-06 21:02:21 UTC (rev 4139)
@@ -39,7 +39,8 @@
! BEGIN DART PREPROCESS USE OF SPECIAL OBS_DEF MODULE
! 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
+! read_radar_ref, get_expected_radar_ref, &
+! get_expected_fall_velocity
! END DART PREPROCESS USE OF SPECIAL OBS_DEF MODULE
!-----------------------------------------------------------------------------
@@ -49,6 +50,8 @@
! 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)
+! case(PRECIPITATION_FALL_SPEED)
+! call get_expected_fall_velocity(state, location, obs_val, istatus)
! END DART PREPROCESS GET_EXPECTED_OBS_FROM_DEF
!-----------------------------------------------------------------------------
@@ -58,6 +61,8 @@
! call read_radial_vel(obs_def%key, ifile, fileformat)
! case(RADAR_REFLECTIVITY)
! call read_radar_ref(obs_val, obs_def%key)
+! case(PRECIPITATION_FALL_SPEED)
+! continue
! END DART PREPROCESS READ_OBS_DEF
!-----------------------------------------------------------------------------
@@ -67,6 +72,8 @@
! call write_radial_vel(obs_def%key, ifile, fileformat)
! case(RADAR_REFLECTIVITY)
! continue
+! case(PRECIPITATION_FALL_SPEED)
+! continue
! END DART PREPROCESS WRITE_OBS_DEF
!-----------------------------------------------------------------------------
@@ -76,6 +83,8 @@
! call interactive_radial_vel(obs_def%key)
! case(RADAR_REFLECTIVITY)
! continue
+! case(PRECIPITATION_FALL_SPEED)
+! continue
! END DART PREPROCESS INTERACTIVE_OBS_DEF
!-----------------------------------------------------------------------------
@@ -103,7 +112,8 @@
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
+ get_expected_radial_vel, get_obs_def_radial_vel, set_radial_vel, &
+ get_expected_fall_velocity
! version controlled file description for error handling, do not edit
character(len=128), parameter :: &
@@ -274,7 +284,32 @@
! 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.
-
+!
+! microphysics_type:
+! Integer. Tells obs_def_radar_mod what microphysical scheme is employed
+! to enable some smarter decisions in handling fall velocity and radar
+! reflectivity.
+! 1 - Kessler scheme
+! 2 - Lin et al. microphysics (default)
+! 3 - User selected scheme where 10 cm reflectivity and power weighted fall
+! velocity are expected in the state vector (failure if not found)
+! 4 - User selected scheme where only power weighted fall velocity is
+! expected (failure if not found)
+! 5 - User selected scheme where only reflectivity is expected (failure
+! if not found)
+! -1 - ASSUME FALL VELOCITY IS ZERO, allows over-riding the failure modes
+! above if reflectivity and/or fall velocity are not available but a
+! result is desired for testing purposes only.
+!
+! allow_dbztowt_conv:
+! Logical. Flag to enable use of the dbztowt routine where reflectivity is
+! available, but not the power-weighted fall velocity. This scheme uses
+! emperical relations between reflectivity and fall velocity, with poor
+! accuracy for highly reflective, low density particles (such as water coated
+! snow aggregates). Expect questionable accuracy in radial velocity from
+! the forward operator with high elevation angles where ice is present in
+! the model state.
+
logical :: apply_ref_limit_to_obs = .false.
real(r8) :: reflectivity_limit_obs = -10.0_r8
real(r8) :: lowest_reflectivity_obs = -10.0_r8
@@ -283,6 +318,8 @@
real(r8) :: lowest_reflectivity_fwd_op = -10.0_r8
integer :: max_radial_vel_obs = 1000000
logical :: allow_wet_graupel = .false.
+integer :: microphysics_type = 2
+logical :: allow_dbztowt_conv = .false.
! Constants which depend on the microphysics scheme used. Should be set
@@ -306,6 +343,8 @@
max_radial_vel_obs, &
allow_wet_graupel, &
dielectric_factor, &
+ microphysics_type, &
+ allow_dbztowt_conv, &
n0_rain, &
n0_graupel, &
n0_snow, &
@@ -822,40 +861,10 @@
return
endif
-! If the model can return the precipitation fall speed directly, let it do
-! so. Otherwise, compute the various fields individually and then do
-! the computation here. Note that the computation here is accurate only
-! for the simple single-moment microphysics schemes (e.g., Kessler or Lin).
-
-call interpolate(state_vector, location, KIND_POWER_WEIGHTED_FALL_SPEED, &
- precip_fall_speed, istatus)
+call get_expected_fall_velocity(state_vector, location, precip_fall_speed, istatus)
if (istatus /= 0) then
- call interpolate(state_vector, location, KIND_RAINWATER_MIXING_RATIO, qr, istatus)
- if (istatus /= 0) then
- radial_vel = missing_r8
- return
- endif
- call interpolate(state_vector, location, KIND_GRAUPEL_MIXING_RATIO, qg, istatus)
- if (istatus /= 0) then
- radial_vel = missing_r8
- return
- endif
- call interpolate(state_vector, location, KIND_SNOW_MIXING_RATIO, qs, istatus)
- if (istatus /= 0) then
- radial_vel = missing_r8
- return
- endif
- call interpolate(state_vector, location, KIND_DENSITY, rho, istatus)
- if (istatus /= 0) then
- radial_vel = missing_r8
- return
- endif
- call interpolate(state_vector, location, KIND_TEMPERATURE, temp, istatus)
- if (istatus /= 0) then
- radial_vel = missing_r8
- return
- endif
- call get_precip_fall_speed(qr, qg, qs, rho, temp, precip_fall_speed)
+ radial_vel = missing_r8
+ return
endif
radial_vel = radial_vel_data(velkey)%beam_direction(1) * u + &
@@ -883,7 +892,115 @@
endif
end subroutine get_expected_radial_vel
+
+!----------------------------------------------------------------------
+!----------------------------------------------------------------------
+! Expected fall velocity section
+!----------------------------------------------------------------------
+!----------------------------------------------------------------------
+subroutine get_expected_fall_velocity(state_vector, location, &
+ precip_fall_speed, istatus)
+
+! This is the main forward operator routine for the expected
+! fall velocity, and it also used as part of computing expected
+! radial velocity.
+
+real(r8), intent(in) :: state_vector(:)
+type(location_type), intent(in) :: location
+real(r8), intent(out) :: precip_fall_speed
+integer, intent(out) :: istatus
+
+! Local vars
+real(r8) :: qr, qg, qs, rho, temp, refl
+logical, save :: first_time = .true.
+
+! If the model can return the precipitation fall speed directly, let it do
+! so. Otherwise, attempt to compute if Kessler or Lin type microphysics,
+! or see if the dbztowt routine is desired. Note that the computation for
+! Lin and Kessler is expected to be reasonably accurate for 10cm radar.
+
+! Should we check microphysics_type var or just go ahead and try to get a value?
+istatus = 0
+precip_fall_speed = 0.0_r8
+call interpolate(state_vector, location, KIND_POWER_WEIGHTED_FALL_SPEED, &
+ precip_fall_speed, istatus)
+
+! If able to get value, return here.
+if (istatus == 0) return
+
+! If the user explicitly wanted to interpolate in the field, try to complain
+! if it could not. Note that the interp could fail for other reasons.
+if (microphysics_type == 3 .or. microphysics_type == 4) then
+ ! Could return a specific istatus code here to indicate this condition.
+ if (first_time) then
+ call error_handler(E_MSG,'get_expected_fall_velocity', &
+ 'interpolate failed. Fall Speed may NOT be in state vector', '', '', '')
+ first_time = .false.
+ endif
+ return
+endif
+
+! If not in the state vector, try to calculate it here based on the
+! setting of the microphysics_type namelist.
+
+! if Kessler or Lin we can compute the fall velocity
+if (microphysics_type == 1 .or. microphysics_type == 2) then
+ call interpolate(state_vector, location, KIND_RAINWATER_MIXING_RATIO, qr, istatus)
+ if (istatus /= 0) then
+ precip_fall_speed = missing_r8
+ return
+ endif
+ if (microphysics_type == 2) then
+ call interpolate(state_vector, location, KIND_GRAUPEL_MIXING_RATIO, qg, istatus)
+ if (istatus /= 0) then
+ precip_fall_speed = missing_r8
+ return
+ endif
+ call interpolate(state_vector, location, KIND_SNOW_MIXING_RATIO, qs, istatus)
+ if (istatus /= 0) then
+ precip_fall_speed = missing_r8
+ return
+ endif
+ endif
+ ! Done with Lin et al specific calls
+ call interpolate(state_vector, location, KIND_DENSITY, rho, istatus)
+ if (istatus /= 0) then
+ precip_fall_speed = missing_r8
+ return
+ endif
+ call interpolate(state_vector, location, KIND_TEMPERATURE, temp, istatus)
+ if (istatus /= 0) then
+ precip_fall_speed = missing_r8
+ return
+ endif
+ call get_LK_precip_fall_speed(qr, qg, qs, rho, temp, precip_fall_speed)
+ ! Done with Lin et al or Kessler -
+else if (microphysics_type == 5 .and. allow_dbztowt_conv) then
+ ! Provided reflectivity field - will estimate fall velocity using empirical relations
+ call get_expected_radar_ref(state_vector, location, refl, istatus)
+ if (istatus /= 0) then
+ precip_fall_speed = missing_r8
+ return
+ endif
+ call interpolate(state_vector, location, KIND_DENSITY, rho, istatus)
+ if (istatus /= 0) then
+ precip_fall_speed = missing_r8
+ return
+ endif
+ precip_fall_speed = dbztowt(refl, rho, missing_r8)
+else if (microphysics_type < 0) then
+ ! User requested setting fall velocity to zero - use with caution
+ precip_fall_speed = 0.0_r8
+else
+ ! Couldn't manage to compute a fall velocity
+ precip_fall_speed = missing_r8
+ istatus = 2
+endif
+
+
+end subroutine get_expected_fall_velocity
+
!----------------------------------------------------------------------
!----------------------------------------------------------------------
! Radar reflectivity
@@ -933,12 +1050,13 @@
! "interpolate()" ultimately calls model_mod routine "model_interpolate()"
! to get model values of qr, qg, qs, rho, and temp at the ob location. Then
-! the routine "get_reflectivity()" is called to compute the radar reflectivity
+! the routine "get_LK_reflectivity()" is called to compute the radar reflectivity
! factor value, Z, that corresponds to the hydrometeor and thermodynamic values.
-real(r8) :: qr, qg, qs, rho, temp
-real(r8) :: debug_location(3)
-logical :: debug = .false. ! set to .true. to enable debug printout
+real(r8) :: qr, qg, qs, rho, temp
+real(r8) :: debug_location(3)
+logical :: debug = .false. ! set to .true. to enable debug printout
+logical, save :: first_time = .true.
if ( .not. module_initialized ) call initialize_module
@@ -954,51 +1072,72 @@
! the computation here. Note that the computation here is accurate only
! for the simple single-moment microphysics schemes (e.g., Kessler or Lin).
+! Try to draw from state vector first
call interpolate(state_vector, location, KIND_RADAR_REFLECTIVITY, ref, istatus)
if (istatus /= 0) then
- call interpolate(state_vector, location, KIND_RAINWATER_MIXING_RATIO, &
- qr, istatus)
- if (istatus /= 0) then
- ref = missing_r8
+ ! If the user explicitly wanted to interpolate in the field, try to complain
+ ! if it could not. Note that the interp could fail for other reasons.
+ if (microphysics_type == 3 .or. microphysics_type == 5) then
+ ! Could return a specific istatus code here to indicate this condition.
+ if (first_time) then
+ call error_handler(E_MSG,'get_expected_radar_ref', &
+ 'interpolate failed. Reflectivity may NOT be in state vector', '', '', '')
+ first_time = .false.
+ endif
return
endif
+
+ if (microphysics_type == 1 .or. microphysics_type == 2) then
+
+ call interpolate(state_vector, location, KIND_RAINWATER_MIXING_RATIO, &
+ qr, istatus)
+ if (istatus /= 0) then
+ ref = missing_r8
+ return
+ endif
+
+ if (microphysics_type == 2) then
+ ! Also need some ice vars
+ call interpolate(state_vector, location, KIND_GRAUPEL_MIXING_RATIO, &
+ qg, istatus)
+ if (istatus /= 0) then
+ ref = missing_r8
+ return
+ endif
+
+ call interpolate(state_vector, location, KIND_SNOW_MIXING_RATIO, &
+ qs, istatus)
+ if (istatus /= 0) then
+ ref = missing_r8
+ return
+ endif
+ endif
+ call interpolate(state_vector, location, KIND_DENSITY, rho, istatus)
+ if (istatus /= 0) then
+ ref = missing_r8
+ return
+ endif
- call interpolate(state_vector, location, KIND_GRAUPEL_MIXING_RATIO, &
- qg, istatus)
- if (istatus /= 0) then
- ref = missing_r8
- return
- endif
+ call interpolate(state_vector, location, KIND_TEMPERATURE, temp, istatus)
+ if (istatus /= 0) then
+ ref = missing_r8
+ return
+ endif
- call interpolate(state_vector, location, KIND_SNOW_MIXING_RATIO, &
- qs, istatus)
- if (istatus /= 0) then
+ call get_LK_reflectivity(qr, qg, qs, rho, temp, ref)
+
+ ! Always convert to dbz. Make sure the value, before taking the logarithm,
+ ! is always slightly positive.
+ ! tiny() is a fortran intrinsic function that is > 0 by a very small amount.
+
+ ref = 10.0_r8 * log10(max(tiny(ref), ref))
+ else
+ ! not in state vector and not Lin et al or Kessler so can't do reflectivity
ref = missing_r8
- return
- endif
-
- call interpolate(state_vector, location, KIND_DENSITY, rho, istatus)
- if (istatus /= 0) then
- ref = missing_r8
- return
- endif
-
- call interpolate(state_vector, location, KIND_TEMPERATURE, temp, istatus)
- if (istatus /= 0) then
- ref = missing_r8
- return
- endif
-
- call get_reflectivity(qr, qg, qs, rho, temp, ref)
-
+ endif
endif
-! Always convert to dbz. Make sure the value, before taking the logarithm,
-! is always slightly positive.
-! tiny() is a fortran intrinsic function that is > 0 by a very small amount.
-ref = 10.0_r8 * log10(max(tiny(ref), ref))
-
if ((apply_ref_limit_to_fwd_op) .and. &
(ref < reflectivity_limit_fwd_op) .and. (ref /= missing_r8)) then
ref = lowest_reflectivity_fwd_op
@@ -1033,7 +1172,7 @@
!----------------------------------------------------------------------
!----------------------------------------------------------------------
-subroutine get_reflectivity(qr, qg, qs, rho, temp, ref)
+subroutine get_LK_reflectivity(qr, qg, qs, rho, temp, ref)
! Computes the equivalent radar reflectivity factor in mm^6 m^-3 for
! simple single-moment microphysics schemes such as the Kessler and
@@ -1090,11 +1229,11 @@
endif
endif
-end subroutine get_reflectivity
+end subroutine get_LK_reflectivity
!----------------------------------------------------------------------
-subroutine get_precip_fall_speed(qr, qg, qs, rho, temp, precip_fall_speed)
+subroutine get_LK_precip_fall_speed(qr, qg, qs, rho, temp, precip_fall_speed)
! Computes power-weighted precipitation fall speed in m s^-1 for simple single-
! moment microphysics schemes such as the Kessler and Lin et al. schemes.
@@ -1154,7 +1293,7 @@
endif
if (precip_fall_speed > 0.0_r8) then
- call get_reflectivity(qr, qg, qs, rho, temp, ref)
+ call get_LK_reflectivity(qr, qg, qs, rho, temp, ref)
if (ref > 0.0_r8) then
precip_fall_speed = precip_fall_speed/ref
else
@@ -1162,10 +1301,36 @@
endif
endif
-end subroutine get_precip_fall_speed
+end subroutine get_LK_precip_fall_speed
!----------------------------------------------------------------------
+function dbztowt(rf, rho, spval)
+
+! Convert reflectivity (in DBZ, not Z) to terminal fall speed.
+! (Code from the pyncommas system - author D. Dowell)
+
+real(r8), intent(in) :: rf ! reflectivity (dBZ)
+real(r8), intent(in) :: rho ! density (km/m**3)
+real(r8), intent(in) :: spval ! bad/missing data flag
+real(r8) :: dbztowt
+
+! Local vars
+real(r8) :: refl ! reflectivity (Z)
+
+if ( (rf == spval) .or. (rho == spval) ) then
+ dbztowt = spval
+else
+ ! Convert back to Z for this calculation.
+ refl = 10.0**(0.1*rf)
+ ! Original code used opposite sign - be careful if updating.
+ dbztowt = 2.6 * refl**0.107 * (1.2/rho)**0.4
+endif
+
+end function dbztowt
+
+!----------------------------------------------------------------------
+
subroutine initialize_constants()
! Initialize module global constants.
@@ -1256,7 +1421,7 @@
! In the expressions for reflectivity for each hydrometeor category, the
! following parameters are computed from the constants that do not vary in time
! and space. Computing these parameters now means that the equations in the
-! get_reflectivity subroutine will have the following simple form:
+! get_LK_reflectivity subroutine will have the following simple form:
!
! ref = param_refl_r * (rho*qr)**1.75.
!
Modified: DART/trunk/obs_def/obs_def_radar_mod.html
===================================================================
--- DART/trunk/obs_def/obs_def_radar_mod.html 2009-11-06 20:13:09 UTC (rev 4138)
+++ DART/trunk/obs_def/obs_def_radar_mod.html 2009-11-06 21:02:21 UTC (rev 4139)
@@ -856,6 +856,8 @@
lowest_reflectivity_fwd_op, &
max_radial_vel_obs, &
allow_wet_graupel, &
+ microphysics_type &
+ allow_dbztowt_conv &
dielectric_factor, &
n0_rain, &
n0_graupel, &
@@ -883,6 +885,8 @@
lowest_reflectivity_fwd_op = -10.0_r8,
max_radial_vel_obs = 1000000,
allow_wet_graupel = .FALSE.,
+ microphysics_type = 2 ,
+ allow_dbztowt_conv = .FALSE.,
dielectric_factor = 0.224_r8,
n0_rain = 8.0e6_r8,
n0_graupel = 4.0e6_r8,
@@ -958,6 +962,48 @@
the freezing level. In the future, it might be better to provide the
option of having a transition layer.
</TD></TR>
+<TR><!--contents--><TD valign=top>microphysics_type</TD>
+ <!-- type --><TD valign=top>integer</TD>
+ <!--descript--><TD>
+If the state vector contains the reflectivity or the power
+weighted fall speed, interpolate directly from those regardless of
+the setting of this item. If the state vector does not
+contain the fields, this value should be set to be
+compatible with whatever microphysical scheme is being
+used by the model. If the model is using a different
+microphysical scheme but has compatible fields to the ones
+listed below, setting this value will select the scheme to use.
+<br>
+1 = Kessler scheme.
+<br>
+2 = Lin et al. microphysics (default)
+<br>
+3 = User selected scheme where 10 cm reflectivity and power weighted fall
+ velocity are expected in the state vector (failure if not found)
+<br>
+4 = User selected scheme where only power weighted fall velocity is
+ expected (failure if not found)
+<br>
+5 = User selected scheme where only reflectivity is expected (failure
+ if not found)
+<br>
+-1 = ASSUME FALL VELOCITY IS ZERO, allows over-riding the failure modes
+ above if reflectivity and/or fall velocity are not available but a
+ result is desired for testing purposes only.
+<br>
+Default value is 2.
+ </TD></TR>
+<TR><!--contents--><TD valign=top>allow_dbztowt_conv</TD>
+ <!-- type --><TD valign=top>logical</TD>
+ <!--descript--><TD>
+Flag to enable use of the dbztowt routine where reflectivity is
+available, but not the power-weighted fall velocity. This scheme uses
+emperical relations between reflectivity and fall velocity, with poor
+accuracy for highly reflective, low density particles (such as water coated
+snow aggregates). Expect questionable accuracy in radial velocity from
+the forward operator with high elevation angles where ice is present in
+the model state. Default value is .FALSE.
+ </TD></TR>
<TR><!--contents--><TD valign=top>dielectric_factor</TD>
<!-- type --><TD valign=top>real(r8)</TD>
<!--descript--><TD>
Modified: DART/trunk/obs_def/obs_def_radar_mod.nml
===================================================================
--- DART/trunk/obs_def/obs_def_radar_mod.nml 2009-11-06 20:13:09 UTC (rev 4138)
+++ DART/trunk/obs_def/obs_def_radar_mod.nml 2009-11-06 21:02:21 UTC (rev 4139)
@@ -13,6 +13,8 @@
lowest_reflectivity_fwd_op = -10.0_r8,
max_radial_vel_obs = 1000000,
allow_wet_graupel = .false.,
+ microphysics_type = 2 ,
+ allow_dbztowt_conv = .false.,
dielectric_factor = 0.224_r8,
n0_rain = 8.0e6_r8,
n0_graupel = 4.0e6_r8,
@@ -21,4 +23,3 @@
rho_graupel = 400.0_r8,
rho_snow = 100.0_r8
/
-
More information about the Dart-dev
mailing list