[Dart-dev] [3997] DART/trunk/observations/gps: Updated COSMIC GPS RO converter.

nancy at ucar.edu nancy at ucar.edu
Mon Aug 10 15:33:28 MDT 2009


An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090810/ae51754c/attachment-0001.html 
-------------- next part --------------
Modified: DART/trunk/observations/gps/convert_cosmic_gps_cdf.f90
===================================================================
--- DART/trunk/observations/gps/convert_cosmic_gps_cdf.f90	2009-08-10 19:13:07 UTC (rev 3996)
+++ DART/trunk/observations/gps/convert_cosmic_gps_cdf.f90	2009-08-10 21:33:27 UTC (rev 3997)
@@ -1,3 +1,8 @@
+! Data Assimilation Research Testbed -- DART
+! Copyright 2004-2008, Data Assimilation Research Section
+! University Corporation for Atmospheric Research
+! Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !
 !   convert_cosmic_gps_cdf - program that reads a COSMIC GPS observation 
@@ -7,70 +12,87 @@
 !     created June 2008 Ryan Torn, NCAR/MMM
 !
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
 program convert_cosmic_gps_cdf
 
 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(-)
-use    utilities_mod, only : initialize_utilities, find_namelist_in_file, &
-                             check_namelist_read, nmlfileunit, do_output
+use time_manager_mod, only : time_type, set_calendar_type, GREGORIAN, set_time,&
+                             increment_time, get_time, set_date, operator(-),  &
+                             print_date
+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, &
+                             nc_check, find_textfile_dims
 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, write_obs_seq, &
-                             append_obs_to_seq, init_obs_sequence, get_num_obs, &
-                             set_copy_meta_data, set_qc_meta_data, set_qc, & 
+use obs_sequence_mod, only : obs_sequence_type, obs_type, read_obs_seq,       &
+                             static_init_obs_sequence, init_obs, destroy_obs, &
+                             write_obs_seq, init_obs_sequence, get_num_obs,   &
+                             insert_obs_in_seq, destroy_obs_sequence,         &
+                             set_copy_meta_data, set_qc_meta_data, set_qc,    & 
                              set_obs_values, set_obs_def, insert_obs_in_seq
 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_error_variance, set_obs_def_location, &
                              set_obs_def_key
 use  obs_def_gps_mod, only : set_gpsro_ref
 use     obs_kind_mod, only : GPSRO_REFRACTIVITY
+
 use           netcdf
 
 implicit none
 
-character(len=19),  parameter :: gpsro_netcdf_file = 'cosmic_gps_input.nc'
-character(len=129), parameter :: gpsro_out_file    = 'obs_seq.gpsro'
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+   source   = "$URL$", &
+   revision = "$Revision$", &
+   revdate  = "$Date$"
 
+
 integer, parameter ::   num_copies = 1,   &   ! number of copies in sequence
                         num_qc     = 1        ! number of QC entries
 
-character (len=129) :: meta_data
+character (len=129) :: meta_data, msgstring, next_infile
 character (len=80)  :: name
 character (len=19)  :: datestr
 character (len=6)   :: subset
-integer :: rcode, ncid, varid, nlevels, k, &
-           aday, asec, dday, dsec, oday, osec, &
-           iyear, imonth, iday, ihour, imin, isec, &
-           glat, glon, zloc, obs_num, io, iunit, nobs
-logical :: file_exist
-real(r8) :: hght_miss, refr_miss, azim_miss, oerr, & 
-            obs_window, qc, lato, lono, hghto, refro, azimo, wght, nx, ny, & 
-            nz, ds, htop, rfict, obsval, phs, obs_val(1), qc_val(1), & 
-            ref_obserr_kuo_percent, excess_obserr_percent
+integer :: rcode, ncid, varid, nlevels, k, nfiles, num_new_obs,  &
+           aday, asec, dday, dsec, 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.
+real(r8) :: hght_miss, refr_miss, azim_miss, oerr,               & 
+            qc, lato, lono, hghto, refro, azimo, wght, nx, ny,   & 
+            nz, ds, htop, rfict, obsval, phs, obs_val(1), qc_val(1)
 
 real(r8), allocatable :: lat(:), lon(:), hght(:), refr(:), azim(:), & 
                          hghtp(:), refrp(:)
 
 type(obs_def_type)      :: obs_def
 type(obs_sequence_type) :: obs_seq
-type(obs_type)          :: obs
+type(obs_type)          :: obs, prev_obs
 type(time_type)         :: time_obs, time_anal
 
 !------------------------------------------------------------------------
 !  Declare namelist parameters
 !------------------------------------------------------------------------
 
-integer, parameter :: nmaxlevels = 200   !  maximum number of observation levels
+integer, parameter :: nmaxlevels = 200   !  max number of observation levels
 
 logical  :: local_operator = .true.
+logical  :: overwrite_time = .false.
 real(r8) :: obs_levels(nmaxlevels) = -1.0_r8
+real(r8) :: obs_window = 12.0     ! accept obs within +/- hours from anal time
 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'
 
 namelist /convert_cosmic_gps_nml/ obs_levels, local_operator, obs_window, &
