!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This program decodes a swath of AIRS L2 standard retrievals, computing RH 
!    from q and writing the resulting soundings into a little_r formatted file 
!    named soundings.little_r. The code is based on that from the example reader
!    programs supplied through the AIRS web site.
!
! Michael Duda -- 26 May 2006
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
program decode_airs

   use read_airs

   implicit none

   real, external :: calc_rh
   integer, external :: iargc
   logical, external :: granule_out_of_timewin

   character (len=69), parameter :: rpt_format  = '(2f20.5,4a40,f20.5,5i10,3L10,2i10,a20,13(f13.5,i7))'
   character (len=15), parameter :: meas_format = '(10(f13.5, i7))'
   character (len=5), parameter  :: end_format  = '(3i7)'
   integer, parameter :: LESS = -1
   integer, parameter :: EQUAL = 0
   integer, parameter :: MORE = 1

   ! Data stored in the header record
   real :: latitude, longitude, elevation
   real :: slp, psfc, refp, grndt, sst, x1, x2, x3, x4, x5, x6, x7, x8
   integer :: islp, ipsfc, n_valid, n_err, n_warning, iseq, n_dups 
   integer :: isut, julian, irefp, igrndt, isst
   integer :: ix1, ix2, ix3, ix4, ix5, ix6, ix7, ix8
   character (len=20) :: date_char
   character (len=40) :: id, name, platform, source
   logical :: is_sound, is_bogus, is_discard

   ! Data stored in the data record
   real :: p, z, t, dewpt, spd, dir, u, v, rh, thk
   integer :: ip, iz, it, idewpt, ispd, idir, iu, iv, irh, ithk

   ! Data stored in the end record
   integer :: num_valid, num_err, num_warn

   type (airs_ret_gran_t) :: ret     ! structure with entire granule 
   integer :: layer                  ! index for cloud or atmospheric layers
   integer :: nargs                  ! number of arguments
   integer :: track                  ! along-track index (scan number)
   integer :: xtrack                 ! across-track index (FOV number)
   character (len=256) :: arg        ! buffer for argument
   character (len=256) :: file_name  ! name of AIRS Level 2 file to read
   character (len=11) :: ref_time
   character (len=2)  :: version

   character (len=14) :: min_time, max_time
   namelist /time_window/ min_time, max_time
   integer :: qual_threshold
   namelist /quality/ qual_threshold

   integer :: delta_time
   integer :: istatus
   integer :: nvalid, nvalid_temp, nvalid_h2o, npolarday_temp, npolarday_h2o
   integer :: fn
   integer :: cmp_min, cmp_max
   integer :: mm, ss
   integer :: qual                ! quality.  0 is best, then 1.  2 is bad.
   integer :: i, j                ! indices 1..3 for 3x3 of IR FOVs per retrieval
   real :: temp                   ! temperature or -9999 if not good enough
   real :: h2o                    ! H2O MMR or -9999 if not good enough
   logical :: polarday            ! true if |lat|>50 and DayNightFlag is not Night

   nargs = iargc()

   if (nargs <= 0) then
     write(6,*) ' '
     write(6,*) 'Usage: decode_airs.exe filename ...'
     write(6,*) '    where filename is the name of an HDFEOS file containing'
     write(6,*) '    a swath of L2 AIRS retreivals.'
     write(6,*) ' '
     stop
   end if


   min_time = '00000000000000'
   max_time = '99999999999999'

   open(11,file='time_window.nl',status='old',iostat=istatus)
   if (istatus == 0) then
      read(11,time_window)
      close(11)
   end if

   qual_threshold = 1   ! 0: best, 1:good, 2:not use

   open(12,file='qual_threshold.nl',status='old',iostat=istatus)
   if (istatus == 0) then
      read(12,quality)
      close(12)
   end if

   ! The time for each FOV is given as the number of seconds
   !   since 1993 Jan 1 00Z
   ref_time = '1993010100 '

   open(10, file='soundings.little_r', form='formatted', status='replace')
   open(20, file='soundings_polarday.little_r', form='formatted', status='replace')

   do fn = 1, nargs
      call getarg (fn, file_name)
   
      ! Read all retrievals from file_name and place them in ret
      call airs_ret_rdr(file_name, ret, version)

      if (granule_out_of_timewin(ret%Time, min_time, max_time)) cycle      
   
      nvalid_temp = 0
      nvalid_h2o  = 0
      npolarday_temp = 0
      npolarday_h2o  = 0

      write(6,'(a,a)') 'DayNightFlag for this granule: ', trim(ret%DayNightFlag)

      do track = 1, AIRS_RET_GEOTRACK
      do xtrack = 1, AIRS_RET_GEOXTRACK

         polarday = .false.
         nvalid = 0
      
         delta_time = nint(ret%Time(xtrack,track))/3600
         date_char(1:6) = '      '
         call geth_newdate(ref_time, delta_time, date_char(7:20))
         mm = mod(nint(ret%Time(xtrack,track)),3600)/60
         ss = mod(nint(ret%Time(xtrack,track)),60)
         write(date_char(17:20),'(i2.2,i2.2)') mm, ss
   
         if (ret%Time(xtrack,track) < 0.) then
            cmp_min = LESS
         else
            call cmp_datestr(date_char(7:20), min_time, cmp_min)
            call cmp_datestr(date_char(7:20), max_time, cmp_max)
         end if


         if (cmp_min /= LESS .and. cmp_max /= MORE) then

            write (6,'(a,f8.4,a5,f9.4,a3)') 'Location: ',                          &
                                            ret%Latitude(xtrack,track),'lat, ',    &
                                            ret%Longitude(xtrack,track),'lon'
            write(6,'(a)') 'Time: '//date_char(7:20)
            ! write(6,'(a,i9)') 'RetQAFlag: ', ret%RetQAFlag(xtrack,track)
            write(6,*) ' '
         
            !
            ! First loop through all levels and find out how many valid obs there are
            !
            if ( version == 'V4' ) then
               do layer=AIRS_RET_STDPRESSURELAY,1,-1
         
                  ! Which Qual flag to check depends on layer
                  if (ret%PressStd(layer) < ret%Press_mid_top_bndry(xtrack,track) ) then
                     qual = ret%Qual_Temp_Profile_Top(xtrack,track)
      
                  else if (ret%PressStd(layer) < ret%Press_bot_mid_bndry(xtrack,track) ) then
                     qual = ret%Qual_Temp_Profile_Mid(xtrack,track)
      
                  else 
                     qual = ret%Qual_Temp_Profile_Bot(xtrack,track)
      
                  endif
               
                  ! Temperature. If quality is bad (2) then put out flag value of -9999.0
                  if (qual <= qual_threshold .and. layer > ret%nSurfStd(xtrack,track)) then
                     temp = ret%TAirStd(layer,xtrack,track)
                     h2o = ret%H2OMMRStd(layer,xtrack,track)

                  ! See p.3 of Level 2 Product Levels and Layers
                  ! else if (qual <= qual_threshold .and. abs(ret%PressStd(layer) - ret%PSurfStd(xtrack,track)) < 5.) then
