[Dart-dev] [3932] DART/trunk/observations/AIRS: Incudes new thinning namelist options, and the

nancy at ucar.edu nancy at ucar.edu
Fri Jun 19 09:27:26 MDT 2009


An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090619/ba6075de/attachment-0001.html 
-------------- next part --------------
Modified: DART/trunk/observations/AIRS/airs_obs_mod.f90
===================================================================
--- DART/trunk/observations/AIRS/airs_obs_mod.f90	2009-06-19 15:26:03 UTC (rev 3931)
+++ DART/trunk/observations/AIRS/airs_obs_mod.f90	2009-06-19 15:27:25 UTC (rev 3932)
@@ -13,20 +13,24 @@
 
 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 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, &
+                             get_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, &
+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(>=)
 
 use    utilities_mod, only : get_unit, open_file, close_file, file_exist, &
-                             register_module, error_handler, &
+                             register_module, error_handler,              &
                              E_ERR, E_MSG, timestamp, is_longitude_between
 
-use     location_mod, only : location_type, set_location, VERTISPRESSURE
+use     location_mod, only : location_type, set_location, VERTISPRESSURE, &
+                             get_location
 
 use     obs_kind_mod, only : get_obs_kind_index, &
                              KIND_TEMPERATURE, KIND_SPECIFIC_HUMIDITY
@@ -34,9 +38,10 @@
 use obs_kind_mod,     only : AIRS_TEMPERATURE, AIRS_SPECIFIC_HUMIDITY
 
 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
+                             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
 
 use airs_JPL_mod   ! need ', only' list here
 
@@ -56,22 +61,18 @@
 
 logical :: DEBUG = .false.
 
-! fixed size storage
-real ::   T    (AIRS_RET_STDPRESSURELEV, AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)
-real ::   T_err(AIRS_RET_STDPRESSURELEV, AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)
-real :: MMR    (AIRS_RET_H2OPRESSURELEV, AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)
-real :: MMR_err(AIRS_RET_H2OPRESSURELEV, AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)
-real ::   Q    (AIRS_RET_H2OPRESSURELEV, AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)
-real ::   Q_err(AIRS_RET_H2OPRESSURELEV, AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)
-real ::       PBest(                     AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)
-real ::       PGood(                     AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)
-integer :: nBestStd(                     AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)
-integer :: nGoodStd(                     AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)
-real ::         lat(                     AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)
-real ::         lon(                     AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)
-real ::         tim(                     AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)
-integer :: qual_h2o(                     AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)
+! the sizes of the Temperature arrays are:
+!   (AIRS_RET_STDPRESSURELAY, AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)
+! the sizes of the MMR arrays are:
+!   (AIRS_RET_H2OPRESSURELAY, AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)
+! the rest of the arrays are:
+!   (AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)
 
+! MMR is converted into Q:  Q = MMR / (1 + MMR)
+real ::   Q    (AIRS_RET_H2OPRESSURELAY, AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)
+real ::   Q_err(AIRS_RET_H2OPRESSURELAY, AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)
+
+
 contains
 
 
@@ -86,41 +87,10 @@
 end subroutine initialize_module
 
 
-subroutine create_output_filename(l2name, ofname)
-!-------------------------------------------------
-! The L2 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)  :: l2name
-character(len=*), intent(OUT) :: ofname
 
-integer :: i, basestart, extstart, strlen
-
-! hardcoded and brittle, but for now...  the first 19 chars
-! of the input filename have the date & granule number, which
-! seems like the bulk of the useful info.  find the last / and
-! copy from there to +19 chars.
-
-strlen = len_trim(l2name)
-
-basestart = 1
-slashloop : do i = strlen-1,1,-1
-   if (l2name(i:i) == '/' ) then
-      basestart = i+1
-      exit slashloop
-   endif
-enddo slashloop
-
-extstart = basestart+19-1
-
-ofname = l2name(basestart:extstart)//'.out'
-print *, 'output filename = ', ofname
-
-end subroutine create_output_filename
-
-
-
-function real_obs_sequence ( granule, lon1, lon2, lat1, lat2 )
+function real_obs_sequence ( granule, lon1, lon2, lat1, lat2, &
+                             min_MMR_threshold, top_pressure_level, &
+                             row_thin, col_thin)
 !------------------------------------------------------------------------------
 !  extract the temperature and humidity observations from a granule
 !  and convert to DART observation format.  allow caller to specify