-                                  ray_ds, ray_htop
+                                  ray_ds, ray_htop, gpsro_netcdf_file,    &
+                                  gpsro_netcdf_filelist, gpsro_out_file 
 
+! could add overwrite_time to namelist
+
 ! initialize some values
 obs_num = 1
 qc = 0.0_r8
@@ -96,6 +118,7 @@
 ! Record the namelist values used for the run 
 if (do_output()) write(nmlfileunit, nml=convert_cosmic_gps_nml)
 
+! namelist checks for sanity
 
 !  count observation levels, make sure observation levels increase from 0
 nlevels = 0
@@ -105,211 +128,271 @@
 end do
 do k = 2, nlevels
   if ( obs_levels(k-1) >= obs_levels(k) ) then
-    write(6,*) 'Observation levels should increase'
-    stop
+    call error_handler(E_ERR, 'convert_cosmic_gps_cdf',       &
+                       'Observation levels should increase',  &
+                       source, revision, revdate)
   end if
 end do
 
 !  should error check the window some
 if (obs_window <= 0.0_r8 .or. obs_window > 24.0_r8) then
-    write(6,*) 'Bad value for obs_window (hours)'
-    stop
+    call error_handler(E_ERR, 'convert_cosmic_gps_cdf',       &
+                       'Bad value for obs_window (hours)',    &
+                       source, revision, revdate)
 else
    obs_window = obs_window * 3600.0_r8
 endif
 
-!  open the occultation profile, check if it is within the window
-rcode = nf90_open(gpsro_netcdf_file, nf90_nowrite, ncid)
-call check( nf90_get_att(ncid, nf90_global, 'year',   iyear)  )
-call check( nf90_get_att(ncid, nf90_global, 'month',  imonth) )
-call check( nf90_get_att(ncid, nf90_global, 'day',    iday)   )
-call check( nf90_get_att(ncid, nf90_global, 'hour',   ihour)  )
-call check( nf90_get_att(ncid, nf90_global, 'minute', imin)   )
-call check( nf90_get_att(ncid, nf90_global, 'second', isec)   )
 
-time_obs = set_date(iyear, imonth, iday, ihour, imin, isec)
-call get_time(time_obs, osec, oday)
+! cannot have both a single filename and a list; the namelist must
+! shut one off.
+if (gpsro_netcdf_file /= '' .and. gpsro_netcdf_filelist /= '') then
+  call error_handler(E_ERR, 'convert_cosmic_gps_cdf',                     &
+                     'One of gpsro_netcdf_file or filelist must be NULL', &
+                     source, revision, revdate)
+endif
+if (gpsro_netcdf_filelist /= '') from_list = .true.
 
-! time1-time2 is always positive no matter the relative magnitudes
-call get_time(time_anal-time_obs,dsec,dday)
-if ( real(dsec+dday*86400) > obs_window ) then
-  write(6,*) 'Observation is outside window'
-  stop
-end if
+! need to know a reasonable max number of obs that could be added here.
+if (from_list) then
+   call find_textfile_dims(gpsro_netcdf_filelist, nfiles, dummy)
+   num_new_obs = nlevels * nfiles
+else
+   num_new_obs = nlevels
+endif
 
-call check( nf90_inq_dimid(ncid, "MSL_alt", varid) )
-call check( nf90_inquire_dimension(ncid, varid, name, nobs) )
-
-call check( nf90_get_att(ncid, nf90_global, 'lon',   glon)  )
-call check( nf90_get_att(ncid, nf90_global, 'lat',   glat)  )
-call check( nf90_get_att(ncid, nf90_global, 'rfict', rfict) )
-rfict = rfict * 1000.0_r8
-
-allocate( lat(nobs))  ;  allocate( lon(nobs))
-allocate(hght(nobs))  ;  allocate(refr(nobs))
-allocate(azim(nobs))
-
-! read the latitude array
-call check( nf90_inq_varid(ncid, "Lat", varid) )
-call check( nf90_get_var(ncid, varid, lat) )
-
-! read the latitude array
-call check( nf90_inq_varid(ncid, "Lon", varid) )
-call check( nf90_get_var(ncid, varid, lon) )
-
-! read the altitude array
-call check( nf90_inq_varid(ncid, "MSL_alt", varid) )
-call check( nf90_get_var(ncid, varid, hght) )
-call check( nf90_get_att(ncid, varid, '_FillValue', hght_miss) )
-
-! read the refractivity
-call check( nf90_inq_varid(ncid, "Ref", varid) )
-call check( nf90_get_var(ncid, varid, refr) )
-call check( nf90_get_att(ncid, varid, '_FillValue', refr_miss) )
-
-! read the dew-point temperature array
-call check( nf90_inq_varid(ncid, "Azim", varid) )
-call check( nf90_get_var(ncid, varid, azim) )
-call check( nf90_get_att(ncid, varid, '_FillValue', azim_miss) )
-
-call check( nf90_close(ncid) )
-
 !  either read existing obs_seq or create a new one
 call static_init_obs_sequence()
 call init_obs(obs, num_copies, num_qc)
