var/da/da_mtgirs/da_get_innov_vector_mtgirs.inc

   1 subroutine da_get_innov_vector_mtgirs (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, k        ! 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    real    :: speed, direction
  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    real, allocatable :: model_q(:,:)  ! Model value q at ob location.
  30
  31    real    :: v_h(kts:kte)      ! Model value h at ob hor. location.
  32    real    :: v_p(kts:kte)      ! Model value p at ob hor. location.
  33
  34
  35    if (trace_use_dull) call da_trace_entry ("da_get_innov_vector_mtgirs")
  36
  37    allocate (model_u(iv%info(mtgirs)%max_lev,iv%info(mtgirs)%n1:iv%info(mtgirs)%n2))
  38    allocate (model_v(iv%info(mtgirs)%max_lev,iv%info(mtgirs)%n1:iv%info(mtgirs)%n2))
  39    allocate (model_t(iv%info(mtgirs)%max_lev,iv%info(mtgirs)%n1:iv%info(mtgirs)%n2))
  40    allocate (model_q(iv%info(mtgirs)%max_lev,iv%info(mtgirs)%n1:iv%info(mtgirs)%n2))
  41
  42    model_u(:,:) = 0.0
  43    model_v(:,:) = 0.0
  44    model_t(:,:) = 0.0
  45    model_q(:,:) = 0.0
  46
  47    if ( it > 1 ) then
  48       do n=iv%info(mtgirs)%n1, iv%info(mtgirs)%n2
  49          do k=1, iv%info(mtgirs)%levels(n)
  50             if (iv%mtgirs(n)%u(k)%qc == fails_error_max) iv%mtgirs(n)%u(k)%qc = 0
  51             if (iv%mtgirs(n)%v(k)%qc == fails_error_max) iv%mtgirs(n)%v(k)%qc = 0
  52             if (iv%mtgirs(n)%t(k)%qc == fails_error_max) iv%mtgirs(n)%t(k)%qc = 0
  53             if (iv%mtgirs(n)%q(k)%qc == fails_error_max) iv%mtgirs(n)%q(k)%qc = 0
  54          end do
  55       end do
  56    end if
  57
  58    do n=iv%info(mtgirs)%n1, iv%info(mtgirs)%n2
  59       if (iv%info(mtgirs)%levels(n) < 1) cycle
  60
  61       ! [1.1] Get horizontal interpolation weights:
  62
  63       if (position_lev_dependant) then
  64          i(:)   = iv%info(mtgirs)%i(:,n)
  65          j(:)   = iv%info(mtgirs)%j(:,n)
  66          dx(:)  = iv%info(mtgirs)%dx(:,n)
  67          dy(:)  = iv%info(mtgirs)%dy(:,n)
  68          dxm(:) = iv%info(mtgirs)%dxm(:,n)
  69          dym(:) = iv%info(mtgirs)%dym(:,n)
  70          do k=kts,kte
  71             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)) &
  72                + 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))
  73             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)) &
  74                + 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))
  75          end do
  76       else
  77          i(1)   = iv%info(mtgirs)%i(1,n)
  78          j(1)   = iv%info(mtgirs)%j(1,n)
  79          dx(1)  = iv%info(mtgirs)%dx(1,n)
  80          dy(1)  = iv%info(mtgirs)%dy(1,n)
  81          dxm(1) = iv%info(mtgirs)%dxm(1,n)
  82          dym(1) = iv%info(mtgirs)%dym(1,n)
  83
  84          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)) &
  85                        + 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))
  86          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)) &
  87                        + 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))
  88       end if
  89
  90       do k=1, iv%info(mtgirs)%levels(n)
  91          if (iv%mtgirs(n)%p(k) > 1.0) then
  92             call da_to_zk (iv%mtgirs(n)%p(k), v_p, v_interp_p, iv%info(mtgirs)%zk(k,n))
  93          else if (iv%mtgirs(n)%h(k) > 0.0) then
  94             call da_to_zk (iv%mtgirs(n)%h(k), v_h, v_interp_h, iv%info(mtgirs)%zk(k,n))
  95          end if
  96       end do
  97
  98    end do
  99
 100    call da_convert_zk (iv%info(mtgirs))
 101
 102    if (.not. anal_type_verify) then
 103       do n=iv%info(mtgirs)%n1,iv%info(mtgirs)%n2
 104          do k=1, iv%info(mtgirs)%levels(n)
 105             if (iv%info(mtgirs)%zk(k,n) < 0.0) then
 106                iv%mtgirs(n)%u(k)%qc = missing_data
 107                iv%mtgirs(n)%v(k)%qc = missing_data
 108                iv%mtgirs(n)%t(k)%qc = missing_data
 109                iv%mtgirs(n)%q(k)%qc = missing_data
 110             end if
 111          end do
 112       end do
 113    end if
 114
 115    ! [1.2] Interpolate horizontally to ob:
 116
 117 #ifdef A2C
 118    call da_interp_lin_3d (grid%xb%u, iv%info(mtgirs), model_u,'u')
 119    call da_interp_lin_3d (grid%xb%v, iv%info(mtgirs), model_v,'v')
 120 #else
 121    call da_interp_lin_3d (grid%xb%u, iv%info(mtgirs), model_u)
 122    call da_interp_lin_3d (grid%xb%v, iv%info(mtgirs), model_v)
 123 #endif
 124    call da_interp_lin_3d (grid%xb%t, iv%info(mtgirs), model_t)
 125    call da_interp_lin_3d (grid%xb%q, iv%info(mtgirs), model_q)
 126
 127    do n=iv%info(mtgirs)%n1, iv%info(mtgirs)%n2
 128       !----------------------------------------------------------------------
 129       ! [2.0] Initialise components of innovation vector:
 130       !----------------------------------------------------------------------
 131
 132       do k = 1, iv%info(mtgirs)%levels(n)
 133          iv%mtgirs(n)%u(k)%inv = 0.0
 134          iv%mtgirs(n)%v(k)%inv = 0.0
 135          iv%mtgirs(n)%t(k)%inv = 0.0
 136          iv%mtgirs(n)%q(k)%inv = 0.0
 137
 138          !-------------------------------------------------------------------
 139          ! [3.0] Interpolation:
 140          !-------------------------------------------------------------------
 141          if (wind_sd_mtgirs) then
 142              call da_ffdduv_model (speed,direction,model_u(k,n), model_v(k,n), convert_uv2fd)
 143
 144               if (ob%mtgirs(n)%u(k) > missing_r .AND. iv%mtgirs(n)%u(k)%qc >= obs_qc_pointer) then
 145                   iv%mtgirs(n)%u(k)%inv = ob%mtgirs(n)%u(k) - speed
 146               end if
 147
 148               if (ob%mtgirs(n)%v(k) > missing_r .AND. iv%mtgirs(n)%v(k)%qc >= obs_qc_pointer) then
 149                   iv%mtgirs(n)%v(k)%inv = ob%mtgirs(n)%v(k) - direction
 150                   if (iv%mtgirs(n)%v(k)%inv > 180.0 ) iv%mtgirs(n)%v(k)%inv = iv%mtgirs(n)%v(k)%inv - 360.0
 151                   if (iv%mtgirs(n)%v(k)%inv < -180.0 ) iv%mtgirs(n)%v(k)%inv = iv%mtgirs (n)%v(k)%inv + 360.0
 152               end if
 153          else
 154               if (ob%mtgirs(n)%u(k) > missing_r .AND. iv%mtgirs(n)%u(k)%qc >= obs_qc_pointer) then
 155                   iv%mtgirs(n)%u(k)%inv = ob%mtgirs(n)%u(k) - model_u(k,n)
 156               end if
 157
 158               if (ob%mtgirs(n)%v(k) > missing_r .AND. iv%mtgirs(n)%v(k)%qc >= obs_qc_pointer) then
 159                   iv%mtgirs(n)%v(k)%inv = ob%mtgirs(n)%v(k) - model_v(k,n)
 160               end if
 161           end if
 162
 163
 164          if (ob%mtgirs(n)%t(k) > missing_r .AND. iv%mtgirs(n)%t(k)%qc >= obs_qc_pointer) then
 165             iv%mtgirs(n)%t(k)%inv = ob%mtgirs(n)%t(k) - model_t(k,n)
 166          end if
 167
 168          if (ob%mtgirs(n)%q(k) > missing_r .AND. iv%mtgirs(n)%q(k)%qc >= obs_qc_pointer) then
 169             iv%mtgirs(n)%q(k)%inv = ob%mtgirs(n)%q(k) - model_q(k,n)
 170          end if
 171       end do
 172    end do
 173
 174    !----------------------------------------------------------------------
 175    ! [5.0] Perform optional maximum error check:
 176    !----------------------------------------------------------------------
 177
 178     if ( check_max_iv ) &
 179        call da_check_max_iv_mtgirs (iv, it, num_qcstat_conv)
 180
 181    deallocate (model_u)
 182    deallocate (model_v)
 183    deallocate (model_t)
 184    deallocate (model_q)
 185
 186    if (trace_use_dull) call da_trace_exit ("da_get_innov_vector_mtgirs")
 187
 188 end subroutine da_get_innov_vector_mtgirs
 189
 190