var/da/da_obs_io/da_read_obs_radar.inc

   1 subroutine da_read_obs_radar (iv, filename, grid)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: Read 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, nn, 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    integer                       :: ntotal,nlocal,ilocal
  29    integer                       :: radar_nlocal
  30    real, allocatable             :: in(:), out(:)
  31
  32
  33    if (trace_use) call da_trace_entry("da_read_obs_radar")
  34
  35    ntotal = iv%info(radar)%ptotal(iv%time-1)
  36    nlocal = iv%info(radar)%plocal(iv%time-1)
  37    ilocal = nlocal
  38    radar_nlocal = nlocal
  39
  40    ! 1. open file
  41
  42    call da_get_unit(iunit)
  43    open(unit   = iunit,     &
  44         FILE   = trim(filename), &
  45         FORM   = 'FORMATTED',  &
  46         ACCESS = 'SEQUENTIAL', &
  47         iostat =  iost,     &
  48         STATUS = 'OLD')
  49
  50    if (iost /= 0) then
  51       ! Missing file does not matter
  52       call da_warning(__FILE__,__LINE__, &
  53          (/"Cannot open radar file "//filename/))
  54       call da_free_unit(iunit)
  55       if (trace_use) call da_trace_exit("da_read_obs_radar")
  56       return
  57    end if
  58
  59    ! 2. read total radar
  60
  61    !  2.1 read first line
  62
  63    read (unit=iunit, fmt = '(A)', iostat = iost) char_total_radar
  64
  65    !  2.2 process error
  66
  67    if (iost /= 0) then
  68      call da_error(__FILE__,__LINE__, &
  69         (/"Cannot read radar file"/))
  70    end if
  71
  72    !  2.3 total radar number
  73
  74    read (unit=char_total_radar (15:17),fmt='(I3)', iostat = iost) total_radar
  75
  76    if (print_detail_radar) write (unit=stdout,fmt='(/,A,I3,/)') &
  77        ' TOTAL RADAR: ', total_radar
  78
  79    !  2.4 skip one line
  80
  81    read (unit=iunit, fmt = '(A)', iostat = iost)
  82
  83    ! 3. read radar data
  84
  85    do nn = 1, total_radar
  86
  87       ! 3.1 skip one blank line
  88
  89       read (unit=iunit, fmt = '(A)', iostat = iost)
  90
  91       ! 3.2 read header
  92
  93       read (unit=iunit, fmt = '(A)', iostat = iost) char_ned
  94
  95       ! 3.3 read header information
  96
  97       read (unit=char_ned (1:5),fmt='(A5)', iostat = iost) platform % stn % platform
  98
  99       if (print_detail_radar) write (unit=stdout,fmt='(A)') 'RADAR Observations information'
 100
 101       read (unit=char_ned (8:19),fmt='(A12)', iostat = iost) platform % stn % name
 102
 103 !      if (print_detail_radar) write (unit=stdout,fmt='(A,A5,A,A12)')  &
 104        write (unit=stdout,fmt='(A,A5,A,A12)')  &
 105                            ' Reading ',platform % stn % platform, &
 106                            ' data at station:', platform % stn % name
 107
 108       read (unit=char_ned(20:27),fmt='(F8.3)', iostat = iost) platform % stn % lon
 109
 110       read (unit=char_ned (30:37),fmt='(F8.3)', iostat = iost) platform % stn % lat
 111
 112       read (unit=char_ned (40:47),fmt='(F8.1)', iostat = iost) platform % stn % elv
 113
 114       if (print_detail_radar) write (unit=stdout,fmt='(A,2(F8.3,2X),F8.1)')  &
 115          'The station longitude, latitude, and altitude are: ', &
 116          platform % stn % lon, &
 117          platform % stn % lat, platform % stn % elv
 118
 119       read (unit=char_ned (50:68),fmt='(A19)', iostat = iost) platform % stn % date_char
 120
 121       if (print_detail_radar) write (unit=stdout,fmt='(A,A19)')   &
 122          'The observation time for this data is ',     &
 123          platform % stn % date_char
 124
 125       read (unit=char_ned (69:74),fmt='(I6)', iostat = iost) platform % stn % numobs
 126
 127       if (print_detail_radar) write (unit=stdout,fmt='(A,I6)')   &
 128          'Total number of Super-observations is ', &
 129          platform % stn % numobs
 130
 131
 132       read (unit=char_ned (75:80),fmt='(I6)', iostat = iost) platform % stn % levels
 133
 134       if (print_detail_radar) write (unit=stdout,fmt='(A,I6)')   &
 135          'Vertical layers for each Super-observation is ', &
 136          platform % stn % levels
 137
 138       ! 3.4 skip two lines
 139
 140       read (unit=iunit, fmt = '(A)', iostat = iost)
 141       read (unit=iunit, fmt = '(A)', iostat = iost)
 142
 143       ! 3.5 loop over records
 144
 145       reports: do j = 1, platform % stn % numobs
 146
 147          ! 3.5.1 read station general info
 148
 149          read (unit = iunit, iostat = iost, &
 150                       fmt = '(A12,3X,A19,2X,2(F12.3,2X),F8.1,2X,I6)') &
 151                       platform % info % platform,  &
 152                       platform % info % date_char, &
 153                       platform % info % lat,       &
 154                       platform % info % lon,       &
 155                       platform % info % elv,       &
 156                       platform % info % levels
 157       if (platform%info%lon == 180.0  ) platform%info%lon =-180.000
 158       ! Fix funny wind direction at Poles
 159       if (platform%info%lat < -89.9999 .or. platform%info%lat > 89.9999) then
 160          platform%info%lon = 0.0
 161       end if
 162
 163          read(platform % info % platform (4:6), '(I3)') fm
 164
 165          ! 3.5.2 read each level
 166
 167          do i = 1, platform % info % levels
 168             ! height
 169             platform%each(i) = radar_each_level_type(missing_r, missing, -1.0,       &
 170                field_type(missing_r, missing, missing_r, missing, missing_r), & ! rv
 171                field_type(missing_r, missing, missing_r, missing, missing_r))   ! rf
 172
 173             read (unit = iunit, fmt = '(3X, F12.1, 2(F12.3,I4,F12.3,2X))') &
 174                              platform % each (i) % height,           &
 175                              platform % each (i) % rv % inv,         &
 176                              platform % each (i) % rv % qc,          &
 177                              platform % each (i) % rv % error,       &
 178                              platform % each (i) % rf % inv,         &
 179                              platform % each (i) % rf % qc,          &
 180                              platform % each (i) % rf % error
 181
 182             if (platform % each (i) % rv % error == 0.0 .or. &
 183                 abs(platform % each (i) % rv % error - missing_r) < 1.0) then
 184                  platform % each (i) % rv % error  = 1.0
 185             end if
 186
 187             if (platform % each (i) % rf % error == 0.0 .or. &
 188                 abs(platform % each (i) % rf % error - missing_r) < 1.0) then
 189                  platform % each (i) % rf % error  = 1.0
 190             end if
 191
 192             if (platform % each (i) % rv % inv   == missing_r .or. &
 193                 platform % each (i) % rv % error == missing_r) then
 194                 platform % each (i) % rv % qc     = missing_data
 195             end if
 196
 197             if (platform % each (i) % rf % inv   == missing_r .or. &
 198                 platform % each (i) % rf % error == missing_r) then
 199                 platform % each (i) % rf % qc     = missing_data
 200             end if
 201          end do
 202
 203          ! Check if outside of the time range:
 204
 205          read (platform%info%date_char,'(i4,4(1x,i2))') &
 206                iyear, imonth, iday, ihour, imin
 207          call da_get_julian_time (iyear,imonth,iday,ihour,imin,obs_time)
 208          ! If you skip this part, you have to skip in da_scan_obs_radar.inc too!
 209          if ( obs_time < time_slots(0) .or. &
 210               obs_time >= time_slots(num_fgat_time) ) then
 211             if (print_detail_radar) then
 212                write(unit=stdout, fmt='(a)') '*** Outside of the time range:'
 213                write(unit=stdout, fmt=fmt_info) &
 214                      platform%info%platform,    &
 215                      platform%info%date_char,   &
 216                      platform%stn%name
 217             end if
 218             cycle reports
 219          endif
 220
 221          call da_llxy (platform%info, platform%loc, outside, outside_all)
 222          if( outside_all .and. multi_inc == 0 ) then
 223             if (print_detail_radar) then
 224                write(unit=stdout, fmt='(a)') '*** Report is outside of domain:'
 225                write(unit=stdout, fmt='(2x,a,2(2x,f8.3),2x,a)') &
 226                      platform%info%platform,    &
 227                      platform%info%lat,   &
 228                      platform%info%lon,   &
 229                      platform%stn%name
 230             end if
 231             cycle reports
 232          end if
 233
 234         read (analysis_date,'(i4,4(1x,i2))') &
 235                                     iyear, imonth, iday, ihour, imin
 236
 237          nlevels = platform%info%levels
 238
 239             iv%info(radar)%max_lev = max(iv%info(radar)%max_lev, platform%info%levels)
 240
 241          if (nlevels > max_ob_levels) then
 242             write(unit=message(1),fmt='(A,2I8)') &
 243                ' radar=> nlevels > max_ob_levels:',nlevels, max_ob_levels
 244             call da_warning(__FILE__,__LINE__,message(1:1))
 245             nlevels = max_ob_levels
 246              platform%info%levels = nlevels
 247          else if (nlevels < 1) then
 248             cycle reports
 249          end if
 250
 251          ! Loop over duplicating obs for global
 252          n_dup = 1
 253          if (global .and. &
 254             (platform%loc%i == ids .or. platform%loc%i == ide)) n_dup= 2
 255          do ndup = 1, n_dup
 256             select case (fm)
 257             case (128)
 258             if (.not.use_radarobs .or. ntotal == max_radar_input) cycle reports
 259
 260                if (ndup==1 ) ntotal = ntotal + 1
 261                if (outside) cycle reports
 262                nlocal = nlocal + 1
 263                ilocal = ilocal + 1
 264                iv % radar (ilocal) % stn_loc % lon = platform % stn % lon
 265                iv % radar (ilocal) % stn_loc % lat = platform % stn % lat
 266                iv % radar (ilocal) % stn_loc % elv = platform % stn % elv
 267
 268                iv%info(radar)%levels(ilocal)    = nlevels
 269                iv%info(radar)%name(ilocal)      = platform%info%name
 270                iv%info(radar)%platform(ilocal)  = platform%info%platform
 271                iv%info(radar)%id(ilocal)        = platform%info%id
 272                iv%info(radar)%date_char(ilocal) = platform%info%date_char
 273                iv%info(radar)%lat(:,ilocal)     = platform%info%lat
 274                iv%info(radar)%lon(:,ilocal)     = platform%info%lon
 275                iv%info(radar)%elv(ilocal)       = platform%info%elv
 276                iv%info(radar)%x(:,ilocal)           = platform%loc%x
 277                iv%info(radar)%y(:,ilocal)           = platform%loc%y
 278                iv%info(radar)%i(:,ilocal)           = platform%loc%i
 279                iv%info(radar)%j(:,ilocal)           = platform%loc%j
 280                iv%info(radar)%dx(:,ilocal)          = platform%loc%dx
 281                iv%info(radar)%dxm(:,ilocal)         = platform%loc%dxm
 282                iv%info(radar)%dy(:,ilocal)          = platform%loc%dy
 283                iv%info(radar)%dym(:,ilocal)         = platform%loc%dym
 284                iv%info(radar)%proc_domain(:,ilocal) = platform%loc%proc_domain
 285                iv%info(radar)%obs_global_index(ilocal) = ntotal
 286                if( nlocal == ilocal) then
 287                   allocate (iv % radar (ilocal) % model_p  (1:iv%info(radar)%max_lev))
 288                   allocate (iv % radar (ilocal) % model_rho(1:iv%info(radar)%max_lev))
 289                   allocate (iv % radar (ilocal) % model_qrn(1:iv%info(radar)%max_lev))
 290                   allocate (iv % radar (ilocal) % model_qcl(1:iv%info(radar)%max_lev))
 291                   allocate (iv % radar (ilocal) % model_qci(1:iv%info(radar)%max_lev))
 292                   allocate (iv % radar (ilocal) % model_qsn(1:iv%info(radar)%max_lev))
 293                   allocate (iv % radar (ilocal) % model_qgr(1:iv%info(radar)%max_lev))
 294                   allocate (iv % radar (ilocal) % model_zmm(1:iv%info(radar)%max_lev))
 295                   allocate (iv % radar (ilocal) % height   (1:iv%info(radar)%max_lev))
 296                   allocate (iv % radar (ilocal) % height_qc(1:iv%info(radar)%max_lev))
 297                   allocate (iv % radar (ilocal) % rv       (1:iv%info(radar)%max_lev))
 298                   allocate (iv % radar (ilocal) % rf       (1:iv%info(radar)%max_lev))
 299                   allocate (iv % radar (ilocal) % zmm      (1:iv%info(radar)%max_lev))
 300                   if ( use_radar_rhv ) then
 301                      allocate (iv % radar (ilocal) % rrn      (1:iv%info(radar)%max_lev))
 302                      allocate (iv % radar (ilocal) % rcl      (1:iv%info(radar)%max_lev))
 303                      allocate (iv % radar (ilocal) % rci      (1:iv%info(radar)%max_lev))
 304                      allocate (iv % radar (ilocal) % rsn      (1:iv%info(radar)%max_lev))
 305                      allocate (iv % radar (ilocal) % rgr      (1:iv%info(radar)%max_lev))
 306                      allocate (iv % radar (ilocal) % rrno     (1:iv%info(radar)%max_lev))
 307                      allocate (iv % radar (ilocal) % rclo     (1:iv%info(radar)%max_lev))
 308                      allocate (iv % radar (ilocal) % rcio     (1:iv%info(radar)%max_lev))
 309                      allocate (iv % radar (ilocal) % rsno     (1:iv%info(radar)%max_lev))
 310                      allocate (iv % radar (ilocal) % rgro     (1:iv%info(radar)%max_lev))
 311                   end if
 312                   if ( use_radar_rqv ) then
 313                      allocate (iv % radar (ilocal) % rqv      (1:iv%info(radar)%max_lev))
 314                      allocate (iv % radar (ilocal) % rqvo     (1:iv%info(radar)%max_lev))
 315                   end if
 316                end if
 317                do i = 1, nlevels
 318                   iv % radar (ilocal) % height(i)    = platform % each(i) % height
 319                   iv % radar (ilocal) % height_qc(i) = platform % each(i) % height_qc
 320                   iv % radar (ilocal) % rv(i)        = platform % each(i) % rv
 321                   iv % radar (ilocal) % rf(i)        = platform % each(i) % rf
 322
 323                   if ( use_radar_rhv ) then
 324                      iv % radar (ilocal) % rrn(i) % inv   = missing_r
 325                      iv % radar (ilocal) % rrn(i) % qc    = missing_data
 326                      iv % radar (ilocal) % rrn(i) % error = missing_r
 327                      iv % radar (ilocal) % rrno(i)        = missing_r
 328
 329                      iv % radar (ilocal) % rcl(i) % inv   = missing_r
 330                      iv % radar (ilocal) % rcl(i) % qc    = missing_data
 331                      iv % radar (ilocal) % rcl(i) % error = missing_r
 332                      iv % radar (ilocal) % rclo(i)        = missing_r
 333
 334                      iv % radar (ilocal) % rci(i) % inv   = missing_r
 335                      iv % radar (ilocal) % rci(i) % qc    = missing_data
 336                      iv % radar (ilocal) % rci(i) % error = missing_r
 337                      iv % radar (ilocal) % rcio(i)        = missing_r
 338
 339                      iv % radar (ilocal) % rsn(i) % inv   = missing_r
 340                      iv % radar (ilocal) % rsn(i) % qc    = missing_data
 341                      iv % radar (ilocal) % rsn(i) % error = missing_r
 342                      iv % radar (ilocal) % rsno(i)        = missing_r
 343
 344                      iv % radar (ilocal) % rgr(i) % inv   = missing_r
 345                      iv % radar (ilocal) % rgr(i) % qc    = missing_data
 346                      iv % radar (ilocal) % rgr(i) % error = missing_r
 347                      iv % radar (ilocal) % rgro(i)        = missing_r
 348                   end if
 349
 350                   if ( use_radar_rqv ) then
 351                      iv % radar (ilocal) % rqv(i) % inv   = missing_r
 352                      iv % radar (ilocal) % rqv(i) % qc    = missing_data
 353                      iv % radar (ilocal) % rqv(i) % error = missing_r
 354                      iv % radar (ilocal) % rqvo(i)        = missing_r
 355                   end if
 356                end do
 357
 358             case default;
 359                write(unit=message(1), fmt='(a)') 'Unsaved obs found:'
 360                write(unit=message(2), fmt='(2a)') &
 361                   'platform % info % platform=', platform % info % platform
 362                write(unit=message(3), fmt='(a, i3)') &
 363                   'platform % info % levels=', platform % info % levels
 364                call da_warning(__FILE__,__LINE__,message(1:3))
 365             end select
 366             if (global .and. ndup == 1) then
 367                if (platform%loc % i >= ide) then
 368                   platform%loc%i = ids
 369                   platform%loc%proc_domain = .false.
 370                else if (platform%loc % i <= ids) then
 371                   platform%loc%i = ide
 372                   platform%loc%proc_domain = .false.
 373                end if
 374             end if
 375          end do        !  loop over duplicate
 376       end do reports
 377
 378                 radar_nlocal = nlocal
 379
 380
 381    end do  ! total_radar
 382
 383    if (print_detail_radar) write (unit=stdout,fmt='(/,A,I3,/)') &
 384        ' Processed TOTAL RADAR: ', total_radar
 385    close(iunit)
 386    call da_free_unit(iunit)
 387    if (trace_use) call da_trace_exit("da_read_obs_radar")
 388 end subroutine da_read_obs_radar
 389