+call init_obs(prev_obs, num_copies, num_qc)
 inquire(file=gpsro_out_file, exist=file_exist)
 if ( file_exist ) then
 
 print *, "found existing obs_seq file, appending to ", trim(gpsro_out_file)
-  call read_obs_seq(gpsro_out_file, 0, 0, nlevels, obs_seq)
+   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 = ", nlevels
-  call init_obs_sequence(obs_seq, num_copies, num_qc, nlevels)
+print *, "max entries = ", num_new_obs
+  call init_obs_sequence(obs_seq, num_copies, num_qc, num_new_obs)
   do k = 1, num_copies
-    meta_data = 'NCEP BUFR observation'
+    meta_data = 'COSMIC GPS observation'
     call set_copy_meta_data(obs_seq, k, meta_data)
   end do
   do k = 1, num_qc
-    meta_data = 'NCEP QC index'
+    meta_data = 'COSMIC QC'
     call set_qc_meta_data(obs_seq, k, meta_data)
   end do
 
 end if
 
 allocate(hghtp(nlevels))  ;  allocate(refrp(nlevels))
-obsloop: do k = 1, nlevels
+did_obs = .false.
 
-  call interp_height_wght(hght, obs_levels(k), nobs, zloc, wght)
-  if ( zloc < 1 )  cycle obsloop
-  hghtp(nlevels-k+1) = obs_levels(k) * 1000.0_r8
-  refrp(nlevels-k+1) = exp( wght * log(refr(zloc)) + & 
-                    (1.0_r8 - wght) * log(refr(zloc+1)) ) * 1.0e-6_r8
+! main loop that does either a single file or a list of files
 
-end do obsloop
+filenum = 1
+fileloop: do      ! until out of files
 
-obsloop2: do k = 1, nlevels
+   ! get the single name, or the next name from a list
+   if (from_list) then 
+      next_infile = get_next_filename(gpsro_netcdf_filelist, filenum)
+   else
+      next_infile = gpsro_netcdf_file
+      if (filenum > 1) next_infile = ''
+   endif
+   if (next_infile == '') exit fileloop
+  
+   !  open the occultation profile, check if it is within the window
+   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')
+   
+   time_obs = set_date(iyear, imonth, iday, ihour, imin, isec)
+   call get_time(time_obs,  osec, oday)
+   
+   ! time1-time2 is always positive no matter the relative magnitudes
+   call get_time(time_anal-time_obs,dsec,dday)
+   if ( real(dsec+dday*86400) > obs_window ) then
+     call print_date(time_obs, 'observation: ')
+     call print_date(time_anal, 'window center: ')
+     write(msgstring, *) 'Window width (hours): ', obs_window / 3600.0_r8
+     call error_handler(E_MSG, 'convert_cosmic_gps_cdf',       &
+                        msgstring, source, revision, revdate)
+     call error_handler(E_ERR, 'convert_cosmic_gps_cdf',       &
+                        'Observation is outside window',       &
+                        source, revision, revdate)
+   end if
+   
+   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_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')
+   rfict = rfict * 1000.0_r8
+   
+   allocate( lat(nobs))  ;  allocate( lon(nobs))
+   allocate(hght(nobs))  ;  allocate(refr(nobs))
+   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')
+   
+   ! 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')
+   
+   ! 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')
+   
+   ! 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')
+   
+   ! 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_close(ncid) , 'close file')
+   
+   obsloop: do k = 1, nlevels
+   
+     call interp_height_wght(hght, obs_levels(k), nobs, zloc, wght)
+     if ( zloc < 1 )  cycle obsloop
+     hghtp(nlevels-k+1) = obs_levels(k) * 1000.0_r8
+     refrp(nlevels-k+1) = exp( wght * log(refr(zloc)) + & 
+                       (1.0_r8 - wght) * log(refr(zloc+1)) ) * 1.0e-6_r8
+   
+   end do obsloop
+   
+   first_obs = .true.
+   
+   obsloop2: do k = 1, nlevels
+   
+     call interp_height_wght(hght, obs_levels(k), nobs, zloc, wght)
+     if ( zloc < 1 )  cycle obsloop2
+   
+     lato  = wght * lat(zloc)  + (1.0_r8 - wght) * lat(zloc+1)
+     lono  = wght * lon(zloc)  + (1.0_r8 - wght) * lon(zloc+1)
+     if ( lono < 0.0_r8 )  lono = lono + 360.0_r8
+     hghto = wght * hght(zloc) + (1.0_r8 - wght) * hght(zloc+1)
+     hghto = hghto * 1000.0_r8
+     refro = wght * refr(zloc) + (1.0_r8 - wght) * refr(zloc+1)
+     azimo = wght * azim(zloc) + (1.0_r8 - wght) * azim(zloc+1)
+   
+     if ( local_operator ) then
+   
+        nx        = 0.0_r8
+        ny        = 0.0_r8
+        nz        = 0.0_r8
+        ray_ds    = 0.0_r8
+        ray_htop  = 0.0_r8
+        rfict     = 0.0_r8
+   
+        obsval = refro
+        oerr   = 0.01_r8 * ref_obserr_kuo_percent(hghto * 0.001_r8) * obsval
+        subset = 'GPSREF'
+   
+     else
+   
+       !  compute tangent unit vector
+       call tanvec01(lono, lato, azimo, nx, ny, nz)
+   
+       !  compute the excess phase
+       call excess(refrp, hghtp, lono, lato, hghto, nx, & 
+                   ny, nz, rfict, ray_ds, ray_htop, phs, nlevels)
+   
+       !  if too high, phs will return as 0.  cycle loop here.
+       if (phs <= 0) cycle obsloop2
+   
+       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
+       subset = 'GPSEXC'
+   
+     end if
+   
+     call set_gpsro_ref(obs_num, nx, ny, nz, rfict, ray_ds, ray_htop, subset)
+     call set_obs_def_location(obs_def,set_location(lono,lato,hghto,VERTISHEIGHT))
+     call set_obs_def_kind(obs_def, GPSRO_REFRACTIVITY)
+     if (overwrite_time) then   ! generally do not want to do this here.
+        call set_obs_def_time(obs_def, set_time(asec, aday))
+     else
+        call set_obs_def_time(obs_def, set_time(osec, oday))
+     endif
+     call set_obs_def_error_variance(obs_def, oerr * oerr)
+     call set_obs_def_key(obs_def, obs_num)
+     call set_obs_def(obs, obs_def)
+   
+     obs_val(1) = obsval
+     call set_obs_values(obs, obs_val)
+     qc_val(1)  = qc
+     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.
+     ! (this code used to call append, but calling insert makes it work even if
+     ! the input files are processed out of strict time order, which one message
+     ! i got seemed to indicate was happening...)
+     if (first_obs) then
+        call insert_obs_in_seq(obs_seq, obs)
+        first_obs = .false.
+   !print *, 'inserting first obs'
+     else
+        call insert_obs_in_seq(obs_seq, obs, prev_obs)
+   !print *, 'inserting other obs'
+     endif
+     obs_num = obs_num+1
+     prev_obs = obs
 