@@ -128,24 +98,28 @@
 
 type(airs_granule_type), intent(in) :: granule
 real(r8), intent(in) :: lon1, lon2, lat1, lat2
+real(r8), intent(in) :: min_MMR_threshold, top_pressure_level
+integer,  intent(in) :: row_thin, col_thin
 
 ! max possible obs from this one granule.
-integer :: max_num=  AIRS_RET_STDPRESSURELEV * AIRS_RET_GEOXTRACK * AIRS_RET_GEOTRACK + &
-                     AIRS_RET_H2OPRESSURELEV * AIRS_RET_GEOXTRACK * AIRS_RET_GEOTRACK 
+integer :: max_num =  &
+   AIRS_RET_STDPRESSURELAY * AIRS_RET_GEOXTRACK * AIRS_RET_GEOTRACK + &
+   AIRS_RET_H2OPRESSURELAY * AIRS_RET_GEOXTRACK * AIRS_RET_GEOTRACK 
 
 type(obs_sequence_type) :: real_obs_sequence
 type(obs_def_type)      :: obs_def
 type(obs_type)          :: obs, prev_obs
+type(location_type)     :: obs_loc
 
 integer :: i, irow, icol, ivert, num_copies, num_qc
 integer :: days, seconds
-integer :: obs_num
+integer :: obs_num, temperature_top_index, humidity_top_index
 integer :: which_vert, tobstype, qobstype
 
 real(r4) :: speed, dir
 real(r8) :: olon, olat, vloc
 real(r8) :: obs_value, obs_var
-real(r8) :: tqc, qqc
+real(r8) :: tqc, qqc, latlon(3)
 real(r8) :: sintheta, costheta, dirvar, speedvar
 
 type(time_type) :: obs_time, base_time, pre_time, time
@@ -184,94 +158,73 @@
 endif
 
 
-
-!filename = '../data/200711/20071101/AIRS.2007.11.01.001.L2.RetStd.v5.2.2.0.G08078150655.hdf'
-
-
 base_time = set_date(1993, 1, 1, 0, 0, 0)   ! reference date: jan 1st, 1993
 
+call debug_print_size_check(granule)
 
-print *, 'size TairStd', size(granule%TairStd)
-T = reshape(granule%TairStd, shape(T))
-print *, 'first  row T', T(:, 1, 1)
-print *, 'second row T', T(:, 2, 1)
-print *, 'second col T', T(:, 1, 2)
+! The file contains 'Retrieved Water Vapor Mass Mixing Ratio'.  Convert
+! to specific humidity here.  Original units: gm/kg - scale to kg/kg first.
 
-print *, 'size TairStdErr', size(granule%TairStdErr)
-T_err = reshape(granule%TairStdErr, shape(T))
-print *, 'first  row T_err', T_err(:, 1, 1)
+where (granule%H2OMMRStd > min_MMR_threshold) 
+   Q = (granule%H2OMMRStd / 1000.) / ( 1.0 + (granule%H2OMMRStd/1000.0))
+elsewhere
+   Q = 0.0
+endwhere
 
-print *, 'PBest', size(granule%PBest)
-PBest = reshape(granule%PBest, shape(PBest))
-print *, 'first  row PBest', PBest(:, 1)
+where (granule%H2OMMRStdErr > 0.0_r8) 
+   Q_err = (granule%H2OMMRStdErr / 1000.) / ( 1.0 + (granule%H2OMMRStdErr/1000.0))
+elsewhere
+   Q_err = 0.0_r8   ! is this really ok?
+endwhere
 
