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