!                 ! temp = ret%TAirStd(layer,xtrack,track)
!                 ! h2o = ret%H2OMMRStd(layer,xtrack,track)

                  else
                     temp = -9999.0
                     h2o = -9999.0
      
                  end if
          
                  if (temp /= -9999.) then
                     nvalid = nvalid + 1
                     if (h2o /= -9999.) nvalid = nvalid + 1
                  end if
      
               end do
            end if  ! end of V4 quality selection

            if ( version == 'V5' ) then
               if ( qual_threshold == 0 ) then       ! Best
                  nvalid = AIRS_RET_STDPRESSURELAY - ret%nBestStd(xtrack,track) + 1
               else if ( qual_threshold == 1 ) then  ! Good and Best
                  if ( ret%nGoodStd(xtrack,track) == 29 ) then
                     nvalid = AIRS_RET_STDPRESSURELAY - ret%nBestStd(xtrack,track) + 1
                  else
                     nvalid = AIRS_RET_STDPRESSURELAY - ret%nGoodStd(xtrack,track) + 1
                  end if
               end if
               if ( (abs(ret%Latitude(xtrack,track)) > 50.0) .and. (trim(ret%DayNightFlag) /= 'Night') ) then
                  polarday = .true.
               end if
            end if  ! end of V5 quality selection
            !
            ! If there are valid obs to be reported, write them out to little_r format
            !
            if (nvalid > 0) then
         
               latitude = ret%Latitude(xtrack,track)
               longitude = ret%Longitude(xtrack,track)
               if ( version == 'V5' ) then
                  id       = 'AIRS L2 Std Retrieval '//version//'                '
               else
                  id       = 'AIRS L2 Std Retrieval                   '
               end if
               name     = 'AIRSRET                                 '
               platform = 'FM-133                                  '
               source   = 'GES DAAC                                '
               elevation = -888888.
               n_valid = nvalid
               n_err = 0
               n_warning = 0
               iseq = 0
               n_dups = 0
               is_sound = .true.
               is_bogus = .false.
               is_discard = .false.
               isut = -888888
               julian = -888888
               slp = -888888.
               islp = -88 
               refp = -888888.
               irefp = -88 
               grndt = -888888.
               igrndt = -88 
               sst = -888888.
               isst = -88 
               psfc = -888888.
               ipsfc = -88 
               x1 = -888888.
               ix1 = -88 
               x2 = -888888.
               ix2 = -88 
               x3 = -888888.
               ix3 = -88 
               x4 = -888888.
               ix4 = -88 
               x5 = -888888.
               ix5 = -88 
               x6 = -888888.
               ix6 = -88 
               x7 = -888888.
               ix7 = -88 
               x8 = -888888.
               ix8 = -88 
         
               if ( .not. polarday ) then
                  write(10, fmt=rpt_format) latitude, longitude, id, name, platform, &
                     source, elevation, n_valid, n_err, n_warning, iseq, n_dups,     &
                     is_sound, is_bogus, is_discard, isut, julian, date_char,        &
                     slp, islp, refp, irefp, grndt, igrndt, sst, isst, psfc, ipsfc,  &
                     x1, ix1, x2, ix2, x3, ix3, x4, ix4, x5, ix5, x6, ix6, x7, ix7,  &
                     x8, ix8
               else
                  write(20, fmt=rpt_format) latitude, longitude, id, name, platform, &
                     source, elevation, n_valid, n_err, n_warning, iseq, n_dups,     &
                     is_sound, is_bogus, is_discard, isut, julian, date_char,        &
                     slp, islp, refp, irefp, grndt, igrndt, sst, isst, psfc, ipsfc,  &
                     x1, ix1, x2, ix2, x3, ix3, x4, ix4, x5, ix5, x6, ix6, x7, ix7,  &
                     x8, ix8
               end if
         
               do layer=AIRS_RET_STDPRESSURELAY,1,-1

                  if ( version == 'V4' ) then
                     ! Which Qual flag to check depends on layer
                     if (ret%PressStd(layer) < ret%Press_mid_top_bndry(xtrack,track) ) then
                        qual = ret%Qual_Temp_Profile_Top(xtrack,track)
      
                     else if (ret%PressStd(layer) < ret%Press_bot_mid_bndry(xtrack,track) ) then
                        qual = ret%Qual_Temp_Profile_Mid(xtrack,track)
      
                     else 
                        qual = ret%Qual_Temp_Profile_Bot(xtrack,track)
      
                     endif
                  
                     ! Temperature. If quality is bad (2) then put out flag value of -888888.0
                     if (qual <= qual_threshold .and. layer > ret%nSurfStd(xtrack,track)) then
                        temp = ret%TAirStd(layer,xtrack,track)
                        if (temp == -9999.) temp = -888888.
                        h2o = ret%H2OMMRStd(layer,xtrack,track)
                        if (h2o == -9999.) h2o = -888888.
      
                     else
                        temp = -888888.0
                        h2o = -888888.0
      
                     end if
                  end if    ! end if V4 check

                  if ( version == 'V5' ) then

                     ! for temperature

                     temp = -888888.0
                     if ( qual_threshold == 0 ) then
                        if ( layer >= ret%nBestStd(xtrack,track) ) then
                           temp = ret%TAirStd(layer,xtrack,track)
                        end if
                     else if ( qual_threshold == 1 ) then
                        if ( layer >= ret%nGoodStd(xtrack,track) .or. layer >= ret%nBestStd(xtrack,track) ) then
                           temp = ret%TAirStd(layer,xtrack,track)
                        end if
                     else
                        temp = -888888.0
                     end if
                     if (temp == -9999.) temp = -888888.

                     ! for moisture

                     h2o = -888888.0
                     if ( qual_threshold == 0 ) then
                        if ( (layer <= AIRS_RET_H2OPRESSURELAY) .and.                &
                             (ret%Qual_H2O(xtrack,track) <= 1 ) .and.    &
                             (ret%pressH2O(layer) <= ret%PBest(xtrack,track)) ) then
                           h2o = ret%H2OMMRStd(layer,xtrack,track)
                        end if
                     else if ( qual_threshold == 1 ) then
                        if ( ret%Qual_H2O(xtrack,track) <= qual_threshold ) then
                           h2o = ret%H2OMMRStd(layer,xtrack,track)
                        end if
                     end if
                     if ( (ret%H2OMMRStdErr(layer,xtrack,track) < 0.0) .or.   &
                          (ret%H2OMMRStdErr(layer,xtrack,track) > 0.5*h2o)  ) then
                        h2o = -9999.0
                     end if
                     if (h2o == -9999.) h2o = -888888.

                  end if    ! end if V5 check

                  if ( temp > 0.0 ) nvalid_temp = nvalid_temp + 1
                  if ( h2o  > 0.0 ) nvalid_h2o  = nvalid_h2o  + 1

                  if ( polarday ) then
                     if ( temp > 0.0 ) then
                        nvalid_temp    = nvalid_temp - 1
                        npolarday_temp = npolarday_temp + 1
                     end if
                     if ( h2o  > 0.0 ) then
                        nvalid_h2o     = nvalid_h2o  - 1
                        npolarday_h2o  = npolarday_h2o  + 1
                     end if
                  end if

                  if (temp /= -888888.) then
             
                     p = ret%pressStd(layer)*100.
                     ip = 0
                     z = ret%GP_Height(layer,xtrack,track)
                     iz = 0
                     t = temp
                     it = 0
                     dewpt = -888888.
                     idewpt = 0
                     spd = -888888.
                     ispd = 0
                     dir = -888888.
                     idir = 0
                     u = -888888.
                     iu = 0
                     v = -888888.
                     iv = 0
                     if (h2o == -888888. .or. t == -888888.) then
                        rh = -888888.
                     else
                        rh = calc_rh(t, h2o/1000., p)
                        if (rh < 0.) rh = 0.
                     end if
                     irh = 0
                     thk = -888888.
                     ithk = 0
                     if ( .not. polarday ) then
                        write(10, fmt=meas_format) p, ip, z, iz, t, it, dewpt, idewpt, &
                              spd, ispd, dir, idir, u, iu, v, iv, rh, irh, thk, ithk
                     else
                        write(20, fmt=meas_format) p, ip, z, iz, t, it, dewpt, idewpt, &
                              spd, ispd, dir, idir, u, iu, v, iv, rh, irh, thk, ithk
                     end if
                  end if
               end do
            
               p = -777777.
               ip = 0
               z = -777777.
               iz = 0
               t = -888888.
               it = 0
               dewpt = -888888.
               idewpt = 0
               spd = -888888.
               ispd = 0
               dir = -888888.
               idir = 0
               u = -888888.
               iu = 0
               v = -888888.
               iv = 0
               rh = -888888.
               irh = 0
               thk = -888888.
               ithk = 0
      
               if ( .not. polarday ) then
                  write(10, fmt=meas_format) p, ip, z, iz, t, it, dewpt, idewpt, &
                        spd, ispd, dir, idir, u, iu, v, iv, rh, irh, thk, ithk
               else
                  write(20, fmt=meas_format) p, ip, z, iz, t, it, dewpt, idewpt, &
                        spd, ispd, dir, idir, u, iu, v, iv, rh, irh, thk, ithk
               end if
      
               num_valid = nvalid
               num_err = 0
               num_warn = 0
               if ( .not. polarday ) then
                  write(10, fmt=end_format) num_valid, num_err, num_warn
               else
                  write(20, fmt=end_format) num_valid, num_err, num_warn
               end if
         
            end if  ! end if nvalid > 0
         
            ! Need to output quality info so users can tell if we have a good
            ! retrieval with no clouds or a bad retrieval.  Both cases will put
            ! all -9999s below.