-print *, 'PGood', size(granule%PGood)
-PGood = reshape(granule%PGood, shape(PGood))
-print *, 'first  row PGood', PGood(:, 1)
+! see what the top pressure allowed will be, and compute the index number
+! that is still below it in each profile.  that will be become the loop end.
+! default to doing the whole column
+temperature_top_index = size(granule%pressStd)
+do i = 1, size(granule%pressStd)
+   if (granule%pressStd(i) < top_pressure_level) then
+      temperature_top_index = i
+      exit
+   endif
+enddo
+if (DEBUG) print *, 'temp_top_index = ', temperature_top_index
 
-print *, 'nBestStd', size(granule%nBestStd)
-nBestStd = reshape(granule%nBestStd, shape(nBestStd))
-print *, 'first  row nBestStd', nBestStd(:, 1)
+humidity_top_index = size(granule%pressH2O)
+do i = 1, size(granule%pressH2O)
+   if (granule%pressH2O(i) < top_pressure_level) then
+      humidity_top_index = i
+      exit
+   endif
+enddo
+if (DEBUG) print *, 'humd_top_index = ', humidity_top_index
 
-print *, 'nGoodStd', size(granule%nGoodStd)
-nGoodStd = reshape(granule%nGoodStd, shape(nGoodStd))
-print *, 'first  row nGoodStd', nGoodStd(:, 1)
-
-print *, 'size water vapor mass mixing ratio', size(granule%H2OMMRStd)
-MMR = reshape(granule%H2OMMRStd, shape(MMR))
-print *, 'first  row MMR', MMR(:, 1, 1)
-print *, 'second col MMR', MMR(:, 2, 1)
-print *, 'second col MMR', MMR(:, 1, 2)
-
-print *, 'size water vapor mass mixing ratio err', size(granule%H2OMMRStdErr)
-MMR_err = reshape(granule%H2OMMRStdErr, shape(MMR_err))
-print *, 'first  row MMR_err', MMR_err(:, 1, 1)
-
-print *, 'water vapor mass mixing ratio quality', size(granule%Qual_H2o)
-qual_h2o = reshape(granule%Qual_H2O, shape(qual_h2o))
-print *, 'first  row qual_h2o', qual_h2o(:, 1)
-
-Q = (MMR / 1000.) / ( 1.0 + (MMR/1000.0))
-Q_err = (MMR_err / 1000.) / ( 1.0 + (MMR_err/1000.0))
-print *, 'first  row Q', Q(:, 1, 1)
-print *, 'second col Q', Q(:, 2, 1)
-print *, 'second col Q', Q(:, 1, 2)
-
-print *, 'first  row Q_err', Q_err(:, 1, 1)
-
-lat = reshape(granule%Latitude, shape(lat))
-lon = reshape(granule%Longitude, shape(lon))
-tim = reshape(granule%Time, shape(tim))
-
-print *, '1st lat, lon, time', lat(1,1), lon(1,1), tim(1,1)
-print *, 'r 2 lat, lon, time', lat(2,1), lon(2,1), tim(2,1)
-print *, 'c 2 lat, lon, time', lat(1,2), lon(1,2), tim(1,2)
-
-! compute time of this obs
-print *, 'int time ', int(tim(1, 1))
-obs_time = base_time + set_time(int(tim(1, 1)))
-call print_date(obs_time)
-print *, 'int time ', int(tim(AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK))
-obs_time = base_time + set_time(int(tim(AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)))
-call print_date(obs_time)
-
-print *, 'pressStd', granule%pressStd
-print *, 'pressH2O', granule%pressH2O
-
 !  loop over all observations within the file
 !------------------------------------------------------------------------------
 
 obs_num    = 1
 which_vert = VERTISPRESSURE
 
+! rows are along-track, stepping in the direction the satellite is moving
 rowloop:  do irow=1,AIRS_RET_GEOTRACK
 
+   ! if we're going to subset, we could cycle here
+   if (row_thin > 0) then
+      if (modulo(irow, row_thin) /= 1) cycle rowloop
+   endif
+
+   ! columns are across-track, varying faster than rows.
    colloop:  do icol=1,AIRS_RET_GEOXTRACK
 
