[Dart-dev] [4126] DART/trunk/observations/WOD: Updated version of this converter, and additional shell scripts
nancy at ucar.edu
nancy at ucar.edu
Thu Nov 5 14:26:32 MST 2009
Revision: 4126
Author: nancy
Date: 2009-11-05 14:26:32 -0700 (Thu, 05 Nov 2009)
Log Message:
-----------
Updated version of this converter, and additional shell scripts
to help automate obs conversions. There are date fixes in this
code which were confirmed to be input data errors or conversion
errors from the data providers. This version includes obs error
values settable by namelist, more output about what incoming QC
values were encountered, and more testing/debugging options.
Also fixed a bug where the temperature QC had been applied to
the salinity obs. Added a README to describe where the data
files came from, and some simple suggestions for the expected
workflow here.
Modified Paths:
--------------
DART/trunk/observations/WOD/WOD.html
DART/trunk/observations/WOD/wod_read_routines.f90
DART/trunk/observations/WOD/wod_to_obs.f90
Added Paths:
-----------
DART/trunk/observations/WOD/README
DART/trunk/observations/WOD/shell_scripts/
DART/trunk/observations/WOD/shell_scripts/input.nml.template
DART/trunk/observations/WOD/shell_scripts/makedaily.sh
DART/trunk/observations/WOD/shell_scripts/runyear.sh
-------------- next part --------------
Added: DART/trunk/observations/WOD/README
===================================================================
--- DART/trunk/observations/WOD/README (rev 0)
+++ DART/trunk/observations/WOD/README 2009-11-05 21:26:32 UTC (rev 4126)
@@ -0,0 +1,26 @@
+
+
+The original data files were yearly, by instrument,
+downloaded from http://dss.ucar.edu/datasets/ds285.0/
+
+We are using the STD files, which are on standard depth
+levels, and not the OBS files, which are raw and in many
+cases too dense in the vertical to be independent measurements.
+
+Look for the yearly_XXX_STD.tar files - there are about
+10 different types of obs, plus documentation.
+
+run wod_to_obs on them (in the shell_scripts dir, see runyear.ksh).
+
+then on either the DASG cluster or bluefire, go into:
+
+/ptmp/dart/Obs_sets/WOD/progs
+
+and see the makedaily.sh file for how to break the full year
+of data into daily obs files. it has code to roll over year
+and month boundaries. the files are daily, centered on 0Z,
+and so run from 12-24Z the previous day and then 0-12Z on the
+current day.
+
+
+
Modified: DART/trunk/observations/WOD/WOD.html
===================================================================
--- DART/trunk/observations/WOD/WOD.html 2009-11-05 21:22:16 UTC (rev 4125)
+++ DART/trunk/observations/WOD/WOD.html 2009-11-05 21:26:32 UTC (rev 4126)
@@ -156,7 +156,14 @@
wod_input_file,
wod_input_filelist,
wod_out_file,
- avg_obs_per_file
+ avg_obs_per_file,
+ debug,
+ max_casts,
+ no_output_file,
+ print_every_nth_cast,
+ print_qc_summary,
+ temperature_error,
+ salinity_error
</em></pre></div>
<H3 class=indent1>Discussion</H3>
<P>This namelist is read in a file called <em class=file>input.nml</em>
@@ -177,8 +184,8 @@
<TR><!--contents--><TD valign=top>wod_input_filelist</TD>
<!-- type --><TD valign=top>character(len=128)</TD>
- <!--descript--><TD>To convert a series of files in a single execution
- create a text file which contains each input file,
+ <!--descript--><TD>To convert one or more files in a single execution
+ create a text file which contains each input filename,
in ascii, one filename per line. Set this item to
the name of that file, and set 'wod_input_file' to
the empty string ('').
@@ -202,20 +209,31 @@
is computed by multiplying this number by the number of input
files. If you get an error because there is no more room to
add observations to the output file, increase this number.
+ Do not make this an unreasonably huge number, however, since
+ the code does preallocate space and will be slow if the number
+ of obs becomes very large.
Default: 500000</TD></TR>
<TR><!--contents--><TD valign=top>print_every_nth_cast</TD>
<!-- type --><TD valign=top>integer</TD>
<!--descript--><TD>If set to a value greater than 0, the program will print a
- message every N casts. This allows the user to monitor the
- progress of the conversion.
+ message after processing every N casts. This allows the user
+ to monitor the progress of the conversion.
Default: -1 (no printing)</TD></TR>
+ <TR><!--contents--><TD valign=top>print_qc_summary</TD>
+ <!-- type --><TD valign=top>logical</TD>
+ <!--descript--><TD>If set to .TRUE., the program will print out a
+ summary of the number of casts which had a non-zero
+ quality control values
+ (current files appear to use values of 1-9).
+ Default: .TRUE. (print QC summary)</TD></TR>
+
<TR><!--contents--><TD valign=top>debug</TD>
<!-- type --><TD valign=top>logical</TD>
<!--descript--><TD>If set to .TRUE., the program will print
out debugging information.
- Default: .FALSE. </TD></TR>
+ Default: .FALSE. (no debug info printed) </TD></TR>
<TR><!--contents--><TD valign=top>max_casts</TD>
<!-- type --><TD valign=top>integer</TD>
@@ -224,6 +242,36 @@
Generally only expected to be useful for debugging.
Default: -1 (convert all data)</TD></TR>
+ <TR><!--contents--><TD valign=top>no_output_file</TD>
+ <!-- type --><TD valign=top>logical</TD>
+ <!--descript--><TD>If set to .TRUE., the converter will do all the work needed
+ to convert the observations, count the number of each category
+ of QC values, etc, but will not create the final obs_seq file.
+ Can be useful if checking an input file for problems, or for
+ getting QC statistics without waiting for a full output file
+ to be constructed, which can be slow for large numbers of obs.
+ Generally only expected to be useful for debugging.
+ Default: .FALSE. (create output file)</TD></TR>
+
+ <TR><!--contents--><TD valign=top>temperature_error</TD>
+ <!-- type --><TD valign=top>real(r8)</TD>
+ <!--descript--><TD>The combined expected error of temperature observations from
+ all sources, including instrument error, model bias, and
+ representativeness error (e.g. larger or smaller grid
+ box sizes affecting expected accuracy), in degrees Centigrade.
+ Default: 0.5 degrees C. Observations in output file store
+ error variance, which will be this value squared.</TD></TR>
+
+ <TR><!--contents--><TD valign=top>salinity_error</TD>
+ <!-- type --><TD valign=top>real(r8)</TD>
+ <!--descript--><TD>The combined expected error of salinity observations from
+ all sources, including instrument error, model bias, and
+ representativeness error (e.g. larger or smaller grid
+ box sizes affecting expected accuracy) in g/kg (psu).
+ Default: 0.5 psu (g/kg). Observations in output file store
+ error variance, and use units of msu (kg/kg), so the values
+ will be this value / 1000.0, squared.</TD></TR>
+
</TABLE>
<P>
</P>
@@ -248,6 +296,9 @@
<BR><HR><BR>
<H2>FUTURE PLANS</H2>
<P>
+This converter is currently being used on WOD05 data, but as soon as the
+WOD09 data is available (data up to 2009), the same program should be
+able to convert those files as well.
</P>
<BR><HR><BR>
Added: DART/trunk/observations/WOD/shell_scripts/input.nml.template
===================================================================
--- DART/trunk/observations/WOD/shell_scripts/input.nml.template (rev 0)
+++ DART/trunk/observations/WOD/shell_scripts/input.nml.template 2009-11-05 21:26:32 UTC (rev 4126)
@@ -0,0 +1,74 @@
+
+
+! flist contains a list of input filenames to convert into a
+! single output file
+&wod_to_obs_nml
+ wod_input_file = '',
+ wod_input_filelist = 'YYYYlist',
+ wod_out_file = '../data/obs_seqYYYY.wod',
+ avg_obs_per_file = 500000,
+ print_every_nth_cast = 1000,
+ /
+
+&preprocess_nml
+ input_obs_kind_mod_file = '../../../obs_kind/DEFAULT_obs_kind_mod.F90',
+ output_obs_kind_mod_file = '../../../obs_kind/obs_kind_mod.f90',
+ input_obs_def_mod_file = '../../../obs_def/DEFAULT_obs_def_mod.F90',
+ output_obs_def_mod_file = '../../../obs_def/obs_def_mod.f90',
+ input_files = '../../../obs_def/obs_def_ocean_mod.f90',
+ /
+
+&obs_sequence_tool_nml
+ filename_seq_list = 'olist',
+ filename_out = '../YYYYMM/obs_seq.0Z.YYYYMMDD',
+ print_only = .false.,
+ gregorian_cal = .true.,
+ first_obs_days = GREG0,
+ first_obs_seconds = 43201,
+ last_obs_days = GREG1,
+ last_obs_seconds = 43200,
+ /
+
+&utilities_nml
+ module_details = .false.,
+ write_nml = 'none',
+ nmlfilename = 'convert.nml',
+ /
+
+&obs_sequence_nml
+ write_binary_obs_sequence = .false.,
+ /
+
+&obs_kind_nml
+ /
+
+&obs_def_gps_nml
+ /
+
+&location_nml
+ /
+
+# roughly monthly bins
+&schedule_nml
+ calendar = 'Gregorian',
+ first_bin_start = YYYY, 1, 1, 0, 0, 0 ,
+ first_bin_end = YYYY, 1,30, 0, 0, 0 ,
+ last_bin_end = YYYY,12,30, 0, 0, 0 ,
+ bin_interval_days = 30,
+ bin_interval_seconds = 0,
+ max_num_bins = 20,
+ print_table = .true.,
+ /
+
+&obs_seq_to_netcdf_nml
+ obs_sequence_name = 'obs_seqYYYY.wod',
+ obs_sequence_list = '',
+ lonlim1 = 0.0,
+ lonlim2 = 360.0,
+ latlim1 = -90.0,
+ latlim2 = 90.0,
+ verbose = .false.,
+ /
+
+
+
Added: DART/trunk/observations/WOD/shell_scripts/makedaily.sh
===================================================================
--- DART/trunk/observations/WOD/shell_scripts/makedaily.sh (rev 0)
+++ DART/trunk/observations/WOD/shell_scripts/makedaily.sh 2009-11-05 21:26:32 UTC (rev 4126)
@@ -0,0 +1,133 @@
+#!/bin/bash
+
+# split the monthly file into "daily" files which start at 12:01Z
+# the previous day and end at 12:00Z on the day that matches the
+# day in the filename.
+
+#BSUB -J splitobs
+#BSUB -W 4:00
+#BSUB -o %J.out
+#BSUB -e %J.err
+#BSUB -q share
+#BSUB -P 93300315
+#BSUB -n 1
+
+
+# set the first and last days to be split. can roll over
+# month and year boundaries now! note that for the first day
+# you need the previous day data available.
+
+let start_year=1998
+let start_month=1
+let start_day=2
+
+let end_year=1999
+let end_month=12
+let end_day=31
+
+EXEDIR='../work'
+DATDIR='../data'
+
+# end of things you should have to set in this script
+
+# convert the start and stop times to gregorian days, so we can
+# compute total number of days including rolling over month and
+# year boundaries. make sure all values have leading 0s if they
+# are < 10. do the end time first so we can use the same values
+# to set the initial day while we are doing the total day calc.
+
+# these outputs from advance time (with the -g flag) are
+# 2 integers: gregorian_day_number seconds
+# and since we don't set hours, minutes, or seconds, the second
+# number is always 0 and uninteresting for us.
+mon2=`printf %02d $end_month`
+day2=`printf %02d $end_day`
+end_d=(`echo ${end_year}${mon2}${day2}00 0 -g | ${EXEDIR}/advance_time`)
+#echo $end_d
+
+mon2=`printf %02d $start_month`
+day2=`printf %02d $start_day`
+start_d=(`echo ${start_year}${mon2}${day2}00 0 -g | ${EXEDIR}/advance_time`)
+#echo $start_d
+
+# these are a string in the format YYYYMMDDHH
+# do them here to prime the loop below which first takes them apart.
+currday=(`echo ${start_year}${mon2}${day2}00 0 | ${EXEDIR}/advance_time`)
+nextday=(`echo ${start_year}${mon2}${day2}00 +1d | ${EXEDIR}/advance_time`)
+prevday=(`echo ${start_year}${mon2}${day2}00 -1d | ${EXEDIR}/advance_time`)
+
+# how many total days are going to be merged (for the loop counter)
+# (pull out the first of the 2 numbers which are output from advance_time)
+let totaldays=${end_d}-${start_d}+1
+#echo $totaldays
+
+# loop over each day
+let d=1
+while (( d <= totaldays)) ; do
+#echo top of loop
+#echo $currday $nextday
+ # parse out the parts from a string which is YYYYMMDDHH
+ # both for the current day and tomorrow
+ cyear=${currday:0:4}
+ cmonth=${currday:4:2}
+ cday=${currday:6:2}
+ nyear=${nextday:0:4}
+ nmonth=${nextday:4:2}
+ nday=${nextday:6:2}
+ pyear=${prevday:0:4}
+ pmonth=${prevday:4:2}
+ pday=${prevday:6:2}
+#echo curr $cyear $cmonth $cday
+#echo next $nyear $nmonth $nday
+#echo prev $pyear $pmonth $pday
+
+ # compute the equivalent gregorian days here.
+ g=(`echo ${cyear}${cmonth}${cday}00 -1d -g | ${EXEDIR}/advance_time`)
+ greg0=${g[0]}
+ g=(`echo ${cyear}${cmonth}${cday}00 0 -g | ${EXEDIR}/advance_time`)
+ greg1=${g[0]}
+ g=(`echo ${cyear}${cmonth}${cday}00 +1d -g | ${EXEDIR}/advance_time`)
+ greg2=${g[0]}
+#echo $greg0 $greg1 $greg2
+
+ echo starting WOD obs for ${cyear}${cmonth}${cday} gregorian= $greg1
+
+ # last 12 hrs yesterdays data plus the first 12 hrs of todays
+ if [[ ${cmonth} == '01' && ${cday} == '01' ]] ; then
+ echo "${DATDIR}/obs_seq${pyear}.wod" > olist
+ echo "${DATDIR}/obs_seq${cyear}.wod" >> olist
+ else
+ echo "${DATDIR}/obs_seq${cyear}.wod" > olist
+ fi
+
+ sed -e "s/YYYY/${cyear}/g" \
+ -e "s/MM/${cmonth}/g" \
+ -e "s/PP/${pmonth}/g" \
+ -e "s/DD/${cday}/g" \
+ -e "s/GREG0/${greg0}/g" \
+ -e "s/GREG1/${greg1}/g" \
+ -e "s/GREG2/${greg2}/g" < ./input.nml.template > input.nml
+
+ # make sure output dir exists
+ if [[ ! -d ../${cyear}${cmonth} ]] ; then
+ mkdir ../${cyear}${cmonth}
+ fi
+
+ # do the extract here
+ ${EXEDIR}/obs_sequence_tool
+
+
+ # advance the day; the output is YYYYMMDD00
+ prevday=$currday
+ currday=$nextday
+ nextday=(`echo ${nyear}${nmonth}${nday}00 +1d | ${EXEDIR}/advance_time`)
+#echo $currday $nextday $prevday
+
+ # advance the loop counter
+ let d=d+1
+#echo d=$d
+
+done
+
+exit 0
+
Property changes on: DART/trunk/observations/WOD/shell_scripts/makedaily.sh
___________________________________________________________________
Added: svn:executable
+ *
Added: DART/trunk/observations/WOD/shell_scripts/runyear.sh
===================================================================
--- DART/trunk/observations/WOD/shell_scripts/runyear.sh (rev 0)
+++ DART/trunk/observations/WOD/shell_scripts/runyear.sh 2009-11-05 21:26:32 UTC (rev 4126)
@@ -0,0 +1,41 @@
+#!/bin/sh
+
+# take a year-long obs_seq file and split it into daily files
+# which start at 12Z the previous day and end at 12Z on the
+# current day.
+
+# set the first and last years.
+let start_year=1998
+let end_year=2005
+
+EXEDIR='../work'
+DATDIR='../data'
+
+# end of things you should have to set in this script
+
+let totalyears=${end_year}-${start_year}+1
+
+
+# loop over each year
+let y=1
+let year=$start_year
+
+while (( y <= totalyears )) ; do
+
+ # status/debug - comment in or out as desired.
+ echo starting processing for ${year}
+
+ ls -1 ${DATDIR}/*/*/*${year} > ${year}list
+
+ sed -e "s/YYYY/$year/" input.nml.template > input.nml
+
+ ${EXEDIR}/wod_to_obs
+
+ # advance the loop counter
+ let y=y+1
+ let year=year+1
+
+done
+
+exit 0
+
Property changes on: DART/trunk/observations/WOD/shell_scripts/runyear.sh
___________________________________________________________________
Added: svn:executable
+ *
Modified: DART/trunk/observations/WOD/wod_read_routines.f90
===================================================================
--- DART/trunk/observations/WOD/wod_read_routines.f90 2009-11-05 21:22:16 UTC (rev 4125)
+++ DART/trunk/observations/WOD/wod_read_routines.f90 2009-11-05 21:26:32 UTC (rev 4126)
@@ -106,7 +106,7 @@
subroutine WODreadDART(nf,iyear,month,iday, &
time,rlat,rlon,levels,istdlev,nvar,ip2,nsecond, &
- bmiss,ieof)
+ bmiss,castid,ieof)
! This subroutine reads in the WOD ASCII format and loads it
! into arrays which are common/shared with the calling program.
@@ -116,7 +116,7 @@
! Passed Variables:
!
! nf - file identification number for input file
-! jj - WOD cast number
+! castid - WOD cast number
! !cc - NODC country code
! !icruise - NODC cruise number
! iyear - year of cast
@@ -280,7 +280,7 @@
!dimension jtottax(maxtcode,maxtax),itaxerr(maxtcode,maxtax)
!dimension itaxorigerr(maxtcode,maxtax)
- integer :: n, nn, n0, n1, n2, i, j, k, ij, jj
+ integer :: n, nn, n0, n1, n2, i, j, k, ij, castid
integer :: nbio, inc, nchar, nlines, istartc, icruise, npinf
integer :: npi, inchad, ica, icn, ns, insec, inbio, nbothtot, itaxtot
@@ -405,7 +405,7 @@
!
! Extract header information from the cast array
!
-! jj - WOD cast number
+! castid - WOD cast number
! cc - NODC country code
! icruise - NODC cruise number
! iyear - year of cast
@@ -418,7 +418,7 @@
istartc=inc+3
read(ichar(istartc:istartc),'(i1)') inc
write(aout(3:3),'(i1)') inc
- read(ichar(istartc+1:istartc+inc),aout) jj
+ read(ichar(istartc+1:istartc+inc),aout) castid
istartc=istartc+inc+1
cc = ichar(istartc:istartc+1)
Modified: DART/trunk/observations/WOD/wod_to_obs.f90
===================================================================
--- DART/trunk/observations/WOD/wod_to_obs.f90 2009-11-05 21:22:16 UTC (rev 4125)
+++ DART/trunk/observations/WOD/wod_to_obs.f90 2009-11-05 21:26:32 UTC (rev 4126)
@@ -18,7 +18,7 @@
use types_mod, only : r8
use time_manager_mod, only : time_type, set_calendar_type, GREGORIAN, set_time,&
increment_time, get_time, set_date, operator(-), &
- print_date, operator(+)
+ print_date, operator(+), leap_year
use utilities_mod, only : initialize_utilities, find_namelist_in_file, &
check_namelist_read, nmlfileunit, do_output, &
get_next_filename, error_handler, E_ERR, E_MSG, &
@@ -100,7 +100,7 @@
character (len=80) :: name
character (len=19) :: datestr
character (len=6) :: subset
-integer :: rcode, ncid, varid, k, nfiles, num_new_obs
+integer :: rcode, ncid, varid, j, k, nfiles, num_new_obs, castid
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
@@ -142,7 +142,8 @@
type(obs_type) :: obs, prev_obs
type(time_type) :: obs_time, delta_time
-integer :: temp_type, salt_type, bad_casts
+integer :: temp_type, salt_type, good_temp, good_salt, bad_temp, bad_salt
+integer :: salt_qc(10), temp_qc(10)
!------------------------------------------------------------------------
! Declare namelist parameters
@@ -153,13 +154,18 @@
character(len=128) :: wod_out_file = 'obs_seq.wod'
integer :: avg_obs_per_file = 500000
logical :: debug = .false.
+logical :: print_qc_summary = .true.
integer :: max_casts = -1
+logical :: no_output_file = .false.
integer :: print_every_nth_cast = -1
+real(r8) :: temperature_error = 0.5 ! degrees C
+real(r8) :: salinity_error = 0.5 ! g/kg
-namelist /wod_to_obs_nml/ wod_input_file, &
- wod_input_filelist, wod_out_file, &
- avg_obs_per_file, debug, max_casts, &
- print_every_nth_cast
+namelist /wod_to_obs_nml/ &
+ 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
! start of executable code
@@ -248,29 +254,24 @@
print *, 'opening file ', trim(next_infile)
! reset anything that's per-file
- bad_casts = 0
+ good_temp = 0
+ bad_temp = 0
+ good_salt = 0
+ bad_salt = 0
+ temp_qc(:) = 0
+ salt_qc(:) = 0
cast = 1
castloop: do ! until out of data
call WODreadDART(funit,obsyear,obsmonth,obsday, &
dtime,lato,lono,levels,istdlev,nvar,ip2,nsecond, &
- bmiss,ieof)
+ bmiss,castid,ieof)
!if (ieof /= 0) print *, 'ieof is ', ieof
! this comes back as 1 when the file has all been read.
if (ieof == 1) exit castloop
- ! if the whole cast is marked with a bad qc, loop here as well.
- ! i'm seeing a few bad dates, e.g. 2/29 on non-leap years,
- ! sept 31. try testing the cast qc before decoding the date,
- ! since it is a fatal error in our set_time() routine to set
- ! an invalid date.
- if (ierror(1) /= 0) then
- bad_casts = bad_casts + 1
- goto 20 ! inc counter, cycle castloop
- endif
-
!print *, 'obsyear, obsmonth, obsday = ', obsyear, obsmonth, obsday
!print *, 'levels = ', levels
!print *, 'lato, lono = ', lato, lono
@@ -282,19 +283,30 @@
! 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.
- ! FIXME: this is the code i wanted to use, but leap_year() isn't
- ! implemented in the current time_manager. hack it for now.
- !if ((.not. leap_year(obsyear)) .and. (obsmonth == 2) .and. (obsday == 29)) then
- if ((obsyear /= 2000) .and. (obsmonth == 2) .and. (obsday == 29)) then
- print *, 'date says: ', obsyear, obsmonth, obsday
- print *, 'which does not exist in this year; skipping obs'
- goto 20 ! inc counter, cycle castloop
+ ! 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; skipping obs'
- goto 20 ! inc counter, cycle castloop
+ 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.
@@ -349,22 +361,41 @@
endif
! FIXME: these are best-guess only.
- terr = 1.0_r8 ! temp error = 1 degrees C
- serr = 0.5_r8 / 1000.0_r8 ! salinity error = 0.5 g/kg, convert to kg/kg
+ terr = temperature_error ! degrees C
+ serr = salinity_error / 1000.0_r8 ! g/kg, convert to kg/kg
if (lono < 0.0_r8) lono = lono + 360.0_r8
obslon = lono
obslat = lato
first_obs = .true.
+
+ if (have_temp) then
+ if (ierror(1) == 0) then
+ good_temp = good_temp + 1
+ else
+ bad_temp = bad_temp + 1
+ temp_qc(ierror(1)) = temp_qc(ierror(1)) + 1
+ endif
+ endif
+ if (have_salt) then
+ if (ierror(2) == 0) then
+ good_salt = good_salt + 1
+ else
+ bad_salt = bad_salt + 1
+ salt_qc(ierror(2)) = salt_qc(ierror(2)) + 1
+ endif
+ endif
+
obsloop: do k = 1, levels
obsdepth = depth(k)
! ierror is whole cast error
if (have_temp .and. ierror(1) == 0 .and. &
- temp_type > 0 .and. temp(k, 1) /= bmiss) then
+ temp_type > 0 .and. temp(k, 1) /= bmiss &
+ .and. (.not. no_output_file)) then
! set location - incoming obs are -180 to 180 in longitude;
! dart wants 0 to 360. (also, r8)
@@ -401,8 +432,9 @@
! ierror is whole cast
- if (have_salt .and. ierror(1) == 0 .and. &
- salt_type > 0 .and. temp(k, 2) /= bmiss) then
+ if (have_salt .and. ierror(2) == 0 .and. &
+ 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))
@@ -452,7 +484,36 @@
call close_file(funit)
filenum = filenum + 1
- print *, 'filename, total casts, bad casts: ', trim(next_infile), cast, bad_casts
+ if (print_qc_summary) then
+ call error_handler(E_MSG, '', '')
+ write(msgstring, *) 'input data filename: ', trim(next_infile)
+ call error_handler(E_MSG, '', msgstring)
+ call error_handler(E_MSG, '', '')
+ write(msgstring, *) 'total casts, total profiles: ', cast, &
+ good_temp+bad_temp+good_salt+bad_salt
+ call error_handler(E_MSG, '', msgstring)
+ write(msgstring, *) 'total/good/bad temperature profiles: ', &
+ good_temp+bad_temp, good_temp, bad_temp
+ call error_handler(E_MSG, '', msgstring)
+ write(msgstring, *) 'total/good/bad salinity profiles: ', &
+ good_salt+bad_salt, good_salt, bad_salt
+ call error_handler(E_MSG, '', msgstring)
+ call error_handler(E_MSG, '', '')
+ do j=1, 10
+ if (temp_qc(j) > 0) then
+ write(msgstring, *) 'temp qc of ', j, ' found ', temp_qc(j), ' times'
+ call error_handler(E_MSG, '', msgstring)
+ endif
+ enddo
+ call error_handler(E_MSG, '', '')
+ do j=1, 10
+ if (salt_qc(j) > 0) then
+ write(msgstring, *) 'salt qc of ', j, ' found ', salt_qc(j), ' times'
+ call error_handler(E_MSG, '', msgstring)
+ endif
+ enddo
+ call error_handler(E_MSG, '', '')
+ endif
end do fileloop
More information about the Dart-dev
mailing list