var/convertor/decode_l2_airs/decode_airs.f90

   1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   2 ! This program decodes a swath of AIRS L2 standard retrievals, computing RH
   3 !    from q and writing the resulting soundings into a little_r formatted file
   4 !    named soundings.little_r. The code is based on that from the example reader
   5 !    programs supplied through the AIRS web site.
   6 !
   7 ! Michael Duda -- 26 May 2006
   8 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   9 program decode_airs
  10
  11    use read_airs
  12
  13    implicit none
  14
  15    real, external :: calc_rh
  16    integer, external :: iargc
  17    logical, external :: granule_out_of_timewin
  18
  19    character (len=69), parameter :: rpt_format  = '(2f20.5,4a40,f20.5,5i10,3L10,2i10,a20,13(f13.5,i7))'
  20    character (len=15), parameter :: meas_format = '(10(f13.5, i7))'
  21    character (len=5), parameter  :: end_format  = '(3i7)'
  22    integer, parameter :: LESS = -1
  23    integer, parameter :: EQUAL = 0
  24    integer, parameter :: MORE = 1
  25
  26    ! Data stored in the header record
  27    real :: latitude, longitude, elevation
  28    real :: slp, psfc, refp, grndt, sst, x1, x2, x3, x4, x5, x6, x7, x8
  29    integer :: islp, ipsfc, n_valid, n_err, n_warning, iseq, n_dups
  30    integer :: isut, julian, irefp, igrndt, isst
  31    integer :: ix1, ix2, ix3, ix4, ix5, ix6, ix7, ix8
  32    character (len=20) :: date_char
  33    character (len=40) :: id, name, platform, source
  34    logical :: is_sound, is_bogus, is_discard
  35
  36    ! Data stored in the data record
  37    real :: p, z, t, dewpt, spd, dir, u, v, rh, thk
  38    integer :: ip, iz, it, idewpt, ispd, idir, iu, iv, irh, ithk
  39
  40    ! Data stored in the end record
  41    integer :: num_valid, num_err, num_warn
  42
  43    type (airs_ret_gran_t) :: ret     ! structure with entire granule
  44    integer :: layer                  ! index for cloud or atmospheric layers
  45    integer :: nargs                  ! number of arguments
  46    integer :: track                  ! along-track index (scan number)
  47    integer :: xtrack                 ! across-track index (FOV number)
  48    character (len=256) :: arg        ! buffer for argument
  49    character (len=256) :: file_name  ! name of AIRS Level 2 file to read
  50    character (len=11) :: ref_time
  51    character (len=2)  :: version
  52
  53    character (len=14) :: min_time, max_time
  54    namelist /time_window/ min_time, max_time
  55    integer :: qual_threshold
  56    namelist /quality/ qual_threshold
  57
  58    integer :: delta_time
  59    integer :: istatus
  60    integer :: nvalid, nvalid_temp, nvalid_h2o, npolarday_temp, npolarday_h2o
  61    integer :: fn
  62    integer :: cmp_min, cmp_max
  63    integer :: mm, ss
  64    integer :: qual                ! quality.  0 is best, then 1.  2 is bad.
  65    integer :: i, j                ! indices 1..3 for 3x3 of IR FOVs per retrieval
  66    real :: temp                   ! temperature or -9999 if not good enough
  67    real :: h2o                    ! H2O MMR or -9999 if not good enough
  68    logical :: polarday            ! true if |lat|>50 and DayNightFlag is not Night
  69
  70    nargs = iargc()
  71
  72    if (nargs <= 0) then
  73      write(6,*) ' '
  74      write(6,*) 'Usage: decode_airs.exe filename ...'
  75      write(6,*) '    where filename is the name of an HDFEOS file containing'
  76      write(6,*) '    a swath of L2 AIRS retreivals.'
  77      write(6,*) ' '
  78      stop
  79    end if
  80
  81
  82    min_time = '00000000000000'
  83    max_time = '99999999999999'
  84
  85    open(11,file='time_window.nl',status='old',iostat=istatus)
  86    if (istatus == 0) then
  87       read(11,time_window)
  88       close(11)
  89    end if
  90
  91    qual_threshold = 1   ! 0: best, 1:good, 2:not use
  92
  93    open(12,file='qual_threshold.nl',status='old',iostat=istatus)
  94    if (istatus == 0) then
  95       read(12,quality)
  96       close(12)
  97    end if
  98
  99    ! The time for each FOV is given as the number of seconds
 100    !   since 1993 Jan 1 00Z
 101    ref_time = '1993010100 '
 102
 103    open(10, file='soundings.little_r', form='formatted', status='replace')
 104    open(20, file='soundings_polarday.little_r', form='formatted', status='replace')
 105
 106    do fn = 1, nargs
 107       call getarg (fn, file_name)
 108
 109       ! Read all retrievals from file_name and place them in ret
 110       call airs_ret_rdr(file_name, ret, version)
 111
 112       if (granule_out_of_timewin(ret%Time, min_time, max_time)) cycle
 113
 114       nvalid_temp = 0
 115       nvalid_h2o  = 0
 116       npolarday_temp = 0
 117       npolarday_h2o  = 0
 118
 119       write(6,'(a,a)') 'DayNightFlag for this granule: ', trim(ret%DayNightFlag)
 120
 121       do track = 1, AIRS_RET_GEOTRACK
 122       do xtrack = 1, AIRS_RET_GEOXTRACK
 123
 124          polarday = .false.
 125          nvalid = 0
 126
 127          delta_time = nint(ret%Time(xtrack,track))/3600
 128          date_char(1:6) = '      '
 129          call geth_newdate(ref_time, delta_time, date_char(7:20))
 130          mm = mod(nint(ret%Time(xtrack,track)),3600)/60
 131          ss = mod(nint(ret%Time(xtrack,track)),60)
 132          write(date_char(17:20),'(i2.2,i2.2)') mm, ss
 133
 134          if (ret%Time(xtrack,track) < 0.) then
 135             cmp_min = LESS
 136          else
 137             call cmp_datestr(date_char(7:20), min_time, cmp_min)
 138             call cmp_datestr(date_char(7:20), max_time, cmp_max)
 139          end if
 140
 141
 142          if (cmp_min /= LESS .and. cmp_max /= MORE) then
 143
 144             write (6,'(a,f8.4,a5,f9.4,a3)') 'Location: ',                          &
 145                                             ret%Latitude(xtrack,track),'lat, ',    &
 146                                             ret%Longitude(xtrack,track),'lon'
 147             write(6,'(a)') 'Time: '//date_char(7:20)
 148             ! write(6,'(a,i9)') 'RetQAFlag: ', ret%RetQAFlag(xtrack,track)
 149             write(6,*) ' '
 150
 151             !
 152             ! First loop through all levels and find out how many valid obs there are
 153             !
 154             if ( version == 'V4' ) then
 155                do layer=AIRS_RET_STDPRESSURELAY,1,-1
 156
 157                   ! Which Qual flag to check depends on layer
 158                   if (ret%PressStd(layer) < ret%Press_mid_top_bndry(xtrack,track) ) then
 159                      qual = ret%Qual_Temp_Profile_Top(xtrack,track)
 160
 161                   else if (ret%PressStd(layer) < ret%Press_bot_mid_bndry(xtrack,track) ) then
 162                      qual = ret%Qual_Temp_Profile_Mid(xtrack,track)
 163
 164                   else
 165                      qual = ret%Qual_Temp_Profile_Bot(xtrack,track)
 166
 167                   endif
 168
 169                   ! Temperature. If quality is bad (2) then put out flag value of -9999.0
 170                   if (qual <= qual_threshold .and. layer > ret%nSurfStd(xtrack,track)) then
 171                      temp = ret%TAirStd(layer,xtrack,track)
 172                      h2o = ret%H2OMMRStd(layer,xtrack,track)
 173
 174                   ! See p.3 of Level 2 Product Levels and Layers
 175                   ! else if (qual <= qual_threshold .and. abs(ret%PressStd(layer) - ret%PSurfStd(xtrack,track)) < 5.) then
 176 !                 ! temp = ret%TAirStd(layer,xtrack,track)
 177 !                 ! h2o = ret%H2OMMRStd(layer,xtrack,track)
 178
 179                   else
 180                      temp = -9999.0
 181                      h2o = -9999.0
 182
 183                   end if
 184
 185                   if (temp /= -9999.) then
 186                      nvalid = nvalid + 1
 187                      if (h2o /= -9999.) nvalid = nvalid + 1
 188                   end if
 189
 190                end do
 191             end if  ! end of V4 quality selection
 192
 193             if ( version == 'V5' ) then
 194                if ( qual_threshold == 0 ) then       ! Best
 195                   nvalid = AIRS_RET_STDPRESSURELAY - ret%nBestStd(xtrack,track) + 1
 196                else if ( qual_threshold == 1 ) then  ! Good and Best
 197                   if ( ret%nGoodStd(xtrack,track) == 29 ) then
 198                      nvalid = AIRS_RET_STDPRESSURELAY - ret%nBestStd(xtrack,track) + 1
 199                   else
 200                      nvalid = AIRS_RET_STDPRESSURELAY - ret%nGoodStd(xtrack,track) + 1
 201                   end if
 202                end if
 203                if ( (abs(ret%Latitude(xtrack,track)) > 50.0) .and. (trim(ret%DayNightFlag) /= 'Night') ) then
 204                   polarday = .true.
 205                end if
 206             end if  ! end of V5 quality selection
 207             !
 208             ! If there are valid obs to be reported, write them out to little_r format
 209             !
 210             if (nvalid > 0) then
 211
 212                latitude = ret%Latitude(xtrack,track)
 213                longitude = ret%Longitude(xtrack,track)
 214                if ( version == 'V5' ) then
 215                   id       = 'AIRS L2 Std Retrieval '//version//'                '
 216                else
 217                   id       = 'AIRS L2 Std Retrieval                   '
 218                end if
 219                name     = 'AIRSRET                                 '
 220                platform = 'FM-133                                  '
 221                source   = 'GES DAAC                                '
 222                elevation = -888888.
 223                n_valid = nvalid
 224                n_err = 0
 225                n_warning = 0
 226                iseq = 0
 227                n_dups = 0
 228                is_sound = .true.
 229                is_bogus = .false.
 230                is_discard = .false.
 231                isut = -888888
 232                julian = -888888
 233                slp = -888888.
 234                islp = -88
 235                refp = -888888.
 236                irefp = -88
 237                grndt = -888888.
 238                igrndt = -88
 239                sst = -888888.
 240                isst = -88
 241                psfc = -888888.
 242                ipsfc = -88
 243                x1 = -888888.
 244                ix1 = -88
 245                x2 = -888888.
 246                ix2 = -88
 247                x3 = -888888.
 248                ix3 = -88
 249                x4 = -888888.
 250                ix4 = -88
 251                x5 = -888888.
 252                ix5 = -88
 253                x6 = -888888.
 254                ix6 = -88
 255                x7 = -888888.
 256                ix7 = -88
 257                x8 = -888888.
 258                ix8 = -88
 259
 260                if ( .not. polarday ) then
 261                   write(10, fmt=rpt_format) latitude, longitude, id, name, platform, &
 262                      source, elevation, n_valid, n_err, n_warning, iseq, n_dups,     &
 263                      is_sound, is_bogus, is_discard, isut, julian, date_char,        &
 264                      slp, islp, refp, irefp, grndt, igrndt, sst, isst, psfc, ipsfc,  &
 265                      x1, ix1, x2, ix2, x3, ix3, x4, ix4, x5, ix5, x6, ix6, x7, ix7,  &
 266                      x8, ix8
 267                else
 268                   write(20, fmt=rpt_format) latitude, longitude, id, name, platform, &
 269                      source, elevation, n_valid, n_err, n_warning, iseq, n_dups,     &
 270                      is_sound, is_bogus, is_discard, isut, julian, date_char,        &
 271                      slp, islp, refp, irefp, grndt, igrndt, sst, isst, psfc, ipsfc,  &
 272                      x1, ix1, x2, ix2, x3, ix3, x4, ix4, x5, ix5, x6, ix6, x7, ix7,  &
 273                      x8, ix8
 274                end if
 275
 276                do layer=AIRS_RET_STDPRESSURELAY,1,-1
 277
 278                   if ( version == 'V4' ) then
 279                      ! Which Qual flag to check depends on layer
 280                      if (ret%PressStd(layer) < ret%Press_mid_top_bndry(xtrack,track) ) then
 281                         qual = ret%Qual_Temp_Profile_Top(xtrack,track)
 282
 283                      else if (ret%PressStd(layer) < ret%Press_bot_mid_bndry(xtrack,track) ) then
 284                         qual = ret%Qual_Temp_Profile_Mid(xtrack,track)
 285
 286                      else
 287                         qual = ret%Qual_Temp_Profile_Bot(xtrack,track)
 288
 289                      endif
 290
 291                      ! Temperature. If quality is bad (2) then put out flag value of -888888.0
 292                      if (qual <= qual_threshold .and. layer > ret%nSurfStd(xtrack,track)) then
 293                         temp = ret%TAirStd(layer,xtrack,track)
 294                         if (temp == -9999.) temp = -888888.
 295                         h2o = ret%H2OMMRStd(layer,xtrack,track)
 296                         if (h2o == -9999.) h2o = -888888.
 297
 298                      else
 299                         temp = -888888.0
 300                         h2o = -888888.0
 301
 302                      end if
 303                   end if    ! end if V4 check
 304
 305                   if ( version == 'V5' ) then
 306
 307                      ! for temperature
 308
 309                      temp = -888888.0
 310                      if ( qual_threshold == 0 ) then
 311                         if ( layer >= ret%nBestStd(xtrack,track) ) then
 312                            temp = ret%TAirStd(layer,xtrack,track)
 313                         end if
 314                      else if ( qual_threshold == 1 ) then
 315                         if ( layer >= ret%nGoodStd(xtrack,track) .or. layer >= ret%nBestStd(xtrack,track) ) then
 316                            temp = ret%TAirStd(layer,xtrack,track)
 317                         end if
 318                      else
 319                         temp = -888888.0
 320                      end if
 321                      if (temp == -9999.) temp = -888888.
 322
 323                      ! for moisture
 324
 325                      h2o = -888888.0
 326                      if ( qual_threshold == 0 ) then
 327                         if ( (layer <= AIRS_RET_H2OPRESSURELAY) .and.                &
 328                              (ret%Qual_H2O(xtrack,track) <= 1 ) .and.    &
 329                              (ret%pressH2O(layer) <= ret%PBest(xtrack,track)) ) then
 330                            h2o = ret%H2OMMRStd(layer,xtrack,track)
 331                         end if
 332                      else if ( qual_threshold == 1 ) then
 333                         if ( ret%Qual_H2O(xtrack,track) <= qual_threshold ) then
 334                            h2o = ret%H2OMMRStd(layer,xtrack,track)
 335                         end if
 336                      end if
 337                      if ( (ret%H2OMMRStdErr(layer,xtrack,track) < 0.0) .or.   &
 338                           (ret%H2OMMRStdErr(layer,xtrack,track) > 0.5*h2o)  ) then
 339                         h2o = -9999.0
 340                      end if
 341                      if (h2o == -9999.) h2o = -888888.
 342
 343                   end if    ! end if V5 check
 344
 345                   if ( temp > 0.0 ) nvalid_temp = nvalid_temp + 1
 346                   if ( h2o  > 0.0 ) nvalid_h2o  = nvalid_h2o  + 1
 347
 348                   if ( polarday ) then
 349                      if ( temp > 0.0 ) then
 350                         nvalid_temp    = nvalid_temp - 1
 351                         npolarday_temp = npolarday_temp + 1
 352                      end if
 353                      if ( h2o  > 0.0 ) then
 354                         nvalid_h2o     = nvalid_h2o  - 1
 355                         npolarday_h2o  = npolarday_h2o  + 1
 356                      end if
 357                   end if
 358
 359                   if (temp /= -888888.) then
 360
 361                      p = ret%pressStd(layer)*100.
 362                      ip = 0
 363                      z = ret%GP_Height(layer,xtrack,track)
 364                      iz = 0
 365                      t = temp
 366                      it = 0
 367                      dewpt = -888888.
 368                      idewpt = 0
 369                      spd = -888888.
 370                      ispd = 0
 371                      dir = -888888.
 372                      idir = 0
 373                      u = -888888.
 374                      iu = 0
 375                      v = -888888.
 376                      iv = 0
 377                      if (h2o == -888888. .or. t == -888888.) then
 378                         rh = -888888.
 379                      else
 380                         rh = calc_rh(t, h2o/1000., p)
 381                         if (rh < 0.) rh = 0.
 382                      end if
 383                      irh = 0
 384                      thk = -888888.
 385                      ithk = 0
 386                      if ( .not. polarday ) then
 387                         write(10, fmt=meas_format) p, ip, z, iz, t, it, dewpt, idewpt, &
 388                               spd, ispd, dir, idir, u, iu, v, iv, rh, irh, thk, ithk
 389                      else
 390                         write(20, fmt=meas_format) p, ip, z, iz, t, it, dewpt, idewpt, &
 391                               spd, ispd, dir, idir, u, iu, v, iv, rh, irh, thk, ithk
 392                      end if
 393                   end if
 394                end do
 395
 396                p = -777777.
 397                ip = 0
 398                z = -777777.
 399                iz = 0
 400                t = -888888.
 401                it = 0
 402                dewpt = -888888.
 403                idewpt = 0
 404                spd = -888888.
 405                ispd = 0
 406                dir = -888888.
 407                idir = 0
 408                u = -888888.
 409                iu = 0
 410                v = -888888.
 411                iv = 0
 412                rh = -888888.
 413                irh = 0
 414                thk = -888888.
 415                ithk = 0
 416
 417                if ( .not. polarday ) then
 418                   write(10, fmt=meas_format) p, ip, z, iz, t, it, dewpt, idewpt, &
 419                         spd, ispd, dir, idir, u, iu, v, iv, rh, irh, thk, ithk
 420                else
 421                   write(20, fmt=meas_format) p, ip, z, iz, t, it, dewpt, idewpt, &
 422                         spd, ispd, dir, idir, u, iu, v, iv, rh, irh, thk, ithk
 423                end if
 424
 425                num_valid = nvalid
 426                num_err = 0
 427                num_warn = 0
 428                if ( .not. polarday ) then
 429                   write(10, fmt=end_format) num_valid, num_err, num_warn
 430                else
 431                   write(20, fmt=end_format) num_valid, num_err, num_warn
 432                end if
 433
 434             end if  ! end if nvalid > 0
 435
 436             ! Need to output quality info so users can tell if we have a good
 437             ! retrieval with no clouds or a bad retrieval.  Both cases will put
 438             ! all -9999s below.
 439 !            write(6,*) '# Cloud quality:'
 440 !            if (ret%Qual_Cloud_OLR(xtrack,track) .eq. 0) then
 441 !            write(6,*) '    0    Very Good'
 442 !            else if (ret%Qual_Cloud_OLR(xtrack,track) .eq. 1) then
 443 !               write(6,*) '    1    Good'
 444 !            else
 445 !               write(6,*) '    2    Do Not Use'
 446 !            end if
 447
 448          end if ! Profile falls within time window
 449
 450       end do
 451       end do
 452
 453       write(6,*) 'nvalid_temp, nvalid_h2o ', nvalid_temp, nvalid_h2o
 454       write(6,*) 'npolarday_temp, npolarday_h2o ', npolarday_temp, npolarday_h2o
 455
 456    end do ! Loop over filenames
 457
 458    close(10)
 459    close(20)
 460
 461    stop
 462
 463 end program decode_airs
 464
 465 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 466 ! Given two date strings of the form YYYYMMDDHHmmss, icmp is set to
 467 !    0 if the two dates are the same,
 468 !   -1 if date1 < date2, and
 469 !   +1 if date1 > date2.
 470 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 471 subroutine cmp_datestr(date1, date2, icmp)
 472
 473    implicit none
 474
 475    ! Arguments
 476    character (len=14), intent(in) :: date1, date2
 477    integer, intent(out) :: icmp
 478
 479    ! Local variables
 480    integer :: yyyy1, yyyy2
 481    integer :: mm1, mm2
 482    integer :: dd1, dd2
 483    integer :: hh1, hh2
 484    integer :: min1, min2
 485    integer :: sec1, sec2
 486
 487    read(date1(1:4), '(i4)') yyyy1
 488    read(date1(5:6), '(i2)') mm1
 489    read(date1(7:8), '(i2)') dd1
 490    read(date1(9:10), '(i2)') hh1
 491    read(date1(11:12), '(i2)') min1
 492    read(date1(13:14), '(i2)') sec1
 493
 494    read(date2(1:4), '(i4)') yyyy2
 495    read(date2(5:6), '(i2)') mm2
 496    read(date2(7:8), '(i2)') dd2
 497    read(date2(9:10), '(i2)') hh2
 498    read(date2(11:12), '(i2)') min2
 499    read(date2(13:14), '(i2)') sec2
 500
 501    if (yyyy1 > yyyy2) then
 502       icmp = 1
 503       return
 504    else if (yyyy1 < yyyy2) then
 505       icmp = -1
 506       return
 507    end if
 508
 509    if (mm1 > mm2) then
 510       icmp = 1
 511       return
 512    else if (mm1 < mm2) then
 513       icmp = -1
 514       return
 515    end if
 516
 517    if (dd1 > dd2) then
 518       icmp = 1
 519       return
 520    else if (dd1 < dd2) then
 521       icmp = -1
 522       return
 523    end if
 524
 525    if (hh1 > hh2) then
 526       icmp = 1
 527       return
 528    else if (hh1 < hh2) then
 529       icmp = -1
 530       return
 531    end if
 532
 533    if (min1 > min2) then
 534       icmp = 1
 535       return
 536    else if (min1 < min2) then
 537       icmp = -1
 538       return
 539    end if
 540
 541    if (sec1 > sec2) then
 542       icmp = 1
 543       return
 544    else if (sec1 < sec2) then
 545       icmp = -1
 546       return
 547    end if
 548
 549    icmp = 0
 550    return
 551
 552 end subroutine cmp_datestr
 553
 554
 555 function granule_out_of_timewin(time_arr, min_time, max_time)
 556
 557    use read_airs
 558
 559    implicit none
 560
 561    ! Parameters
 562    integer, parameter :: LESS = -1
 563    integer, parameter :: EQUAL = 0
 564    integer, parameter :: MORE = 1
 565
 566    ! Arguments
 567    double precision, dimension(AIRS_RET_GEOXTRACK,AIRS_RET_GEOTRACK), intent(in) :: time_arr
 568    character (len=14), intent(in) :: min_time, max_time
 569
 570    ! Local variables
 571    integer :: i, j
 572    integer :: delta_time, cmp_min, cmp_max, mm, ss
 573    character (len=11) :: ref_time
 574    character (len=14) :: date_char
 575
 576    ! Return value
 577    logical :: granule_out_of_timewin
 578
 579    granule_out_of_timewin = .true.
 580    ref_time = '1993010100 '
 581
 582    do i=1,AIRS_RET_GEOXTRACK,AIRS_RET_GEOXTRACK-1
 583    do j=1,AIRS_RET_GEOTRACK,AIRS_RET_GEOTRACK-1
 584
 585       delta_time = nint(time_arr(i,j))/3600
 586       call geth_newdate(ref_time, delta_time, date_char)
 587       mm = mod(nint(time_arr(i,j)),3600)/60
 588       ss = mod(nint(time_arr(i,j)),60)
 589       write(date_char(11:14),'(i2.2,i2.2)') mm, ss
 590
 591       if (time_arr(i,j) < 0.) then
 592          cmp_min = LESS
 593       else
 594          call cmp_datestr(date_char, min_time, cmp_min)
 595          call cmp_datestr(date_char, max_time, cmp_max)
 596       end if
 597
 598       if (cmp_min /= LESS .and. cmp_max /= MORE) then
 599          granule_out_of_timewin = .false.
 600          return
 601       end if
 602
 603    end do
 604    end do
 605
 606 end function granule_out_of_timewin