-      olat  = lat(icol,irow) ! valid range [ -90.00,  90.00]
-      olon  = lon(icol,irow) ! valid range [-180.00, 180.00]
+      ! if we're going to subset, ditto
+      if (col_thin > 0) then
+         if (modulo(icol, col_thin) /= 1) cycle colloop
+      endif
 
-      ! verify the location is not outside valid limits
+      ! observation lat, lon:
+      olat  = granule%Latitude (icol,irow) ! valid range [ -90.00,  90.00]
+      olon  = granule%Longitude(icol,irow) ! valid range [-180.00, 180.00]
+
+      ! verify the location is not outside valid limits.  AIRS  uses -180/180
       if((olon > 180.0_r8) .or. (olon < -180.0_r8) .or.  &
          (olat >  90.0_r8) .or. (olat <  -90.0_r8)) then
          write(*,*)'WARNING : invalid location.  col,row,lon,lat = ', icol,irow,olon,olat
@@ -285,21 +238,28 @@
       ! make sure lon is between 0 and 360
       if (olon < 0) olon = olon + 360.0_r8
 
-      obs_time = base_time + set_time(int(tim(icol, irow)))
+      obs_time = base_time + set_time(int(granule%Time(icol, irow)))
       call get_time(obs_time, seconds, days)
      
 
-      vert_T_loop:  do ivert=nBestStd(icol, irow), AIRS_RET_STDPRESSURELEV
+      vert_T_loop: do ivert=granule%nBestStd(icol, irow), temperature_top_index
 
-         tqc = 0   ! FIXME
+         tqc = 0   ! if we get here, the quality control is 'best' == 0
 
          ! create the obs_def for this observation, add to sequence
          !---------------------------------------------------------------
       
+
+         ! apparently -9999 is missing data, outside of qc mechanism
+         obs_value = granule%TAirStd(ivert, icol, irow)
+         if (obs_value == -9999.0_r8) cycle vert_T_loop
+
+
+         obs_var = granule%TAirStdErr(ivert, icol, irow) * &
+                   granule%TAirStdErr(ivert, icol, irow)
+
          vloc = granule%pressStd(ivert)
-   
-         obs_value = T(ivert, icol, irow)
-         obs_var = T_err(ivert, icol, irow) * T_err(ivert, icol, irow)
+
          call real_obs(num_copies, num_qc, obs, olon, olat, vloc, obs_value, &
                        obs_var, tqc, AIRS_TEMPERATURE, which_vert, seconds, days)
       
@@ -323,19 +283,22 @@
  
       enddo vert_T_loop
 
-      vert_Q_loop:  do ivert=1,AIRS_RET_H2OPRESSURELEV
+      vert_Q_loop:  do ivert=1,humidity_top_index
 
-         if (qual_h2o(icol, irow) > 0) exit vert_Q_loop
+         if (granule%Qual_H2O(icol, irow) > 0) exit vert_Q_loop
 
-         qqc = 0   ! FIXME
+         qqc = 0   ! if we get here, the quality control is 'Best' == 0
 
          ! create the obs_def for this observation, add to sequence
          !---------------------------------------------------------------
-      
-         vloc = granule%pressH2O(ivert)
    
+         ! if original MMR data was -9999, that is apparently a missing val
+         if (granule%H2OMMRStd(ivert, irow, icol) == -9999.0_r8) cycle vert_Q_loop
+
          obs_value = Q(ivert, icol, irow)
          obs_var = Q_err(ivert, icol, irow) * Q_err(ivert, icol, irow)
+         vloc = granule%pressH2O(ivert)
+
          call real_obs(num_copies, num_qc, obs, olon, olat, vloc, obs_value, &
                        obs_var, qqc, AIRS_SPECIFIC_HUMIDITY, which_vert, seconds, days)
       
@@ -363,21 +326,30 @@
 enddo rowloop
 
 ! Print a little summary
-print*, 'obs used = ', obs_num, ' obs skipped = ', max_num - obs_num
+write(msgstring,*) 'obs used = ', obs_num, ' obs skipped = ', max_num - obs_num
+call error_handler(E_MSG, ' ', msgstring)
 
 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 ')
