da_get_innov_vector_radar.inc
References to this file elsewhere.
1 subroutine da_get_innov_vector_radar (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, k ! Index dimension.
16 real :: dx, dxm ! Interpolation weights.
17 real :: dy, dym ! Interpolation weights.
18
19 real, allocatable :: model_p(:,:)
20 real, allocatable :: model_u(:,:)
21 real, allocatable :: model_v(:,:)
22 real, allocatable :: model_w(:,:)
23 real, allocatable :: model_rho(:,:)
24 real, allocatable :: model_qrn(:,:)
25 real, allocatable :: model_rv(:,:)
26 real, allocatable :: model_rf(:,:)
27 real, allocatable :: model_ps(:)
28
29 real :: v_h(kms:kme) ! Model value h at ob hor. location.
30
31 real :: xr,yr,zr
32 integer :: irv, irvf
33 integer :: irf, irff
34
35 if (trace_use) call da_trace_entry("da_get_innov_vector_radar")
36
37 irv = 0; irvf = 0; irf = 0; irff = 0
38
39 allocate (model_p(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
40 allocate (model_u(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
41 allocate (model_v(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
42 allocate (model_w(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
43 allocate (model_rho(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
44 allocate (model_qrn(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
45 allocate (model_rv(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
46 allocate (model_rf(iv%info(radar)%max_lev,iv%info(radar)%n1:iv%info(radar)%n2))
47 allocate (model_ps(iv%info(radar)%n1:iv%info(radar)%n2))
48
49 model_p(:,:) = 0
50 model_u(:,:) = 0
51 model_v(:,:) = 0
52 model_w(:,:) = 0
53 model_rho(:,:) = 0
54 model_qrn(:,:) = 0
55 model_rv(:,:) = 0
56 model_rf(:,:) = 0
57 model_ps(:) = 0
58
59 do n=iv%info(radar)%n1,iv%info(radar)%n2
60 if (iv%info(radar)%levels(n) < 1) cycle
61
62 ! [1.0] Get horizontal interpolation weights:
63
64 i = iv%info(radar)%i(1,n)
65 j = iv%info(radar)%j(1,n)
66 dx = iv%info(radar)%dx(1,n)
67 dy = iv%info(radar)%dy(1,n)
68 dxm = iv%info(radar)%dxm(1,n)
69 dym = iv%info(radar)%dym(1,n)
70
71 do k=kts,kte
72 v_h(k) = dym*(dxm*grid%xb%h(i,j,k)+dx*grid%xb%h(i+1,j,k)) + dy*(dxm*grid%xb%h(i,j+1,k)+dx*grid%xb%h(i+1,j+1,k))
73 end do
74
75 do k=1, iv%info(radar)%levels(n)
76 call da_to_zk(iv%radar(n)%height(k), v_h, v_interp_h, iv%info(radar)%zk(k,n))
77
78 if (iv%info(radar)%zk(k,n) < 1.0) then
79 iv%radar(n)%height_qc(k) = below_model_surface
80 else if (iv%info(radar)%zk(k,n) > mkz) then
81 iv%radar(n)%height_qc(k) = above_model_lid
82 end if
83 end do
84 end do
85
86 call da_convert_zk (iv%info(radar))
87
88 ! [2.0] Interpolate horizontally to ob points:
89
90 call da_interp_lin_3d (grid%xb % p, iv%info(radar), model_p)
91 call da_interp_lin_3d (grid%xb % u, iv%info(radar), model_u)
92 call da_interp_lin_3d (grid%xb % v, iv%info(radar), model_v)
93 call da_interp_lin_3d (grid%xb % wh, iv%info(radar), model_w)
94 call da_interp_lin_3d (grid%xb % rho, iv%info(radar), model_rho)
95 call da_interp_lin_3d (grid%xb % qrn, iv%info(radar), model_qrn)
96
97 ! Test 5.0E-5 as critical value. It can not be smaller.
98
99 do n=iv%info(radar)%n1,iv%info(radar)%n2
100 do k=1,iv%info(radar)%levels(n)
101 model_qrn(k,n)=amax1(5.0E-5,model_qrn(k,n))
102 end do
103
104 i = iv%info(radar)%i(1,n)
105 j = iv%info(radar)%j(1,n)
106 dx = iv%info(radar)%dx(1,n)
107 dy = iv%info(radar)%dy(1,n)
108 dxm = iv%info(radar)%dxm(1,n)
109 dym = iv%info(radar)%dym(1,n)
110
111
112 model_ps(n) = dxm *(dym * grid%xb % psac(i, j) + dy * grid%xb%psac(i+1, j)) + &
113 dx *(dym * grid%xb % psac(i,j+1) + dy * grid%xb%psac(i+1,j+1)) + &
114 grid%xb % ptop
115
116 iv%radar(n)%model_p(:) = model_p(1:iv%info(radar)%levels(n),n)
117 iv%radar(n)%model_rho(:) = model_rho(1:iv%info(radar)%levels(n),n)
118 iv%radar(n)%model_qrn(:) = model_qrn(1:iv%info(radar)%levels(n),n)
119 iv%radar(n)%model_ps = model_ps(n)
120
121 ! [3.0] Calculate rv, rf at OBS location and initialise components of &
122 ! innovation vector:
123
124 if (fg_format == fg_format_wrf) then
125 call da_llxy_wrf(map_info, &
126 iv%radar(n)%stn_loc%lat, iv%radar(n)%stn_loc%lon, &
127 iv%radar(n)%stn_loc%x, iv%radar(n)%stn_loc%y)
128 else
129 call da_llxy_default( iv%radar(n)%stn_loc%lat, iv%radar(n)%stn_loc%lon, &
130 iv%radar(n)%stn_loc%x, iv%radar(n)%stn_loc%y)
131 end if
132
133 xr = grid%xb%ds *(iv%info(radar)%x(1,n) - iv%radar(n)%stn_loc%x)
134 yr = grid%xb%ds *(iv%info(radar)%y(1,n) - iv%radar(n)%stn_loc%y)
135
136 do k=1, iv%info(radar)%levels(n)
137 iv % radar(n) % rv(k) % inv = 0.0
138 iv % radar(n) % rf(k) % inv = 0.0
139
140 if (iv % radar(n) % height_qc(k) /= below_model_surface .and. &
141 iv % radar(n) % height_qc(k) /= above_model_lid) then
142
143 if (use_radar_rv) then
144 if (abs(iv % radar(n) % rv(k) % qc - missing_data) > 1) then
145 if (abs(ob % radar(n) % rv(k) - missing_r) > 1.0 .AND. &
146 iv % radar(n) % rv(k) % qc >= obs_qc_pointer) then
147
148 zr=iv%radar(n)%height(k) - iv%radar(n)%stn_loc%elv
149
150 call da_radial_velocity(model_rv(k,n), model_p(k,n), &
151 model_u(k,n), model_v(k,n), model_w(k,n), &
152 model_qrn(k,n), model_ps(n), xr, yr, zr)
153
154 iv % radar(n) % rv(k) % inv = ob % radar(n) % rv(k) - model_rv(k,n)
155 end if
156 end if
157 end if
158
159 if (use_radar_rf) then
160 if (abs(iv % radar(n) % rf(k) % qc - missing_data) > 1) then
161 if (abs(ob % radar(n) % rf(k) - missing_r) > 1.0 .AND. &
162 iv % radar(n) % rf(k) % qc >= obs_qc_pointer) then
163 model_rf(k,n) = leh1 + leh2 * alog10(model_rho(k,n) * model_qrn(k,n) * 1.0e+3)
164 iv % radar(n) % rf(k) % inv = ob % radar(n) % rf(k) - model_rf(k,n)
165 end if
166 end if
167 end if
168 end if
169 end do
170 end do
171
172 !------------------------------------------------------------------------
173 ! [4.0] Perform optional maximum error check:
174 !------------------------------------------------------------------------
175
176 if (check_max_iv) then
177 call da_check_max_iv_radar(iv, it, irv, irf, irvf, irff)
178 end if
179
180 if (rootproc .and. check_max_iv_print) then
181 write(unit = check_max_iv_unit, fmt ='(/,A,i5,A)') &
182 'For outer iteration ', it, ', Total Rejections for radar follows:'
183
184 if (use_radar_rv) then
185 write( unit = check_max_iv_unit, fmt = '(/,2(A,I6))') &
186 'Number of failed rv observations: ',irvf, ' on ',irv
187 end if
188
189 if (use_radar_rf) then
190 write( unit = check_max_iv_unit, fmt = '(/,2(A,I6))') &
191 'Number of failed rf observations: ',irff, ' on ',irf
192 end if
193 end if
194
195 deallocate (model_p)
196 deallocate (model_u)
197 deallocate (model_v)
198 deallocate (model_w)
199 deallocate (model_qrn)
200 deallocate (model_rv)
201 deallocate (model_rf)
202 deallocate (model_ps)
203
204 if (trace_use) call da_trace_exit("da_get_innov_vector_radar")
205
206 end subroutine da_get_innov_vector_radar
207
208