var/da/da_gpspw/da_transform_xtoy_gpspw.inc

   1 subroutine da_transform_xtoy_gpspw (grid, iv, y)
   2
   3    !----------------------------------------------------------------
   4    ! Purpose: Calculate the difference in the elevation between model surface and GPS site
   5    !----------------------------------------------------------------
   6
   7    implicit none
   8
   9    type (domain),  intent(in)    :: grid
  10    type (iv_type), intent(in)    :: iv       ! Innovation vector (O-B).
  11    type (y_type),  intent(inout) :: y        ! y = h (grid%xa)
  12
  13    integer :: n        ! Loop counter.
  14    integer :: i, j     ! Index dimension.
  15    real    :: dx, dxm  !
  16    real    :: dy, dym  !
  17
  18    integer :: k          ! Index dimension.
  19    real    :: dpw, ddpw  ! adjustment pw [mm]
  20    real    :: obs_terr   ! real terrain height at GPS site [m]
  21    real    :: model_terr ! model terrain height at GPS site[m]
  22    real    :: model_q(kts:kte)    ! model q at GPS site [kg/kg]
  23    real    :: model_z(kts:kte)    ! model z at GPS site [m]
  24    real    :: model_rho(kts:kte)  ! model rho at GPS site [kg/m^3]
  25    real    :: model_dq(kts:kte)   ! model dq at GPS site [kg/kg]
  26    real    :: model_drho(kts:kte) ! model drho at GPS site [kg/m^3]
  27
  28    if (trace_use_dull) call da_trace_entry("da_transform_xtoy_gpspw")
  29
  30    do n=iv%info(gpspw)%n1,iv%info(gpspw)%n2
  31       i   = iv%info(gpspw)%i(1,n)
  32       dy  = iv%info(gpspw)%dy(1,n)
  33       dym = iv%info(gpspw)%dym(1,n)
  34       j   = iv%info(gpspw)%j(1,n)
  35       dx  = iv%info(gpspw)%dx(1,n)
  36       dxm = iv%info(gpspw)%dxm(1,n)
  37
  38       ! Mathematical transformation:
  39
  40       ! kusaka 2004-04-08
  41       ! y % gpspw(n)% tpw = dym* (dxm * grid%xa%tpw(i,j) + &
  42       !    dx  * grid%xa%tpw(i+1,j)) + &
  43       !    dy * (dxm * grid%xa%tpw(i,j+1) + &
  44       !    dx  * grid%xa%tpw(i+1,j+1))
  45       dpw = 0.0
  46       obs_terr   = iv%info(gpspw)%elv(n)
  47       model_terr = dym*(dxm*grid%xb%terr(i,j)   + dx*grid%xb%terr(i+1,j)) + &
  48                    dy *(dxm*grid%xb%terr(i,j+1) + dx*grid%xb%terr(i+1,j+1))
  49
  50       if (obs_terr <= model_terr) then
  51          model_q(1) = dym*(dxm*grid%xb%q(i,j,1)   + dx*grid%xb%q(i+1,j,1)) + &
  52              dy *(dxm*grid%xb%q(i,j+1,1) + dx*grid%xb%q(i+1,j+1,1))
  53          model_rho(1) = dym*(dxm*grid%xb%rho(i,j,1)   + dx*grid%xb%rho(i+1,j,1)) + &
  54             dy *(dxm*grid%xb%rho(i,j+1,1) + dx*grid%xb%rho(i+1,j+1,1))
  55
  56          model_dq(1) = dym*(dxm*grid%xa%q(i,j,1)   + dx*grid%xa%q(i+1,j,1)) + &
  57             dy *(dxm*grid%xa%q(i,j+1,1) + dx*grid%xa%q(i+1,j+1,1))
  58          model_drho(1) = dym*(dxm*grid%xa%rho(i,j,1)   + dx*grid%xa%rho(i+1,j,1)) + &
  59             dy *(dxm*grid%xa%rho(i,j+1,1) + dx*grid%xa%rho(i+1,j+1,1))
  60
  61          dpw = (model_rho(1)*model_dq(1) + model_drho(1)*model_q(1)) &
  62             * (obs_terr - model_terr)
  63       else
  64          model_z(1) = dym*(dxm*grid%xb%hf(i,j,1)   + dx*grid%xb%hf(i+1,j,1)) + &
  65             dy *(dxm*grid%xb%hf(i,j+1,1) + dx*grid%xb%hf(i+1,j+1,1))
  66          do k = kts, kte
  67             if (model_z(k) >= obs_terr) exit
  68
  69             model_z(k+1) = dym*(dxm*grid%xb%hf(i,j,k+1)   + dx*grid%xb%hf(i+1,j,k+1)) + &
  70                dy *(dxm*grid%xb%hf(i,j+1,k+1) + dx*grid%xb%hf(i+1,j+1,k+1))
  71
  72             model_q(k) = dym*(dxm*grid%xb%q(i,j,k)   + dx*grid%xb%q(i+1,j,k)) + &
  73                dy *(dxm*grid%xb%q(i,j+1,k) + dx*grid%xb%q(i+1,j+1,k))
  74             model_rho(k) = dym*(dxm*grid%xb%rho(i,j,k)   + dx*grid%xb%rho(i+1,j,k)) + &
  75                dy *(dxm*grid%xb%rho(i,j+1,k) + dx*grid%xb%rho(i+1,j+1,k))
  76
  77             model_dq(k) = dym*(dxm*grid%xa%q(i,j,k)   + dx*grid%xa%q(i+1,j,k)) + &
  78                dy *(dxm*grid%xa%q(i,j+1,k) + dx*grid%xa%q(i+1,j+1,k))
  79             model_drho(k) = dym*(dxm*grid%xa%rho(i,j,k)   + dx*grid%xa%rho(i+1,j,k)) + &
  80                dy *(dxm*grid%xa%rho(i,j+1,k) + dx*grid%xa%rho(i+1,j+1,k))
  81
  82             ! Here assumed that "model_z" is constant, i.e. grid%xa%hf=0.0. With MM5,
  83             ! this is true, but with WRF, grid%xa%hf may not be zero (?). In the WRF
  84             ! model space (x,y,znu), only the "znu" is constant, but all variables
  85             ! including height could be changed at the "znu" levels. So here is only
  86             ! an approximation for WRF. The following few lines of code is written
  87             ! by Y.-R. Guo 07/16/2004.
  88
  89             if (model_z(k+1) <= obs_terr) then
  90                ddpw = (model_rho(k)*model_dq(k) + model_drho(k)*model_q(k)) * (model_z(k+1)-model_z(k))
  91             else
  92                ddpw = (model_rho(k)*model_dq(k) + model_drho(k)*model_q(k)) * (obs_terr-model_z(k))
  93             end if
  94
  95             dpw = dpw + ddpw
  96          end do
  97       end if
  98       y % gpspw(n)% tpw = dym* (dxm * grid%xa%tpw(i,j) + &
  99                                 dx  * grid%xa%tpw(i+1,j)) + &
 100                          dy * (dxm * grid%xa%tpw(i,j+1) + &
 101                                 dx  * grid%xa%tpw(i+1,j+1)) &
 102                          + 0.1*dpw
 103    end do
 104
 105    if (trace_use_dull) call da_trace_exit("da_transform_xtoy_gpspw")
 106
 107 end subroutine da_transform_xtoy_gpspw
 108
 109