-  call interp_height_wght(hght, obs_levels(k), nobs, zloc, wght)
-  if ( zloc < 1 )  cycle obsloop2
+     if (.not. did_obs) did_obs = .true.
+   
+   end do obsloop2
 
-  lato  = wght * lat(zloc)  + (1.0_r8 - wght) * lat(zloc+1)
-  lono  = wght * lon(zloc)  + (1.0_r8 - wght) * lon(zloc+1)
-  if ( lono < 0.0_r8 )  lono = lono + 360.0_r8
-  hghto = wght * hght(zloc) + (1.0_r8 - wght) * hght(zloc+1)
-  hghto = hghto * 1000.0_r8
-  refro = wght * refr(zloc) + (1.0_r8 - wght) * refr(zloc+1)
-  azimo = wght * azim(zloc) + (1.0_r8 - wght) * azim(zloc+1)
+  ! clean up and loop if there is another input file
+  deallocate( lat, lon, hght, refr, azim )
 
-  if ( local_operator ) then
+  filenum = filenum + 1
 
-     nx    = 0.0_r8
-     ny    = 0.0_r8
-     nz    = 0.0_r8
-     ds    = 0.0_r8
-     htop  = 0.0_r8
-     rfict = 0.0_r8
+end do fileloop
 
-     obsval = refro
-     oerr   = 0.01_r8 * ref_obserr_kuo_percent(hghto * 0.001_r8) * obsval
-     subset = 'GPSREF'
+! 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)
 
-  else
+   ! 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.
+!print *, 'calling destroy_obs'
+   call destroy_obs(obs)
+   call destroy_obs(prev_obs)
+print *, 'skipping destroy_seq'
+   ! get core dumps here, not sure why?
+   !if (get_num_obs(obs_seq) > 0) call destroy_obs_sequence(obs_seq)
+endif
 
-    !  compute tangent unit vector
-    call tanvec01(lono, lato, azimo, nx, ny, nz)
+! END OF MAIN ROUTINE
 
-    !  compute the excess phase
-    call excess(refrp, hghtp, lono, lato, hghto, nx, & 
-                ny, nz, rfict, ray_ds, ray_htop, phs, nlevels)
+contains
 
-    !  if too high, phs will return as 0.  cycle loop here.
-    if (phs <= 0) cycle
+! local subroutines/functions follow
 
