var/da/da_geoamv/da_get_innov_vector_geoamv.inc

   1 subroutine da_get_innov_vector_geoamv( it,num_qcstat_conv,  grid, ob, iv)
   2
   3    !-------------------------------------------------------------------------
   4    ! Purpose: Calculates innovation vector does QC for geoamv
   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(in)    :: 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 :: k              ! Index dimension.
  19    integer :: num_levs       ! Number of obs levels.
  20    real    :: speed, direction
  21
  22    integer :: i  (kms:kme)
  23    integer :: j  (kms:kme)
  24    real    :: dx (kms:kme)
  25    real    :: dxm(kms:kme)
  26    real    :: dy (kms:kme)
  27    real    :: dym(kms:kme)
  28    real,allocatable :: model_u(:,:)
  29    real,allocatable :: model_v(:,:)
  30
  31    real    :: v_p(kts:kte)      ! Model value p at ob hor. location.
  32
  33    integer :: itu,ituf,itvv,itvvf
  34
  35    if (trace_use_dull) call da_trace_entry ("da_get_innov_vector_geoamv")
  36
  37    allocate (model_u(iv%info(geoamv)%max_lev,iv%info(geoamv)%n1:iv%info(geoamv)%n2))
  38    allocate (model_v(iv%info(geoamv)%max_lev,iv%info(geoamv)%n1:iv%info(geoamv)%n2))
  39
  40    model_u(:,:) = 0.0
  41    model_v(:,:) = 0.0
  42
  43    if ( it > 1 ) then
  44       do n = iv%info(geoamv)%n1, iv%info(geoamv)%n2
  45          do k = 1, iv%info(geoamv)%levels(n)
  46             if (iv%geoamv(n)%u(k)%qc == fails_error_max) iv%geoamv(n)%u(k)%qc = 0
  47             if (iv%geoamv(n)%v(k)%qc == fails_error_max) iv%geoamv(n)%v(k)%qc = 0
  48          end do
  49       end do
  50    end if
  51
  52    do n = iv%info(geoamv)%n1, iv%info(geoamv)%n2
  53       ! [1.3] Get horizontal interpolation weights:
  54
  55       num_levs = iv%info(geoamv)%levels(n)
  56       if (num_levs < 1) cycle
  57
  58       ! slower
  59       ! i(:)   = iv%info(geoamv)%i(:,n)
  60       ! j(:)   = iv%info(geoamv)%j(:,n)
  61       ! dx(:)  = iv%info(geoamv)%dx(:,n)
  62       ! dy(:)  = iv%info(geoamv)%dy(:,n)
  63       ! dxm(:) = iv%info(geoamv)%dxm(:,n)
  64       ! dym(:) = iv%info(geoamv)%dym(:,n)
  65
  66       ! faster
  67       i(1)   = iv%info(geoamv)%i(1,n)
  68       j(1)   = iv%info(geoamv)%j(1,n)
  69       dx(1)  = iv%info(geoamv)%dx(1,n)
  70       dy(1)  = iv%info(geoamv)%dy(1,n)
  71       dxm(1) = iv%info(geoamv)%dxm(1,n)
  72       dym(1) = iv%info(geoamv)%dym(1,n)
  73
  74       ! if position varies with height, slower
  75       ! do k=kts,kte
  76       !    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)) &
  77       !       + 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))
  78       ! end do
  79
  80       ! If position does not, faster
  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
  84       do k=1, iv%info(geoamv)%levels(n)
  85          if (iv%geoamv(n)%p(k) > 1.0) then
  86             call da_to_zk (iv%geoamv(n)%p(k), v_p, v_interp_p, iv%info(geoamv)%zk(k,n))
  87          end if
  88       end do
  89
  90    end do
  91
  92    call da_convert_zk (iv%info(geoamv))
  93
  94    if (.not. anal_type_verify) then
  95       do n = iv%info(geoamv)%n1, iv%info(geoamv)%n2
  96          do k=1, iv%info(geoamv)%levels(n)
  97             if (iv%info(geoamv)%zk(k,n) < 0.0) then
  98                iv%geoamv(n)%u(k)% qc = missing_data
  99                iv%geoamv(n)%v(k)% qc = missing_data
 100             end if
 101          end do
 102       end do
 103    end if
 104
 105 #ifdef A2C
 106    call da_interp_lin_3d (grid%xb%u, iv%info(geoamv), model_u,'u')
 107    call da_interp_lin_3d (grid%xb%v, iv%info(geoamv), model_v,'v')
 108 #else
 109    call da_interp_lin_3d (grid%xb%u, iv%info(geoamv), model_u)
 110    call da_interp_lin_3d (grid%xb%v, iv%info(geoamv), model_v)
 111 #endif
 112
 113    do n = iv%info(geoamv)%n1, iv%info(geoamv)%n2
 114       do k = 1, iv%info(geoamv)%levels(n)
 115          iv%geoamv(n)%u(k)%inv = 0.0
 116          iv%geoamv(n)%v(k)%inv = 0.0
 117
 118          if (wind_sd_geoamv) then
 119              call da_ffdduv_model (speed,direction,model_u(k,n), model_v(k,n), convert_uv2fd)
 120
 121              if (ob%geoamv(n)%u(k) > missing_r .AND. iv%geoamv(n)%u(k)%qc >= obs_qc_pointer) then
 122                  iv%geoamv(n)%u(k)%inv = ob%geoamv(n)%u(k) - speed
 123              end if
 124
 125              if (ob%geoamv(n)%v(k) > missing_r .AND. iv%geoamv(n)%v(k)%qc >= obs_qc_pointer) then
 126                  iv%geoamv(n)%v(k)%inv = ob%geoamv(n)%v(k) - direction
 127                  if (iv%geoamv(n)%v(k)%inv > 180.0 ) iv%geoamv(n)%v(k)%inv = iv%geoamv(n)%v(k)%inv - 360.0
 128                  if (iv%geoamv(n)%v(k)%inv < -180.0 ) iv%geoamv(n)%v(k)%inv = iv%geoamv(n)%v(k)%inv + 360.0
 129              end if
 130          else
 131              if (ob%geoamv(n)%u(k) > missing_r .AND. iv%geoamv(n)%u(k)%qc >= obs_qc_pointer) then
 132                  iv%geoamv(n)%u(k)%inv = ob%geoamv(n)%u(k) - model_u(k,n)
 133              end if
 134
 135              if (ob%geoamv(n)%v(k) > missing_r .AND. iv%geoamv(n)%v(k)%qc >= obs_qc_pointer) then
 136                  iv%geoamv(n)%v(k)%inv = ob%geoamv(n)%v(k) - model_v(k,n)
 137              end if
 138          end if
 139
 140       end do
 141    end do
 142
 143    !------------------------------------------------------------------------
 144    ! Perform optional maximum error check:
 145    !------------------------------------------------------------------------
 146    if ( check_max_iv ) &
 147       call da_check_max_iv_geoamv(iv,it,num_qcstat_conv)
 148
 149    deallocate (model_u)
 150    deallocate (model_v)
 151
 152    if (trace_use_dull) call da_trace_exit ("da_get_innov_vector_geoamv")
 153
 154 end subroutine da_get_innov_vector_geoamv