var/da/da_obs_io/da_scan_obs_rain.inc

   1 subroutine da_scan_obs_rain (iv, filename, ifgat)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Scan the rain observation file
   5    !---------------------------------------------------------------------------
   6
   7    implicit none
   8
   9    type (iv_type),    intent(inout) :: iv
  10    character(len=*),  intent(in)    :: filename
  11    integer,           intent(in)    :: ifgat
  12    integer                       :: i, j, n, iost, nlevels, fm
  13    integer                       :: file_rain
  14    integer                       :: iunit
  15
  16    type (rain_single_level_type) :: platform
  17
  18    character (len = 120)         :: char_file_rain
  19    character (len = 120)         :: fmt_name
  20    character (len = 160)         :: info_string
  21    logical                       :: outside, outside_all
  22    integer                       :: n_dup, ndup
  23
  24    real*8                        :: obs_time
  25    integer                       :: iyear, imonth, iday, ihour, imin
  26
  27    ! for thinning
  28    real                          :: dlat_earth,dlon_earth,crit
  29    integer                       :: itt,itx,iout
  30    logical                       :: iuse
  31    integer                       :: tp, nlocal
  32    real, allocatable             :: in(:), out(:)
  33
  34    integer                      :: num_outside_all, num_outside_time, num_thinned, num_report
  35
  36    if (trace_use) call da_trace_entry("da_scan_obs_rain")
  37
  38    num_report       = 0
  39    num_outside_all  = 0
  40    num_outside_time = 0
  41    num_thinned      = 0
  42    tp               = iv%info(rain)%plocal(ifgat-1)
  43    nlocal           = 0
  44
  45    ! 1. open file
  46    ! ============
  47    call da_get_unit(iunit)
  48    open(unit   = iunit,     &
  49         FILE   = trim(filename), &
  50         FORM   = 'FORMATTED',  &
  51         ACCESS = 'SEQUENTIAL', &
  52         iostat =  iost,     &
  53         STATUS = 'OLD')
  54
  55    if (iost /= 0) then
  56       write(unit=message(1),fmt='(A,I5,A)') &
  57          "Error",iost," opening rainfall obs file "//trim(filename)
  58       call da_warning(__FILE__,__LINE__,message(1:1))
  59       call da_free_unit(iunit)
  60       if (trace_use) call da_trace_exit("da_scan_obs_rain")
  61       return
  62    end if
  63
  64    ! 2. read rainfall data
  65    ! ===================
  66
  67    ! 2.1 read first line
  68    !     ---------------
  69
  70    read (unit=iunit, fmt = '(A)', iostat = iost) char_file_rain
  71
  72    ! 2.2 process error
  73
  74    if (iost /= 0) then
  75       ! Does matter if present and unreadable
  76       call da_error(__FILE__,__LINE__, &
  77          (/"Cannot read rainfall file"/))
  78    end if
  79
  80    ! 2.3read header info
  81
  82    head_info: do
  83
  84       read (unit=iunit, fmt = '(A)', iostat = iost) info_string
  85
  86       if (iost /= 0) then
  87          write(unit=message(1),fmt='(A,I3,A,I3)') &
  88             "Error",iost,"reading rainfall obs header on unit",iunit
  89          call da_warning(__FILE__,__LINE__,message(1:1))
  90       if (trace_use) call da_trace_exit("da_scan_obs_rain")
  91          return
  92       end if
  93
  94       if (info_string(1:6) == 'EACH  ') exit
  95
  96    end do head_info
  97
  98    ! 2.3 total rainfall data info
  99
 100    read (unit=char_file_rain (8:14),fmt='(I7)', iostat = iost) file_rain
 101
 102    ! 2.4 skip one lines
 103
 104    read (unit=iunit, fmt = '(A)', iostat = iost)
 105
 106    ! 3. read rain data
 107
 108    reports:   do n = 1, file_rain
 109
 110       ! 3.1 read station general info
 111
 112       read (unit = iunit, iostat = iost, &
 113          fmt = '(A12,1X,A19,1X,I6,2(F12.3,2X),F8.1,1X,A5)') &
 114          platform % info % platform,  &
 115          platform % info % date_char, &
 116          platform % info % levels,    &
 117          platform % info % lat,       &
 118          platform % info % lon,       &
 119          platform % info % elv,       &
 120          platform % info % id
 121
 122       if (print_detail_rain) then
 123          write(unit=stdout, fmt = '(A12,1X,A19,1X,I6,2(F12.3,2X),F8.1,1X,A5)') &
 124             platform % info % platform,  &
 125             platform % info % date_char, &
 126             platform % info % levels,    &
 127             platform % info % lat,       &
 128             platform % info % lon,       &
 129             platform % info % elv,       &
 130             platform % info % id
 131       end if
 132
 133       read(unit=platform % info % platform (4:6), fmt='(I3)') fm
 134
 135       num_report = num_report+1
 136
 137       ! 3.2 read rainfall data
 138
 139       platform%each (1) = rain_each_type(missing_r, missing, -1.0,&
 140          field_type(missing_r, missing, missing_r, missing, missing_r))
 141
 142       read (unit = iunit, fmt = '(F12.3,F12.3,I4,F12.3)') &
 143          platform % each (1) % height,             &
 144          platform % each (1) % rain % inv,         &
 145          platform % each (1) % rain % qc,          &
 146          platform % each (1) % rain % error
 147
 148      ! 3.3 Check if outside of the time range:
 149
 150       read (platform%info%date_char,'(i4,4(1x,i2))') &
 151             iyear, imonth, iday, ihour, imin
 152       call da_get_julian_time (iyear,imonth,iday,ihour,imin,obs_time)
 153       if ( obs_time < time_slots(0) .or. &
 154            obs_time >= time_slots(num_fgat_time) ) then
 155            num_outside_time = num_outside_time + 1
 156          cycle reports
 157       endif
 158
 159       call da_llxy (platform%info, platform%loc, outside, outside_all)
 160
 161       nlevels = platform%info%levels
 162
 163       if (outside_all) then
 164          num_outside_all = num_outside_all + 1
 165          cycle reports
 166       end if
 167
 168       dlat_earth = platform%info%lat
 169       dlon_earth = platform%info%lon
 170       if (dlon_earth < 0.0) dlon_earth = dlon_earth + 360.0
 171       if (dlon_earth >= 360.0) dlon_earth = dlon_earth - 360.0
 172       dlat_earth = dlat_earth * deg2rad
 173       dlon_earth = dlon_earth * deg2rad
 174
 175       ! 3.4  Loop over duplicating obs for global
 176       n_dup = 1
 177       if (global .and. &
 178          (platform%loc%i == ids .or. platform%loc%i == ide)) n_dup= 2
 179
 180       if (test_transforms) ndup = 1
 181
 182       do n_dup = 1, n_dup
 183          select case (fm)
 184
 185          case (129)
 186             if (iv%info(rain)%nlocal > max_rain_input) then
 187                write(unit=message(1),fmt='(A,I6,A,I6)') &
 188                   ' rain #= ',iv%info(rain)%nlocal, ' > max_rain_input = ', max_rain_input
 189                call da_error(__FILE__,__LINE__,message(1:1))
 190             end if
 191             iv%info(rain)%ntotal = iv%info(rain)%ntotal + 1
 192             if (outside) cycle reports
 193             if ( thin_rainobs ) then
 194                crit = 1.0
 195                call map2grids(rain,ifgat,dlat_earth,dlon_earth,crit,iv%info(rain)%nlocal,itx,1,itt,iout,iuse)
 196                if ( .not. iuse ) then
 197                   num_thinned = num_thinned + 1
 198                   cycle reports
 199                end if
 200             else
 201                iv%info(rain)%nlocal = iv%info(rain)%nlocal + 1
 202             endif
 203
 204          case default;
 205             write(unit=message(1), fmt='(a)') 'unsaved obs found:'
 206             write(unit=message(2), fmt='(2a)') &
 207                'platform%info%platform=', platform%info%platform
 208             write(unit=message(3), fmt='(a, i3)') &
 209                  'platform%info%levels=', platform%info%levels
 210             call da_warning(__FILE__,__LINE__,message(1:3))
 211             cycle reports
 212          end select
 213
 214          iv%info(rain)%max_lev = 1
 215
 216       end do        !  loop over duplicate
 217    end do reports
 218
 219    ! thinning check
 220    if ( thin_rainobs ) then
 221      if ( iv%info(rain)%ntotal > 0 ) then
 222
 223 #ifdef DM_PARALLEL
 224          ! Get minimum crit and associated processor index.
 225          allocate ( in  (thinning_grid(rain,ifgat)%itxmax) )
 226          allocate ( out (thinning_grid(rain,ifgat)%itxmax) )
 227          do i = 1, thinning_grid(rain,ifgat)%itxmax
 228             in(i) = thinning_grid(rain,ifgat)%score_crit(i)
 229          end do
 230
 231          call mpi_reduce(in, out, thinning_grid(rain,ifgat)%itxmax, true_mpi_real, mpi_min, root, comm, ierr)
 232          call wrf_dm_bcast_real (out, thinning_grid(rain,ifgat)%itxmax)
 233
 234          do i = 1, thinning_grid(rain,ifgat)%itxmax
 235             if ( abs(out(i)-thinning_grid(rain,ifgat)%score_crit(i)) > 1.0E-10 ) then
 236                thinning_grid(rain,ifgat)%ibest_obs(i) = 0
 237             end if
 238          end do
 239
 240          thinning_grid(rain,ifgat)%score_crit(:) = out(:)
 241
 242          deallocate( in  )
 243          deallocate( out )
 244 #endif
 245
 246          do j = (1+tp), iv%info(rain)%nlocal
 247             do i = 1, thinning_grid(rain,ifgat)%itxmax
 248                if ( thinning_grid(rain,ifgat)%ibest_obs(i) == j .and.         &
 249                     thinning_grid(rain,ifgat)%score_crit(i) < 9.99e6 ) then
 250                   nlocal = nlocal + 1
 251                   exit
 252                end if
 253             end do
 254          end do
 255
 256       num_thinned = num_thinned + iv%info(rain)%nlocal - nlocal
 257       iv%info(rain)%nlocal = tp + nlocal
 258       end if
 259    end if  ! thin_rainobs
 260
 261    write(unit=message(1),fmt='(A,4(1x,i7))') &
 262       'da_scan_obs_rain: num_report, num_outside_all, num_outside_time, num_thinned: ', &
 263       num_report, num_outside_all, num_outside_time, num_thinned
 264    call da_message(message(1:1))
 265
 266    close (iunit)
 267    call da_free_unit(iunit)
 268
 269    if (trace_use) call da_trace_exit("da_scan_obs_rain")
 270
 271
 272 end subroutine da_scan_obs_rain
 273
 274