+   call print_date(pre_time,' which is gregorian date   ')
+   obs_loc = get_obs_def_location(obs_def)
+   latlon = get_location(obs_loc)
+   write(msgstring,*)'lat,lon,pres =', latlon(1), latlon(2), latlon(3)
+   call error_handler(E_MSG, 'location: ', msgstring)
 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 ')
+   call print_date(time,' which is gregorian date   ')
+   obs_loc = get_obs_def_location(obs_def)
+   latlon = get_location(obs_loc)
+   write(msgstring,*)'lat,lon,pres =', latlon(1), latlon(2), latlon(3)
+   call error_handler(E_MSG, 'location: ', msgstring)
 endif
-print*, ''
+call error_handler(E_MSG, ' ', ' ')
 
 end function real_obs_sequence
 
@@ -441,4 +413,126 @@
 end subroutine real_obs_def
 
 
+subroutine check_size(size1, size2, varlabel, subrlabel)
+!----------------------------------------------------------------------
+integer,          intent(in) :: size1, size2
+character(len=*), intent(in) :: varlabel, subrlabel
+
+if (size1 /= size2) then
+   write(msgstring, '(A,I6,A,I6)') 'unexpected size '//trim(varlabel)//': ', &
+         size1, ' /= ', size2
+   call error_handler(E_ERR, trim(subrlabel), msgstring,source,revision,revdate)
+endif
+
+end subroutine check_size
+
+
+subroutine create_output_filename(l2name, ofname)
+!-------------------------------------------------
+! The L2 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)  :: l2name
+character(len=*), intent(OUT) :: ofname
+
+integer :: i, basestart, extstart, strlen
+
+! hardcoded and brittle, but for now...  the first 19 chars
+! of the input filename have the date & granule number, which
+! seems like the bulk of the useful info.  find the last / and
+! copy from there to +19 chars.
+
+strlen = len_trim(l2name)
+
+basestart = 1
+slashloop : do i = strlen-1,1,-1
+   if (l2name(i:i) == '/' ) then
+      basestart = i+1
+      exit slashloop
+   endif
+enddo slashloop
+
+extstart = basestart+19-1
+
+ofname = l2name(basestart:extstart)//'.out'
+if (DEBUG) print *, 'output filename = ', ofname
+
+end subroutine create_output_filename
+
+
+subroutine debug_print_size_check(granule)
+!-------------------------------------------------
+! bit of sanity checking - before we loop over these arrays, make sure
+! they are the size we expect them to be.
+type(airs_granule_type), intent(in) :: granule
+
+type(time_type) :: obs_time, base_time
+
+base_time = set_date(1993, 1, 1, 0, 0, 0)   ! reference date: jan 1st, 1993
+
+call check_size(size(granule%TAirStd),                           &
+   AIRS_RET_STDPRESSURELAY*AIRS_RET_GEOXTRACK*AIRS_RET_GEOTRACK, &
+   'TAirStd (T)','real_obs_sequence')
+
+if (DEBUG) print *, 'first  row T', granule%TAirStd(:, 1, 1)
+if (DEBUG) print *, 'second row T', granule%TAirStd(:, 2, 1)
+if (DEBUG) print *, 'second col T', granule%TAirStd(:, 1, 2)
+
+call check_size(size(granule%TAirStdErr), size(granule%TAirStd), &
+   'TAirStdErr (T_err)','real_obs_sequence')
+
+if (DEBUG) print *, 'first  row T err', granule%TAirStdErr(:, 1, 1)
+
+! First (lowest) good pressure level number
+call check_size(size(granule%nBestStd), AIRS_RET_GEOXTRACK*AIRS_RET_GEOTRACK, &
+   'nBestStd (T QC)','real_obs_sequence')
+
+if (DEBUG) print *, 'first  row nBestStd', granule%nBestStd(:, 1)
+
+call check_size(size(granule%H2OMMRStd),                         &
+   AIRS_RET_H2OPRESSURELAY*AIRS_RET_GEOXTRACK*AIRS_RET_GEOTRACK, &
+   'H2OMMRStd (MMR)','real_obs_sequence')
+
+if (DEBUG) print *, 'first  row MMR', granule%H2OMMRStd(:, 1, 1)
+if (DEBUG) print *, 'second col MMR', granule%H2OMMRStd(:, 2, 1)
+if (DEBUG) print *, 'second col MMR', granule%H2OMMRStd(:, 1, 2)
+
+call check_size(size(granule%H2OMMRStdErr), size(granule%H2OMMRStd), &
+   'H2OMMRStdErr (MMR_err)','real_obs_sequence')
+
+if (DEBUG) print *, 'first  row MMR err', granule%H2OMMRStdErr(:, 1, 1)
+
+call check_size(size(granule%Qual_H2O), AIRS_RET_GEOXTRACK*AIRS_RET_GEOTRACK, &
+   'Qual_H2O (Q QC)','real_obs_sequence')
+
+if (DEBUG) print *, 'first  row qual_h2o', granule%Qual_H2O(:, 1)
+
+if (DEBUG) print *, 'first  row Q', Q(:, 1, 1)
+if (DEBUG) print *, 'second col Q', Q(:, 2, 1)
+if (DEBUG) print *, 'second col Q', Q(:, 1, 2)
+
+if (DEBUG) print *, 'first  row Q_err', Q_err(:, 1, 1)
+
+
+if (DEBUG) print *, '1st lat, lon, time', &
+           granule%Latitude(1,1), granule%Longitude(1,1), granule%Time(1,1)
+if (DEBUG) print *, 'r 2 lat, lon, time', &
+           granule%Latitude(2,1), granule%Longitude(2,1), granule%Time(2,1)
+if (DEBUG) print *, 'c 2 lat, lon, time', &
+           granule%Latitude(1,2), granule%Longitude(1,2), granule%Time(1,2)
+
+! debugging for first/last time of obs in this file
+if (DEBUG) print *, 'int start time ', int(granule%Time(1, 1))
+if (DEBUG) obs_time = base_time + set_time(int(granule%Time(1, 1)))
+if (DEBUG) call print_date(obs_time)
+
+if (DEBUG) print *, 'int end   time ', int(granule%Time(AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK))
+if (DEBUG) obs_time = base_time + set_time(int(granule%Time(AIRS_RET_GEOXTRACK, AIRS_RET_GEOTRACK)))
+if (DEBUG) call print_date(obs_time)
+
+if (DEBUG) print *, 'pressStd', granule%pressStd
+if (DEBUG) print *, 'pressH2O', granule%pressH2O
+
+end subroutine debug_print_size_check
+
 end module airs_obs_mod

