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
  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       gmi   = trim(rttov_inst_name(rtminit_sensor(i))) == 'gmi'
  70       if (hirs) then
  71          ! 1.0 QC for HIRS
  72          call da_qc_hirs(it, i,nchan,ob,iv)
  73       else if (airs) then
  74          call da_qc_airs(it, i,nchan,ob,iv)
  75       else if ( hsb ) then
  76          ! call da_qc_hsb(it, i,nchan,ob,iv)
  77          call da_warning(__FILE__,__LINE__,(/'QC Not implemented for HSB'/))
  78       else if (amsua) then
  79          call da_qc_amsua(it,i,nchan,ob,iv)
  80       else if ( amsub ) then
  81          call da_qc_amsub(it,i,nchan,ob,iv)
  82       else if (msu) then
  83          ! call da_qc_msu(it, i,nchan, ob,iv)
  84          call da_warning(__FILE__,__LINE__,(/'QC Not implemented for MSU'/))
  85       else if (ssmis) then
  86          call da_qc_ssmis(it, i,nchan,ob,iv)
  87       else if (mhs) then
  88          call da_qc_mhs(it,i,nchan,ob,iv)
  89       else if (iasi) then
  90          call da_qc_iasi(it,i,nchan,ob,iv)
  91       else if (mwhs) then
  92          call da_qc_mwhs(it,i,nchan,ob,iv)
  93       else if (mwhs2) then
  94          call da_qc_mwhs2(it,i,nchan,ob,iv)
  95       else if (mwts) then
  96          call da_qc_mwts(it,i,nchan,ob,iv)
  97       else if (atms) then
  98          call da_qc_atms(it,i,nchan,ob,iv)
  99       else if (seviri) then
 100          call da_qc_seviri(it,i,nchan,ob,iv)
 101       else if (amsr2) then
 102          call da_qc_amsr2(it,i,nchan,ob,iv)
 103       else if (ahi) then
 104          call da_qc_ahi(it,i,nchan,ob,iv)
 105       else if (imager) then
 106          call da_qc_goesimg(it,i,nchan,ob,iv)
 107       else if (gmi) then
 108          call da_qc_gmi(it,i,nchan,ob,iv)
 109       else
 110          write(unit=message(1),fmt='(A,A)') &
 111             "Unrecognized instrument",trim(rttov_inst_name(rtminit_sensor(i)))
 112          call da_error(__FILE__,__LINE__,message(1:1))
 113       end if
 114
 115       ! Report number of observations to other processors via rootproc
 116
 117       num_tovs_before(i,:) = 0
 118       num_tovs_before(i,myproc+1)=iv%instid(i)%num_rad
 119       temp(:)= num_tovs_before(i,:)
 120       call da_proc_sum_ints(temp(:))
 121
 122 #ifdef DM_PARALLEL
 123       call wrf_dm_bcast_integer(temp(:),num_procs)
 124 #endif
 125       num_tovs_before(i,:) = temp(:)
 126
 127       num_tovs_after(i,:) = num_tovs_before(i,:)
 128
 129       if (rootproc .and. print_detail_rad) then
 130          write(unit=message(1),fmt='(A,I1,A)') "Instrument ",i, &
 131             " initial tovs distribution"
 132          write(unit=message(2),fmt=*) num_tovs_before(i,:)
 133          call da_message(message(1:2))
 134       end if
 135
 136       ! Decide how to reallocate observations
 137
 138       num_tovs_avg=sum(num_tovs_before(i,:))/num_procs
 139
 140       call da_trace_int_sort(num_tovs_before(i,:),num_procs,index)
 141
 142       do p=1,num_procs
 143          excess_count(p)=num_tovs_before(i,index(p))-num_tovs_avg
 144          spare_count(p)=num_tovs_avg-num_tovs_before(i,index(p))
 145       end do
 146
 147       tovs_copy_count(i) = 0
 148       tovs_send_start(i,:) = 0
 149       tovs_send_count(i,:) = 0
 150
 151       do
 152          copy_found = .false.
 153          do p=1,num_procs
 154             if (spare_count(p) > tovs_min_transfer) then
 155                do j=num_procs,1,-1
 156                   if (excess_count(j) > tovs_min_transfer) then
 157                      copy_found = .true.
 158                      tovs_copy_count(i)=tovs_copy_count(i)+1
 159                      tovs_send_pe(i,tovs_copy_count(i)) = index(j)-1
 160                      tovs_recv_pe(i,tovs_copy_count(i)) = index(p)-1
 161                      transfer=min(spare_count(p),excess_count(j))
 162                      tovs_send_count(i,tovs_copy_count(i)) = transfer
 163                      tovs_recv_start(i,tovs_copy_count(i)) = num_tovs_after(i,index(p))+1
 164                      num_tovs_after(i,index(p))=num_tovs_after(i,index(p))+transfer
 165                      num_tovs_after(i,index(j))=num_tovs_after(i,index(j))-transfer
 166                      tovs_send_start(i,tovs_copy_count(i)) = num_tovs_after(i,index(j))+1
 167                      spare_count(p)=spare_count(p)-transfer
 168                      excess_count(j)=excess_count(j)-transfer
 169                      exit
 170                   end if
 171                end do
 172             end if
 173          end do
 174          if (.not. copy_found) exit
 175       end do
 176
 177       if (print_detail_rad) then
 178          write(unit=message(1),fmt='(A,I1,A)') "Instrument ",i," final tovs distribution"
 179          write(unit=message(2),fmt=*) num_tovs_after(i,:)
 180          call da_message(message(1:2))
 181       end if
 182
 183       iv % instid(i) % num_rad_glo = sum(num_tovs_after(i,:))
 184    end do
 185
 186    deallocate (index)
 187    deallocate (excess_count)
 188    deallocate (spare_count)
 189
 190    if (trace_use) call da_trace_exit("da_qc_rad")
 191
 192 end subroutine da_qc_rad
 193
 194