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