Modified: DART/trunk/observations/AIRS/convert_airs_L2.f90
===================================================================
--- DART/trunk/observations/AIRS/convert_airs_L2.f90	2009-06-19 15:26:03 UTC (rev 3931)
+++ DART/trunk/observations/AIRS/convert_airs_L2.f90	2009-06-19 15:27:25 UTC (rev 3932)
@@ -22,9 +22,9 @@
                              error_handler, timestamp, E_ERR, E_MSG, &
                              find_namelist_in_file, check_namelist_read, &
                              do_nml_file, do_nml_term, &
-                             logfileunit, nmlfileunit
+                             logfileunit, nmlfileunit, get_next_filename
 
-use airs_obs_mod   ! FIXME: need to add ,only :   
+use airs_obs_mod,     only : real_obs_sequence, create_output_filename
 
 implicit none
 
@@ -32,11 +32,11 @@
 ! Declare local parameters
 ! ----------------------------------------------------------------------
 
-character(len=256)      :: datafile(1), output_name, dartfile, string1
+character(len=256)      :: datafile(1), output_name, dartfile, msgstring
 type(airs_granule_type) :: granule
 type(obs_sequence_type) :: seq
 
-integer :: io, iunit
+integer :: io, iunit, f, nfiles, index
 
 ! version controlled file description for error handling, do not edit
 character(len=128), parameter :: &
