var/da/da_gpspw/da_transform_xtoy_gpspw_adj.inc

   1 subroutine da_transform_xtoy_gpspw_adj(grid, iv, jo_grad_y, jo_grad_x)
   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          ! obs. inc vector (o-b).
  11    type (y_type) , intent(in)    :: jo_grad_y   ! grad_y(jo)
  12    type (x_type) , intent(inout) :: jo_grad_x   ! grad_x(jo)
  13
  14    integer :: n        ! Loop counter.
  15    integer :: i, j     ! Index dimension.
  16    real    :: dx, dxm  !
  17    real    :: dy, dym  !
  18
  19    integer :: k        ! Index dimension.
  20    real    :: dpw,ddpw     ! adjustment pw [mm]
  21    real    :: obs_terr     ! real terrain height at GPS site [m]
  22    real    :: model_terr   ! model terrain height at GPS site[m]
  23    real    :: model_q(kts:kte)   ! model q at GPS site [kg/kg]
  24    real    :: model_z(kts:kte)   ! model z at GPS site [m]
  25    real    :: model_rho(kts:kte) ! model rho at GPS site [kg/m^3]
  26    real    :: model_dq(kts:kte)   ! model dq at GPS site [kg/kg]
  27    real    :: model_drho(kts:kte) ! model drho at GPS site [kg/m^3]
  28
  29    if (trace_use_dull) call da_trace_entry("da_transform_xtoy_gpspw_adj")
  30
  31    do n=iv%info(gpspw)%n1,iv%info(gpspw)%n2
  32       i   = iv%info(gpspw)%i(1,n)
  33       dy  = iv%info(gpspw)%dy(1,n)
  34       dym = iv%info(gpspw)%dym(1,n)
  35       j   = iv%info(gpspw)%j(1,n)
  36       dx  = iv%info(gpspw)%dx(1,n)
  37       dxm = iv%info(gpspw)%dxm(1,n)
  38
  39       !  Initialise the varibles
  40       dpw           = 0.0
  41       model_q(:)    = 0.0
  42       model_dq(:)   = 0.0
  43       model_rho(:)  = 0.0
  44       model_drho(:) = 0.0
  45
  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       dpw = 0.1 *jo_grad_y%gpspw(n)%tpw
  51
  52       jo_grad_x%tpw(i  ,j )  = jo_grad_x%tpw(i  ,j )  + dxm*dym*jo_grad_y%gpspw(n)%tpw
  53       jo_grad_x%tpw(i+1,j )  = jo_grad_x%tpw(i+1,j )  + dym*dx *jo_grad_y%gpspw(n)%tpw
  54       jo_grad_x%tpw(i  ,j+1) = jo_grad_x%tpw(i,j+1)   + dxm *dy*jo_grad_y%gpspw(n)%tpw
  55       jo_grad_x%tpw(i+1,j+1) = jo_grad_x%tpw(i+1,j+1) + dx *dy *jo_grad_y%gpspw(n)%tpw
  56
  57       if (obs_terr <= model_terr) then
  58          model_q(1)   = dym*(dxm*grid%xb%q(i,j,1)   + dx*grid%xb%q(i+1,j,1)) + &
  59             dy *(dxm*grid%xb%q(i,j+1,1) + dx*grid%xb%q(i+1,j+1,1))
  60          model_rho(1) = dym*(dxm*grid%xb%rho(i,j,1)   + dx*grid%xb%rho(i+1,j,1)) + &
  61             dy *(dxm*grid%xb%rho(i,j+1,1) + dx*grid%xb%rho(i+1,j+1,1))
  62
  63          model_dq(1)   =  model_rho(1)*(obs_terr - model_terr)*dpw
  64          model_drho(1) =  model_q(1)  *(obs_terr - model_terr)*dpw
  65
  66          jo_grad_x%q(i,j,1)       = jo_grad_x%q(i,j,1)       + dym*dxm*model_dq(1)
  67          jo_grad_x%q(i+1,j,1)     = jo_grad_x%q(i+1,j,1)     + dym*dx*model_dq(1)
  68          jo_grad_x%q(i,j+1,1)     = jo_grad_x%q(i,j+1,1)     + dy*dxm*model_dq(1)
  69          jo_grad_x%q(i+1,j+1,1)   = jo_grad_x%q(i+1,j+1,1)   + dy*dx*model_dq(1)
  70
  71          jo_grad_x%rho(i,j,1)     = jo_grad_x%rho(i,j,1)     + dym*dxm*model_drho(1)
  72          jo_grad_x%rho(i+1,j,1)   = jo_grad_x%rho(i+1,j,1)   + dym*dx*model_drho(1)
  73          jo_grad_x%rho(i,j+1,1)   = jo_grad_x%rho(i,j+1,1)   + dy*dxm*model_drho(1)
  74          jo_grad_x%rho(i+1,j+1,1) = jo_grad_x%rho(i+1,j+1,1) + dy*dx*model_drho(1)
  75       else
  76          model_z(1) = dym*(dxm*grid%xb%hf(i,j,1)   + dx*grid%xb%hf(i+1,j,1)) + &
  77             dy *(dxm*grid%xb%hf(i,j+1,1) + dx*grid%xb%hf(i+1,j+1,1))
  78
  79          do k = kts, kte
  80            if (model_z(k) >= obs_terr) exit
  81
  82            model_z(k+1) = dym*(dxm*grid%xb%hf(i,j,k+1)   + dx*grid%xb%hf(i+1,j,k+1)) + &
  83               dy *(dxm*grid%xb%hf(i,j+1,k+1) + dx*grid%xb%hf(i+1,j+1,k+1))
  84            model_q(k) = dym*(dxm*grid%xb%q(i,j,k)   + dx*grid%xb%q(i+1,j,k)) + &
  85               dy *(dxm*grid%xb%q(i,j+1,k) + dx*grid%xb%q(i+1,j+1,k))
  86            model_rho(k) = dym*(dxm*grid%xb%rho(i,j,k)   + dx*grid%xb%rho(i+1,j,k)) + &
  87               dy *(dxm*grid%xb%rho(i,j+1,k) + dx*grid%xb%rho(i+1,j+1,k))
  88
  89            ddpw = dpw
  90
  91            if (model_z(k+1) <= obs_terr) then
  92              model_dq  (k) = model_rho(k) *(model_z(k+1) - model_z(k))*ddpw
  93              model_drho(k) = model_q(k)   *(model_z(k+1) - model_z(k))*ddpw
  94            else
  95              model_dq  (k) = model_rho(k) *(obs_terr-model_z(k))*ddpw
  96              model_drho(k) = model_q(k)   *(obs_terr-model_z(k))*ddpw
  97            end if
  98
  99            ! No feedback to x%hf was considered here (Refer to comments in
 100            ! da_transform_xtoy_gpspw.inc).       Y.-R. Guo  07/15/2002
 101            !.........................................................................
 102
 103            jo_grad_x%q(i,j,k)       = jo_grad_x%q(i,j,k)       + dym*dxm*model_dq(k)
 104            jo_grad_x%q(i+1,j,k)     = jo_grad_x%q(i+1,j,k)     + dym*dx*model_dq(k)
 105            jo_grad_x%q(i,j+1,k)     = jo_grad_x%q(i,j+1,k)     + dy*dxm*model_dq(k)
 106            jo_grad_x%q(i+1,j+1,k)   = jo_grad_x%q(i+1,j+1,k)   + dy*dx*model_dq(k)
 107
 108            jo_grad_x%rho(i,j,k)     = jo_grad_x%rho(i,j,k)     + dym*dxm*model_drho(k)
 109            jo_grad_x%rho(i+1,j,k)   = jo_grad_x%rho(i+1,j,k)   + dym*dx*model_drho(k)
 110            jo_grad_x%rho(i,j+1,k)   = jo_grad_x%rho(i,j+1,k)   + dy*dxm*model_drho(k)
 111            jo_grad_x%rho(i+1,j+1,k) = jo_grad_x%rho(i+1,j+1,k) + dy*dx*model_drho(k)
 112         end do
 113       end if
 114    end do
 115
 116    if (trace_use_dull) call da_trace_exit("da_transform_xtoy_gpspw_adj")
 117
 118 end subroutine da_transform_xtoy_gpspw_adj
 119
 120