da_get_innov_vector_rttov.inc
References to this file elsewhere.
1 subroutine da_get_innov_vector_rttov (it, xb, xp, ob, iv)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Calculate innovation vector for radiance data.
5 !
6 ! METHOD: d = y - H(x)
7 ! 1. interpolate xb to obs location
8 ! 2. call forward RTM to get simulated bright temperature
9 ! 3. obs BT - simulated BT
10 !---------------------------------------------------------------------------
11
12 implicit none
13
14 !#include "rttov_q2v.interface"
15
16 integer, intent(in) :: it ! External iteration.
17 type (xb_type), intent(in) :: xb ! first guess state.
18 type (xpose_type), intent(in) :: xp ! Domain decomposition vars.
19 type (y_type), intent(inout) :: ob ! Observation structure.
20 type (ob_type), intent(inout) :: iv ! O-B structure.
21
22 integer :: n,i ! Loop counter.
23 integer :: k ! Index dimension.
24 integer :: nlevels ! Number of obs levels.
25
26 real :: seap, icep, lndp, snop
27 real, allocatable :: v_p(:,:)
28
29 integer :: inst, nchan
30 real, allocatable :: polar(:)
31 real, allocatable :: pres(:)
32 logical :: amsua, amsub, ssmi
33 ! logical :: mhs
34
35 integer :: kn,knchpf,sensor_type,n1,n2,n1n2
36 integer, allocatable :: kprof(:), kchan(:),kochan(:),knchan(:),indexn(:,:)
37 real , allocatable :: pemsv(:), pemsh(:), pems5(:)
38 real , allocatable :: zasat(:), zlsat(:)
39
40 real,allocatable :: temp(:)
41
42 type(con_vars_type), allocatable :: con_vars(:)
43 type(aux_vars_type), allocatable :: aux_vars(:)
44
45 ! variables for computing clwp
46 real, allocatable :: dpf(:,:), clw(:,:), pf(:,:)
47
48 if ( iv%num_inst < 1 ) return
49
50 #if !defined(RTTOV)
51 call da_error(__FILE__,__LINE__, &
52 (/"Must compile with $RTTOV option for radiances"/))
53 #else
54
55 if (trace_use) call da_trace_entry("da_get_innov_vector_rttov")
56
57 !------------------------------------------------------
58 ! [1.0] calculate the background bright temperature
59 !-------------------------------------------------------
60
61 do inst = 1, iv%num_inst ! loop for sensor
62 if ( iv%instid(inst)%num_rad < 1 ) cycle
63 nlevels = iv%instid(inst)%nlevels
64 nchan = iv%instid(inst)%nchan
65
66 amsua = trim( inst_name(rtminit_sensor(inst))) == 'amsua'
67 amsub = trim( inst_name(rtminit_sensor(inst))) == 'amsub'
68 ! mhs = trim( inst_name(rtminit_sensor(inst))) == 'mhs'
69 ssmi = trim( inst_name(rtminit_sensor(inst))) == 'ssmi'
70 sensor_type = coefs(inst)%id_sensor - 1 ! 0:IR 1:MW
71
72 if (iv%ob_numb(iv%current_ob_time)%radiance(inst) <= &
73 iv%ob_numb(iv%current_ob_time-1)%radiance(inst) ) cycle
74 n1 = iv%ob_numb(iv%current_ob_time-1)%radiance(inst)+1
75 n2 = iv%ob_numb(iv%current_ob_time)%radiance(inst)
76
77 n1n2=n2-n1+1
78
79 allocate (pres(1:nlevels))
80 allocate (con_vars(n1:n2))
81 allocate (aux_vars(n1:n2))
82 allocate (zasat(n1:n2))
83 allocate (zlsat(n1:n2))
84 allocate (polar(nchan))
85
86 knchpf = n1n2*nchan
87 ! knchpf = 1*nchan
88
89 allocate (knchan (n1:n2))
90 ! allocate (knchan (1))
91 allocate (kprof (knchpf))
92 allocate (kchan (knchpf))
93 allocate (kochan (knchpf))
94 allocate (indexn (1:nchan,n1:n2))
95 allocate (pemsv (knchpf))
96 allocate (pemsh (knchpf))
97 allocate (pems5 (knchpf))
98
99 pres(1:nlevels) = coefs(inst) % ref_prfl_p(1:nlevels)
100
101 allocate(v_p(xp%kms:xp%kme,n1:n2))
102 v_p(:,:)=0.0
103 allocate(clw(xp%kms:xp%kme,n1:n2))
104 allocate(dpf(xp%kms:xp%kme,n1:n2))
105 allocate(pf(xp%kms:xp%kme+1,n1:n2))
106
107 ! horizontal interpolate xb pressure to ob position for every xb levels
108 do n=n1,n2
109 do k=xp%kts,xp%kte
110 ! convert to mb
111 v_p(k,n) = 0.01*(iv%instid(inst)%loc_dym(n)*( &
112 iv%instid(inst)%loc_dxm(n)*xb%p(iv%instid(inst)%loc_i(n), iv%instid(inst)%loc_j(n),k) + &
113 iv%instid(inst)%loc_dx(n) *xb%p(iv%instid(inst)%loc_i(n)+1,iv%instid(inst)%loc_j(n),k)) + &
114 iv%instid(inst)%loc_dy(n) *( &
115 iv%instid(inst)%loc_dxm(n)*xb%p(iv%instid(inst)%loc_i(n), iv%instid(inst)%loc_j(n)+1,k) + &
116 iv%instid(inst)%loc_dx(n) *xb%p(iv%instid(inst)%loc_i(n)+1,iv%instid(inst)%loc_j(n)+1,k)))
117 end do
118 end do
119
120 call da_to_zk_new(pres, v_p(:,n1:n2), xp, v_interp_p, &
121 iv%instid(inst)%zk(:,n1:n2), n1n2,nlevels )
122
123 call da_togrid_new (iv%instid(inst)%zk(:,n1:n2), xp%kts, xp%kte, &
124 iv%instid(inst)%loc_k(:,n1:n2), &
125 iv%instid(inst)%loc_dz(:,n1:n2), &
126 iv%instid(inst)%loc_dzm(:,n1:n2))
127 ! [1.2] Interpolate horizontally to ob:
128
129 call da_interp_lin_3d_new( xb%t, xp, &
130 iv%instid(inst)%loc_i(n1:n2), &
131 iv%instid(inst)%loc_j(n1:n2), &
132 iv%instid(inst)%loc_k(:,n1:n2), &
133 iv%instid(inst)%loc_dx(n1:n2), &
134 iv%instid(inst)%loc_dy(n1:n2), &
135 iv%instid(inst)%loc_dz(:,n1:n2), &
136 iv%instid(inst)%loc_dxm(n1:n2), &
137 iv%instid(inst)%loc_dym(n1:n2), &
138 iv%instid(inst)%loc_dzm(:,n1:n2), &
139 iv%instid(inst)%t(:,n1:n2), nlevels, n1n2 )
140
141 call da_interp_lin_3d_new( xb%q, xp, &
142 iv%instid(inst)%loc_i(n1:n2), &
143 iv%instid(inst)%loc_j(n1:n2), &
144 iv%instid(inst)%loc_k(:,n1:n2), &
145 iv%instid(inst)%loc_dx(n1:n2), &
146 iv%instid(inst)%loc_dy(n1:n2), &
147 iv%instid(inst)%loc_dz(:,n1:n2), &
148 iv%instid(inst)%loc_dxm(n1:n2), &
149 iv%instid(inst)%loc_dym(n1:n2), &
150 iv%instid(inst)%loc_dzm(:,n1:n2), &
151 iv%instid(inst)%mr(:,n1:n2), nlevels, n1n2 )
152
153 do n= n1,n2
154 do k=1, nlevels
155 if (iv%instid(inst)%zk(k,n) <= 0.0) then
156 iv%instid(inst)%t(k,n) = coefs(inst) % ref_prfl_t(k,gas_id_watervapour) ! outside model level
157 iv%instid(inst)%mr(k,n) = coefs(inst) % ref_prfl_mr(k,gas_id_watervapour)
158 else
159 iv%instid(inst)%mr(k,n) = iv%instid(inst)%mr(k,n) * q2ppmv
160 ! call rttov_q2v (&
161 ! gas_unit_specconc, &! in
162 ! iv%instid(inst)%q(k,n), &! in
163 ! gas_id_watervapour, &! in
164 ! iv%instid(inst)%q(k,n), &! in
165 ! iv%instid(inst)%mr(k,n) ) ! inout
166 end if
167 if (pres(k) < 100.0) iv%instid(inst)%mr(k,n) = coefs(inst) % ref_prfl_mr(k,gas_id_watervapour)
168 end do
169
170 ! determine surface type of obs location
171 !-----------------------------------------
172 call da_detsurtyp( xb%snow, xb%xice, xb%landmask, &
173 xb%ivgtyp, xb%isltyp, &
174 xp%ims, xp%ime, xp%jms, xp%jme, &
175 iv%instid(inst)%loc_i(n), iv%instid(inst)%loc_j(n), &
176 iv%instid(inst)%loc_dx(n), iv%instid(inst)%loc_dy(n), &
177 iv%instid(inst)%loc_dxm(n), iv%instid(inst)%loc_dym(n), &
178 iv%instid(inst)%isflg(n),iv%instid(inst)%vegtyp(n), iv%instid(inst)%soiltyp(n), &
179 seap, icep, lndp, snop )
180 end do
181
182 call da_interp_lin_2d_new (xb % u10, xp%ims, xp%ime, xp%jms, xp%jme, &
183 iv%instid(inst)%loc_i(n1:n2), &
184 iv%instid(inst)%loc_j(n1:n2), &
185 iv%instid(inst)%loc_dx(n1:n2), &
186 iv%instid(inst)%loc_dy(n1:n2), &
187 iv%instid(inst)%loc_dxm(n1:n2), &
188 iv%instid(inst)%loc_dym(n1:n2), &
189 iv%instid(inst)%u10(n1:n2),n1n2)
190
191 call da_interp_lin_2d_new (xb % v10, xp%ims, xp%ime, xp%jms, xp%jme, &
192 iv%instid(inst)%loc_i(n1:n2), &
193 iv%instid(inst)%loc_j(n1:n2), &
194 iv%instid(inst)%loc_dx(n1:n2), &
195 iv%instid(inst)%loc_dy(n1:n2), &
196 iv%instid(inst)%loc_dxm(n1:n2), &
197 iv%instid(inst)%loc_dym(n1:n2), &
198 iv%instid(inst)%v10(n1:n2),n1n2)
199
200 call da_interp_lin_2d_new (xb % t2, xp%ims, xp%ime, xp%jms, xp%jme, &
201 iv%instid(inst)%loc_i(n1:n2), &
202 iv%instid(inst)%loc_j(n1:n2), &
203 iv%instid(inst)%loc_dx(n1:n2), &
204 iv%instid(inst)%loc_dy(n1:n2), &
205 iv%instid(inst)%loc_dxm(n1:n2), &
206 iv%instid(inst)%loc_dym(n1:n2), &
207 iv%instid(inst)%t2m(n1:n2) ,n1n2)
208
209 call da_interp_lin_2d_new (xb % q2, xp%ims, xp%ime, xp%jms, xp%jme, &
210 iv%instid(inst)%loc_i(n1:n2), &
211 iv%instid(inst)%loc_j(n1:n2), &
212 iv%instid(inst)%loc_dx(n1:n2), &
213 iv%instid(inst)%loc_dy(n1:n2), &
214 iv%instid(inst)%loc_dxm(n1:n2), &
215 iv%instid(inst)%loc_dym(n1:n2), &
216 iv%instid(inst)%q2m(n1:n2),n1n2)
217
218 call da_interp_lin_2d_new(xb % psfc, xp%ims, xp%ime, xp%jms, xp%jme, &
219 iv%instid(inst)%loc_i(n1:n2), &
220 iv%instid(inst)%loc_j(n1:n2), &
221 iv%instid(inst)%loc_dx(n1:n2), &
222 iv%instid(inst)%loc_dy(n1:n2), &
223 iv%instid(inst)%loc_dxm(n1:n2), &
224 iv%instid(inst)%loc_dym(n1:n2), &
225 iv%instid(inst)%ps(n1:n2),n1n2)
226
227 do n=n1,n2
228 if (iv%instid(inst)%isflg(n) == 0) then ! over sea using SST
229 call da_interp_lin_2d( xb % tgrn, &
230 xp%ims, xp%ime, xp%jms, xp%jme, &
231 iv%instid(inst)%loc_i(n), &
232 iv%instid(inst)%loc_j(n), &
233 iv%instid(inst)%loc_dx(n), &
234 iv%instid(inst)%loc_dy(n), &
235 iv%instid(inst)%loc_dxm(n), &
236 iv%instid(inst)%loc_dym(n), &
237 iv%instid(inst)%ts(n))
238 else
239 call da_interp_lin_2d( xb % tsk, &
240 xp%ims, xp%ime, xp%jms, xp%jme, &
241 iv%instid(inst)%loc_i(n), &
242 iv%instid(inst)%loc_j(n), &
243 iv%instid(inst)%loc_dx(n), &
244 iv%instid(inst)%loc_dy(n), &
245 iv%instid(inst)%loc_dxm(n), &
246 iv%instid(inst)%loc_dym(n), &
247 iv%instid(inst)%ts(n))
248 end if
249
250 end do
251 call da_interp_lin_2d_new(xb % terr, &
252 xp%ims, xp%ime, xp%jms, xp%jme, &
253 iv%instid(inst)%loc_i(n1:n2), &
254 iv%instid(inst)%loc_j(n1:n2), &
255 iv%instid(inst)%loc_dx(n1:n2), &
256 iv%instid(inst)%loc_dy(n1:n2), &
257 iv%instid(inst)%loc_dxm(n1:n2), &
258 iv%instid(inst)%loc_dym(n1:n2), &
259 iv%instid(inst)%elevation(n1:n2),n1n2 )
260
261 ! variables for emissivity calculations
262 !----------------------------------------
263 call da_interp_lin_2d_new( xb % smois, xp%ims, xp%ime, xp%jms, xp%jme, &
264 iv%instid(inst)%loc_i(n1:n2), &
265 iv%instid(inst)%loc_j(n1:n2), &
266 iv%instid(inst)%loc_dx(n1:n2), &
267 iv%instid(inst)%loc_dy(n1:n2), &
268 iv%instid(inst)%loc_dxm(n1:n2), &
269 iv%instid(inst)%loc_dym(n1:n2), &
270 iv%instid(inst)%smois(n1:n2),n1n2)
271
272 call da_interp_lin_2d_new( xb % tslb, xp%ims, xp%ime, xp%jms, xp%jme, &
273 iv%instid(inst)%loc_i(n1:n2), &
274 iv%instid(inst)%loc_j(n1:n2), &
275 iv%instid(inst)%loc_dx(n1:n2), &
276 iv%instid(inst)%loc_dy(n1:n2), &
277 iv%instid(inst)%loc_dxm(n1:n2), &
278 iv%instid(inst)%loc_dym(n1:n2), &
279 iv%instid(inst)%tslb(n1:n2),n1n2)
280
281 call da_interp_lin_2d_new( xb % snowh, xp%ims, xp%ime, xp%jms, xp%jme, &
282 iv%instid(inst)%loc_i(n1:n2), &
283 iv%instid(inst)%loc_j(n1:n2), &
284 iv%instid(inst)%loc_dx(n1:n2), &
285 iv%instid(inst)%loc_dy(n1:n2), &
286 iv%instid(inst)%loc_dxm(n1:n2), &
287 iv%instid(inst)%loc_dym(n1:n2), &
288 iv%instid(inst)%snowh(n1:n2),n1n2)
289
290 call da_interp_lin_2d_new( xb % vegfra, xp%ims, xp%ime, xp%jms, xp%jme, &
291 iv%instid(inst)%loc_i(n1:n2), &
292 iv%instid(inst)%loc_j(n1:n2), &
293 iv%instid(inst)%loc_dx(n1:n2), &
294 iv%instid(inst)%loc_dy(n1:n2), &
295 iv%instid(inst)%loc_dxm(n1:n2), &
296 iv%instid(inst)%loc_dym(n1:n2), &
297 iv%instid(inst)%vegfra(n1:n2),n1n2)
298
299 ! iv%instid(inst)%snowh(n1:n2) = iv%instid(inst)%snowh(n1:n2)*100.0 ! convert from m to mm
300 iv%instid(inst)%vegfra(n1:n2) = 0.01*iv%instid(inst)%vegfra(n1:n2) ! convert range to 0~1
301
302 ! variables for calculation of cloud affected radiance
303 !-------------------------------------------------------
304 do k=xp%kts,xp%kte
305 call da_interp_lin_2d_new( xb%t(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
306 iv%instid(inst)%loc_i(n1:n2), &
307 iv%instid(inst)%loc_j(n1:n2), &
308 iv%instid(inst)%loc_dx(n1:n2), &
309 iv%instid(inst)%loc_dy(n1:n2), &
310 iv%instid(inst)%loc_dxm(n1:n2), &
311 iv%instid(inst)%loc_dym(n1:n2), &
312 iv%instid(inst)%tm(k,n1:n2),n1n2)
313
314 call da_interp_lin_2d_new( xb%q(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
315 iv%instid(inst)%loc_i(n1:n2), &
316 iv%instid(inst)%loc_j(n1:n2), &
317 iv%instid(inst)%loc_dx(n1:n2), &
318 iv%instid(inst)%loc_dy(n1:n2), &
319 iv%instid(inst)%loc_dxm(n1:n2), &
320 iv%instid(inst)%loc_dym(n1:n2), &
321 iv%instid(inst)%qm(k,n1:n2),n1n2)
322
323 call da_interp_lin_2d_new( xb%qrn(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
324 iv%instid(inst)%loc_i(n1:n2), &
325 iv%instid(inst)%loc_j(n1:n2), &
326 iv%instid(inst)%loc_dx(n1:n2), &
327 iv%instid(inst)%loc_dy(n1:n2), &
328 iv%instid(inst)%loc_dxm(n1:n2), &
329 iv%instid(inst)%loc_dym(n1:n2), &
330 iv%instid(inst)%qrn(k,n1:n2),n1n2)
331
332 call da_interp_lin_2d_new( xb%qcw(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
333 iv%instid(inst)%loc_i(n1:n2), &
334 iv%instid(inst)%loc_j(n1:n2), &
335 iv%instid(inst)%loc_dx(n1:n2), &
336 iv%instid(inst)%loc_dy(n1:n2), &
337 iv%instid(inst)%loc_dxm(n1:n2), &
338 iv%instid(inst)%loc_dym(n1:n2), &
339 iv%instid(inst)%qcw(k,n1:n2),n1n2)
340
341 call da_interp_lin_2d_new( xb%qci(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
342 iv%instid(inst)%loc_i(n1:n2), &
343 iv%instid(inst)%loc_j(n1:n2), &
344 iv%instid(inst)%loc_dx(n1:n2), &
345 iv%instid(inst)%loc_dy(n1:n2), &
346 iv%instid(inst)%loc_dxm(n1:n2), &
347 iv%instid(inst)%loc_dym(n1:n2), &
348 iv%instid(inst)%qci(k,n1:n2),n1n2)
349 call da_interp_lin_2d_new( xb%qsn(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
350 iv%instid(inst)%loc_i(n1:n2), &
351 iv%instid(inst)%loc_j(n1:n2), &
352 iv%instid(inst)%loc_dx(n1:n2), &
353 iv%instid(inst)%loc_dy(n1:n2), &
354 iv%instid(inst)%loc_dxm(n1:n2), &
355 iv%instid(inst)%loc_dym(n1:n2), &
356 iv%instid(inst)%qsn(k,n1:n2),n1n2)
357
358 call da_interp_lin_2d_new( xb%qgr(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
359 iv%instid(inst)%loc_i(n1:n2), &
360 iv%instid(inst)%loc_j(n1:n2), &
361 iv%instid(inst)%loc_dx(n1:n2), &
362 iv%instid(inst)%loc_dy(n1:n2), &
363 iv%instid(inst)%loc_dxm(n1:n2), &
364 iv%instid(inst)%loc_dym(n1:n2), &
365 iv%instid(inst)%qgr(k,n1:n2),n1n2)
366
367 iv%instid(inst)%pm(k,n1:n2) = v_p(k,n1:n2)
368
369 end do
370
371 ! convert to hPa
372 iv%instid(inst)%ps(n1:n2) = 0.01 * iv%instid(inst)%ps(n1:n2)
373 iv%instid(inst)%mr2m(n1:n2) = iv%instid(inst)%q2m(n1:n2) * q2ppmv
374
375 ! ADD for computing cloud liquid water path (mm) from guess
376
377 pf(xp%kts,n1:n2) = 100.0*iv%instid(inst)%ps(n1:n2)
378 do k = 2,xp%kte
379 pf(k,n1:n2) = 50.0*(v_p(k-1,n1:n2)+v_p(k,n1:n2))
380 end do
381 pf(xp%kte+1,n1:n2)= 50.0*v_p(xp%kte,n1:n2)
382
383 iv%instid(inst)%clwp(n1:n2) = 0.0
384 do k = xp%kts,xp%kte
385 dpf(k,n1:n2) = pf(k,n1:n2) - pf(k+1,n1:n2)
386 clw(k,n1:n2) = iv%instid(inst)%qcw(k,n1:n2)*dpf(k,n1:n2)/gravity
387 where (v_p(k,n1:n2)<100.0) clw (k,n1:n2) = 0.0
388 iv%instid(inst)%clwp(n1:n2) =iv%instid(inst)%clwp(n1:n2) + clw(k,n1:n2)
389 end do
390
391 ! call rttov_q2v (&
392 ! gas_unit_specconc, &! in
393 ! iv%instid(inst)%q2m(n1:n2), &! in
394 ! gas_id_watervapour, &! in
395 ! iv%instid(inst)%q2m(n1:n2), &! in
396 ! iv%instid(inst)%mr2m(n1:n2) ! inout
397
398 ! Compute surface emissivity for observation
399 !-------------------------------------------
400
401 zasat(n1:n2) = iv%instid(inst)%satzen(n1:n2) * deg2rad
402 zlsat(n1:n2) = iv%instid(inst)%satazi(n1:n2) * deg2rad
403
404 ! do i = n1,n2
405 ! knchan(i) = nchan
406 ! do kn=1,nchan
407 ! kprof((i-1)*nchan+kn) = i
408 ! kchan((i-1)*nchan+kn) = kn
409 ! kochan((i-1)*nchan+kn) = kn
410 ! indexn(kn,i) = (-1)*nchan+kn
411 ! end do
412 ! end do
413
414 pemsv(:) = 0.0
415 pemsh(:) = 0.0
416
417 if (amsua) then
418 polar(1:4) = 0.0 ! polarizition: 0:vertical; 1:horizontal
419 polar(5:6) = 1.0
420 polar(7) = 0.0
421 polar(8:14)= 1.0
422 polar(15) = 0.0
423 end if
424 if (amsub) then
425 polar(1:5) = 0.0
426 end if
427
428 do n=n1,n2
429 knchan(n) = nchan
430 do kn=1,nchan
431 kprof((n-1)*nchan+kn) = n
432 kchan((n-1)*nchan+kn) = kn
433 kochan((n-1)*nchan+kn) = kn
434 indexn(kn,n) = (n-1)*nchan+kn
435 end do
436 end do
437
438 if (coefs(inst)%id_sensor == 2) then ! microwave
439 call da_trace("da_get_innov_vector_rttov",message="Calling gsi_emiss")
440 call gsi_emiss(inst,knchpf,kprof,kchan,knchan,indexn, &
441 iv%instid(inst)%satzen * deg2rad, &
442 iv%instid(inst)%satazi * deg2rad, polar, &
443 iv%instid(inst)%isflg,pemsv,pemsh,pems5, &
444 iv%instid(inst)%ts, &
445 iv%instid(inst)%soiltyp, &
446 iv%instid(inst)%tslb, &
447 iv%instid(inst)%smois, &
448 iv%instid(inst)%vegtyp, &
449 iv%instid(inst)%vegfra, &
450 iv%instid(inst)%snowh, &
451 sensor_type,n1n2,nchan, &
452 iv%instid(inst)%tb_inv(:,:),amsua,amsub,ssmi, &
453 iv%instid(inst)%u10, &
454 iv%instid(inst)%v10)
455 end if
456
457 do n=n1,n2
458 con_vars(n) % nlevels = nlevels
459 allocate (con_vars(n) % t(nlevels))
460 allocate (con_vars(n) % q(nlevels))
461
462 con_vars(n) % t(1:nlevels) = iv%instid(inst)%t(1:nlevels,n)
463 con_vars(n) % q(1:nlevels) = iv%instid(inst)%mr(1:nlevels,n)
464 con_vars(n) % ps = iv%instid(inst)%ps(n)
465
466 aux_vars(n) % t2m = iv%instid(inst)%t2m(n)
467 aux_vars(n) % q2m = iv%instid(inst)%mr2m(n)
468 aux_vars(n) % u10 = iv%instid(inst)%u10(n)
469 aux_vars(n) % v10 = iv%instid(inst)%v10(n)
470 aux_vars(n) % surftype = iv%instid(inst)%landsea_mask(n)
471 aux_vars(n) % surft = iv%instid(inst)%ts(n)
472 aux_vars(n) % satzen = iv%instid(inst)%satzen(n)
473 aux_vars(n) % satazi = iv%instid(inst)%satazi(n)
474
475 ! [1.3] Call RTM forward model
476
477 ! da_rttov_direct nominally an array version, but doesn't handle arrays
478 ! of surface flags properly
479 allocate(temp(nchan))
480 call da_rttov_direct (inst, iv%instid(inst)%isflg(n), nchan, 1, &
481 con_vars(n:n), aux_vars(n:n), &
482 ! iv%instid(inst)%tb_xb(:,n), pemsv, pemsh, &
483 temp(:), pemsv((n-1)*nchan+1:n*nchan), pemsh((n-1)*nchan+1:n*nchan), &
484 iv%instid(inst)%emiss(:,n))
485 iv%instid(inst)%tb_xb(:,n)=temp(:)
486 deallocate(temp)
487 do k=1,nchan
488 ! if (iv%instid(inst)%tb_inv(k,n) > 15.0) then
489 ! write (0,*) "n,k,tb,xb",n,k, ob%instid(inst)%tb(k,n),iv%instid(inst)%tb_xb(k,n)
490 ! end if
491 end do
492 end do
493
494 if (iv%instid(inst)%nchannels > nchan) then
495 iv%instid(inst)%emiss(nchan+1:iv%instid(inst)%nchannels,n1:n2) = 0.0
496 end if
497
498 if (coefs(inst)%id_sensor == 2) then ! microwave
499 do n=n1,n2
500 iv%instid(inst)%emiss(1:nchan,n) = pems5(1+(n-n1)*nchan:(n-n1+1)*nchan)
501 end do
502 end if
503
504 do n=n1,n2
505 deallocate (con_vars(n) % t)
506 deallocate (con_vars(n) % q)
507 end do
508
509 !----------------------------------------------------------------
510 ! [2.0] calculate components of innovation vector:
511 !----------------------------------------------------------------
512
513 where (iv%instid(inst)%tb_inv(1:nchan,n1:n2) > missing_r)
514 iv%instid(inst)%tb_inv(1:nchan,n1:n2) = &
515 ob%instid(inst)%tb(1:nchan,n1:n2) - iv%instid(inst)%tb_xb(1:nchan,n1:n2)
516 elsewhere
517 iv%instid(inst)%tb_inv(1:nchan,n1:n2) = missing_r
518 endwhere
519 !do n=n1,n2
520 ! do k=1,nchan
521 ! if (iv%instid(inst)%tb_inv(k,n) > 15.0) then
522 ! write (0,*) "n,k,tb,xb,inv",n,k, ob%instid(inst)%tb(k,n),iv%instid(inst)%tb_xb(k,n),iv%instid(inst)%tb_inv(k,n)
523 ! end if
524 ! end do
525 !end do
526
527 deallocate (v_p)
528 deallocate (clw)
529 deallocate (dpf)
530 deallocate (pf)
531 deallocate (pres)
532 deallocate (con_vars)
533 deallocate (aux_vars)
534 deallocate (zasat)
535 deallocate (zlsat)
536 deallocate (polar)
537 deallocate (knchan)
538 deallocate (kprof)
539 deallocate (kchan)
540 deallocate (kochan)
541 deallocate (indexn)
542 deallocate (pemsv)
543 deallocate (pemsh)
544 deallocate (pems5)
545
546 !------------------------------------------------------------------------
547 ! [3.0] Perform bias correction
548 !------------------------------------------------------------------------
549
550 if (biascorr) then
551 write(unit=stdout,fmt='(A,A)') 'Performing bias correction for ', &
552 trim(iv%instid(inst)%rttovid_string)
553 call da_biascorr(inst,ob,iv)
554 end if
555
556 call da_togrid_new (iv%instid(inst)%zk(:,n1:n2), xp%kts, xp%kte, &
557 iv%instid(inst)%loc_k(:,n1:n2), &
558 iv%instid(inst)%loc_dz(:,n1:n2), &
559 iv%instid(inst)%loc_dzm(:,n1:n2))
560 end do ! end loop for sensor
561
562 !------------------------------------------------------------------------
563 ! [4.0] Perform QC check
564 !------------------------------------------------------------------------
565
566 if (qc_rad) then
567 call da_qc_rad(ob, iv)
568 end if
569
570 !------------------------------------------
571 ! 5.0 prepare bias statistics files
572 !------------------------------------------
573
574 if (biasprep) then
575 do inst = 1, iv%num_inst
576 write(unit=stdout,fmt='(A,A)') 'Preparing bias statistics files for ', &
577 trim(iv%instid(inst)%rttovid_string)
578 call da_biasprep(inst,ob,iv)
579 end do
580 end if
581
582 if (trace_use) call da_trace_exit("da_get_innov_vector_rttov")
583 #endif
584
585 end subroutine da_get_innov_vector_rttov
586
587