wrfv2_fire/var/da/da_airep/da_get_innov_vector_airep.inc

   1 subroutine da_get_innov_vector_airep( it,num_qcstat_conv, grid, ob, iv)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: TBD
   5    !    Updated for Analysis on Arakawa-C grid
   6    !    Author: Syed RH Rizvi,  MMM/ESSL/NCAR,  Date: 10/22/2008
   7    !-----------------------------------------------------------------------
   8
   9    implicit none
  10
  11    integer,          intent(in)    :: it       ! External iteration.
  12    type(domain),     intent(in)    :: grid     ! first guess state.
  13    type(y_type),     intent(inout) :: ob       ! Observation structure.
  14    type(iv_type),    intent(inout) :: iv       ! O-B structure.
  15    integer,          intent(inout) :: num_qcstat_conv(:,:,:,:)
  16
  17    integer :: n                         ! Loop counter.
  18    integer :: i  (kms:kme)
  19    integer :: j  (kms:kme)
  20    real    :: dx (kms:kme)
  21    real    :: dxm(kms:kme)
  22    real    :: dy (kms:kme)
  23    real    :: dym(kms:kme)
  24    integer :: k                   ! Index dimension.
  25
  26    real, allocatable :: model_u(:,:)  ! Model value u at ob location.
  27    real, allocatable :: model_v(:,:)  ! Model value v at ob location.
  28    real, allocatable :: model_t(:,:)  ! Model value t at ob location.
  29
  30    real    :: v_h(kms:kme)      ! Model value h at ob hor. location.
  31    real    :: v_p(kms:kme)      ! Model value p at ob hor. location.
  32
  33    if (trace_use_dull) call da_trace_entry("da_get_innov_vector_airep")
  34
  35    allocate (model_u(iv%info(airep)%max_lev,iv%info(airep)%n1:iv%info(airep)%n2))
  36    allocate (model_v(iv%info(airep)%max_lev,iv%info(airep)%n1:iv%info(airep)%n2))
  37    allocate (model_t(iv%info(airep)%max_lev,iv%info(airep)%n1:iv%info(airep)%n2))
  38
  39    model_u(:,:) = 0.0
  40    model_v(:,:) = 0.0
  41    model_t(:,:) = 0.0
  42
  43    if ( it > 1) then
  44       do n=iv%info(airep)%n1, iv%info(airep)%n2
  45          do k=1, iv%info(airep)%levels(n)
  46             if (iv%airep(n)%u(k)%qc == fails_error_max) iv%airep(n)%u(k)%qc = 0
  47             if (iv%airep(n)%v(k)%qc == fails_error_max) iv%airep(n)%v(k)%qc = 0
  48             if (iv%airep(n)%t(k)%qc == fails_error_max) iv%airep(n)%t(k)%qc = 0
  49          end do
  50       end do
  51    end if
  52
  53    do n=iv%info(airep)%n1, iv%info(airep)%n2
  54       if (iv%info(airep)%levels(n) < 1) cycle
  55
  56       ! [1.1] Get horizontal interpolation weights:
  57
  58       if (position_lev_dependant) then
  59          i(:)   = iv%info(airep)%i(:,n)
  60          j(:)   = iv%info(airep)%j(:,n)
  61          dx(:)  = iv%info(airep)%dx(:,n)
  62          dy(:)  = iv%info(airep)%dy(:,n)
  63          dxm(:) = iv%info(airep)%dxm(:,n)
  64          dym(:) = iv%info(airep)%dym(:,n)
  65          do k=kts,kte
  66             v_h(k) = dym(k)*(dxm(k)*grid%xb%h(i(k),j(k),  k) + dx(k)*grid%xb%h(i(k)+1,j(k),  k)) &
  67                     + dy(k)*(dxm(k)*grid%xb%h(i(k),j(k)+1,k) + dx(k)*grid%xb%h(i(k)+1,j(k)+1,k))
  68             v_p(k) = dym(k)*(dxm(k)*grid%xb%p(i(k),j(k),  k) + dx(k)*grid%xb%p(i(k)+1,j(k),  k)) &
  69                     + dy(k)*(dxm(k)*grid%xb%p(i(k),j(k)+1,k) + dx(k)*grid%xb%p(i(k)+1,j(k)+1,k))
  70          end do
  71       else
  72          i(1)   = iv%info(airep)%i(1,n)
  73          j(1)   = iv%info(airep)%j(1,n)
  74          dx(1)  = iv%info(airep)%dx(1,n)
  75          dy(1)  = iv%info(airep)%dy(1,n)
  76          dxm(1) = iv%info(airep)%dxm(1,n)
  77          dym(1) = iv%info(airep)%dym(1,n)
  78
  79          v_h(kts:kte) = dym(1)*(dxm(1)*grid%xb%h(i(1),j(1),  kts:kte) + dx(1)*grid%xb%h(i(1)+1,j(1),  kts:kte)) &
  80                        + dy(1)*(dxm(1)*grid%xb%h(i(1),j(1)+1,kts:kte) + dx(1)*grid%xb%h(i(1)+1,j(1)+1,kts:kte))
  81          v_p(kts:kte) = dym(1)*(dxm(1)*grid%xb%p(i(1),j(1),  kts:kte) + dx(1)*grid%xb%p(i(1)+1,j(1),  kts:kte)) &
  82                        + dy(1)*(dxm(1)*grid%xb%p(i(1),j(1)+1,kts:kte) + dx(1)*grid%xb%p(i(1)+1,j(1)+1,kts:kte))
  83       end if
  84
  85       do k=1, iv%info(airep)%levels(n)
  86          if (iv%airep(n)%p(k) > 1.0) then
  87             call da_to_zk (iv%airep(n)%p(k), v_p, v_interp_p, iv%info(airep)%zk(k,n))
  88          else if (iv%airep(n)%h(k) > 0.0) then
  89             call da_to_zk (iv%airep(n)%h(k), v_h, v_interp_h, iv%info(airep)%zk(k,n))
  90          end if
  91       end do
  92    end do
  93
  94    call da_convert_zk (iv%info(airep))
  95
  96    if (.not. anal_type_verify) then
  97       do n=iv%info(airep)%n1,iv%info(airep)%n2
  98          do k=1, iv%info(airep)%levels(n)
  99             if (iv%info(airep)%zk(k,n) < 0.0) then
 100                iv%airep(n)%u(k)%qc = missing_data
 101                iv%airep(n)%v(k)%qc = missing_data
 102                iv%airep(n)%t(k)%qc = missing_data
 103             end if
 104          end do
 105       end do
 106    end if
 107
 108 #ifdef A2C
 109    call da_interp_lin_3d (grid%xb%u, iv%info(airep), model_u,'u')
 110    call da_interp_lin_3d (grid%xb%v, iv%info(airep), model_v,'v')
 111 #else
 112    call da_interp_lin_3d (grid%xb%u, iv%info(airep), model_u)
 113    call da_interp_lin_3d (grid%xb%v, iv%info(airep), model_v)
 114 #endif
 115    call da_interp_lin_3d (grid%xb%t, iv%info(airep), model_t)
 116
 117    do n=iv%info(airep)%n1, iv%info(airep)%n2
 118
 119       !-------------------------------------------------------------------
 120       ! [2.0] Initialise components of innovation vector:
 121       !-------------------------------------------------------------------
 122
 123       do k = 1, iv%info(airep)%levels(n)
 124          iv % airep(n) % u(k) % inv = 0.0
 125          iv % airep(n) % v(k) % inv = 0.0
 126          iv % airep(n) % t(k) % inv = 0.0
 127
 128          !----------------------------------------------------------------
 129          ! [3.0] Fast interpolation:
 130          !----------------------------------------------------------------
 131
 132          if (ob % airep(n) % u(k) > missing_r .AND. iv % airep(n) % u(k) % qc >= obs_qc_pointer) then
 133             iv % airep(n) % u(k) % inv = ob % airep(n) % u(k) - model_u(k,n)
 134          end if
 135
 136          if (ob % airep(n) % v(k) > missing_r .AND. iv % airep(n) % v(k) % qc >= obs_qc_pointer) then
 137             iv % airep(n) % v(k) % inv = ob % airep(n) % v(k) - model_v(k,n)
 138          end if
 139
 140          if (ob % airep(n) % t(k) > missing_r .AND. iv % airep(n) % t(k) % qc >= obs_qc_pointer) then
 141             iv % airep(n) % t(k) % inv = ob % airep(n) % t(k) - model_t(k,n)
 142          end if
 143       end do
 144    end do
 145
 146    !-------------------------------------------------------------------
 147    ! [5.0] Perform optional maximum error check:
 148    !-------------------------------------------------------------------
 149
 150    if ( check_max_iv ) &
 151       call da_check_max_iv_airep (iv, it,num_qcstat_conv)
 152
 153    deallocate (model_u)
 154    deallocate (model_v)
 155    deallocate (model_t)
 156
 157    if (trace_use_dull) call da_trace_exit("da_get_innov_vector_airep")
 158
 159 end subroutine da_get_innov_vector_airep
 160
 161