var/da/da_obs_io/da_scan_obs_radar.inc

   1 subroutine da_scan_obs_radar (iv, filename, grid)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Scan the radar observation file
   5    !---------------------------------------------------------------------------
   6
   7    implicit none
   8
   9    type (iv_type),    intent(inout) :: iv
  10    character(len=*),  intent(in)    :: filename
  11    type(domain),     intent(in)     :: grid     ! first guess state.
  12
  13    integer                       :: i, j, n, iost, nlevels, fm
  14    integer                       :: total_radar
  15    integer                       :: iunit
  16
  17    type (radar_multi_level_type) :: platform
  18
  19    character (LEN = 120)         :: char_total_radar
  20    character (LEN = 120)         :: char_ned
  21
  22    logical                       :: outside, outside_all
  23    integer                       :: n_dup, ndup
  24
  25    real*8                        :: obs_time
  26    integer                       :: iyear, imonth, iday, ihour, imin
  27
  28
  29    if (trace_use) call da_trace_entry("da_scan_obs_radar")
  30
  31    ! 1. open file
  32    ! ============
  33
  34    call da_get_unit(iunit)
  35    open(unit   = iunit,     &
  36         FILE   = trim(filename), &
  37         FORM   = 'FORMATTED',  &
  38         ACCESS = 'SEQUENTIAL', &
  39         iostat =  iost,     &
  40         STATUS = 'OLD')
  41
  42    if (iost /= 0) then
  43       ! Does not matter of radar file missing
  44       call da_warning(__FILE__,__LINE__, &
  45          (/"Cannot open radar file "//filename/))
  46       call da_free_unit(iunit)
  47       if (trace_use) call da_trace_exit("da_scan_obs_radar")
  48       return
  49    end if
  50    ! 1.1  Initialize
  51    ! ============
  52
  53
  54    ! 2. read total radar
  55    ! ===================
  56
  57    ! 2.1 read first line
  58    !     ---------------
  59
  60    read (unit=iunit, fmt = '(A)', iostat = iost) char_total_radar
  61    if (iost /= 0) then
  62       ! Does matter if present and unreadable
  63       call da_error(__FILE__,__LINE__, &
  64          (/"Cannot read radar file"/))
  65    end if
  66
  67    ! 2.3 total radar number
  68
  69    read (unit=char_total_radar (15:17),fmt='(I3)', iostat = iost) total_radar
  70
  71    ! 2.4 skip one lines
  72
  73    read (unit=iunit, fmt = '(A)', iostat = iost)
  74
  75    ! 3. read radar data
  76
  77    do n = 1, total_radar
  78
  79       ! 3.1 skip one blank line
  80
  81       read (unit=iunit, fmt = '(A)', iostat = iost)
  82
  83       ! 3.2 read header
  84
  85       read (unit=iunit, fmt = '(A)', iostat = iost) char_ned
  86       if (iost /= 0) then
  87          write(unit=message(1),fmt='(A,I3)')'Error reading header information for station number ',n
  88          call da_error(__FILE__,__LINE__,message(1:1))
  89       end if
  90
  91       ! 3.3 read header information
  92
  93       read (unit=char_ned (69:74), fmt='(I6)', iostat = iost) platform % stn % numobs
  94
  95       if(char_ned (1:5).ne."RADAR") then
  96          write(unit=message(1),fmt='(A,I3)')'Error reading header information for station number ',n
  97          call da_error(__FILE__,__LINE__,message(1:1))
  98       end if
  99
 100       ! 3.4 skip two lines
 101
 102       read (unit=iunit, fmt = '(A)', iostat = iost)
 103       read (unit=iunit, fmt = '(A)', iostat = iost)
 104
 105       ! 3.5 loop over records
 106
 107       reports: do j = 1, platform % stn % numobs
 108
 109          ! 3.5.1 read station general info
 110
 111          read (unit = iunit, iostat = iost, &
 112                       fmt = '(A12,3X,A19,2X,2(F12.3,2X),F8.1,2X,I6)') &
 113                       platform % info % platform,  &
 114                       platform % info % date_char, &
 115                       platform % info % lat,       &
 116                       platform % info % lon,       &
 117                       platform % info % elv,       &
 118                       platform % info % levels
 119
 120          if (iost /= 0) then
 121             write(unit=message(1),fmt='(A)')'Error reading next observation for '
 122             write(unit=message(2),fmt='(A)')char_ned
 123             write(unit=message(3),fmt='(A,I6,A,I6)')'Observation number ',j,' out of ',platform % stn % numobs
 124             call da_error(__FILE__,__LINE__,message(1:3))
 125          end if
 126
 127          if (platform%info%lon == 180.0  ) platform%info%lon =-180.000
 128          ! Fix funny wind direction at Poles
 129          if (platform%info%lat < -89.9999 .or. platform%info%lat > 89.9999) then
 130             platform%info%lon = 0.0
 131          end if
 132
 133          read(unit=platform % info % platform (4:6), fmt='(I3)') fm
 134
 135          !     3.5.2 read each level
 136
 137          do i = 1, platform % info % levels
 138             ! height
 139             platform%each (i) = radar_each_level_type(missing_r, missing, -1.0,&
 140                field_type(missing_r, missing, missing_r, missing, missing_r), & ! rv
 141                field_type(missing_r, missing, missing_r, missing, missing_r))   ! rf
 142
 143             read (unit = iunit, fmt = '(3X, F12.1, 2(F12.3,I4,F12.3,2X))', iostat = iost) &
 144                              platform % each (i) % height,           &
 145                              platform % each (i) % rv % inv,         &
 146                              platform % each (i) % rv % qc,          &
 147                              platform % each (i) % rv % error,       &
 148                              platform % each (i) % rf % inv,         &
 149                              platform % each (i) % rf % qc,          &
 150                              platform % each (i) % rf % error
 151
 152          if (iost /= 0) then
 153             write(unit=message(1),fmt='(A)')'Error reading next level for '
 154             write(unit=message(2),fmt='(A)')char_ned
 155             write(unit=message(3),fmt='(A,I6,A,I6)')'Observation number ',j,' out of ',platform % stn % numobs
 156             call da_error(__FILE__,__LINE__,message(1:3))
 157          end if
 158
 159             if (platform % each (i) % rv % error == 0.0) then
 160                  platform % each (i) % rv % error  = 1.0
 161             end if
 162
 163             if (platform % each (i) % rf % error == 0.0) then
 164                  platform % each (i) % rf % error  = 1.0
 165             end if
 166
 167             if (platform % each (i) % rv % inv   == missing_r .or. &
 168                 platform % each (i) % rv % error == missing_r) then
 169                 platform % each (i) % rv % qc     = missing_data
 170             end if
 171
 172             if (platform % each (i) % rf % inv   == missing_r .or. &
 173                 platform % each (i) % rf % error == missing_r) then
 174                 platform % each (i) % rf % qc     = missing_data
 175             end if
 176
 177          end do
 178
 179          ! Check if outside of the time range:
 180
 181          read (platform%info%date_char,'(i4,4(1x,i2))', iostat = iost) &
 182                iyear, imonth, iday, ihour, imin
 183          if (iost /= 0) then
 184             write(unit=message(1),fmt='(A,A,A)')'Bad date: ',platform % info % date_char,' found for'
 185             write(unit=message(2),fmt='(A)')char_ned
 186             write(unit=message(3),fmt='(A,I6,A,I6)')'Observation number ',j,' out of ',platform % stn % numobs
 187             write(unit=message(4),fmt='(A)')'This is usually due to an incorrect number of vertical levels in the previous observation'
 188             call da_error(__FILE__,__LINE__,message(1:4))
 189          end if
 190          call da_get_julian_time (iyear,imonth,iday,ihour,imin,obs_time)
 191          if ( obs_time < time_slots(0) .or. &
 192               obs_time >= time_slots(num_fgat_time) ) then
 193             cycle reports
 194          endif
 195
 196          call da_llxy (platform%info, platform%loc, outside, outside_all)
 197          if( outside_all .and. multi_inc == 0 ) cycle reports
 198
 199          nlevels = platform%info%levels
 200
 201          if (nlevels > max_ob_levels) then
 202              write(unit=message(1),fmt='(A,2I8)') &
 203                 ' radar=> nlevels > max_ob_levels:',nlevels, max_ob_levels
 204              call da_warning(__FILE__,__LINE__,message(1:1))
 205
 206              nlevels = max_ob_levels
 207              platform%info%levels = nlevels
 208          else if (nlevels < 1) then
 209             cycle reports
 210          end if
 211
 212
 213          ! Loop over duplicating obs for global
 214          n_dup = 1
 215          if (global .and. (platform%loc%i == ids .or. platform%loc%i == ide)) n_dup= 2
 216
 217          do ndup = 1, n_dup
 218             select case (fm)
 219
 220             case (128)
 221                if (.not.use_radarobs .or. iv%info(radar)%ntotal == max_radar_input) cycle reports
 222                if (ndup==1) iv%info(radar)%ntotal = iv%info(radar)%ntotal + 1
 223                if (outside) cycle reports
 224                iv%info(radar)%nlocal = iv%info(radar)%nlocal + 1
 225
 226             case default;
 227                write(unit=stdout, fmt='(a)') 'Warning: unsaved obs found:'
 228
 229                write(unit=stdout, fmt='(2a)') &
 230                   'platform % info % platform=', platform % info % platform
 231
 232                write(unit=stdout, fmt='(a, i3)') &
 233                   'platform % info % levels=', platform % info % levels
 234             end select
 235
 236             iv%info(radar)%max_lev = max(iv%info(radar)%max_lev, platform%info%levels)
 237          end do        !  loop over duplicate
 238       end do reports
 239
 240    end do ! total_radar
 241
 242    close (iunit)
 243    call da_free_unit(iunit)
 244
 245    if (trace_use) call da_trace_exit("da_scan_obs_radar")
 246
 247 end subroutine da_scan_obs_radar