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
41 itu = 0; itvv = 0; itp = 0; itt = 0; itqv = 0;
42 ituf = 0; itvvf = 0; itpf = 0; ittf = 0; itqvf = 0;
43
44 do n=iv%ob_numb(iv%current_ob_time-1)%buoy + 1, &
45 iv%ob_numb(iv%current_ob_time)%buoy
46 ! [1.1] Get horizontal interpolation weights:
47
48 i = iv%buoy(n)%loc%i
49 j = iv%buoy(n)%loc%j
50 dx = iv%buoy(n)%loc%dx
51 dy = iv%buoy(n)%loc%dy
52 dxm = iv%buoy(n)%loc%dxm
53 dym = iv%buoy(n)%loc%dym
54
55 ! Surface correction
56
57 iv%buoy(n)%p%inv = ob%buoy(n)%p
58 iv%buoy(n)%t%inv = ob%buoy(n)%t
59 iv%buoy(n)%q%inv = ob%buoy(n)%q
60 iv%buoy(n)%u%inv = ob%buoy(n)%u
61 iv%buoy(n)%v%inv = ob%buoy(n)%v
62
63 if (sfc_assi_options == 1) then
64 iv%buoy(n)%zk=missing_r
65
66 if (iv % buoy(n) % h > missing_r) then
67 do k=xp%kts,xp%kte
68 v_h(k) = dym*(dxm*xb%h(i,j ,k) + dx*xb%h(i+1,j ,k)) &
69 + dy *(dxm*xb%h(i,j+1,k) + dx*xb%h(i+1,j+1,k))
70 end do
71
72 hd = v_h(xp%kts) - iv % buoy(n) % h
73
74 if (abs(hd) <= Max_StHeight_Diff .or. anal_type_verify) then
75 if (iv % buoy(n) % h < v_h(xp%kts)) then
76 iv%buoy(n)%zk = 1.0+1.0e-6
77 call da_obs_sfc_correction(iv%buoy(n), xb)
78
79 ! To keep the original "ob" with no change for multiple
80 ! outer-loops use:
81
82 ! ob%buoy(n)%p = iv%buoy(n)%p%inv
83 ! ob%buoy(n)%t = iv%buoy(n)%t%inv
84 ! ob%buoy(n)%q = iv%buoy(n)%q%inv
85 ! ob%buoy(n)%u = iv%buoy(n)%u%inv
86 ! ob%buoy(n)%v = iv%buoy(n)%v%inv
87 else
88 call da_to_zk(iv % buoy(n) % h, v_h, xp, v_interp_h, &
89 iv%buoy(n)%zk)
90 end if
91 else
92 iv%buoy(n)%zk = missing_r
93 end if
94 else if (ob % buoy(n) % p > 1.0) then
95 do k=xp%kts,xp%kte
96 v_p(k) = dym*(dxm*xb%p(i,j ,k) + dx*xb%p(i+1,j ,k)) &
97 + dy *(dxm*xb%p(i,j+1,k) + dx*xb%p(i+1,j+1,k))
98 end do
99
100 call da_to_zk(ob % buoy(n) % p, v_p, xp, v_interp_p, &
101 iv%buoy(n)%zk)
102
103 if (iv%buoy(n)%zk < 0.0 .and. .not.anal_type_verify) then
104 iv % buoy(n) % p % inv = missing_r
105 iv % buoy(n) % p % qc = missing
106
107 iv%buoy(n)%zk = 1.0+1.0e-6
108 end if
109 end if
110
111 !------------------------------------------------------------------
112 ! [2.0] Initialise components of innovation vector:
113 !------------------------------------------------------------------
114
115 if (iv%buoy(n)%zk < 0.0 .and. .not.anal_type_verify) then
116 iv % buoy(n) % u % qc = missing
117 iv % buoy(n) % v % qc = missing
118 iv % buoy(n) % t % qc = missing
119 iv % buoy(n) % q % qc = missing
120 iv % buoy(n) % p % qc = missing
121 else
122 ! [1.2] Interpolate horizontally:
123 call da_interp_obs_lin_2d( xb % u, xp, i, j, dx, dy, dxm, dym, &
124 model_u, iv%buoy(n)%zk)
125 call da_interp_obs_lin_2d( xb % v, xp, i, j, dx, dy, dxm, dym, &
126 model_v, iv%buoy(n)%zk)
127 call da_interp_obs_lin_2d( xb % t, xp, i, j, dx, dy, dxm, dym, &
128 model_t, iv%buoy(n)%zk)
129 call da_interp_obs_lin_2d( xb % q, xp, i, j, dx, dy, dxm, dym, &
130 model_q, iv%buoy(n)%zk)
131 call da_interp_obs_lin_2d( xb % p, xp, i, j, dx, dy, dxm, dym, &
132 model_p, iv%buoy(n)%zk)
133 end if
134 else if (sfc_assi_options == 2) then
135
136 ! Surface data assimilation approca 2
137 ! -----------------------------------
138
139 ! 1.2.1 Surface assimilation approach 2(10-m u, v, 2-m t, q,
140 ! and sfc_p)
141
142 call da_interp_lin_2d( xb % u10, xp%ims, xp%ime, xp%jms, xp%jme, &
143 i, j, dx, dy, dxm, dym, &
144 model_u)
145 call da_interp_lin_2d( xb % v10, xp%ims, xp%ime, xp%jms, xp%jme, &
146 i, j, dx, dy, dxm, dym, &
147 model_v)
148 call da_interp_lin_2d( xb % t2, xp%ims, xp%ime, xp%jms, xp%jme, &
149 i, j, dx, dy, dxm, dym, &
150 model_t)
151 call da_interp_lin_2d( xb % q2, xp%ims, xp%ime, xp%jms, xp%jme, &
152 i, j, dx, dy, dxm, dym, &
153 model_q)
154 call da_interp_lin_2d( xb % psfc, xp%ims, xp%ime, xp%jms, xp%jme, &
155 i, j, dx, dy, dxm, dym, &
156 model_p)
157
158 if (iv%buoy(n)%p%qc >= 0) then
159 ! model surface p, t, q, h at observed site:
160
161 call da_interp_lin_2d( xb % terr, xp%ims, xp%ime, xp%jms,xp%jme,&
162 i, j, dx, dy, dxm, dym, &
163 hsm)
164
165 ho = iv%buoy(n)%h
166 to = -888888.0
167 qo = -888888.0
168
169 if (iv%buoy(n)%t%qc >= 0 .and. iv%buoy(n)%q%qc >= 0) then
170 to = ob%buoy(n)%t
171 qo = ob%buoy(n)%q
172 call da_sfc_pre(psfcm, model_p, model_t, model_q, &
173 hsm, ho, to, qo)
174 else if (iv%buoy(n)%t%qc >= 0 .and. iv%buoy(n)%q%qc < 0) then
175 to = ob%buoy(n)%t
176 call da_sfc_pre(psfcm, model_p, model_t, model_q, hsm, ho, to)
177 else
178 call da_sfc_pre(psfcm, model_p, model_t, model_q, hsm, ho)
179 end if
180
181 ! Pressure at the observed height:
182 model_p = psfcm
183 end if
184 end if
185
186 !-----------------------------------------------------------------------
187 ! [3.0] Fast interpolation:
188 !-----------------------------------------------------------------------
189
190 if (ob % buoy(n) % u > missing_r .AND. &
191 iv % buoy(n) % u % qc >= obs_qc_pointer) then
192 iv % buoy(n) % u % inv = iv%buoy(n)%u%inv - model_u
193 else
194 iv % buoy(n) % u % inv = 0.0
195 end if
196
197 if (ob % buoy(n) % v > missing_r .AND. &
198 iv % buoy(n) % v % qc >= obs_qc_pointer) then
199 iv % buoy(n) % v % inv = iv%buoy(n)%v%inv - model_v
200 else
201 iv % buoy(n) % v % inv = 0.0
202 end if
203
204 if (ob % buoy(n) % p > 0.0 .AND. &
205 iv % buoy(n) % p % qc >= obs_qc_pointer) then
206 iv % buoy(n) % p % inv = iv%buoy(n)%p%inv - model_p
207 else
208 iv % buoy(n) % p % inv = 0.0
209 end if
210
211 if (ob % buoy(n) % t > 0.0 .AND. &
212 iv % buoy(n) % t % qc >= obs_qc_pointer) then
213 iv % buoy(n) % t % inv = iv%buoy(n)%t%inv - model_t
214 else
215 iv % buoy(n) % t % inv = 0.0
216 end if
217
218 if (ob % buoy(n) % q > 0.0 .AND. &
219 iv % buoy(n) % q % qc >= obs_qc_pointer) then
220 iv % buoy(n) % q % inv = iv%buoy(n)%q%inv - model_q
221 else
222 iv % buoy(n) % q % inv = 0.0
223 end if
224
225 ! -----------------------------------------------------------------------
226 ! [5.0] Perform optional maximum error check:
227 !-----------------------------------------------------------------------
228
229 if (check_max_iv) then
230 call da_check_max_iv_buoy(it, iv % buoy(n), &
231 itu,ituf,itvv,itvvf,itp,itpf,itt,ittf,itqv,itqvf)
232 end if
233 end do
234
235 if (rootproc .and. check_max_iv_print) then
236 write(unit= check_max_iv_unit, fmt ='(A,i5,A)')'For outer iteration ',it, &
237 ', Total Rejections for Buoy follows:'
238
239 write(unit = check_max_iv_unit, fmt = '(/,10(2(A,I6),/))') &
240 'Number of failed u-wind observations: ',ituf, ' on ',itu, &
241 'Number of failed v-wind observations: ',itvvf,' on ',itvv, &
242 'Number of failed pressure observations: ',itpf, ' on ',itp, &
243 'Number of failed temperature observations:',ittf, ' on ',itt, &
244 'Number of failed mixing ratio observations:',itqvf,' on ',itqv, &
245 'Finally Total Buoy rejections ',ituf+itvvf+itpf+ittf+itqvf,' on ',&
246 itu +itvv +itp +itt +itqv
247 end if
248
249 end subroutine da_get_innov_vector_buoy
250
251