[Dart-dev] [3801] DART/trunk/observations: Adding the first cut at U, V support of JPL hdf-format QuikSCAT winds.
nancy at ucar.edu
nancy at ucar.edu
Wed Apr 8 22:38:14 MDT 2009
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090408/939bb289/attachment-0001.html
-------------- next part --------------
Added: DART/trunk/observations/quikscat/Makefile
===================================================================
--- DART/trunk/observations/quikscat/Makefile (rev 0)
+++ DART/trunk/observations/quikscat/Makefile 2009-04-09 04:38:13 UTC (rev 3801)
@@ -0,0 +1,52 @@
+# Makefile for Fortran QuikSCAT HDF readers
+# Ted Lungu
+#
+# HDF libraries used to link the sample program
+HDFDIR= /usr/local/hdf
+HDFDIR= /contrib/hdf5-1.8.2_intel-10.1-64
+HDFINC = $(HDFDIR)/include
+
+LDFLAGS = -L$(HDFDIR)/lib -lmfhdf -ldf -ljpeg -lz
+#******************************************************************
+# Uncomment next 4 lines if using SUN solaris machine
+#******************************************************************
+#MACHINE = SUN
+#SYS_DEF = -lm -lnsl
+#F77 = f77
+#FFLAGS = -w -e -cg92 -Nl99 -I$(HDFINC)
+
+#******************************************************************
+# Uncomment next 4 lines if using Intel x86 Linux machine & g77
+#******************************************************************
+MACHINE = i586
+SYS_DEF = -lm
+F77 = /contrib/hdf-4.2r4_gnu-4.1.2-64/bin/h4fc
+FFLAGS = -ffixed-line-length-132 -O2 -I$(HDFINC)
+
+#******************************************************************
+# Uncomment next 4 lines if using SUN OS
+#******************************************************************
+# MACHINE = SUN
+# SYS_DEF = -lm
+# F77 = f77
+# FFLAGS = -w -I$(HDFINC)
+
+#******************************************************************
+# Uncomment next 3 lines if using SGI machine
+#******************************************************************
+# MACHINE = IRIX
+# SYS_DEF = -lm
+# F77 = f77
+# FFLAGS = -w -n32 -static -I$(HDFINC)
+
+.KEEP_STATE:
+
+read_qscat2b: read_qscat2b.o
+ $(F77) read_qscat2b.o $(FFLAGS) $(LDFLAGS) $(SYS_DEF) -o $@
+
+read_qscat2b.o: read_qscat2b.f
+ $(F77) -c $(FFLAGS) read_qscat2b.f
+
+clean:
+ rm -f *.o
+
Property changes on: DART/trunk/observations/quikscat/Makefile
___________________________________________________________________
Name: svn:mime-type
+ text/plain
Name: svn:eol-style
+ native
Added: DART/trunk/observations/quikscat/convert_L2b.f90
===================================================================
--- DART/trunk/observations/quikscat/convert_L2b.f90 (rev 0)
+++ DART/trunk/observations/quikscat/convert_L2b.f90 2009-04-09 04:38:13 UTC (rev 3801)
@@ -0,0 +1,100 @@
+! Data Assimilation Research Testbed -- DART
+! Copyright 2004-2007, Data Assimilation Research Section
+! University Corporation for Atmospheric Research
+! Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+
+program convert_L2b
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
+! Initial program to read the raw ocean observations and insert them
+! into an observation sequence. To make things easy ... we will mandate
+! an assimilatin interval of 1 day - so all observations will be
+! redefined to occur at NOON on the day they were observed.
+
+use types_mod, only : r8, deg2rad, PI
+use obs_sequence_mod, only : obs_sequence_type, write_obs_seq, &
+ static_init_obs_sequence, destroy_obs_sequence
+use quikscat_JPL_mod, only : real_obs_sequence, read_qscat2b, orbit_type, &
+ create_output_filename
+use utilities_mod, only : initialize_utilities, register_module, &
+ do_output, logfileunit, nmlfileunit, &
+ error_handler, timestamp, E_ERR, E_MSG, &
+ find_namelist_in_file, check_namelist_read
+
+implicit none
+
+! ----------------------------------------------------------------------
+! Declare local parameters
+! ----------------------------------------------------------------------
+
+character(len=256) :: datafile, output_name, dartfile, string1
+type(orbit_type) :: orbit
+type(obs_sequence_type) :: seq
+
+integer :: io, iunit
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+ source = "$URL$", &
+ revision = "$Revision$", &
+ revdate = "$Date$"
+
+! ----------------------------------------------------------------------
+! Declare namelist parameters
+! ----------------------------------------------------------------------
+
+character(len=128) :: l2b_file = ''
+character(len=128) :: datadir = '.'
+character(len=128) :: outputdir = '.'
+
+real(r8) :: lon1 = 0.0_r8, & ! lower longitude bound
+ lon2 = 360.0_r8, & ! upper longitude bound
+ lat1 = -90.0_r8, & ! lower latitude bound
+ lat2 = 90.0_r8 ! upper latitude bound
+
+namelist /convert_L2b_nml/ l2b_file, datadir, outputdir, &
+ lon1, lon2, lat1, lat2
+
+! ----------------------------------------------------------------------
+! start of executable program code
+! ----------------------------------------------------------------------
+
+call initialize_utilities('convert_L2b')
+call register_module(source,revision,revdate)
+
+! Initialize the obs_sequence module ...
+
+call static_init_obs_sequence()
+
+!----------------------------------------------------------------------
+! Read the namelist
+!----------------------------------------------------------------------
+
+call find_namelist_in_file('input.nml', 'convert_L2b_nml', iunit)
+read(iunit, nml = convert_L2b_nml, iostat = io)
+call check_namelist_read(iunit, io, 'convert_L2b_nml')
+
+! Record the namelist values used for the run ...
+write(nmlfileunit, nml=convert_L2b_nml)
+write( * , nml=convert_L2b_nml)
+
+lon1 = min(max(lon1,0.0_r8),360.0_r8)
+lon2 = min(max(lon2,0.0_r8),360.0_r8)
+if ( lon1 > lon2 ) lon2 = lon2 + 360.0_r8
+
+call create_output_filename(l2b_file, output_name)
+datafile = trim( datadir)//'/'//trim(l2b_file)
+dartfile = trim(outputdir)//'/'//trim(output_name)
+
+call read_qscat2b(datafile, orbit) ! read from HDF file into a structure
+seq = real_obs_sequence(orbit, lon1, lon2, lat1, lat2) ! convert structure to a sequence
+call write_obs_seq(seq, dartfile)
+call destroy_obs_sequence(seq) ! release the memory of the seq
+call timestamp(source,revision,revdate,'end') ! close the log file
+
+end program convert_L2b
Property changes on: DART/trunk/observations/quikscat/convert_L2b.f90
___________________________________________________________________
Name: svn:keywords
+ Date Revision Author HeadURL Id
Name: svn:eol-style
+ native
Added: DART/trunk/observations/quikscat/quikscat_JPL_mod.f90
===================================================================
--- DART/trunk/observations/quikscat/quikscat_JPL_mod.f90 (rev 0)
+++ DART/trunk/observations/quikscat/quikscat_JPL_mod.f90 2009-04-09 04:38:13 UTC (rev 3801)
@@ -0,0 +1,869 @@
+! Data Assimilation Research Testbed -- DART
+! Copyright 2004-2007, Data Assimilation Research Section
+! University Corporation for Atmospheric Research
+! Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+
+module quikscat_JPL_mod
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
+use types_mod, only : r4, r8, digits12, deg2rad, rad2deg
+
+use obs_def_mod, only : obs_def_type, get_obs_def_time, read_obs_def, &
+ write_obs_def, destroy_obs_def, interactive_obs_def, &
+ copy_obs_def, set_obs_def_time, set_obs_def_kind, &
+ set_obs_def_error_variance, set_obs_def_location
+
+use time_manager_mod, only : time_type, get_date, set_date, get_time, set_time, &
+ set_calendar_type, GREGORIAN, print_date, print_time, &
+ operator(==), operator(>), operator(<), operator(>=), &
+ operator(/=), operator(+), operator(-)
+
+use utilities_mod, only : get_unit, open_file, close_file, file_exist, &
+ register_module, error_handler, &
+ E_ERR, E_MSG, timestamp
+
+use location_mod, only : location_type, set_location, VERTISHEIGHT
+
+use obs_kind_mod, only : get_obs_kind_index, &
+ KIND_U_WIND_COMPONENT, KIND_V_WIND_COMPONENT
+
+use obs_kind_mod, only : QKSWND_U_WIND_COMPONENT, QKSWND_V_WIND_COMPONENT
+
+use obs_sequence_mod, only : init_obs_sequence, init_obs, insert_obs_in_seq, &
+ set_obs_values, set_qc, obs_sequence_type, obs_type, &
+ copy_obs, set_copy_meta_data, set_qc_meta_data, set_obs_def, &
+ get_first_obs, get_last_obs, get_obs_def
+
+implicit none
+private
+
+public :: real_obs_sequence, read_qscat2b, orbit_type, create_output_filename
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+ source = "$URL$", &
+ revision = "$Revision$", &
+ revdate = "$Date$"
+
+logical, save :: module_initialized = .false.
+character(len=128) :: msgstring
+
+logical :: DEBUG = .true.
+
+! For the 25.0km resolution WVC ... MAX_ROWS = 1624 MAX_WVC = 76
+! For the 12.5km resolution WVC ... MAX_ROWS = 3248 MAX_WVC = 152
+
+integer :: MAX_ROWS, MAX_WVC, MAX_AMBIG
+parameter (MAX_ROWS = 1624)
+parameter (MAX_WVC = 76)
+parameter (MAX_AMBIG = 4)
+
+type orbit_type
+ private
+ character(len=32) :: solution
+
+ type(time_type),dimension(MAX_ROWS) :: row_time
+ real(r4), dimension(MAX_ROWS) :: wvc_row
+
+ ! real(r4), dimension(MAX_WVC,MAX_ROWS) :: wvc_index
+ ! real(r4), dimension(MAX_WVC,MAX_ROWS) :: num_in_fore
+ ! real(r4), dimension(MAX_WVC,MAX_ROWS) :: num_in_aft
+ ! real(r4), dimension(MAX_WVC,MAX_ROWS) :: num_out_fore
+ ! real(r4), dimension(MAX_WVC,MAX_ROWS) :: num_out_aft
+ ! real(r4), dimension(MAX_WVC,MAX_ROWS) :: num_ambigs
+ ! real(r4), dimension(MAX_WVC,MAX_ROWS) :: wvc_selection
+ ! real(r4), dimension(MAX_WVC,MAX_ROWS) :: nof_rain_index
+ integer, dimension(MAX_WVC,MAX_ROWS) :: wvc_index
+ integer, dimension(MAX_WVC,MAX_ROWS) :: num_in_fore
+ integer, dimension(MAX_WVC,MAX_ROWS) :: num_in_aft
+ integer, dimension(MAX_WVC,MAX_ROWS) :: num_out_fore
+ integer, dimension(MAX_WVC,MAX_ROWS) :: num_out_aft
+ integer, dimension(MAX_WVC,MAX_ROWS) :: num_ambigs
+ integer, dimension(MAX_WVC,MAX_ROWS) :: wvc_selection
+ integer, dimension(MAX_WVC,MAX_ROWS) :: nof_rain_index
+
+ real(r4), dimension(MAX_WVC,MAX_ROWS) :: wvc_lat, wvc_lon
+ real(r4), dimension(MAX_WVC,MAX_ROWS) :: wvc_quality_flag
+ real(r4), dimension(MAX_WVC,MAX_ROWS) :: atten_corr
+ real(r4), dimension(MAX_WVC,MAX_ROWS) :: mp_rain_probability
+ real(r4), dimension(MAX_WVC,MAX_ROWS) :: srad_rain_rate
+
+ real(r4), dimension(MAX_AMBIG,MAX_WVC,MAX_ROWS) :: wind_speed
+ real(r4), dimension(MAX_AMBIG,MAX_WVC,MAX_ROWS) :: wind_dir
+ real(r4), dimension(MAX_AMBIG,MAX_WVC,MAX_ROWS) :: wind_speed_err
+ real(r4), dimension(MAX_AMBIG,MAX_WVC,MAX_ROWS) :: wind_dir_err
+ real(r4), dimension(MAX_AMBIG,MAX_WVC,MAX_ROWS) :: max_likelihood_est
+
+ real(r4), dimension(MAX_WVC,MAX_ROWS) :: model_speed ! NWP Wind Vector
+ real(r4), dimension(MAX_WVC,MAX_ROWS) :: model_dir ! NWP Wind Vector
+ real(r4), dimension(MAX_WVC,MAX_ROWS) :: wind_speed_selection ! DRE Wind Vector
+ real(r4), dimension(MAX_WVC,MAX_ROWS) :: wind_dir_selection ! DRE Wind Vector
+
+end type orbit_type
+
+contains
+
+
+subroutine initialize_module
+!-------------------------------------------------
+call register_module(source, revision, revdate)
+
+call set_calendar_type(GREGORIAN)
+
+module_initialized = .true.
+
+end subroutine initialize_module
+
+
+subroutine create_output_filename(l2bname, ofname)
+!-------------------------------------------------
+! The L2b filenames have a very long extension that
+! records when the data was published - not very interesting
+! for our purposes. replace with something DART-y.
+character(len=*), intent(IN) :: l2bname
+character(len=*), intent(OUT) :: ofname
+
+integer :: i, strlen, extstart
+
+strlen = len_trim(l2bname)
+
+dotloop : do i = strlen,1,-1
+ if (l2bname(i:i) == '.' ) then
+ extstart = i
+ exit dotloop
+ endif
+enddo dotloop
+
+ofname = l2bname(1:extstart)//'obs_seq_out'
+
+end subroutine create_output_filename
+
+
+
+function real_obs_sequence ( orbit, lon1, lon2, lat1, lat2 )
+!------------------------------------------------------------------------------
+! this function is to prepare data to DART sequence format
+!
+type(orbit_type), intent(in) :: orbit
+real(r8), intent(in) :: lon1, lon2, lat1, lat2
+
+integer :: max_num=MAX_WVC*MAX_ROWS*2
+
+type(obs_sequence_type) :: real_obs_sequence
+type(obs_def_type) :: obs_def
+type(obs_type) :: obs, prev_obs
+
+integer :: i, irow, iwvc, iamb, num_copies, num_qc
+integer :: days, seconds
+integer :: obs_num
+integer :: which_vert, uobstype, vobstype
+
+real(r4) :: speed, dir
+real(r8) :: lonc, lon, lat, vloc
+real(r8) :: u_obs, v_obs, u_var, v_var
+real(r8) :: aqc
+real(r8) :: sintheta, costheta, dirvar, speedvar
+
+type(time_type) :: time, pre_time
+
+character(len = 32) :: obs_kind_name
+character(len = 129) :: meta_data
+
+if ( .not. module_initialized ) call initialize_module
+
+num_copies = 1
+num_qc = 1
+
+! Initialize an obs_sequence
+call init_obs_sequence(real_obs_sequence, num_copies, num_qc, max_num)
+
+! set meta data of obs_seq
+do i = 1, num_copies
+ meta_data = 'observation'
+ call set_copy_meta_data(real_obs_sequence, i, meta_data)
+end do
+
+do i = 1, num_qc
+ meta_data = 'QC index - wvc_quality_flag'
+ call set_qc_meta_data(real_obs_sequence, i, meta_data)
+end do
+
+! Initialize the obs variables
+call init_obs( obs, num_copies, num_qc)
+call init_obs(prev_obs, num_copies, num_qc)
+
+! assign each observation the correct observation type
+uobstype = get_obs_kind_index('QKSWND_U_WIND_COMPONENT')
+vobstype = get_obs_kind_index('QKSWND_V_WIND_COMPONENT')
+if(( uobstype < 1) .or. ( vobstype < 1)) then
+ msgstring = 'unknown observation type [QKSWND_U_WIND_COMPONENT] ... dying ...'
+ call error_handler(E_ERR,'real_obs_sequence',msgstring,source,revision,revdate)
+endif
+
+! loop over all observations within the file
+!------------------------------------------------------------------------------
+
+obs_num = 0
+vloc = 10 ! all L2b data is 10m height
+which_vert = VERTISHEIGHT
+
+!rowloop: do irow=400,403
+rowloop: do irow=1,MAX_ROWS
+
+ time = orbit%row_time(irow)
+ call get_time(time, seconds, days)
+
+ wvcloop: do iwvc=1,MAX_WVC
+
+ ! no ambiguities means no retrieval
+
+ if (orbit%num_ambigs(iwvc,irow) < 1) cycle wvcloop
+
+ ! only use data flagged as 'best'
+ ! this should be namelist controlled at some point
+
+ aqc = orbit%wvc_quality_flag(iwvc,irow)
+ if ( aqc /= 0 ) cycle wvcloop
+
+ lat = orbit%wvc_lat(iwvc,irow) ! valid range [-90.00, 90.00]
+ lon = orbit%wvc_lon(iwvc,irow) ! valid range [ 0.00, 359.99]
+
+ ! verify the location is not outside valid limits
+ if((lon > 360.0_r8) .or. (lon < 0.0_r8) .or. &
+ (lat > 90.0_r8) .or. (lat < -90.0_r8)) then
+ write(*,*)'WARNING : invalid location. wvc,row,lon,lat = ', iwvc,irow,lon,lat
+ cycle wvcloop
+ endif
+
+ ! reject observations outside the bounding box (allowing wrapping)
+ lonc = lon
+ if (lon2 > 360.0_r8 .and. lon < 180.0_r8) lonc = lon + 360.0_r8
+
+ if(( lat < lat1) .or. ( lat > lat2) .or. &
+ (lonc < lon1) .or. (lonc > lon2)) cycle wvcloop
+
+ ! QuikSCAT uses the oceanographic/flow convention ...
+ ! 0.0 is TOWARD the north - in direct contradiction to
+ ! atmospheric convention. Convert by adding 180 modulo 360
+
+ ! using the 'selected wind' ... aka highest-ranked ambiguity,
+ ! Not DIRTH nor NWP ... should also be namelist selected
+
+ iamb = orbit%wvc_selection( iwvc,irow)
+ speed = orbit%wind_speed( iamb,iwvc,irow)
+ dir = mod(orbit%wind_dir( iamb,iwvc,irow)+180.0_r4, 360.0_r4)
+
+ if ( speed < 1.0_r4 ) cycle wvcloop ! everything unreliable
+
+ ! The requirements for QuikSCAT were 2 m/s speed (rms) over 3-20m/s
+ ! or 10% of the speed from 20-30 m/s and 20 degrees (rms) 3-30 m/s
+
+ dirvar = 20.0_r4*deg2rad ! 20 degree (rms) by spec
+ speedvar = max(2.0_r4, speed*0.1_r4)
+
+ sintheta = sin(dir*deg2rad)
+ costheta = cos(dir*deg2rad)
+
+ ! converting the speed and direction and uncertainties to U,V
+ ! is a bit tedious for the uncertainties. Really could/should
+ ! have a native type for speed and direction and a forward
+ ! operator that takes U,V and calculates speed/dir ...(trivial)
+ ! left for a future upgrade TJH
+
+ u_obs = speed * sintheta
+ v_obs = speed * costheta
+
+ u_var = abs(speedvar*sintheta) ! This is THE WRONG FORMULA ...
+ v_var = abs(speedvar*costheta) ! This is THE WRONG FORMULA ...
+
+ ! verify the location is not outside valid limits
+ if((lon > 360.0_r8) .or. (lon < 0.0_r8) .or. &
+ (lat > 90.0_r8) .or. (lat < -90.0_r8)) then
+ write(*,*) 'invalid location. lon,lat = ', lon, lat
+ cycle wvcloop
+ endif
+
+ obs_num = obs_num + 1
+
+ ! create the obs_def for this observation, add to sequence
+ !---------------------------------------------------------------
+
+ call real_obs(num_copies, num_qc, obs, lon, lat, vloc, u_obs, &
+ u_var, aqc, uobstype, which_vert, seconds, days)
+
+ if(obs_num == 1) then ! for the first observation
+ call insert_obs_in_seq(real_obs_sequence, obs)
+ call copy_obs(prev_obs, obs)
+ pre_time = time
+ else ! not the first observation
+ if(time == pre_time) then ! same time as previous observation
+ call insert_obs_in_seq(real_obs_sequence, obs, prev_obs)
+ call copy_obs(prev_obs, obs)
+ else ! not the same time
+ call insert_obs_in_seq(real_obs_sequence, obs)
+ call copy_obs(prev_obs, obs)
+ pre_time = time
+ endif
+ endif
+
+ ! The (paired) V observation must be at the same time,
+ ! and cannot be first, so it is easier.
+ call real_obs(num_copies, num_qc, obs, lon, lat, vloc, v_obs, &
+ v_var, aqc, vobstype, which_vert, seconds, days)
+ call insert_obs_in_seq(real_obs_sequence, obs, prev_obs)
+ call copy_obs(prev_obs, obs)
+
+ obs_num = obs_num + 1
+
+ if ( DEBUG ) then
+ write(*,*)irow, iwvc, aqc, lat, lon, &
+ sqrt(u_obs**2+v_obs**2), atan(u_obs/v_obs)*rad2deg
+ endif
+
+ enddo wvcloop
+enddo rowloop
+
+! Print a little summary
+print*, 'obs used = ', obs_num, ' obs skipped = ', 2*MAX_WVC*MAX_ROWS - obs_num
+
+if ( get_first_obs(real_obs_sequence, obs) ) then
+ call get_obs_def(obs, obs_def)
+ pre_time = get_obs_def_time(obs_def)
+ call print_time(pre_time,' first time in sequence is ')
+ call print_date(pre_time,' first date in sequence is ')
+endif
+if( get_last_obs(real_obs_sequence, obs)) then
+ call get_obs_def(obs, obs_def)
+ time = get_obs_def_time(obs_def)
+ call print_time(time,' last time in sequence is ')
+ call print_date(time,' last date in sequence is ')
+endif
+print*, ''
+
+end function real_obs_sequence
+
+
+
+subroutine real_obs(num_copies, num_qc, obs, lon, lat, vloc, obs_value, &
+ var2, aqc, obs_kind, which_vert, seconds, days)
+!------------------------------------------------------------------------------
+integer, intent(in) :: num_copies, num_qc
+type(obs_type), intent(inout) :: obs
+real(r8), intent(in) :: lon, lat, vloc, obs_value, var2, aqc
+integer, intent(in) :: obs_kind, which_vert, seconds, days
+
+integer :: i
+real(r8) :: aqc01(1), obs_value01(1)
+type(obs_def_type) :: obsdef0
+
+if ( .not. module_initialized ) call initialize_module
+
+! Does real initialization of an observation type
+
+call real_obs_def(obsdef0, lon, lat, vloc, &
+ var2, obs_kind, which_vert, seconds, days)
+call set_obs_def(obs, obsdef0)
+
+do i = 1, num_copies
+ obs_value01(1) = obs_value
+ call set_obs_values(obs, obs_value01(1:1) )
+end do
+
+do i = 1, num_qc
+ aqc01(1) = aqc
+ call set_qc(obs, aqc01(1:1))
+end do
+
+end subroutine real_obs
+
+
+
+subroutine real_obs_def(obs_def, lon, lat, vloc, &
+ var2, obs_kind, which_vert, seconds, days)
+!----------------------------------------------------------------------
+type(obs_def_type), intent(inout) :: obs_def
+real(r8),intent(in) :: lon, lat, vloc, var2
+integer, intent(in) :: obs_kind, which_vert, seconds, days
+
+type(location_type) :: loc0
+
+if ( .not. module_initialized ) call initialize_module
+
+! set obs location
+loc0 = set_location(lon, lat, vloc, which_vert )
+call set_obs_def_location(obs_def, loc0)
+
+! set obs kind
+call set_obs_def_kind(obs_def, obs_kind)
+
+call set_obs_def_time(obs_def, set_time(seconds, days) )
+call set_obs_def_error_variance(obs_def, var2)
+
+end subroutine real_obs_def
+
+
+
+subroutine read_qscat2b(l2b_file, orbit)
+!====================================================================
+!
+! Description:
+!
+! This file contains 3 subroutines in order to read the
+! QuikSCAT Level 2B data in Hierarchical Data Format (HDF).
+! The subroutines are as follows.
+!
+! 1. read_attrib_byname(): a subroutine to read the name
+! and value(s) of a global attribute
+! referenced by its name.
+!
+! 2. read_timetags(): a subroutine to read the timetag info
+! contained in the HDF VDATA
+!
+! 3. extract_sds(): a subroutine to read the contents of an
+! SDS from an HDF file
+!
+! NOTES:
+! 1. Please refer all questions concerning this program and
+! QuikSCAT data obtained from the JPL PO.DAAC to
+! qscat at podaac.jpl.nasa.gov.
+!
+! 2. The HDF library must be installed before this program will
+! work properly. The HDF library and further information
+! about HDF may be obtained from the National Center for
+! Supercomputing Applications (NCSA) at http://hdf.ncsa.uiuc.edu.
+!
+! 3. The L2B data are read in their entirety. Examples of reading
+! the QuikSCAT data by slabs can be found in read_qscat1b.f
+! read_qscat2a.f.
+!
+!====================================================================
+
+character(len=128), intent(in) :: l2b_file
+type(orbit_type), intent(out) :: orbit
+
+! Set Parameters
+
+integer :: DFACC_RDONLY
+parameter (DFACC_RDONLY = 1)
+
+! Define Variables
+
+character(len=1) :: product(8)
+character(len=21) :: TimeTags(MAX_ROWS)
+
+real(r4), dimension(MAX_WVC,MAX_ROWS) :: datmat
+
+integer :: sd_id, retn, sfstart, sfend
+integer :: irec1, irec2, itmp, irow, iwvc, iamb
+integer :: ntype, nval
+
+integer :: year, doy, hh, mm
+real(r4) :: ss
+type(time_type) :: basetime, offset
+
+if ( .not. module_initialized ) call initialize_module
+
+call error_handler(E_MSG,'read_qscat2b','FILENAME: '//trim(l2b_file))
+
+if ( .not. file_exist(l2b_file)) then
+ call error_handler(E_ERR,'read_qscat2b', &
+ 'QuikSCAT L2B file does not exist.', source, revision, revdate)
+endif
+
+! Open the HDF input file and initiate the SD interface
+sd_id = sfstart(l2b_file,DFACC_RDONLY)
+
+! Make sure that the file is a QuikSCAT Level 2B file
+! I've never been able to make this work with F90 ... crazy games are afoot.
+!----------------------------------------------------------------------
+! call read_attrib_byname(sd_id,'ShortName',ntype,nval,product)
+! write(*,*) product
+! write(*,*) product .ne. 'QSCATL2B'
+!if ( product .ne. 'QSCATL2B') then
+! print *,'The input file is not a QuikSCAT Level 2B file'
+! print *,'*** Aborting program ***'
+! stop
+!endif
+
+! Read the timetag info contained in the HDF VDATA
+call read_timetags(l2b_file, TimeTags)
+
+! Read each SDS in its entirety. For an example of reading
+! the QuikSCAT SDS data in slabs, please refer to read_qscat2a.f.
+
+irow=1
+call extract_sds(sd_id,'wvc_row', irow, MAX_ROWS, orbit%wvc_row)
+call extract_sds(sd_id,'wvc_lat', irow, MAX_ROWS, orbit%wvc_lat)
+call extract_sds(sd_id,'wvc_lon', irow, MAX_ROWS, orbit%wvc_lon)
+call extract_sds(sd_id,'atten_corr', irow, MAX_ROWS, orbit%atten_corr)
+call extract_sds(sd_id,'model_speed', irow, MAX_ROWS, orbit%model_speed)
+call extract_sds(sd_id,'model_dir', irow, MAX_ROWS, orbit%model_dir)
+call extract_sds(sd_id,'wind_speed', irow, MAX_ROWS, orbit%wind_speed)
+call extract_sds(sd_id,'wind_dir', irow, MAX_ROWS, orbit%wind_dir)
+call extract_sds(sd_id,'wind_speed_err', irow, MAX_ROWS, orbit%wind_speed_err)
+call extract_sds(sd_id,'wind_dir_err', irow, MAX_ROWS, orbit%wind_dir_err)
+call extract_sds(sd_id,'max_likelihood_est', irow, MAX_ROWS, orbit%max_likelihood_est)
+call extract_sds(sd_id,'wind_speed_selection',irow, MAX_ROWS, orbit%wind_speed_selection)
+call extract_sds(sd_id,'wind_dir_selection', irow, MAX_ROWS, orbit%wind_dir_selection)
+call extract_sds(sd_id,'mp_rain_probability', irow, MAX_ROWS, orbit%mp_rain_probability)
+call extract_sds(sd_id,'srad_rain_rate', irow, MAX_ROWS, orbit%srad_rain_rate)
+call extract_sds(sd_id,'wvc_index', irow, MAX_ROWS, datmat)
+ orbit%wvc_index = nint(datmat)
+call extract_sds(sd_id,'num_in_fore', irow, MAX_ROWS, datmat)
+ orbit%num_in_fore = nint(datmat)
+call extract_sds(sd_id,'num_in_aft', irow, MAX_ROWS, datmat)
+ orbit%num_in_aft = nint(datmat)
+call extract_sds(sd_id,'num_out_fore', irow, MAX_ROWS, datmat)
+ orbit%num_out_fore = nint(datmat)
+call extract_sds(sd_id,'num_out_aft', irow, MAX_ROWS, datmat)
+ orbit%num_out_aft = nint(datmat)
+call extract_sds(sd_id,'wvc_quality_flag', irow, MAX_ROWS, datmat)
+ orbit%wvc_quality_flag = int(datmat)
+call extract_sds(sd_id,'num_ambigs', irow, MAX_ROWS, datmat)
+ orbit%num_ambigs = nint(datmat)
+call extract_sds(sd_id,'wvc_selection', irow, MAX_ROWS, datmat)
+ orbit%wvc_selection = nint(datmat)
+call extract_sds(sd_id,'nof_rain_index', irow, MAX_ROWS, datmat)
+ orbit%nof_rain_index = nint(datmat)
+
+! Convert time tag to a dart timetype.
+!00000000011111111112
+!12345678901234567890
+!2008-001T01:16:24.584
+
+do irow = 1,MAX_ROWS
+ read(TimeTags(irow),'( i4)')year
+ read(TimeTags(irow),'( 5x, i3)')doy
+ read(TimeTags(irow),'( 9x, i2)')hh
+ read(TimeTags(irow),'(12x, i2)')mm
+ read(TimeTags(irow),'(15x,f7.3)')ss
+
+ basetime = set_date(year,1,1,0,0,0) ! start of the (right) year
+ offset = set_time((hh*60 + mm)*60 + nint(ss), doy-1)
+ orbit%row_time(irow) = basetime + offset
+
+! if (DEBUG) call print_date(orbit%row_time(irow), TimeTags(irow))
+enddo
+
+! Print results to screen
+if ( DEBUG ) then
+RECORDLOOP : do irow=400,410 ! or 1,MAX_ROWS
+
+ write(*,*) ' '
+ write(*,*) 'TIME: ', TimeTags(irow)
+ write(*,'(''WVC ROW: '',f5.0)') orbit%wvc_row(irow)
+
+ write(*,105)
+ 105 format('WVC#',' WVC_Qual', &
+ ' WVC Latitude/Longitude', &
+ ' Selected Wind Vector', &
+ ' NWP Wind Vector', &
+ ' Num/Sel ambig', &
+ ' DRE Wind Vector', &
+ ' MUDH Prob', &
+ ' NOF Index', &
+ ' SRR')
+
+ WVCLOOP : do iwvc = 1,MAX_WVC
+ if (orbit%num_ambigs(iwvc,irow).gt.0) then
+ iamb=orbit%wvc_selection(iwvc,irow)
+ write(*,110) orbit%wvc_index(iwvc,irow), & ! WVC#
+ int(orbit%wvc_quality_flag(iwvc,irow)), & ! WVC_Qual
+ orbit%wvc_lat(iwvc,irow), orbit%wvc_lon(iwvc,irow), &
+ orbit%wind_speed(iamb,iwvc,irow), orbit%wind_dir(iamb,iwvc,irow), & ! Selected Wind Vector
+ orbit%model_speed(iwvc,irow), orbit%model_dir(iwvc,irow), & ! NWP Wind Vector
+ orbit%num_ambigs(iwvc,irow), orbit%wvc_selection(iwvc,irow), & ! Num/Sel ambig
+ orbit%wind_speed_selection(iwvc,irow), orbit%wind_dir_selection(iwvc,irow), & ! DRE Wind Vector
+ orbit%mp_rain_probability(iwvc,irow), & ! MUDH Prob
+ orbit%nof_rain_index(iwvc,irow), & ! NOF Index
+ orbit%srad_rain_rate(iwvc,irow) ! SRR
+ endif
+ enddo WVCLOOP
+enddo RECORDLOOP
+endif
+
+ 110 format(i3,4x,"0X",z4.4,8x,f6.2,3x,f6.2,6x,f6.2,3x,f6.2, &
+ 4x,f6.2,3x,f6.2,7x,i2,3x,i2,3x,f6.2,3x,f6.2, &
+ 3x,f6.2,3x,i6,2x,f6.2)
+
+retn = sfend(sd_id)
+end subroutine read_qscat2b
+
+!====================================================================
+! READ_ATTRIB_BYNAME: a subroutine to read the name and
+! value(s) of a global attribute
+! referenced by its name.
+!
+! 5/14/1998 R.S. Dunbar
+!====================================================================
+
+ subroutine read_attrib_byname(sd_id, in_attr_name, &
+ num_type, n_values, fvalues)
+
+ character*(*) fvalues(*)
+
+ integer :: sd_id,num_type,n_values
+ integer :: attr_index,count,retn,n,oldn
+ integer :: sffattr,sfgainfo,sfrattr
+ integer :: ival, i
+
+ integer :: MAX_NC_NAME
+ parameter (MAX_NC_NAME=256)
+
+ character*(*) in_attr_name
+ character attr_name*(MAX_NC_NAME),attr_data*512
+ character*(MAX_NC_NAME) values(20)
+ character :: cr
+
+ if ( .not. module_initialized ) call initialize_module
+
+! Find the attribute assigned to in_attr_name
+ attr_index = sffattr(sd_id,in_attr_name)
+
+! Get information about the file attribute
+ retn = sfgainfo(sd_id,attr_index,attr_name,num_type,count)
+
+! Read the attribute data
+ retn = sfrattr(sd_id,attr_index,attr_data)
+
+ cr = char(10)
+ ival = 0
+ oldn = 1
+ 5 continue
+
+! QuikSCAT attributes have atleast three lines:
+! metadata type, array size and metadata contents
+! Use "blank spaces" to identify the end of a line
+
+ n = index(attr_data(oldn:(count-1)),cr)
+
+! Read all of the metadata lines
+ if (n .eq. 0) then
+ ival=ival+1
+ values(ival) = attr_data(oldn:(count-1))
+ goto 99
+ else
+ ival=ival+1
+ values(ival) = attr_data(oldn:(oldn+n-2))
+ endif
+ oldn=n+oldn
+ goto 5
+
+ 99 continue
+ n_values = ival - 2
+ do i=1,n_values
+ fvalues(i) = values(i+2)
+ enddo
+ return
+ end subroutine read_attrib_byname
+
+!====================================================================
+! READ_TIMETAGS: a subroutine to read the timetag info
+! contained in the HDF VDATA
+!
+! 5/1998 R.S. Dunbar
+!
+! Revisions:
+! 7/1999 Code adapted to read timetags in their entirety.
+! Commenter were also added. K.L. Perry
+!====================================================================
+ subroutine read_timetags(filename,timetags)
+
+ character(len=80) :: filename
+ character(len=21) :: timetags(*)
+ character(len=60) :: fields
+ character vdata_name*30
+ integer :: file_id,vdata_ref,vdata_id
+ integer :: n_records,interlace,vdata_size
+ integer :: hopen,vsfgid,vsfatch,vsfinq,vsfread,vfsdtch,hclose
+ integer :: retn
+
+ integer :: DFACC_RDONLY,FULL_INTERLACE
+ parameter(DFACC_RDONLY=1)
+ parameter(FULL_INTERLACE=0)
+
+ if ( .not. module_initialized ) call initialize_module
+
+! Open the HDF file
+ file_id = hopen(filename,DFACC_RDONLY,0)
+
+! Initialize the VS interface
+ call vfstart(file_id)
+
+! Get the reference number for the first vdata in the file
+ vdata_ref = -1
+ vdata_ref = vsfgid(file_id,vdata_ref)
+
+! Attach to the vdata for reading if it is found, otherwise
+! exit the program.
+ if (vdata_ref.eq.0) then
+ print *,'No Timetags were found in the HDF VDATA'
+ print *,'*** Aborting program ***'
+ stop
+ endif
+
+ vdata_id = vsfatch(file_id,vdata_ref,'r')
+
+! Get n_records
+ retn=vsfinq(vdata_id,n_records,interlace,fields,vdata_size,vdata_name)
+
+! Read the timetags
+ retn = vsfread(vdata_id,timetags,n_records,FULL_INTERLACE)
+
+! Terminate access to the vdata and to the VS interface,
+! then close the HDF file.
+
+! retn = vsfdtch(vdata_id)
+ call vsfdtch(vdata_id)
+ call vfend(file_id)
+ retn = hclose(file_id)
+
+ return
+ end subroutine read_timetags
+
+
+
+!====================================================================
+! EXTRACT_SDS: a subroutine to read the contents of an
+! SDS from an HDF file
+!
+! 5/12/1998 R.S. Dunbar
+!
+! Revisions:
+! 7/1999 Code adapted to read input in bytes as well as ints
+! and floats. Comments were also added. K.L. Perry
+!
+! 3/2000 Corrected code for 8-bit unsigned integers. "buffer"
+! was used instead of "buffer2". K.L. Perry
+!
+! 5/2000 Changed MAX_BUF_SIZE from 1000000 to 10000000.
+! Corrected code for 32-bit unsigned integers. Created
+! "buffer3" array of int*4 to correctly read in uint32.
+! K.L. Perry
+!
+!====================================================================
+ subroutine extract_sds(sd_id,in_var,irec,slab_size,out_var)
+
+ integer :: MAX_BUF_SIZE
+ parameter (MAX_BUF_SIZE=10000000)
+ integer :: sd_id,sds_index,sds_id,retn
+ integer :: rank,dim_sizes(3),data_type,nattrs,num_type
+ integer :: edge(3),stride(3),start(3),irec,slab_size
+ real(digits12) :: cal, cal_err, off, off_err
+ integer :: iprod,i,itmp
+
+ character*(*) in_var
+ character(len=256) :: name
+ integer :: sfn2index,sfselect,sfginfo,sfrdata,sfgcal,sfendacc
+
+ integer*2 :: buffer(MAX_BUF_SIZE)
+ byte buffer2(MAX_BUF_SIZE)
+ integer*4 :: buffer3(MAX_BUF_SIZE)
+ real(r4) :: out_var(MAX_BUF_SIZE)
+
+ if ( .not. module_initialized ) call initialize_module
+
+! Search for the index of "in_var"
+ sds_index = sfn2index(sd_id, in_var)
+
+! Select data set corresponding to the returned index
+ sds_id = sfselect(sd_id,sds_index)
+ retn = sfginfo(sds_id,name,rank,dim_sizes,data_type,nattrs)
+
+ do i=1,rank
+ edge(i)=dim_sizes(i)
+ start(i)=0
+ stride(i)=1
+ enddo
+ edge(rank)=slab_size
+ start(rank)=irec-1
+
+ iprod=1
+ do i=1,rank
+ iprod=iprod*edge(i)
+ enddo
+
+! Get the calibration and offset values of input
+ retn = sfgcal(sds_id,cal,cal_err,off,off_err,num_type)
+
+! Read Arrays which are not float32 or int8 or uint8 or uint32
+ if ((data_type.ne. 5).and.(data_type.ne.20).and. &
+ (data_type.ne.21).and.(data_type.ne.25)) then
+
+! Read the data set into the "buffer" array
+ retn=sfrdata(sds_id,start,stride,edge,buffer)
+
+! Calibrate the output
+ do i=1,iprod
+! Correct for 16-bit unsigned integers
+ if ((data_type.eq.23).and.(buffer(i).lt.0)) then
+ out_var(i)=buffer(i)+65536.0
+
+! No correction needed for signed or positive unsigned integers
+ else
+ out_var(i)=buffer(i)
+ endif
+
+ out_var(i)=out_var(i)*cal
+ enddo
+
+! Read uint32 arrays.
+ else if (data_type.eq.25) then
+
+! Read the data set into the "buffer3" uint32 array
+ retn=sfrdata(sds_id,start,stride,edge,buffer3)
+
+! Calibrate the output
+ do i=1,iprod
+! Correct for 32-bit unsigned integers
+ if ((data_type.eq.25).and.(buffer3(i).lt.0)) then
+ out_var(i)=buffer3(i)+4294967296.0
+ else
+ out_var(i)=buffer3(i)
+ endif
+ out_var(i)=out_var(i)*cal
+ enddo
+
+! Read int8 and uint8 arrays.
+ else if ((data_type.eq.20).or.(data_type.eq.21)) then
+
+! Read the data set into the "buffer2" byte-array
+ retn=sfrdata(sds_id,start,stride,edge,buffer2)
+
+! Calibrate the output
+ do i=1,iprod
+
+! Correct for 8-bit unsigned integers
+ itmp=buffer2(i)
+ if ((data_type.eq.21).and.(buffer2(i).lt.0)) then
+ itmp=itmp+256
+ out_var(i)=itmp
+
+! No correction needed for signed or positive unsigned integers
+ else
+ out_var(i)=itmp
+ endif
+
+ out_var(i)=out_var(i)*cal
+ enddo
+
+ else
+! Read float32 arrays directly into the "out_var" array
+ retn=sfrdata(sds_id,start,stride,edge,out_var)
+
+! Calibrate the output
+ do i=1,iprod
+ out_var(i)=out_var(i)*cal
+ enddo
+ endif
+
+! Terminate access to the array data set.
+ retn = sfendacc(sds_id)
+ end subroutine extract_sds
+
+end module quikscat_JPL_mod
Property changes on: DART/trunk/observations/quikscat/quikscat_JPL_mod.f90
___________________________________________________________________
Name: svn:keywords
+ Date Revision Author HeadURL Id
Name: svn:eol-style
+ native
Added: DART/trunk/observations/quikscat/read_qscat2b.f
===================================================================
--- DART/trunk/observations/quikscat/read_qscat2b.f (rev 0)
+++ DART/trunk/observations/quikscat/read_qscat2b.f 2009-04-09 04:38:13 UTC (rev 3801)
@@ -0,0 +1,489 @@
+ program read_qscat2b
+
+c====================================================================
+c Filename: read_qscat2b.f
+c
+c Usage:
+c
+c To run this program, use the following command:
+c
+c your_computer% read_qscat2b <filename>
+c
+c where "<filename>" is the name of the QuikSCAT Level 2B
+c input file
+c
+c Description:
+c
+c This file contains 3 subroutines in order to read the
+c QuikSCAT Level 2B data in Hierarchical Data Format (HDF).
+c The subroutines are as follows.
+c
+c 1. read_attrib_byname(): a subroutine to read the name
+c and value(s) of a global attribute
+c referenced by its name.
+c
+c 2. read_timetags(): a subroutine to read the timetag info
+c contained in the HDF VDATA
+c
+c 3. extract_sds(): a subroutine to read the contents of an
+c SDS from an HDF file
+c
+c NOTES:
+c 1. Please refer all questions concerning this program and
+c QuikSCAT data obtained from the JPL PO.DAAC to
+c qscat at podaac.jpl.nasa.gov.
+c
+c 2. The HDF library must be installed before this program will
+c work properly. The HDF library and further information
+c about HDF may be obtained from the National Center for
+c Supercomputing Applications (NCSA) at http://hdf.ncsa.uiuc.edu.
+c
+c 3. The L2B data are read in their entirety. Examples of reading
+c the QuikSCAT data by slabs can be found in read_qscat1b.f
+c read_qscat2a.f.
+c
+c 7/1/1999 R.S. Dunbar, K.L. Perry
+c
+c Modifications:
+c
+c 11/18/1999 Added capability to read new SDSs,
+c wind_speed_selection and wind_dir_selection.
+c K.L. Perry
+c
+c 1/12/2000 Added capability to read new SDSs,
+c mp_rain_probability and nof_rain_index. K.L. Perry
+c
+c 3/10/2000 Corrected handling of 8-bit unsigned integers
+c in extract_sds subroutine. K.L. Perry
+c
+c 5/15/2000 Corrected handling of 32-bit unsigned integers
+c in extract_sds subroutine. This change should
+c not affect the Level 2A or 2B software. K.L. Perry
+c
+c 7/13/2006 sdrad_rain_rate and adapted code for 12.5 km product reading
+C Ted Lungu
+c Copyright 1999-2006, California Institute of Technology
+c====================================================================
+
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list