da_get_innov_vector_gpsref.inc
References to this file elsewhere.
1 subroutine da_get_innov_vector_gpsref( it, xb, xp, ob, iv)
2
3 !-----------------------------------------------------------------------
4 ! Purpose: TBD
5 !-----------------------------------------------------------------------
6
7 ! ------------------------------------------------------------------------
8 ! Called by da_obs/da_get_innov_vector.inc
9 ! ------------------------------------------------------------------------
10
11 implicit none
12
13 integer, intent(in) :: it ! External iteration.
14 type(xb_type), intent(in) :: xb ! first guess state.
15 type(xpose_type), intent(in) :: xp ! Domain decomposition vars.
16 type(y_type), intent(inout) :: ob ! Observation structure.
17 type(ob_type), intent(inout) :: iv ! O-B structure.
18
19 integer :: n ! Loop counter.
20 integer :: i, j, k ! Index dimension.
21 integer :: num_levs ! Number of obs levels.
22 real :: dx, dxm ! Interpolation weights.
23 real :: dy, dym ! Interpolation weights.
24 real, dimension(1:max_ob_levels) :: model_ref !Model gpsref at ob loc
25 real, dimension(xp%kms:xp%kme) :: v_h ! Model value h at ob
26 ! hor. location.
27 integer :: pref, fref
28
29 ! For quality control
30
31 real , parameter :: h_1 = 7000.0, h_2 = 25000.0
32 ! Lidia Cucurull values:
33 real , parameter :: pcnt1 = 0.05, pcnt2 = 0.04, pcnt3 = 0.10
34 ! testing values:
35 ! real , parameter :: pcnt1 = 0.02, pcnt2 = 0.01, pcnt3 = 0.03
36 integer, parameter :: qc_below = -31, qc_middle = -32, qc_above = -33
37
38 integer :: nn, n1, ntotal, nqc0, nqc1, nqc2, nqc3
39 real :: percnt
40 character(len=40), dimension(5000) :: name_qc
41 real , dimension(5000) :: height_below
42
43
44 ! GPS REF Pseudo OBS test:
45
46 if (pseudo_var(1:3) == 'ref' .and. num_pseudo > 0) then
47
48 ! Deallocate:
49 if (iv%num_gpsref > 0) then
50 do n = 1, iv%num_gpsref
51 deallocate(iv % gpsref(n) % h)
52 deallocate(iv % gpsref(n) % zk)
53 deallocate(iv % gpsref(n) % ref)
54 deallocate(iv % gpsref(n) % p)
55 deallocate(iv % gpsref(n) % t)
56 deallocate(iv % gpsref(n) % q)
57 deallocate(ob % gpsref(n) % ref)
58 end do
59 deallocate(iv % gpsref)
60 deallocate(ob % gpsref)
61 end if
62
63 use_gpsrefobs = .true.
64
65 ! Allocate:
66 iv%num_gpsref = num_pseudo
67 iv%ob_numb(1)%gpsref = num_pseudo
68 iv%num_pseudo = 0
69
70 allocate(iv % gpsref(1:num_pseudo))
71 allocate(iv%gpsref(num_pseudo)%zk(1:1))
72 allocate(iv%gpsref(num_pseudo)%ref(1:1))
73 allocate(ob%gpsref(1:num_pseudo))
74 allocate(ob%gpsref(num_pseudo)%ref(1:1))
75
76 write(stdout,'(a,i2)') '==> GPS REF pseudo OBS test: num_pseudo=',num_pseudo
77
78 iv%gpsref(1)%info%levels = 1
79
80 iv % gpsref(1)%loc%x = pseudo_x
81 iv % gpsref(1)%loc%y = pseudo_y
82
83 iv%gpsref(1)%loc%i = int(pseudo_x)
84 iv%gpsref(1)%loc%j = int(pseudo_y)
85
86 iv%gpsref(1)%loc%dx = pseudo_x-real(iv%gpsref(1)%loc%i)
87 iv%gpsref(1)%loc%dy = pseudo_y-real(iv%gpsref(1)%loc%j)
88 iv%gpsref(1)%loc%dxm=1.0-iv%gpsref(1)%loc%dx
89 iv%gpsref(1)%loc%dym=1.0-iv%gpsref(1)%loc%dy
90
91 iv % gpsref(1) %ref(1) % inv = pseudo_val
92 iv % gpsref(1) %ref(1) % error = pseudo_err
93 iv % gpsref(1) %ref(1) % qc = 0
94
95 ! Set halo:
96 if ((iv%gpsref(1)%loc%i < xp%its-1) .or.(iv%gpsref(1)%loc%i > xp%ite) .or. &
97 (iv%gpsref(1)%loc%j < xp%jts-1) .or.(iv%gpsref(1)%loc%j > xp%jte)) then
98 iv%gpsref(1)%loc%proc_domain = .false.
99 else
100 iv%gpsref(1)%loc%proc_domain = .false.
101
102 if (iv%gpsref(1)%loc%i >= xp%its .and. iv%gpsref(1)%loc%i <= xp%ite .and. &
103 iv%gpsref(1)%loc%j >= xp%jts .and. iv%gpsref(1)%loc%j <= xp%jte) then
104 iv%gpsref(1)%loc%proc_domain = .true.
105 end if
106 end if
107
108 write(stdout,'(a4,2f15.5)') pseudo_var, pseudo_val, pseudo_err
109 write(stdout,'(3f15.5)') pseudo_x, pseudo_y, pseudo_z
110 end if
111
112 if (iv%num_gpsref < 1) return
113
114 ntotal = 0
115 pref = 0 ; fref = 0
116
117 do n=iv%ob_numb(iv%current_ob_time-1)%gpsref + 1, iv%ob_numb(iv%current_ob_time)%gpsref
118 num_levs = iv%gpsref(n)%info%levels
119
120 if (num_levs < 1) cycle
121
122 model_ref(:) = 0.0
123
124 ! Get cross pt. horizontal interpolation weights:
125
126 i = iv%gpsref(n)%loc%i
127 j = iv%gpsref(n)%loc%j
128 dx = iv%gpsref(n)%loc%dx
129 dy = iv%gpsref(n)%loc%dy
130 dxm = iv%gpsref(n)%loc%dxm
131 dym = iv%gpsref(n)%loc%dym
132
133 if (.not.(pseudo_var(1:3) == 'ref' .and. num_pseudo > 0)) then
134
135 ! Get the gpsref%zk from gpsref%h:
136
137 do k=xp%kts,xp%kte
138 v_h(k) = dym*(dxm*xb%h(i,j ,k) + dx*xb%h(i+1,j ,k)) &
139 + dy *(dxm*xb%h(i,j+1,k) + dx*xb%h(i+1,j+1,k))
140 end do
141 do k=1, num_levs
142 iv%gpsref(n)%zk(k)=missing_r
143 if (iv%gpsref(n)%h(k) > 0.0) &
144 call da_to_zk(iv%gpsref(n)%h(k), v_h, xp, v_interp_h, iv%gpsref(n)%zk(k))
145 if (iv%gpsref(n)%zk(k) < 0.0 .and. .not. anal_type_verify) then
146 iv%gpsref(n)%ref(k)%qc = missing
147 end if
148 end do
149 else
150 iv % gpsref(1)%zk(1) = pseudo_z
151 end if
152
153 call da_interp_lin_3d( xb%ref, xp, i, j, dx, dy, dxm, dym, &
154 model_ref, max_ob_levels, iv%gpsref(n)%zk, num_levs)
155
156 if (.not.(pseudo_var(1:3) == 'ref' .and. num_pseudo > 0)) then
157 do k = 1, num_levs
158 iv%gpsref(n)%ref(k)%inv = 0.0
159
160 if (ob%gpsref(n)%ref(k) > missing_r .AND. &
161 iv%gpsref(n)%ref(k)%qc >= obs_qc_pointer) then
162 iv%gpsref(n)%ref(k)%inv = ob%gpsref(n)%ref(k) - model_ref(k)
163 end if
164 end do
165 else
166 ob % gpsref(1)%ref(1) = model_ref(1) + iv %gpsref(1)%ref(1)%inv
167 end if
168
169 ! Quality check 1: Gross error(departure from the background) check
170
171 if (check_max_iv) &
172 call da_check_max_iv_gpsref(it, iv%gpsref(n), pref, fref)
173
174 ! Quality check 2: Error percentage check.
175
176 if (.not. anal_type_verify) then
177 if (.not.(pseudo_var(1:3) == 'ref' .and. num_pseudo > 0)) then
178 do k=1, num_levs
179
180 ! incremetal refractivity or the relative error:
181 ! abs[(O-B)/{(O+B)/2}] (Lidia Cucurull 2005)
182
183 ntotal = ntotal + 1
184 percnt = 2.0 * abs(iv%gpsref(n)%ref(k)%inv / &
185 (ob%gpsref(n)%ref(k) + model_ref(k)))
186
187 if (iv%gpsref(n)%ref(k)%qc >= obs_qc_pointer) then
188
189 if (iv%gpsref(n)%h(k) < h_1) then
190 if (percnt > pcnt1) iv%gpsref(n)%ref(k)%qc = qc_below
191 else if (iv%gpsref(n)%h(k) > h_2) then
192 if (percnt > pcnt3) iv%gpsref(n)%ref(k)%qc = qc_above
193 else
194 if (percnt > pcnt2) iv%gpsref(n)%ref(k)%qc = qc_middle
195 end if
196 end if
197 end do
198 end if
199 end if ! end of if verify check
200 end do
201
202 ! Quality check 3: Low levels quality control
203
204 if (.not. anal_type_verify) then
205 if (.not.(pseudo_var(1:3) == 'ref' .and. num_pseudo > 0)) then
206 ! Search for the GPS RO's name with the 'qc_below':
207
208 nn = 0
209 height_below = 0.0
210 name_qc = ' '
211
212 do n=iv%ob_numb(iv%current_ob_time-1)%gpsref + 1, iv%ob_numb(iv%current_ob_time)%gpsref
213 nn = nn + 1
214 num_levs = iv%gpsref(n)%info%levels
215 do k=1, num_levs
216 if (iv%gpsref(n)%ref(k)%qc == qc_below) then
217 name_qc(nn) = iv%gpsref(n)%info%name
218 height_below(nn) = max(iv%gpsref(n)%h(k),height_below(nn))
219 end if
220 end do
221 if (height_below(nn) == 0.0) nn = nn - 1
222 end do
223
224 ! Set the flag qc_below to the levels below percnt < pcnt1::
225
226 ntotal = 0
227 nqc0 = 0
228 nqc1 = 0
229 nqc2 = 0
230 nqc3 = 0
231
232 do n=iv%ob_numb(iv%current_ob_time-1)%gpsref + 1, iv%ob_numb(iv%current_ob_time)%gpsref
233 num_levs = iv%gpsref(n)%info%levels
234
235 do n1 = 1,nn
236 if (iv%gpsref(n)%info%name == name_qc(n1)) then
237 do k=1, num_levs
238 if (iv%gpsref(n)%h(k) < height_below(n1) .and. &
239 iv%gpsref(n)%ref(k)%qc >= 0) iv%gpsref(n)%ref(k)%qc = qc_below
240 end do
241 exit
242 end if
243 end do
244
245 do k=1, num_levs
246 ntotal = ntotal + 1
247 if (iv%gpsref(n)%ref(k)%qc == fails_error_max) nqc0 = nqc0 + 1
248 if (iv%gpsref(n)%ref(k)%qc == qc_middle) nqc1 = nqc1 + 1
249 if (iv%gpsref(n)%ref(k)%qc == qc_below) nqc2 = nqc2 + 1
250 if (iv%gpsref(n)%ref(k)%qc == qc_above) nqc3 = nqc3 + 1
251 end do
252 end do
253 end if
254 end if ! end of if verify check
255
256 if (rootproc .and. check_max_iv_print) then
257 write(unit = check_max_iv_unit, fmt ='(A,i5,A)') &
258 'For outer iteration ',it, &
259 ', Total Rejections for GPSRef follows:'
260 write(unit = check_max_iv_unit, fmt = '(/,10(2(A,I6),/))') &
261 'Number of failed GPSRef observations:', &
262 fref, ' on ',pref
263 end if
264
265 end subroutine da_get_innov_vector_gpsref
266
267