wrf svn trunk commit r4103
[wrffire.git] / wrfv2_fire / var / da / da_qscat / da_get_innov_vector_qscat.inc
blob33acf75720aba833ad3ea2c49df54b744fb9ebea
1 subroutine da_get_innov_vector_qscat (it,num_qcstat_conv, grid, ob, iv)
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    !-----------------------------------------------------------------------
9    implicit none
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(:,:,:,:)
17    integer :: n        ! Loop counter.
18    integer :: i, j, k  ! Index dimension.
19    real    :: dx, dxm  ! Interpolation weights.
20    real    :: dy, dym  ! Interpolation weights.
22    real, allocatable :: model_u(:,:)  ! Model value u at ob location.
23    real, allocatable :: model_v(:,:)  ! Model value v at ob location.
25    real    :: v_h(kms:kme)      ! Model value h at ob hor. location.
27    if (trace_use_dull) call da_trace_entry("da_get_innov_vector_qscat")
29    allocate (model_u(iv%info(qscat)%max_lev,iv%info(qscat)%n1:iv%info(qscat)%n2))
30    allocate (model_v(iv%info(qscat)%max_lev,iv%info(qscat)%n1:iv%info(qscat)%n2))
32    if ( it > 1 ) then
33       do n=iv%info(qscat)%n1,iv%info(qscat)%n2
34          if (iv%qscat(n)%u%qc == fails_error_max) iv%qscat(n)%u%qc = 0
35          if (iv%qscat(n)%v%qc == fails_error_max) iv%qscat(n)%v%qc = 0
36       end do
37    end if
39    do n=iv%info(qscat)%n1,iv%info(qscat)%n2
41       ! [1.1] Get horizontal interpolation weights:
43       i   = iv%info(qscat)%i(1,n)
44       j   = iv%info(qscat)%j(1,n)
45       dx  = iv%info(qscat)%dx(1,n)
46       dy  = iv%info(qscat)%dy(1,n)
47       dxm = iv%info(qscat)%dxm(1,n)
48       dym = iv%info(qscat)%dym(1,n)
50       do k=kts,kte
51          v_h(k) = dym*(dxm*grid%xb%h(i,j,k)+dx*grid%xb%h(i+1,j,k)) + dy*(dxm*grid%xb%h(i,j+1,k)+dx*grid%xb%h(i+1,j+1,k))
52       end do
54       if (iv % qscat(n) % h > missing_r) then
55          call da_to_zk(iv % qscat(n) % h, v_h, v_interp_h, iv%info(qscat)%zk(1,n))
56          if (iv%info(qscat)%zk(1,n) < 1.0) then
57             iv%info(qscat)%zk(1,n) = 1.0
58          end if
59       end if
60    end do
62    call da_convert_zk (iv%info(qscat))
64    if (.not. anal_type_verify) then
65       do n=iv%info(qscat)%n1,iv%info(qscat)%n2
66          if (iv%info(qscat)%zk(1,n) < 0.0) then
67             iv%qscat(n)%u%qc = missing_data
68             iv%qscat(n)%v%qc = missing_data
69          end if
70       end do
71    end if
73 #ifdef A2C
74    call da_interp_lin_3d (grid%xb%u, iv%info(qscat), model_u,'u')
75    call da_interp_lin_3d (grid%xb%v, iv%info(qscat), model_v,'v')
76 #else
77    call da_interp_lin_3d (grid%xb%u, iv%info(qscat), model_u)
78    call da_interp_lin_3d (grid%xb%v, iv%info(qscat), model_v)
79 #endif
81    do n=iv%info(qscat)%n1,iv%info(qscat)%n2
83       !------------------------------------------------------------------------
84       ! [2.0] Initialise components of innovation vector:
85       !------------------------------------------------------------------------
87       !------------------------------------------------------------------------
88       ! [3.0] Fast interpolation:
89       !------------------------------------------------------------------------
91       if (ob % qscat(n) % u > missing_r .AND. &
92           iv % qscat(n) % u % qc >= obs_qc_pointer) then
93          iv % qscat(n) % u % inv = ob % qscat(n) % u - model_u(1,n)
94       end if
96       if (ob % qscat(n) % v > missing_r .AND. &
97           iv % qscat(n) % v % qc >= obs_qc_pointer) then
98          iv % qscat(n) % v % inv = ob % qscat(n) % v - model_v(1,n)
99       end if
100    end do
102    !------------------------------------------------------------------------
103    ! [5.0] Perform optional maximum error check:
104    !------------------------------------------------------------------------
106    if ( check_max_iv ) &
107       call da_check_max_iv_qscat(iv, it, num_qcstat_conv)       
109    deallocate (model_u)
110    deallocate (model_v)
112    if (trace_use_dull) call da_trace_exit("da_get_innov_vector_qscat")
114 end subroutine da_get_innov_vector_qscat