-    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
-    subset = 'GPSEXC'
-
-  end if
-
-  call set_gpsro_ref(obs_num, nx, ny, nz, rfict, ray_ds, ray_htop, subset)
-  call set_obs_def_location(obs_def,set_location(lono,lato,hghto,VERTISHEIGHT))
-  call set_obs_def_kind(obs_def, GPSRO_REFRACTIVITY)
-  call set_obs_def_time(obs_def, set_time(osec, oday))
-  call set_obs_def_error_variance(obs_def, oerr * oerr)
-  call set_obs_def_key(obs_def, obs_num)
-  call set_obs_def(obs, obs_def)
-
-  obs_val(1) = obsval
-  call set_obs_values(obs, obs_val)
-  qc_val(1)  = qc
-  call set_qc(obs, qc_val)
-
-  ! if this gives an error, use insert instead.  slower, but times do not
-  ! need to be monotonically increasing.
-  call append_obs_to_seq(obs_seq, obs)
-  ! this will speed up if we keep track of the previous time, and pass in an
-  ! existing obs to insert after.
-  !call insert_obs_in_seq(obs_seq, obs)
-  obs_num = obs_num+1
-
-end do obsloop2
-
-if ( get_num_obs(obs_seq) > 0 ) call write_obs_seq(obs_seq, gpsro_out_file)
-
-stop
-end
-
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !
-!   check - subroutine that checks the flag from a netCDF function.  If
-!           there is an error, the message is displayed.
-!
-!    istatus - netCDF output flag
-!
-!     created Dec. 2007 Ryan Torn, NCAR/MMM
-!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-subroutine check( istatus ) 
-
-use netcdf
-
-implicit none
-
-integer, intent (in) :: istatus
-
-if(istatus /= nf90_noerr) then 
-  print*,'Netcdf error: ',trim(nf90_strerror(istatus))
-  stop
-end if
-end subroutine check
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!
 !   excess - subroutine that computes the excess phase based on the 
 !            method of Sergey et al. (2004), eq. 1 and 2
 !
@@ -815,3 +898,5 @@
 
 return
 end subroutine vprod
+
+end program

Modified: DART/trunk/observations/gps/gps.html
===================================================================
--- DART/trunk/observations/gps/gps.html	2009-08-10 19:13:07 UTC (rev 3996)
+++ DART/trunk/observations/gps/gps.html	2009-08-10 21:33:27 UTC (rev 3997)
@@ -35,6 +35,8 @@
 <A HREF="#Overview">OVERVIEW</A> / 
 <A HREF="#DataSources">DATA SOURCES</A> / 
 <A HREF="#Programs">PROGRAMS</A> /
+<A HREF="#Modules">MODULES USES</A> /
+<A HREF="#Namelist">NAMELIST</A> /
 <A HREF="#KnownBugs">KNOWN BUGS</A> /
 <A HREF="#FuturePlans">FUTURE PLANS</A> 
 </DIV>
@@ -52,7 +54,7 @@
 <P>
 GPS Radio Occultation data is being returned
 from a series of satellites as part of the
-<a href="http://www.cosmic.ucar.edu">Cosmic</a> 
+<a href="http://www.cosmic.ucar.edu">COSMIC</a> 
 project.
 The programs in this directory help to extract the
 data from the distribution files and put them into
@@ -80,16 +82,24 @@
 
 <P>
 Data from 
-<a href="http://www.cosmic.ucar.edu">Cosmic</a> 
+<a href="http://www.cosmic.ucar.edu">COSMIC</a> 
 is available by signing up on the
 <a href="http://cosmic-io.cosmic.ucar.edu/cdaac">data access</a>
 web page.
 It is delivered in
 <a href="http://www.unidata.ucar.edu/software/netcdf">netCDF</a>
 file format.
+The files we use as input to these conversion programs are
+the Level 2 data, Atmospheric Profiles (filenames include
+the string 'atmPrf').
 </P>
 
 <P>
+Currently each vertical profile is stored in a separate file,
+and there are between 1000-2000 profiles/day, so converting a day's
+worth of observations involves downloading many individual files.
+There are plans in place to bundle these profiles together in
+a tar file to make it easier to download the raw data.
 </P>
 
 <!--==================================================================-->
@@ -99,26 +109,176 @@
 <H2>PROGRAMS</H2>
 The data is distributed in
 <a href="http://www.unidata.ucar.edu/software/netcdf">netCDF</a>
-file format.
-Both Fortran and
-<a href="http://www.mathworks.com/">Matlab</a>
-programs are provided here to read the data and output
-DART obs_seq file formats.
+file format.  DART requires all observations to be in a proprietary
+format often called DART "obs_seq" format.
+The files in this directory, a combination
+of C shell scripts and a Fortran source executable,
+do this data conversion.
 <P>
+Optional namelist interface
+<A HREF="#Namelist"> <em class=code>&amp;convert_cosmic_gps_nml</em> </A>
+may be read from file <em class=file>input.nml</em>.
+</P> 
+<P>
+The work directory contains several scripts, including one which downloads
+the raw data files a day at a time (cosmic_download.csh), and one 
+which executes the conversion program (convert_script.csh).  These
+scripts make 6 hour files by default, but have options for other
+times.  The input files are downloaded a day at a time.  Be aware
+that each profile is stored in a separate netcdf file, and there
+are between 1000-2000 files/day, so the download process can be
+lengthy.  You probably want to download as a separate preprocess step
+and do not use the options to automatically delete the input files.
+Keep the files around until you are sure you are satisified with the
+output files and then delete them by hand.
+</P> 
+<P>
+The conversion executable, convert_cosmic_gps_cdf, reads the namelist
+from the file 'input.nml', but also reads an analysis time from the
+terminal or the standard input unit.  This makes it simpler to convert
+multiple files without editing the namelist.  The program prompts for
+the time string with the required time format, 
+<em class=code>yyyy-mm-dd_hh:mm:ss</em> .
 </P>
 
 <!--==================================================================-->
