wrfv2_fire/var/da/da_ssmi/da_get_innov_vector_ssmt1.inc

   1 subroutine da_get_innov_vector_ssmt1(it,num_qcstat_conv,grid, ob, iv)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: TBD
   5    !-----------------------------------------------------------------------
   6
   7    implicit none
   8
   9    integer,          intent(in)    :: it       ! External iteration.
  10    type(domain),     intent(in)    :: grid       ! first guess state.
  11    type(y_type),     intent(inout) :: ob       ! Observation structure.
  12    type(iv_type),    intent(inout) :: iv       ! O-B structure.
  13    integer,          intent(inout) :: num_qcstat_conv(:,:,:,:)
  14
  15
  16    integer :: n        ! Loop counter.
  17    integer :: i, j, k  ! Index dimension.
  18    integer :: num_levs ! Number of obs levels.
  19    real    :: dx, dxm  ! Interpolation weights.
  20    real    :: dy, dym  ! Interpolation weights.
  21    real, allocatable :: model_t(:,:)  ! Model value t at ob location.
  22
  23    real   :: v_h(kms:kme)      ! Model value h at ob hor. location.
  24    real   :: v_p(kms:kme)      ! Model value p at ob hor. location.
  25
  26    if (trace_use_dull) call da_trace_entry("da_get_innov_vector_ssmt1")
  27
  28    allocate (model_t(1:iv%info(ssmt1)%max_lev,iv%info(ssmt1)%n1:iv%info(ssmt1)%n2))
  29    model_t(:,:) = 0.0
  30
  31    if ( it > 1 )then
  32       do n=iv%info(ssmt1)%n1,iv%info(ssmt1)%n2
  33          do k = 1, iv%info(ssmt1)%levels(n)
  34             if(iv % ssmt1(n) % t(k) % qc == fails_error_max)iv % ssmt1(n) % t(k) % qc = 0
  35          end do
  36       end do
  37    end if
  38
  39    do n=iv%info(ssmt1)%n1,iv%info(ssmt1)%n2
  40
  41       num_levs = iv%info(ssmt1)%levels(n)
  42
  43       if (num_levs < 1) cycle
  44
  45       ! [1.1] Get horizontal interpolation weights:
  46
  47       i   = iv%info(ssmt1)%i(1,n)
  48       j   = iv%info(ssmt1)%j(1,n)
  49       dx  = iv%info(ssmt1)%dx(1,n)
  50       dy  = iv%info(ssmt1)%dy(1,n)
  51       dxm = iv%info(ssmt1)%dxm(1,n)
  52       dym = iv%info(ssmt1)%dym(1,n)
  53
  54       do k=kts,kte
  55          v_h(k) = dym*(dxm*grid%xb%h(i,j,k)+dx*grid%xb%h(i+1,j,k)) + dy*(dxm*grid%xb%h(i,j+1,k)+dx*grid%xb%h(i+1,j+1,k))
  56          v_p(k) = dym*(dxm*grid%xb%p(i,j,k)+dx*grid%xb%p(i+1,j,k)) + dy*(dxm*grid%xb%p(i,j+1,k)+dx*grid%xb%p(i+1,j+1,k))
  57       end do
  58
  59       num_levs=0
  60
  61       do k=1, iv%info(ssmt1)%levels(n)
  62          if (iv % ssmt1(n) % h(k) > 0.0) then
  63             call da_to_zk(iv % ssmt1(n) % h(k), v_h, v_interp_h, iv%info(ssmt1)%zk(k,n))
  64          else if (iv % ssmt1(n) % p(k) > 1.0) then
  65             call da_to_zk(iv % ssmt1(n) % p(k), v_p, v_interp_p, iv%info(ssmt1)%zk(k,n))
  66          end if
  67
  68          if ( iv%info(ssmt1)%zk(k,n) > 0.0) then
  69             num_levs=num_levs+1
  70             iv%info(ssmt1)%zk(num_levs,n)=iv%info(ssmt1)%zk(k,n)
  71
  72             ob % ssmt1(n) % t(num_levs)         = ob % ssmt1(n) % t(k)
  73             iv % ssmt1(n) % t(num_levs) % qc    = iv % ssmt1(n) % t(k) % qc
  74             iv % ssmt1(n) % t(num_levs) % error = iv % ssmt1(n) % t(k) % error
  75          end if
  76       end do
  77
  78       iv%info(ssmt1)%levels(n) = num_levs
  79    end do
  80
  81    call da_convert_zk (iv%info(ssmt1))
  82
  83    ! [1.2] Interpolate horizontally to ob:
  84
  85    call da_interp_lin_3d (grid%xb%t, iv%info(ssmt1), model_t)
  86
  87    do n=iv%info(ssmt1)%n1,iv%info(ssmt1)%n2
  88
  89       !---------------------------------------------------------------------
  90       ! [2.0] Initialise components of innovation vector:
  91       !---------------------------------------------------------------------
  92
  93       do k = 1, iv%info(ssmt1)%levels(n)
  94          iv % ssmt1(n) % t(k) % inv = 0.0
  95
  96          !----------------------------------------------------------------
  97          ! [3.0] Interpolation:
  98          !----------------------------------------------------------------
  99
 100          if (ob % ssmt1(n) % t(k) > missing_r .AND. &
 101              iv % ssmt1(n) % t(k) % qc >= obs_qc_pointer) then
 102             iv % ssmt1(n) % t(k) % inv = ob % ssmt1(n) % t(k) - model_t(k,n)
 103          end if
 104       end do
 105    end do
 106
 107    !------------------------------------------------------------------
 108    ! [5.0] Perform optional maximum error check:
 109    !------------------------------------------------------------------
 110
 111    if (check_max_iv) call da_check_max_iv_ssmt1(iv, it, num_qcstat_conv)
 112
 113    deallocate (model_t)
 114
 115    if (trace_use_dull) call da_trace_exit("da_get_innov_vector_ssmt1")
 116
 117 end subroutine da_get_innov_vector_ssmt1
 118
 119