[Dart-dev] [4171] DART/trunk/models/wrf: Observation preprocessor which is WRF aware, from Ryan Torn.
nancy at ucar.edu
nancy at ucar.edu
Mon Nov 30 13:58:52 MST 2009
Revision: 4171
Author: nancy
Date: 2009-11-30 13:58:52 -0700 (Mon, 30 Nov 2009)
Log Message:
-----------
Observation preprocessor which is WRF aware, from Ryan Torn.
Will select obs only within the wrf domain, will superob,
will select only particular obs types based on the namelist.
Added Paths:
-----------
DART/trunk/models/wrf/WRF_DART_utilities/wrf_dart_obs_preprocess.f90
DART/trunk/models/wrf/WRF_DART_utilities/wrf_dart_obs_preprocess.nml
DART/trunk/models/wrf/work/mkmf_wrf_dart_obs_preprocess
DART/trunk/models/wrf/work/path_names_wrf_dart_obs_preprocess
-------------- next part --------------
Added: DART/trunk/models/wrf/WRF_DART_utilities/wrf_dart_obs_preprocess.f90
===================================================================
--- DART/trunk/models/wrf/WRF_DART_utilities/wrf_dart_obs_preprocess.f90 (rev 0)
+++ DART/trunk/models/wrf/WRF_DART_utilities/wrf_dart_obs_preprocess.f90 2009-11-30 20:58:52 UTC (rev 4171)
@@ -0,0 +1,1715 @@
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! wrf_dart_obs_preprocess - WRF-DART utility program that at a
+! minimum removes all observations outside
+! of a WRF domain and add observations from
+! supplimental obs sequences. The program
+! assumes all data is from one observation
+! time. In addition, this program allows
+! the user to do the following functions:
+!
+! - remove observations near the lateral boundaries
+! - increase observation error near lateral boundaries
+! - remove observations above certain pressure/height levels
+! - remove observations where the model and obs topography are large
+! - remove significant level rawinsonde data
+! - remove rawinsonde observations near TC core
+! - superob aircraft and satellite wind data
+!
+! created Oct. 2007 Ryan Torn, NCAR/MMM
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+program wrf_dart_obs_preprocess
+
+use types_mod, only : r8
+use obs_sequence_mod, only : obs_sequence_type, static_init_obs_sequence, &
+ read_obs_seq_header, destroy_obs_sequence, &
+ get_num_obs, write_obs_seq
+use utilities_mod, only : find_namelist_in_file, check_namelist_read, nc_check
+use obs_kind_mod, only : RADIOSONDE_U_WIND_COMPONENT, ACARS_U_WIND_COMPONENT, &
+ MARINE_SFC_U_WIND_COMPONENT, LAND_SFC_U_WIND_COMPONENT, &
+ SAT_U_WIND_COMPONENT, VORTEX_LAT
+use model_mod, only : static_init_model
+use netcdf
+
+implicit none
+
+! ----------------------------------------------------------------------
+! Declare namelist parameters
+! ----------------------------------------------------------------------
+
+! Generic parameters
+character(len=129) :: file_name_input = 'obs_seq.old', &
+ file_name_output = 'obs_seq.new', &
+ sonde_extra = 'obs_seq.rawin', &
+ acars_extra = 'obs_seq.acars', &
+ land_sfc_extra = 'obs_seq.land_sfc', &
+ marine_sfc_extra = 'obs_seq.marine_sfc', &
+ sat_wind_extra = 'obs_seq.satwnd', &
+ trop_cyclone_extra = 'obs_seq.tc'
+integer :: max_num_obs = 100000 ! Largest number of obs in one sequence
+
+! boundary-specific parameters
+real(r8) :: obs_boundary = 0.0_r8 ! number of grid points to remove obs near boundary
+logical :: increase_bdy_error = .false. ! true to increase obs error near boundary
+real(r8) :: maxobsfac = 2.5_r8 ! maximum increase in obs error near boundary
+real(r8) :: obsdistbdy = 15.0_r8 ! number of grid points to increase obs err.
+
+! parameters used to reduce observations
+logical :: sfc_elevation_check = .false. ! remove obs where model-obs topography is large
+real(r8) :: sfc_elevation_tol = 300.0_r8 ! largest difference between model and obs. topo.
+real(r8) :: obs_pressure_top = 0.0_r8 ! remove all obs at lower pressure
+real(r8) :: obs_height_top = 2.0e10_r8 ! remove all obs at higher height
+
+! Rawinsonde-specific parameters
+logical :: include_sig_data = .true. ! include significant-level data
+real(r8) :: tc_sonde_radii = -1.0_r8 ! remove sonde obs closer than this to TC
+
+! aircraft-specific parameters
+logical :: superob_aircraft = .false. ! super-ob aircraft data
+real(r8) :: aircraft_horiz_int = 36000.0_r8 ! horizontal interval for super-ob
+real(r8) :: aircraft_pres_int = 2500.0_r8 ! pressure interval for super-ob
+
+! sat wind specific parameters
+logical :: superob_sat_winds = .false. ! super-ob sat wind data
+real(r8) :: sat_wind_horiz_int = 100.0_r8 ! horizontal interval for super-ob
+real(r8) :: sat_wind_pres_int = 2500.0_r8 ! pressure interval for super-ob
+logical :: overwrite_ncep_satwnd_qc = .false. ! true to overwrite NCEP QC (see instructions)
+
+! surface obs. specific parameters
+logical :: overwrite_ncep_sfc_qc = .false. ! true to overwrite NCEP QC (see instructions)
+
+namelist /wrf_obs_preproc_nml/file_name_input, file_name_output, &
+ include_sig_data, superob_aircraft, superob_sat_winds, &
+ sfc_elevation_check, overwrite_ncep_sfc_qc, overwrite_ncep_satwnd_qc, &
+ aircraft_pres_int, sat_wind_pres_int, sfc_elevation_tol, &
+ obs_pressure_top, obs_height_top, obs_boundary, sonde_extra, &
+ acars_extra, land_sfc_extra, marine_sfc_extra, sat_wind_extra, &
+ trop_cyclone_extra, tc_sonde_radii, increase_bdy_error, &
+ maxobsfac, obsdistbdy, sat_wind_horiz_int, aircraft_horiz_int
+
+! ----------------------------------------------------------------------
+! Declare other variables
+! ----------------------------------------------------------------------
+
+character(len=129) :: obs_seq_read_format
+character(len=80) :: name
+
+integer :: io, iunit, fid, var_id, obs_seq_file_id, num_copies, &
+ num_qc, num_obs, max_obs_seq, nx, ny
+
+logical :: file_exist, pre_I_format
+
+type(obs_sequence_type) :: seq_all, seq_rawin, seq_sfc, seq_acars, seq_satwnd, &
+ seq_tc, seq_other
+
+call static_init_obs_sequence()
+call static_init_model()
+
+call find_namelist_in_file("input.nml", "wrf_obs_preproc_nml", iunit)
+read(iunit, nml = wrf_obs_preproc_nml, iostat = io)
+call check_namelist_read(iunit, io, "wrf_obs_preproc_nml")
+
+! open a wrfinput file, which is on this domain
+call nc_check( nf90_open(path = "wrfinput_d01", mode = nf90_nowrite, ncid = fid), &
+ 'main', 'open wrfinput_d01' )
+call nc_check( nf90_inq_dimid(fid, "west_east", var_id), &
+ 'main', 'inq. dimid west_east' )
+call nc_check( nf90_inquire_dimension(fid, var_id, name, nx), &
+ 'main', 'inquire dimension west_east' )
+call nc_check( nf90_inq_dimid(fid, "south_north", var_id), &
+ 'main', 'inq. dimid south_north' )
+call nc_check( nf90_inquire_dimension(fid, var_id, name, ny), &
+ 'main', 'inquire dimension south_north' )
+call nc_check( nf90_close(fid), 'main', 'close wrfinput_d01' )
+
+! if obs_seq file exists, read in the data, otherwise, create a blank one.
+inquire(file = trim(adjustl(file_name_input)), exist = file_exist)
+if ( file_exist ) then
+
+ call read_obs_seq_header(file_name_input, num_copies, num_qc, num_obs, max_obs_seq, &
+ obs_seq_file_id, obs_seq_read_format, pre_I_format, close_the_file = .true.)
+
+else
+
+ num_copies = 1 ; num_qc = 1 ; max_obs_seq = max_num_obs * 3
+ call create_new_obs_seq(num_copies, num_qc, max_obs_seq, seq_all)
+
+end if
+
+! create obs sequences for different obs types
+call create_new_obs_seq(num_copies, num_qc, max_num_obs, seq_rawin)
+call create_new_obs_seq(num_copies, num_qc, max_num_obs, seq_sfc)
+call create_new_obs_seq(num_copies, num_qc, max_num_obs, seq_acars)
+call create_new_obs_seq(num_copies, num_qc, max_num_obs, seq_satwnd)
+call create_new_obs_seq(num_copies, num_qc, 100, seq_tc)
+call create_new_obs_seq(num_copies, num_qc, max_num_obs, seq_other)
+
+! read input obs_seq file, divide into platforms
+call read_and_parse_input_seq(file_name_input, dble(nx), dble(ny), obs_boundary, &
+include_sig_data, obs_pressure_top, obs_height_top, sfc_elevation_check, &
+sfc_elevation_tol, overwrite_ncep_sfc_qc, overwrite_ncep_satwnd_qc, seq_rawin, &
+seq_sfc, seq_acars, seq_satwnd, seq_tc, seq_other)
+
+! add supplimental rawinsonde observations from file
+call add_supplimental_obs(sonde_extra, seq_rawin, max_obs_seq, &
+RADIOSONDE_U_WIND_COMPONENT, nx, ny, obs_boundary, include_sig_data, &
+obs_pressure_top, sfc_elevation_check, sfc_elevation_tol)
+
+! add supplimental ACARS observations from file
+call add_supplimental_obs(acars_extra, seq_acars, max_obs_seq, &
+ACARS_U_WIND_COMPONENT, nx, ny, obs_boundary, include_sig_data, &
+obs_pressure_top, sfc_elevation_check, sfc_elevation_tol)
+
+! add supplimental marine observations from file
+call add_supplimental_obs(marine_sfc_extra, seq_sfc, max_obs_seq, &
+MARINE_SFC_U_WIND_COMPONENT, nx, ny, obs_boundary, include_sig_data, &
+obs_pressure_top, sfc_elevation_check, sfc_elevation_tol)
+
+! add supplimental land surface observations from file
+call add_supplimental_obs(land_sfc_extra, seq_sfc, max_obs_seq, &
+LAND_SFC_U_WIND_COMPONENT, nx, ny, obs_boundary, include_sig_data, &
+obs_pressure_top, sfc_elevation_check, sfc_elevation_tol)
+
+! add supplimental satellite wind observations from file
+call add_supplimental_obs(sat_wind_extra, seq_satwnd, max_obs_seq, &
+SAT_U_WIND_COMPONENT, nx, ny, obs_boundary, include_sig_data, &
+obs_pressure_top, sfc_elevation_check, sfc_elevation_tol)
+
+! add supplimental tropical cyclone vortex observations from file
+call add_supplimental_obs(trop_cyclone_extra, seq_tc, max_obs_seq, &
+VORTEX_LAT, nx, ny, obs_boundary, include_sig_data, &
+obs_pressure_top, sfc_elevation_check, sfc_elevation_tol)
+
+! remove all sonde observations within radius of TC if desired
+if ( tc_sonde_radii > 0.0_r8 ) call remove_sondes_near_tc(seq_tc, &
+ seq_rawin, tc_sonde_radii)
+
+! super-ob ACARS data
+if ( superob_aircraft ) call superob_aircraft_data(seq_acars, &
+ aircraft_horiz_int, aircraft_pres_int)
+
+! super-ob satellite wind data
+if ( superob_sat_winds ) call superob_sat_wind_data(seq_satwnd, &
+ sat_wind_horiz_int, sat_wind_pres_int)
+
+max_obs_seq = get_num_obs(seq_tc) + get_num_obs(seq_rawin) + &
+ get_num_obs(seq_sfc) + get_num_obs(seq_acars) + &
+ get_num_obs(seq_satwnd) + get_num_obs(seq_other)
+
+call create_new_obs_seq(num_copies, num_qc, max_obs_seq, seq_all)
+
+call build_master_sequence(seq_tc, seq_all)
+call destroy_obs_sequence(seq_tc)
+
+call build_master_sequence(seq_rawin, seq_all)
+call destroy_obs_sequence(seq_rawin)
+
+call build_master_sequence(seq_sfc, seq_all)
+call destroy_obs_sequence(seq_sfc)
+
+call build_master_sequence(seq_acars, seq_all)
+call destroy_obs_sequence(seq_acars)
+
+call build_master_sequence(seq_satwnd, seq_all)
+call destroy_obs_sequence(seq_satwnd)
+
+call build_master_sequence(seq_other, seq_all)
+call destroy_obs_sequence(seq_other)
+
+write(6,*) 'Total number of observations:', get_num_obs(seq_all)
+
+! increase the observation error along the lateral boundary
+if ( increase_bdy_error ) call increase_obs_err_bdy(seq_all, &
+ obsdistbdy, maxobsfac, dble(nx), dble(ny))
+
+! write the observation sequence to file
+call write_obs_seq(seq_all, file_name_output)
+call destroy_obs_sequence(seq_all)
+
+stop
+end
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! aircraft_obs_check - function that determines whether to include an
+! aircraft observation in the sequence. For now,
+! this function is a placeholder and returns true.
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+function aircraft_obs_check()
+
+use types_mod, only : r8
+
+implicit none
+
+logical :: aircraft_obs_check
+
+aircraft_obs_check = .true.
+
+return
+end function aircraft_obs_check
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! add_supplimental_obs - subroutine that reads observation data from
+! a supplimental obs sequence file, performs
+! validation checks and adds it to the
+! platform-specific obs sequence.
+!
+! filename - name of supplimental obs sequence file
+! obs_seq - platform-specific obs sequence
+! max_obs_seq - maximum number of observations in sequence
+! plat_kind - integer kind of platform (used for print statements)
+! nx - number of grid points in x direction
+! ny - number of grid points in y direction
+! obs_bdy - grid point buffer to remove observations
+! siglevel - true to include sonde significant level data
+! ptop - lowest pressure to include in sequence
+! htop - highest height level to include in sequence
+! sfcelev - true to perform surface obs. elevation check
+! elev_max - maximum difference between model and obs. height
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine add_supplimental_obs(filename, obs_seq, max_obs_seq, plat_kind, &
+ nx, ny, obs_bdy, siglevel, ptop, htop, &
+ sfcelev, elev_max)
+
+use types_mod, only : r8
+use location_mod, only : location_type, get_location, vert_is_pressure, &
+ vert_is_height
+use obs_sequence_mod, only : obs_sequence_type, obs_type, init_obs, &
+ get_num_copies, get_num_qc, read_obs_seq, &
+ get_first_obs, get_obs_def, get_next_obs, &
+ copy_obs, append_obs_to_seq, destroy_obs_sequence
+use obs_def_mod, only : obs_def_type, get_obs_kind, &
+ get_obs_def_location
+use obs_kind_mod, only : RADIOSONDE_U_WIND_COMPONENT, ACARS_U_WIND_COMPONENT, &
+ LAND_SFC_U_WIND_COMPONENT, MARINE_SFC_U_WIND_COMPONENT, &
+ SAT_U_WIND_COMPONENT, VORTEX_LAT
+use model_mod, only : get_domain_info
+
+implicit none
+
+character(len=129), intent(in) :: filename
+type(obs_sequence_type), intent(inout) :: obs_seq
+integer, intent(in) :: max_obs_seq, plat_kind, nx, ny
+logical, intent(in) :: siglevel, sfcelev
+real(r8), intent(in) :: obs_bdy, ptop, htop, elev_max
+
+integer :: nloc, okind, dom_id
+logical :: file_exist, last_obs, pass_checks, original_observation, &
+ rawinsonde_obs_check, aircraft_obs_check, surface_obs_check, &
+ sat_wind_obs_check
+real(r8) :: xyz_loc(3), xloc, yloc
+
+type(location_type) :: obs_loc_list(max_obs_seq), obs_loc
+type(obs_def_type) :: obs_def
+type(obs_sequence_type) :: supp_obs_seq
+type(obs_type) :: obs_in, prev_obs, obs
+
+inquire(file = trim(adjustl(filename)), exist = file_exist)
+if ( .not. file_exist ) return
+
+select case (plat_kind)
+
+ case (RADIOSONDE_U_WIND_COMPONENT)
+ write(6,*) 'Adding Supplimental Rawinsonde Data'
+ case (ACARS_U_WIND_COMPONENT)
+ write(6,*) 'Adding Supplimental ACARS Data'
+ case (MARINE_SFC_U_WIND_COMPONENT)
+ write(6,*) 'Adding Supplimental Marine Surface Data'
+ case (LAND_SFC_U_WIND_COMPONENT)
+ write(6,*) 'Adding Supplimental Land Surface Data'
+ case (SAT_U_WIND_COMPONENT)
+ write(6,*) 'Adding Supplimental Satellite Wind Data'
+ case (VORTEX_LAT)
+ write(6,*) 'Adding Supplimental Tropical Cyclone Data'
+
+end select
+
+call init_obs(obs_in, get_num_copies(obs_seq), get_num_qc(obs_seq))
+call init_obs(obs, get_num_copies(obs_seq), get_num_qc(obs_seq))
+call init_obs(prev_obs, get_num_copies(obs_seq), get_num_qc(obs_seq))
+
+! create list of observations in plaform sequence
+call build_obs_loc_list(obs_seq, max_obs_seq, nloc, obs_loc_list)
+
+last_obs = .false.
+call read_obs_seq(trim(adjustl(filename)), 0, 0, 0, supp_obs_seq)
+if ( .not. get_first_obs(supp_obs_seq, obs_in) ) last_obs = .true.
+
+ObsLoop: do while ( .not. last_obs ) ! loop over all observations in a sequence
+
+ ! read data from observation
+ call get_obs_def(obs_in, obs_def)
+ okind = get_obs_kind(obs_def)
+ obs_loc = get_obs_def_location(obs_def)
+ xyz_loc = get_location(obs_loc)
+ call get_domain_info(xyz_loc(1),xyz_loc(2),dom_id,xloc,yloc)
+
+ ! check if the observation is within the domain
+ if ( ((xloc < (obs_bdy+1.0_r8) .or. xloc > (dble(nx)-obs_bdy-1.0_r8) .or. &
+ yloc < (obs_bdy+1.0_r8) .or. yloc > (dble(ny)-obs_bdy-1.0_r8)) .and. &
+ (dom_id == 1)) .or. dom_id < 1 ) then
+
+ prev_obs = obs_in
+ call get_next_obs(supp_obs_seq, prev_obs, obs_in, last_obs)
+ cycle ObsLoop
+
+ end if
+
+ ! check if the observation is within vertical bounds of domain
+ if ( (vert_is_pressure(obs_loc) .and. xyz_loc(3) < ptop) .or. &
+ (vert_is_height(obs_loc) .and. xyz_loc(3) > htop) ) then
+
+ prev_obs = obs_in
+ call get_next_obs(supp_obs_seq, prev_obs, obs_in, last_obs)
+ cycle ObsLoop
+
+ end if
+
+ ! check if the observation already exists
+ if ( .not. original_observation(obs_loc, obs_loc_list, nloc) ) then
+
+ prev_obs = obs_in
+ call get_next_obs(supp_obs_seq, prev_obs, obs_in, last_obs)
+ cycle ObsLoop
+
+ end if
+
+ ! perform platform-specific checks
+ select case (plat_kind)
+
+ case (RADIOSONDE_U_WIND_COMPONENT)
+ pass_checks = rawinsonde_obs_check(obs_loc, okind, siglevel, &
+ sfcelev, elev_max)
+ case (ACARS_U_WIND_COMPONENT)
+ pass_checks = aircraft_obs_check()
+ case (MARINE_SFC_U_WIND_COMPONENT)
+ pass_checks = surface_obs_check(sfcelev, elev_max, xyz_loc)
+ case (LAND_SFC_U_WIND_COMPONENT)
+ pass_checks = surface_obs_check(sfcelev, elev_max, xyz_loc)
+ case (SAT_U_WIND_COMPONENT)
+ pass_checks = sat_wind_obs_check()
+ case default
+ pass_checks = .true.
+
+ end select
+
+ if ( pass_checks ) then
+
+ call copy_obs(obs, obs_in)
+ call append_obs_to_seq(obs_seq, obs)
+
+ end if
+
+ prev_obs = obs_in
+ call get_next_obs(supp_obs_seq, prev_obs, obs_in, last_obs)
+
+end do ObsLoop
+
+call destroy_obs_sequence(supp_obs_seq)
+
+return
+end subroutine add_supplimental_obs
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! create_new_obs_seq - subroutine that is used to create a new
+! observation sequence.
+!
+! num_copies - number of copies associated with each observation
+! num_qc - number of quality control reports in each obs.
+! max_num - maximum number of observations in sequence
+! seq - observation sequence
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine create_new_obs_seq(num_copies, num_qc, max_num, seq)
+
+use obs_sequence_mod, only : obs_sequence_type, init_obs_sequence, &
+ set_copy_meta_data, set_qc_meta_data
+
+implicit none
+
+integer, intent(in) :: num_copies, num_qc, max_num
+type(obs_sequence_type), intent(out) :: seq
+
+character(len=129) :: copy_meta_data, qc_meta_data
+integer :: i
+
+call init_obs_sequence(seq, num_copies, num_qc, max_num)
+do i = 1, num_copies
+ copy_meta_data = 'NCEP BUFR observation'
+ call set_copy_meta_data(seq, i, copy_meta_data)
+end do
+do i = 1, num_qc
+ qc_meta_data = 'NCEP QC index'
+ call set_qc_meta_data(seq, i, qc_meta_data)
+end do
+
+return
+end subroutine create_new_obs_seq
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! build_master_sequence - subroutine used to take observations from
+! a smaller observation sequence and appends
+! them to a larger observation sequence.
+! Note that this routine only works if the
+! observations are at the same time.
+!
+! seq_type - observation sequence with one observation type
+! seq_all - observation sequence with more observations
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine build_master_sequence(seq_type, seq_all)
+
+use obs_sequence_mod, only : obs_type, obs_sequence_type, init_obs, &
+ get_first_obs, copy_obs, append_obs_to_seq, &
+ get_next_obs, get_obs_def, get_num_copies, &
+ get_num_qc
+
+implicit none
+
+type(obs_sequence_type), intent(in) :: seq_type
+type(obs_sequence_type), intent(inout) :: seq_all
+
+logical :: last_obs
+type(obs_type) :: obs_in, obs, prev_obs
+
+last_obs = .false.
+call init_obs(obs_in, get_num_copies(seq_type), get_num_qc(seq_type))
+call init_obs(obs, get_num_copies(seq_type), get_num_qc(seq_type))
+call init_obs(prev_obs, get_num_copies(seq_type), get_num_qc(seq_type))
+
+if ( .not. get_first_obs(seq_type, obs_in) ) return
+
+do while ( .not. last_obs )
+
+ call copy_obs(obs, obs_in)
+ call append_obs_to_seq(seq_all, obs)
+
+ prev_obs = obs_in
+ call get_next_obs(seq_type, prev_obs, obs_in, last_obs)
+
+end do
+
+return
+end subroutine build_master_sequence
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! build_obs_loc_list - subroutine that creates an array of locations
+! of the observations in a sequence.
+!
+! obs_seq - observation sequence to read locations from
+! maxobs - maximum number of observations in a sequence
+! nloc - number of individual locations
+! obs_loc_list - array of observation locations
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine build_obs_loc_list(seq, maxobs, nloc, obs_loc_list)
+
+use location_mod, only : location_type
+use obs_sequence_mod, only : obs_sequence_type, get_num_copies, &
+ get_num_qc, init_obs, get_first_obs, &
+ get_obs_def, get_next_obs, obs_type
+use obs_def_mod, only : obs_def_type, get_obs_def_location
+
+implicit none
+
+integer, intent(in) :: maxobs
+type(obs_sequence_type), intent(in) :: seq
+integer, intent(out) :: nloc
+type(location_type), intent(out) :: obs_loc_list(maxobs)
+
+logical :: last_obs, original_observation
+type(obs_type) :: obs, prev_obs
+type(obs_def_type) :: obs_def
+type(location_type) :: obs_loc
+
+call init_obs(obs, get_num_copies(seq), get_num_qc(seq))
+call init_obs(prev_obs, get_num_copies(seq), get_num_qc(seq))
+
+last_obs = .false. ; nloc = 0
+if ( .not. get_first_obs(seq, obs) ) last_obs = .true.
+
+do while ( .not. last_obs ) ! loop over all observations in a sequence
+
+ call get_obs_def(obs, obs_def)
+ obs_loc = get_obs_def_location(obs_def)
+
+ ! construct a list of observation locations
+ if ( original_observation(obs_loc, obs_loc_list, nloc) ) then
+
+ nloc = nloc + 1
+ obs_loc_list(nloc) = obs_loc
+
+ end if
+ prev_obs = obs
+ call get_next_obs(seq, prev_obs, obs, last_obs)
+
+end do
+
+return
+end subroutine build_obs_loc_list
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! create_obs_type - subroutine that is used to create an observation
+! type from observation data.
+!
+! lat - latitude of observation
+! lon - longitude of observation
+! vloc - vertical location of observation
+! vcord - DART vertical coordinate integer
+! obsv - observation value
+! okind - observation kind
+! oerr - observation error
+! day - gregorian day of the observation
+! sec - gregorian second of the observation
+! qc - integer quality control value
+! obs - observation type that includes the observation information
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine create_obs_type(lat, lon, vloc, vcord, obsv, okind, oerr, qc, otime, obs)
+
+use types_mod, only : r8
+use obs_sequence_mod, only : obs_type, set_obs_values, set_qc, set_obs_def
+use obs_def_mod, only : obs_def_type, set_obs_def_time, set_obs_def_kind, &
+ set_obs_def_error_variance, set_obs_def_location
+use location_mod, only : location_type, set_location
+use time_manager_mod, only : time_type
+
+implicit none
+
+integer, intent(in) :: okind, vcord
+real(r8), intent(in) :: lat, lon, vloc, obsv, oerr, qc
+type(time_type), intent(in) :: otime
+type(obs_type), intent(inout) :: obs
+
+real(r8) :: obs_val(1), qc_val(1)
+type(obs_def_type) :: obs_def
+
+call set_obs_def_location(obs_def, set_location(lon, lat, vloc, vcord))
+call set_obs_def_kind(obs_def, okind)
+call set_obs_def_time(obs_def, otime)
+call set_obs_def_error_variance(obs_def, oerr)
+call set_obs_def(obs, obs_def)
+
+obs_val(1) = obsv
+call set_obs_values(obs, obs_val)
+qc_val(1) = qc
+call set_qc(obs, qc_val)
+
+return
+end subroutine create_obs_type
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! increase_obs_err_bdy - subroutine that increases the observation
+! error based on proximity to the lateral
+! boundary.
+!
+! seq - observation sequence
+! obsbdy - number of grid points near boundary to increase error
+! maxfac - factor to increase observation error at boundary
+! nx - number of grid points in the x direction
+! ny - number of grid points in the y direction
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine increase_obs_err_bdy(seq, obsbdy, maxfac, nx, ny)
+
+use types_mod, only : r8
+use location_mod, only : location_type, get_location
+use obs_sequence_mod, only : obs_sequence_type, obs_type, init_obs, &
+ get_num_copies, get_num_qc, get_first_obs, &
+ get_obs_def, set_obs_def, set_obs, &
+ get_next_obs, get_obs_key
+use obs_def_mod, only : obs_def_type, get_obs_def_error_variance, &
+ set_obs_def_error_variance, get_obs_def_location
+use model_mod, only : get_domain_info
+
+implicit none
+
+type(obs_sequence_type), intent(inout) :: seq
+real(r8), intent(in) :: obsbdy, maxfac, nx, ny
+
+integer :: dom_id
+logical :: last_obs
+real(r8) :: mobse, bobse, xyz_loc(3), xloc, yloc, bdydist, obsfac
+
+type(obs_def_type) :: obs_def
+type(obs_type) :: obs, prev_obs
+
+write(6,*) 'Increasing the Observation Error Near the Lateral Boundary'
+
+call init_obs(obs, get_num_copies(seq), get_num_qc(seq))
+call init_obs(prev_obs, get_num_copies(seq), get_num_qc(seq))
+
+! compute slope and intercept for error increase factor
+mobse = (maxfac - 1.0_r8) / (1.0_r8 - obsbdy)
+bobse = maxfac - mobse
+
+last_obs = .false.
+if ( .not. get_first_obs(seq, obs) ) last_obs = .true.
+
+do while ( .not. last_obs )
+
+ ! get location information
+ call get_obs_def(obs, obs_def)
+ xyz_loc = get_location(get_obs_def_location(obs_def))
+ call get_domain_info(xyz_loc(1),xyz_loc(2),dom_id,xloc,yloc)
+
+ ! compute distance to boundary, increase based on this distance
+ bdydist = min(xloc-1.0_r8, yloc-1.0_r8, nx-xloc, ny-yloc)
+ if ( bdydist <= obsbdy ) then
+
+ obsfac = mobse * bdydist + bobse
+ call set_obs_def_error_variance(obs_def, &
+ get_obs_def_error_variance(obs_def) * obsfac * obsfac)
+ call set_obs_def(obs, obs_def)
+ call set_obs(seq, obs, get_obs_key(obs))
+
+ end if
+ prev_obs = obs
+ call get_next_obs(seq, prev_obs, obs, last_obs)
+
+end do
+
+return
+end subroutine increase_obs_err_bdy
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! isManLevel - function that returns a logical true if the input
+! pressure level is a mandatory rawinsonde level.
+!
+! plevel - pressure level to check (Pa)
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+function isManLevel(plevel)
+
+use types_mod, only : r8
+
+implicit none
+
+real(r8), intent(in) :: plevel
+
+integer, parameter :: nman = 16
+integer :: kk
+logical :: isManLevel
+real (r8) raw_man_levels(nman) &
+ / 100000.0_r8, 92500.0_r8, 85000.0_r8, 70000.0_r8, 50000.0_r8, 40000.0_r8, &
+ 30000.0_r8, 25000.0_r8, 20000.0_r8, 15000.0_r8, 10000.0_r8, 7000.0_r8, &
+ 5000.0_r8, 3000.0_r8, 2000.0_r8, 1000.0_r8 /
+
+isManLevel = .false.
+do kk = 1, nman
+ if ( plevel == raw_man_levels(kk) ) then
+ isManLevel = .true.
+ return
+ end if
+end do
+
+return
+end function isManLevel
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! original_observation - function that returns true if the location
+! is not within an array of locations
+!
+! obsloc - location to check
+! obsloc_list - array of locations to look through
+! nloc - number of locations in array
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+function original_observation(obsloc, obsloc_list, nloc)
+
+use types_mod, only : r8
+use location_mod, only : location_type, get_dist
+
+real(r8), parameter :: dist_epsilon = 0.00001_r8
+
+integer, intent(in) :: nloc
+type(location_type), intent(in) :: obsloc, obsloc_list(nloc)
+
+integer :: n
+logical :: original_observation
+
+original_observation = .true.
+
+do n = 1, nloc
+
+ if ( get_dist(obsloc, obsloc_list(n), 1, 1, .true.) <= dist_epsilon ) then
+ original_observation = .false.
+ return
+ end if
+
+end do
+
+return
+end function original_observation
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! rawinsonde_obs_check - function that performs obsrvation checks
+! specific to rawinsonde observations.
+!
+! obs_loc - observation location
+! obs_kind - DART observation kind
+! siglevel - true to include significant level data
+! elev_check - true to check differene between model and obs elev.
+! elev_max - maximum difference between model and obs elevation
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+function rawinsonde_obs_check(obs_loc, obs_kind, siglevel, &
+ elev_check, elev_max)
+
+use types_mod, only : r8
+use obs_kind_mod, only : RADIOSONDE_SURFACE_ALTIMETER, KIND_SURFACE_ELEVATION
+use model_mod, only : model_interpolate
+use location_mod, only : location_type, set_location, get_location, &
+ vert_is_pressure
+
+implicit none
+
+type(location_type), intent(in) :: obs_loc
+integer, intent(in) :: obs_kind
+logical, intent(in) :: siglevel, elev_check
+real(r8), intent(in) :: elev_max
+
+integer :: istatus
+logical :: rawinsonde_obs_check, isManLevel
+real(r8) :: xyz_loc(3), xmod(1), hsfc
+
+rawinsonde_obs_check = .true.
+xyz_loc = get_location(obs_loc)
+
+if ( obs_kind /= RADIOSONDE_SURFACE_ALTIMETER ) then
+
+ ! check if vertical level is mandatory level
+ if ( (.not. siglevel) .and. (.not. isManLevel(xyz_loc(3))) ) then
+ rawinsonde_obs_check = .false.
+ return
+ end if
+
+else
+
+ ! perform elevation check for altimeter
+ if ( elev_check ) then
+
+ call model_interpolate(xmod, obs_loc, KIND_SURFACE_ELEVATION, hsfc, istatus)
+ if ( abs(hsfc - xyz_loc(3)) > elev_max ) rawinsonde_obs_check = .false.
+
+ end if
+
+end if
+
+return
+end function rawinsonde_obs_check
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! read_and_parse_input_seq - subroutine that reads a generic
+! observation sequence and divides the
+! obs into sequences for each platform.
+!
+! filename - name of input obs sequence
+! nx - number of grid points in x direction
+! ny - number of grid points in y direction
+! obs_bdy - grid point buffer to remove observations
+! siglevel - true to include sonde significant level data
+! ptop - lowest pressure to include in sequence
+! htop - highest height level to include in sequence
+! sfcelev - true to perform surface obs. elevation check
+! elev_max - maximum difference between model and obs. height
+! new_sfc_qc - true to replace NCEP surface QC
+! new_satwnd_qc - true to replace NCEP sat wind QC over ocean
+! rawin_seq - rawinsonde sequence
+! sfc_seq - surface sequence
+! acars_seq - aircraft sequence
+! satwnd_seq - satellite wind sequence
+! tc_seq - TC data sequence
+! other_seq - remaining observation sequence
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine read_and_parse_input_seq(filename, nx, ny, obs_bdy, siglevel, ptop, &
+ htop, sfcelev, elev_max, new_sfc_qc, &
+ new_satwnd_qc, rawin_seq, sfc_seq, &
+ acars_seq, satwnd_seq, tc_seq, other_seq)
+
+use types_mod, only : r8
+use utilities_mod, only : nc_check
+use location_mod, only : location_type, get_location, vert_is_pressure, &
+ vert_is_height
+use obs_sequence_mod, only : obs_sequence_type, obs_type, init_obs, &
+ get_num_copies, get_num_qc, get_qc_meta_data, &
+ get_first_obs, get_obs_def, copy_obs, get_num_qc, &
+ append_obs_to_seq, get_next_obs, get_qc, set_qc, &
+ destroy_obs_sequence, read_obs_seq
+use obs_def_mod, only : obs_def_type, get_obs_kind, get_obs_def_location
+use obs_kind_mod, only : RADIOSONDE_U_WIND_COMPONENT, RADIOSONDE_V_WIND_COMPONENT, &
+ RADIOSONDE_SURFACE_ALTIMETER, RADIOSONDE_TEMPERATURE, &
+ RADIOSONDE_SPECIFIC_HUMIDITY, &
+ AIRCRAFT_U_WIND_COMPONENT, AIRCRAFT_V_WIND_COMPONENT, &
+ AIRCRAFT_TEMPERATURE, AIRCRAFT_SPECIFIC_HUMIDITY, &
+ ACARS_U_WIND_COMPONENT, ACARS_V_WIND_COMPONENT, &
+ ACARS_TEMPERATURE, ACARS_SPECIFIC_HUMIDITY, &
+ MARINE_SFC_U_WIND_COMPONENT, MARINE_SFC_V_WIND_COMPONENT, &
+ MARINE_SFC_TEMPERATURE, MARINE_SFC_SPECIFIC_HUMIDITY, &
+ LAND_SFC_U_WIND_COMPONENT, LAND_SFC_V_WIND_COMPONENT, &
+ LAND_SFC_TEMPERATURE, LAND_SFC_SPECIFIC_HUMIDITY, &
+ MARINE_SFC_ALTIMETER, LAND_SFC_ALTIMETER, &
+ SAT_U_WIND_COMPONENT, SAT_V_WIND_COMPONENT, &
+ VORTEX_LAT, VORTEX_LON, VORTEX_PMIN, VORTEX_WMAX
+use model_mod, only : get_domain_info
+use netcdf
+
+implicit none
+
+real(r8), parameter :: satwnd_qc_ok = 15.0_r8
+real(r8), parameter :: sfc_qc_ok = 9.0_r8
+real(r8), parameter :: new_qc_value = 2.0_r8
+
+character(len=129), intent(in) :: filename
+real(r8), intent(in) :: nx, ny, obs_bdy, ptop, htop, elev_max
+logical, intent(in) :: siglevel, sfcelev, new_sfc_qc, &
+ new_satwnd_qc
+type(obs_sequence_type), intent(inout) :: rawin_seq, sfc_seq, acars_seq, &
+ satwnd_seq, tc_seq, other_seq
+
+character(len=129) :: qcmeta
+integer :: fid, var_id, okind, dom_id, i, j
+logical :: file_exist, last_obs, input_ncep_qc, rawinsonde_obs_check, &
+ surface_obs_check, aircraft_obs_check, sat_wind_obs_check
+real(r8), allocatable :: xland(:,:), qc(:)
+real(r8) :: xyz_loc(3), xloc, yloc
+
+type(location_type) :: obs_loc
+type(obs_def_type) :: obs_def
+type(obs_sequence_type) :: seq
+type(obs_type) :: obs, obs_in, prev_obs
+
+inquire(file = trim(adjustl(filename)), exist = file_exist)
+if ( .not. file_exist ) return
+
+call read_obs_seq(filename, 0, 0, 0, seq)
+
+call init_obs(obs, get_num_copies(seq), get_num_qc(seq))
+call init_obs(obs_in, get_num_copies(seq), get_num_qc(seq))
+call init_obs(prev_obs, get_num_copies(seq), get_num_qc(seq))
+allocate(qc(get_num_qc(seq)))
+
+! read land distribution
+allocate(xland(nint(nx),nint(ny)))
+call nc_check( nf90_open(path = "wrfinput_d01", mode = nf90_nowrite, ncid = fid), &
+ 'read_and_parse_input_seq', 'open wrfinput_d01')
+call nc_check( nf90_inq_varid(fid, "XLAND", var_id), &
+ 'read_and_parse_input_seq', 'inquire XLAND ID')
+call nc_check( nf90_get_var(fid, var_id, xland), &
+ 'read_and_parse_input_seq', 'read XLAND')
+call nc_check( nf90_close(fid), 'read_and_parse_input_seq', 'close wrfinput_d01')
+
+input_ncep_qc = .false.
+qcmeta = get_qc_meta_data(seq, 1)
+if ( trim(adjustl(qcmeta)) == 'NCEP QC index' ) input_ncep_qc = .true.
+
+last_obs = .false.
+if ( .not. get_first_obs(seq, obs_in) ) last_obs = .true.
+
+InputObsLoop: do while ( .not. last_obs ) ! loop over all observations in a sequence
+
+ ! Get the observation information, check if it is in the domain
+ call get_obs_def(obs_in, obs_def)
+ okind = get_obs_kind(obs_def)
+ obs_loc = get_obs_def_location(obs_def)
+ xyz_loc = get_location(obs_loc)
+ call get_domain_info(xyz_loc(1),xyz_loc(2),dom_id,xloc,yloc)
+ i = nint(xloc); j = nint(yloc)
+
+ ! check horizontal location
+ if ( ((xloc < (obs_bdy+1.0_r8) .or. xloc > (nx-obs_bdy-1.0_r8) .or. &
+ yloc < (obs_bdy+1.0_r8) .or. yloc > (ny-obs_bdy-1.0_r8)) .and. &
+ (dom_id == 1)) .or. dom_id < 1 ) then
+
+ prev_obs = obs_in
+ call get_next_obs(seq, prev_obs, obs_in, last_obs)
+ cycle InputObsLoop
+
+ end if
+
+ ! check vertical location
+ if ( (vert_is_pressure(obs_loc) .and. xyz_loc(3) < ptop) .or. &
+ (vert_is_height(obs_loc) .and. xyz_loc(3) > htop) ) then
+
+ prev_obs = obs_in
+ call get_next_obs(seq, prev_obs, obs_in, last_obs)
+ cycle InputObsLoop
+
+ end if
+
+ ! perform platform-specific checks
+ select case (okind)
+
+ case ( RADIOSONDE_U_WIND_COMPONENT, RADIOSONDE_V_WIND_COMPONENT, &
+ RADIOSONDE_TEMPERATURE, RADIOSONDE_SPECIFIC_HUMIDITY, &
+ RADIOSONDE_SURFACE_ALTIMETER)
+
+ if ( rawinsonde_obs_check(obs_loc, okind, siglevel, sfcelev, elev_max) ) then
+
+ call copy_obs(obs, obs_in)
+ call append_obs_to_seq(rawin_seq, obs)
+
+ end if
+
+ case ( LAND_SFC_U_WIND_COMPONENT, LAND_SFC_V_WIND_COMPONENT, &
+ LAND_SFC_TEMPERATURE, LAND_SFC_SPECIFIC_HUMIDITY, &
+ MARINE_SFC_U_WIND_COMPONENT, MARINE_SFC_V_WIND_COMPONENT, &
+ MARINE_SFC_TEMPERATURE, MARINE_SFC_SPECIFIC_HUMIDITY, &
+ LAND_SFC_ALTIMETER, MARINE_SFC_ALTIMETER )
+
+ if ( surface_obs_check(sfcelev, elev_max, xyz_loc) ) then
+
+ call copy_obs(obs, obs_in)
+ if ( new_sfc_qc .and. okind /= LAND_SFC_ALTIMETER .and. &
+ okind /= MARINE_SFC_ALTIMETER ) then
+
+ call get_qc(obs, qc)
+ if ( qc(1) == sfc_qc_ok .and. input_ncep_qc ) then
+ qc(1) = new_qc_value
+ call set_qc(obs, qc)
+ end if
+
+ end if
+ call append_obs_to_seq(sfc_seq, obs)
+
+ endif
+
+ case ( AIRCRAFT_U_WIND_COMPONENT, AIRCRAFT_V_WIND_COMPONENT, &
+ AIRCRAFT_TEMPERATURE, AIRCRAFT_SPECIFIC_HUMIDITY, &
+ ACARS_U_WIND_COMPONENT, ACARS_V_WIND_COMPONENT, &
+ ACARS_TEMPERATURE, ACARS_SPECIFIC_HUMIDITY )
+
+ if ( aircraft_obs_check() ) then
+
+ call copy_obs(obs, obs_in)
+ call append_obs_to_seq(acars_seq, obs)
+
+ end if
+
+ case ( SAT_U_WIND_COMPONENT, SAT_V_WIND_COMPONENT )
+
+ if ( sat_wind_obs_check() ) then
+
+ call copy_obs(obs, obs_in)
+ if ( new_satwnd_qc ) then
+
+ call get_qc(obs, qc)
+ if ( qc(1) == satwnd_qc_ok .and. input_ncep_qc .and. &
+ xland(i,j) > 1.0_r8 ) then
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list