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