-<!-- Describe the bugs.                                               -->
+
+<A NAME="Modules"></A>
+<HR>
+<H2>MODULES USED</H2>
+<PRE>
+types_mod
+time_manager_mod
+utilities_mod
+location_mod
+obs_sequence_mod
+obs_def_mod
+obs_def_gps_mod
+obs_kind_mod
+netcdf
+</PRE>
+
 <!--==================================================================-->
 
+<A NAME="Namelist"></A>
+<BR><HR><BR>
+<H2>NAMELIST</H2>
+ <P>We adhere to the F90 standard of starting a namelist with an ampersand
+ '&amp;' and terminating with a slash '/'.
+ <div class=namelist><pre>
+ <em class=call>namelist / convert_cosmic_gps_nml / </em>
+    obs_levels, local_operator, obs_window, &
+    ray_ds, ray_htop, gpsro_netcdf_file,    &
+    gpsro_netcdf_filelist, gpsro_out_file
+ </em></pre></div> 
+ <H3 class=indent1>Discussion</H3>
+ <P>This namelist is read in a file called <em class=file>input.nml</em>
+ </P>
+ <TABLE border=0 cellpadding=3 width=100%>
+ <TR><TH align=left>Contents    </TH>
+     <TH align=left>Type        </TH>
+     <TH align=left>Description </TH></TR>
+   
+ <TR><!--contents--><TD valign=top>obs_levels</TD>
+     <!--  type  --><TD valign=top>integer(200)</TD>
+     <!--descript--><TD>A series of heights, in kilometers, where observations
+                        from this profile should be interpolated.  (Note that
+                        the other distances and heights in the namelist are
+                        specified in meters.)  The values should be listed in
+                        increasing height order.
+                        Default: none</TD></TR>
+   
+ <TR><!--contents--><TD valign=top>local_operator</TD>
+     <!--  type  --><TD valign=top>logical</TD>
+     <!--descript--><TD>If .true., compute the observation using a method
+                        which assumes all effects occur at the tangent point.
+                        If .false., integrate along the tangent line and do
+                        ray-path reconstruction.
+                        Default: .true.</TD></TR>
+   
+ <TR><!--contents--><TD valign=top>obs_window</TD>
+     <!--  type  --><TD valign=top>real(r8)</TD>
+     <!--descript--><TD>Accept and convert observations if they are within
+                        plus or minus this number of hours from the analysis
+                        time (which is read as an input from the terminal or
+                        the standard input unit).
+                        Default: 12.0</TD></TR>
+   
+ <TR><!--contents--><TD valign=top>ray_ds</TD>
+     <!--  type  --><TD valign=top>real(r8)</TD>
+     <!--descript--><TD>For the non-local operator only, the delta stepsize, 
+                        in meters, to use for the along-path integration in
+                        each direction out from the tangent point.
+                        Default: 5000.0</TD></TR>
+   
+ <TR><!--contents--><TD valign=top>ray_htop</TD>
+     <!--  type  --><TD valign=top>real(r8)</TD>
+     <!--descript--><TD>For the non-local operator only, stop the integration
+                        when one of the endpoints of the next integration step
+                        goes above this height.  Specified in meters.
+                        Default: 15000.0</TD></TR>
+   
+ <TR><!--contents--><TD valign=top>gpsro_netcdf_file</TD>
+     <!--  type  --><TD valign=top>character(len=128)</TD>
+     <!--descript--><TD>The input filename when converting a single profile.  
+                        Only one of the 2 filenames can have a valid value, 
+                        so to use the single filename set the list name
+                        ('gpsro_netcdf_filelist') to the empty string ('').
+                        Default: 'cosmic_gps_input.nc'</TD></TR>
+   
+ <TR><!--contents--><TD valign=top>gpsro_netcdf_filelist</TD>
+     <!--  type  --><TD valign=top>character(len=128)</TD>
+     <!--descript--><TD>To convert a series of profiles in a single execution
+                        create a text file which contains each input file,
+                        in ascii, one filename per line.  Set this item to
+                        the name of that file, and set 'gpsro_netcdf_file' to
+                        the empty string ('').
+                        Default: 'cosmic_gps_input_list'</TD></TR>
+   
+ <TR><!--contents--><TD valign=top>gpsro_out_file</TD>
+     <!--  type  --><TD valign=top>character(len=128)</TD>
+     <!--descript--><TD>The output file to be created.  Note that to be
+                        compatible with earlier versions of this program, if
+                        this file already exists it will be read in and the
+                        new data will be inserted into that file.
+                        Default: 'obs_seq.gpsro'</TD></TR>
+   
+ <TR><!--contents--><TD valign=top>overwrite_time</TD>
+     <!--  type  --><TD valign=top>logical</TD>
+     <!--descript--><TD>This item is NOT in the standard namelist, but code
+                        exists in the program to support the functionality
+                        if the source is edited and this variable is added 
+                        to the namelist.  If set to
+                        .true., the output file will be created with all
+                        observations marked with the specified analysis time
+                        instead of the actual observation collection time.
+                        This is generally not what you want, but might be
+                        useful if you are binning observations by time to do 
+                        thinning or super-ob'ing data.
+                        Default: .false.</TD></TR>
+   
+ </TABLE>
+<P>
+</P>
+
+<!--==================================================================-->
+
 <A NAME="KnownBugs"></A>
 <BR><HR><BR>
 <H2>KNOWN BUGS</H2>
 <P>
