var/da/da_sound/da_check_buddy_sound.inc

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