[Dart-dev] [4590] DART/trunk/observations/WOD/wod_to_obs.f90: This has the changes to avoid bad times in the input files,
nancy at ucar.edu
nancy at ucar.edu
Fri Dec 3 15:09:34 MST 2010
Revision: 4590
Author: nancy
Date: 2010-12-03 15:09:34 -0700 (Fri, 03 Dec 2010)
Log Message:
-----------
This has the changes to avoid bad times in the input files,
and does the time conversion from day/fractional hours to
days and seconds correctly. Also a name change from CODAR
to HFRADAR to match the change in obs_def_ocean_mod.f90.
Modified Paths:
--------------
DART/trunk/observations/WOD/wod_to_obs.f90
-------------- next part --------------
Modified: DART/trunk/observations/WOD/wod_to_obs.f90
===================================================================
--- DART/trunk/observations/WOD/wod_to_obs.f90 2010-12-03 18:30:08 UTC (rev 4589)
+++ DART/trunk/observations/WOD/wod_to_obs.f90 2010-12-03 22:09:34 UTC (rev 4590)
@@ -39,6 +39,7 @@
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, &
set_obs_def_key
+use obs_kind_mod, only : get_obs_kind_name
! FIXME: i don't have all the actual kinds yet (bottle? underway? )
! but it's closer than before. must have obs_def_ocean_mod.f90 in preprocess list.
@@ -79,13 +80,13 @@
DOPPLER_W_CURRENT_COMPONENT, KIND_W_CURRENT_COMPONENT, &
SATELLITE_MICROWAVE_SST, SATELLITE_INFRARED_SST, &
SATELLITE_SSH, SATELLITE_SSS, &
- CODAR_RADIAL_VELOCITY, KIND_VELOCITY
+ HFRADAR_RADIAL_VELOCITY, KIND_VELOCITY
! not clean interface; going for working code first.
use WOD_read_routines_mod, only : WODREADDART, depth, temp, ierror, iderror, &
- maxlevel, maxcalc, kdim, maxtcode, maxtax, maxpinf, bmiss, &
- isec, sechead
+ maxlevel, maxcalc, kdim, maxtcode, maxtax, maxpinf, &
+ bmiss, isec, sechead
@@ -101,11 +102,11 @@
integer, parameter :: num_copies = 1, & ! number of copies in sequence
num_qc = 1 ! number of QC entries
-character (len=129) :: meta_data, msgstring, next_infile
+character (len=129) :: meta_data, msgstring, next_infile, cdummy
character (len=80) :: name
character (len=19) :: datestr
character (len=6) :: subset
-integer :: rcode, ncid, varid, j, k, nfiles, num_new_obs, castid
+integer :: rcode, ncid, varid, j, k, nfiles, num_new_obs, castid, l
integer :: aday, asec, dday, dsec, oday, osec
integer :: obsyear, obsmonth, obsday, obshour, obsmin, obssec
integer :: zloc, obs_num, io, iunit, filenum, dummy, i_qc, nc_rc
@@ -120,8 +121,8 @@
real :: dtime, lato, lono
real(r8) :: obslat, obslon, obsdepth
-! platform codes
-integer :: ptype(2, 16) = (/ &
+! platform codes - reshape required by latest XLF compiler.
+integer :: ptype(2, 16) = reshape( (/ &
MBT_TEMPERATURE, MBT_SALINITY, & ! ptype(1) = MBT
XBT_TEMPERATURE, XBT_SALINITY, & ! ptype(2) = XBT
DBT_TEMPERATURE, DBT_SALINITY, & ! ptype(3) = DBT
@@ -137,16 +138,18 @@
APB_TEMPERATURE, APB_SALINITY, & ! ptype(13) = animal
0, 0, & ! ptype(14) = bucket
GLIDER_TEMPERATURE, GLIDER_SALINITY, & ! ptype(15) = glider
- MBT_TEMPERATURE, MBT_SALINITY /) ! ptype(16) = microBT
+ MBT_TEMPERATURE, MBT_SALINITY /), & ! ptype(16) = microBT
+ (/ 2, 16 /) ) ! reshape 1d array into (2,16)
real(r8), allocatable :: lat(:), lon(:), dep(:), err(:) !, d_qc(:)
-type(obs_def_type) :: obs_def
type(obs_sequence_type) :: obs_seq
type(obs_type) :: obs, prev_obs
type(time_type) :: obs_time, delta_time, prev_time
+integer :: histbin(-1:26) = 0 ! for time debug
+integer :: histcount = 0 ! for time debug
integer :: temp_type, salt_type, good_temp, good_salt, bad_temp, bad_salt
integer :: salt_qc(10), temp_qc(10)
@@ -159,6 +162,7 @@
character(len=128) :: wod_out_file = 'obs_seq.wod'
integer :: avg_obs_per_file = 500000
logical :: debug = .false.
+logical :: timedebug = .false.
logical :: print_qc_summary = .true.
integer :: max_casts = -1
logical :: no_output_file = .false.
@@ -170,7 +174,7 @@
wod_input_file, wod_input_filelist, wod_out_file, &
avg_obs_per_file, debug, max_casts, no_output_file, &
print_every_nth_cast, print_qc_summary, &
- temperature_error, salinity_error
+ temperature_error, salinity_error, timedebug
! start of executable code
@@ -278,63 +282,40 @@
! this comes back as 1 when the file has all been read.
if (ieof == 1) exit castloop
+ ! see what data we have
+ if (nvar == 0) then
+ goto 20 ! inc counter, cycle castloop
+ endif
+
!print *, 'obsyear, obsmonth, obsday = ', obsyear, obsmonth, obsday
!print *, 'levels = ', levels
!print *, 'lato, lono = ', lato, lono
- ! FIXME: need to gather:
+ ! need to gather:
! time (year, month, day, time), lat, lon, depth, obs type, ?.
! start out with converting the real time.
+ ! do error check first. if fail, cycle here
+ ! this routine alters bad days and times to be valid, or
+ ! returns a failure and we loop here.
+ if (.not. date_ok(obsyear, obsmonth, obsday, &
+ dtime, bmiss, castid)) goto 20
- ! at least 2 files have dates of 2/29 on non-leap years. test for
- ! that and reject those obs. also found a 9/31, apparently without
- ! a bad cast qc. the 2/29s were examined and determined to be a bad
- ! conversion - should be mar 1 the 9/31s were recording errors and
- ! based on the other obs, should have been 9/30.
-
- ! don't do the year test unless this is leap day.
- if ((obsmonth == 2) .and. (obsday == 29)) then
- ! jan 1st always exists. set this so we can test for leap year.
- obs_time = set_date(obsyear, 1, 1)
- if (.not. leap_year(obs_time)) then
- print *, 'cast number: ', castid
- print *, 'date says: ', obsyear, obsmonth, obsday
- print *, 'which does not exist in this year; setting to mar 1'
- ! convert this to mar 1 (date conversion error)
- obsmonth = 3
- obsday = 1
- endif
- endif
- if ((obsmonth == 9) .and. (obsday == 31)) then
- print *, 'cast number: ', castid
- print *, 'date says: ', obsyear, obsmonth, obsday
- print *, 'which does not exist in any year; setting to sept 30'
- ! convert this to sept 30 (apparent date input error)
- obsmonth = 9
- obsday = 30
- endif
-
- ! convert to integer days and seconds, and add on to reference time.
- if (dtime == bmiss) then
- obssec = 0
- else
- obssec = int(dtime * 86400.0)
- endif
- ! fixme?
- if (obsday == 0) then
- obsday = 1
- endif
+ ! convert fractional hours (0-24) to integer days and seconds
+ ! and add on to reference time. dtime was set to 0 if it was
+ ! missing in the date_ok() routine.
+ obssec = int(dtime * 3600.0)
delta_time = set_time(obssec)
obs_time = set_date(obsyear, obsmonth, obsday) + delta_time
call get_time(obs_time, osec, oday)
!print *, 'oday, osec = ', oday, osec
!call print_date(obs_time)
-
- ! see what data we have
- if (nvar == 0) then
- goto 20 ! inc counter, cycle castloop
+ if (debug) then
+ print *, ' --- '
+ print *, 'cast number: ', castid
+ print *, 'obsyear,month,day,time = ', obsyear, obsmonth, obsday, dtime
+ print *, 'lato, lono = ', lato, lono
endif
have_temp = .false.
@@ -362,7 +343,14 @@
if (debug) then
print *, 'num levels: ', levels
print *, 'has temp, salinity: ', have_temp, have_salt
- print *, 'temp, salinity type: ', temp_type, salt_type
+ if (have_temp) then
+ cdummy = get_obs_kind_name(temp_type)
+ print *, 'temp type/name: ', temp_type, trim(cdummy)
+ endif
+ if (have_salt) then
+ cdummy = get_obs_kind_name(salt_type)
+ print *, 'salt type/name: ', salt_type, trim(cdummy)
+ endif
call print_date(obs_time, 'obs time')
endif
@@ -392,10 +380,20 @@
endif
endif
+ if (debug) then
+ if (ierror(1) /= 0) print *, 'whole temp cast discarded, ierror == ', ierror(1)
+ if (ierror(2) /= 0) print *, 'whole salt cast discarded, ierror == ', ierror(2)
+ endif
+
obsloop: do k = 1, levels
obsdepth = depth(k)
+ if (debug) then
+ write(*, "(A,F6.0,F8.3,F4.0,F8.3,F4.0)") "depth, temp/salt value,qc: ",&
+ depth(k), temp(k,1), iderror(k,1), temp(k, 2), iderror(k,2)
+ endif
+
! ierror is whole cast error
if (have_temp .and. ierror(1) == 0 .and. &
temp_type > 0 .and. temp(k, 1) /= bmiss &
@@ -403,34 +401,11 @@
! set location - incoming obs are -180 to 180 in longitude;
! dart wants 0 to 360. (also, r8)
- call set_obs_def_location(obs_def, &
- set_location(obslon, obslat, obsdepth,VERTISHEIGHT))
- call set_obs_def_kind(obs_def, temp_type)
- call set_obs_def_time(obs_def, obs_time)
-
- call set_obs_def_error_variance(obs_def, terr * terr)
- call set_obs_def_key(obs_def, obs_num)
- call set_obs_def(obs, obs_def)
-
-!print *, 'temp = ', temp(k, 1)
- obs_val(1) = temp(k, 1)
- call set_obs_values(obs, obs_val)
- qc_val(1) = iderror(k, 1) ! iderror is per-depth error
- call set_qc(obs, qc_val)
-
- ! first one, insert with no prev. otherwise, since all times will be the
- ! same for this column, insert with the prev obs as the starting point.
- ! (the first insert with no prev means it will search for the right
- ! time ordered starting point.)
- if (prev_time > obs_time) then
- call insert_obs_in_seq(obs_seq, obs)
- else
- call insert_obs_in_seq(obs_seq, obs, prev_obs)
- endif
- obs_num = obs_num+1
- prev_obs = obs
- prev_time = obs_time
-
+ call fill_obs(obs, obs_num, obslon, obslat, obsdepth, temp_type, &
+ obs_time, terr, real(temp(k, 1), r8), real(iderror(k, 1), r8))
+
+ call add_obs(obs_seq, obs, obs_time, prev_obs, prev_time)
+
if (.not. did_obs) did_obs = .true.
endif
@@ -440,34 +415,11 @@
salt_type > 0 .and. temp(k, 2) /= bmiss &
.and. (.not. no_output_file)) then
- call set_obs_def_location(obs_def, &
- set_location(obslon, obslat, obsdepth,VERTISHEIGHT))
- call set_obs_def_kind(obs_def, salt_type)
- call set_obs_def_time(obs_def, obs_time)
-
- call set_obs_def_error_variance(obs_def, serr * serr)
- call set_obs_def_key(obs_def, obs_num)
- call set_obs_def(obs, obs_def)
-
-!print *, 'salt = ', temp(k, 2) / 1000.0_r8
- obs_val(1) = temp(k, 2) / 1000.0_r8
- call set_obs_values(obs, obs_val)
- qc_val(1) = iderror(k, 2)
- call set_qc(obs, qc_val)
-
- ! first one, insert with no prev. otherwise, since all times will be the
- ! same for this column, insert with the prev obs as the starting point.
- ! (the first insert with no prev means it will search for the right
- ! time ordered starting point.)
- if (prev_time > obs_time) then
- call insert_obs_in_seq(obs_seq, obs)
- else
- call insert_obs_in_seq(obs_seq, obs, prev_obs)
- endif
- obs_num = obs_num+1
- prev_obs = obs
- prev_time = obs_time
-
+ call fill_obs(obs, obs_num, obslon, obslat, obsdepth, salt_type, &
+ obs_time, serr, temp(k, 2) / 1000.0_r8, real(iderror(k, 2), r8))
+
+ call add_obs(obs_seq, obs, obs_time, prev_obs, prev_time)
+
if (.not. did_obs) did_obs = .true.
endif
@@ -536,17 +488,227 @@
if (get_num_obs(obs_seq) > 0) call destroy_obs_sequence(obs_seq)
endif
+! histogram of times, debug only
+if (timedebug) then
+ print *, "total obs examined = ", histcount
+ print *, "bin -1 is < 0"
+ print *, "bin 0 is == 0"
+ print *, "bin 25 is > 24 but != 99.9"
+ print *, "bin 26 is == 99.9"
+ print *, " "
+ print *, 'bin count'
+ do l=lbound(histbin, 1), ubound(histbin, 1)
+ print *, l, histbin(l)
+ enddo
+endif
+
call timestamp(source,revision,revdate,'end')
call finalize_utilities()
! END OF MAIN ROUTINE
-! contains
+contains
! local subroutines/functions follow
-! something to set the err var - make it a subroutine so we can muck with it
+! subroutine to fill an obs. assumes loc is 3d, and vert is in meters.
+subroutine fill_obs(obs, onum, olon, olat, odepth, otype, otime, oerr, &
+ oval, oqc)
+ type(obs_type), intent(inout) :: obs
+ integer, intent(inout) :: onum
+ real(r8), intent(in) :: olon, olat, odepth, oerr, oval, oqc
+ integer, intent(in) :: otype
+ type(time_type), intent(in) :: otime
-! subroutine to fill an obs?
+type(obs_def_type) :: obs_def
+real(r8) :: valarr(1), qcarr(1)
+call set_obs_def_location(obs_def, &
+ set_location(olon, olat, odepth, VERTISHEIGHT))
+call set_obs_def_kind(obs_def, otype)
+call set_obs_def_time(obs_def, otime)
+
+call set_obs_def_error_variance(obs_def, oerr * oerr)
+call set_obs_def_key(obs_def, onum)
+call set_obs_def(obs, obs_def)
+
+valarr(1) = oval
+call set_obs_values(obs, valarr)
+qcarr(1) = oqc
+call set_qc(obs, qcarr)
+
+onum = onum + 1
+
+end subroutine fill_obs
+
+! add an obs to the sequence. if prev_time same or earlier
+! than obs time, insert with search starting from prev obs.
+subroutine add_obs(seq, obs, obs_time, prev_obs, prev_time)
+ type(obs_sequence_type), intent(inout) :: seq
+ type(obs_type), intent(inout) :: obs, prev_obs
+ type(time_type), intent(in) :: obs_time
+ type(time_type), intent(inout) :: prev_time
+
+logical, save :: first_obs = .true.
+
+if (first_obs .or. prev_time > obs_time) then
+ call insert_obs_in_seq(obs_seq, obs)
+ first_obs = .false.
+else
+ call insert_obs_in_seq(obs_seq, obs, prev_obs)
+endif
+
+prev_obs = obs
+prev_time = obs_time
+
+end subroutine add_obs
+
+! date check, since a lot of these obs seem to have bad times/dates.
+! modify testmonth, testday, dtime if we know what the right answer is.
+! return .TRUE. if date ok, .FALSE. if bad and not fixable.
+function date_ok(testyear, testmonth, testday, dtime, bmiss, castid)
+ integer, intent(inout) :: testyear, testmonth, testday
+ real, intent(inout) :: dtime ! note real, not real(r8), to match
+ real, intent(in) :: bmiss ! WOD-supplied code which uses reals
+ integer, intent(in) :: castid
+
+ logical :: date_ok
+
+integer :: days_per_month(12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)
+type(time_type) :: test_time
+
+! at least 2 files have dates of 2/29 on non-leap years. test for
+! that and fix those dates. also found a 9/31, apparently without
+! a bad cast qc. the 2/29s were examined and determined to be a bad
+! conversion - should be mar 1 the 9/31s were recording errors and
+! based on the other obs, should have been 9/30.
+
+10 format (A,I9,A,I6,I4,I4,F6.2)
+
+! bad month - return fail
+if (testmonth <= 0 .or. testmonth > 12) then
+ write(*, 10) 'cast number: ', castid, ' date, dtime: ', testyear, testmonth, testday, dtime
+ write(*, *) ' bad month number, discarding observation'
+ date_ok = .FALSE.
+ return
+endif
+
+! try to estimate the time shifting error impacts from old code.
+if (timedebug) then
+ histcount = histcount + 1
+ if (dtime == bmiss) then
+ ! no count. it is ok since we discarded it before
+ else if (dtime < 0.0) then
+ histbin(-1) = histbin(-1) + 1
+ else if (dtime == 0.0) then
+ histbin(0) = histbin(0) + 1
+ else if (dtime == 99.9) then
+ histbin(26) = histbin(26) + 1
+ else if (dtime > 24.0) then
+ histbin(25) = histbin(25) + 1
+ else
+ histbin(int(dtime)) = histbin(int(dtime)) + 1
+ endif
+endif
+
+! check time - if standard missing value, don't squawk.
+! otherwise, note the cast number and original value, and
+! then set to 0, or fail.
+if (dtime == bmiss) then
+ dtime = 0.0
+else if (dtime == 99.90) then
+ ! seems to be an alternative missing value in some few cases.
+ ! got confirmation from WOD providers that this was an incorrect
+ ! missing value and future versions of the WOD will use the standard
+ ! 'missing_value' read out of the file (bmiss in this subroutine,
+ ! appears to be something like -999.99)
+ write(*, 10) 'cast number: ', castid, ' date, dtime: ', testyear, testmonth, testday, dtime
+ write(*, *) ' setting to 0 and continuing'
+ dtime = 0.0
+else if (dtime < 0) then
+ write(*, 10) 'cast number: ', castid, ' date, dtime: ', testyear, testmonth, testday, dtime
+ write(*, *) ' setting to 0 and continuing'
+ dtime = 0.0
+else if (dtime > 24.0) then
+ write(*, 10) 'cast number: ', castid, ' date, dtime: ', testyear, testmonth, testday, dtime
+ ! some obs seem to have more than 24 hours in the time variable.
+ ! got confirmation from WOD providers that some appear to be keypunch
+ ! entry errors. this field should never be > 24. for now, set to 0
+ ! because the date is believed to be good.
+ !write(*, *) ' discarding obs'
+ !date_ok = .FALSE.
+ !return
+ write(*, *) ' setting to 0 and continuing'
+ dtime = 0.0
+endif
+
+! good day - return success
+if (testday > 0 .and. testday <= days_per_month(testmonth)) then
+ date_ok = .TRUE.
+ return
+endif
+
+! don't do the year test unless this is leap day.
+if ((testmonth == 2) .and. (testday == 29)) then
+ ! jan 1st always exists. set this so we can test for leap year.
+ test_time = set_date(testyear, 1, 1)
+ if (.not. leap_year(test_time)) then
+ write(*, 10) 'cast number: ', castid, ' date, dtime: ', testyear, testmonth, testday, dtime
+ write(*, *) ' does not exist in this non-leap year; setting to mar 1'
+ ! convert this to mar 1 (date conversion error)
+ testmonth = 3
+ testday = 1
+ endif
+ date_ok = .TRUE.
+ return
+endif
+
+! day 0 is possibly a processing error, but actual day is currently unknown.
+! i have been advised to discard these, which i am reluctantly going to do.
+if (testday == 0) then
+ write(*, 10) 'cast number: ', castid, ' date, dtime: ', testyear, testmonth, testday, dtime
+ write(*, *) ' does not exist; discarding obs'
+ date_ok = .FALSE.
+ return
+ ! alternative code - keep obs and set day to 1 and time to 0.
+ !write(*, *) ' setting day to 1, dtime to 0'
+ !testday = 1
+ !dtime = 0.0
+ !date_ok = .TRUE.
+ return
+endif
+
+! if only 1 > end of month, set to end of month and succeed.
+! if more than 1 > end of month, fail.
+! found other files with 4/31 and 6/39 as dates.
+if (testday > days_per_month(testmonth)+1) then
+ write(*, 10) 'cast number: ', castid, ' date, dtime: ', testyear, testmonth, testday, dtime
+ write(*, *) ' does not exist; discarding obs'
+ date_ok = .FALSE.
+ return
+endif
+
+! not sure this is a good idea; maybe should discard these obs?
+if (testday == days_per_month(testmonth)+1) then
+ write(*, 10) 'cast number: ', castid, ' date, dtime: ', testyear, testmonth, testday, dtime
+ write(*, *) ' does not exist; assuming rollover to day 1, month+1, dtime to 0'
+ dtime = 0.0
+ testday = 1
+ testmonth = testmonth + 1
+ if (testmonth > 12) then
+ testmonth = 1
+ testyear = testyear + 1
+ endif
+ date_ok = .TRUE.
+ return
+endif
+
+! shouldn't get here. assume bad date.
+write(*, 10) 'cast number: ', castid, ' date, dtime: ', testyear, testmonth, testday, dtime
+write(*, *) ' bad date; discarding obs'
+date_ok = .FALSE.
+return
+
+end function date_ok
+
end program
More information about the Dart-dev
mailing list