var/da/da_metar/da_get_innov_vector_metar.inc

   1 subroutine da_get_innov_vector_metar(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    real  :: speed, direction
  17    integer :: n        ! Loop counter.
  18    integer :: i, j, k  ! Index dimension.
  19    real    :: dx, dxm  ! Interpolation weights.
  20    real    :: dy, dym  ! Interpolation weights.
  21    real, allocatable :: model_u(:,:)  ! Model value u at oblocation.
  22    real, allocatable :: model_v(:,:)  ! Model value v at oblocation.
  23    real, allocatable :: model_t(:,:)  ! Model value t at oblocation.
  24    real, allocatable :: model_p(:,:)  ! Model value p at oblocation.
  25    real, allocatable :: model_q(:,:)  ! Model value q at oblocation.
  26    real, allocatable :: model_hsm(:,:)
  27
  28    real    :: v_h(kms:kme)      ! Model value h at ob hor. location.
  29    real    :: v_p(kms:kme)      ! Model value p at ob hor. location.
  30
  31    real    :: hd, psfcm
  32
  33    real    :: ho, to, qo
  34
  35    if (trace_use_dull) call da_trace_entry("da_get_innov_vector_metar")
  36
  37    allocate (model_u  (1,iv%info(metar)%n1:iv%info(metar)%n2))
  38    allocate (model_v  (1,iv%info(metar)%n1:iv%info(metar)%n2))
  39    allocate (model_t  (1,iv%info(metar)%n1:iv%info(metar)%n2))
  40    allocate (model_q  (1,iv%info(metar)%n1:iv%info(metar)%n2))
  41    allocate (model_p  (1,iv%info(metar)%n1:iv%info(metar)%n2))
  42    allocate (model_hsm(1,iv%info(metar)%n1:iv%info(metar)%n2))
  43
  44    if ( it > 1 ) then
  45       do n=iv%info(metar)%n1,iv%info(metar)%n2
  46          if (iv%metar(n)%u%qc == fails_error_max) iv%metar(n)%u%qc = 0
  47          if (iv%metar(n)%v%qc == fails_error_max) iv%metar(n)%v%qc = 0
  48          if (iv%metar(n)%t%qc == fails_error_max) iv%metar(n)%t%qc = 0
  49          if (iv%metar(n)%p%qc == fails_error_max) iv%metar(n)%p%qc = 0
  50          if (iv%metar(n)%q%qc == fails_error_max) iv%metar(n)%q%qc = 0
  51       end do
  52    end if
  53
  54    if (sfc_assi_options == sfc_assi_options_1) then
  55
  56       do n=iv%info(metar)%n1,iv%info(metar)%n2
  57          ! [1.1] Get horizontal interpolation weights:
  58
  59          i   = iv%info(metar)%i(1,n)
  60          j   = iv%info(metar)%j(1,n)
  61          dx  = iv%info(metar)%dx(1,n)
  62          dy  = iv%info(metar)%dy(1,n)
  63          dxm = iv%info(metar)%dxm(1,n)
  64          dym = iv%info(metar)%dym(1,n)
  65
  66          ! Surface correction
  67
  68          iv%metar(n)%p%inv = ob%metar(n)%p
  69          iv%metar(n)%t%inv = ob%metar(n)%t
  70          iv%metar(n)%q%inv = ob%metar(n)%q
  71          iv%metar(n)%u%inv = ob%metar(n)%u
  72          iv%metar(n)%v%inv = ob%metar(n)%v
  73
  74          if (iv % metar(n) % h > missing_r) then
  75             do k=kts,kte
  76                v_h(k) = dym*(dxm*grid%xb%h(i,j  ,k) + dx*grid%xb%h(i+1,j  ,k)) &
  77                   + dy *(dxm*grid%xb%h(i,j+1,k) + dx*grid%xb%h(i+1,j+1,k))
  78             end do
  79
  80             hd = v_h(kts) - iv % metar(n) % h
  81
  82             if (abs(hd) <= Max_StHeight_Diff .or. anal_type_verify) then
  83                if (iv % metar(n) % h < v_h(kts)) then
  84                   iv%info(metar)%zk(:,n) = 1.0+1.0e-6
  85                   call da_obs_sfc_correction(iv%info(metar), iv%metar(n), n, grid%xb)
  86                else
  87                   call da_to_zk(iv % metar(n) % h, v_h, v_interp_h, iv%info(metar)%zk(1,n))
  88                end if
  89             end if
  90          else if (ob % metar(n) % p > 1.0) then
  91             do k=kts,kte
  92                v_p(k) = dym*(dxm*grid%xb%p(i,j  ,k) + dx*grid%xb%p(i+1,j  ,k)) &
  93                   + dy *(dxm*grid%xb%p(i,j+1,k) + dx*grid%xb%p(i+1,j+1,k))
  94             end do
  95
  96             call da_to_zk(ob % metar(n) % p, v_p, v_interp_p, iv%info(metar)%zk(1,n))
  97
  98             if (iv%info(metar)%zk(1,n) < 0.0 .and.  .not.anal_type_verify) then
  99                iv % metar(n) % p % inv = missing_r
 100                iv % metar(n) % p % qc  = missing_data
 101                iv%info(metar)%zk(:,n) = 1.0+1.0e-6
 102             end if
 103          end if
 104       end do
 105
 106       call da_convert_zk(iv%info(metar))
 107
 108       !------------------------------------------------------------------------
 109       ! [2.0] Initialise components of innovation vector:
 110       !------------------------------------------------------------------------
 111
 112       if (.not.anal_type_verify ) then
 113          do n=iv%info(metar)%n1,iv%info(metar)%n2
 114             if (iv%info(metar)%zk(1,n) < 0.0) then
 115                iv % metar(n) % u % qc = missing_data
 116                iv % metar(n) % v % qc = missing_data
 117                iv % metar(n) % t % qc = missing_data
 118                iv % metar(n) % q % qc = missing_data
 119                iv % metar(n) % p % qc = missing_data
 120             end if
 121          end do
 122       end if
 123
 124       ! [1.2] Interpolate horizontally:
 125 #ifdef A2C
 126       call da_interp_lin_3d (grid%xb%u, iv%info(metar),model_u,'u')
 127       call da_interp_lin_3d (grid%xb%v, iv%info(metar),model_v,'v')
 128 #else
 129       call da_interp_lin_3d (grid%xb%u, iv%info(metar),model_u)
 130       call da_interp_lin_3d (grid%xb%v, iv%info(metar),model_v)
 131 #endif
 132       call da_interp_lin_3d (grid%xb%t, iv%info(metar),model_t)
 133       call da_interp_lin_3d (grid%xb%q, iv%info(metar),model_q)
 134       call da_interp_lin_3d (grid%xb%p, iv%info(metar),model_p)
 135
 136    else if (sfc_assi_options == sfc_assi_options_2) then
 137       ! 1.2.1 Surface assimilation approach 2
 138       !(10-m u, v, 2-m t, q, and sfc_p)
 139
 140       call da_interp_lin_2d (grid%xb%u10,  iv%info(metar), 1, model_u)
 141       call da_interp_lin_2d (grid%xb%v10,  iv%info(metar), 1, model_v)
 142       call da_interp_lin_2d (grid%xb%t2,   iv%info(metar), 1, model_t)
 143       call da_interp_lin_2d (grid%xb%q2,   iv%info(metar), 1, model_q)
 144       call da_interp_lin_2d (grid%xb%psfc, iv%info(metar), 1, model_p)
 145       call da_interp_lin_2d (grid%xb%terr, iv%info(metar), 1, model_hsm)
 146
 147       do n=iv%info(metar)%n1,iv%info(metar)%n2
 148
 149          iv%metar(n)%p%inv = ob%metar(n)%p
 150          iv%metar(n)%t%inv = ob%metar(n)%t
 151          iv%metar(n)%q%inv = ob%metar(n)%q
 152          iv%metar(n)%u%inv = ob%metar(n)%u
 153          iv%metar(n)%v%inv = ob%metar(n)%v
 154
 155          if (iv%metar(n)%p%qc >= 0) then
 156             ! model surface p, t, q, h at observed site:
 157
 158             ho = iv%metar(n)%h
 159             to = -888888.0
 160             qo = -888888.0
 161
 162             if (iv%metar(n)%t%qc >= 0 .and. iv%metar(n)%q%qc >= 0) then
 163                to = ob%metar(n)%t
 164                qo = ob%metar(n)%q
 165                call da_sfc_pre(psfcm, model_p(1,n), model_t(1,n), model_q(1,n), model_hsm(1,n), ho, to, qo)
 166             else if (iv%metar(n)%t%qc >= 0 .and. iv%metar(n)%q%qc < 0) then
 167                to = ob%metar(n)%t
 168                call da_sfc_pre(psfcm, model_p(1,n), model_t(1,n), model_q(1,n), model_hsm(1,n), ho, to)
 169             else
 170                call da_sfc_pre(psfcm, model_p(1,n), model_t(1,n), model_q(1,n), model_hsm(1,n), ho)
 171             end if
 172
 173             ! Pressure at the observed height:
 174             model_p(1,n) = psfcm
 175          end if
 176       end do
 177    end if
 178
 179    do n=iv%info(metar)%n1,iv%info(metar)%n2
 180
 181       !-----------------------------------------------------------------------
 182       ! [3.0] Fast interpolation:
 183       !-----------------------------------------------------------------------
 184       if (wind_sd_metar) then
 185          call da_ffdduv_model (speed,direction,model_u(1,n), model_v(1,n), convert_uv2fd)
 186
 187          if (ob%metar(n)%u > missing_r .AND. iv%metar(n)%u%qc >= obs_qc_pointer) then
 188              iv%metar(n)%u%inv = iv%metar(n)%u%inv - speed
 189          else
 190              iv % metar(n) % u % inv = 0.0
 191          end if
 192
 193          if (ob%metar(n)%v > missing_r .AND. iv%metar(n)%v%qc >= obs_qc_pointer) then
 194              iv%metar(n)%v%inv = iv%metar(n)%v%inv - direction
 195              if (iv%metar(n)%v%inv > 180.0 ) iv%metar(n)%v%inv = iv%metar(n)%v%inv - 360.0
 196              if (iv%metar(n)%v%inv < -180.0 ) iv%metar(n)%v%inv = iv%metar(n)%v%inv + 360.0
 197          else
 198              iv % metar(n) % v % inv = 0.0
 199          end if
 200       else
 201          if (ob % metar(n) % u > missing_r .AND. iv % metar(n) % u % qc >= obs_qc_pointer) then
 202              iv % metar(n) % u % inv = iv%metar(n)%u%inv - model_u(1,n)
 203          else
 204              iv % metar(n) % u % inv = 0.0
 205          end if
 206
 207          if (ob % metar(n) % v > missing_r .AND. iv % metar(n) % v % qc >= obs_qc_pointer) then
 208              iv % metar(n) % v % inv = iv%metar(n)%v%inv - model_v(1,n)
 209          else
 210              iv % metar(n) % v % inv = 0.0
 211          end if
 212       end if
 213
 214       !if (ob % metar(n) % p > 0.0 .AND. iv % metar(n) % p % qc >= obs_qc_pointer) then
 215       if ( iv % metar(n) % p % qc >= obs_qc_pointer ) then
 216          iv % metar(n) % p % inv = iv%metar(n)%p%inv - model_p(1,n)
 217       else
 218          iv % metar(n) % p % inv = 0.0
 219       end if
 220
 221       if (ob % metar(n) % t > 0.0 .AND. &
 222          iv % metar(n) % t % qc >= obs_qc_pointer) then
 223          iv % metar(n) % t % inv = iv%metar(n)%t%inv - model_t(1,n)
 224       else
 225          iv % metar(n) % t % inv = 0.0
 226       end if
 227
 228       if (ob % metar(n) % q > 0.0 .AND. &
 229          iv % metar(n) % q % qc >= obs_qc_pointer) then
 230          iv % metar(n) % q % inv = iv%metar(n)%q%inv - model_q(1,n)
 231       else
 232          iv % metar(n) % q % inv = 0.0
 233       end if
 234    end do
 235
 236    !-----------------------------------------------------------------------
 237    ! [5.0] Perform optional maximum error check:
 238    !-----------------------------------------------------------------------
 239
 240    if ( check_max_iv ) &
 241       call da_check_max_iv_metar(iv,ob, it, num_qcstat_conv)
 242
 243    deallocate (model_u)
 244    deallocate (model_v)
 245    deallocate (model_t)
 246    deallocate (model_q)
 247    deallocate (model_p)
 248    deallocate (model_hsm)
 249
 250    if (trace_use_dull) call da_trace_exit("da_get_innov_vector_metar")
 251
 252 end subroutine da_get_innov_vector_metar
 253
 254