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