da_get_innov_vector_buoy.inc
References to this file elsewhere.
1 subroutine da_get_innov_vector_buoy( it, grid, ob, iv)
2
3 !-----------------------------------------------------------------------
4 ! Purpose: TBD
5 !-----------------------------------------------------------------------
6
7 implicit none
8
9 integer, intent(in) :: it ! External iteration.
10 type(domain), intent(in) :: grid ! first guess state.
11 type(y_type), intent(inout) :: ob ! Observation structure.
12 type(iv_type), intent(inout) :: iv ! O-B structure.
13
14 integer :: n ! Loop counter.
15 integer :: i, j, k ! Index dimension.
16 real :: dx, dxm ! Interpolation weights.
17 real :: dy, dym ! Interpolation weights.
18 real, allocatable :: model_u(:,:) ! Model value u at oblocation.
19 real, allocatable :: model_v(:,:) ! Model value v at oblocation.
20 real, allocatable :: model_t(:,:) ! Model value t at oblocation.
21 real, allocatable :: model_p(:,:) ! Model value p at oblocation.
22 real, allocatable :: model_q(:,:) ! Model value q at oblocation.
23 real, allocatable :: model_hsm(:,:)
24
25 real :: v_h(kms:kme) ! Model value h at ob hor. location.
26 real :: v_p(kms:kme) ! Model value p at ob hor. location.
27
28 real :: hd, psfcm
29
30 real :: ho, to, qo
31
32 integer :: itu,ituf,itvv,itvvf,itp,itpf,itt,ittf,itqv,itqvf
33
34 if (trace_use_dull) call da_trace_entry("da_get_innov_vector_buoy")
35
36 allocate (model_u(1,iv%info(buoy)%n1:iv%info(buoy)%n2))
37 allocate (model_v(1,iv%info(buoy)%n1:iv%info(buoy)%n2))
38 allocate (model_t(1,iv%info(buoy)%n1:iv%info(buoy)%n2))
39 allocate (model_p(1,iv%info(buoy)%n1:iv%info(buoy)%n2))
40 allocate (model_q(1,iv%info(buoy)%n1:iv%info(buoy)%n2))
41 allocate (model_hsm(1,iv%info(buoy)%n1:iv%info(buoy)%n2))
42
43 itu = 0; itvv = 0; itp = 0; itt = 0; itqv = 0;
44 ituf = 0; itvvf = 0; itpf = 0; ittf = 0; itqvf = 0;
45
46 if (sfc_assi_options == sfc_assi_options_1) then
47 do n=iv%info(buoy)%n1,iv%info(buoy)%n2
48 ! [1.1] Get horizontal interpolation weights:
49
50 i = iv%info(buoy)%i(1,n)
51 j = iv%info(buoy)%j(1,n)
52 dx = iv%info(buoy)%dx(1,n)
53 dy = iv%info(buoy)%dy(1,n)
54 dxm = iv%info(buoy)%dxm(1,n)
55 dym = iv%info(buoy)%dym(1,n)
56
57 ! Surface correction
58
59 iv%buoy(n)%p%inv = ob%buoy(n)%p
60 iv%buoy(n)%t%inv = ob%buoy(n)%t
61 iv%buoy(n)%q%inv = ob%buoy(n)%q
62 iv%buoy(n)%u%inv = ob%buoy(n)%u
63 iv%buoy(n)%v%inv = ob%buoy(n)%v
64
65 if (iv % buoy(n) % h > missing_r) then
66 do k=kts,kte
67 v_h(k) = dym*(dxm*grid%xb%h(i,j ,k) + dx*grid%xb%h(i+1,j ,k)) &
68 + dy *(dxm*grid%xb%h(i,j+1,k) + dx*grid%xb%h(i+1,j+1,k))
69 end do
70
71 hd = v_h(kts) - iv % buoy(n) % h
72
73 if (abs(hd) <= Max_StHeight_Diff .or. anal_type_verify) then
74 if (iv % buoy(n) % h < v_h(kts)) then
75 iv%info(buoy)%zk(:,n) = 1.0+1.0e-6
76 call da_obs_sfc_correction(iv%info(buoy), iv%buoy(n), n, grid%xb)
77
78 ! To keep the original "ob" with no change for multiple
79 ! outer-loops use:
80
81 ! ob%buoy(n)%p = iv%buoy(n)%p%inv
82 ! ob%buoy(n)%t = iv%buoy(n)%t%inv
83 ! ob%buoy(n)%q = iv%buoy(n)%q%inv
84 ! ob%buoy(n)%u = iv%buoy(n)%u%inv
85 ! ob%buoy(n)%v = iv%buoy(n)%v%inv
86 else
87 call da_to_zk(iv % buoy(n) % h, v_h, v_interp_h, iv%info(buoy)%zk(1,n))
88 end if
89 end if
90 else if (ob % buoy(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 % buoy(n) % p, v_p, v_interp_p, iv%info(buoy)%zk(1,n))
97
98 if (iv%info(buoy)%zk(1,n) < 0.0 .and. .not.anal_type_verify) then
99 iv % buoy(n) % p % inv = missing_r
100 iv % buoy(n) % p % qc = missing
101 iv%info(buoy)%zk(:,n) = 1.0+1.0e-6
102 end if
103 end if
104 end do
105
106 call da_convert_zk (iv%info(buoy))
107
108 if (.not.anal_type_verify) then
109 do n=iv%info(buoy)%n1,iv%info(buoy)%n2
110 if (iv%info(buoy)%zk(1,n) < 0.0) then
111 iv % buoy(n) % u % qc = missing
112 iv % buoy(n) % v % qc = missing
113 iv % buoy(n) % t % qc = missing
114 iv % buoy(n) % q % qc = missing
115 iv % buoy(n) % p % qc = missing
116 end if
117 end do
118 end if
119
120 ! [1.2] Interpolate horizontally:
121 call da_interp_lin_3d (grid%xb%u, iv%info(buoy),model_u)
122 call da_interp_lin_3d (grid%xb%v, iv%info(buoy),model_v)
123 call da_interp_lin_3d (grid%xb%t, iv%info(buoy),model_t)
124 call da_interp_lin_3d (grid%xb%q, iv%info(buoy),model_q)
125 call da_interp_lin_3d (grid%xb%p, iv%info(buoy),model_p)
126 else if (sfc_assi_options == 2) then
127
128 ! Surface data assmiilation approach 2
129 ! -----------------------------------
130
131 ! 1.2.1 Surface assmiilation approach 2(10-m u, v, 2-m t, q,
132 ! and sfc_p)
133
134 call da_interp_lin_2d (grid%xb%u10, iv%info(buoy), 1,model_u)
135 call da_interp_lin_2d (grid%xb%v10, iv%info(buoy), 1,model_v)
136 call da_interp_lin_2d (grid%xb%t2, iv%info(buoy), 1,model_t)
137 call da_interp_lin_2d (grid%xb%q2, iv%info(buoy), 1,model_q)
138 call da_interp_lin_2d (grid%xb%psfc, iv%info(buoy), 1,model_p)
139
140 do n=iv%info(buoy)%n1,iv%info(buoy)%n2
141 if (iv%buoy(n)%p%qc >= 0) then
142 ! model surface p, t, q, h at observed site:
143
144 call da_interp_lin_2d_partial (grid%xb%terr, iv%info(buoy), 1, n, n, model_hsm(:,n))
145
146 ho = iv%buoy(n)%h
147 to = -888888.0
148 qo = -888888.0
149
150 if (iv%buoy(n)%t%qc >= 0 .and. iv%buoy(n)%q%qc >= 0) then
151 to = ob%buoy(n)%t
152 qo = ob%buoy(n)%q
153 call da_sfc_pre(psfcm, model_p(1,n), model_t(1,n), model_q(1,n), model_hsm(1,n), ho, to, qo)
154 else if (iv%buoy(n)%t%qc >= 0 .and. iv%buoy(n)%q%qc < 0) then
155 to = ob%buoy(n)%t
156 call da_sfc_pre(psfcm, model_p(1,n), model_t(1,n), model_q(1,n), model_hsm(1,n), ho, to)
157 else
158 call da_sfc_pre(psfcm, model_p(1,n), model_t(1,n), model_q(1,n), model_hsm(1,n), ho)
159 end if
160
161 ! Pressure at the observed height:
162 model_p(1,n) = psfcm
163 end if
164 end do
165 end if
166
167 do n=iv%info(buoy)%n1,iv%info(buoy)%n2
168
169 !-----------------------------------------------------------------------
170 ! [3.0] Fast interpolation:
171 !-----------------------------------------------------------------------
172
173 if (ob % buoy(n) % u > missing_r .AND. iv % buoy(n) % u % qc >= obs_qc_pointer) then
174 iv % buoy(n) % u % inv = iv%buoy(n)%u%inv - model_u(1,n)
175 else
176 iv % buoy(n) % u % inv = 0.0
177 end if
178
179 if (ob % buoy(n) % v > missing_r .AND. iv % buoy(n) % v % qc >= obs_qc_pointer) then
180 iv % buoy(n) % v % inv = iv%buoy(n)%v%inv - model_v(1,n)
181 else
182 iv % buoy(n) % v % inv = 0.0
183 end if
184
185 if (ob % buoy(n) % p > 0.0 .AND. iv % buoy(n) % p % qc >= obs_qc_pointer) then
186 iv % buoy(n) % p % inv = iv%buoy(n)%p%inv - model_p(1,n)
187 else
188 iv % buoy(n) % p % inv = 0.0
189 end if
190
191 if (ob % buoy(n) % t > 0.0 .AND. iv % buoy(n) % t % qc >= obs_qc_pointer) then
192 iv % buoy(n) % t % inv = iv%buoy(n)%t%inv - model_t(1,n)
193 else
194 iv % buoy(n) % t % inv = 0.0
195 end if
196
197 if (ob % buoy(n) % q > 0.0 .AND. iv % buoy(n) % q % qc >= obs_qc_pointer) then
198 iv % buoy(n) % q % inv = iv%buoy(n)%q%inv - model_q(1,n)
199 else
200 iv % buoy(n) % q % inv = 0.0
201 end if
202 end do
203
204 ! -----------------------------------------------------------------------
205 ! [5.0] Perform optional maximum error check:
206 !-----------------------------------------------------------------------
207
208 if (check_max_iv) call da_check_max_iv_buoy(iv, it, itu,ituf,itvv,itvvf,itp,itpf,itt,ittf,itqv,itqvf)
209
210 if (rootproc .and. check_max_iv_print) then
211 write(unit= check_max_iv_unit, fmt ='(A,i5,A)')'For outer iteration ',it, &
212 ', Total Rejections for Buoy follows:'
213
214 write(unit = check_max_iv_unit, fmt = '(/,10(2(A,I6),/))') &
215 'Number of failed u-wind observations: ',ituf, ' on ',itu, &
216 'Number of failed v-wind observations: ',itvvf,' on ',itvv, &
217 'Number of failed pressure observations: ',itpf, ' on ',itp, &
218 'Number of failed temperature observations:',ittf, ' on ',itt, &
219 'Number of failed mixing ratio observations:',itqvf,' on ',itqv, &
220 'Finally Total Buoy rejections ',ituf+itvvf+itpf+ittf+itqvf,' on ',&
221 itu +itvv +itp +itt +itqv
222 end if
223
224 deallocate (model_u)
225 deallocate (model_v)
226 deallocate (model_t)
227 deallocate (model_p)
228 deallocate (model_q)
229 deallocate (model_hsm)
230
231 if (trace_use_dull) call da_trace_exit("da_get_innov_vector_buoy")
232
233 end subroutine da_get_innov_vector_buoy
234
235