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 :: vtm(:,:), dlnpf(:,:), 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(vtm(xp%kms:xp%kme,n1:n2))
104 allocate(clw(xp%kms:xp%kme,n1:n2))
105 allocate(dlnpf(xp%kms:xp%kme,n1:n2))
106 allocate(pf(xp%kms:xp%kme+1,n1:n2))
107
108 ! horizontal interpolate xb pressure to ob position for every xb levels
109 do n=n1,n2
110 do k=xp%kts,xp%kte
111 ! convert to mb
112 v_p(k,n) = 0.01*(iv%instid(inst)%loc_dym(n)*( &
113 iv%instid(inst)%loc_dxm(n)*xb%p(iv%instid(inst)%loc_i(n), iv%instid(inst)%loc_j(n),k) + &
114 iv%instid(inst)%loc_dx(n) *xb%p(iv%instid(inst)%loc_i(n)+1,iv%instid(inst)%loc_j(n),k)) + &
115 iv%instid(inst)%loc_dy(n) *( &
116 iv%instid(inst)%loc_dxm(n)*xb%p(iv%instid(inst)%loc_i(n), iv%instid(inst)%loc_j(n)+1,k) + &
117 iv%instid(inst)%loc_dx(n) *xb%p(iv%instid(inst)%loc_i(n)+1,iv%instid(inst)%loc_j(n)+1,k)))
118 end do
119 end do
120
121 call da_to_zk_new(pres, v_p(:,n1:n2), xp, v_interp_p, &
122 iv%instid(inst)%zk(:,n1:n2), n1n2,nlevels )
123
124 call da_togrid_new (iv%instid(inst)%zk(:,n1:n2), xp%kts, xp%kte, &
125 iv%instid(inst)%loc_k(:,n1:n2), &
126 iv%instid(inst)%loc_dz(:,n1:n2), &
127 iv%instid(inst)%loc_dzm(:,n1:n2))
128 ! [1.2] Interpolate horizontally to ob:
129
130 call da_interp_lin_3d_new( xb%t, xp, &
131 iv%instid(inst)%loc_i(n1:n2), &
132 iv%instid(inst)%loc_j(n1:n2), &
133 iv%instid(inst)%loc_k(:,n1:n2), &
134 iv%instid(inst)%loc_dx(n1:n2), &
135 iv%instid(inst)%loc_dy(n1:n2), &
136 iv%instid(inst)%loc_dz(:,n1:n2), &
137 iv%instid(inst)%loc_dxm(n1:n2), &
138 iv%instid(inst)%loc_dym(n1:n2), &
139 iv%instid(inst)%loc_dzm(:,n1:n2), &
140 iv%instid(inst)%t(:,n1:n2), nlevels, n1n2 )
141
142 call da_interp_lin_3d_new( xb%q, xp, &
143 iv%instid(inst)%loc_i(n1:n2), &
144 iv%instid(inst)%loc_j(n1:n2), &
145 iv%instid(inst)%loc_k(:,n1:n2), &
146 iv%instid(inst)%loc_dx(n1:n2), &
147 iv%instid(inst)%loc_dy(n1:n2), &
148 iv%instid(inst)%loc_dz(:,n1:n2), &
149 iv%instid(inst)%loc_dxm(n1:n2), &
150 iv%instid(inst)%loc_dym(n1:n2), &
151 iv%instid(inst)%loc_dzm(:,n1:n2), &
152 iv%instid(inst)%mr(:,n1:n2), nlevels, n1n2 )
153
154 do n= n1,n2
155 do k=1, nlevels
156 if (iv%instid(inst)%zk(k,n) <= 0.0) then
157 iv%instid(inst)%t(k,n) = coefs(inst) % ref_prfl_t(k,gas_id_watervapour) ! outside model level
158 iv%instid(inst)%mr(k,n) = coefs(inst) % ref_prfl_mr(k,gas_id_watervapour)
159 else
160 iv%instid(inst)%mr(k,n) = iv%instid(inst)%mr(k,n) * q2ppmv
161 ! call rttov_q2v (&
162 ! gas_unit_specconc, &! in
163 ! iv%instid(inst)%q(k,n), &! in
164 ! gas_id_watervapour, &! in
165 ! iv%instid(inst)%q(k,n), &! in
166 ! iv%instid(inst)%mr(k,n) ) ! inout
167 end if
168 if (pres(k) < 100.0) iv%instid(inst)%mr(k,n) = coefs(inst) % ref_prfl_mr(k,gas_id_watervapour)
169 end do
170
171 ! determine surface type of obs location
172 !-----------------------------------------
173 call da_detsurtyp( xb%snow, xb%xice, xb%landmask, &
174 xb%ivgtyp, xb%isltyp, &
175 xp%ims, xp%ime, xp%jms, xp%jme, &
176 iv%instid(inst)%loc_i(n), iv%instid(inst)%loc_j(n), &
177 iv%instid(inst)%loc_dx(n), iv%instid(inst)%loc_dy(n), &
178 iv%instid(inst)%loc_dxm(n), iv%instid(inst)%loc_dym(n), &
179 iv%instid(inst)%isflg(n),iv%instid(inst)%vegtyp(n), iv%instid(inst)%soiltyp(n), &
180 seap, icep, lndp, snop )
181 end do
182
183 call da_interp_lin_2d_new (xb % u10, xp%ims, xp%ime, xp%jms, xp%jme, &
184 iv%instid(inst)%loc_i(n1:n2), &
185 iv%instid(inst)%loc_j(n1:n2), &
186 iv%instid(inst)%loc_dx(n1:n2), &
187 iv%instid(inst)%loc_dy(n1:n2), &
188 iv%instid(inst)%loc_dxm(n1:n2), &
189 iv%instid(inst)%loc_dym(n1:n2), &
190 iv%instid(inst)%u10(n1:n2),n1n2)
191
192 call da_interp_lin_2d_new (xb % v10, xp%ims, xp%ime, xp%jms, xp%jme, &
193 iv%instid(inst)%loc_i(n1:n2), &
194 iv%instid(inst)%loc_j(n1:n2), &
195 iv%instid(inst)%loc_dx(n1:n2), &
196 iv%instid(inst)%loc_dy(n1:n2), &
197 iv%instid(inst)%loc_dxm(n1:n2), &
198 iv%instid(inst)%loc_dym(n1:n2), &
199 iv%instid(inst)%v10(n1:n2),n1n2)
200
201 call da_interp_lin_2d_new (xb % t2, xp%ims, xp%ime, xp%jms, xp%jme, &
202 iv%instid(inst)%loc_i(n1:n2), &
203 iv%instid(inst)%loc_j(n1:n2), &
204 iv%instid(inst)%loc_dx(n1:n2), &
205 iv%instid(inst)%loc_dy(n1:n2), &
206 iv%instid(inst)%loc_dxm(n1:n2), &
207 iv%instid(inst)%loc_dym(n1:n2), &
208 iv%instid(inst)%t2m(n1:n2) ,n1n2)
209
210 call da_interp_lin_2d_new (xb % q2, xp%ims, xp%ime, xp%jms, xp%jme, &
211 iv%instid(inst)%loc_i(n1:n2), &
212 iv%instid(inst)%loc_j(n1:n2), &
213 iv%instid(inst)%loc_dx(n1:n2), &
214 iv%instid(inst)%loc_dy(n1:n2), &
215 iv%instid(inst)%loc_dxm(n1:n2), &
216 iv%instid(inst)%loc_dym(n1:n2), &
217 iv%instid(inst)%q2m(n1:n2),n1n2)
218
219 call da_interp_lin_2d_new(xb % psfc, xp%ims, xp%ime, xp%jms, xp%jme, &
220 iv%instid(inst)%loc_i(n1:n2), &
221 iv%instid(inst)%loc_j(n1:n2), &
222 iv%instid(inst)%loc_dx(n1:n2), &
223 iv%instid(inst)%loc_dy(n1:n2), &
224 iv%instid(inst)%loc_dxm(n1:n2), &
225 iv%instid(inst)%loc_dym(n1:n2), &
226 iv%instid(inst)%ps(n1:n2),n1n2)
227
228 do n=n1,n2
229 if (iv%instid(inst)%isflg(n) == 0) then ! over sea using SST
230 call da_interp_lin_2d( xb % tgrn, &
231 xp%ims, xp%ime, xp%jms, xp%jme, &
232 iv%instid(inst)%loc_i(n), &
233 iv%instid(inst)%loc_j(n), &
234 iv%instid(inst)%loc_dx(n), &
235 iv%instid(inst)%loc_dy(n), &
236 iv%instid(inst)%loc_dxm(n), &
237 iv%instid(inst)%loc_dym(n), &
238 iv%instid(inst)%ts(n))
239 else
240 call da_interp_lin_2d( xb % tsk, &
241 xp%ims, xp%ime, xp%jms, xp%jme, &
242 iv%instid(inst)%loc_i(n), &
243 iv%instid(inst)%loc_j(n), &
244 iv%instid(inst)%loc_dx(n), &
245 iv%instid(inst)%loc_dy(n), &
246 iv%instid(inst)%loc_dxm(n), &
247 iv%instid(inst)%loc_dym(n), &
248 iv%instid(inst)%ts(n))
249 end if
250
251 end do
252 call da_interp_lin_2d_new(xb % terr, &
253 xp%ims, xp%ime, xp%jms, xp%jme, &
254 iv%instid(inst)%loc_i(n1:n2), &
255 iv%instid(inst)%loc_j(n1:n2), &
256 iv%instid(inst)%loc_dx(n1:n2), &
257 iv%instid(inst)%loc_dy(n1:n2), &
258 iv%instid(inst)%loc_dxm(n1:n2), &
259 iv%instid(inst)%loc_dym(n1:n2), &
260 iv%instid(inst)%elevation(n1:n2),n1n2 )
261
262 ! variables for emissivity calculations
263 !----------------------------------------
264 call da_interp_lin_2d_new( xb % smois, xp%ims, xp%ime, xp%jms, xp%jme, &
265 iv%instid(inst)%loc_i(n1:n2), &
266 iv%instid(inst)%loc_j(n1:n2), &
267 iv%instid(inst)%loc_dx(n1:n2), &
268 iv%instid(inst)%loc_dy(n1:n2), &
269 iv%instid(inst)%loc_dxm(n1:n2), &
270 iv%instid(inst)%loc_dym(n1:n2), &
271 iv%instid(inst)%smois(n1:n2),n1n2)
272
273 call da_interp_lin_2d_new( xb % tslb, xp%ims, xp%ime, xp%jms, xp%jme, &
274 iv%instid(inst)%loc_i(n1:n2), &
275 iv%instid(inst)%loc_j(n1:n2), &
276 iv%instid(inst)%loc_dx(n1:n2), &
277 iv%instid(inst)%loc_dy(n1:n2), &
278 iv%instid(inst)%loc_dxm(n1:n2), &
279 iv%instid(inst)%loc_dym(n1:n2), &
280 iv%instid(inst)%tslb(n1:n2),n1n2)
281
282 call da_interp_lin_2d_new( xb % snowh, xp%ims, xp%ime, xp%jms, xp%jme, &
283 iv%instid(inst)%loc_i(n1:n2), &
284 iv%instid(inst)%loc_j(n1:n2), &
285 iv%instid(inst)%loc_dx(n1:n2), &
286 iv%instid(inst)%loc_dy(n1:n2), &
287 iv%instid(inst)%loc_dxm(n1:n2), &
288 iv%instid(inst)%loc_dym(n1:n2), &
289 iv%instid(inst)%snowh(n1:n2),n1n2)
290
291 call da_interp_lin_2d_new( xb % vegfra, xp%ims, xp%ime, xp%jms, xp%jme, &
292 iv%instid(inst)%loc_i(n1:n2), &
293 iv%instid(inst)%loc_j(n1:n2), &
294 iv%instid(inst)%loc_dx(n1:n2), &
295 iv%instid(inst)%loc_dy(n1:n2), &
296 iv%instid(inst)%loc_dxm(n1:n2), &
297 iv%instid(inst)%loc_dym(n1:n2), &
298 iv%instid(inst)%vegfra(n1:n2),n1n2)
299
300 ! iv%instid(inst)%snowh(n1:n2) = iv%instid(inst)%snowh(n1:n2)*100.0 ! convert from m to mm
301 iv%instid(inst)%vegfra(n1:n2) = 0.01*iv%instid(inst)%vegfra(n1:n2) ! convert range to 0~1
302
303 ! variables for calculation of cloud affected radiance
304 !-------------------------------------------------------
305 do k=xp%kts,xp%kte
306 call da_interp_lin_2d_new( xb%t(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
307 iv%instid(inst)%loc_i(n1:n2), &
308 iv%instid(inst)%loc_j(n1:n2), &
309 iv%instid(inst)%loc_dx(n1:n2), &
310 iv%instid(inst)%loc_dy(n1:n2), &
311 iv%instid(inst)%loc_dxm(n1:n2), &
312 iv%instid(inst)%loc_dym(n1:n2), &
313 iv%instid(inst)%tm(k,n1:n2),n1n2)
314
315 call da_interp_lin_2d_new( xb%q(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
316 iv%instid(inst)%loc_i(n1:n2), &
317 iv%instid(inst)%loc_j(n1:n2), &
318 iv%instid(inst)%loc_dx(n1:n2), &
319 iv%instid(inst)%loc_dy(n1:n2), &
320 iv%instid(inst)%loc_dxm(n1:n2), &
321 iv%instid(inst)%loc_dym(n1:n2), &
322 iv%instid(inst)%qm(k,n1:n2),n1n2)
323
324 call da_interp_lin_2d_new( xb%qrn(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
325 iv%instid(inst)%loc_i(n1:n2), &
326 iv%instid(inst)%loc_j(n1:n2), &
327 iv%instid(inst)%loc_dx(n1:n2), &
328 iv%instid(inst)%loc_dy(n1:n2), &
329 iv%instid(inst)%loc_dxm(n1:n2), &
330 iv%instid(inst)%loc_dym(n1:n2), &
331 iv%instid(inst)%qrn(k,n1:n2),n1n2)
332
333 call da_interp_lin_2d_new( xb%qcw(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
334 iv%instid(inst)%loc_i(n1:n2), &
335 iv%instid(inst)%loc_j(n1:n2), &
336 iv%instid(inst)%loc_dx(n1:n2), &
337 iv%instid(inst)%loc_dy(n1:n2), &
338 iv%instid(inst)%loc_dxm(n1:n2), &
339 iv%instid(inst)%loc_dym(n1:n2), &
340 iv%instid(inst)%qcw(k,n1:n2),n1n2)
341
342 call da_interp_lin_2d_new( xb%qci(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
343 iv%instid(inst)%loc_i(n1:n2), &
344 iv%instid(inst)%loc_j(n1:n2), &
345 iv%instid(inst)%loc_dx(n1:n2), &
346 iv%instid(inst)%loc_dy(n1:n2), &
347 iv%instid(inst)%loc_dxm(n1:n2), &
348 iv%instid(inst)%loc_dym(n1:n2), &
349 iv%instid(inst)%qci(k,n1:n2),n1n2)
350 call da_interp_lin_2d_new( xb%qsn(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
351 iv%instid(inst)%loc_i(n1:n2), &
352 iv%instid(inst)%loc_j(n1:n2), &
353 iv%instid(inst)%loc_dx(n1:n2), &
354 iv%instid(inst)%loc_dy(n1:n2), &
355 iv%instid(inst)%loc_dxm(n1:n2), &
356 iv%instid(inst)%loc_dym(n1:n2), &
357 iv%instid(inst)%qsn(k,n1:n2),n1n2)
358
359 call da_interp_lin_2d_new( xb%qgr(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
360 iv%instid(inst)%loc_i(n1:n2), &
361 iv%instid(inst)%loc_j(n1:n2), &
362 iv%instid(inst)%loc_dx(n1:n2), &
363 iv%instid(inst)%loc_dy(n1:n2), &
364 iv%instid(inst)%loc_dxm(n1:n2), &
365 iv%instid(inst)%loc_dym(n1:n2), &
366 iv%instid(inst)%qgr(k,n1:n2),n1n2)
367
368 iv%instid(inst)%pm(:,n1:n2) = v_p(:,n1:n2)
369
370 ! ADD for computing cloud liquid water path
371 ! virtual T
372 vtm (k,n1:n2) = (1.+0.608*iv%instid(inst)%qm(k,n1:n2))*iv%instid(inst)%tm(k,n1:n2)
373 ! convert kg/kg to g/m3
374 ! clw (k,n1:n2) = iv%instid(inst)%qcw(k,n1:n2)*v_p(k,n1:n2)*100.*0.0289644/8.31451/iv%instid(inst)%tm(k,n1:n2)
375 clw(k,n1:n2) = 0.34836*iv%instid(inst)%qcw(k,n1:n2)*v_p(k,n1:n2)/iv%instid(inst)%tm(k,n1:n2)
376 where (v_p(k,n1:n2)<100.0) clw (k,n1:n2) = 0.0
377
378 end do
379
380 ! convert to hPa
381 iv%instid(inst)%ps(n1:n2) = 0.01 * iv%instid(inst)%ps(n1:n2)
382 iv%instid(inst)%mr2m(n1:n2) = iv%instid(inst)%q2m(n1:n2) * q2ppmv
383
384 ! ADD for computing cloud liquid water path (mm) from guess
385
386 pf(xp%kts,n1:n2) = 100.0*iv%instid(inst)%ps(n1:n2)
387 do k = 2,xp%kte
388 ! pf(k,n1:n2) = 100.0*(v_p(k-1,n1:n2)+v_p(k,n1:n2))*0.5 ! full level pressure
389 pf(k,n1:n2) = 50.0*(v_p(k-1,n1:n2)+v_p(k,n1:n2))
390 end do
391 ! pf(xp%kte+1,n1:n2)= 100.0*v_p(xp%kte,n1:n2)*0.5
392 pf(xp%kte+1,n1:n2)= 50.0*v_p(xp%kte,n1:n2)
393
394 iv%instid(inst)%clwp(n1:n2) = 0.0
395 do k = xp%kts,xp%kte
396 dlnpf(k,n1:n2) = log(pf(k,n1:n2)) - log(pf(k+1,n1:n2))
397 ! clw (k,n1:n2) = clw(k,n1:n2)*vtm(k,n1:n2)*287.05/9.80665*dlnpf(k,n1:n2)
398 clw(k,n1:n2) = 29.27095*clw(k,n1:n2)*vtm(k,n1:n2)*dlnpf(k,n1:n2)
399 iv%instid(inst)%clwp(n1:n2) =iv%instid(inst)%clwp(n1:n2) + clw(k,n1:n2)
400 end do
401 ! iv%instid(inst)%clwp = 0.001*iv%instid(inst)%clwp ! kg/m2 = mm
402
403 ! call rttov_q2v (&
404 ! gas_unit_specconc, &! in
405 ! iv%instid(inst)%q2m(n1:n2), &! in
406 ! gas_id_watervapour, &! in
407 ! iv%instid(inst)%q2m(n1:n2), &! in
408 ! iv%instid(inst)%mr2m(n1:n2) ! inout
409
410 ! Compute surface emissivity for observation
411 !-------------------------------------------
412
413 zasat(n1:n2) = iv%instid(inst)%satzen(n1:n2) * deg2rad
414 zlsat(n1:n2) = iv%instid(inst)%satazi(n1:n2) * deg2rad
415
416 ! do i = n1,n2
417 ! knchan(i) = nchan
418 ! do kn=1,nchan
419 ! kprof((i-1)*nchan+kn) = i
420 ! kchan((i-1)*nchan+kn) = kn
421 ! kochan((i-1)*nchan+kn) = kn
422 ! indexn(kn,i) = (-1)*nchan+kn
423 ! end do
424 ! end do
425
426 pemsv(:) = 0.0
427 pemsh(:) = 0.0
428
429 if (amsua) then
430 polar(1:4) = 0.0 ! polarizition: 0:vertical; 1:horizontal
431 polar(5:6) = 1.0
432 polar(7) = 0.0
433 polar(8:14)= 1.0
434 polar(15) = 0.0
435 end if
436 if (amsub) then
437 polar(1:5) = 0.0
438 end if
439
440 do n=n1,n2
441 knchan(n) = nchan
442 do kn=1,nchan
443 kprof((n-1)*nchan+kn) = n
444 kchan((n-1)*nchan+kn) = kn
445 kochan((n-1)*nchan+kn) = kn
446 indexn(kn,n) = (n-1)*nchan+kn
447 end do
448 end do
449
450 if (coefs(inst)%id_sensor == 2) then ! microwave
451 call da_trace("da_get_innov_vector_rttov",message="Calling gsi_emiss")
452 call gsi_emiss(inst,knchpf,kprof,kchan,knchan,indexn, &
453 iv%instid(inst)%satzen * deg2rad, &
454 iv%instid(inst)%satazi * deg2rad, polar, &
455 iv%instid(inst)%isflg,pemsv,pemsh,pems5, &
456 iv%instid(inst)%ts, &
457 iv%instid(inst)%soiltyp, &
458 iv%instid(inst)%tslb, &
459 iv%instid(inst)%smois, &
460 iv%instid(inst)%vegtyp, &
461 iv%instid(inst)%vegfra, &
462 iv%instid(inst)%snowh, &
463 sensor_type,n1n2,nchan, &
464 iv%instid(inst)%tb_inv(:,:),amsua,amsub,ssmi, &
465 iv%instid(inst)%u10, &
466 iv%instid(inst)%v10)
467 end if
468
469 do n=n1,n2
470 con_vars(n) % nlevels = nlevels
471 allocate (con_vars(n) % t(nlevels))
472 allocate (con_vars(n) % q(nlevels))
473
474 con_vars(n) % t(1:nlevels) = iv%instid(inst)%t(1:nlevels,n)
475 con_vars(n) % q(1:nlevels) = iv%instid(inst)%mr(1:nlevels,n)
476 con_vars(n) % ps = iv%instid(inst)%ps(n)
477
478 aux_vars(n) % t2m = iv%instid(inst)%t2m(n)
479 aux_vars(n) % q2m = iv%instid(inst)%mr2m(n)
480 aux_vars(n) % u10 = iv%instid(inst)%u10(n)
481 aux_vars(n) % v10 = iv%instid(inst)%v10(n)
482 aux_vars(n) % surftype = iv%instid(inst)%landsea_mask(n)
483 aux_vars(n) % surft = iv%instid(inst)%ts(n)
484 aux_vars(n) % satzen = iv%instid(inst)%satzen(n)
485 aux_vars(n) % satazi = iv%instid(inst)%satazi(n)
486
487 ! [1.3] Call RTM forward model
488
489 ! da_rttov_direct nominally an array version, but doesn't handle arrays
490 ! of surface flags properly
491 allocate(temp(nchan))
492 call da_rttov_direct (inst, iv%instid(inst)%isflg(n), nchan, 1, &
493 con_vars(n:n), aux_vars(n:n), &
494 ! iv%instid(inst)%tb_xb(:,n), pemsv, pemsh, &
495 temp(:), pemsv((n-1)*nchan+1:n*nchan), pemsh((n-1)*nchan+1:n*nchan), &
496 iv%instid(inst)%emiss(:,n))
497 iv%instid(inst)%tb_xb(:,n)=temp(:)
498 deallocate(temp)
499 do k=1,nchan
500 ! if (iv%instid(inst)%tb_inv(k,n) > 15.0) then
501 ! write (0,*) "n,k,tb,xb",n,k, ob%instid(inst)%tb(k,n),iv%instid(inst)%tb_xb(k,n)
502 ! end if
503 end do
504 end do
505
506 if (iv%instid(inst)%nchannels > nchan) then
507 iv%instid(inst)%emiss(nchan+1:iv%instid(inst)%nchannels,n1:n2) = 0.0
508 end if
509
510 if (coefs(inst)%id_sensor == 2) then ! microwave
511 do n=n1,n2
512 iv%instid(inst)%emiss(1:nchan,n) = pems5(1+(n-n1)*nchan:(n-n1+1)*nchan)
513 end do
514 end if
515
516 do n=n1,n2
517 deallocate (con_vars(n) % t)
518 deallocate (con_vars(n) % q)
519 end do
520
521 !----------------------------------------------------------------
522 ! [2.0] calculate components of innovation vector:
523 !----------------------------------------------------------------
524
525 where (iv%instid(inst)%tb_inv(1:nchan,n1:n2) > missing_r)
526 iv%instid(inst)%tb_inv(1:nchan,n1:n2) = &
527 ob%instid(inst)%tb(1:nchan,n1:n2) - iv%instid(inst)%tb_xb(1:nchan,n1:n2)
528 elsewhere
529 iv%instid(inst)%tb_inv(1:nchan,n1:n2) = missing_r
530 endwhere
531 !do n=n1,n2
532 ! do k=1,nchan
533 ! if (iv%instid(inst)%tb_inv(k,n) > 15.0) then
534 ! 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)
535 ! end if
536 ! end do
537 !end do
538
539 deallocate (v_p)
540 deallocate (vtm)
541 deallocate (clw)
542 deallocate (dlnpf)
543 deallocate (pf)
544 deallocate (pres)
545 deallocate (con_vars)
546 deallocate (aux_vars)
547 deallocate (zasat)
548 deallocate (zlsat)
549 deallocate (polar)
550 deallocate (knchan)
551 deallocate (kprof)
552 deallocate (kchan)
553 deallocate (kochan)
554 deallocate (indexn)
555 deallocate (pemsv)
556 deallocate (pemsh)
557 deallocate (pems5)
558 end do ! end loop for sensor
559
560
561 !------------------------------------------------------------------------
562 ! [3.0] Perform bias correction
563 !------------------------------------------------------------------------
564
565 if (biascorr) then
566 do inst = 1, iv%num_inst
567 write(unit=stdout,fmt='(A,A)') 'Performing bias correction for ', &
568 trim(iv%instid(inst)%rttovid_string)
569 call da_biascorr(inst,ob,iv)
570 end do
571 end if
572
573
574 !------------------------------------------------------------------------
575 ! [4.0] Perform QC check
576 !------------------------------------------------------------------------
577
578 if (qc_rad) then
579 call da_qc_rad(ob, iv)
580 end if
581
582 do inst = 1, iv%num_inst
583 if (iv%instid(inst)%num_rad < 1) cycle
584 call da_togrid_new (iv%instid(inst)%zk(:,n1:n2), xp%kts, xp%kte, &
585 iv%instid(inst)%loc_k(:,n1:n2), &
586 iv%instid(inst)%loc_dz(:,n1:n2), &
587 iv%instid(inst)%loc_dzm(:,n1:n2))
588 end do
589
590 !------------------------------------------
591 ! 5.0 prepare bias statistics files
592 !------------------------------------------
593
594 if (biasprep) then
595 do inst = 1, iv%num_inst
596 write(unit=stdout,fmt='(A,A)') 'Preparing bias statistics files for ', &
597 trim(iv%instid(inst)%rttovid_string)
598 call da_biasprep(inst,ob,iv)
599 end do
600 end if
601
602 if (print_detail_radiance) call da_status_rad(iv,__FILE__,__LINE__)
603
604 if (trace_use) call da_trace_exit("da_get_innov_vector_rttov")
605 #endif
606
607 end subroutine da_get_innov_vector_rttov
608
609