[Dart-dev] [8147] DART/trunk/observations/gps/convert_cosmic_gps_cdf.f90: add code to look for new 'bad' global attribute.
nancy at ucar.edu
nancy at ucar.edu
Tue Jun 30 13:29:00 MDT 2015
Revision: 8147
Author: nancy
Date: 2015-06-30 13:28:59 -0600 (Tue, 30 Jun 2015)
Log Message:
-----------
add code to look for new 'bad' global attribute. if found
and if the character(len=1) value is not "0", the profile
failed the quality control process and should not be used.
also lengthened some character arrays, removed the reset
of 'first_obs' inside the loop (should only be set once at
the start of the program execution), fixed some comments to
be more accurate, and removed a redundant test for whether
we'd added new obs to the output obs_seq file.
Modified Paths:
--------------
DART/trunk/observations/gps/convert_cosmic_gps_cdf.f90
-------------- next part --------------
Modified: DART/trunk/observations/gps/convert_cosmic_gps_cdf.f90
===================================================================
--- DART/trunk/observations/gps/convert_cosmic_gps_cdf.f90 2015-06-29 20:11:35 UTC (rev 8146)
+++ DART/trunk/observations/gps/convert_cosmic_gps_cdf.f90 2015-06-30 19:28:59 UTC (rev 8147)
@@ -23,7 +23,7 @@
use utilities_mod, only : initialize_utilities, find_namelist_in_file, &
check_namelist_read, nmlfileunit, do_nml_file, &
get_next_filename, error_handler, E_ERR, E_MSG, &
- nc_check, find_textfile_dims, do_nml_term
+ nc_check, find_textfile_dims, do_nml_term, finalize_utilities
use location_mod, only : VERTISHEIGHT, set_location
use obs_sequence_mod, only : obs_sequence_type, obs_type, read_obs_seq, &
static_init_obs_sequence, init_obs, destroy_obs, &
@@ -52,13 +52,14 @@
integer, parameter :: num_copies = 1, & ! number of copies in sequence
num_qc = 1 ! number of QC entries
-character (len=129) :: msgstring, next_infile
-character (len=80) :: name
+character (len=512) :: msgstring, msgstring2, next_infile
+character (len=NF90_MAX_NAME) :: name
character (len=6) :: subset
integer :: ncid, varid, nlevels, k, nfiles, num_new_obs, oday, osec, &
iyear, imonth, iday, ihour, imin, isec, glat, glon, zloc, obs_num, &
- io, iunit, nobs, filenum, dummy
-logical :: file_exist, first_obs, did_obs, from_list = .false.
+ io, iunit, nobs, filenum, dummy, numrejected
+character (len=1) :: badqc
+logical :: file_exist, first_obs, from_list = .false.
real(r8) :: hght_miss, refr_miss, azim_miss, oerr, &
qc, lato, lono, hghto, refro, azimo, wght, nx, ny, &
nz, rfict, obsval, phs, obs_val(1), qc_val(1)
@@ -75,15 +76,15 @@
! Declare namelist parameters
!------------------------------------------------------------------------
-integer, parameter :: nmaxlevels = 200 ! max number of observation levels
+integer, parameter :: NMAXLEVELS = 200 ! max number of observation levels
-logical :: local_operator = .true. ! see html file for more on non/local
-real(r8) :: obs_levels(nmaxlevels) = -1.0_r8
-real(r8) :: ray_ds = 5000.0_r8 ! delta stepsize (m) along ray, nonlocal op
-real(r8) :: ray_htop = 15000.0_r8 ! max height (m) for nonlocal op
-character(len=128) :: gpsro_netcdf_file = 'cosmic_gps_input.nc'
-character(len=128) :: gpsro_netcdf_filelist = 'cosmic_gps_input_list'
-character(len=128) :: gpsro_out_file = 'obs_seq.gpsro'
+logical :: local_operator = .true. ! see html file for more on non/local
+real(r8) :: obs_levels(NMAXLEVELS) = -1.0_r8
+real(r8) :: ray_ds = 5000.0_r8 ! delta stepsize (m) along ray, nonlocal op
+real(r8) :: ray_htop = 15000.0_r8 ! max height (m) for nonlocal op
+character(len=256) :: gpsro_netcdf_file = 'cosmic_gps_input.nc'
+character(len=256) :: gpsro_netcdf_filelist = 'cosmic_gps_input_list'
+character(len=256) :: gpsro_out_file = 'obs_seq.gpsro'
namelist /convert_cosmic_gps_nml/ obs_levels, local_operator, ray_ds, &
ray_htop, gpsro_netcdf_file, &
@@ -110,10 +111,10 @@
! count observation levels, make sure observation levels increase from 0
nlevels = 0
-do k = 1, nmaxlevels
- if ( obs_levels(k) == -1.0_r8 ) exit
+LEVELLOOP: do k = 1, NMAXLEVELS
+ if ( obs_levels(k) == -1.0_r8 ) exit LEVELLOOP
nlevels = k
-end do
+end do LEVELLOOP
do k = 2, nlevels
if ( obs_levels(k-1) >= obs_levels(k) ) then
call error_handler(E_ERR, 'convert_cosmic_gps_cdf', &
@@ -146,30 +147,32 @@
inquire(file=gpsro_out_file, exist=file_exist)
if ( file_exist ) then
- print *, "found existing obs_seq file, appending to ", trim(gpsro_out_file)
+ write(msgstring, *) "found existing obs_seq file, appending to ", trim(gpsro_out_file)
+ write(msgstring2, *) "adding up to a maximum of ", num_new_obs, " new observations"
+ call error_handler(E_MSG, 'convert_cosmic_gps_cdf', msgstring, &
+ source, revision, revdate, text2=msgstring2)
call read_obs_seq(gpsro_out_file, 0, 0, num_new_obs, obs_seq)
else
- print *, "no existing obs_seq file, creating ", trim(gpsro_out_file)
- print *, "max entries = ", num_new_obs
- call init_obs_sequence(obs_seq, num_copies, num_qc, num_new_obs)
+ write(msgstring, *) "no previously existing obs_seq file, creating ", trim(gpsro_out_file)
+ write(msgstring2, *) "with up to a maximum of ", num_new_obs, " observations"
+ call error_handler(E_MSG, 'convert_cosmic_gps_cdf', msgstring, &
+ source, revision, revdate, text2=msgstring2)
- do k = 1, num_copies
- call set_copy_meta_data(obs_seq, k, 'COSMIC GPS observation')
- end do
- do k = 1, num_qc
- call set_qc_meta_data(obs_seq, k, 'COSMIC QC')
- end do
+ call init_obs_sequence(obs_seq, num_copies, num_qc, num_new_obs)
+ call set_copy_meta_data(obs_seq, 1, 'COSMIC GPS observation')
+ call set_qc_meta_data(obs_seq, 1, 'COSMIC QC')
end if
-did_obs = .false.
+! main loop that does either a single file or a list of files.
+! the data is currently distributed as a single profile per file.
-! main loop that does either a single file or a list of files
+filenum = 1
+numrejected = 0
-filenum = 1
fileloop: do ! until out of files
! get the single name, or the next name from a list
@@ -181,24 +184,41 @@
endif
if (next_infile == '') exit fileloop
- ! open the occultation profile, check if it is within the window
+ ! open the next occultation profile
call nc_check( nf90_open(next_infile, nf90_nowrite, ncid), 'file open', next_infile)
- call nc_check( nf90_get_att(ncid,nf90_global,'year', iyear) ,'get_att year')
- call nc_check( nf90_get_att(ncid,nf90_global,'month', imonth),'get_att month')
- call nc_check( nf90_get_att(ncid,nf90_global,'day', iday) ,'get_att day')
- call nc_check( nf90_get_att(ncid,nf90_global,'hour', ihour) ,'get_att hour')
- call nc_check( nf90_get_att(ncid,nf90_global,'minute',imin) ,'get_att minute')
- call nc_check( nf90_get_att(ncid,nf90_global,'second',isec) ,'get_att second')
+
+ ! this is new with the cosmic 2013 reprocessed obs. they have included profiles
+ ! which have failed their QC checking but added a global attribute 'bad' which
+ ! is char "1" for bad profiles. if there is no 'bad' attribute, or if there is one
+ ! and the value is "0", then it is ok to proceed. otherwise close the file and
+ ! cycle the loop to pick up the next file.
+
+ io = nf90_get_att(ncid,nf90_global,'bad', badqc)
+ if ((io == nf90_noerr) .and. (badqc /= "0")) then
+ ! profile with a bad QC
+ numrejected = numrejected + 1
+ call nc_check( nf90_close(ncid) , 'close file', next_infile)
+ filenum = filenum + 1
+ cycle fileloop
+ endif
+
+ ! process the profile
+ call nc_check( nf90_get_att(ncid,nf90_global,'year', iyear) ,'get_att year', next_infile)
+ call nc_check( nf90_get_att(ncid,nf90_global,'month', imonth),'get_att month', next_infile)
+ call nc_check( nf90_get_att(ncid,nf90_global,'day', iday) ,'get_att day', next_infile)
+ call nc_check( nf90_get_att(ncid,nf90_global,'hour', ihour) ,'get_att hour', next_infile)
+ call nc_check( nf90_get_att(ncid,nf90_global,'minute',imin) ,'get_att minute', next_infile)
+ call nc_check( nf90_get_att(ncid,nf90_global,'second',isec) ,'get_att second', next_infile)
time_obs = set_date(iyear, imonth, iday, ihour, imin, isec)
call get_time(time_obs, osec, oday)
- call nc_check( nf90_inq_dimid(ncid, "MSL_alt", varid), 'inq dimid MSL_alt')
- call nc_check( nf90_inquire_dimension(ncid, varid, name, nobs), 'inq dim MSL_alt')
+ call nc_check( nf90_inq_dimid(ncid, "MSL_alt", varid), 'inq dimid MSL_alt', next_infile)
+ call nc_check( nf90_inquire_dimension(ncid, varid, name, nobs), 'inq dim MSL_alt', next_infile)
- call nc_check( nf90_get_att(ncid,nf90_global,'lon', glon) ,'get_att lon' )
- call nc_check( nf90_get_att(ncid,nf90_global,'lat', glat) ,'get_att lat' )
- call nc_check( nf90_get_att(ncid,nf90_global,'rfict',rfict),'get_att rfict')
+ call nc_check( nf90_get_att(ncid,nf90_global,'lon', glon) ,'get_att lon', next_infile)
+ call nc_check( nf90_get_att(ncid,nf90_global,'lat', glat) ,'get_att lat', next_infile)
+ call nc_check( nf90_get_att(ncid,nf90_global,'rfict',rfict),'get_att rfict', next_infile)
rfict = rfict * 1000.0_r8
allocate(hghtp(nobs)) ; allocate(refrp(nobs))
@@ -207,36 +227,34 @@
allocate(azim(nobs))
! read the latitude array
- call nc_check( nf90_inq_varid(ncid, "Lat", varid) ,'inq varid Lat')
- call nc_check( nf90_get_var(ncid, varid, lat) ,'get var Lat')
+ call nc_check( nf90_inq_varid(ncid, "Lat", varid) ,'inq varid Lat', next_infile)
+ call nc_check( nf90_get_var(ncid, varid, lat) ,'get var Lat', next_infile)
! read the latitude array
- call nc_check( nf90_inq_varid(ncid, "Lon", varid) ,'inq varid Lon')
- call nc_check( nf90_get_var(ncid, varid, lon) ,'get var Lon')
+ call nc_check( nf90_inq_varid(ncid, "Lon", varid) ,'inq varid Lon', next_infile)
+ call nc_check( nf90_get_var(ncid, varid, lon) ,'get var Lon', next_infile)
! read the altitude array
- call nc_check( nf90_inq_varid(ncid, "MSL_alt", varid) ,'inq varid MSL_alt')
- call nc_check( nf90_get_var(ncid, varid, hght) ,'get_var MSL_alt')
- call nc_check( nf90_get_att(ncid, varid, '_FillValue', hght_miss) ,'get_att _FillValue MSL_alt')
+ call nc_check( nf90_inq_varid(ncid, "MSL_alt", varid) ,'inq varid MSL_alt', next_infile)
+ call nc_check( nf90_get_var(ncid, varid, hght) ,'get_var MSL_alt', next_infile)
+ call nc_check( nf90_get_att(ncid, varid, '_FillValue', hght_miss) ,'get_att _FillValue MSL_alt', next_infile)
! read the refractivity
- call nc_check( nf90_inq_varid(ncid, "Ref", varid) ,'inq varid Ref')
- call nc_check( nf90_get_var(ncid, varid, refr) ,'get var Ref')
- call nc_check( nf90_get_att(ncid, varid, '_FillValue', refr_miss) ,'get_att _FillValue Ref')
+ call nc_check( nf90_inq_varid(ncid, "Ref", varid) ,'inq varid Ref', next_infile)
+ call nc_check( nf90_get_var(ncid, varid, refr) ,'get var Ref', next_infile)
+ call nc_check( nf90_get_att(ncid, varid, '_FillValue', refr_miss) ,'get_att _FillValue Ref', next_infile)
! read the dew-point temperature array
- call nc_check( nf90_inq_varid(ncid, "Azim", varid) ,'inq varid Azim')
- call nc_check( nf90_get_var(ncid, varid, azim) ,'get var Azim')
- call nc_check( nf90_get_att(ncid, varid, '_FillValue', azim_miss) ,'get_att _FillValue Azim')
+ call nc_check( nf90_inq_varid(ncid, "Azim", varid) ,'inq varid Azim', next_infile)
+ call nc_check( nf90_get_var(ncid, varid, azim) ,'get var Azim', next_infile)
+ call nc_check( nf90_get_att(ncid, varid, '_FillValue', azim_miss) ,'get_att _FillValue Azim', next_infile)
- call nc_check( nf90_close(ncid) , 'close file')
+ call nc_check( nf90_close(ncid) , 'close file', next_infile)
! convert units here.
hghtp(:) = hght(:) * 1000.0_r8
refrp(:) = refr(:) * 1.0e-6_r8
- first_obs = .true.
-
obsloop2: do k = 1, nlevels
call interp_height_wght(hght, obs_levels(k), nobs, zloc, wght)
@@ -278,8 +296,8 @@
obsval = phs
oerr = 0.01_r8 * excess_obserr_percent(hghto * 0.001_r8) * obsval
- !print *, 'hghto,obsval,perc,err=', hghto, obsval, &
- ! excess_obserr_percent(hghto * 0.001_r8), oerr
+ !DEBUG print *, 'hghto,obsval,perc,err=', hghto, obsval, &
+ !DEBUG excess_obserr_percent(hghto * 0.001_r8), oerr
subset = 'GPSEXC'
end if
@@ -301,8 +319,6 @@
obs_num = obs_num+1
- if (.not. did_obs) did_obs = .true.
-
end do obsloop2
! clean up and loop if there is another input file
@@ -312,23 +328,23 @@
end do fileloop
-! done with main loop. if we added any obs to the sequence, write it out.
-if (did_obs) then
-!print *, 'ready to write, nobs = ', get_num_obs(obs_seq)
- if (get_num_obs(obs_seq) > 0) &
- call write_obs_seq(obs_seq, gpsro_out_file)
+! done with main loop. if we added any new obs to the sequence, write it out.
+if (obs_num > 1) call write_obs_seq(obs_seq, gpsro_out_file)
- ! minor stab at cleanup, in the off chance this will someday get turned
- ! into a subroutine in a module. probably not all that needs to be done,
- ! but a start.
- call destroy_obs(obs)
-! if obs == prev_obs then you can't delete the same obs twice.
-! buf if they differ, then it's a leak. for now, don't delete prev
-! since the program is exiting here anyway.
- !call destroy_obs(prev_obs)
- if (get_num_obs(obs_seq) > 0) call destroy_obs_sequence(obs_seq)
+! cleanup memory
+call destroy_obs_sequence(obs_seq)
+call destroy_obs(obs) ! do not destroy prev_obs, which is same as obs
+
+write(msgstring, *) 'processed ', filenum, ' total profiles'
+call error_handler(E_MSG, 'convert_cosmic_gps_obs', msgstring, source, revision, revdate)
+
+if (numrejected > 0) then
+ write(msgstring, *) numrejected, ' profiles rejected for bad incoming quality control'
+ call error_handler(E_MSG, 'convert_cosmic_gps_obs', msgstring, source, revision, revdate)
endif
+call finalize_utilities()
+
! END OF MAIN ROUTINE
contains
More information about the Dart-dev
mailing list