var/da/da_pilot/da_get_innov_vector_pilot.inc

   1 subroutine da_get_innov_vector_pilot( 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, j, k  ! Index dimension.
  19
  20    real    :: dx, dxm  ! Interpolation weights.
  21    real    :: dy, dym  ! Interpolation weights.
  22
  23    real, allocatable :: model_u(:,:)  ! Model value u at ob location.
  24    real, allocatable :: model_v(:,:)  ! Model value v at ob location.
  25
  26    real    :: v_h(kms:kme)      ! Model value h at ob hor. location.
  27    real    :: v_p(kms:kme)      ! Model value p at ob hor. location.
  28    real    :: speed, direction
  29
  30    if (trace_use_dull) call da_trace_entry("da_get_innov_vector_pilot")
  31
  32    allocate (model_u(iv%info(pilot)%max_lev,iv%info(pilot)%n1:iv%info(pilot)%n2))
  33    allocate (model_v(iv%info(pilot)%max_lev,iv%info(pilot)%n1:iv%info(pilot)%n2))
  34
  35    model_u(:,:) = 0.0
  36    model_v(:,:) = 0.0
  37
  38    if ( it > 1 ) then
  39       do n=iv%info(pilot)%n1,iv%info(pilot)%n2
  40          do k=1, iv%info(pilot)%levels(n)
  41             if (iv%pilot(n)%u(k)%qc == fails_error_max) iv%pilot(n)%u(k)%qc = 0
  42             if (iv%pilot(n)%v(k)%qc == fails_error_max) iv%pilot(n)%v(k)%qc = 0
  43          end do
  44       end do
  45    end if
  46
  47    do n=iv%info(pilot)%n1,iv%info(pilot)%n2
  48       ! [1.3] Get horizontal interpolation weights:
  49
  50       i   = iv%info(pilot)%i(1,n)
  51       j   = iv%info(pilot)%j(1,n)
  52       dx  = iv%info(pilot)%dx(1,n)
  53       dy  = iv%info(pilot)%dy(1,n)
  54       dxm = iv%info(pilot)%dxm(1,n)
  55       dym = iv%info(pilot)%dym(1,n)
  56
  57       do k=kts,kte
  58          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))
  59          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))
  60       end do
  61
  62       do k=1, iv%info(pilot)%levels(n)
  63          if (iv % pilot(n) % p(k) > 1.0) then
  64             call da_to_zk(iv % pilot(n) % p(k), v_p, v_interp_p, iv%info(pilot)%zk(k,n))
  65          else if (iv%pilot(n) % h(k) > missing_r) then
  66             call da_to_zk(iv % pilot(n) % h(k), v_h, v_interp_h, iv%info(pilot)%zk(k,n))
  67          end if
  68
  69          if (iv%info(pilot)%zk(k,n) < 0.0 .and.  .not.anal_type_verify) then
  70             iv % pilot(n) % u(k) % qc = missing_data
  71             iv % pilot(n) % v(k) % qc = missing_data
  72          end if
  73       end do
  74    end do
  75
  76    call da_convert_zk (iv%info(pilot))
  77
  78    ! [1.4] Interpolate horizontally:
  79 #ifdef A2C
  80    call da_interp_lin_3d (grid%xb%u, iv%info(pilot), model_u,'u')
  81    call da_interp_lin_3d (grid%xb%v, iv%info(pilot), model_v,'v')
  82 #else
  83    call da_interp_lin_3d (grid%xb%u, iv%info(pilot), model_u)
  84    call da_interp_lin_3d (grid%xb%v, iv%info(pilot), model_v)
  85 #endif
  86
  87    do n=iv%info(pilot)%n1,iv%info(pilot)%n2
  88       !------------------------------------------------------------------------
  89       ! [2.0] Initialise components of innovation vector:
  90       !------------------------------------------------------------------------
  91
  92       do k = 1, iv%info(pilot)%levels(n)
  93          iv % pilot(n) % u(k) % inv = 0.0
  94          iv % pilot(n) % v(k) % inv = 0.0
  95
  96          !------------------------------------------------------------------------
  97          ! [4.0] Fast interpolation:
  98          !------------------------------------------------------------------------
  99
 100
 101           if (wind_sd_pilot) then
 102               call da_ffdduv_model (speed,direction,model_u(k,n), model_v(k,n), convert_uv2fd)
 103
 104               if (ob%pilot(n)%u(k) > missing_r .AND. iv%pilot(n)%u(k)%qc >= obs_qc_pointer) then
 105                   iv%pilot(n)%u(k)%inv = ob%pilot(n)%u(k) - speed
 106               end if
 107
 108               if (ob%pilot(n)%v(k) > missing_r .AND. iv%pilot(n)%v(k)%qc >= obs_qc_pointer) then
 109                   iv%pilot(n)%v(k)%inv = ob%pilot(n)%v(k) - direction
 110                   if (iv%pilot(n)%v(k)%inv > 180.0 ) iv%pilot(n)%v(k)%inv = iv%pilot(n)%v(k)%inv - 360.0
 111                   if (iv%pilot(n)%v(k)%inv < -180.0 ) iv%pilot(n)%v(k)%inv = iv%pilot(n)%v(k)%inv + 360.0
 112               end if
 113           else
 114              if (ob%pilot(n)%u(k) > missing_r .AND. iv%pilot(n)%u(k)%qc >= obs_qc_pointer) then
 115                  iv%pilot(n)%u(k)%inv = ob%pilot(n)%u(k) - model_u(k,n)
 116              end if
 117
 118              if (ob%pilot(n)%v(k) > missing_r .AND. iv%pilot(n)%v(k)%qc >= obs_qc_pointer) then
 119                  iv%pilot(n)%v(k)%inv = ob%pilot(n)%v(k) - model_v(k,n)
 120              end if
 121           end if
 122
 123       end do
 124    end do
 125
 126    !------------------------------------------------------------------------
 127    ! [5.0] Perform optional maximum error check:
 128    !------------------------------------------------------------------------
 129
 130    if ( check_max_iv ) &
 131       call da_check_max_iv_pilot(iv, it, num_qcstat_conv)
 132
 133    deallocate (model_u)
 134    deallocate (model_v)
 135
 136    if (trace_use_dull) call da_trace_exit("da_get_innov_vector_pilot")
 137
 138 end subroutine da_get_innov_vector_pilot
 139
 140