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