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