var/da/da_varbc_tamdar/da_varbc_tamdar_precond.inc

   1    subroutine da_varbc_tamdar_precond (iv)
   2
   3    !----------------------------------------------------!
   4    !  Estimate Hessian for Preconditioning              !
   5    !                                                    !
   6    !  Preconditioning = inverse square root of Hessian  !
   7    !                                                    !
   8    !  Hermite real matrix:                              !
   9    !                 A = Q*D*Q^T                        !
  10    !           A^(1/k) = Q*D^(1/K)*Q^T                  !
  11    !----------------------------------------------------!
  12
  13    implicit none
  14
  15    type (iv_type), intent (inout) :: iv
  16
  17    integer                        :: i,j,iflt,iobs,isn,ipred,iphase,npred
  18    real                           :: lhessian,lbgerr,predi,predj,verr
  19
  20    real, allocatable              :: hessian(:,:)
  21    real*8, allocatable            :: eignvec(:,:), eignval(:)
  22
  23
  24    if (trace_use) call da_trace_entry("da_varbc_tamdar_precond")
  25
  26    if (rootproc) &
  27        write(unit=varbc_tamdar_unit,fmt='(//5X,A)') 'Estimating hessian for preconditioning'
  28
  29    npred = iv%varbc_tamdar%npred
  30
  31    allocate ( hessian(npred, npred) )
  32    allocate ( eignvec(npred, npred) )
  33    allocate ( eignval(npred)        )
  34
  35    do iflt = 1, iv%varbc_tamdar%nair
  36       do iphase = 1, iv%varbc_tamdar%nphase
  37          if (iv%varbc_tamdar%nobs_sum(iphase,iflt) >= varbc_tamdar_nobsmin) then
  38
  39              hessian(:,:) = 0.
  40              eignvec(:,:) = 0.
  41              eignval(:)   = 0.
  42
  43              lbgerr = 0.0
  44              if (iv%varbc_tamdar%nobs(iphase,iflt) > 0 .and. iv%varbc_tamdar%ifuse(iphase,iflt) > 0) then
  45                  do iobs = 1, iv%varbc_tamdar%nobs(iphase,iflt)
  46                     isn = iv%varbc_tamdar%obs_sn(iobs,iphase,iflt)
  47                     verr = 0.0
  48                     if (iv%tamdar(isn)%t(1)%qc >= 0) verr = iv%tamdar(isn)%t(1)%error
  49                     lbgerr = lbgerr + verr**2.0/varbc_tamdar_nbgerr
  50                  end do
  51              end if
  52              lbgerr = wrf_dm_sum_real(lbgerr)
  53
  54              do i = 1, npred
  55                 do j = i, npred
  56                    lhessian = 0.0
  57
  58                    if (iv%varbc_tamdar%nobs(iphase,iflt) > 0 .and. iv%varbc_tamdar%ifuse(iphase,iflt) > 0) then
  59                        do iobs = 1, iv%varbc_tamdar%nobs(iphase,iflt)
  60                           isn = iv%varbc_tamdar%obs_sn(iobs,iphase,iflt)
  61                           predi = iv%varbc_tamdar%pred(i,iphase,iflt)
  62                           predj = iv%varbc_tamdar%pred(j,iphase,iflt)
  63
  64                           verr = 0.0
  65                           if (iv%tamdar(isn)%t(1)%qc >= 0) verr = iv%tamdar(isn)%t(1)%error
  66                           if (verr > 0.) lhessian = lhessian + predi*predj/verr**2.0
  67                        end do
  68                    end if
  69                    hessian(i,j) = wrf_dm_sum_real(lhessian)
  70                    hessian(j,i) = hessian(i,j)
  71                 end do
  72
  73                 iv%varbc_tamdar%bgerr(i,iphase,iflt) = lbgerr/iv%varbc_tamdar%nobs_sum(iphase,iflt)
  74                 if (iv%varbc_tamdar%bgerr(i,iphase,iflt) > 0.) &
  75                     hessian(i,i)=hessian(i,i)+1.0/iv%varbc_tamdar%bgerr(i,iphase,iflt)
  76              end do
  77
  78              if (iv%varbc_tamdar%nobs(iphase,iflt) > 0 .and. iv%varbc_tamdar%ifuse(iphase,iflt) > 0) then
  79
  80                  call da_eof_decomposition(npred, hessian(1:npred,1:npred), &
  81                                            eignvec(1:npred,1:npred),eignval(1:npred))
  82
  83                  if (ANY( eignval(1:npred) <= 0 )) then
  84                      write(unit=stdout,fmt='(A,I8,A,2F12.5)') &
  85                           'VARBC_TAMDAR: non-positive Hessian for tail_id ', iv%varbc_tamdar%tail_id(iflt), &
  86                           '. Eigenvalues =',eignval(1:npred)
  87
  88                      do i = 1, npred
  89                         if (hessian(i,i) > 0) iv%varbc_tamdar%vtox(i,i,iphase,iflt) = 1.0/sqrt(hessian(i,i))
  90                      end do
  91                  else
  92                      do i = 1, npred
  93                      do j = i, npred
  94                         iv%varbc_tamdar%vtox(i,j,iphase,iflt) = SUM( eignvec(i,1:npred)* &
  95                                                                      sqrt(1.0/eignval(1:npred)) * &
  96                                                                      eignvec(j,1:npred) )
  97                         iv%varbc_tamdar%vtox(j,i,iphase,iflt) = iv%varbc_tamdar%vtox(i,j,iphase,iflt)
  98                      end do
  99                      end do
 100                  end if
 101              end if
 102          end if
 103       end do
 104    end do
 105
 106    deallocate(hessian, eignvec, eignval)
 107
 108    if (rootproc) &
 109        write(unit=varbc_tamdar_unit,fmt='(/5X,A)') 'End estimating hessian for preconditioning'
 110
 111    if (trace_use) call da_trace_exit("da_varbc_tamdar_precond")
 112
 113    end subroutine da_varbc_tamdar_precond