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