!            write(6,*) '# Cloud quality:'
!            if (ret%Qual_Cloud_OLR(xtrack,track) .eq. 0) then
!            write(6,*) '    0    Very Good'
!            else if (ret%Qual_Cloud_OLR(xtrack,track) .eq. 1) then
!               write(6,*) '    1    Good'
!            else
!               write(6,*) '    2    Do Not Use'
!            end if

         end if ! Profile falls within time window

      end do
      end do

      write(6,*) 'nvalid_temp, nvalid_h2o ', nvalid_temp, nvalid_h2o
      write(6,*) 'npolarday_temp, npolarday_h2o ', npolarday_temp, npolarday_h2o

   end do ! Loop over filenames

   close(10)
   close(20)

   stop
    
end program decode_airs

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Given two date strings of the form YYYYMMDDHHmmss, icmp is set to 
!    0 if the two dates are the same, 
!   -1 if date1 < date2, and 
!   +1 if date1 > date2.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine cmp_datestr(date1, date2, icmp)

   implicit none

   ! Arguments
   character (len=14), intent(in) :: date1, date2
   integer, intent(out) :: icmp

   ! Local variables
   integer :: yyyy1, yyyy2
   integer :: mm1, mm2
   integer :: dd1, dd2
   integer :: hh1, hh2
   integer :: min1, min2
   integer :: sec1, sec2

   read(date1(1:4), '(i4)') yyyy1
   read(date1(5:6), '(i2)') mm1
   read(date1(7:8), '(i2)') dd1
   read(date1(9:10), '(i2)') hh1
   read(date1(11:12), '(i2)') min1
   read(date1(13:14), '(i2)') sec1

   read(date2(1:4), '(i4)') yyyy2
   read(date2(5:6), '(i2)') mm2
   read(date2(7:8), '(i2)') dd2
   read(date2(9:10), '(i2)') hh2
   read(date2(11:12), '(i2)') min2
   read(date2(13:14), '(i2)') sec2

   if (yyyy1 > yyyy2) then
      icmp = 1
      return
   else if (yyyy1 < yyyy2) then
      icmp = -1
      return
   end if

   if (mm1 > mm2) then
      icmp = 1
      return
   else if (mm1 < mm2) then
      icmp = -1
      return
   end if

   if (dd1 > dd2) then
      icmp = 1
      return
   else if (dd1 < dd2) then
      icmp = -1
      return
   end if

   if (hh1 > hh2) then
      icmp = 1
      return
   else if (hh1 < hh2) then
      icmp = -1
      return
   end if

   if (min1 > min2) then
      icmp = 1
      return
   else if (min1 < min2) then
      icmp = -1
      return
   end if

   if (sec1 > sec2) then
      icmp = 1
      return
   else if (sec1 < sec2) then
      icmp = -1
      return
   end if

   icmp = 0
   return

