var/da/da_synop/da_check_buddy_synop.inc

   1 subroutine da_check_buddy_synop(iv, ob, dx, it)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: Buddy check for SYNOP observations.
   5    !
   6    ! For SYNOP, there may not need the binning procedure before going
   7    ! into the da_buddy_qc. So  bottom_pressure = 30000.0 nad num_bins_p = 1.
   8    ! If you want to do binning, minor modifications needed.
   9    !
  10    !                       Yong-Run Guo, 10/10/2008
  11    !-----------------------------------------------------------------------
  12
  13    implicit none
  14
  15    type(iv_type), intent(inout) :: iv
  16    type(y_type),  intent(in)    :: ob      ! Observation structure
  17    integer,       intent(in)    :: it      ! Outer iteration
  18    real,          intent(in)    :: dx
  19
  20    integer :: k, n, bin, i, j, m_max, kk, nn, numobs
  21    real    :: dx_km, Press_mb
  22
  23 ! All data in one bin:
  24    integer, parameter               :: num_bins_p = 1
  25    real, parameter                  :: bottom_pressure = 30000.0
  26 !
  27 ! Total of 13 bins used:
  28 !   integer, parameter               :: num_bins_p = 13
  29 !   real, parameter                  :: bottom_pressure = 100000.0
  30
  31    real, parameter                  :: bin_width_p = 10000.0
  32    real   , dimension(0:num_bins_p) :: bin_start_p, pressure, bin_width
  33    integer, dimension(0:num_bins_p) :: num
  34 !
  35    integer, allocatable, dimension(:,:) :: n_iv
  36
  37    integer,          allocatable, dimension(:) :: qc_flag_small
  38    real,             allocatable, dimension(:) :: xob, yob, obs
  39    character(len=5), allocatable, dimension(:) :: station_id
  40 !-----------------------------------------------------------------------------
  41
  42 !   if (trace_use_dull) call da_trace_entry("da_check_buddy_synop")
  43
  44    !---------------------------------------------------------------------------
  45    ! [1.0] Open diagnostic file:
  46    !---------------------------------------------------------------------------
  47
  48    if (rootproc .and. check_buddy_print) then
  49       write (check_buddy_unit,'(/A)')  &
  50          '================================================================'
  51       write (unit = check_buddy_unit, fmt = '(A,i4,A,i4/)') &
  52             'SYNOP BUDDY TEST QC:  no_buddies_qc=',no_buddies,&
  53             '  fails_buddy_check_qc=',fails_buddy_check
  54    end if
  55
  56    !---------------------------------------------------------------------------
  57    ! [2.0] Bin the data vertically based on the obs p::
  58    !---------------------------------------------------------------------------
  59
  60 !   print*,'==> Synop Buddy check: num_bins_p = ',num_bins_p
  61    dx_km = dx / 1000.0
  62 !
  63    bin_start_p(0) = bottom_pressure
  64    pressure   (0) = bin_start_p(0)
  65    bin_width      (0) = 0.0
  66    do n = 1, num_bins_p
  67       bin_start_p(n) = bin_start_p(n-1) - bin_width(n-1)
  68       if (bin_start_p(n) > 30000.0) then
  69          bin_width(n) = bin_width_p
  70       else
  71          bin_width(n) = bin_width_p / 2.0
  72       endif
  73       pressure(n) = bin_start_p(n) - bin_width(n)/2.0
  74    enddo
  75    bin_start_p(0) = bottom_pressure + 10.0
  76 !
  77 ! Only 1 bin=0 used, if you want to do the normal binning, comment out
  78 ! the line below:
  79    pressure   (0) = 100000.0
  80
  81 !   print '(I3,2x,"start_p=",f10.1," mid-pressure=",f10.1," width=",f10.1)', &
  82 !        (n, bin_start_p(n), pressure(n), bin_width(n), n=0, num_bins_p)
  83 !
  84 ! 2.1 Get the maximum dimension for all the bins:
  85 !
  86    num = 0
  87    do n = iv%info(synop)%n1,iv%info(synop)%n2
  88          if (ob%synop(n)%p > missing_r) then
  89            do i = 0, num_bins_p - 1
  90               if (iv%synop(n)%p%qc >=0 .and.             &
  91                   (ob%synop(n)%p <= bin_start_p(i) .and. &
  92                    ob%synop(n)%p >  bin_start_p(i+1)) ) then
  93                  bin = i
  94                  exit
  95               endif
  96            enddo
  97 !           bin = int( (bottom_pressure - ob%synop(n)%p)/bin_width(n) ) + 1
  98            if (ob%synop(n)%p > bottom_pressure) bin = 0
  99            if (ob%synop(n)%p <=  bin_start_p(num_bins_p)) bin = num_bins_p
 100            num(bin) = num(bin) + 1
 101          endif
 102    enddo
 103    m_max = maxval(num)
 104 !   print *,(i,num(i),i=0,num_bins_p)
 105 !   print *,"m_max=", m_max
 106 !
 107 ! 2.2 Save the location indices (n,k) for each of bins:
 108 !
 109 !   print '("Synop n1=",i5,"  n2=",i5)',iv%info(synop)%n1,iv%info(synop)%n2
 110    allocate ( n_iv( 0: num_bins_p,1:m_max+10 ) )
 111
 112    num = 0
 113    do n = iv%info(synop)%n1,iv%info(synop)%n2
 114          if (ob%synop(n)%p > missing_r) then
 115            do i = 0, num_bins_p - 1
 116               if (iv%synop(n)%p%qc >=0 .and.             &
 117                   (ob%synop(n)%p <= bin_start_p(i) .and. &
 118                    ob%synop(n)%p >  bin_start_p(i+1)) ) then
 119                  bin = i
 120                  exit
 121               endif
 122            enddo
 123 !           bin = int( (bottom_pressure - ob%synop(n)%p)/bin_width(n) ) + 1
 124            if (ob%synop(n)%p > bottom_pressure) bin = 0
 125            if (ob%synop(n)%p <=  bin_start_p(num_bins_p)) bin = num_bins_p
 126
 127            num(bin) = num(bin) + 1
 128            n_iv(bin,num(bin)) = n
 129
 130          endif
 131    end do
 132 !
 133 ! 2.3 Print out the binned results:
 134 !
 135 !   do i = 0, num_bins_p
 136 !      print '("bin:",I2,"  start_p=",f8.1," num=",i5)', &
 137 !                      i, bin_start_p(i), num(i)
 138 !      do j = 1, num(i)
 139 !         n = n_iv(i,j)
 140 !         print '("j, n:",2i5,2x,"p=",f10.1)', &
 141 !                  j, n, ob%synop(n)%p
 142 !      enddo
 143 !   enddo
 144    !---------------------------------------------------------------------------
 145    ! [3.0] Buddy check for each of the pressure-bins::
 146    !---------------------------------------------------------------------------
 147
 148    do i = 0, num_bins_p
 149
 150      if (num(i) <= 1) cycle
 151
 152 ! 3.1 Get the Station locations:
 153
 154      ! Pressure level:
 155      Press_mb = pressure(i) / 100.0
 156      numobs = num(i)
 157
 158      allocate(xob(1:numobs))
 159      allocate(yob(1:numobs))
 160      allocate(obs(1:numobs))
 161      allocate(qc_flag_small(1:numobs))
 162      allocate(station_id   (1:numobs))
 163
 164      if (rootproc .and. check_buddy_print) then
 165          write (check_buddy_unit,'(5X,A)')  &
 166          '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
 167       write (unit = check_buddy_unit, fmt = '(5X,A,I3,2X,A,I6)') &
 168              'BIN =', i, 'NUMOBS =', numobs
 169      end if
 170 !     print *,'SYNOP: BIN=', i, '  numobs=',num(i)
 171
 172      ! Station locations
 173
 174      do n = 1, numobs
 175         nn = n_iv(i,n)
 176 !
 177         station_id(n)        = iv%info(synop)%id(nn)
 178         xob(n)               = iv%info(synop)%x(1,nn)
 179         yob(n)               = iv%info(synop)%y(1,nn)
 180      enddo
 181
 182 ! 3.2 U-component buddy check:
 183
 184      if (rootproc .and. check_buddy_print) &
 185          write (check_buddy_unit,'(8X,A,A,f10.1,3(A,f6.1))')  &
 186                 'UU      ', ' Pressure(mb)=',Press_mb, ' ds(km)=',dx_km,&
 187                 '  buddy_weight=', buddy_weight , &
 188                 '  max_buddy_uv=', max_buddy_uv
 189
 190      obs = 0.0
 191      qc_flag_small(n) = missing
 192      do n = 1, numobs
 193         nn = n_iv(i,n)
 194         obs(n)               = iv%synop(nn)%u%inv
 195         qc_flag_small(n)     = iv%synop(nn)%u%qc
 196      enddo
 197
 198      call da_buddy_qc (numobs, m_max, station_id, xob, yob, obs, qc_flag_small,&
 199                        'UU      ', Press_mb, dx_km, buddy_weight , &
 200                        max_buddy_uv , check_buddy_unit, check_buddy_print )
 201
 202    !  Put the qc_flag back into the permanent space.
 203
 204      do n = 1, numobs
 205         nn = n_iv(i,n)
 206         iv%synop(nn)%u%qc = qc_flag_small(n)
 207      enddo
 208
 209 ! 3.2 V-component buddy check:
 210
 211      if (rootproc .and. check_buddy_print) &
 212          write (check_buddy_unit,'(8X,A,A,f10.1,3(A,f6.1))')  &
 213                 'VV      ', ' Pressure(mb)=',Press_mb, ' ds(km)=',dx_km,&
 214                 '  buddy_weight=', buddy_weight , &
 215                 '  max_buddy_uv=', max_buddy_uv
 216
 217
 218      obs = 0.0
 219      qc_flag_small(n) = missing
 220      do n = 1, numobs
 221         nn = n_iv(i,n)
 222         obs(n)               = iv%synop(nn)%v%inv
 223         qc_flag_small(n)     = iv%synop(nn)%v%qc
 224      enddo
 225
 226      call da_buddy_qc (numobs, m_max, station_id, xob, yob, obs, qc_flag_small,&
 227                        'VV      ', Press_mb, dx_km, buddy_weight , &
 228                        max_buddy_uv , check_buddy_unit, check_buddy_print )
 229
 230    !  Put the qc_flag back into the permanent space.
 231
 232      do n = 1, numobs
 233         nn = n_iv(i,n)
 234         iv%synop(nn)%v%qc = qc_flag_small(n)
 235      enddo
 236
 237 ! 3.3 Temperature buddy check:
 238
 239      if (rootproc .and. check_buddy_print) &
 240          write (check_buddy_unit,'(8X,A,A,f10.1,3(A,f6.1))')  &
 241                 'TT      ', ' Pressure(mb)=',Press_mb, ' ds(km)=',dx_km,&
 242                 '  buddy_weight=', buddy_weight , &
 243                 '  max_buddy_t=', max_buddy_t
 244
 245      obs = 0.0
 246      qc_flag_small(n) = missing
 247      do n = 1, numobs
 248         nn = n_iv(i,n)
 249         obs(n)               = iv%synop(nn)%t%inv
 250         qc_flag_small(n)     = iv%synop(nn)%t%qc
 251      enddo
 252
 253      call da_buddy_qc (numobs, m_max, station_id, xob, yob, obs, qc_flag_small,&
 254                        'TT      ', Press_mb, dx_km, buddy_weight , &
 255                        max_buddy_t , check_buddy_unit, check_buddy_print )
 256
 257    !  Put the qc_flag back into the permanent space.
 258
 259      do n = 1, numobs
 260         nn = n_iv(i,n)
 261         iv%synop(nn)%t%qc = qc_flag_small(n)
 262      enddo
 263
 264 ! 3.3 Specific humidity buddy check:
 265
 266      if (rootproc .and. check_buddy_print) &
 267          write (check_buddy_unit,'(8X,A,A,f10.1,3(A,f6.1))')  &
 268                 'QQ      ', ' Pressure(mb)=',Press_mb, ' ds(km)=',dx_km,&
 269                 '  buddy_weight=', buddy_weight , &
 270                 '  max_buddy_rh=', max_buddy_rh
 271
 272      obs = 0.0
 273      qc_flag_small(n) = missing
 274      do n = 1, numobs
 275         nn = n_iv(i,n)
 276         obs(n)               = iv%synop(nn)%q%inv
 277         qc_flag_small(n)     = iv%synop(nn)%q%qc
 278      enddo
 279
 280      call da_buddy_qc (numobs, m_max, station_id, xob, yob, obs, qc_flag_small,&
 281                        'QQ      ', Press_mb, dx_km, buddy_weight , &
 282                        max_buddy_rh , check_buddy_unit, check_buddy_print )
 283
 284    !  Put the qc_flag back into the permanent space.
 285
 286      do n = 1, numobs
 287         nn = n_iv(i,n)
 288         iv%synop(nn)%q%qc = qc_flag_small(n)
 289      enddo
 290
 291 ! 3.4 Pressure buddy check:
 292
 293      if (rootproc .and. check_buddy_print) &
 294          write (check_buddy_unit,'(8X,A,A,f10.1,3(A,f6.1))')  &
 295                 'PMSL    ', ' Pressure(mb)=',Press_mb, ' ds(km)=',dx_km,&
 296                 '  buddy_weight=', buddy_weight , &
 297                 '  max_buddy_p=', max_buddy_p
 298
 299      obs = 0.0
 300      qc_flag_small(n) = missing
 301      do n = 1, numobs
 302         nn = n_iv(i,n)
 303         obs(n)               = iv%synop(nn)%p%inv
 304         qc_flag_small(n)     = iv%synop(nn)%p%qc
 305      enddo
 306
 307      call da_buddy_qc (numobs, m_max, station_id, xob, yob, obs, qc_flag_small,&
 308                        'PMSL    ', Press_mb, dx_km, buddy_weight , &
 309                        max_buddy_p , check_buddy_unit, check_buddy_print )
 310
 311    !  Put the qc_flag back into the permanent space.
 312
 313      do n = 1, numobs
 314         nn = n_iv(i,n)
 315         iv%synop(nn)%p%qc = qc_flag_small(n)
 316      enddo
 317
 318 ! 3.5 Deallocate arrays
 319
 320 1234 continue
 321
 322      deallocate(xob)
 323      deallocate(yob)
 324      deallocate(obs)
 325      deallocate(qc_flag_small)
 326      deallocate(station_id   )
 327
 328    enddo
 329
 330    deallocate ( n_iv )
 331
 332    if (trace_use_dull) call da_trace_exit("da_check_buddy_synop")
 333
 334 end subroutine da_check_buddy_synop