var/da/da_ssmi/da_scan_obs_ssmi.inc

   1 subroutine da_scan_obs_ssmi (iv, filename)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Read the header of a SSMI 2.0 GTS observation file
   5    !---------------------------------------------------------------------------
   6
   7    implicit none
   8
   9    type(iv_type),    intent(inout) :: iv
  10    character(len=*), intent(in)    :: filename
  11
  12    character (len =  10)        :: fmt_name
  13    character (len = 160)        :: fmt_info, fmt_loc
  14    character (len = 120)        :: char_ned
  15
  16    integer                      :: iost, fm,iunit
  17
  18    type (model_loc_type)        :: loc
  19    type (info_type)             :: info
  20    type (field_type)            :: speed, tpw
  21
  22    type (field_type)            :: tb19v, tb19h, tb22v
  23    type (field_type)            :: tb37v, tb37h, tb85v, tb85h
  24
  25    logical                      :: isfilter,ipass
  26    logical                      :: outside
  27    integer                      :: irain, iprecip
  28
  29    if (trace_use) call da_trace_entry("da_scan_obs_ssmi")
  30
  31    isfilter = .true. ! filter out rain points
  32    irain = 0
  33
  34    ! open FILE
  35    call da_get_unit(iunit)
  36    open(unit   = iunit,     &
  37       FILE   = trim(filename), &
  38       FORM   = 'FORMATTED',  &
  39       ACCESS = 'SEQUENTIAL', &
  40       iostat =  iost,     &
  41       STATUS = 'OLD')
  42
  43    if (iost /= 0) then
  44       call da_warning(__FILE__,__LINE__, (/"Cannot open SSMI file "//filename/))
  45       call da_free_unit(iunit)
  46       return
  47    end if
  48
  49    rewind (unit = iunit)
  50
  51    ! 2.  read header
  52    ! ===============
  53
  54    ! 2.1 read first line
  55    !     ---------------
  56
  57    read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
  58    if (iost /= 0) then
  59       call da_error(__FILE__,__LINE__, (/"Cannot read SSMI file "//filename/))
  60    end if
  61
  62    ! 2.3 read number of reports
  63    do
  64       read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
  65       if (iost /= 0) then
  66          call da_error(__FILE__,__LINE__, (/"Cannot read SSMI file "//filename/))
  67       end if
  68       if (char_ned(1:6) == 'NESTIX') exit
  69    end do
  70
  71    do
  72       read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
  73       if (char_ned(1:6) == 'INFO  ') exit
  74    end do
  75
  76    read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
  77
  78    ! read formats
  79    read (unit=iunit, fmt = '(A,1X,A)') &
  80       fmt_name, fmt_info, &
  81       fmt_name, fmt_loc
  82
  83    ! skip 1 line
  84    read (unit=iunit, fmt = '(A)') fmt_name
  85
  86    ! loop over records
  87    reports: do
  88       ! read station general info
  89       read (unit=iunit, fmt = fmt_info, iostat = iost) &
  90          info % platform,    &
  91          info % date_char,   &
  92          info % name,        &
  93          info % levels,      &
  94          info % lat,         &
  95          info % lon,         &
  96          info % elv,         &
  97          info % id
  98
  99       read(unit=info % platform (4:6),fmt='(I3)') fm
 100       if (iost /= 0) exit reports
 101
 102       select case(fm)
 103       case (125)    ;
 104          ! read surface wind speed and precipitable water
 105          read (unit=iunit, fmt = fmt_loc) &
 106                speed%inv, speed%qc, speed%error, &
 107                tpw%inv, tpw%qc, tpw%error
 108       case (126)    ;
 109          read (unit=iunit, fmt = fmt_loc)  &
 110                tb19v%inv, tb19v%qc, tb19v%error, &
 111                tb19h%inv, tb19h%qc, tb19h%error, &
 112                tb22v%inv, tb22v%qc, tb22v%error, &
 113                tb37v%inv, tb37v%qc, tb37v%error, &
 114                tb37h%inv, tb37h%qc, tb37h%error, &
 115                tb85v%inv, tb85v%qc, tb85v%error, &
 116                tb85h%inv, tb85h%qc, tb85h%error
 117       case default;
 118          write(unit=message(1), fmt='(a, i6)') 'unsaved ssmi obs found, fm=', fm
 119          write(unit=message(2), fmt='(a, 2f12.6)') 'info%(lon,lat)=', info%lon, info%lat
 120          call da_warning(__FILE__,__LINE__,message(1:2))
 121       end select
 122
 123       ! check if obs is in horizontal domain
 124
 125       ! Compute the model horizontal coordinate x, y
 126       ! Check if obs is wihin horizontal domain
 127
 128       call da_llxy (info, loc, outside)
 129
 130       select case(fm)
 131       case (125) ;
 132          if (.not. use_ssmiretrievalobs .or. iv%info(ssmi_rv)%ntotal == max_ssmi_rv_input) cycle reports
 133
 134          ! Check if at least one field is present
 135          if ((tpw % qc == missing) .AND. (speed % qc == missing)) then
 136             cycle reports
 137          end if
 138          iv%info(ssmi_rv)%ntotal = iv%info(ssmi_rv)%ntotal + 1
 139          if (outside) cycle reports
 140
 141          iv%info(ssmi_rv)%nlocal = iv%info(ssmi_rv)%nlocal + 1
 142
 143       case (126) ;
 144          if (.not. use_ssmitbobs .or. iv%info(ssmi_tb)%ntotal == max_ssmi_tb_input) cycle reports
 145
 146          ! Check if at least one field is present
 147
 148          if ((tb19v % qc == missing) .AND. (tb19h % qc == missing)  .AND. &
 149              (tb22v % qc == missing)                                .AND. &
 150              (tb37v % qc == missing) .AND. (tb37h % qc == missing)  .AND. &
 151              (tb85v % qc == missing) .AND. (tb85h % qc == missing)) then
 152             cycle reports
 153          end if
 154
 155          ! filter rain pixels
 156          ! ====================================
 157
 158          if (isfilter) then
 159             ipass = .false.
 160             iprecip = 0
 161             call filter(tb19v%inv, tb19h%inv, tb22v%inv, tb37v%inv, &
 162                tb37h%inv, tb85v%inv, tb85h%inv, ipass, iprecip, &
 163                info%lat)
 164             if (iprecip == 1) then
 165                irain = irain + 1
 166                cycle reports
 167             end if
 168          end if
 169
 170          iv%info(ssmi_tb)%ntotal = iv%info(ssmi_tb)%ntotal + 1
 171          if (outside) cycle reports
 172          iv%info(ssmi_tb)%nlocal = iv%info(ssmi_tb)%nlocal + 1
 173
 174       case default;
 175          ! Do nothing.
 176       end select
 177
 178    end do reports
 179
 180    iv%info(ssmt1)%max_lev   = 1
 181    iv%info(ssmt2)%max_lev   = 1
 182    iv%info(ssmi_tb)%max_lev = 1
 183    iv%info(ssmi_rv)%max_lev = 1
 184
 185    close(unit=iunit)
 186    call da_free_unit(iunit)
 187
 188    if (irain > 0) then
 189       write(unit=stdout, fmt='(/,5x,a,i6/)') 'Rejected rain contaminated ssmi_tb =', irain
 190       write(unit=stdout, fmt='(A)') ' '
 191    end if
 192
 193    if (trace_use) call da_trace_exit("da_scan_obs_ssmi")
 194
 195 end subroutine da_scan_obs_ssmi
 196
 197