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