+Some COSMIC files seem to have internal times which differ from the
+times encoded in the filenames by as much as 2-3 minutes.  If it is
+important to get all the observations within a particular time window
+files with filenames from a few minutes before and after the window 
+should be converted.
+Times really outside the window can be excluded either by setting
+the 'obs_window' namelist variable, or can be trimmed later with the
+<a href="../../obs_sequence/obs_sequence_tool.html">obs_sequence_tool</a>.
 </P>
 
 <!--==================================================================-->
-<!-- Descibe Future Plans.                                            -->
+<!-- Describe Future Plans.                                           -->
 <!--==================================================================-->
 
 <A NAME="FuturePlans"></A>

Modified: DART/trunk/observations/gps/work/convert_script.csh
===================================================================
--- DART/trunk/observations/gps/work/convert_script.csh	2009-08-10 19:13:07 UTC (rev 3996)
+++ DART/trunk/observations/gps/work/convert_script.csh	2009-08-10 21:33:27 UTC (rev 3997)
@@ -3,20 +3,26 @@
 # Main script:
 # generate multiple days of gps observations
 #
-# calls the cosmic_to_obsseq script with a date, the working
-# directory location, and whether to download the data automatically
-# from the cosmic web site (downloading data requires signing up 
-# for a username/password to access the site, and then setting 
-# the username and password here before running this script.)
+# calls the cosmic_to_obsseq script with 4 args:
 #
+#  - the date in YYYYMMDD format
+#  - the working directory location
+#  - whether to download the data automatically from the cosmic web site 
+#     (downloading data requires signing up for a username/password to 
+#     access the site, and then setting the username and password here 
+#     before running this script.)  set to no if the data has already been
+#     downloaded separately before now.
+#  - whether to delete the data automatically from the local disk after the
+#     conversion is done.
+#
 
 setenv cosmic_user xxx
 setenv cosmic_pw   yyy
 
-./cosmic_to_obsseq.csh 20061101 .. yes
-./cosmic_to_obsseq.csh 20061102 .. yes
-./cosmic_to_obsseq.csh 20061103 .. yes
-./cosmic_to_obsseq.csh 20061104 .. yes
-./cosmic_to_obsseq.csh 20061105 .. yes
-./cosmic_to_obsseq.csh 20061106 .. yes
-./cosmic_to_obsseq.csh 20061107 .. yes
+./cosmic_to_obsseq.csh 20061101 .. no no
+./cosmic_to_obsseq.csh 20061102 .. no no
+./cosmic_to_obsseq.csh 20061103 .. no no
+./cosmic_to_obsseq.csh 20061104 .. no no
+./cosmic_to_obsseq.csh 20061105 .. no no
+./cosmic_to_obsseq.csh 20061106 .. no no
+./cosmic_to_obsseq.csh 20061107 .. no no

