var/da/da_sound/da_get_innov_vector_sonde_sfc.inc

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