var/da/da_radiance/da_qc_rad.inc

   1 subroutine da_qc_rad (it, ob, iv)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: perform quality control for radiance data.
   5    !
   6    ! METHOD:  separated QC for each sensor
   7    !---------------------------------------------------------------------------
   8
   9    implicit none
  10
  11    integer      ,  intent(in)      :: it         ! outer loop count
  12    type (y_type),  intent(in)      :: ob         ! Observation structure.
  13    type (iv_type), intent(inout)   :: iv         ! O-B structure.
  14
  15    integer :: i, nchan,p,j
  16    logical   :: amsua, amsub, hirs, msu,airs, hsb, ssmis, mhs, iasi, seviri
  17    logical   :: mwts, mwhs, atms, amsr2, imager, ahi, mwhs2, gmi, abi
  18
  19    integer, allocatable :: index(:)
  20    integer :: num_tovs_avg
  21    integer, allocatable :: excess_count(:)
  22    integer, allocatable :: spare_count(:)
  23    integer :: transfer
  24    logical :: copy_found
  25    integer :: temp(num_procs)
  26
  27    if (trace_use) call da_trace_entry("da_qc_rad")
  28
  29    if ( .not. allocated(num_tovs_before) )  allocate (num_tovs_before(iv%num_inst,num_procs))
  30    if ( .not. allocated(num_tovs_after) )   allocate (num_tovs_after(iv%num_inst,num_procs))
  31
  32    ! Cannot be more total send,receives than combination of processors
  33    if ( .not. allocated(tovs_copy_count) )  allocate (tovs_copy_count(iv%num_inst))
  34    if ( .not. allocated(tovs_send_pe) )     allocate (tovs_send_pe(iv%num_inst,num_procs*num_procs))
  35    if ( .not. allocated(tovs_recv_pe) )     allocate (tovs_recv_pe(iv%num_inst,num_procs*num_procs))
  36    if ( .not. allocated(tovs_send_start) )  allocate (tovs_send_start(iv%num_inst,num_procs*num_procs))
  37    if ( .not. allocated(tovs_send_count) )  allocate (tovs_send_count(iv%num_inst,num_procs*num_procs))
  38    if ( .not. allocated(tovs_recv_start) )  allocate (tovs_recv_start(iv%num_inst,num_procs*num_procs))
  39
  40    call da_trace("da_qc_rad", message="allocated tovs redistibution arrays")
  41
  42    if ( .not. allocated(index) )            allocate (index(num_procs))
  43    if ( .not. allocated(excess_count) )     allocate (excess_count(num_procs))
  44    if ( .not. allocated(spare_count) )      allocate (spare_count(num_procs))
  45
  46    do i = 1, iv%num_inst
  47
  48       !if (iv%instid(i)%info%n2 < iv%instid(i)%info%n1) cycle
  49
  50       nchan    = iv%instid(i)%nchan
  51
  52       amsua = trim(rttov_inst_name(rtminit_sensor(i))) == 'amsua'
  53       amsub = trim(rttov_inst_name(rtminit_sensor(i))) == 'amsub'
  54       hirs  = trim(rttov_inst_name(rtminit_sensor(i))) == 'hirs'
  55       msu   = trim(rttov_inst_name(rtminit_sensor(i))) == 'msu'
  56       airs  = trim(rttov_inst_name(rtminit_sensor(i))) == 'airs'
  57       hsb   = trim(rttov_inst_name(rtminit_sensor(i))) == 'hsb'
  58       ssmis = trim(rttov_inst_name(rtminit_sensor(i))) == 'ssmis'
  59       mhs   = trim(rttov_inst_name(rtminit_sensor(i))) == 'mhs'
  60       iasi  = trim(rttov_inst_name(rtminit_sensor(i))) == 'iasi'
  61       mwts  = trim(rttov_inst_name(rtminit_sensor(i))) == 'mwts'
  62       mwhs  = trim(rttov_inst_name(rtminit_sensor(i))) == 'mwhs'
  63       mwhs2 = trim(rttov_inst_name(rtminit_sensor(i))) == 'mwhs2'
  64       atms  = trim(rttov_inst_name(rtminit_sensor(i))) == 'atms'
  65       seviri = trim(rttov_inst_name(rtminit_sensor(i))) == 'seviri'
  66       amsr2 = trim(rttov_inst_name(rtminit_sensor(i))) == 'amsr2'
  67       imager = trim(rttov_inst_name(rtminit_sensor(i))) == 'imager'
  68          ahi = trim(rttov_inst_name(rtminit_sensor(i))) == 'ahi'
  69          abi = trim(rttov_inst_name(rtminit_sensor(i))) == 'abi'
  70       gmi   = trim(rttov_inst_name(rtminit_sensor(i))) == 'gmi'
  71       if (hirs) then
  72          ! 1.0 QC for HIRS
  73          call da_qc_hirs(it, i,nchan,ob,iv)
  74       else if (airs) then
  75          call da_qc_airs(it, i,nchan,ob,iv)
  76       else if ( hsb ) then
  77          ! call da_qc_hsb(it, i,nchan,ob,iv)
  78          call da_warning(__FILE__,__LINE__,(/'QC Not implemented for HSB'/))
  79       else if (amsua) then
  80          call da_qc_amsua(it,i,nchan,ob,iv)
  81       else if ( amsub ) then
  82          call da_qc_amsub(it,i,nchan,ob,iv)
  83       else if (msu) then
  84          ! call da_qc_msu(it, i,nchan, ob,iv)
  85          call da_warning(__FILE__,__LINE__,(/'QC Not implemented for MSU'/))
  86       else if (ssmis) then
  87          call da_qc_ssmis(it, i,nchan,ob,iv)
  88       else if (mhs) then
  89          call da_qc_mhs(it,i,nchan,ob,iv)
  90       else if (iasi) then
  91          call da_qc_iasi(it,i,nchan,ob,iv)
  92       else if (mwhs) then
  93          call da_qc_mwhs(it,i,nchan,ob,iv)
  94       else if (mwhs2) then
  95          call da_qc_mwhs2(it,i,nchan,ob,iv)
  96       else if (mwts) then
  97          call da_qc_mwts(it,i,nchan,ob,iv)
  98       else if (atms) then
  99          call da_qc_atms(it,i,nchan,ob,iv)
 100       else if (seviri) then
 101          call da_qc_seviri(it,i,nchan,ob,iv)
 102       else if (amsr2) then
 103          call da_qc_amsr2(it,i,nchan,ob,iv)
 104       else if (ahi) then
 105          call da_qc_ahi(it,i,nchan,ob,iv)
 106       else if (imager) then
 107          call da_qc_goesimg(it,i,nchan,ob,iv)
 108       else if (abi) then
 109          call da_qc_goesabi(it,i,nchan,ob,iv)
 110       else if (gmi) then
 111          call da_qc_gmi(it,i,nchan,ob,iv)
 112       else
 113          write(unit=message(1),fmt='(A,A)') &
 114             "Unrecognized instrument",trim(rttov_inst_name(rtminit_sensor(i)))
 115          call da_error(__FILE__,__LINE__,message(1:1))
 116       end if
 117
 118       ! Report number of observations to other processors via rootproc
 119
 120       num_tovs_before(i,:) = 0
 121       num_tovs_before(i,myproc+1)=iv%instid(i)%num_rad
 122       temp(:)= num_tovs_before(i,:)
 123       call da_proc_sum_ints(temp(:))
 124
 125 #ifdef DM_PARALLEL
 126       call wrf_dm_bcast_integer(temp(:),num_procs)
 127 #endif
 128       num_tovs_before(i,:) = temp(:)
 129
 130       num_tovs_after(i,:) = num_tovs_before(i,:)
 131
 132       if (rootproc .and. print_detail_rad) then
 133          write(unit=message(1),fmt='(A,I1,A)') "Instrument ",i, &
 134             " initial tovs distribution"
 135          write(unit=message(2),fmt=*) num_tovs_before(i,:)
 136          call da_message(message(1:2))
 137       end if
 138
 139       ! Decide how to reallocate observations
 140
 141       num_tovs_avg=sum(num_tovs_before(i,:))/num_procs
 142
 143       call da_trace_int_sort(num_tovs_before(i,:),num_procs,index)
 144
 145       do p=1,num_procs
 146          excess_count(p)=num_tovs_before(i,index(p))-num_tovs_avg
 147          spare_count(p)=num_tovs_avg-num_tovs_before(i,index(p))
 148       end do
 149
 150       tovs_copy_count(i) = 0
 151       tovs_send_start(i,:) = 0
 152       tovs_send_count(i,:) = 0
 153
 154       do
 155          copy_found = .false.
 156          do p=1,num_procs
 157             if (spare_count(p) > tovs_min_transfer) then
 158                do j=num_procs,1,-1
 159                   if (excess_count(j) > tovs_min_transfer) then
 160                      copy_found = .true.
 161                      tovs_copy_count(i)=tovs_copy_count(i)+1
 162                      tovs_send_pe(i,tovs_copy_count(i)) = index(j)-1
 163                      tovs_recv_pe(i,tovs_copy_count(i)) = index(p)-1
 164                      transfer=min(spare_count(p),excess_count(j))
 165                      tovs_send_count(i,tovs_copy_count(i)) = transfer
 166                      tovs_recv_start(i,tovs_copy_count(i)) = num_tovs_after(i,index(p))+1
 167                      num_tovs_after(i,index(p))=num_tovs_after(i,index(p))+transfer
 168                      num_tovs_after(i,index(j))=num_tovs_after(i,index(j))-transfer
 169                      tovs_send_start(i,tovs_copy_count(i)) = num_tovs_after(i,index(j))+1
 170                      spare_count(p)=spare_count(p)-transfer
 171                      excess_count(j)=excess_count(j)-transfer
 172                      exit
 173                   end if
 174                end do
 175             end if
 176          end do
 177          if (.not. copy_found) exit
 178       end do
 179
 180       if (print_detail_rad) then
 181          write(unit=message(1),fmt='(A,I1,A)') "Instrument ",i," final tovs distribution"
 182          write(unit=message(2),fmt=*) num_tovs_after(i,:)
 183          call da_message(message(1:2))
 184       end if
 185
 186       iv % instid(i) % num_rad_glo = sum(num_tovs_after(i,:))
 187    end do
 188
 189    deallocate (index)
 190    deallocate (excess_count)
 191    deallocate (spare_count)
 192
 193    if (trace_use) call da_trace_exit("da_qc_rad")
 194
 195 end subroutine da_qc_rad
 196
 197