var/da/da_gpseph/da_get_innov_vector_gpseph.inc

   1 subroutine da_get_innov_vector_gpseph(it, ob, iv)
   2 !==========================================================================
   3 !  Purpose: Calculate invovation (O-B) for gpseph
   4 !  Author:  Shu-Ya Chen/Shu-Hua Chen        Date: 01/27/2011
   5 !  History: December 2015: re-implementation (Jamie Bresch)
   6 !=========================================================================
   7    implicit none
   8
   9    integer,          intent(in)    :: it       ! External iteration
  10    type(y_type),     intent(inout) :: ob       ! Observation structure
  11    type(iv_type),    intent(inout) :: iv       ! O-B structure
  12
  13    integer :: n        ! Loop counter
  14    integer :: i, j, k  ! Index dimension
  15    real, dimension(:,:), allocatable :: model_eph !Model gpseph at ob loc
  16    real, dimension(:,:), allocatable :: model_ref !Model gpsref at ob loc
  17    real, dimension(kds:kde)          :: obs_ref   !obs gpsref at ob loc
  18    real, dimension(kds:kde)          :: obs_eph   !obs gpseph at ob loc
  19    real, dimension(kds:kde)          :: mean_h
  20    real, dimension(:), allocatable   :: err_height,err_percent,intp_err_percent
  21    real, dimension(:,:), allocatable :: cc
  22    real    :: tmp_height, tmp_percent
  23    integer :: kk, count_lev, err_lev
  24    integer :: nlevels
  25    integer :: n1, n2
  26    integer :: iunit, junit, iost
  27    character(len=256) :: fname
  28    logical :: write_iv_gpseph = .true. !.false.
  29
  30    if (trace_use_dull) call da_trace_entry("da_get_innov_vector_gpseph")
  31 !=============================================================================
  32
  33    mean_h(:) = global_h_mean(:) !km
  34
  35    if (iv%info(gpseph)%nlocal < 1) return
  36
  37    nlevels = kde-kds+1
  38    allocate (model_eph(nlevels,iv%info(gpseph)%n1:iv%info(gpseph)%n2))
  39    model_eph(:,:) = 0.0
  40
  41    do n=iv%info(gpseph)%n1,iv%info(gpseph)%n2
  42       do k=1, iv%info(gpseph)%levels(n)
  43          if( iv%gpseph(n)%eph(k)%qc == fails_error_max .and. it > 1 ) &
  44              iv%gpseph(n)%eph(k)%qc = 0
  45       end do
  46    end do
  47
  48 ! Before gross error check, fill in eph data into iv and ob
  49
  50    if ( it > 1 ) then
  51       if ( allocated (gps_rays ) ) then
  52          n1 = Lbound ( gps_rays, 1 )
  53          n2 = ubound ( gps_rays, 1 )
  54          do n = n1, n2
  55             if ( allocated (gps_rays(n)%je2) ) deallocate (gps_rays(n)%je2)
  56             if ( gps_rays(n)%nbot == 0 .and. gps_rays(n)%ntop == 0 ) cycle
  57             do k = gps_rays(n)%nbot, gps_rays(n)%ntop
  58                deallocate ( gps_rays(n)%ip123(k)%i1 )
  59                deallocate ( gps_rays(n)%ip123(k)%i2 )
  60                deallocate ( gps_rays(n)%ip123(k)%i3 )
  61                deallocate ( gps_rays(n)%ip123(k)%w1 )
  62                deallocate ( gps_rays(n)%ip123(k)%w2 )
  63             enddo
  64             deallocate (gps_rays(n)%ip123)
  65          enddo
  66          deallocate (gps_rays)
  67       endif
  68    endif
  69
  70    if ( .not. allocated ( gps_rays ) ) then
  71       allocate ( gps_rays(iv%info(gpseph)%n1:iv%info(gpseph)%n2) )
  72       gps_rays%nbot = 0 !initialize
  73       gps_rays%ntop = 0 !initialize
  74    end if
  75
  76    obs_loop: do n=iv%info(gpseph)%n1,iv%info(gpseph)%n2
  77
  78       if ( iv%info(gpseph)%levels(n) < 1 ) cycle
  79
  80       call da_gpseph_rays(gps_rays(n), iv%gpseph(n))
  81
  82       obs_ref(:) = iv%gpseph(n)%ref(:)%inv
  83       call da_obs_ref_to_eph(gps_rays(n), obs_ref, obs_eph)
  84
  85       call da_mdl_ref_to_eph(gps_rays(n), global_ref_mean_h, model_eph(:,n))
  86
  87       do k=kds,kde
  88          ob%gpseph(n)%eph(k)=obs_eph(k)
  89       enddo
  90
  91       if (iv%gpseph(n)%level1 .gt. iv%gpseph(n)%level2) cycle obs_loop
  92
  93       ! set quality control (QC) ---------------------------------------------------
  94       do k = kds, kde
  95          if ((k >= iv%gpseph(n)%level1) .and. (k <= iv%gpseph(n)%level2)) then
  96             iv%gpseph(n)%eph(k)%qc = 0
  97          else
  98             iv%gpseph(n)%eph(k)%qc = missing_data
  99          endif
 100       enddo
 101
 102       do k = kds, kde
 103          iv%gpseph(n)%eph(k)%inv = 0.0
 104          if (ob%gpseph(n)%eph(k) > missing_r .AND. &
 105              iv%gpseph(n)%eph(k)%qc >= obs_qc_pointer) then
 106             iv%gpseph(n)%eph(k)%inv = ob%gpseph(n)%eph(k) - model_eph(k,n)
 107          end if
 108       end do
 109
 110    end do obs_loop
 111
 112    !----------------------------------------------------------------
 113    ! Reading observation errors from an external table.
 114    ! The percentage values in the table were generated using
 115    ! Hollingsworth-Lonngberg method from one-month
 116    ! (2003/08/15~2003/09/15) statistics
 117    !----------------------------------------------------------------
 118    ! OBSERR = percentage(%) * Excess Phase
 119    call da_get_unit(iunit)
 120    open (unit=iunit, file='OBSERROR_GPSEPH.TBL', form='formatted', &
 121       status='old', iostat=iost)
 122    if ( iost /= 0 ) then
 123       call da_free_unit(iunit)
 124       call da_error(__FILE__,__LINE__, (/"Error opening OBSERROR_GPSEPH.TBL"/))
 125    else
 126       allocate(err_height(nlevels))
 127       allocate(err_percent(nlevels))
 128       err_percent(:) = missing_r
 129       read(iunit,*)  !the header line
 130       count_lev = 0
 131       read_loop: do
 132          read(iunit,*,iostat=iost) err_lev, tmp_height, tmp_percent
 133          if ( iost /=0 ) exit read_loop
 134          count_lev = count_lev + 1
 135          err_height(count_lev)  = tmp_height
 136          err_percent(count_lev) = tmp_percent*0.01
 137       end do read_loop
 138       call da_free_unit(iunit)
 139       close(iunit)
 140       if ( count_lev < 1 ) then
 141          call da_error(__FILE__,__LINE__, (/"Error reading OBSERROR_GPSEPH.TBL"/))
 142       end if
 143    end if
 144
 145    allocate(intp_err_percent(kde-kds+1))
 146    allocate(cc(3,count_lev))
 147    call da_splinx(count_lev,err_height(1:count_lev),err_percent(1:count_lev),kde-kds+1,mean_h(kds:kde),cc,intp_err_percent)
 148
 149    do n=iv%info(gpseph)%n1,iv%info(gpseph)%n2
 150       if ( iv%info(gpseph)%levels(n) < 2 ) cycle
 151       iv%gpseph(n)%eph(:)%error = missing_r
 152       do kk = iv%gpseph(n)%level1, iv%gpseph(n)%level2
 153          iv%gpseph(n)%eph(kk)%error = intp_err_percent(kk)*ob%gpseph(n)%eph(kk)
 154       enddo
 155    enddo
 156    deallocate(intp_err_percent,err_height,err_percent,cc)
 157
 158    ! Quality check: Gross error(departure from the background) check
 159
 160    if ( check_max_iv ) then
 161       call da_check_max_iv_gpseph(iv, it)
 162    end if
 163
 164    if ( write_iv_gpseph ) then
 165       if (.not. anal_type_verify) then
 166          ! Write out GPS EPH data:
 167          call da_get_unit(junit)
 168          do n=iv%info(gpseph)%n1,iv%info(gpseph)%n2
 169             if ( iv%info(gpseph)%levels(n) < 2 ) cycle
 170             write(fname,'(a,i2.2)') 'RO_Innov_'//iv%info(gpseph)%date_char(n)//'_', it
 171             open (unit=junit, file=trim(fname), form='formatted', status='unknown')
 172             write(unit=junit, fmt='(/i5,2x,a,2x,a,2x,4f10.3,i5)') n, &
 173                iv%info(gpseph)%date_char(n), iv%info(gpseph)%id(n),  &
 174                iv%info(gpseph)%lat(1,n)  , iv%info(gpseph)%lon(1,n), &
 175                iv%info(gpseph)%x(1,n)  , iv%info(gpseph)%y(1,n), &
 176                iv%info(gpseph)%levels(n)
 177             write(unit=junit, fmt='(a5,3x,6a14)') 'level','     height   ', &
 178                '    Obs_eph   ','  model_eph   ','  Innov_eph   ',          &
 179                '  error_eph   ',' qc_eph       '
 180             do k = kds,kde
 181                write(unit=junit, fmt='(i5,1x,5f14.3,i10)')  k, &
 182                   mean_h(k)*1000., ob%gpseph(n)%eph(k),        &
 183                   model_eph(k,n), iv%gpseph(n)%eph(k)%inv,     &
 184                   iv%gpseph(n)%eph(k)%error, iv%gpseph(n)%eph(k)%qc
 185             end do
 186          end do
 187          call da_free_unit(junit)
 188          close(junit)
 189       end if  ! end of verify check
 190    end if ! write_iv_gpseph
 191
 192    deallocate (model_eph)
 193
 194    if (trace_use_dull) call da_trace_exit("da_get_innov_vector_gpseph")
 195
 196 end subroutine da_get_innov_vector_gpseph