da_get_innov_vector_gpspw.inc
References to this file elsewhere.
1 subroutine da_get_innov_vector_gpspw( it, xb, xp, ob, iv)
2
3 !----------------------------------------------------------------
4 ! Purpose: TBD
5 !----------------------------------------------------------------
6
7 implicit none
8
9 integer, intent(in) :: it ! External iteration.
10 type(xb_type), intent(in) :: xb ! first guess state.
11 type(xpose_type), intent(in) :: xp ! Domain decomposition vars.
12 type(y_type), intent(inout) :: ob ! Observation structure.
13 type(ob_type), intent(inout) :: iv ! O-B structure.
14
15 integer :: n ! Loop counter.
16 integer :: i, j ! Index dimension.
17 real :: dx, dxm ! Interpolation weights.
18 real :: dy, dym ! Interpolation weights.
19 real :: model_tpw! Model value u at oblocation.
20 integer :: ittpw,ittpwf
21
22 !--------------------------------------------------------------------------
23 integer :: k ! Index dimension
24 real :: dpw, ddpw ! adjustment pw [mm]
25 real :: obs_terr ! real terrain height at GPS site [m]
26 real :: model_terr ! model terrain height at GPS site[m]
27 real,dimension(xp%kts:xp%kte):: model_q ! model q at GPS site [kg/kg]
28 real,dimension(xp%kts:xp%kte):: model_z ! model z at GPS site [m]
29 real,dimension(xp%kts:xp%kte):: model_rho ! model rho at GPS site [kg/m^3]
30
31 ! unit_gps = myproc + 140
32 !---------------------------------------------------------------------------
33
34 !
35 ! GPS TPW Pseudo OBS test:
36 !
37 if (pseudo_var(1:3) == 'tpw' .and. num_pseudo > 0) then
38 ! Deallocate:
39 if (iv%num_gpspw > 0) then
40 write(unit=stdout, fmt='(a,i4)') &
41 'iv%num_gpspw =', iv%num_gpspw
42 deallocate(iv % gpspw)
43 deallocate(ob % gpspw)
44 end if
45
46 Use_GpspwObs = .true.
47
48 ! Allocate:
49 iv % num_gpspw = num_pseudo
50 iv % ob_numb(1)%gpspw = num_pseudo
51 iv % num_pseudo = 0
52
53 allocate(iv % gpspw(1:num_pseudo))
54 allocate(ob % gpspw(1:num_pseudo))
55
56 write(unit=stdout,fmt='(a,i2)') &
57 '==> GPS TPW pseudo OBS test: num_pseudo=',num_pseudo
58
59 iv % gpspw(1) % loc % x = pseudo_x
60 iv % gpspw(1) % loc % y = pseudo_y
61 iv % gpspw(1) % loc % i = int(pseudo_x)
62 iv % gpspw(1) % loc % j = int(pseudo_y)
63 iv % gpspw(1) % loc % dx = pseudo_x-real(iv % gpspw(1) % loc % i)
64 iv % gpspw(1) % loc % dy = pseudo_y-real(iv % gpspw(1) % loc % j)
65 iv % gpspw(1) % loc % dxm = 1.0 - iv % gpspw(1) % loc % dx
66 iv % gpspw(1) % loc % dym = 1.0 - iv % gpspw(1) % loc % dy
67
68 iv % gpspw(1) % tpw % inv = pseudo_val
69 iv % gpspw(1) % tpw % qc = 0
70 iv % gpspw(1) % tpw % error = pseudo_err
71
72 ! To consider the site elevation effect, set the model terrain height
73 ! to elevation for pseudo OBS:
74
75 i = iv%gpspw(1)%loc%i
76 j = iv%gpspw(1)%loc%j
77 dx = iv%gpspw(1)%loc%dx
78 dy = iv%gpspw(1)%loc%dy
79 dxm = iv%gpspw(1)%loc%dxm
80 dym = iv%gpspw(1)%loc%dym
81
82 iv%gpspw(1)%info%elv = dym*(dxm*xb%terr(i,j) + dx*xb%terr(i+1,j)) + &
83 dy *(dxm*xb%terr(i,j+1) + dx*xb%terr(i+1,j+1))
84
85 ! Set halo:
86 if ((iv%gpspw(1)%loc%i < xp%its-1) .or.(iv%gpspw(1)%loc%i > xp%ite) .or. &
87 (iv%gpspw(1)%loc%j < xp%jts-1) .or.(iv%gpspw(1)%loc%j > xp%jte)) then
88 call da_error(__FILE__,__LINE__,(/"Should never have obs outside halo by now"/))
89 ! iv%gpspw(1)%loc%proc_domain_with_halo = .false.
90 iv%gpspw(1)%loc%proc_domain = .false.
91 else
92 ! iv%gpspw(1)%loc%proc_domain_with_halo = .true.
93 iv%gpspw(1)%loc%proc_domain = .false.
94
95 if (iv%gpspw(1)%loc%i >= xp%its .and. iv%gpspw(1)%loc%i <= xp%ite .and. &
96 iv%gpspw(1)%loc%j >= xp%jts .and. iv%gpspw(1)%loc%j <= xp%jte) then
97 iv%gpspw(1)%loc%proc_domain = .true.
98 end if
99 end if
100
101 write(unit=stdout,fmt='(a4,2f15.5)') pseudo_var, pseudo_val, pseudo_err
102 write(unit=stdout,fmt='(3f15.5)') pseudo_x, pseudo_y, pseudo_z
103 write(unit=stdout,fmt='(a,f8.2)') 'iv%gpspw: elv=',iv%gpspw(1)%info%elv
104 end if
105
106 if (iv % num_gpspw > 0) then
107
108 ittpw = 0
109 ittpwf = 0
110
111 ! write(unit=unit_gps,fmt='(3x,a3,10a10)') ' n ',' lat ', &
112 ! ' lon ', ' obs ght ', ' mdl ght ', &
113 ! ' obsh-mdlh', ' obs pw ', ' model pw', &
114 ! ' O-B pw ', ' Dpw ', ' O-B+Dpw '
115
116 do n=iv%ob_numb(iv%current_ob_time-1)%gpspw + 1, iv%ob_numb(iv%current_ob_time)%gpspw
117
118 ! [1.1] Get horizontal interpolation weights:
119
120 i = iv%gpspw(n)%loc%i
121 j = iv%gpspw(n)%loc%j
122 dx = iv%gpspw(n)%loc%dx
123 dy = iv%gpspw(n)%loc%dy
124 dxm = iv%gpspw(n)%loc%dxm
125 dym = iv%gpspw(n)%loc%dym
126
127 model_tpw = dym*(dxm*xb%tpw(i,j) + dx*xb%tpw(i+1,j)) + &
128 dy *(dxm*xb%tpw(i,j+1) + dx*xb%tpw(i+1,j+1))
129
130 if (pseudo_var(1:3) == 'tpw' .and. num_pseudo > 0) then
131 ! To compute the 'ob':
132 ob % gpspw(n) % tpw = iv % gpspw(n) % tpw % inv + model_tpw
133 else
134 ! To compute the 'inv':
135 iv % gpspw(n) % tpw % inv = 0.0
136 if (ob % gpspw(n) % tpw > missing_r .AND. &
137 iv % gpspw(n) % tpw % qc >= obs_qc_pointer) then
138 dpw = 0.0
139 obs_terr = iv%gpspw(n)%info%elv
140 model_terr = dym*(dxm*xb%terr(i,j) + dx*xb%terr(i+1,j)) + &
141 dy *(dxm*xb%terr(i,j+1) + dx*xb%terr(i+1,j+1))
142
143 if (obs_terr <= model_terr) then
144 model_q(1) = dym*(dxm*xb%q(i,j,1) + dx*xb%q(i+1,j,1)) + &
145 dy *(dxm*xb%q(i,j+1,1) + dx*xb%q(i+1,j+1,1))
146 model_rho(1) = dym*(dxm*xb%rho(i,j,1) + dx*xb%rho(i+1,j,1)) + &
147 dy *(dxm*xb%rho(i,j+1,1) + dx*xb%rho(i+1,j+1,1))
148 dpw = model_rho(1) * model_q(1) *( obs_terr - model_terr)
149 else
150 model_z(1) = dym*(dxm*xb%hf(i,j,1) + dx*xb%hf(i+1,j,1)) + &
151 dy *(dxm*xb%hf(i,j+1,1) + dx*xb%hf(i+1,j+1,1))
152 do k = xp%kts, xp%kte
153 if (model_z(k) >= obs_terr) exit
154
155 model_z(k+1) = dym*(dxm*xb%hf(i,j,k+1) + dx*xb%hf(i+1,j,k+1)) + &
156 dy *(dxm*xb%hf(i,j+1,k+1) + dx*xb%hf(i+1,j+1,k+1))
157 model_q(k) = dym*(dxm*xb%q(i,j,k) + dx*xb%q(i+1,j,k)) + &
158 dy *(dxm*xb%q(i,j+1,k) + dx*xb%q(i+1,j+1,k))
159 model_rho(k) = dym*(dxm*xb%rho(i,j,k) + dx*xb%rho(i+1,j,k)) + &
160 dy *(dxm*xb%rho(i,j+1,k) + dx*xb%rho(i+1,j+1,k))
161
162 if (model_z(k+1) <= obs_terr) then
163 ddpw = model_rho(k) * model_q(k) *( model_z(k+1) - model_z(k))
164 else
165 ddpw = model_rho(k) * model_q(k) *( obs_terr - model_z(k))
166 end if
167
168 dpw = dpw + ddpw
169 end do
170 end if
171
172 iv % gpspw(n) % tpw % inv = ob % gpspw(n) % tpw - model_tpw &
173 + 0.1*dpw
174 end if
175 end if
176
177 !------------------------------------------------------------------------
178 ! [5.0] Perform optional maximum error check:
179 !------------------------------------------------------------------------
180
181 if ( check_max_iv) then
182 call da_check_max_iv_gpspw(it, iv % gpspw(n), ittpw,ittpwf)
183 end if
184
185 ! write(unit=unit_gps, fmt='(i4,10f10.3)') n, &
186 ! iv%gpspw(n)%info%lat, iv%gpspw(n)%info%lon, obs_terr, &
187 ! model_terr, obs_terr - model_terr, ob%gpspw(n)%tpw, &
188 ! model_tpw , ob%gpspw(n)%tpw-model_tpw, 0.1*dpw, &
189 ! ob%gpspw(n)%tpw-model_tpw+0.1*dpw
190 end do
191
192 if (rootproc .and. check_max_iv_print) then
193 write(unit = check_max_iv_unit, fmt ='(A,i5,A)')'For outer iteration ',it, &
194 ', Total Rejections for Gpspw follows:'
195 write(unit = check_max_iv_unit, fmt = '(/,10(2(A,I6),/))') &
196 'Number of failed TPW-observations: ',ittpwf, ' on ',ittpw, &
197 'Finally all Gpspw rejections ',ittpwf,' on ',ittpw
198 end if
199 end if
200
201 end subroutine da_get_innov_vector_gpspw
202
203