[Dart-dev] [4977] DART/trunk: This is a much more mature version of the obs-space forecast
nancy at ucar.edu
nancy at ucar.edu
Thu Jun 9 15:15:38 MDT 2011
Revision: 4977
Author: thoar
Date: 2011-06-09 15:15:38 -0600 (Thu, 09 Jun 2011)
Log Message:
-----------
This is a much more mature version of the obs-space forecast
evaluation tools and documentation. The R diagnostic scripts
are not ready for distribution - yet.
Modified Paths:
--------------
DART/trunk/models/wrf/work/input.nml
DART/trunk/obs_sequence/obs_seq_coverage.f90
DART/trunk/obs_sequence/obs_seq_coverage.nml
Added Paths:
-----------
DART/trunk/doc/images/forecasting_diagram.png
DART/trunk/doc/images/obs_seq_coverage_diagram.png
DART/trunk/doc/images/obs_seq_verify_diagram.png
DART/trunk/doc/images/simple_forecast.png
DART/trunk/doc/images/verification_48hrX6hr.png
DART/trunk/doc/images/verification_time_icon.png
DART/trunk/models/wrf/work/mkmf_obs_seq_verify
DART/trunk/models/wrf/work/path_names_obs_seq_verify
DART/trunk/obs_sequence/obs_seq_coverage.html
DART/trunk/obs_sequence/obs_seq_verify.f90
DART/trunk/obs_sequence/obs_seq_verify.html
DART/trunk/obs_sequence/obs_seq_verify.nml
-------------- next part --------------
Added: DART/trunk/doc/images/forecasting_diagram.png
===================================================================
(Binary files differ)
Property changes on: DART/trunk/doc/images/forecasting_diagram.png
___________________________________________________________________
Added: svn:mime-type
+ image/png
Added: DART/trunk/doc/images/obs_seq_coverage_diagram.png
===================================================================
(Binary files differ)
Property changes on: DART/trunk/doc/images/obs_seq_coverage_diagram.png
___________________________________________________________________
Added: svn:mime-type
+ image/png
Added: DART/trunk/doc/images/obs_seq_verify_diagram.png
===================================================================
(Binary files differ)
Property changes on: DART/trunk/doc/images/obs_seq_verify_diagram.png
___________________________________________________________________
Added: svn:mime-type
+ image/png
Added: DART/trunk/doc/images/simple_forecast.png
===================================================================
(Binary files differ)
Property changes on: DART/trunk/doc/images/simple_forecast.png
___________________________________________________________________
Added: svn:mime-type
+ image/png
Added: DART/trunk/doc/images/verification_48hrX6hr.png
===================================================================
(Binary files differ)
Property changes on: DART/trunk/doc/images/verification_48hrX6hr.png
___________________________________________________________________
Added: svn:mime-type
+ image/png
Added: DART/trunk/doc/images/verification_time_icon.png
===================================================================
(Binary files differ)
Property changes on: DART/trunk/doc/images/verification_time_icon.png
___________________________________________________________________
Added: svn:mime-type
+ image/png
Modified: DART/trunk/models/wrf/work/input.nml
===================================================================
--- DART/trunk/models/wrf/work/input.nml 2011-06-09 20:02:27 UTC (rev 4976)
+++ DART/trunk/models/wrf/work/input.nml 2011-06-09 21:15:38 UTC (rev 4977)
@@ -19,7 +19,7 @@
output_forward_op_errors = .false.,
print_every_nth_obs = -1,
silence = .false.,
- /
+ /
&filter_nml
async = 2,
@@ -154,7 +154,7 @@
&assim_model_nml
write_binary_restart_files = .true.,
netCDF_large_file_support = .false.,
- /
+ /
# Notes for model_nml:
# (1) vert_localization_coord must be one of:
@@ -333,8 +333,18 @@
selections_is_obs_seq = .false.,
print_only = .false.,
calendar = 'gregorian',
-/
+ /
+&obs_seq_verify_nml
+ obs_sequence_list = '',
+ obs_sequence_name = 'obs_seq.forecast',
+ station_template = 'obsdef_mask.nc',
+ netcdf_out = 'forecast.nc',
+ obtype_string = 'METAR_U_10_METER_WIND',
+ verbose = .true.,
+ debug = .false.,
+ /
+
&restart_file_tool_nml
input_file_name = "filter_restart",
output_file_name = "filter_updated_restart",
Added: DART/trunk/models/wrf/work/mkmf_obs_seq_verify
===================================================================
--- DART/trunk/models/wrf/work/mkmf_obs_seq_verify (rev 0)
+++ DART/trunk/models/wrf/work/mkmf_obs_seq_verify 2011-06-09 21:15:38 UTC (rev 4977)
@@ -0,0 +1,18 @@
+#!/bin/csh
+#
+# DART software - Copyright 2004 - 2011 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$
+
+../../../mkmf/mkmf -p obs_seq_verify -t ../../../mkmf/mkmf.template -c"-Duse_netCDF" \
+ -a "../../.." $* path_names_obs_seq_verify
+
+exit $status
+
+# <next few lines under version control, do not edit>
+# $URL$
+# $Revision$
+# $Date$
+
Property changes on: DART/trunk/models/wrf/work/mkmf_obs_seq_verify
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Added: DART/trunk/models/wrf/work/path_names_obs_seq_verify
===================================================================
--- DART/trunk/models/wrf/work/path_names_obs_seq_verify (rev 0)
+++ DART/trunk/models/wrf/work/path_names_obs_seq_verify 2011-06-09 21:15:38 UTC (rev 4977)
@@ -0,0 +1,15 @@
+assim_model/assim_model_mod.f90
+common/types_mod.f90
+location/threed_sphere/location_mod.f90
+models/wrf/model_mod.f90
+models/wrf/module_map_utils.f90
+mpi_utilities/null_mpi_utilities_mod.f90
+obs_def/obs_def_mod.f90
+obs_kind/obs_kind_mod.f90
+obs_sequence/obs_sequence_mod.f90
+obs_sequence/obs_seq_verify.f90
+random_nr/random_nr_mod.f90
+random_seq/random_seq_mod.f90
+sort/sort_mod.f90
+time_manager/time_manager_mod.f90
+utilities/utilities_mod.f90
Modified: DART/trunk/obs_sequence/obs_seq_coverage.f90
===================================================================
--- DART/trunk/obs_sequence/obs_seq_coverage.f90 2011-06-09 20:02:27 UTC (rev 4976)
+++ DART/trunk/obs_sequence/obs_seq_coverage.f90 2011-06-09 21:15:38 UTC (rev 4977)
@@ -18,9 +18,9 @@
! The observation sequence file only contains lat/lon/level/which_vert,
! so this is all we have to work with.
!
-! From: Soyoung Ha <syha at ucar.edu>
+! From: Soyoung Ha
! Date: November 5, 2010 12:00:22 PM MDT
-! To: Tim Hoar <thoar at ucar.edu>, Nancy Collins <nancy at ucar.edu>
+! To: Tim Hoar, Nancy Collins
! Subject: obs_seq.out on bluefire
!
! Hi Tim and Nancy,
@@ -39,13 +39,14 @@
!-----------------------------------------------------------------------
-use types_mod, only : r4, r8, digits12, MISSING_R8, MISSING_R4
+use types_mod, only : r4, r8, digits12, MISSING_R8, MISSING_R4, &
+ metadatalength, obstypelength
use obs_sequence_mod, only : read_obs_seq, obs_type, obs_sequence_type, get_first_obs, &
get_obs_def, get_copy_meta_data, &
get_next_obs, init_obs, init_obs_sequence, &
- assignment(=), get_num_copies, static_init_obs_sequence, &
- get_qc, destroy_obs_sequence, read_obs_seq_header, &
- destroy_obs, get_qc_meta_data
+ assignment(=), get_num_copies, get_num_qc, get_qc, &
+ static_init_obs_sequence, destroy_obs_sequence, destroy_obs, &
+ read_obs_seq_header, get_qc_meta_data
use obs_def_mod, only : obs_def_type, get_obs_def_time, get_obs_kind, write_obs_def, &
get_obs_def_location, set_obs_def_time, &
set_obs_def_location, set_obs_def_kind, set_obs_def_error_variance
@@ -56,9 +57,12 @@
set_location, is_location_in_region, query_location, &
nc_write_location_atts, nc_get_location_varids, &
nc_write_location, LocationDims
-use time_manager_mod, only : time_type, set_date, set_time, get_time, print_time, &
- set_time_missing, operator(>), operator(<), operator(==), &
- operator(<=), operator(-), operator(+), operator(/=)
+use time_manager_mod, only : time_type, set_date, set_time, get_time, &
+ set_calendar_type, get_calendar_string, &
+ print_time, print_date, &
+ operator(+), operator(-), operator(<), operator(>), &
+ operator(==), operator(/=), operator(<=), operator(/), &
+ operator(>=), operator(*)
use utilities_mod, only : get_unit, close_file, register_module, &
file_exist, error_handler, E_ERR, E_WARN, E_MSG, &
initialize_utilities, nmlfileunit, timestamp, &
@@ -96,9 +100,7 @@
integer :: station_id ! the index (into stations) of an existing location
integer :: timeindex ! the index (into the time array of a station)
integer :: num_output ! total number of desired locations and times found
-integer :: num_max ! most number of desired times found at any location
-integer, parameter :: STRINGLENGTH = 32
integer, parameter :: MAX_OBS_INPUT_TYPES = 500 ! lazy, just going big
!---------------------------------------------------------------------
@@ -109,7 +111,6 @@
type(obs_type) :: obs1, obs2
type(obs_def_type) :: obs_def
type(location_type) :: obs_loc, minl, maxl
-real(r8), dimension(LocationDims) :: locarray
character(len = 129) :: obs_seq_in_file_name
character(len = 129), allocatable, dimension(:) :: obs_seq_filenames
@@ -125,53 +126,72 @@
! Namelist with (some scalar) default values
!-----------------------------------------------------------------------
+character(len = 129) :: obs_sequence_list = 'obs_coverage_list.txt'
+character(len = 129) :: obs_sequence_name = ''
+character(len = obstypelength) :: obs_of_interest(MAX_OBS_INPUT_TYPES) = ''
character(len = 129) :: textfile_out = 'obsdef_mask.txt'
character(len = 129) :: netcdf_out = 'obsdef_mask.nc'
-character(len = 129) :: obs_sequence_name = 'obs_seq.final'
-character(len = 129) :: obs_sequence_list = ''
-character(len = STRINGLENGTH) :: obs_of_interest(MAX_OBS_INPUT_TYPES) = ''
+character(len = 129) :: calendar = 'Gregorian'
-real(r8) :: lonlim1= MISSING_R8, lonlim2= MISSING_R8
-real(r8) :: latlim1= MISSING_R8, latlim2= MISSING_R8
-integer :: nTmin = 0 ! minimum number of times required
-integer :: nTmax = 0 ! maximum number of times required
+integer, dimension(6) :: first_analysis = (/ 2003, 1, 1, 0, 0, 0 /)
+integer, dimension(6) :: last_analysis = (/ 2003, 1, 2, 0, 0, 0 /)
+integer :: forecast_length_days = 1
+integer :: forecast_length_seconds = 0
+integer :: verification_interval_seconds = 21600 ! 6 hours
+real(r8) :: temporal_coverage_percent = 100.0 ! all times required
+real(r8) :: lonlim1 = MISSING_R8
+real(r8) :: lonlim2 = MISSING_R8
+real(r8) :: latlim1 = MISSING_R8
+real(r8) :: latlim2 = MISSING_R8
+logical :: verbose = .false.
+logical :: debug = .false. ! undocumented ... on purpose
-logical :: debug = .false. ! undocumented ... on purpose
-logical :: verbose = .false.
+namelist /obs_seq_coverage_nml/ obs_sequence_list, obs_sequence_name, &
+ obs_of_interest, textfile_out, netcdf_out, calendar, &
+ first_analysis, last_analysis, forecast_length_days, &
+ forecast_length_seconds, verification_interval_seconds, &
+ temporal_coverage_percent, lonlim1, lonlim2, latlim1, latlim2, &
+ verbose, debug
-namelist /obs_seq_coverage_nml/ obs_sequence_name, obs_sequence_list, &
- lonlim1, lonlim2, latlim1, latlim2, &
- nTmin, nTmax, obs_of_interest, &
- verbose, debug, textfile_out, netcdf_out
-
!-----------------------------------------------------------------------
! Quantities of interest
!-----------------------------------------------------------------------
-integer, parameter :: Ncopies = 1
-integer :: allNcopies
-character(len=STRINGLENGTH), dimension(Ncopies) :: copy_names = &
- (/ 'observation error variance' /)
+integer :: qc_index ! copy index of the original qc value
+integer :: dart_qc_index ! copy index of the DART qc value
-character(len=STRINGLENGTH), allocatable, dimension(:) :: module_obs_copy_names
-character(len=STRINGLENGTH), allocatable, dimension(:) :: obs_copy_names
-character(len=STRINGLENGTH), allocatable, dimension(:) :: module_qc_copy_names
-character(len=STRINGLENGTH), allocatable, dimension(:) :: qc_copy_names
+character(len=metadatalength), allocatable, dimension(:) :: module_obs_copy_names
+character(len=metadatalength), allocatable, dimension(:) :: obs_copy_names
+character(len=metadatalength), allocatable, dimension(:) :: module_qc_copy_names
+character(len=metadatalength), allocatable, dimension(:) :: qc_copy_names
-real(r8), allocatable, dimension(:) :: qc
integer, dimension(max_obs_kinds) :: obs_kinds_inds = 0
+real(r8), allocatable, dimension(:) :: qc_values
+type(time_type), allocatable, dimension(:) :: all_verif_times
+type(time_type), allocatable, dimension(:,:) :: verification_times
+real(digits12), allocatable, dimension(:,:) :: experiment_Tr8
+type(time_type) :: verification_stride, half_stride
+integer :: num_analyses ! # of fcsts from first_analysis to last_analysis
+integer :: num_verify_per_fcst
+integer :: num_verification_times ! number of verification times - total
+integer :: nT_minimum ! will settle for this many verif times - total
+
!-----------------------------------------------------------------------
! General purpose variables
!-----------------------------------------------------------------------
integer :: ifile, nread, ngood
-integer :: i, io, ncunit
+integer :: i, j, io, ncunit
-type(time_type) :: obs_time
+type(time_type) :: obs_time, no_time, last_possible_time
character(len = 129) :: ncName, string1, string2, string3
+! ~# of degrees for 1/2 meter at Earth equator
+! 360 deg-earth/(40000 km-earth * 1000m-km)
+real(r8), parameter :: HALF_METER = 180.0_r8 / (40000.0_r8 * 1000.0_r8)
+
!=======================================================================
! Get the party started
!=======================================================================
@@ -183,14 +203,8 @@
call init_obs(obs1, 0, 0)
call init_obs(obs2, 0, 0)
call init_obs_sequence(seq,0,0,0)
+no_time = set_time(0,0)
-! Allocate a hunk of stations. If we fill this up, we will
-! have to create temporary storage, copy, deallocate, reallocate ...
-
-num_stations = 0
-max_stations = 4000
-call initialize_stations(max_stations, stations)
-
!----------------------------------------------------------------------
! Read the namelist
!----------------------------------------------------------------------
@@ -204,14 +218,12 @@
if (do_nml_term()) write( * , nml=obs_seq_coverage_nml)
! Check the user input for sanity
-if (nTmin > nTmax) then
- write(string1,*)'namelist: nTmin (',nTmin,') must be <= nTmax (',nTmax,')'
+if (temporal_coverage_percent < 100.0_r8) then
+ write(string1,*)'namelist: temporal_coverage_percent (',temporal_coverage_percent,&
+ ') must be == 100.0 for now.)'
call error_handler(E_ERR, 'obs_seq_coverage', string1, source, revision, revdate)
endif
-if (nTmin < 0) then
- write(string1,*)'nTmin must be > 0, was read as ',nTmin
- call error_handler(E_ERR, 'obs_seq_coverage', string1, source, revision, revdate)
-endif
+
if ((obs_sequence_name /= '') .and. (obs_sequence_list /= '')) then
write(string1,*)'specify "obs_sequence_name" or "obs_sequence_list"'
write(string2,*)'set other to an empty string ... i.e. ""'
@@ -219,9 +231,13 @@
revdate, text2=string2)
endif
-!----------------------------------------------------------------------
+call set_calendar_type(calendar)
+call get_calendar_string(calendar)
+
+minl = set_location( (/ lonlim1, latlim1, 0.0_r8, 1.0_r8 /)) ! vertical unimportant
+maxl = set_location( (/ lonlim2, latlim2, 0.0_r8, 1.0_r8 /)) ! vertical unimportant
+
! Determine if the desired observation types exist
-!----------------------------------------------------------------------
TypeLoop : do i = 1,MAX_OBS_INPUT_TYPES
@@ -243,11 +259,34 @@
enddo TypeLoop
+! Set the verification time array (global storage)
+
+call set_required_times(first_analysis, last_analysis, &
+ forecast_length_days, forecast_length_seconds, &
+ verification_interval_seconds, temporal_coverage_percent)
+
+if (verbose) then
+ write(*,*) ! whitespace
+ write(*,*)'At least',nT_minimum,' observations times are required at:'
+ do i=1,num_verification_times
+ write(string1,*)'verification # ',i,' at '
+ call print_date(all_verif_times(i),trim(string1))
+ enddo
+ write(*,*) ! whitespace
+endif
+
+last_possible_time = all_verif_times(num_verification_times) + half_stride
+
+! Allocate a hunk of stations. If we fill this up, we will
+! have to create temporary storage, copy, deallocate, reallocate ...
+
+num_stations = 0
+max_stations = 4000
+call initialize_stations(max_stations, stations)
+
!====================================================================
!====================================================================
-minl = set_location( (/ lonlim1, latlim1, 0.0_r8, 1.0_r8 /)) ! vertical unimportant
-maxl = set_location( (/ lonlim2, latlim2, 0.0_r8, 1.0_r8 /)) ! vertical unimportant
!----------------------------------------------------------------------
! Prepare the variables
@@ -266,7 +305,7 @@
call destroy_obs(obs2)
call destroy_obs_sequence(seq)
- if (allocated(qc)) deallocate(qc)
+ if (allocated(qc_values)) deallocate(qc_values)
if (allocated(qc_copy_names)) deallocate(qc_copy_names)
if (allocated(obs_copy_names)) deallocate(obs_copy_names)
@@ -286,6 +325,7 @@
write(string1,*)trim(obs_seq_in_file_name),&
' does not exist. Finishing up.'
call error_handler(E_MSG,'obs_seq_coverage',string1,source,revision,revdate)
+ write(*,*) ! whitespace
exit ObsFileLoop
endif
@@ -307,14 +347,12 @@
! I am taking the observational error variance and making it one of the copies
- allNcopies = num_copies + Ncopies
-
if ((num_qc <= 0) .or. (num_copies <=0)) then
write(string1,*)'need at least 1 qc and 1 observation copy'
call error_handler(E_ERR,'obs_seq_coverage',string1,source,revision,revdate)
endif
- allocate( obs_copy_names(allNcopies), qc_copy_names(num_qc), qc(num_qc))
+ allocate( obs_copy_names(num_copies), qc_copy_names(num_qc), qc_values(num_qc))
if ( debug ) then
write(*,*)
@@ -341,37 +379,37 @@
call read_obs_seq(obs_seq_in_file_name, 0, 0, 0, seq)
do i=1, num_copies
- string1 = trim(get_copy_meta_data(seq,i))//' '
- obs_copy_names(i) = string1(1:STRINGLENGTH)
+ string1 = trim(get_copy_meta_data(seq,i))
+ obs_copy_names(i) = adjustl(string1)
enddo
- do i=1, Ncopies
- obs_copy_names(num_copies+i) = trim(copy_names(i))
- enddo
do i=1, num_qc
- string1 = trim(get_qc_meta_data(seq,i))//' '
- qc_copy_names(i) = string1(1:STRINGLENGTH)
+ string1 = trim(get_qc_meta_data(seq,i))
+ qc_copy_names(i) = adjustl(string1)
enddo
+ call find_our_copies(seq, qc_index, dart_qc_index)
+
if ( ifile == 1 ) then ! record the metadata for comparison
- allocate(module_obs_copy_names(allNcopies), &
+ allocate(module_obs_copy_names(num_copies), &
module_qc_copy_names(num_qc) )
do i=1, num_copies
- string1 = trim(get_copy_meta_data(seq,i))//' '
- module_obs_copy_names(i) = string1(1:STRINGLENGTH)
+ module_obs_copy_names(i) = obs_copy_names(i)
enddo
- do i=1, Ncopies
- module_obs_copy_names(num_copies+i) = trim(copy_names(i))
- enddo
do i=1, num_qc
- string1 = trim(get_qc_meta_data(seq,i))//' '
- module_qc_copy_names(i) = string1(1:STRINGLENGTH)
+ module_qc_copy_names(i) = qc_copy_names(i)
enddo
else ! Compare all subsequent files' metadata to the first one
- do i = 1,allNcopies
+ if (num_copies /= size(module_obs_copy_names)) then
+ write(string1,'(''num_copies '',i3,'' does not match '',i3)') &
+ num_copies, size(module_obs_copy_names)
+ call error_handler(E_ERR,'obs_seq_coverage',string1,source,revision,revdate)
+ endif
+
+ do i = 1,num_copies
if (trim(obs_copy_names(i)) /= trim(module_obs_copy_names(i))) then
write(string1,'(''obs copy '',i3,'' from '',a)') i,trim(obs_seq_in_file_name)
call error_handler(E_MSG,'obs_seq_coverage',string1,source,revision,revdate)
@@ -414,43 +452,51 @@
ObservationLoop : do nread = 1,num_obs
!--------------------------------------------------------------------
- if ( verbose .and. (mod(nread,1000) == 0) ) &
+ if ( verbose .and. (mod(nread,10000) == 0) ) &
write(*,*)'Processing obs ',nread,' of ',num_obs
- call get_obs_def(obs1, obs_def)
- call get_qc( obs1, qc)
-
- flavor = get_obs_kind(obs_def)
- obs_time = get_obs_def_time(obs_def)
+ call get_obs_def(obs1, obs_def)
+ flavor = get_obs_kind( obs_def)
+ obs_time = get_obs_def_time( obs_def)
obs_loc = get_obs_def_location(obs_def)
- if (verbose .and. (nread == 1)) call print_time(obs_time,'First observation time')
+ call get_qc( obs1, qc_values)
+ if (verbose .and. (nread == 1)) then
+ call print_time(obs_time,'First observation time')
+ call print_date(obs_time,'First observation date')
+ write(*,*) ! whitespace
+ endif
+
+ if (obs_time > last_possible_time) exit ObsFileLoop
+
!-----------------------------------------------------------------
- ! determine if obs is a new location or time at an existing loc
- ! first : reject if not in desired region or not a type we want
+ ! * reject if not in desired region
+ ! * reject if not a type we want [tracked in obs_kinds_inds(:)]
+ ! * reject if dart_qc exists and is a 4 ...
!-----------------------------------------------------------------
- if ( is_location_in_region(obs_loc,minl,maxl) .and. &
- (obs_kinds_inds(flavor) > 0) ) then
+ if ( .not. is_location_in_region(obs_loc,minl,maxl) ) goto 100
+ if ( obs_kinds_inds(flavor) <= 0) goto 100
+ if ( dart_qc_index > 0 ) then
+ if (qc_values(dart_qc_index) == 4) goto 100
+ endif
- ngood = ngood + 1
- station_id = find_station_location(flavor, obs_loc, stations)
+ ngood = ngood + 1
- if ( station_id < 1 ) then
- station_id = add_new_station(flavor, obs_loc, stations)
- ! else
- ! if (verbose) write(*,*)'obs(',nread,') matches station ',station_id
- endif
+ ! determine if obs is a new location or time at an existing loc
- if ( is_time_new( obs_time, station_id, stations, timeindex) ) &
- call update_time( obs_time, station_id, stations, timeindex)
+ station_id = find_station_location(flavor, obs_loc, stations)
- else
- locarray = get_location(obs_loc)
- if (debug) write(*,*)'obs(',nread,') type ',flavor,'is not wanted',locarray
+ if ( station_id < 1 ) then
+ station_id = add_new_station(flavor, obs_loc, stations)
endif
+ if ( is_time_wanted( obs_time, station_id, stations, timeindex) ) &
+ call update_time( obs_time, station_id, stations, timeindex)
+
+ 100 continue
+
call get_next_obs(seq, obs1, obs2, last_ob_flag)
if (.not. last_ob_flag) obs1 = obs2
@@ -467,20 +513,25 @@
allocate(DesiredStations(num_stations))
DesiredStations = .FALSE.
num_output = 0
-num_max = stations(i)%ntimes
-TimeLoop : do i = 1,num_stations
+do i = 1,num_stations
- if (stations(i)%ntimes > num_max) num_max = stations(i)%ntimes
+ stations(i)%ntimes = 0
- if ( (stations(i)%ntimes >= nTmin) .and. &
- (stations(i)%ntimes <= nTmax) ) then
+ do j = 1,num_verification_times
+ if (stations(i)%times(j) /= no_time) &
+ stations(i)%ntimes = stations(i)%ntimes + 1
+ enddo
+
+ if (stations(i)%ntimes >= nT_minimum) then
DesiredStations(i) = .TRUE.
- num_output = num_output + stations(i)%ntimes
+ num_output = num_output + 1
endif
-enddo TimeLoop
+enddo
+if (verbose) write(*,*)'There were ',num_output,' stations matching the input criterion.'
+
! Output a netCDF file of 'all' observations locations and times.
! Used to explore what is available.
@@ -492,15 +543,14 @@
! if no stations are selected, do something.
if (num_output < 1) then
- write(string1,*)'No location had at least ',nTmin,' reporting times.'
- write(string2,*)'Most was ',num_max
+ write(string1,*)'No location had at least ',nT_minimum,' reporting times.'
call error_handler(E_ERR, 'obs_seq_coverage', string1, source, revision, &
revdate, text2=string2)
endif
! Output the file of desired observation locations and times.
! Used to subset the observation sequence files.
-call print_summary
+call write_obsdefs
!-----------------------------------------------------------------------
! Really, really, done.
@@ -511,7 +561,7 @@
call destroy_obs_sequence(seq)
call destroy_stations(stations)
-if (allocated(qc)) deallocate(qc)
+if (allocated(qc_values)) deallocate(qc_values)
if (allocated(qc_copy_names)) deallocate(qc_copy_names)
if (allocated(obs_copy_names)) deallocate(obs_copy_names)
if (allocated(module_obs_copy_names)) deallocate(module_obs_copy_names)
@@ -525,8 +575,14 @@
CONTAINS
!======================================================================
+
function find_station_location(ObsType, ObsLocation, stationlist) result(station_id)
! Simply try to find a matching lat/lon for an observation type
+! The lons/lats get yanked around "a lot" - being converted from ASCII radians
+! to r8 degrees to r8 radians to r8 degrees and then checked for "equality".
+! So - we're actually just checking to see if the lat/lon is within something
+! like 500 cm either direction. Seems like a reasonable definition of 'match'.
+!
integer, intent(in) :: ObsType
type(location_type), intent(in) :: ObsLocation
type(station), dimension(:), intent(inout) :: stationlist
@@ -546,13 +602,10 @@
londiff = abs(obslocarray(1) - stnlocarray(1))
latdiff = abs(obslocarray(2) - stnlocarray(2))
- if ( (londiff <= epsilon(londiff)) .and. &
- (latdiff <= epsilon(latdiff)) .and. &
+ if ( (londiff <= HALF_METER) .and. &
+ (latdiff <= HALF_METER) .and. &
(ObsType == stationlist(i)%obs_type) ) then
- ! if ( (ObsLocation == stationlist(i)%location) .and. &
- ! (ObsType == stationlist(i)%obs_type) ) then
-
station_id = i
exit FindLoop
endif
@@ -562,7 +615,9 @@
end function find_station_location
+!============================================================================
+
function add_new_station(ObsType, ObsLocation, stationlist) result(station_id)
! Ugh ... if a new location is found, add it. If the stationlist does not have
@@ -655,105 +710,68 @@
end function add_new_station
+!============================================================================
-function is_time_new(ObsTime, stationid, stationlist, timeindex)
-! Determine if the observation time is not already in the registry
-! of the times for the particular station.
+function is_time_wanted(ObsTime, stationid, stationlist, timeindex)
+! The station has a list of the observation times closest to the
+! verification times. Determine if the observation time is closer to
+! the verification time than what we already have.
+
type(time_type), intent(in) :: ObsTime
integer, intent(in) :: stationid
type(station), dimension(:), intent(in) :: stationlist
integer, intent(out) :: timeindex
-logical :: is_time_new
+logical :: is_time_wanted
-type(time_type) :: stnhour, obhour, stndelta, obdelta
+type(time_type) :: stndelta, obdelta
integer :: i
-logical :: have_this_hour
-have_this_hour = .FALSE.
-timeindex = 0
+timeindex = 0
+is_time_wanted = .FALSE.
-if ( stationlist(stationid)%ntimes == 0 ) then
- is_time_new = .TRUE.
- timeindex = 1
- return
-endif
+! the time_minus function always returns a positive difference
-TimeLoop : do i = 1,stationlist(stationid)%ntimes
+TimeLoop : do i = 1,num_verification_times
- stnhour = nearest_hour(stationlist(stationid)%times(i))
- obhour = nearest_hour(ObsTime)
+ obdelta = ObsTime - all_verif_times(i)
- ! Make sure we only compare observations to the same hour
- if (stnhour /= obhour) cycle TimeLoop
+ ! If observation is not within half a verification step,
+ ! try the next one.
+ if (obdelta >= half_stride) cycle TimeLoop
- have_this_hour = .TRUE.
+ stndelta = stationlist(stationid)%times(i) - all_verif_times(i)
- ! the time_minus function always returns a positive difference
- stndelta = stationlist(stationid)%times(i) - stnhour
- obdelta = ObsTime - obhour
-
- ! Check to see if the increment is smaller
- ! if it is, then we know which one to overwrite
+ ! Check to see if the observation is closer to the verification time
+ ! than the one we have.
if (obdelta < stndelta) then
if (debug) call print_time(stationlist(stationid)%times(i),'replacing ')
if (debug) call print_time(ObsTime,'with this observation time')
timeindex = i
+ is_time_wanted = .TRUE.
exit TimeLoop
endif
enddo TimeLoop
-if (have_this_hour .and. (timeindex == 0) ) then
- ! the time we already have is closer than the candidate
- is_time_new = .FALSE.
-elseif ( have_this_hour ) then ! candidate is closer
- is_time_new = .TRUE.
-else ! must be a new observation hour
- is_time_new = .TRUE.
- timeindex = stationlist(stationid)%ntimes + 1
-endif
+end function is_time_wanted
-end function is_time_new
+!============================================================================
-function nearest_hour(sometime)
-! Return the hour nearest to the input time
-type(time_type), intent(in) :: sometime
-type(time_type) :: nearest_hour
-
-type(time_type) :: thirty, tplus30
-integer :: days, hours, secs
-
-thirty = set_time(60*30-1,0) ! almost thirty minutes
-tplus30 = sometime + thirty
-call get_time(tplus30, secs, days)
-hours = secs/(60*60)
-
-nearest_hour = set_time(hours*60*60, days)
-
-if (debug) call print_time(sometime, 'input time')
-if (debug) call print_time(nearest_hour,'top of the hour')
-
-end function nearest_hour
-
-
-
subroutine update_time(ObsTime, stationid, stationlist, timeindex)
+! The station has a list of the observation times closest to the
+! verification times.
! Add a new time to the station registry.
-! If there is no additional space, must take action.
type(time_type), intent(in) :: ObsTime
integer, intent(in) :: stationid
type(station), dimension(:), intent(inout) :: stationlist
integer, intent(in) :: timeindex
-type(time_type), allocatable, dimension(:) :: temptimes
-integer :: ntimes
-
! Update stuff that seems like a good idea,
! but I don't really know if I'll use it ...
if ( stationlist(stationid)%ntimes == 0 ) then
@@ -766,55 +784,30 @@
if ( stationlist(stationid)%last_time < ObsTime ) &
stationlist(stationid)%last_time = ObsTime
-! Do we need to make room for the new time
+if (debug) write(*,*)'Stuffing time into station ',stationid,' at timestep ', timeindex
-ntimes = size(stationlist(stationid)%times)
+! as long as ntimes /= 0 we are OK.
+! When the stations get written to the netCDF file, count the
+! number of non-zero times in the times array for a real count.
+stationlist(stationid)%ntimes = timeindex
-if ( (stationlist(stationid)%ntimes >= ntimes) .and. &
- (timeindex > stationlist(stationid)%ntimes) ) then
-
- allocate(temptimes(ntimes))
- temptimes = stationlist(stationid)%times
-
- if (associated(stationlist(stationid)%times)) then
- deallocate( stationlist(stationid)%times )
- nullify( stationlist(stationid)%times )
- endif
-
- allocate( stationlist(stationid)%times(2*ntimes) )
- stationlist(stationid)%times(1:ntimes) = temptimes(1:ntimes)
-
- deallocate(temptimes)
-
-endif
-
-! If the time is new one, increment counter
-if (timeindex > stationlist(stationid)%ntimes) then
- stationlist(stationid)%ntimes = stationlist(stationid)%ntimes + 1
- if (debug) write(*,*)'Adding a new time to ',stationid, &
- ' count ', stationlist(stationid)%ntimes, &
- ' of ', size(stationlist(stationid)%times)
-endif
-
! Stuff the time in the appropriate slot ... finally.
stationlist(stationid)%times(timeindex) = ObsTime
-if (debug) write(*,*)'Stuffing time into ',stationid, &
- ' at ', timeindex, &
- ' of ', size(stationlist(stationid)%times)
-
end subroutine update_time
+!============================================================================
+
Function InitNetCDF(fname)
character(len=*), intent(in) :: fname
integer :: InitNetCDF
integer :: ncid, i, indx1, nlines, linelen
integer :: LineLenDimID, nlinesDimID, stringDimID
-integer :: TimeDimID, StationsDimID
-integer :: VarID
+integer :: TimeDimID, StationsDimID, FcstDimID, VerifyDimID
+integer :: VarID, FcstVarID, VerifVarID, ExperimentVarID
character(len=8) :: crdate ! needed by F90 DATE_AND_TIME intrinsic
character(len=10) :: crtime ! needed by F90 DATE_AND_TIME intrinsic
@@ -822,9 +815,13 @@
integer, dimension(8) :: values ! needed by F90 DATE_AND_TIME intrinsic
character(len=129), allocatable, dimension(:) :: textblock
+real(digits12), allocatable, dimension(:) :: mytimes
+integer, allocatable, dimension(:) :: forecast_length
-integer :: nmost, ntypes
+integer :: ntypes, secs, days, ndims, mylen
+integer, dimension(nf90_max_var_dims) :: dimIDs
+
if(.not. byteSizesOK()) then
call error_handler(E_ERR,'InitNetCDF', &
'Compiler does not support required kinds of variables.',source,revision,revdate)
@@ -832,10 +829,10 @@
InitNetCDF = 0
-call nc_check(nf90_create(path = trim(fname), cmode = nf90_share, &
+call nc_check(nf90_create(path = trim(fname), cmode = nf90_clobber, &
ncid = ncid), 'obs_seq_coverage:InitNetCDF', 'create '//trim(fname))
-write(string1,*)trim(fname), ' is fortran unit ',ncid
+if (debug) write(string1,*)trim(fname), ' is fortran unit ',ncid
call error_handler(E_MSG,'InitNetCDF',string1,source,revision,revdate)
!----------------------------------------------------------------------------
@@ -854,10 +851,14 @@
'InitNetCDF', 'put_att obs_seq_coverage_revision '//trim(fname))
call nc_check(nf90_put_att(ncid, NF90_GLOBAL, 'obs_seq_coverage_revdate', revdate ), &
'InitNetCDF', 'put_att obs_seq_coverage_revdate '//trim(fname))
-call nc_check(nf90_put_att(ncid, NF90_GLOBAL, 'nTmin', nTmin ), &
- 'InitNetCDF', 'put_att nTmin '//trim(fname))
-call nc_check(nf90_put_att(ncid, NF90_GLOBAL, 'nTmax', nTmax ), &
- 'InitNetCDF', 'put_att nTmax '//trim(fname))
+call nc_check(nf90_put_att(ncid, NF90_GLOBAL, 'min_steps_required', nT_minimum ), &
+ 'InitNetCDF', 'put_att min_steps_required '//trim(fname))
+call nc_check(nf90_put_att(ncid, NF90_GLOBAL, 'forecast_length_days', &
+ forecast_length_days ), 'InitNetCDF', 'put_att forecast days '//trim(fname))
+call nc_check(nf90_put_att(ncid, NF90_GLOBAL, 'forecast_length_seconds', &
+ forecast_length_seconds ), 'InitNetCDF', 'put_att forecast seconds '//trim(fname))
+call nc_check(nf90_put_att(ncid, NF90_GLOBAL, 'verification_interval_seconds', &
+ verification_interval_seconds ), 'InitNetCDF', 'put_att verif interval '//trim(fname))
! Write all desired observation types.
! As a sanity check - do it from our working array.
@@ -898,30 +899,82 @@
!----------------------------------------------------------------------------
call nc_check(nf90_set_fill(ncid, NF90_NOFILL, i), &
- 'obs_seq_coverage:InitNetCDF', 'set_nofill '//trim(fname))
+ 'InitNetCDF', 'set_fill '//trim(fname))
+! the number of stations
+
call nc_check(nf90_def_dim(ncid=ncid, &
- name='stations', len = NF90_UNLIMITED, dimid = StationsDimID), &
- 'InitNetCDF', 'def_dim:stations '//trim(fname))
+ name='station', len = NF90_UNLIMITED, dimid = StationsDimID), &
+ 'InitNetCDF', 'def_dim:station '//trim(fname))
-call nc_check(nf90_def_var(ncid=ncid, name="stations", xtype=nf90_int, &
+call nc_check(nf90_def_var(ncid=ncid, name='station', xtype=nf90_int, &
dimids = (/ StationsDimID /), varid=VarID), &
- 'InitNetCDF', 'stations:def_var')
-call nc_check(nf90_put_att(ncid, VarID, "long_name", "desired station flag"), &
- 'InitNetCDF', 'stations:long_name')
-call nc_check(nf90_put_att(ncid, VarID, "description", "1 == good station"), &
- 'InitNetCDF', 'stations:description')
+ 'InitNetCDF', 'station:def_var')
+call nc_check(nf90_put_att(ncid, VarID, 'long_name', 'desired station flag'), &
+ 'InitNetCDF', 'station:long_name')
+call nc_check(nf90_put_att(ncid, VarID, 'description', '1 == good station'), &
+ 'InitNetCDF', 'station:description')
-! Find the station with the longest time array and define a dimension.
-nmost = 0
-do i = 1,num_stations
- if (stations(i)%ntimes > nmost) nmost = stations(i)%ntimes
-enddo
+! the number of verification times
call nc_check(nf90_def_dim(ncid=ncid, &
- name='time', len = nmost, dimid = TimeDimID), &
+ name='time', len = num_verification_times, dimid = TimeDimID), &
'InitNetCDF', 'def_dim:time '//trim(fname))
+call nc_check(nf90_def_var(ncid=ncid, name='time', xtype=nf90_double, &
+ dimids = (/ TimeDimID /), varid=VarID), &
+ 'InitNetCDF', 'time:def_var')
+call nc_check(nf90_put_att(ncid, VarID, 'long_name', 'verification time'), &
+ 'InitNetCDF', 'time:long_name')
+call nc_check(nf90_put_att(ncid, VarID, 'units', 'days since 1601-1-1'), &
+ 'InitNetCDF', 'time:put_att units')
+call nc_check(nf90_put_att(ncid, VarID, 'calendar', trim(calendar)), &
+ 'InitNetCDF', 'time:put_att calendar')
+
+! the number of supported forecasts
+
+call nc_check(nf90_def_dim(ncid=ncid, &
+ name='analysisT', len = num_analyses, dimid = FcstDimID), &
+ 'InitNetCDF', 'def_dim:analysisT '//trim(fname))
+call nc_check(nf90_def_var(ncid=ncid, name='analysisT', xtype=nf90_double, &
+ dimids = (/ FcstDimID /), varid=FcstVarID), &
+ 'InitNetCDF', 'analysisT:def_var')
+call nc_check(nf90_put_att(ncid, FcstVarID, 'long_name', 'analysis (start) time of each forecast'), &
+ 'InitNetCDF', 'analysisT:long_name')
+call nc_check(nf90_put_att(ncid, FcstVarID, 'units', 'days since 1601-1-1'), &
+ 'InitNetCDF', 'analysisT:put_att units')
+call nc_check(nf90_put_att(ncid, FcstVarID, 'calendar', trim(calendar)), &
+ 'InitNetCDF', 'analysisT:put_att calendar')
+
+! the number of verification times per forecast
+
+call nc_check(nf90_def_dim(ncid=ncid, &
+ name='forecast_lead',len= num_verify_per_fcst,dimid=VerifyDimID), &
+ 'InitNetCDF', 'def_dim:forecast_lead '//trim(fname))
+call nc_check(nf90_def_var(ncid=ncid, name='forecast_lead', xtype=nf90_int, &
+ dimids = (/ VerifyDimID /), varid=VerifVarID), &
+ 'InitNetCDF', 'forecast_lead:def_var')
+call nc_check(nf90_put_att(ncid, VerifVarID, 'long_name', 'current forecast length'), &
+ 'InitNetCDF', 'forecast_lead:long_name')
+call nc_check(nf90_put_att(ncid, VerifVarID, 'units', 'seconds'), &
+ 'InitNetCDF', 'forecast_lead:put_att units')
+
+! the verification times for each forecast
+
+call nc_check(nf90_def_var(ncid=ncid, name='verification_times', xtype=nf90_double, &
+ dimids = (/ VerifyDimID, FcstDimID /), varid=ExperimentVarID), &
+ 'InitNetCDF', 'experiment:def_var')
+call nc_check(nf90_put_att(ncid, ExperimentVarID, 'long_name', 'verification times during each forecast run'), &
+ 'InitNetCDF', 'experiment:long_name')
+call nc_check(nf90_put_att(ncid, ExperimentVarID, 'units', 'days since 1601-1-1'), &
+ 'InitNetCDF', 'experiment:put_att units')
+call nc_check(nf90_put_att(ncid, ExperimentVarID, 'calendar', trim(calendar)), &
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list