var/da/da_gpspw/da_transform_xtoy_gpsztd.inc

   1 subroutine da_transform_xtoy_gpsztd ( grid, iv, y )
   2 !----------------------------------------------------------------
   3 ! Purpose: To calculate the observed ZTD with the height
   4 !          correction in considering the differexbnce between
   5 !          the model terrain height and receiver height.
   6 !
   7 !    The logic is similar to the Gpspw:
   8 !
   9 !      when the receiver below the model surface, a correction,
  10 !      "dzd" should be subtructed from the observed ZTD;
  11 !      when the receiver higher than the model suface, a
  12 !      correction "ztd" should be added to the observed ZTD.
  13 !
  14 !----------------------------------------------------------------
  15
  16    IMPLICIT NONE
  17
  18    type (domain),  intent(in)    :: grid
  19    type (iv_type), intent(in)    :: iv       ! Innovation vector (O-B).
  20    type (y_type),  intent(INOUT) :: y        ! y = h (xa)
  21
  22    INTEGER                      :: n        ! Loop counter.
  23    INTEGER                      :: i, j     ! Index dimension.
  24    REAL                         :: dx, dxm  !
  25    REAL                         :: dy, dym  !
  26 !--
  27    INTEGER                :: k          ! Index dimension.
  28    REAL                   :: dzd, ddzd  ! adjustment pw [mm]
  29    REAL                   :: obs_terr   ! real terrain height at GPS site [m]
  30    REAL                   :: model_terr ! model terrain height at GPS site[m]
  31    REAL,DIMENSION(kts:kte):: model_ztd   ! model ztd at GPS site [m]
  32    REAL,DIMENSION(kts:kte):: model_z     ! model height at GPS site [m]
  33 !--
  34    if (trace_use_dull) call da_trace_entry("da_transform_xtoy_gpsztd")
  35
  36       !bugfix-for-4dvar y % gpspw(:)% tpw = 0.0
  37
  38       do n=iv%info(gpspw)%n1,iv%info(gpspw)%n2
  39
  40          i   = iv%info(gpspw)%i(1,n)
  41          dy  = iv%info(gpspw)%dy(1,n)
  42          dym = iv%info(gpspw)%dym(1,n)
  43          j   = iv%info(gpspw)%j(1,n)
  44          dx  = iv%info(gpspw)%dx(1,n)
  45          dxm = iv%info(gpspw)%dxm(1,n)
  46
  47 ! Mathematical transformation:
  48
  49          dzd = 0.0
  50          obs_terr   = iv%info(gpspw)%elv(n)
  51          model_terr = dym*(dxm*grid%xb%terr(i,j)   + dx*grid%xb%terr(i+1,j)) + &
  52                       dy *(dxm*grid%xb%terr(i,j+1) + dx*grid%xb%terr(i+1,j+1))
  53
  54          if ( obs_terr <= model_terr ) then
  55
  56             model_ztd(1) = dym*(dxm*grid%xa%ref(i,j,1)   + dx*grid%xa%ref(i+1,j,1)) + &
  57                            dy *(dxm*grid%xa%ref(i,j+1,1) + dx*grid%xa%ref(i+1,j+1,1))
  58
  59             dzd =  model_ztd(1) * ( obs_terr - model_terr )
  60
  61          else
  62
  63             model_z(1) = dym*(dxm*grid%xb%hf(i,j,1)   + dx*grid%xb%hf(i+1,j,1)) + &
  64                          dy *(dxm*grid%xb%hf(i,j+1,1) + dx*grid%xb%hf(i+1,j+1,1))
  65
  66             do k = kts, kte
  67
  68               if (model_z(k) >= obs_terr ) exit
  69
  70               model_z(k+1) = dym*(dxm*grid%xb%hf(i,j,k+1)   + dx*grid%xb%hf(i+1,j,k+1)) + &
  71                              dy *(dxm*grid%xb%hf(i,j+1,k+1) + dx*grid%xb%hf(i+1,j+1,k+1))
  72
  73               model_ztd(k) = dym*(dxm*grid%xa%ref(i,j,k)   + dx*grid%xa%ref(i+1,j,k)) + &
  74                              dy *(dxm*grid%xa%ref(i,j+1,k) + dx*grid%xa%ref(i+1,j+1,k))
  75 !
  76 ! Here assumed that "model_z" is constant, i.e. grid%xa%hf=0.0. With MM5,
  77 ! this is true, but with WRF, grid%xa%hf may not be zero (?). In the WRF
  78 ! model space (x,y,znu), only the "znu" is constant, but all variables
  79 ! including height could be changed at the "znu" levels. So here is only
  80 ! an approximation for WRF. The following few lines of code is written
  81 ! by Y.-R. Guo 07/16/2004.
  82 !
  83               if ( model_z(k+1) <= obs_terr ) then
  84                  ddzd = model_ztd(k) * ( model_z(k+1) - model_z(k) )
  85               else
  86                  ddzd = model_ztd(k) * ( obs_terr -  model_z(k) )
  87               endif
  88
  89               dzd = dzd + ddzd
  90
  91             end do
  92          end if
  93
  94          y % gpspw(n)% tpw = dym* ( dxm * grid%xa%ztd(i,j) + &
  95                                     dx  * grid%xa%ztd(i+1,j) ) + &
  96                              dy * ( dxm * grid%xa%ztd(i,j+1) + &
  97                                     dx  * grid%xa%ztd(i+1,j+1) ) &
  98                              + 1.e-4 * dzd
  99
 100       end do
 101
 102    if (trace_use_dull) call da_trace_exit("da_transform_xtoy_gpsztd")
 103
 104 end subroutine da_transform_xtoy_gpsztd
 105