Added: DART/trunk/observations/gps/work/cosmic_download.csh
===================================================================
--- DART/trunk/observations/gps/work/cosmic_download.csh	                        (rev 0)
+++ DART/trunk/observations/gps/work/cosmic_download.csh	2009-08-10 21:33:27 UTC (rev 3997)
@@ -0,0 +1,153 @@
+#!/bin/csh 
+########################################################################
+#
+#   cosmic_download.csh - script that downloads COSMIC observations.
+#         then you can run cosmic_to_obsseq.csh with a third argument of
+#         'no' so it does not re-download the same files.
+#
+# requires 2 args:
+#    $1 - analysis date (yyyymmdd format)
+#    $2 - base observation directory
+#
+# and update the 3 env settings at the top to match your system.
+#
+#     created May  2009, nancy collins ncar/cisl
+#                        converted from cosmic_to_obsseq.csh
+#     updated Aug  2009, nancy collins ncar/cisl
+#
+# from the cosmic web site about the use of 'wget' to download
+# the many files needed to do this process:
+# ------- 
+# Hints for using wget for fetching CDAAC files from CDAAC:
+# 
+# Here is one recipe for fetching all cosmic real time atmPrf files for one day:
+# 
+# wget -nd -np -r -l 10 -w 2 --http-user=xxxx --http-passwd=xxxx http://cosmic-io.cosmic.ucar.edu/cdaac/login/cosmicrt/level2/atmPrf/2009.007/
+# 
+# The option -np (no parents) is important. Without it, all manner of 
+# files from throughout the site will be loaded, I think due to the 
+# links back to the main page which are everywhere.
+# 
+# The option -r (recursive fetch) is necessary unless there is just 
+# one file you want to fetch.
+# 
+# The option -l 10 (limit of recursive depth to 10 levels) is necessary 
+# in order to get around the default 5 level depth.
+# 
+# The option -nd dumps all fetched files into your current directory. 
+# Otherwise a directory hierarchy will be created: 
+#   cosmic-io.cosmic.ucar.edu/cdaac/login/cosmic/level2/atmPrf/2006.207
+# 
+# The option -w 2 tells wget to wait two seconds between each file fetch 
+# so as to have pity on our poor web server.
+# ------- 
+# 
+########################################################################
+
+# should only have to set the DART_DIR and the rest should be in the
+# right place relative to it.
+setenv DART_DIR      /home/user/dart
+setenv cosmic_user   xxx
+setenv cosmic_pw     yyy
+
+setenv DART_WORK_DIR  ${DART_DIR}/observations/gps/work
+setenv CONV_PROG      convert_cosmic_gps_cdf
+setenv DATE_PROG      advance_time
+
+set chatty=yes
+set downld=yes
+
+set datea   = ${1}     # target date, YYYYMMDD
+set datadir = ${2}     # where to process the files
+
+set assim_freq          = 6  # hours, sets centers of windows.
+set download_window     = 3  # window half-width (some users choose 2 hours)
+set gps_repository_path = 'http://cosmic-io.cosmic.ucar.edu/cdaac/login'
+set wget_cmd            = 'wget -q -nd -np -r -l 10 -w 1'
+
+# i've done this wrong enough times and wasted a lot of download time,
+# so do a bunch of bullet-proofing here before going on.
+
+# verify the dirs all exist, the input.nml is in place.
+if ( ! -d ${DART_WORK_DIR} ) then
+  echo 'work directory not found: ' ${DART_WORK_DIR}
+  exit
+endif
+
+echo 'current dir is ' `pwd`
+if ( `pwd` != ${DART_WORK_DIR} ) then
+   echo 'if not already there, changing directory to work dir.'
+   cd ${DART_WORK_DIR}
+   echo 'current dir now ' `pwd`
+endif
+
+if ( ! -d ${datadir} ) then
+  echo 'data processing directory not found: ' ${datadir}
+  echo 'creating now.'
+  mkdir ${datadir}
+  ls -ld ${datadir}
+endif
+
+if ( ! -e ${datadir}/${DATE_PROG} ) then
+  echo 'data processing directory does not contain the time converter'
+  echo 'copying from work dir to data proc dir'
+  echo `pwd`/${DATE_PROG} '->' ${datadir}/${DATE_PROG}
+  cp -f ./${DATE_PROG} ${datadir}
+else
+  echo 'using time conversion program already found in data proc dir'
+endif
+
+echo 'changing dir to data proc directory'
+cd ${datadir}
+echo 'current dir now ' `pwd`
+
+ if ( ! $?cosmic_user || ! $?cosmic_pw ) then
+    echo "You must setenv cosmic_user to your username for the cosmic web site"
+    echo "and setenv cosmic_pw to your password. (or export cosmic_user=val for"
+    echo "ksh/bash users) then rerun this script. "
+    exit -1
+ endif
+
+if ( $chatty == 'yes' ) then
+   echo 'starting raw file download at' `date`
+endif
+
+set get = "${wget_cmd} --http-user=${cosmic_user} --http-passwd=${cosmic_pw}" 
+
+set yyyy   = `echo $datea | cut -b1-4`
+
+if ( ! -d ${datea} ) then
+  echo 'year/month/day directory not found: ' ${datea}
+  echo 'creating now.'
+  mkdir ${datea}
+endif
+
+cd ${datea}
+
+echo $datea
+set jyyyydd = `echo $datea 0 -j | ../${DATE_PROG}` 
+echo $jyyyydd
+@ mday = $jyyyydd[2] + 1000  ;  set mday = `echo $mday | cut -b2-4`
+${get} ${gps_repository_path}/cosmic/level2/atmPrf/${yyyy}.${mday}/
+rm -f *.html *.txt
+${get} ${gps_repository_path}/champ/level2/atmPrf/${yyyy}.${mday}/
+rm -f *.html *.txt
+
+set jyyyydd = `echo $datea 24 -j | ../${DATE_PROG}`
+@ mday = $jyyyydd[2] + 1000  ;  set mday = `echo $mday | cut -b2-4` 
+${get} ${gps_repository_path}/cosmic/level2/atmPrf/${yyyy}.${mday}/
+rm -f *.html *.txt
+${get} ${gps_repository_path}/champ/level2/atmPrf/${yyyy}.${mday}/
+rm -f *.html *.txt
+
+if ( $chatty == 'yes' ) then
+   # the ls arg list line gets too long in some cases
+   echo `/bin/ls | grep _nc | wc -l` 'raw files'
+   echo 'all raw files download at ' `date`
+endif
+
+cd ..
+
+exit 0
+
+


Property changes on: DART/trunk/observations/gps/work/cosmic_download.csh
___________________________________________________________________
Name: svn:executable
   + *

Modified: DART/trunk/observations/gps/work/cosmic_to_obsseq.csh
===================================================================

@@ Diff output truncated at 40000 characters. @@


More information about the Dart-dev mailing list