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