var/da/da_gpspw/da_get_innov_vector_gpspw.inc

   1 subroutine da_get_innov_vector_gpspw (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
  18    integer :: n        ! Loop counter.
  19    integer :: i, j     ! Index dimension.
  20    real    :: dx, dxm  ! Interpolation weights.
  21    real    :: dy, dym  ! Interpolation weights.
  22    real    :: model_tpw! Model value u at oblocation.
  23
  24    !--------------------------------------------------------------------------
  25    integer :: k            ! Index dimension
  26    real    :: dpw, ddpw    ! adjustment pw [mm]
  27    real    :: obs_terr     ! real terrain height at GPS site [m]
  28    real    :: model_terr   ! model terrain height at GPS site[m]
  29    real    :: model_q(kts:kte) ! model q at GPS site [kg/kg]
  30    real    :: model_z(kts:kte) ! model z at GPS site [m]
  31    real    :: model_rho(kts:kte) ! model rho at GPS site [kg/m^3]
  32
  33    if (trace_use_dull) call da_trace_entry("da_get_innov_vector_gpspw")
  34
  35    ! unit_gps = myproc + 140
  36    !---------------------------------------------------------------------------
  37
  38    if (iv%info(gpspw)%nlocal > 0) then
  39
  40
  41       do n=iv%info(gpspw)%n1,iv%info(gpspw)%n2
  42
  43          if( iv % gpspw(n) % tpw % qc == fails_error_max .and. it > 1) &
  44              iv % gpspw(n) % tpw % qc = 0
  45
  46          ! [1.1] Get horizontal interpolation weights:
  47
  48          i   = iv%info(gpspw)%i(1,n)
  49          j   = iv%info(gpspw)%j(1,n)
  50          dx  = iv%info(gpspw)%dx(1,n)
  51          dy  = iv%info(gpspw)%dy(1,n)
  52          dxm = iv%info(gpspw)%dxm(1,n)
  53          dym = iv%info(gpspw)%dym(1,n)
  54
  55          model_tpw = dym*(dxm*grid%xb%tpw(i,j)   + dx*grid%xb%tpw(i+1,j)) + &
  56                      dy *(dxm*grid%xb%tpw(i,j+1) + dx*grid%xb%tpw(i+1,j+1))
  57
  58             ! To compute the 'inv':
  59             if ( .not. pseudo_tpw ) iv % gpspw(n) % tpw % inv = 0.0
  60             if (ob % gpspw(n) % tpw > missing_r .AND. &
  61                 iv % gpspw(n) % tpw % qc >= obs_qc_pointer) then
  62                dpw = 0.0
  63                obs_terr   = iv%info(gpspw)%elv(n)
  64                model_terr = dym*(dxm*grid%xb%terr(i,j)   + dx*grid%xb%terr(i+1,j)) + &
  65                           dy *(dxm*grid%xb%terr(i,j+1) + dx*grid%xb%terr(i+1,j+1))
  66                if (obs_terr <= model_terr) then
  67                   model_q(1) = dym*(dxm*grid%xb%q(i,j,1)   + dx*grid%xb%q(i+1,j,1)) + &
  68                              dy *(dxm*grid%xb%q(i,j+1,1) + dx*grid%xb%q(i+1,j+1,1))
  69                   model_rho(1) = dym*(dxm*grid%xb%rho(i,j,1) + dx*grid%xb%rho(i+1,j,1)) + &
  70                              dy *(dxm*grid%xb%rho(i,j+1,1) + dx*grid%xb%rho(i+1,j+1,1))
  71                   dpw = model_rho(1) * model_q(1) *( obs_terr - model_terr)
  72                else
  73                   model_z(1) = dym*(dxm*grid%xb%hf(i,j,1)   + dx*grid%xb%hf(i+1,j,1)) + &
  74                             dy *(dxm*grid%xb%hf(i,j+1,1) + dx*grid%xb%hf(i+1,j+1,1))
  75                   do k = kts, kte
  76                      if (model_z(k) >= obs_terr) exit
  77
  78                      model_z(k+1) = dym*(dxm*grid%xb%hf(i,j,k+1)   + dx*grid%xb%hf(i+1,j,k+1)) + &
  79                                 dy *(dxm*grid%xb%hf(i,j+1,k+1) + dx*grid%xb%hf(i+1,j+1,k+1))
  80                      model_q(k) = dym*(dxm*grid%xb%q(i,j,k)   + dx*grid%xb%q(i+1,j,k)) + &
  81                               dy *(dxm*grid%xb%q(i,j+1,k) + dx*grid%xb%q(i+1,j+1,k))
  82                      model_rho(k) = dym*(dxm*grid%xb%rho(i,j,k)   + dx*grid%xb%rho(i+1,j,k)) + &
  83                                 dy *(dxm*grid%xb%rho(i,j+1,k) + dx*grid%xb%rho(i+1,j+1,k))
  84
  85                      if (model_z(k+1) <= obs_terr) then
  86                         ddpw = model_rho(k) * model_q(k) *( model_z(k+1) - model_z(k))
  87                      else
  88                         ddpw = model_rho(k) * model_q(k) *( obs_terr - model_z(k))
  89                      end if
  90
  91                      dpw = dpw + ddpw
  92                   end do
  93                end if
  94                if ( pseudo_tpw .and. it == 1 ) then
  95                ! To compute the 'ob':
  96                  ob % gpspw(n) % tpw = iv % gpspw(n) % tpw % inv + model_tpw - 0.1*dpw
  97                else
  98                  iv % gpspw(n) % tpw % inv = ob % gpspw(n) % tpw - model_tpw + 0.1*dpw
  99                end if
 100             end if
 101
 102       end do
 103
 104       !------------------------------------------------------------------------
 105       ! [5.0] Perform optional maximum error check:
 106       !------------------------------------------------------------------------
 107       if ( .not. pseudo_tpw .and. check_max_iv ) &
 108          call da_check_max_iv_gpspw(iv, it, num_qcstat_conv)
 109    end if
 110    if (trace_use_dull) call da_trace_exit("da_get_innov_vector_gpspw")
 111
 112 end subroutine da_get_innov_vector_gpspw
 113
 114