var/da/da_gpspw/da_get_innov_vector_gpsztd.inc

   1 SUBROUTINE da_get_innov_vector_gpsztd ( it, num_qcstat_conv, grid, ob, iv )
   2 !----------------------------------------------------------------
   3 ! Innovation for Ground-based ZTD.
   4 !
   5 ! Because we can ONLY assimilate either GPS PW or GPS ZTD,
   6 ! never assimilate both of them simultaneously, here we
   7 ! used the PW structure for ZTD to avoid declaration of the
   8 ! another structure.
   9 !                                 Y.-R. Guo           05/21/2008
  10 ! History:
  11 !    2017-06: Jamie Bresch
  12 !             (1) reject obs-model height difference larger than Max_StHeight_Diff_ztd
  13 !             (2) properly write out ztd innov info
  14 !             (3) minor clean-up
  15 !----------------------------------------------------------------
  16    IMPLICIT NONE
  17
  18    integer,       intent(in)    :: it      ! External iteration.
  19    type(domain),  intent(in)    :: grid    ! first guess state
  20    type(y_type),  intent(inout) :: ob      ! Observation structure.
  21    type(iv_type), intent(inout) :: iv      ! O-B structure.
  22    integer,       intent(inout) :: num_qcstat_conv(:,:,:,:)
  23
  24    INTEGER                      :: n        ! Loop counter.
  25    INTEGER                      :: i, j     ! Index dimension.
  26    REAL                         :: dx, dxm  ! Interpolation weights.
  27    REAL                         :: dy, dym  ! Interpolation weights.
  28    REAL                         :: mdl_ztd  ! Model ztd at ob location.
  29
  30 !--------------------------------------------------------------------------
  31    INTEGER                :: k            ! Index dimension
  32    REAL                   :: dzd, ddzd    ! adjustment ztd (ref*dz*10**6)
  33    REAL                   :: obs_terr     ! real terrain height at GPS site [m]
  34    REAL                   :: model_terr   ! model terrain height at GPS site[m]
  35    REAL,DIMENSION(kts:kte):: model_ztd    ! model ref at GPS site [kg/kg]
  36    REAL,DIMENSION(kts:kte):: model_z      ! model z at GPS site [m]
  37    INTEGER                :: iunit, ierr
  38    character(len=256)     :: fname
  39
  40
  41    if ( iv%info(gpspw)%nlocal <= 0 ) return
  42
  43    if (trace_use_dull) call da_trace_entry("da_get_innov_vector_gpsztd")
  44
  45    if ( write_iv_gpsztd ) then
  46       if ( num_fgat_time > 1 ) then
  47          write(unit=fname, fmt='(i2.2,a,i2.2,a,i4.4)') &
  48             it,'_inv_gpsztd_t', iv%time, '.', myproc
  49       else
  50          write(unit=fname, fmt='(i2.2,a,i4.4)') it,'_inv_gpsztd.', myproc
  51       end if
  52       call da_get_unit(iunit)
  53       open(unit=iunit,file=trim(fname),form='formatted',iostat=ierr)
  54       write(unit=iunit,fmt='(a4,1x,a4,12a10)')  'ztd:', '   n',   &
  55                        '       lat', '       lon', '   obs_ght',  &
  56                        '   mdl_ght', ' obsh-mdlh', '   obs_ztd',  &
  57                        '   mdl_ztd', '   O-B ztd', '      Dztd',  &
  58                        '  O-B+Dztd', '   Obs_err', '        qc'
  59    end if
  60
  61    do n = iv%info(gpspw)%n1,iv%info(gpspw)%n2
  62
  63       if ( iv % gpspw(n) % tpw % qc == fails_error_max .and. it > 1 ) then
  64          iv % gpspw(n) % tpw % qc = 0
  65       end if
  66
  67       ! Get horizontal interpolation weights:
  68
  69       i   = iv%info(gpspw)%i(1,n)
  70       j   = iv%info(gpspw)%j(1,n)
  71       dx  = iv%info(gpspw)%dx(1,n)
  72       dy  = iv%info(gpspw)%dy(1,n)
  73       dxm = iv%info(gpspw)%dxm(1,n)
  74       dym = iv%info(gpspw)%dym(1,n)
  75
  76       ! xb%ztd (cm) is computed in da_transfer_wrftoxb by calling da_transform_xtoztd
  77       mdl_ztd   = dym*(dxm*grid%xb%ztd(i,j)   + dx*grid%xb%ztd(i+1,j)) + &
  78                   dy *(dxm*grid%xb%ztd(i,j+1) + dx*grid%xb%ztd(i+1,j+1))
  79
  80       if ( .not. pseudo_ztd ) iv % gpspw(n) % tpw % inv = 0.0
  81       if ( ob % gpspw(n) % tpw > missing_r .and. &
  82            iv % gpspw(n) % tpw % qc >= obs_qc_pointer ) then
  83
  84          dzd = 0.0
  85          obs_terr   = iv%info(gpspw)%elv(n)
  86          model_terr = dym*(dxm*grid%xb%terr(i,j)   + dx*grid%xb%terr(i+1,j)) + &
  87                       dy *(dxm*grid%xb%terr(i,j+1) + dx*grid%xb%terr(i+1,j+1))
  88
  89          if ( obs_terr <= model_terr ) then
  90
  91             model_ztd(1) = dym*(dxm*grid%xb%ref(i,j,1)   + dx*grid%xb%ref(i+1,j,1)) + &
  92                            dy *(dxm*grid%xb%ref(i,j+1,1) + dx*grid%xb%ref(i+1,j+1,1))
  93             dzd = model_ztd(1) * ( obs_terr - model_terr )
  94
  95          else
  96
  97             model_z(1) = dym*(dxm*grid%xb%hf(i,j,1)   + dx*grid%xb%hf(i+1,j,1)) + &
  98                          dy *(dxm*grid%xb%hf(i,j+1,1) + dx*grid%xb%hf(i+1,j+1,1))
  99
 100             do k = kts, kte
 101
 102                if (model_z(k) >= obs_terr ) exit
 103
 104                model_z(k+1) = dym*(dxm*grid%xb%hf(i,j,k+1)   + dx*grid%xb%hf(i+1,j,k+1)) + &
 105                               dy *(dxm*grid%xb%hf(i,j+1,k+1) + dx*grid%xb%hf(i+1,j+1,k+1))
 106                model_ztd(k) = dym*(dxm*grid%xb%ref(i,j,k)   + dx*grid%xb%ref(i+1,j,k)) + &
 107                               dy *(dxm*grid%xb%ref(i,j+1,k) + dx*grid%xb%ref(i+1,j+1,k))
 108
 109                if ( model_z(k+1) <= obs_terr ) then
 110                   ddzd = model_ztd(k) * ( model_z(k+1) - model_z(k) )
 111                else
 112                   ddzd = model_ztd(k) * ( obs_terr - model_z(k) )
 113                endif
 114
 115                dzd = dzd + ddzd
 116             end do
 117          end if
 118
 119          if ( pseudo_ztd .and. it == 1 ) then
 120
 121             ! compute ob from input inv for pseudo_ztd
 122             ob % gpspw(n) % tpw = iv % gpspw(n) % tpw % inv + mdl_ztd - 1.e-4 * dzd
 123
 124          else
 125
 126             ! compute inv for ztd
 127             iv % gpspw(n) % tpw % inv = ob % gpspw(n) % tpw - mdl_ztd &
 128                                           + 1.e-4 * dzd
 129 !
 130 ! Overwrite the observation error specification (YRG):
 131 !
 132 !           iv % gpspw(n) % tpw % error = 1.0 + 0.02*(ob%gpspw(n)%tpw-200.)
 133
 134          end if !pseudo_ztd
 135       end if ! valid obs
 136
 137       if ( abs(obs_terr - model_terr) > Max_StHeight_Diff_ztd ) then
 138          iv%gpspw(n)%tpw%qc = -66
 139       end if
 140
 141       if ( write_iv_gpsztd ) then
 142          write(unit=iunit, fmt='(a4,1x,i4,11f10.3,i10)') 'ztd:', n, &
 143               iv%info(gpspw)%lat(1,n), iv%info(gpspw)%lon(1,n), obs_terr, &
 144               model_terr, obs_terr - model_terr, ob%gpspw(n)%tpw,   &
 145               mdl_ztd , ob%gpspw(n)%tpw-mdl_ztd, 1.e-4*dzd,       &
 146               ob%gpspw(n)%tpw-mdl_ztd+1.e-4*dzd, iv%gpspw(n)%tpw%error,&
 147               iv%gpspw(n)%tpw%qc
 148       end if
 149
 150    end do ! n1-n2 loop
 151
 152    ! Perform optional maximum error check:
 153    if ( .not. pseudo_ztd .and. check_max_iv ) &
 154          call da_check_max_iv_gpspw(iv, it, num_qcstat_conv)
 155
 156    if ( write_iv_gpsztd ) then
 157       call da_free_unit(iunit)
 158       close(iunit)
 159    end if
 160
 161    if (trace_use_dull) call da_trace_exit("da_get_innov_vector_gpsztd")
 162
 163 end subroutine da_get_innov_vector_gpsztd
 164