var/da/da_synop/da_get_innov_vector_synop.inc

   1 subroutine da_get_innov_vector_synop( 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    ! model 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
  23    real    :: v_h(kms:kme)      ! Model value h at ob hor. loc
  24    real    :: v_p(kms:kme)      ! Model value p at ob hor. loc
  25
  26    real    :: hd, psfcm
  27
  28    real    :: ho, to, qo
  29    real    :: speed, direction
  30    real    :: rh_error, es, qs, rho
  31    real    :: xcorr
  32    real    :: err_inflate_fac
  33
  34    real, allocatable :: model_u(:,:)
  35    real, allocatable :: model_v(:,:)
  36    real, allocatable :: model_t(:,:)
  37    real, allocatable :: model_q(:,:)
  38    real, allocatable :: model_p(:,:)
  39    real, allocatable :: model_psfc(:,:)
  40    real, allocatable :: model_hsm(:,:)
  41    real, allocatable :: model_qs(:)
  42
  43    if (trace_use_dull) call da_trace_entry("da_get_innov_vector_synop")
  44
  45    allocate (model_u  (1,iv%info(synop)%n1:iv%info(synop)%n2))
  46    allocate (model_v  (1,iv%info(synop)%n1:iv%info(synop)%n2))
  47    allocate (model_t  (1,iv%info(synop)%n1:iv%info(synop)%n2))
  48    allocate (model_q  (1,iv%info(synop)%n1:iv%info(synop)%n2))
  49    allocate (model_p  (1,iv%info(synop)%n1:iv%info(synop)%n2))
  50    allocate (model_psfc(1,iv%info(synop)%n1:iv%info(synop)%n2))
  51    allocate (model_hsm(1,iv%info(synop)%n1:iv%info(synop)%n2))
  52
  53    if ( it > 1 ) then
  54       do n=iv%info(synop)%n1,iv%info(synop)%n2
  55          if (iv%synop(n)%u%qc == fails_error_max) iv%synop(n)%u%qc = 0
  56          if (iv%synop(n)%v%qc == fails_error_max) iv%synop(n)%v%qc = 0
  57          if (iv%synop(n)%t%qc == fails_error_max) iv%synop(n)%t%qc = 0
  58          if (iv%synop(n)%p%qc == fails_error_max) iv%synop(n)%p%qc = 0
  59          if (iv%synop(n)%q%qc == fails_error_max) iv%synop(n)%q%qc = 0
  60       end do
  61    end if
  62
  63    if (sfc_assi_options == sfc_assi_options_1) then
  64
  65       do n=iv%info(synop)%n1,iv%info(synop)%n2
  66
  67          if ( it == 1 .and. sfc_hori_intp_options == 2 ) then
  68             ! choose the nearest model point with matching land/ocean type and
  69             ! smallest height difference with the ob
  70             ! dx, dy, dxm, dym are reset in da_sfc_hori_interp_weights
  71             call da_sfc_hori_interp_weights(n, iv%info(synop), iv%synop(n), grid%xb, 1)
  72          end if
  73
  74          ! [1.1] Get horizontal interpolation weights:
  75          i   = iv%info(synop)%i(1,n)
  76          j   = iv%info(synop)%j(1,n)
  77          dx  = iv%info(synop)%dx(1,n)
  78          dy  = iv%info(synop)%dy(1,n)
  79          dxm = iv%info(synop)%dxm(1,n)
  80          dym = iv%info(synop)%dym(1,n)
  81
  82          ! Surface correction
  83
  84          iv%synop(n)%p%inv = ob%synop(n)%p
  85          iv%synop(n)%t%inv = ob%synop(n)%t
  86          iv%synop(n)%q%inv = ob%synop(n)%q
  87          iv%synop(n)%u%inv = ob%synop(n)%u
  88          iv%synop(n)%v%inv = ob%synop(n)%v
  89
  90          if ( it == 1 .and. sfc_hori_intp_options == 2 ) then
  91             ! some obs might have been rejected in da_sfc_hori_interp_weights
  92             ! due to mismatched land/ocean type
  93             if ( iv%synop(n)%p%qc < 0 ) iv%synop(n)%p%inv = missing_r
  94             if ( iv%synop(n)%t%qc < 0 ) iv%synop(n)%t%inv = missing_r
  95             if ( iv%synop(n)%q%qc < 0 ) iv%synop(n)%q%inv = missing_r
  96             if ( iv%synop(n)%u%qc < 0 ) iv%synop(n)%u%inv = missing_r
  97             if ( iv%synop(n)%v%qc < 0 ) iv%synop(n)%v%inv = missing_r
  98          end if
  99
 100          if (iv % synop(n) % h > missing_r) then
 101             do k=kts,kte
 102                v_h(k) = dym*(dxm*grid%xb%h(i,j  ,k) + dx*grid%xb%h(i+1,j  ,k)) &
 103                        + dy *(dxm*grid%xb%h(i,j+1,k) + dx*grid%xb%h(i+1,j+1,k))
 104             end do
 105
 106             hd = v_h(kts) - iv % synop(n) % h
 107             if (abs(hd) <= Max_StHeight_Diff .or.  anal_type_verify ) then
 108                if (iv % synop(n) % h < v_h(kts)) then
 109                   iv%info(synop)%zk(:,n) = 1.0+1.0e-6
 110                   call da_obs_sfc_correction(iv%info(synop), iv%synop(n), n, grid%xb)
 111
 112                   if ( sfcht_adjust_q ) then
 113                      ! calculate RH with original t, p, q
 114                      if ( abs(ob%synop(n)%q-missing_r) > 1.0 ) then
 115                         if ( abs(ob%synop(n)%p-missing_r) > 1.0 .and. &
 116                              abs(ob%synop(n)%t-missing_r) > 1.0 ) then
 117                            call da_tpq_to_rh(ob%synop(n)%t, ob%synop(n)%p, ob%synop(n)%q, es, qs, rho)
 118                         end if
 119                      end if
 120                      ! calculate corrected q with corrected t, p and original RH
 121                      if ( abs(ob%synop(n)%q-missing_r) > 1.0 ) then
 122                         if ( abs(iv%synop(n)%p%inv-missing_r) > 1.0 .and. &
 123                              abs(iv%synop(n)%t%inv-missing_r) > 1.0 ) then
 124                            call da_tp_to_qs(iv%synop(n)%t%inv, iv%synop(n)%p%inv, es, qs)
 125                            iv%synop(n)%q%inv = qs * rho * 0.01
 126                            iv%synop(n)%q%qc  = surface_correction
 127                         end if
 128                     end if
 129                   end if !sfcht_adjust_q
 130
 131                else
 132                   call da_to_zk(iv % synop(n) % h, v_h, v_interp_h, iv%info(synop)%zk(1,n))
 133                end if
 134             end if
 135
 136             if ( (it == 1) .and. (obs_err_inflate) .and. (abs(hd) <= Max_StHeight_Diff) ) then
 137                err_inflate_fac = MIN(exp(abs(hd)/stn_ht_diff_scale),100.)
 138                !iv%synop(n)%u%error = iv%synop(n)%u%error * err_inflate_fac
 139                !iv%synop(n)%v%error = iv%synop(n)%v%error * err_inflate_fac
 140                iv%synop(n)%t%error = iv%synop(n)%t%error * err_inflate_fac
 141                iv%synop(n)%p%error = iv%synop(n)%p%error * err_inflate_fac
 142                iv%synop(n)%q%error = iv%synop(n)%q%error * err_inflate_fac
 143             end if
 144
 145          else if (ob % synop(n) % p > 1.0) then
 146             do k=kts,kte
 147               v_p(k) = dym*(dxm*grid%xb%p(i,j  ,k) + dx*grid%xb%p(i+1,j  ,k)) &
 148                        + dy *(dxm*grid%xb%p(i,j+1,k) + dx*grid%xb%p(i+1,j+1,k))
 149             end do
 150
 151             call da_to_zk(ob % synop(n) % p, v_p, v_interp_p, iv%info(synop)%zk(1,n))
 152             if (iv%info(synop)%zk(1,n) < 0.0 .and. .not.anal_type_verify ) then
 153                iv % synop(n) % p % inv = missing_r
 154                iv % synop(n) % p % qc  = missing_data
 155                iv%info(synop)%zk(:,n) = 1.0+1.0e-6
 156             end if
 157          end if
 158       end do
 159
 160       call da_convert_zk (iv%info(synop))
 161
 162       if (.not.anal_type_verify ) then
 163          do n=iv%info(synop)%n1,iv%info(synop)%n2
 164             if (iv%info(synop)%zk(1,n) < 0.0) then
 165                iv % synop(n) % u % qc = missing_data
 166                iv % synop(n) % v % qc = missing_data
 167                iv % synop(n) % t % qc = missing_data
 168                iv % synop(n) % q % qc = missing_data
 169                iv % synop(n) % p % qc = missing_data
 170             end if
 171          end do
 172       end if
 173
 174       ! [1.2] Interpolate horizontally:
 175 #ifdef A2C
 176       call da_interp_lin_3d (grid%xb%u, iv%info(synop), model_u,'u')
 177       call da_interp_lin_3d (grid%xb%v, iv%info(synop), model_v,'v')
 178 #else
 179       call da_interp_lin_3d (grid%xb%u, iv%info(synop), model_u)
 180       call da_interp_lin_3d (grid%xb%v, iv%info(synop), model_v)
 181 #endif
 182       call da_interp_lin_3d (grid%xb%t, iv%info(synop), model_t)
 183       call da_interp_lin_3d (grid%xb%q, iv%info(synop), model_q)
 184       call da_interp_lin_3d (grid%xb%p, iv%info(synop), model_p)
 185       call da_interp_lin_2d (grid%xb%psfc, iv%info(synop), 1, model_psfc(1,:))
 186       call da_interp_lin_2d (grid%xb%terr, iv%info(synop), 1, model_hsm(1,:))
 187
 188       if ( it == 1 .and. q_error_options == 2 ) then
 189          allocate (model_qs(iv%info(synop)%n1:iv%info(synop)%n2))
 190          do n=iv%info(synop)%n1,iv%info(synop)%n2
 191            if ( abs(ob%synop(n)%q-missing_r) > 1.0 ) then
 192                call da_tp_to_qs(model_t(1,n), model_p(1,n), es, model_qs(n))
 193                rh_error = iv%synop(n)%q%error !q error is rh at this stage
 194                iv%synop(n)%q%error = model_qs(n)*rh_error*0.01
 195            else
 196                iv%synop(n)%q%error = missing_r
 197                iv%synop(n)%q%qc    = missing_data
 198            end if
 199          end do
 200          deallocate (model_qs)
 201       end if
 202
 203    else if (sfc_assi_options == sfc_assi_options_2) then
 204       ! Surface data assimilation approach 2
 205       !------------------------------------
 206
 207       ! 1.2.1 Surface assimilation approach 2(10-m u, v, 2-m t, q, and sfc_p)
 208
 209       call da_interp_lin_2d (grid%xb%u10,  iv%info(synop), 1, model_u(1,:))
 210       call da_interp_lin_2d (grid%xb%v10,  iv%info(synop), 1, model_v(1,:))
 211       call da_interp_lin_2d (grid%xb%t2,   iv%info(synop), 1, model_t(1,:))
 212       call da_interp_lin_2d (grid%xb%q2,   iv%info(synop), 1, model_q(1,:))
 213       call da_interp_lin_2d (grid%xb%psfc, iv%info(synop), 1, model_p(1,:))
 214       call da_interp_lin_2d (grid%xb%terr, iv%info(synop), 1, model_hsm(1,:))
 215
 216       do n=iv%info(synop)%n1,iv%info(synop)%n2
 217
 218          iv%synop(n)%p%inv = ob%synop(n)%p
 219          iv%synop(n)%t%inv = ob%synop(n)%t
 220          iv%synop(n)%q%inv = ob%synop(n)%q
 221          iv%synop(n)%u%inv = ob%synop(n)%u
 222          iv%synop(n)%v%inv = ob%synop(n)%v
 223
 224          if (iv%synop(n)%p%qc >= 0) then
 225             ho = iv%synop(n)%h
 226             to = -888888.0
 227             qo = -888888.0
 228
 229             if (iv%synop(n)%t%qc >= 0 .and. iv%synop(n)%q%qc >= 0) then
 230                to = ob%synop(n)%t
 231                qo = ob%synop(n)%q
 232                call da_sfc_pre (psfcm, model_p(1,n), model_t(1,n), model_q(1,n), model_hsm(1,n), ho, to, qo)
 233             else if (iv%synop(n)%t%qc >= 0 .and. iv%synop(n)%q%qc < 0) then
 234                to = ob%synop(n)%t
 235                call da_sfc_pre (psfcm, model_p(1,n), model_t(1,n), model_q(1,n), model_hsm(1,n), ho, to)
 236             else
 237                call da_sfc_pre (psfcm, model_p(1,n), model_t(1,n), model_q(1,n), model_hsm(1,n), ho)
 238             end if
 239
 240             ! Pressure at the observed height:
 241             model_p(1,n) = psfcm
 242          end if
 243       end do
 244    end if
 245
 246    do n=iv%info(synop)%n1,iv%info(synop)%n2
 247
 248       !--------------------------------------------------------------------
 249       !     [3.0] Fast interpolation:
 250       !--------------------------------------------------------------------
 251       if(wind_sd_synop)then
 252          call da_ffdduv_model (speed,direction,model_u(1,n), model_v(1,n), convert_uv2fd)
 253
 254          if (ob%synop(n)%u > missing_r .AND. iv%synop(n)%u%qc >= obs_qc_pointer) then
 255              iv%synop(n)%u%inv = iv%synop(n)%u%inv - speed
 256          else
 257              iv % synop(n) % u % inv = 0.0
 258          end if
 259
 260          if (ob%synop(n)%v > missing_r .AND. iv%synop(n)%v%qc >= obs_qc_pointer) then
 261              iv%synop(n)%v%inv = iv%synop(n)%v%inv - direction
 262              if (iv%synop(n)%v%inv > 180.0 ) iv%synop(n)%v%inv = iv%synop(n)%v%inv - 360.0
 263              if (iv%synop(n)%v%inv < -180.0 ) iv%synop(n)%v%inv = iv%synop(n)%v%inv + 360.0
 264          else
 265              iv % synop(n) % v % inv = 0.0
 266          end if
 267       else
 268          if (ob % synop(n) % u > missing_r .AND. iv % synop(n) % u % qc >= obs_qc_pointer) then
 269              iv % synop(n) % u % inv = iv%synop(n)%u%inv - model_u(1,n)
 270          else
 271              iv % synop(n) % u % inv = 0.0
 272          end if
 273
 274          if (ob % synop(n) % v > missing_r .AND. iv % synop(n) % v % qc >= obs_qc_pointer) then
 275              iv % synop(n) % v % inv = iv%synop(n)%v%inv - model_v(1,n)
 276          else
 277              iv % synop(n) % v % inv = 0.0
 278          end if
 279       end if
 280
 281       !if (ob % synop(n) % p > 0.0 .AND. iv % synop(n) % p % qc >= obs_qc_pointer) then
 282       if ( iv % synop(n) % p % qc >= obs_qc_pointer ) then
 283          iv % synop(n) % p % inv = iv%synop(n)%p%inv - model_p(1,n)
 284       else
 285          iv % synop(n) % p % inv = 0.0
 286       end if
 287
 288       if (ob % synop(n) % t > 0.0 .AND. iv % synop(n) % t % qc >= obs_qc_pointer) then
 289          iv % synop(n) % t % inv = iv%synop(n)%t%inv - model_t(1,n)
 290       else
 291          iv % synop(n) % t % inv = 0.0
 292       end if
 293
 294       if (ob % synop(n) % q > 0.0 .AND. iv % synop(n) % q % qc >= obs_qc_pointer) then
 295          if ( sfcht_adjust_q ) then
 296             iv % synop(n) % q % inv = iv%synop(n)%q%inv - model_q(1,n)
 297          else
 298             iv % synop(n) % q % inv = ob%synop(n)%q - model_q(1,n)
 299          end if
 300       else
 301          iv % synop(n) % q % inv = 0.0
 302       end if
 303    end do
 304
 305    !--------------------------------------------------------------------
 306    ! [5.0] Perform optional maximum error check:
 307    !--------------------------------------------------------------------
 308
 309    if ( check_max_iv ) &
 310       call da_check_max_iv_synop(iv,ob, it,num_qcstat_conv)
 311
 312    if (check_buddy) call da_check_buddy_synop(iv, ob, grid%dx, it)
 313  !
 314    deallocate (model_u)
 315    deallocate (model_v)
 316    deallocate (model_t)
 317    deallocate (model_q)
 318    deallocate (model_p)
 319    deallocate (model_psfc)
 320    deallocate (model_hsm)
 321
 322    if (trace_use_dull) call da_trace_exit("da_get_innov_vector_synop")
 323
 324 end subroutine da_get_innov_vector_synop
 325
 326