[Dart-dev] [7234] DART/trunk/models/mpas_atm: first commit for new code from soyoung.
nancy at ucar.edu
nancy at ucar.edu
Thu Oct 30 16:26:23 MDT 2014
Revision: 7234
Author: nancy
Date: 2014-10-30 16:26:22 -0600 (Thu, 30 Oct 2014)
Log Message:
-----------
first commit for new code from soyoung. the program
mpas_art_obs_preprocess is based on the wrf_dart_obs_preprocess
donated by ryan torn. this version does some of the same
types of preprocessing, including superob'ing (averaging)
dense obs like satellite winds and aircraft obs. one input
to this program is an mpas grid file. the superobs will be
done on all obs inside a single mpas grid cell. to adjust
the density of the obs, specify a coarser or finer mpas grid.
the obs are also averaged in the vertical; the size of the
vertical bins is set via namelist value.
this code is new and has not been rigoriously tested; please
use with caution, examine the outputs, and report any problems
to us? thanks.
Modified Paths:
--------------
DART/trunk/models/mpas_atm/model_mod.f90
DART/trunk/models/mpas_atm/work/input.nml
Added Paths:
-----------
DART/trunk/models/mpas_atm/mpas_dart_obs_preprocess.f90
DART/trunk/models/mpas_atm/mpas_dart_obs_preprocess.html
DART/trunk/models/mpas_atm/mpas_dart_obs_preprocess.nml
DART/trunk/models/mpas_atm/work/mkmf_mpas_dart_obs_preprocess
DART/trunk/models/mpas_atm/work/path_names_mpas_dart_obs_preprocess
-------------- next part --------------
Modified: DART/trunk/models/mpas_atm/model_mod.f90
===================================================================
--- DART/trunk/models/mpas_atm/model_mod.f90 2014-10-30 21:02:23 UTC (rev 7233)
+++ DART/trunk/models/mpas_atm/model_mod.f90 2014-10-30 22:26:22 UTC (rev 7234)
@@ -120,7 +120,8 @@
get_analysis_time, &
write_model_time, &
get_grid_dims, &
- print_variable_ranges
+ print_variable_ranges, &
+ find_closest_cell_center
! version controlled file description for error handling, do not edit
character(len=256), parameter :: source = &
Added: DART/trunk/models/mpas_atm/mpas_dart_obs_preprocess.f90
===================================================================
--- DART/trunk/models/mpas_atm/mpas_dart_obs_preprocess.f90 (rev 0)
+++ DART/trunk/models/mpas_atm/mpas_dart_obs_preprocess.f90 2014-10-30 22:26:22 UTC (rev 7234)
@@ -0,0 +1,2078 @@
+! DART software - Copyright 2004 - 2013 UCAR. This open source software is
+! provided by UCAR, "as is", without charge, subject to all terms of use at
+! http://www.image.ucar.edu/DAReS/DART/DART_download
+!
+! $Id$
+
+program mpas_dart_obs_preprocess
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! mpas_dart_obs_preprocess - MPAS-DART utility program that adds
+! observations from supplimental obs sequences.
+! The program assumes all data is from one
+! obsevation time. In addition, this program allows
+! the user to do the following functions:
+!
+! - 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 based on wrf_dart_obs_preprocess by Soyoung Ha, NCAR/MMM (Oct. 2014)
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+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, &
+ METAR_U_10_METER_WIND, GPSRO_REFRACTIVITY, &
+ SAT_U_WIND_COMPONENT, PROFILER_U_WIND_COMPONENT, VORTEX_LAT
+use time_manager_mod, only : time_type, set_calendar_type, GREGORIAN, set_time
+use model_mod, only : static_init_model
+use netcdf
+
+implicit none
+
+! ----------------------------------------------------------------------
+! Declare namelist parameters
+! ----------------------------------------------------------------------
+
+! Generic parameters
+character(len=256) :: static_mpas_filename = 'mpas_init.nc'
+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', &
+ metar_extra = 'obs_seq.metar', &
+ marine_sfc_extra = 'obs_seq.marine', &
+ sat_wind_extra = 'obs_seq.satwnd', &
+ profiler_extra = 'obs_seq.profiler', &
+ gpsro_extra = 'obs_seq.gpsro', &
+ trop_cyclone_extra = 'obs_seq.tc'
+integer :: max_num_obs = 1000000 ! Largest number of obs in one sequence
+logical :: overwrite_obs_time = .false. ! true to overwrite all observation times
+
+! 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_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_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 /mpas_obs_preproc_nml/ static_mpas_filename, &
+ file_name_input, file_name_output, max_num_obs, &
+ 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, sonde_extra, metar_extra, &
+ acars_extra, land_sfc_extra, marine_sfc_extra, sat_wind_extra, profiler_extra, &
+ trop_cyclone_extra, gpsro_extra, tc_sonde_radii, overwrite_obs_time
+
+!----------------------------------------------------------------------
+! 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, gday, gsec
+
+logical :: file_exist, pre_I_format
+
+type(obs_sequence_type) :: seq_all, seq_rawin, seq_sfc, seq_acars, seq_satwnd, &
+ seq_prof, seq_tc, seq_gpsro, seq_other
+
+type(time_type) :: anal_time
+
+integer :: nCells = -1 ! Total number of cells making up the grid
+integer :: nTimes = -1 ! Total number of times
+integer :: dimid, ncid, VarID
+real(r8), allocatable :: xland(:,:) ! indicator for land (1.0_r8) or ocean (> 1.0_r8)
+
+print*,'Enter target assimilation time (gregorian day, second): '
+read*,gday,gsec
+call set_calendar_type(GREGORIAN)
+anal_time = set_time(gsec, gday)
+
+call static_init_obs_sequence()
+
+call static_init_model()
+call nc_check(nf90_open(trim(static_mpas_filename), NF90_NOWRITE, ncid), &
+ 'mpas_dart_obs_preprocess','open '//trim(static_mpas_filename))
+call nc_check(nf90_inq_dimid(ncid, 'Time', dimid), &
+ 'mpas_dart_obs_preprocess','inq_dimid Time '//trim(static_mpas_filename))
+call nc_check(nf90_inquire_dimension(ncid, dimid, len=nTimes), &
+ 'mpas_dart_obs_preprocess','inquire_dimension nCells '//trim(static_mpas_filename))
+call nc_check(nf90_inq_dimid(ncid, 'nCells', dimid), &
+ 'mpas_dart_obs_preprocess','inq_dimid nCells '//trim(static_mpas_filename))
+call nc_check(nf90_inquire_dimension(ncid, dimid, len=nCells), &
+ 'mpas_dart_obs_preprocess','inquire_dimension nCells '//trim(static_mpas_filename))
+
+! FIXME: how can I count or remove time dimension in [Time|1] x [nCell]?
+allocate(xland(nTimes,nCells))
+call nc_check(nf90_inq_varid(ncid, 'xland', VarID), &
+ 'mpas_dart_obs_preprocess', 'inq_varid xland '//trim(static_mpas_filename))
+call nc_check(nf90_get_var( ncid, VarID, xland), &
+ 'mpas_dart_obs_preprocess', 'get_var xland '//trim(static_mpas_filename))
+
+call find_namelist_in_file("input.nml", "mpas_obs_preproc_nml", iunit)
+read(iunit, nml = mpas_obs_preproc_nml, iostat = io)
+call check_namelist_read(iunit, io, "mpas_obs_preproc_nml")
+
+! 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_obs_seq, seq_rawin)
+call create_new_obs_seq(num_copies, num_qc, max_obs_seq, seq_sfc)
+call create_new_obs_seq(num_copies, num_qc, max_obs_seq, seq_acars)
+call create_new_obs_seq(num_copies, num_qc, max_obs_seq, seq_satwnd)
+call create_new_obs_seq(num_copies, num_qc, max_obs_seq, seq_prof)
+call create_new_obs_seq(num_copies, num_qc, max_obs_seq, seq_gpsro)
+call create_new_obs_seq(num_copies, num_qc, 100, seq_tc)
+call create_new_obs_seq(num_copies, num_qc, max_obs_seq, seq_other)
+
+! read input obs_seq file, divide into platforms
+call read_and_parse_input_seq(file_name_input, nCells, xland(1,:), &
+include_sig_data, obs_pressure_top, obs_height_top, sfc_elevation_check, &
+sfc_elevation_tol, overwrite_ncep_sfc_qc, overwrite_ncep_satwnd_qc, &
+overwrite_obs_time, anal_time, seq_rawin, seq_sfc, seq_acars, seq_satwnd, &
+seq_tc, seq_gpsro, seq_other)
+
+! add supplimental rawinsonde observations from file
+call add_supplimental_obs(sonde_extra, seq_rawin, max_obs_seq, &
+RADIOSONDE_U_WIND_COMPONENT, include_sig_data, &
+obs_pressure_top, obs_height_top, sfc_elevation_check, sfc_elevation_tol, &
+overwrite_obs_time, anal_time)
+
+! add supplimental ACARS observations from file
+call add_supplimental_obs(acars_extra, seq_acars, max_obs_seq, &
+ACARS_U_WIND_COMPONENT, include_sig_data, &
+obs_pressure_top, obs_height_top, sfc_elevation_check, sfc_elevation_tol, &
+overwrite_obs_time, anal_time)
+
+! add supplimental marine observations from file
+call add_supplimental_obs(marine_sfc_extra, seq_sfc, max_obs_seq, &
+MARINE_SFC_U_WIND_COMPONENT, include_sig_data, &
+obs_pressure_top, obs_height_top, sfc_elevation_check, sfc_elevation_tol, &
+overwrite_obs_time, anal_time)
+
+! add supplimental land surface observations from file
+call add_supplimental_obs(land_sfc_extra, seq_sfc, max_obs_seq, &
+LAND_SFC_U_WIND_COMPONENT, include_sig_data, &
+obs_pressure_top, obs_height_top, sfc_elevation_check, sfc_elevation_tol, &
+overwrite_obs_time, anal_time)
+
+! add supplimental metar observations from file
+call add_supplimental_obs(metar_extra, seq_sfc, max_obs_seq, &
+METAR_U_10_METER_WIND, include_sig_data, &
+obs_pressure_top, obs_height_top, sfc_elevation_check, sfc_elevation_tol, &
+overwrite_obs_time, anal_time)
+
+! add supplimental satellite wind observations from file
+call add_supplimental_obs(sat_wind_extra, seq_satwnd, max_obs_seq, &
+SAT_U_WIND_COMPONENT, include_sig_data, &
+obs_pressure_top, obs_height_top, sfc_elevation_check, sfc_elevation_tol, &
+overwrite_obs_time, anal_time)
+
+! add supplimental profiler observations from file
+call add_supplimental_obs(profiler_extra, seq_prof, max_obs_seq, &
+PROFILER_U_WIND_COMPONENT, include_sig_data, &
+obs_pressure_top, obs_height_top, sfc_elevation_check, sfc_elevation_tol, &
+overwrite_obs_time, anal_time)
+
+! add supplimental GPSRO observations from file
+call add_supplimental_obs(gpsro_extra, seq_gpsro, max_obs_seq, &
+GPSRO_REFRACTIVITY, include_sig_data, &
+obs_pressure_top, obs_height_top, sfc_elevation_check, sfc_elevation_tol, &
+overwrite_obs_time, anal_time)
+
+! add supplimental tropical cyclone vortex observations from file
+call add_supplimental_obs(trop_cyclone_extra, seq_tc, max_obs_seq, &
+VORTEX_LAT, include_sig_data, &
+obs_pressure_top, obs_height_top, sfc_elevation_check, sfc_elevation_tol, &
+overwrite_obs_time, anal_time)
+
+! 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, nCells, anal_time, &
+ aircraft_pres_int, obs_pressure_top)
+
+! super-ob satellite wind data
+if ( superob_sat_winds ) call superob_sat_wind_data(seq_satwnd, nCells, anal_time, &
+ sat_wind_pres_int, obs_pressure_top)
+
+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_prof) + &
+ get_num_obs(seq_gpsro) + 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_gpsro, seq_all)
+call destroy_obs_sequence(seq_gpsro)
+
+call build_master_sequence(seq_satwnd, seq_all)
+call destroy_obs_sequence(seq_satwnd)
+
+call build_master_sequence(seq_prof, seq_all)
+call destroy_obs_sequence(seq_prof)
+
+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)
+
+! 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)
+! 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, &
+ siglevel, ptop, htop, &
+ sfcelev, elev_max, overwrite_time, atime)
+
+use types_mod, only : r8
+use time_manager_mod, only : time_type, operator(>=)
+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, set_obs_def, &
+ get_num_copies, get_num_qc, read_obs_seq, copy_obs, &
+ get_first_obs, get_obs_def, get_next_obs, &
+ get_last_obs, insert_obs_in_seq, destroy_obs_sequence
+use obs_def_mod, only : obs_def_type, get_obs_kind, set_obs_def_time, &
+ get_obs_def_location, get_obs_def_time
+use obs_kind_mod, only : RADIOSONDE_U_WIND_COMPONENT, ACARS_U_WIND_COMPONENT, &
+ LAND_SFC_U_WIND_COMPONENT, MARINE_SFC_U_WIND_COMPONENT, &
+ METAR_U_10_METER_WIND, GPSRO_REFRACTIVITY, &
+ SAT_U_WIND_COMPONENT, VORTEX_LAT
+
+implicit none
+
+character(len=129), intent(in) :: filename
+type(time_type), intent(in) :: atime
+type(obs_sequence_type), intent(inout) :: obs_seq
+integer, intent(in) :: max_obs_seq, plat_kind
+logical, intent(in) :: siglevel, sfcelev, overwrite_time
+real(r8), intent(in) :: 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, first_obs
+real(r8) :: llv_loc(3)
+
+
+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_obsi, prev_obso, obs
+type(time_type) :: obs_time, prev_time
+
+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 (METAR_U_10_METER_WIND)
+ write(6,*) 'Adding Supplimental METAR Data'
+ case (SAT_U_WIND_COMPONENT)
+ write(6,*) 'Adding Supplimental Satellite Wind Data'
+ case (VORTEX_LAT)
+ write(6,*) 'Adding Supplimental Tropical Cyclone Data'
+ case (GPSRO_REFRACTIVITY)
+ write(6,*) 'Adding Supplimental GPS RO 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_obsi, get_num_copies(obs_seq), get_num_qc(obs_seq))
+call init_obs(prev_obso, 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)
+
+! find the last observation in the sequence
+if ( get_last_obs(obs_seq, prev_obso) ) then
+
+ first_obs = .false.
+ call get_obs_def(prev_obso, obs_def)
+ prev_time = get_obs_def_time(obs_def)
+
+else
+
+ first_obs = .true.
+
+end if
+
+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)
+ llv_loc = get_location(obs_loc)
+
+ prev_obsi = obs_in
+ call get_next_obs(supp_obs_seq, prev_obsi, obs_in, last_obs)
+ cycle ObsLoop
+
+ ! check if the observation is within vertical bounds of domain
+ if ( (vert_is_pressure(obs_loc) .and. llv_loc(3) < ptop) .or. &
+ (vert_is_height(obs_loc) .and. llv_loc(3) > htop) ) then
+
+ prev_obsi = obs_in
+ call get_next_obs(supp_obs_seq, prev_obsi, 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_obsi = obs_in
+ call get_next_obs(supp_obs_seq, prev_obsi, obs_in, last_obs)
+ cycle ObsLoop
+
+ end if
+
+ ! overwrite the observation time with the analysis time if desired
+ if ( overwrite_time ) then
+
+ call set_obs_def_time(obs_def, atime)
+ call set_obs_def(obs_in, obs_def)
+
+ 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, llv_loc)
+ case (LAND_SFC_U_WIND_COMPONENT)
+ pass_checks = surface_obs_check(sfcelev, elev_max, llv_loc)
+ case (METAR_U_10_METER_WIND)
+ pass_checks = surface_obs_check(sfcelev, elev_max, llv_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 get_obs_def(obs, obs_def)
+ obs_time = get_obs_def_time(obs_def)
+
+ if (obs_time >= prev_time .and. (.not. first_obs)) then ! same time or later than previous obs
+ call insert_obs_in_seq(obs_seq, obs, prev_obso)
+ else ! earlier, search from start of seq
+ call insert_obs_in_seq(obs_seq, obs)
+ end if
+
+ first_obs = .false.
+ prev_obso = obs
+ prev_time = obs_time
+
+ end if
+
+ prev_obsi = obs_in
+ call get_next_obs(supp_obs_seq, prev_obsi, 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 time_manager_mod, only : time_type, operator(>=)
+use obs_sequence_mod, only : obs_type, obs_sequence_type, init_obs, &
+ get_first_obs, copy_obs, insert_obs_in_seq, &
+ get_next_obs, get_obs_def, get_num_copies, &
+ get_num_qc, get_obs_def
+use obs_def_mod, only : obs_def_type, get_obs_def_time
+
+implicit none
+
+type(obs_sequence_type), intent(in) :: seq_type
+type(obs_sequence_type), intent(inout) :: seq_all
+
+logical :: last_obs, first_obs
+type(obs_def_type) :: obs_def
+type(obs_type) :: obs_in, obs, prev_obsi, prev_obsa
+type(time_type) :: obs_time, prev_time
+
+last_obs = .false. ; first_obs = .true.
+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_obsi, get_num_copies(seq_type), get_num_qc(seq_type))
+call init_obs(prev_obsa, 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 get_obs_def(obs, obs_def)
+ obs_time = get_obs_def_time(obs_def)
+
+ if (obs_time >= prev_time .and. (.not. first_obs)) then ! same time or later than previous obs
+ call insert_obs_in_seq(seq_all, obs, prev_obsa)
+ else ! earlier, search from start of seq
+ call insert_obs_in_seq(seq_all, obs)
+ end if
+
+ first_obs = .false.
+ prev_obsi = obs_in
+ prev_obsa = obs
+ prev_time = obs_time
+
+ call get_next_obs(seq_type, prev_obsi, 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
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! 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) :: llv_loc(3), xmod(1), hsfc
+
+rawinsonde_obs_check = .true.
+llv_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(llv_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 - llv_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
+! 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, nCell, xland, siglevel, ptop, &
+ htop, sfcelev, elev_max, new_sfc_qc, &
+ new_satwnd_qc, overwrite_time, atime, &
+ rawin_seq, sfc_seq, acars_seq, satwnd_seq, &
+ tc_seq, gpsro_seq, other_seq)
+
+use types_mod, only : r8
+use utilities_mod, only : nc_check
+use time_manager_mod, only : time_type
+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, set_obs_def
+use obs_def_mod, only : obs_def_type, get_obs_kind, get_obs_def_location, &
+ set_obs_def_time
+use obs_kind_mod, only : RADIOSONDE_U_WIND_COMPONENT, RADIOSONDE_V_WIND_COMPONENT, &
+ RADIOSONDE_SURFACE_ALTIMETER, RADIOSONDE_TEMPERATURE, &
+ RADIOSONDE_SPECIFIC_HUMIDITY, RADIOSONDE_DEWPOINT, &
+ RADIOSONDE_RELATIVE_HUMIDITY, GPSRO_REFRACTIVITY, &
+ AIRCRAFT_U_WIND_COMPONENT, AIRCRAFT_V_WIND_COMPONENT, &
+ AIRCRAFT_TEMPERATURE, AIRCRAFT_SPECIFIC_HUMIDITY, &
+ ACARS_DEWPOINT, ACARS_RELATIVE_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, &
+ MARINE_SFC_RELATIVE_HUMIDITY, MARINE_SFC_DEWPOINT, &
+ LAND_SFC_U_WIND_COMPONENT, LAND_SFC_V_WIND_COMPONENT, &
+ LAND_SFC_TEMPERATURE, LAND_SFC_SPECIFIC_HUMIDITY, &
+ LAND_SFC_RELATIVE_HUMIDITY, LAND_SFC_DEWPOINT, &
+ METAR_U_10_METER_WIND, METAR_V_10_METER_WIND, &
+ METAR_TEMPERATURE_2_METER, METAR_SPECIFIC_HUMIDITY_2_METER, &
+ METAR_DEWPOINT_2_METER, METAR_RELATIVE_HUMIDITY_2_METER, &
+ METAR_ALTIMETER, 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 : find_closest_cell_center
+use netcdf
+
+implicit none
+
+real(r8), parameter :: satwnd_qc_ok = 15.0_r8
+real(r8), parameter :: sfc_qc_ok1 = 9.0_r8
+real(r8), parameter :: sfc_qc_ok2 = 15.0_r8
+real(r8), parameter :: new_qc_value = 2.0_r8
+
+character(len=129), intent(in) :: filename
+integer, intent(in) :: nCell
+real(r8), intent(in) :: xland(nCell)
+real(r8), intent(in) :: ptop, htop, elev_max
+logical, intent(in) :: siglevel, sfcelev, new_sfc_qc, &
+ new_satwnd_qc, overwrite_time
+type(time_type), intent(in) :: atime
+type(obs_sequence_type), intent(inout) :: rawin_seq, sfc_seq, acars_seq, &
+ satwnd_seq, tc_seq, gpsro_seq, other_seq
+
+character(len=129) :: qcmeta
+integer :: fid, var_id, okind, dom_id, cellid
+logical :: file_exist, last_obs, input_ncep_qc, rawinsonde_obs_check, &
+ surface_obs_check, aircraft_obs_check, sat_wind_obs_check
+real(r8), allocatable :: qc(:)
+real(r8) :: llv_loc(3)
+
+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)))
+
+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)
+ llv_loc = get_location(obs_loc)
+ cellid = find_closest_cell_center(llv_loc(2), llv_loc(1))
+
+ prev_obs = obs_in
+ call get_next_obs(seq, prev_obs, obs_in, last_obs)
+ cycle InputObsLoop
+
+ ! check vertical location
+ if ( (vert_is_pressure(obs_loc) .and. llv_loc(3) < ptop) .or. &
+ (vert_is_height(obs_loc) .and. llv_loc(3) > htop) ) then
+
+ prev_obs = obs_in
+ call get_next_obs(seq, prev_obs, obs_in, last_obs)
+ cycle InputObsLoop
+
+ end if
+
+ ! overwrite the observation time with the analysis time if desired
+ if ( overwrite_time ) then
+
+ call set_obs_def_time(obs_def, atime)
+ call set_obs_def(obs_in, obs_def)
+
+ end if
+
+ ! perform platform-specific checks
+ select case (okind)
+
+ case ( RADIOSONDE_U_WIND_COMPONENT, RADIOSONDE_V_WIND_COMPONENT, &
+ RADIOSONDE_TEMPERATURE, RADIOSONDE_SPECIFIC_HUMIDITY, &
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list