var/da/da_varbc/da_varbc_precond.inc

   1   subroutine da_varbc_precond (iv)
   2
   3    !---------------------------------------------------------------------------
   4    !  PURPOSE:  Calculate covariance matrix b/w bias predictors
   5    !            for VarBC preconditioning
   6    !
   7    !  Called from da_get_innov_vector_radiance
   8    !
   9    !  HISTORY: 11/07/2007 - Creation                     Tom Auligne
  10    !---------------------------------------------------------------------------
  11
  12    implicit none
  13
  14    type (iv_type), intent (inout)   :: iv             ! Innovation
  15
  16    integer                          :: inst, n, i, j, k, ii, jj
  17    integer                          :: npred, npredmax, num_rad, num_rad_active
  18    real, allocatable                :: hessian(:,:)
  19    real*8, allocatable              :: eignvec(:,:), eignval(:)
  20    real                             :: hessian_local, bgerr_local, pred_i, pred_j
  21
  22    if ( iv%num_inst < 1 ) RETURN
  23
  24    if (trace_use) call da_trace_entry("da_varbc_precond")
  25
  26    write(unit=stdout,fmt='(A)') 'VARBC: Estimate Hessian for preconditioning'
  27
  28    do inst = 1, iv%num_inst      ! loop for sensors
  29       npredmax = iv%instid(inst)%varbc_info%npredmax
  30
  31       allocate ( hessian(npredmax, npredmax) )
  32       allocate ( eignvec(npredmax, npredmax) )
  33       allocate ( eignval(npredmax)           )
  34
  35       do k = 1, iv%instid(inst)%nchan         ! loop for channels
  36          npred = iv%instid(inst)%varbc(k)%npred
  37          if (npred   <= 0) cycle              ! VarBC channel only
  38
  39          num_rad        = iv%instid(inst)%num_rad
  40          if (num_rad > 0) then
  41             num_rad_active = COUNT( (iv%instid(inst)%info%proc_domain(1,1:num_rad)) &
  42                                .AND.(iv%instid(inst)%tb_qc(k,1:num_rad) > qc_varbc_bad) )
  43          else
  44             num_rad_active = 0
  45          end if
  46          iv%instid(inst)%varbc(k)%nobs  = wrf_dm_sum_integer(num_rad_active)
  47
  48          if ( satinfo(inst)%iuse(k) == 1 ) &
  49             write(unit=stdout,fmt='(A,I6,3A,I5)') &
  50             'VARBC:',iv%instid(inst)%varbc(k)%nobs,' active observations for ', &
  51             trim(adjustl(iv%instid(inst)%rttovid_string)),' channel',           &
  52             iv%instid(inst)%ichan(k)
  53
  54          if (iv%instid(inst)%varbc(k)%nobs == 0) cycle
  55
  56         !---------------------------------------------------------
  57         ! Calculate estimation of the Hessian for preconditioning
  58         !---------------------------------------------------------
  59          do i = 1, npred
  60             ii = iv%instid(inst)%varbc(k)%ipred(i)
  61
  62            ! Observation term
  63            !------------------
  64             do j = i, npred
  65                jj = iv%instid(inst)%varbc(k)%ipred(j)
  66                hessian_local = 0.0
  67
  68                do n= 1, num_rad      ! loop for pixel
  69                   if (iv%instid(inst)%tb_qc(k,n) <= qc_varbc_bad)   cycle  ! good obs only
  70                   if (.not. iv%instid(inst)%info%proc_domain(1,n))  cycle  ! do not sum up HALO data
  71
  72                   if (ii == iv%instid(inst)%varbc_info%gammapred) then
  73                      pred_i = iv%instid(inst)%gamma_jacobian(k,n)
  74                   else
  75                      pred_i = iv%instid(inst)%varbc_info%pred(ii,n)
  76                   end if
  77
  78                   if (jj == iv%instid(inst)%varbc_info%gammapred) then
  79                      pred_j = iv%instid(inst)%gamma_jacobian(k,n)
  80                   else
  81                      pred_j = iv%instid(inst)%varbc_info%pred(jj,n)
  82                   end if
  83
  84                   hessian_local = hessian_local + pred_i * pred_j / &
  85                                   iv%instid(inst)%tb_error(k,n)**2
  86                end do                                               !  end loop for pixel
  87
  88                ! Sum hessian preconditioning across processors
  89                hessian(i,j) = wrf_dm_sum_real(hessian_local)
  90                hessian(j,i) = hessian(i,j)
  91             end do
  92
  93            ! Background term
  94            !-----------------
  95             if (iv%instid(inst)%varbc_info%nbgerr(ii) <= 0) cycle
  96             bgerr_local = 0.0
  97             do n= 1, num_rad
  98                if (iv%instid(inst)%tb_qc(k,n) <= qc_varbc_bad)   cycle  ! good obs only
  99                if (.not. iv%instid(inst)%info%proc_domain(1,n))  cycle  ! do not sum up HALO data
 100
 101                bgerr_local = bgerr_local + iv%instid(inst)%tb_error(k,n)**2  / &
 102                                              varbc_nbgerr
 103                                            ! iv%instid(inst)%varbc_info%nbgerr(ii)
 104             end do
 105
 106             iv%instid(inst)%varbc(k)%bgerr(i) = wrf_dm_sum_real(bgerr_local) / &
 107                                                 iv%instid(inst)%varbc(k)%nobs
 108
 109             if (iv%instid(inst)%varbc(k)%bgerr(i) > 0) &
 110                hessian(i,i) = hessian(i,i) + 1/iv%instid(inst)%varbc(k)%bgerr(i)
 111          end do
 112
 113         !--------------------------------------------------
 114         ! Preconditioning = inverse square root of Hessian
 115         !--------------------------------------------------
 116          hessian = hessian / varbc_factor**2
 117
 118          if (npred == 1) then
 119             iv%instid(inst)%varbc(k)%vtox(1,1)     = 1.0 / sqrt(hessian(1,1))
 120          else
 121             call da_eof_decomposition(npred, hessian(1:npred,1:npred), &
 122                             eignvec(1:npred,1:npred),eignval(1:npred))
 123
 124             if (ANY(eignval(1:npred) <= 0)) then
 125                write(unit=stdout,fmt='(3A,I4,A,10F18.5)') &
 126                'VARBC: non-positive Hessian for ', iv%instid(inst)%rttovid_string, &
 127                       ' ,channel ',k,'--> Eigenvalues =',eignval(1:npred)
 128                do i = 1, npred
 129                   if (hessian(i,i) > 0) &
 130                      iv%instid(inst)%varbc(k)%vtox(i,i) = 1.0 / sqrt(hessian(i,i))
 131                end do
 132             else
 133                do i = 1, npred
 134                   do j = i, npred
 135                      iv%instid(inst)%varbc(k)%vtox(i,j) = sum(          &
 136                                            eignvec(i,1:npred)         * &
 137                                            sqrt(1.0/eignval(1:npred)) * &
 138                                            eignvec(j,1:npred) )
 139
 140                      iv%instid(inst)%varbc(k)%vtox(j,i) = &
 141                      iv%instid(inst)%varbc(k)%vtox(i,j)
 142                   end do
 143                end do
 144             end if
 145          end if
 146       end do                     !  end loop for channels
 147       deallocate(hessian, eignvec, eignval)
 148    end do                        !  end loop for sensor
 149
 150
 151    if (trace_use) call da_trace_exit("da_varbc_precond")
 152
 153  end subroutine da_varbc_precond