1 subroutine da_get_innov_vector_sound (it,num_qcstat_conv, grid, ob, iv)
3 !-----------------------------------------------------------------------
5 ! Updated for Analysis on Arakawa-C grid
6 ! Author: Syed RH Rizvi, MMM/ESSL/NCAR, Date: 10/22/2008
7 !-----------------------------------------------------------------------
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, k ! Loop counter.
18 integer :: i (kms:kme)
19 integer :: j (kms:kme)
25 real, allocatable :: model_u(:,:) ! Model value u at ob location.
26 real, allocatable :: model_v(:,:) ! Model value v at ob location.
27 real, allocatable :: model_t(:,:) ! Model value t at ob location.
28 real, allocatable :: model_q(:,:) ! Model value q at ob location.
30 real :: v_h(kts:kte) ! Model value h at ob hor. location.
31 real :: v_p(kts:kte) ! Model value p at ob hor. location.
33 if (trace_use_dull) call da_trace_entry ("da_get_innov_vector_sound")
35 allocate (model_u(iv%info(sound)%max_lev,iv%info(sound)%n1:iv%info(sound)%n2))
36 allocate (model_v(iv%info(sound)%max_lev,iv%info(sound)%n1:iv%info(sound)%n2))
37 allocate (model_t(iv%info(sound)%max_lev,iv%info(sound)%n1:iv%info(sound)%n2))
38 allocate (model_q(iv%info(sound)%max_lev,iv%info(sound)%n1:iv%info(sound)%n2))
46 do n=iv%info(sound)%n1,iv%info(sound)%n2
47 do k=1, iv%info(sound)%levels(n)
48 if (iv%sound(n)%u(k)%qc == fails_error_max) iv%sound(n)%u(k)%qc = 0
49 if (iv%sound(n)%v(k)%qc == fails_error_max) iv%sound(n)%v(k)%qc = 0
50 if (iv%sound(n)%t(k)%qc == fails_error_max) iv%sound(n)%t(k)%qc = 0
51 if (iv%sound(n)%q(k)%qc == fails_error_max) iv%sound(n)%q(k)%qc = 0
56 do n=iv%info(sound)%n1, iv%info(sound)%n2
57 if (iv%info(sound)%levels(n) < 1) cycle
59 ! [1.1] Get horizontal interpolation weights:
61 if (position_lev_dependant) then
62 i(:) = iv%info(sound)%i(:,n)
63 j(:) = iv%info(sound)%j(:,n)
64 dx(:) = iv%info(sound)%dx(:,n)
65 dy(:) = iv%info(sound)%dy(:,n)
66 dxm(:) = iv%info(sound)%dxm(:,n)
67 dym(:) = iv%info(sound)%dym(:,n)
69 v_h(k) = dym(k)*(dxm(k)*grid%xb%h(i(k),j(k),k) + dx(k)*grid%xb%h(i(k)+1,j(k),k)) &
70 + dy(k) *(dxm(k)*grid%xb%h(i(k),j(k)+1,k) + dx(k)*grid%xb%h(i(k)+1,j(k)+1,k))
71 v_p(k) = dym(k)*(dxm(k)*grid%xb%p(i(k),j(k),k) + dx(k)*grid%xb%p(i(k)+1,j(k),k)) &
72 + dy(k) *(dxm(k)*grid%xb%p(i(k),j(k)+1,k) + dx(k)*grid%xb%p(i(k)+1,j(k)+1,k))
75 i(1) = iv%info(sound)%i(1,n)
76 j(1) = iv%info(sound)%j(1,n)
77 dx(1) = iv%info(sound)%dx(1,n)
78 dy(1) = iv%info(sound)%dy(1,n)
79 dxm(1) = iv%info(sound)%dxm(1,n)
80 dym(1) = iv%info(sound)%dym(1,n)
82 v_h(kts:kte) = dym(1) * (dxm(1)*grid%xb%h(i(1),j(1),kts:kte) + dx(1)*grid%xb%h(i(1)+1,j(1),kts:kte)) &
83 + dy(1) * (dxm(1)*grid%xb%h(i(1),j(1)+1,kts:kte) + dx(1)*grid%xb%h(i(1)+1,j(1)+1,kts:kte))
84 v_p(kts:kte) = dym(1) * (dxm(1)*grid%xb%p(i(1),j(1),kts:kte) + dx(1)*grid%xb%p(i(1)+1,j(1),kts:kte)) &
85 + dy(1) * (dxm(1)*grid%xb%p(i(1),j(1)+1,kts:kte) + dx(1)*grid%xb%p(i(1)+1,j(1)+1,kts:kte))
88 do k=1, iv%info(sound)%levels(n)
89 if (iv%sound(n)%p(k) > 1.0) then
90 call da_to_zk (iv%sound(n)%p(k), v_p, v_interp_p, iv%info(sound)%zk(k,n))
91 else if (iv%sound(n)%h(k) > 0.0) then
92 call da_to_zk (iv%sound(n)%h(k), v_h, v_interp_h, iv%info(sound)%zk(k,n))
99 call da_convert_zk (iv%info(sound))
101 if (.not. anal_type_verify) then
102 do n=iv%info(sound)%n1,iv%info(sound)%n2
103 do k=1, iv%info(sound)%levels(n)
104 if (iv%info(sound)%zk(k,n) < 0.0) then
105 iv%sound(n)%u(k)%qc = missing_data
106 iv%sound(n)%v(k)%qc = missing_data
107 iv%sound(n)%t(k)%qc = missing_data
108 iv%sound(n)%q(k)%qc = missing_data
114 ! [1.2] Interpolate horizontally to ob:
117 call da_interp_lin_3d (grid%xb%u, iv%info(sound), model_u,'u')
118 call da_interp_lin_3d (grid%xb%v, iv%info(sound), model_v,'v')
120 call da_interp_lin_3d (grid%xb%u, iv%info(sound), model_u)
121 call da_interp_lin_3d (grid%xb%v, iv%info(sound), model_v)
123 call da_interp_lin_3d (grid%xb%t, iv%info(sound), model_t)
124 call da_interp_lin_3d (grid%xb%q, iv%info(sound), model_q)
126 do n=iv%info(sound)%n1, iv%info(sound)%n2
127 !----------------------------------------------------------------------
128 ! [2.0] Initialise components of innovation vector:
129 !----------------------------------------------------------------------
131 do k = 1, iv%info(sound)%levels(n)
132 iv%sound(n)%u(k)%inv = 0.0
133 iv%sound(n)%v(k)%inv = 0.0
134 iv%sound(n)%t(k)%inv = 0.0
135 iv%sound(n)%q(k)%inv = 0.0
136 !-------------------------------------------------------------------
137 ! [3.0] Interpolation:
138 !-------------------------------------------------------------------
140 if (ob%sound(n)%u(k) > missing_r .AND. iv%sound(n)%u(k)%qc >= obs_qc_pointer) then
141 iv%sound(n)%u(k)%inv = ob%sound(n)%u(k) - model_u(k,n)
144 if (ob%sound(n)%v(k) > missing_r .AND. iv%sound(n)%v(k)%qc >= obs_qc_pointer) then
145 iv%sound(n)%v(k)%inv = ob%sound(n)%v(k) - model_v(k,n)
148 if (ob%sound(n)%t(k) > missing_r .AND. iv%sound(n)%t(k)%qc >= obs_qc_pointer) then
149 iv%sound(n)%t(k)%inv = ob%sound(n)%t(k) - model_t(k,n)
152 if (ob%sound(n)%q(k) > missing_r .AND. iv%sound(n)%q(k)%qc >= obs_qc_pointer) then
153 iv%sound(n)%q(k)%inv = ob%sound(n)%q(k) - model_q(k,n)
158 !----------------------------------------------------------------------
159 ! [5.0] Perform optional maximum error check:
160 !----------------------------------------------------------------------
162 if ( check_max_iv ) &
163 call da_check_max_iv_sound (iv, it, num_qcstat_conv)
165 if (check_buddy) call da_check_buddy_sound(iv, grid%dx, it)
172 if (trace_use_dull) call da_trace_exit ("da_get_innov_vector_sound")
174 end subroutine da_get_innov_vector_sound