end subroutine cmp_datestr


function granule_out_of_timewin(time_arr, min_time, max_time)

   use read_airs

   implicit none

   ! Parameters
   integer, parameter :: LESS = -1
   integer, parameter :: EQUAL = 0
   integer, parameter :: MORE = 1

   ! Arguments
   double precision, dimension(AIRS_RET_GEOXTRACK,AIRS_RET_GEOTRACK), intent(in) :: time_arr 
   character (len=14), intent(in) :: min_time, max_time

   ! Local variables
   integer :: i, j
   integer :: delta_time, cmp_min, cmp_max, mm, ss
   character (len=11) :: ref_time
   character (len=14) :: date_char

   ! Return value
   logical :: granule_out_of_timewin

   granule_out_of_timewin = .true.
   ref_time = '1993010100 '

   do i=1,AIRS_RET_GEOXTRACK,AIRS_RET_GEOXTRACK-1
   do j=1,AIRS_RET_GEOTRACK,AIRS_RET_GEOTRACK-1

      delta_time = nint(time_arr(i,j))/3600
      call geth_newdate(ref_time, delta_time, date_char)
      mm = mod(nint(time_arr(i,j)),3600)/60
      ss = mod(nint(time_arr(i,j)),60)
      write(date_char(11:14),'(i2.2,i2.2)') mm, ss
   
      if (time_arr(i,j) < 0.) then
         cmp_min = LESS
      else
         call cmp_datestr(date_char, min_time, cmp_min)
         call cmp_datestr(date_char, max_time, cmp_max)
      end if

      if (cmp_min /= LESS .and. cmp_max /= MORE) then
         granule_out_of_timewin = .false.
         return
      end if

   end do
   end do

end function granule_out_of_timewin