@@ -48,8 +48,12 @@
 ! Declare namelist parameters
 ! ----------------------------------------------------------------------
         
-character(len=128) ::   l2_file = ''
-character(len=128) ::   datadir = '.'
+integer, parameter :: MAXFILES = 256
+character(len=128) :: nextfile
+
+character(len=128) :: l2_files(MAXFILES) = ''
+character(len=128) :: l2_file_list       = ''
+character(len=128) :: datadir   = '.'
 character(len=128) :: outputdir = '.'
 
 real(r8) :: lon1 =   0.0_r8,  &   !  lower longitude bound
@@ -57,9 +61,17 @@
             lat1 = -90.0_r8,  &   !  lower latitude bound
             lat2 =  90.0_r8       !  upper latitude bound
 
-namelist /convert_airs_L2_nml/ l2_file, datadir, outputdir, &
-                           lon1, lon2, lat1, lat2
+real(r8) :: min_MMR_threshold = 1.0e-30
+real(r8) :: top_pressure_level = 0.0001    ! no obs higher than this
+integer  :: cross_track_thin = 0
+integer  :: along_track_thin = 0
 
+namelist /convert_airs_L2_nml/ l2_files, l2_file_list, &
+                               datadir, outputdir, &
+                               lon1, lon2, lat1, lat2, &
+                               min_MMR_threshold, top_pressure_level, &
+                               cross_track_thin, along_track_thin
+
 ! ----------------------------------------------------------------------
 ! start of executable program code
 ! ----------------------------------------------------------------------
@@ -83,15 +95,54 @@
 if (do_nml_file()) write(nmlfileunit, nml=convert_airs_L2_nml)
 if (do_nml_term()) write(    *      , nml=convert_airs_L2_nml)
 
-call create_output_filename(l2_file, output_name)
-datafile(1) = trim(datadir) // '/' // trim(l2_file)
-dartfile = trim(outputdir) // '/' // trim(output_name)
-!dartfile = trim(outputdir) // '/test.out'
+if ((l2_files(1) /= '') .and. (l2_file_list /= '')) then
+   write(msgstring,*)'cannot specify both an input file and an input file list'
+   call error_handler(E_ERR, 'convert_airs_L2', msgstring, &
+                      source, revision, revdate)
+endif
 
-call airs_ret_rdr(datafile, granule)   ! read from HDF file into a structure
-seq = real_obs_sequence(granule, 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
+
+index = 0
+
+! do loop without an index.  will loop until exit called.
+do
+   index = index + 1
+   if (l2_files(1) /= '') then
+      if (index > size(l2_files)) then
+         write(msgstring,*)'cannot specify more than ', size(l2_files), ' files'
+         call error_handler(E_ERR, 'convert_airs_L2', msgstring, &
+                            source, revision, revdate)
+      endif
+      nextfile = l2_files(index)
+   else
+      ! this is the new routine
+      ! it opens the listfile, returns the index-th one
+      nextfile = get_next_filename(l2_file_list, index)
+   endif
+
+   if (nextfile == '') exit
+
+   ! construct an appropriate output filename
+   call create_output_filename(nextfile, output_name)
+   datafile(1) = trim(datadir)   // '/' // trim(nextfile)
+   dartfile    = trim(outputdir) // '/' // trim(output_name)
+   
+   ! read from HDF file into a derived type that holds all the information
+   call airs_ret_rdr(datafile, granule)   
+
+   ! convert derived type information to DART sequence
+   seq = real_obs_sequence(granule, lon1, lon2, lat1, lat2, &
+                           min_MMR_threshold, top_pressure_level, &
+                           along_track_thin, cross_track_thin) 
+
+   ! write the sequence to a disk file
+   call write_obs_seq(seq, dartfile) 
+ 
+   ! release the sequence memory
+   call destroy_obs_sequence(seq)
+
+enddo
+
 call timestamp(source,revision,revdate,'end') ! close the log file
 
 end program convert_airs_L2


More information about the Dart-dev mailing list