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