var/da/da_gpsref/da_get_innov_vector_gpsref.inc

   1 subroutine da_get_innov_vector_gpsref(it, grid, ob, iv)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: calculate innovation vectors of gpsro refractivity
   5    ! History:
   6    !   2020-06: move QC out of da_get_innov_vector_gpsref
   7    !-----------------------------------------------------------------------
   8
   9    implicit none
  10
  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.
  15
  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
  22
  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.
  31
  32    character(len=2) :: c_it
  33    character(len=7) :: c_nt
  34
  35
  36    if (trace_use_dull) call da_trace_entry("da_get_innov_vector_gpsref")
  37
  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
  43
  44    allocate (model_ref(iv%info(gpsref)%max_lev,iv%info(gpsref)%n1:iv%info(gpsref)%n2))
  45
  46    model_ref(:,:) = 0.0
  47
  48    do n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
  49
  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
  61
  62       ! Get cross pt. horizontal interpolation weights:
  63
  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)
  70
  71       if ( .not. pseudo_ref ) then
  72
  73          ! Get the zk from gpsref%h:
  74
  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))
  78 !
  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
  92
  93 !
  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 !..............................................................................
 106
 107    end do
 108
 109    call da_convert_zk (iv%info(gpsref))
 110
 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)
 113
 114    do n=iv%info(gpsref)%n1,iv%info(gpsref)%n2
 115
 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
 120
 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
 124
 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
 135
 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
 142
 143    if ( .not. anal_type_verify ) then
 144       call da_qc_gpsref(it, grid, ob, iv)
 145    end if
 146
 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
 153
 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)
 163
 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
 179
 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
 184
 185           do l=1, iv%info(gpsref)%levels(n)
 186           write(533,*) n,l,iv%gpsref(n)%ref(l)%qc
 187           enddo
 188
 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
 207
 208 ! Write out GPS Ref data:
 209
 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
 245
 246    deallocate (model_ref)
 247
 248    if (trace_use_dull) call da_trace_exit("da_get_innov_vector_gpsref")
 249
 250 end subroutine da_get_innov_vector_gpsref