Merge remote-tracking branch 'origin/release-v4.5.2'
[WRF.git] / var / da / da_gpsref / da_get_innov_vector_gpsref.inc
blob581c619815f99b37aef918e8f0b5ced6066c50d5
1 subroutine da_get_innov_vector_gpsref(it, grid, ob, iv)
3    !-----------------------------------------------------------------------
4    ! Purpose: calculate innovation vectors of gpsro refractivity
5    ! History:
6    !   2020-06: move QC out of da_get_innov_vector_gpsref
7    !-----------------------------------------------------------------------
9    implicit none
11    integer,          intent(in)    :: it       ! External iteration.
12    type(domain),     intent(in)    :: grid     ! first guess state.
13    type(y_type),     intent(inout) :: ob       ! Observation structure.
14    type(iv_type),    intent(inout) :: iv       ! O-B structure.
16    integer :: n        ! Loop counter.
17    integer :: i, j, k, kk  ! Index dimension.
18    real    :: dx, dxm, dz, dzm  ! Interpolation weights.
19    real    :: dy, dym  ! Interpolation weights.
20    real,allocatable :: model_ref(:,:) !Model gpsref at ob loc
21    real    :: v_h(kms:kme), v_p(kms:kme)     ! Model value h at ob
23    integer :: Iu_ref, l
24    integer,allocatable :: used_lev(:,:) ! for obs. data thinning
25                                         ! record the closest level with model
26    integer,allocatable :: qc(:)         ! record iv%gpsref(n)%ref(k)%qc
27                                ! hor. location.
28    real    :: distance_h       ! cal. h-difference between obs and model
29    real,allocatable :: min_dis(:)   ! minimum difference
30                                                ! hor. location.
32    character(len=2) :: c_it
33    character(len=7) :: c_nt
36    if (trace_use_dull) call da_trace_entry("da_get_innov_vector_gpsref")
38 ! all processors (even when nlocal=0) have to call da_qc_gpsref later in this subroutine
39 !   if (iv%info(gpsref)%nlocal < 1) then
40 !      if (trace_use_dull) call da_trace_exit("da_get_innov_vector_gpsref")
41 !      return
42 !   end if
44    allocate (model_ref(iv%info(gpsref)%max_lev,iv%info(gpsref)%n1:iv%info(gpsref)%n2))
46    model_ref(:,:) = 0.0
48    do n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
50       do k=1, iv%info(gpsref)%levels(n)
51          if ( ( iv%gpsref(n)%ref(k)%qc == fails_error_max )    .or. &
52               ( iv%gpsref(n)%ref(k)%qc == qcflag_pcnt_below )  .or. &
53               ( iv%gpsref(n)%ref(k)%qc == qcflag_pcnt_middle ) .or. &
54               ( iv%gpsref(n)%ref(k)%qc == qcflag_pcnt_above )  .or. &
55               ( iv%gpsref(n)%ref(k)%qc == qcflag_dndz )        .or. &
56               ( iv%gpsref(n)%ref(k)%qc == qcflag_dndz2 )       .or. &
57               ( iv%gpsref(n)%ref(k)%qc == qcflag_cutoff ) ) then
58             if( it > 1 ) iv%gpsref(n)%ref(k)%qc = 0
59          end if
60       end do
62       ! Get cross pt. horizontal interpolation weights:
64       i   = iv%info(gpsref)%i(1,n)
65       j   = iv%info(gpsref)%j(1,n)
66       dx  = iv%info(gpsref)%dx(1,n)
67       dy  = iv%info(gpsref)%dy(1,n)
68       dxm = iv%info(gpsref)%dxm(1,n)
69       dym = iv%info(gpsref)%dym(1,n)
71       if ( .not. pseudo_ref ) then
73          ! Get the zk from gpsref%h:
75          do k=kts,kte
76             v_h(k) = dym*(dxm*grid%xb%h(i,j  ,k) + dx*grid%xb%h(i+1,j  ,k)) &
77                    + dy *(dxm*grid%xb%h(i,j+1,k) + dx*grid%xb%h(i+1,j+1,k))
79             v_p(k) = dym*(dxm*grid%xb%p(i,j  ,k) + dx*grid%xb%p(i+1,j  ,k)) &
80                    + dy *(dxm*grid%xb%p(i,j+1,k) + dx*grid%xb%p(i+1,j+1,k))
81          end do
82          do k=1, iv%info(gpsref)%levels(n)
83             if (iv%gpsref(n)%h(k) > 0.0) &
84                call da_to_zk(iv%gpsref(n)%h(k), v_h, v_interp_h, iv%info(gpsref)%zk(k,n))
85             if (iv%info(gpsref)%zk(k,n) < 0.0 .and. .not. anal_type_verify) then
86                iv%gpsref(n)%ref(k)%qc = missing_data
87             end if
88          end do
89       else
90          iv%info(gpsref)%zk(:,n) = pseudo_z
91       end if
94 ! To assign the retrieved pressure to GPSREF data (YRG, 06/15/2011)............
95       do k=1, iv%info(gpsref)%levels(n)
96          kk = int (iv%info(gpsref)%zk(k,n))
97          if (kk >= kts .and. kk+1 <= kte) then
98            dz  = iv%info(gpsref)%zk(k,n) - real(kk)
99            dzm = 1.0 - dz
100            iv%gpsref(n)%p(k)%inv = v_p(kk  ) * dzm + v_p(kk+1) * dz 
101            ob%gpsref(n)%p(k) = iv%gpsref(n)%p(k)%inv
102            iv%gpsref(n)%p(k)%qc = -5
103          endif
104       end do
105 !..............................................................................
107    end do
109    call da_convert_zk (iv%info(gpsref))
111 ! t_iwabuchi 20111216 (hmn's update) linear interpolation of log(n) 
112    call da_interp_lin_3d (grid%xb%reflog, iv%info(gpsref), model_ref)
114    do n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
116 ! t_iwabuchi 20111216 compute exp of log (N)
117       do k = 1, iv%info(gpsref)%levels(n)
118         model_ref(k,n)=exp(model_ref(k,n))
119       end do
121       if ( ( .not. pseudo_ref ) .or. it > 1 ) then
122          do k = 1, iv%info(gpsref)%levels(n)
123             iv%gpsref(n)%ref(k)%inv = 0.0
125             if (ob%gpsref(n)%ref(k) > missing_r .AND. &
126                  !iv%gpsref(n)%ref(k)%qc >= obs_qc_pointer) then
127                  iv%gpsref(n)%ref(k)%qc /= missing_data) then
128                  iv%gpsref(n)%ref(k)%inv = ob%gpsref(n)%ref(k) - model_ref(k,n)
129             end if
130          end do
131       else
132          ob % gpsref(1)%ref(1) = model_ref(1,n) + iv %gpsref(1)%ref(1)%inv 
133       end if
134    end do
136    if ( pseudo_ref ) then
137       ! Done for pseudo_ref after getting iv%gpsref(n)%ref(k)%inv
138       deallocate (model_ref)
139       if (trace_use_dull) call da_trace_exit("da_get_innov_vector_gpsref")
140       return
141    end if
143    if ( .not. anal_type_verify ) then
144       call da_qc_gpsref(it, grid, ob, iv)
145    end if
147 ! ------------------------------------------------------------------------------
148 !   GPSRO thinning  (Shu-Ya Chen 20090701)
149    if (.not. anal_type_verify) then
150     IF ( gpsref_thinning ) THEN
151       allocate (used_lev(kms:kme,iv%info(gpsref)%n1:iv%info(gpsref)%n2))
152       used_lev(:,:) = missing_data
154        DO n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
155        allocate(min_dis(kms:kme))
156        allocate(qc(iv%info(gpsref)%levels(n)))
157        i   = iv%info(gpsref)%i(1,n)
158        j   = iv%info(gpsref)%j(1,n)
159        dx  = iv%info(gpsref)%dx(1,n)
160        dy  = iv%info(gpsref)%dy(1,n)
161        dxm = iv%info(gpsref)%dxm(1,n)
162        dym = iv%info(gpsref)%dym(1,n)
164        ! Get the zk from gpsref%h:
165           do k=kts,kte
166              v_h(k) = dym*(dxm*grid%xb%h(i,j  ,k) + dx*grid%xb%h(i+1,j  ,k)) &
167                     + dy *(dxm*grid%xb%h(i,j+1,k) + dx*grid%xb%h(i+1,j+1,k))
168           end do
169           do k=kts,kte 
170           min_dis(k)=1.0E10
171              do l=1, iv%info(gpsref)%levels(n)
172                 if ( iv%gpsref(n)%ref(l)%qc >= obs_qc_pointer ) then
173                 distance_h=abs(iv%gpsref(n)%h(l)-v_h(k))
174                 min_dis(k)=min(min_dis(k),distance_h)
175                 if ( min_dis(k) == distance_h ) used_lev(k,n)=l
176                 end if
177              end do
178           end do
180           write(533,*) 'obs_qc_pointer=',obs_qc_pointer,'missing_data=',missing_data
181           do k=kts,kte
182           write(533,*) n,k,'used_lev=',used_lev(k,n)
183           enddo
185           do l=1, iv%info(gpsref)%levels(n)
186           write(533,*) n,l,iv%gpsref(n)%ref(l)%qc
187           enddo
189           do l=1, iv%info(gpsref)%levels(n)
190           qc(l)=iv%gpsref(n)%ref(l)%qc
191           end do
192           do k=kts,kte
193            qc(used_lev(k,n))=1   ! which level is closest to model level
194           end do
195        !  data thinning (set thinned levels to be -99)
196           do l=1, iv%info(gpsref)%levels(n)
197           if ( iv%gpsref(n)%ref(l)%qc >= obs_qc_pointer &
198                .and. qc(l) /= 1 ) then
199           iv%gpsref(n)%ref(l)%qc = -99
200           end if
201           end do
202        deallocate(min_dis)
203        deallocate(qc)
204        END DO
205        deallocate (used_lev)
206     END IF
208 ! Write out GPS Ref data:
210    if ( write_iv_gpsref ) then
211      DO n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
212      Iu_ref = 336
213      write(c_it,'(I2.2)') it
214      c_nt=iv%info(gpsref)%name(n)(8:11)//iv%info(gpsref)%name(n)(28:30)
215      open (unit=Iu_ref, file='RO_Innov_'//iv%info(gpsref)%date_char(n)//'_'//c_nt//'.'//c_it, &
216            form='formatted')
217            write(unit=Iu_ref, fmt='(/i5,2x,a,2x,a,2x,4f10.3,i5)') n, &
218               iv%info(gpsref)%date_char(n), iv%info(gpsref)%id(n),  &
219               iv%info(gpsref)%lat(1,n)  , iv%info(gpsref)%lon(1,n), &
220               iv%info(gpsref)%x(1,n)  , iv%info(gpsref)%y(1,n), &
221               iv%info(gpsref)%levels(n)
222            write(unit=Iu_ref, fmt='(a5,3x,6a14)') 'level','     height   ', &
223                        '    Obs_ref   ','  model_ref   ','  Innov_ref   ', &
224                        '  error_ref   ',' qc_ref       '
225            do k = 1, iv%info(gpsref)%levels(n)
226 !             if ( gpsref_thinning ) then
227 !               if ( iv%gpsref(n)%ref(l)%qc >= obs_qc_pointer ) then
228 !                  write(unit=Iu_ref, fmt='(i3,1x,5f14.3,i10)')  k, &
229 !                  iv%gpsref(n)%h(k), ob%gpsref(n)%ref(k),   &
230 !                  model_ref(k,n), iv%gpsref(n)%ref(k)%inv,   &
231 !                  iv%gpsref(n)%ref(k)%error, iv%gpsref(n)%ref(k)%qc
232 !               end if
233 !             else
234                write(unit=Iu_ref, fmt='(i3,1x,5f14.3,i10)')  k, &
235                   iv%gpsref(n)%h(k), ob%gpsref(n)%ref(k),   &
236                   model_ref(k,n), iv%gpsref(n)%ref(k)%inv,   &
237                   iv%gpsref(n)%ref(k)%error, iv%gpsref(n)%ref(k)%qc
238 !             end if
239            end do
240      close(Iu_ref)
241      END DO
242    end if ! write_iv_gpsref
243    ! .........................................................................
244    end if  ! end of verify check
246    deallocate (model_ref)
248    if (trace_use_dull) call da_trace_exit("da_get_innov_vector_gpsref")
250 end subroutine da_get_innov_vector_gpsref