[Dart-dev] [3476] DART/trunk/obs_def/obs_def_gps_mod.f90:
Updated gps observation module.
nancy at ucar.edu
nancy at ucar.edu
Mon Jul 28 11:50:03 MDT 2008
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20080728/175a6480/attachment-0001.html
-------------- next part --------------
Modified: DART/trunk/obs_def/obs_def_gps_mod.f90
===================================================================
--- DART/trunk/obs_def/obs_def_gps_mod.f90 2008-07-25 21:20:18 UTC (rev 3475)
+++ DART/trunk/obs_def/obs_def_gps_mod.f90 2008-07-28 17:50:02 UTC (rev 3476)
@@ -3,6 +3,12 @@
! University Corporation for Atmospheric Research
! Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+! Note: This version has a namelist item for the max number of
+! gps observations that can be read in, but it is currently commented out.
+! Search below for the 'NAMELIST' string to find where to comment the
+! code back in. There are 2 places. This maximum must be the total of
+! all gps obs in all input files, which if you are reading multiple obs_seq
+! files (e.g. for the obs_diag program) might be a larger number than 100K.
! <next few lines under version control, do not edit>
! $URL$
@@ -19,7 +25,8 @@
! BEGIN DART PREPROCESS USE OF SPECIAL OBS_DEF MODULE
-! use obs_def_gps_mod, only : get_expected_gpsro_ref, read_gpsro_ref, write_gpsro_ref
+! use obs_def_gps_mod, only : get_expected_gpsro_ref, interactive_gpsro_ref, &
+! read_gpsro_ref, write_gpsro_ref
! END DART PREPROCESS USE OF SPECIAL OBS_DEF MODULE
@@ -31,20 +38,19 @@
! BEGIN DART PREPROCESS READ_OBS_DEF
! case(GPSRO_REFRACTIVITY)
-! call read_gpsro_ref(obs_def%key, ifile, fileformat)
-! continue
+! call read_gpsro_ref(obs_def%key, ifile, fileformat)
! END DART PREPROCESS READ_OBS_DEF
! BEGIN DART PREPROCESS WRITE_OBS_DEF
! case(GPSRO_REFRACTIVITY)
-! call write_gpsro_ref(obs_def%key, ifile, fileformat)
+! call write_gpsro_ref(obs_def%key, ifile, fileformat)
! END DART PREPROCESS WRITE_OBS_DEF
! BEGIN DART PREPROCESS INTERACTIVE_OBS_DEF
! case(GPSRO_REFRACTIVITY)
-! ! call interactive_gpsro_ref(obs_def%key)
+! call interactive_gpsro_ref(obs_def%key)
! END DART PREPROCESS INTERACTIVE_OBS_DEF
@@ -52,23 +58,27 @@
module obs_def_gps_mod
use types_mod, only : r8, missing_r8, RAD2DEG, DEG2RAD, PI
-use utilities_mod, only : register_module, error_handler, E_ERR, E_MSG, file_exist, &
- open_file, close_file
-use location_mod, only : location_type, set_location, get_location , write_location, &
- read_location
-use time_manager_mod, only : time_type, read_time, write_time, set_time, set_time_missing, &
- interactive_time
+use utilities_mod, only : register_module, error_handler, E_ERR, E_MSG, &
+ file_exist, open_file, close_file, nmlfileunit, &
+ check_namelist_read, find_namelist_in_file, &
+ do_output
+use location_mod, only : location_type, set_location, get_location, &
+ write_location, read_location, vert_is_height, &
+ VERTISHEIGHT
+use time_manager_mod, only : time_type, read_time, write_time, &
+ set_time, set_time_missing, interactive_time
use assim_model_mod, only : interpolate
use obs_kind_mod, only : KIND_U_WIND_COMPONENT, &
KIND_V_WIND_COMPONENT, KIND_SURFACE_PRESSURE, &
- KIND_TEMPERATURE, KIND_SPECIFIC_HUMIDITY, KIND_PRESSURE, &
- KIND_GPSRO
+ KIND_TEMPERATURE, KIND_SPECIFIC_HUMIDITY, &
+ KIND_PRESSURE, KIND_GPSRO
implicit none
private
-public :: set_gpsro_ref, write_gpsro_ref, read_gpsro_ref, get_expected_gpsro_ref
+public :: set_gpsro_ref, get_gpsro_ref, write_gpsro_ref, read_gpsro_ref, &
+ get_expected_gpsro_ref, interactive_gpsro_ref
! version controlled file description for error handling, do not edit
character(len=128), parameter :: &
@@ -80,63 +90,148 @@
! Storage for the special information required for GPS RO observations
!
-integer, parameter :: max_gpsro_obs = 100000
-character(len=6) :: gpsro_ref_form(max_gpsro_obs)
-real(r8) :: ray_direction(3, max_gpsro_obs)
-real(r8) :: rfict(max_gpsro_obs) ! local curvature radius
-real(r8) :: step_size(max_gpsro_obs)
-real(r8) :: ray_top(max_gpsro_obs)
+!integer, parameter :: max_gpsro_obs = 100000
+!character(len=6) :: gpsro_ref_form(max_gpsro_obs)
+!real(r8) :: ray_direction(3, max_gpsro_obs)
+!real(r8) :: rfict(max_gpsro_obs) ! local curvature radius
+!real(r8) :: step_size(max_gpsro_obs)
+!real(r8) :: ray_top(max_gpsro_obs)
+
+! Because we are currently only generating one observation type
+! (GPSRO_REFRACTIVITY), there must be enough of these to cover all gps
+! obs in all obs_seq files that are read in (e.g. for obs_diag if you
+! cover multiple days or weeks, you must have enough room for all of them.)
+! the local operator needs none of this additional info; the best approach
+! would be to keep a single KIND_GPSRO, but make 2 observation types.
+! the local has no additional metadata; the nonlocal needs one of these
+! allocated and filled in.
+integer :: max_gpsro_obs = 100000
+type gps_nonlocal_type
+ private
+ character(len=6) :: gpsro_ref_form
+ real(r8) :: ray_direction(3)
+ real(r8) :: rfict
+ real(r8) :: step_size
+ real(r8) :: ray_top
+end type gps_nonlocal_type
+
+type(gps_nonlocal_type), allocatable :: gps_data(:)
+
+! NAMELIST - comment in this plus one more place below to enable run-time limit.
+!namelist /obs_def_gps_nml/ max_gpsro_obs
+
+character(len=129) :: msgstring
integer :: ii
-integer :: keycount = 0
+integer :: keycount
contains
-!---------------------------------------------------------------------------------
+!------------------------------------------------------------------------------
subroutine initialize_module
-!---------------------------------------------------------------------------------
-! subroutine initialize_module
+!------------------------------------------------------------------------------
!
-! pretty simple for this module.
+! initialize global gps private key number and allocate space for obs data
+integer :: rc, iunit
+
call register_module(source, revision, revdate)
module_initialized = .true.
+! global count of all gps observations from any input file
+keycount = 0
+
+! NAMELIST - comment in this block to enable runtime setting of the max obs.
+!! Read the namelist entry
+!! find max number of gps obs which can be stored, and initialize type
+!call find_namelist_in_file("input.nml", "obs_def_gps_nml", iunit)
+!read(iunit, nml = obs_def_gps_nml, iostat = rc)
+!call check_namelist_read(iunit, rc, "obs_def_gps_nml")
+!
+!! Record the namelist values used for the run ...
+!if (do_output()) write(nmlfileunit, nml=obs_def_gps_nml)
+!if (do_output()) write( * , nml=obs_def_gps_nml)
+
+allocate(gps_data(max_gpsro_obs), stat = rc)
+if (rc /= 0) then
+ write(msgstring, *) 'initial allocation failed for gps observation data,', &
+ 'itemcount = ', max_gpsro_obs
+ call error_handler(E_ERR,'initialize_module', msgstring, &
+ source, revision, revdate)
+endif
+
end subroutine initialize_module
- subroutine set_gpsro_ref(key, nx, ny, nz, rfict0, ds, htop, subset0)
-!---------------------------------------------------------------------------------
-!subroutine set_gpsro_ref(key, nx, ny, nz, rfict0, ds, htop, subset0)
+ subroutine set_gpsro_ref(gpskey, nx, ny, nz, rfict0, ds, htop, subset0)
+!------------------------------------------------------------------------------
!
+! increment key and set all private data for this observation
-integer, intent(in) :: key
-real(r8), intent(in) :: nx, ny, nz, rfict0, ds, htop
-character(len=6), intent(in) :: subset0
+integer, intent(out) :: gpskey
+real(r8), intent(in) :: nx, ny, nz, rfict0, ds, htop
+character(len=6), intent(in) :: subset0
if ( .not. module_initialized ) call initialize_module
-ray_direction(1, key) = nx
-ray_direction(2, key) = ny
-ray_direction(3, key) = nz
-gpsro_ref_form(key) = subset0
+keycount = keycount + 1
+gpskey = keycount
- rfict(key) = rfict0
-step_size(key) = ds
- ray_top(key) = htop
+if(gpskey > max_gpsro_obs) then
+ write(msgstring, *) 'key (',gpskey,') exceeds max_gpsro_obs (',max_gpsro_obs,')'
+ call error_handler(E_ERR,'read_gpsro_ref', msgstring, &
+ source, revision, revdate)
+endif
+gps_data(gpskey)%ray_direction(1) = nx
+gps_data(gpskey)%ray_direction(2) = ny
+gps_data(gpskey)%ray_direction(3) = nz
+gps_data(gpskey)%gpsro_ref_form = subset0
+
+gps_data(gpskey)%rfict = rfict0
+gps_data(gpskey)%step_size = ds
+gps_data(gpskey)%ray_top = htop
+
end subroutine set_gpsro_ref
+ subroutine get_gpsro_ref(gpskey, nx, ny, nz, rfict0, ds, htop, subset0)
+!------------------------------------------------------------------------------
+!
+! return all private data for this observation
- subroutine write_gpsro_ref(key, ifile, fform)
-!---------------------------------------------------------------------------------
-!subroutine write_gpsro_ref(key, ifile, fform)
+integer, intent(in) :: gpskey
+real(r8), intent(out) :: nx, ny, nz, rfict0, ds, htop
+character(len=6), intent(out) :: subset0
+
+if ( .not. module_initialized ) call initialize_module
+
+if (gpskey < 1 .or. gpskey > keycount) then
+ write(msgstring, *) 'key (',gpskey,') out of valid range (1<=key<=',keycount,')'
+ call error_handler(E_ERR,'get_gpsro_ref', msgstring, &
+ source, revision, revdate)
+endif
+
+nx = gps_data(gpskey)%ray_direction(1)
+ny = gps_data(gpskey)%ray_direction(2)
+nz = gps_data(gpskey)%ray_direction(3)
+subset0 = gps_data(gpskey)%gpsro_ref_form
+
+rfict0 = gps_data(gpskey)%rfict
+ds = gps_data(gpskey)%step_size
+htop = gps_data(gpskey)%ray_top
+
+end subroutine get_gpsro_ref
+
+
+
+ subroutine write_gpsro_ref(gpskey, ifile, fform)
+!------------------------------------------------------------------------------
!
-integer, intent(in) :: key, ifile
+integer, intent(in) :: gpskey, ifile
character(len=*), intent(in), optional :: fform
character(len=32) :: fileformat
@@ -151,16 +246,20 @@
SELECT CASE (fileformat)
CASE ("unf", "UNF", "unformatted", "UNFORMATTED")
- write(ifile) key
- write(ifile) rfict(key), step_size(key), ray_top(key), &
- (ray_direction(ii, key), ii=1, 3), gpsro_ref_form(key)
+ write(ifile) gpskey
+ write(ifile) gps_data(gpskey)%rfict, gps_data(gpskey)%step_size, &
+ gps_data(gpskey)%ray_top, &
+ (gps_data(gpskey)%ray_direction(ii), ii=1, 3), &
+ gps_data(gpskey)%gpsro_ref_form
continue
CASE DEFAULT
- write(ifile,11) key
- write(ifile, *) rfict(key), step_size(key), ray_top(key), &
- (ray_direction(ii, key), ii=1, 3), gpsro_ref_form(key)
+ write(ifile,11) gpskey
+ write(ifile, *) gps_data(gpskey)%rfict, gps_data(gpskey)%step_size, &
+ gps_data(gpskey)%ray_top, &
+ (gps_data(gpskey)%ray_direction(ii), ii=1, 3), &
+ gps_data(gpskey)%gpsro_ref_form
END SELECT
11 format('gpsroref', i8)
@@ -168,88 +267,164 @@
- subroutine read_gpsro_ref(key, ifile, fform)
-!---------------------------------------------------------------------------------
-!subroutine read_gpsro_ref(key, ifile, fform)
+ subroutine read_gpsro_ref(gpskey, ifile, fform)
+!------------------------------------------------------------------------------
!
-! Every GPS observation has its own (metadata) key.
-! When you read two gps observation sequence files, it is necessary to track the
-! total number of metadata keys read ... not just the number in the current file.
+! Every GPS observation has its own (metadata) gpskey.
+! When you read multiple gps observation sequence files, it is necessary
+! to track the total number of metadata gpskeys read, not just the number
+! in the current file.
!
-integer, intent(out) :: key
+integer, intent(out) :: gpskey
integer, intent(in) :: ifile
character(len=*), intent(in), optional :: fform
-integer :: keyin ! the metadata key in the current obs sequence
+integer :: keyin ! the metadata key in the current obs sequence
-character(len=8) :: header
-character(len=32) :: fileformat
+real(r8) :: nx, ny, nz, rfict0, ds, htop
+character(len=6) :: subset0
+character(len=8) :: header
+character(len=32) :: fileformat
if ( .not. module_initialized ) call initialize_module
fileformat = "ascii" ! supply default
if(present(fform)) fileformat = trim(adjustl(fform))
-keycount = keycount + 1 ! the total metadata key count from all sequences
-key = keycount ! copied to the output variable
-
-if(key > max_gpsro_obs) then
- write(*, *) 'key (',key,') exceeds max_gpsro_obs (',max_gpsro_obs,')'
- call error_handler(E_ERR,'read_gpsro_ref', &
- 'Increase max_gpsro_obs.', source, revision, revdate)
-endif
-
! Read the character identifier for verbose formatted output
SELECT CASE (fileformat)
CASE ("unf", "UNF", "unformatted", "UNFORMATTED")
- read(ifile) keyin ! read and throw away
- read(ifile) rfict(key), step_size(key), ray_top(key), &
- (ray_direction(ii, key), ii=1, 3), gpsro_ref_form(key)
- continue
+ read(ifile) keyin ! read and throw away
+ read(ifile) rfict0, ds, htop, nx, ny, nz, subset0
CASE DEFAULT
read(ifile, FMT='(a8, i8)') header, keyin ! throw away keyin
if(header /= 'gpsroref') then
- call error_handler(E_ERR,'read_gpsro_ref', &
- 'Expected header "gpsroref" in input file', source, revision, revdate)
+ call error_handler(E_ERR,'read_gpsro_ref', &
+ 'Expected header "gpsroref" in input file', source, revision, revdate)
endif
- read(ifile, *) rfict(key), step_size(key), ray_top(key), &
- (ray_direction(ii, key), ii=1, 3), gpsro_ref_form(key)
+ read(ifile, *) rfict0, ds, htop, nx, ny, nz, subset0
END SELECT
+
+! increment key and set all private data for this observation
+call set_gpsro_ref(gpskey, nx, ny, nz, rfict0, ds, htop, subset0)
+
end subroutine read_gpsro_ref
+subroutine interactive_gpsro_ref(gpskey)
+!----------------------------------------------------------------------
+!
+! Interactively prompt for the info needed to create a gps refractivity
+! observation. Increments the key number and returns it.
- subroutine get_expected_gpsro_ref(state_vector, location, key, ro_ref, istatus)
-!---------------------------------------------------------------------------------
-!subroutine get_expected_gpsro_ref(state_vector, location, key, ro_ref, istatus)
+integer, intent(out) :: gpskey
+
+real(r8) :: nx, ny, nz, rfict0, ds, htop
+character(len=6) :: subset0
+type(location_type) :: location
+integer :: gpstype
+
+
+if ( .not. module_initialized ) call initialize_module
+
+!Now interactively obtain reflectivity type information
+! valid choices are local or non-local
+
+write(*, *)
+write(*, *) 'Beginning to inquire information on reflectivity type.'
+write(*, *)
+
+100 continue
+write(*, *) 'Enter 1 for local refractivity (GPSREF)'
+write(*, *) 'Enter 2 for non-local refractivity/excess phase delay (GPSEXC)'
+write(*, *)
+
+read(*,*) gpstype
+
+select case (gpstype)
+ case (1)
+ subset0 = 'GPSREF'
+ case (2)
+ subset0 = 'GPSEXC'
+ case default
+ write(*,*) 'Bad value, must enter 1 or 2'
+ goto 100
+end select
+
+if (gpstype == 2) then
+ ! FIXME: i have no idea what valid values are for any
+ ! of the following items, so i cannot add any error checking or
+ ! guidance for the user.
+
+ write(*, *)
+ write(*, *) 'Enter X, Y, Z value for ray direction'
+ write(*, *)
+ read(*,*) nx, ny, nz
+
+ write(*, *)
+ write(*, *) 'Enter local curvature radius'
+ write(*, *)
+ read(*,*) rfict0
+
+ write(*, *)
+ write(*, *) 'Enter step size'
+ write(*, *)
+ read(*,*) ds
+
+ write(*, *)
+ write(*, *) 'Enter ray top'
+ write(*, *)
+ read(*,*) htop
+else
+ nx = 0.0
+ ny = 0.0
+ nz = 0.0
+ rfict0 = 0.0
+ ds = 0.0
+ htop = 0.0
+endif
+
+! increment key and set all private data for this observation
+call set_gpsro_ref(gpskey, nx, ny, nz, rfict0, ds, htop, subset0)
+
+write(*, *)
+write(*, *) 'End of specialized section for gps observation data.'
+write(*, *) 'You will now have to enter the regular obs information.'
+write(*, *)
+
+end subroutine interactive_gpsro_ref
+
+ subroutine get_expected_gpsro_ref(state_vector, location, gpskey, ro_ref, istatus)
+!------------------------------------------------------------------------------
!
! Purpose: Calculate GPS RO local refractivity or non_local (integrated)
! refractivity (excess phase, Sergey Sokolovskiy et al., 2005)
-!---------------------------------------------------------------------------------
+!------------------------------------------------------------------------------
!
! inputs:
! state_vector: DART state vector
!
! output parameters:
-! ro_ref: modeled local refractivity (N-1)*1.0e6 or non_local refractivity (excess phase, m)
+! ro_ref: modeled local refractivity (N-1)*1.0e6 or non_local
+! refractivity (excess phase, m)
! (according to the input data parameter subset)
! istatus: =0 normal; =1 outside of domain.
-!---------------------------------------------------------------------------------
+!------------------------------------------------------------------------------
! Author: Hui Liu
! Version 1.1: June 15, 2004: Initial version CAM
!
! Version 1.2: July 29, 2005: revised for new obs_def and WRF
-!---------------------------------------------------------------------------------
+!------------------------------------------------------------------------------
implicit none
real(r8), intent(in) :: state_vector(:)
-integer, intent(in) :: key
type(location_type), intent(in) :: location
+integer, intent(in) :: gpskey
real(r8), intent(out) :: ro_ref
integer, intent(out) :: istatus
@@ -267,8 +442,15 @@
if ( .not. module_initialized ) call initialize_module
+if ( .not. vert_is_height(location)) then
+ write(msgstring, *) 'vertical location must be height; gps obs key ', gpskey
+ call error_handler(E_ERR,'get_expected_gpsro_ref', msgstring, &
+ source, revision, revdate)
+endif
+
obsloc = get_location(location)
+
lon = obsloc(1) ! degree: 0 to 360
lat = obsloc(2) ! degree: -90 to 90
height = obsloc(3) ! (m)
@@ -276,34 +458,32 @@
! calculate refractivity at perigee
call ref_local(state_vector, location, height, lat, lon, ref_perigee, istatus0)
-istatus = istatus0
+! if istatus > 0, the interpolation failed and we should return failure now.
+if(istatus0 > 0) then
+ istatus = istatus0
+ ro_ref = missing_r8
+endif
-choose: if(gpsro_ref_form(key) == 'GPSREF') then ! use local refractivity
+choose: if(gps_data(gpskey)%gpsro_ref_form == 'GPSREF') then
+ ! use local refractivity
ro_ref = ref_perigee * 1.0e6 ! in (N-1)*1.0e6, same with obs
- if(istatus > 1) then
- istatus = 1
- ro_ref = missing_r8
- endif
+else ! gps_data(gpskey)%gpsro_ref_form == 'GPSEXC'
- return
-
-else ! gpsro_ref_form(key) == 'GPSEXC'
-
! otherwise, use non_local refractivity(excess phase delay)
- ! Initilization
+ ! Initialization
phase = 0.0_r8
dist_to_perigee = 0.0_r8 ! distance to perigee from a point of the ray
- nx = ray_direction(1, key)
- ny = ray_direction(2, key)
- nz = ray_direction(3, key)
+ nx = gps_data(gpskey)%ray_direction(1)
+ ny = gps_data(gpskey)%ray_direction(2)
+ nz = gps_data(gpskey)%ray_direction(3)
! convert location of the perigee from geodetic to Cartesian coordinate
- call geo2carte (height, lat, lon, xo, yo, zo, rfict(key) )
+ call geo2carte (height, lat, lon, xo, yo, zo, gps_data(gpskey)%rfict )
! currently, use a straight line passing the perigee point as ray model.
! later, more sophisticated ray models can be used.
@@ -322,7 +502,7 @@
100 continue
iter = iter + 1
- dist_to_perigee = dist_to_perigee + step_size(key)
+ dist_to_perigee = dist_to_perigee + gps_data(gpskey)%step_size
! integrate to one direction of the ray for one step
xx = xo + dist_to_perigee * nx
@@ -332,14 +512,19 @@
! convert the location of the point to geodetic coordinates
! height(m), lat, lon(deg)
- call carte2geo(xx, yy, zz, height1, lat1, lon1, rfict(key) )
+ call carte2geo(xx, yy, zz, height1, lat1, lon1, gps_data(gpskey)%rfict )
! get the refractivity at this ray point(ref00)
call ref_local(state_vector, location, height1, lat1, lon1, ref00, istatus0)
- istatus = istatus + istatus0
+ ! when any point of the ray is problematic, return failure
+ if(istatus0 > 0) then
+ istatus = istatus0
+ ro_ref = missing_r8
+ return
+ endif
! get the excess phase due to this ray interval
- delta_phase1 = (ref1 + ref00) * step_size(key) * 0.5_r8
+ delta_phase1 = (ref1 + ref00) * gps_data(gpskey)%step_size * 0.5_r8
! save the refractivity for integration of next ray interval
ref1 = ref00
@@ -349,14 +534,19 @@
yy = yo - dist_to_perigee * ny
zz = zo - dist_to_perigee * nz
- call carte2geo (xx, yy, zz, height1, lat1, lon1, rfict(key) )
+ call carte2geo (xx, yy, zz, height1, lat1, lon1, gps_data(gpskey)%rfict )
! get the refractivity at this ray point(ref00)
call ref_local(state_vector, location, height1, lat1, lon1, ref00, istatus0)
- istatus = istatus + istatus0
+ ! when any point of the ray is problematic, return failure
+ if(istatus0 > 0) then
+ istatus = istatus0
+ ro_ref = missing_r8
+ return
+ endif
! get the excess phase due to this ray interval
- delta_phase2 = (ref2 + ref00) * step_size(key) * 0.5_r8
+ delta_phase2 = (ref2 + ref00) * gps_data(gpskey)%step_size * 0.5_r8
! save the refractivity for integration of next ray interval
ref2 = ref00
@@ -364,39 +554,35 @@
phase = phase + delta_phase1 + delta_phase2
! print*, 'phase= ', phase, delta_phase1, delta_phase2
- if(height1 .lt. ray_top(key) ) go to 100 !! do one more step
+ if(height1 .lt. gps_data(gpskey)%ray_top ) go to 100 !! do one more step
! finish the integration of the excess phase along the ray
ro_ref = phase ! in m
- ! when any point of the ray is problematic
- if(istatus > 1) then
- istatus = 1
- ro_ref = missing_r8
- endif
-
! print*, 'xx = ', lon, lat, height, ro_ref
endif choose
+! ended ok, return local or non-local accumulated value
+istatus = 0
+
end subroutine get_expected_gpsro_ref
subroutine ref_local(state_vector, location, height, lat, lon, ref00, istatus0)
-!---------------------------------------------------------------------------------
-!subroutine ref_local(state_vector, location, height, lat, lon, ref00, istatus0)
+!------------------------------------------------------------------------------
!
-! Calculate local refractivity at any GPS ray point (height, lat, lon)
+! Calculate local refractivity at any GPS ray point (height, lat, lon)
!
-! inputs:
-! height, lat, lon: GPS observation location (units: m, degree)
+! inputs:
+! height, lat, lon: GPS observation location (units: m, degree)
!
-! output:
-! ref00: modeled local refractivity at ray point(unit: N-1, ~1.0e-4 to e-6)
+! output:
+! ref00: modeled local refractivity at ray point(unit: N-1, ~1.0e-4 to e-6)
!
-!---------------------------------------------------------------------------------
+!------------------------------------------------------------------------------
implicit none
real(r8), intent(in) :: state_vector(:)
@@ -418,13 +604,20 @@
if(lon > 360.0_r8 ) lon2 = lon - 360.0_r8
if(lon < 0.0_r8 ) lon2 = lon + 360.0_r8
-which_vert = 3
+which_vert = VERTISHEIGHT
location2 = set_location(lon2, lat, height, which_vert)
-istatus0 = 0
+! set return values assuming failure, so we can simply return if any
+! of the interpolation calls below fail.
+istatus0 = 1
+ref00 = missing_r8
+
call interpolate(state_vector, location2, KIND_TEMPERATURE, t, istatus0)
+if (istatus0 > 0) return
call interpolate(state_vector, location2, KIND_SPECIFIC_HUMIDITY, q, istatus0)
+if (istatus0 > 0) return
call interpolate(state_vector, location2, KIND_PRESSURE, p, istatus0)
+if (istatus0 > 0) return
! required variable units for calculation of GPS refractivity
! t : Kelvin, from top to bottom
@@ -437,22 +630,24 @@
ew = q * p/(rdorv + (1.0_r8-rdorv)*q )
ref00 = c1*p/t + c2*ew/(t**2) ! (N-1)
+! now we have succeeded, set istatus to good
+istatus0 = 0
+
end subroutine ref_local
subroutine geo2carte (s1, s2, s3, x1, x2, x3, rfict0)
-!---------------------------------------------------------------------------------
-!subroutine geo2carte (s1, s2, s3, x1, x2, x3, rfict0)
+!------------------------------------------------------------------------------
!
! Converts geodetical coordinates to cartesian with a reference sphere
-!---------------------------------------------------------------------------------
+!------------------------------------------------------------------------------
! input parameters:
! s - geodetical coordinates
! (height (m), latitude (degree), longitude (degree))
! -90 to 90 0 to 360
! output parameters:
! x - cartesian coordinates (m) connected with the earth(x, y, z-coordinate)
-!---------------------------------------------------------------------------------
+!------------------------------------------------------------------------------
implicit none
real(r8), intent(in) :: s1, s2, s3, rfict0 ! units: m
real(r8), intent(out) :: x1, x2 ,x3
@@ -470,8 +665,7 @@
subroutine carte2geo (x1, x2, x3, s1, s2, s3, rfict0)
-!---------------------------------------------------------------------------------
-!subroutine carte2geo (x1, x2, x3, s1, s2, s3, rfict0)
+!------------------------------------------------------------------------------
!
! Converts cartesian coordinates to geodetical.
!
@@ -482,7 +676,7 @@
! s - geodetical coordinates
! (height (m), latitude (deg), longitude (deg))
! -90 to 90 0 to 360
-!---------------------------------------------------------------------------------
+!------------------------------------------------------------------------------
implicit none
real(r8), intent(in) :: x1, x2, x3, rfict0
real(r8), intent(out) :: s1, s2, s3
@@ -507,4 +701,5 @@
end subroutine carte2geo
end module obs_def_gps_mod
+
! END DART PREPROCESS MODULE CODE
More information about the Dart-dev
mailing list