da_get_innov_vector_rttov.inc
References to this file elsewhere.
1 subroutine da_get_innov_vector_rttov (it, grid, ob, iv)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Calculate innovation vector for radiance data.
5 !
6 ! METHOD: d = y - H(x)
7 ! 1. interpolate grid%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 (domain), intent(in) :: grid ! first guess state.
18 type (y_type), intent(inout) :: ob ! Observation structure.
19 type (iv_type), intent(inout) :: iv ! O-B structure.
20
21 #if defined(RTTOV)
22
23 integer :: n ! Loop counter.
24 integer :: k ! Index dimension.
25 integer :: nlevels ! Number of obs levels.
26
27 real :: seap, icep, lndp, snop
28 real, allocatable :: v_p(:,:)
29
30 integer :: inst, nchan
31 real, allocatable :: polar(:)
32 real, allocatable :: pres(:)
33 logical :: amsua, amsub, ssmi
34 ! logical :: mhs
35
36 integer :: kn,knchpf,sensor_type,n1,n2,n1n2
37 integer, allocatable :: kprof(:), kchan(:),kochan(:),knchan(:),indexn(:,:)
38 real , allocatable :: pemsv(:), pemsh(:), pems5(:)
39 real , allocatable :: zasat(:), zlsat(:)
40
41 ! FIX?
42 real,allocatable :: temp(:)
43
44 real, allocatable :: tgrn(:)
45 real, allocatable :: tsk(:)
46
47 type(con_vars_type), allocatable :: con_vars(:)
48 type(aux_vars_type), allocatable :: aux_vars(:)
49
50 ! variables for computing clwp
51 real, allocatable :: dpf(:,:), clw(:,:), pf(:,:)
52
53 if (trace_use) call da_trace_entry("da_get_innov_vector_rttov")
54
55 !------------------------------------------------------
56 ! [1.0] calculate the background bright temperature
57 !-------------------------------------------------------
58
59 do inst = 1, iv%num_inst ! loop for sensor
60 if ( iv%instid(inst)%num_rad < 1 ) cycle
61 nlevels = iv%instid(inst)%nlevels
62 nchan = iv%instid(inst)%nchan
63
64 amsua = trim( inst_name(rtminit_sensor(inst))) == 'amsua'
65 amsub = trim( inst_name(rtminit_sensor(inst))) == 'amsub'
66 ! mhs = trim( inst_name(rtminit_sensor(inst))) == 'mhs'
67 ssmi = trim( inst_name(rtminit_sensor(inst))) == 'ssmi'
68 sensor_type = coefs(inst)%id_sensor - 1 ! 0:IR 1:MW
69
70 if (iv%instid(inst)%info%n2 <= iv%instid(inst)%info%n1) cycle
71 n1 = iv%instid(inst)%info%n1
72 n2 = iv%instid(inst)%info%n2
73
74 n1n2=n2-n1+1
75
76 allocate (pres(1:nlevels))
77 allocate (con_vars(n1:n2))
78 allocate (aux_vars(n1:n2))
79 allocate (zasat(n1:n2))
80 allocate (zlsat(n1:n2))
81 allocate (polar(nchan))
82
83 knchpf = n1n2*nchan
84 ! knchpf = 1*nchan
85
86 allocate (knchan (n1:n2))
87 ! allocate (knchan (1))
88 allocate (kprof (knchpf))
89 allocate (kchan (knchpf))
90 allocate (kochan (knchpf))
91 allocate (indexn (1:nchan,n1:n2))
92 allocate (pemsv (knchpf))
93 allocate (pemsh (knchpf))
94 allocate (pems5 (knchpf))
95
96 pres(1:nlevels) = coefs(inst) % ref_prfl_p(1:nlevels)
97
98 allocate(v_p(kms:kme,n1:n2))
99 v_p(:,:)=0.0
100 allocate(clw(kms:kme,n1:n2))
101 allocate(dpf(kms:kme,n1:n2))
102 allocate(pf(kms:kme+1,n1:n2))
103
104 ! horizontal interpolate grid%xb pressure to ob position for every grid%xb levels
105 do n=n1,n2
106 do k=kts,kte
107 ! convert to mb
108 v_p(k,n) = 0.01*(iv%instid(inst)%info%dym(k,n)*( &
109 iv%instid(inst)%info%dxm(k,n)*grid%xb%p(iv%instid(inst)%info%i(k,n), iv%instid(inst)%info%j(k,n),k) + &
110 iv%instid(inst)%info%dx(k,n) *grid%xb%p(iv%instid(inst)%info%i(k,n)+1,iv%instid(inst)%info%j(k,n),k)) + &
111 iv%instid(inst)%info%dy(k,n) *( &
112 iv%instid(inst)%info%dxm(k,n)*grid%xb%p(iv%instid(inst)%info%i(k,n), iv%instid(inst)%info%j(k,n)+1,k) + &
113 iv%instid(inst)%info%dx(k,n) *grid%xb%p(iv%instid(inst)%info%i(k,n)+1,iv%instid(inst)%info%j(k,n)+1,k)))
114 end do
115 end do
116
117 call da_to_zk_new(pres, v_p(:,n1:n2), v_interp_p, n1n2,nlevels,iv%instid(inst)%info%zk(:,n1:n2))
118
119 call da_convert_zk (iv%instid(inst)%info)
120
121 ! [1.2] Interpolate horizontally to ob:
122 call da_interp_lin_3d (grid%xb%t, iv%instid(inst)%info, iv%instid(inst)%t (:,n1:n2))
123 call da_interp_lin_3d (grid%xb%q, iv%instid(inst)%info, iv%instid(inst)%mr(:,n1:n2))
124
125 do n= n1,n2
126 do k=1, nlevels
127 if (iv%instid(inst)%info%zk(k,n) <= 0.0) then
128 iv%instid(inst)%t(k,n) = coefs(inst) % ref_prfl_t(k,gas_id_watervapour) ! outside model level
129 iv%instid(inst)%mr(k,n) = coefs(inst) % ref_prfl_mr(k,gas_id_watervapour)
130 else
131 iv%instid(inst)%mr(k,n) = iv%instid(inst)%mr(k,n) * q2ppmv
132 ! call rttov_q2v (&
133 ! gas_unit_specconc, &! in
134 ! iv%instid(inst)%q(k,n), &! in
135 ! gas_id_watervapour, &! in
136 ! iv%instid(inst)%q(k,n), &! in
137 ! iv%instid(inst)%mr(k,n) ) ! inout
138 end if
139 if (pres(k) < 100.0) iv%instid(inst)%mr(k,n) = coefs(inst) % ref_prfl_mr(k,gas_id_watervapour)
140 end do
141
142 ! determine surface type of obs location
143 !-----------------------------------------
144 call da_detsurtyp( grid%xb%snow, grid%xb%xice, grid%xb%landmask, &
145 grid%xb%ivgtyp, grid%xb%isltyp, &
146 ims, ime, jms, jme, &
147 iv%instid(inst)%info%i(1,n), iv%instid(inst)%info%j(1,n), &
148 iv%instid(inst)%info%dx(1,n), iv%instid(inst)%info%dy(1,n), &
149 iv%instid(inst)%info%dxm(1,n), iv%instid(inst)%info%dym(1,n), &
150 iv%instid(inst)%isflg(n),iv%instid(inst)%vegtyp(n), iv%instid(inst)%soiltyp(n), &
151 seap, icep, lndp, snop )
152 end do
153
154 call da_interp_lin_2d (grid%xb % u10, iv%instid(inst)%info, 1,iv%instid(inst)%u10(n1:n2))
155 call da_interp_lin_2d (grid%xb % v10, iv%instid(inst)%info, 1,iv%instid(inst)%v10(n1:n2))
156 call da_interp_lin_2d (grid%xb % t2, iv%instid(inst)%info, 1,iv%instid(inst)%t2m(n1:n2))
157 call da_interp_lin_2d (grid%xb % q2, iv%instid(inst)%info, 1,iv%instid(inst)%q2m(n1:n2))
158 call da_interp_lin_2d (grid%xb % psfc, iv%instid(inst)%info, 1,iv%instid(inst)%ps (n1:n2))
159
160 allocate (tgrn(n1:n2))
161 allocate (tsk(n1:n2))
162 call da_interp_lin_2d (grid%xb % tgrn, iv%instid(inst)%info, 1, tgrn)
163 call da_interp_lin_2d (grid%xb % tsk, iv%instid(inst)%info, 1, tsk)
164
165 do n=n1,n2
166 if (iv%instid(inst)%isflg(n) == 0) then ! over sea using SST
167 iv%instid(inst)%ts(n) = tgrn(n)
168 else
169 iv%instid(inst)%ts(n) = tsk(n)
170 end if
171 end do
172 deallocate (tgrn)
173 deallocate (tsk)
174
175 call da_interp_lin_2d (grid%xb%terr, iv%instid(inst)%info, 1, iv%instid(inst)%elevation(n1:n2))
176
177 ! variables for emissivity calculations
178 !----------------------------------------
179
180 call da_interp_lin_2d (grid%xb%smois, iv%instid(inst)%info, 1, iv%instid(inst)%smois(n1:n2))
181 call da_interp_lin_2d (grid%xb%tslb, iv%instid(inst)%info, 1, iv%instid(inst)%tslb(n1:n2))
182 call da_interp_lin_2d (grid%xb%snowh, iv%instid(inst)%info, 1, iv%instid(inst)%snowh(n1:n2))
183 call da_interp_lin_2d (grid%xb%vegfra, iv%instid(inst)%info, 1, iv%instid(inst)%vegfra(n1:n2))
184
185 ! iv%instid(inst)%snowh(n1:n2) = iv%instid(inst)%snowh(n1:n2)*100.0 ! convert from m to mm
186 iv%instid(inst)%vegfra(n1:n2) = 0.01*iv%instid(inst)%vegfra(n1:n2) ! convert range to 0~1
187
188 ! variables for calculation of cloud affected radiance
189 !-------------------------------------------------------
190 do k=kts,kte
191 call da_interp_lin_2d (grid%xb%t (:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%tm(k,:))
192 call da_interp_lin_2d (grid%xb%q (:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%qm(k,:))
193 call da_interp_lin_2d (grid%xb%qrn(:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%qrn(k,:))
194 call da_interp_lin_2d (grid%xb%qcw(:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%qcw(k,:))
195 call da_interp_lin_2d (grid%xb%qci(:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%qci(k,:))
196 call da_interp_lin_2d (grid%xb%qsn(:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%qsn(k,:))
197 call da_interp_lin_2d (grid%xb%qgr(:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%qgr(k,:))
198 end do
199
200 iv%instid(inst)%pm(:,n1:n2) = v_p(:,n1:n2)
201
202 ! convert to hPa
203 iv%instid(inst)%ps(n1:n2) = 0.01 * iv%instid(inst)%ps(n1:n2)
204 iv%instid(inst)%mr2m(n1:n2) = iv%instid(inst)%q2m(n1:n2) * q2ppmv
205
206 ! ADD for computing cloud liquid water path (mm) from guess
207
208 pf(kts,n1:n2) = 100.0*iv%instid(inst)%ps(n1:n2)
209 do k = 2,kte
210 pf(k,n1:n2) = 50.0*(v_p(k-1,n1:n2)+v_p(k,n1:n2))
211 end do
212 pf(kte+1,n1:n2)= 50.0*v_p(kte,n1:n2)
213
214 iv%instid(inst)%clwp(n1:n2) = 0.0
215 do k = kts,kte
216 dpf(k,n1:n2) = pf(k,n1:n2) - pf(k+1,n1:n2)
217 clw(k,n1:n2) = iv%instid(inst)%qcw(k,n1:n2)*dpf(k,n1:n2)/gravity
218 where (v_p(k,n1:n2)<100.0) clw (k,n1:n2) = 0.0
219 iv%instid(inst)%clwp(n1:n2) =iv%instid(inst)%clwp(n1:n2) + clw(k,n1:n2)
220 end do
221
222 ! call rttov_q2v (&
223 ! gas_unit_specconc, &! in
224 ! iv%instid(inst)%q2m(n1:n2), &! in
225 ! gas_id_watervapour, &! in
226 ! iv%instid(inst)%q2m(n1:n2), &! in
227 ! iv%instid(inst)%mr2m(n1:n2) ! inout
228
229 ! Compute surface emissivity for observation
230 !-------------------------------------------
231
232 zasat(n1:n2) = iv%instid(inst)%satzen(n1:n2) * deg2rad
233 zlsat(n1:n2) = iv%instid(inst)%satazi(n1:n2) * deg2rad
234
235 ! do i = n1,n2
236 ! knchan(i) = nchan
237 ! do kn=1,nchan
238 ! kprof((i-1)*nchan+kn) = i
239 ! kchan((i-1)*nchan+kn) = kn
240 ! kochan((i-1)*nchan+kn) = kn
241 ! indexn(kn,i) = (-1)*nchan+kn
242 ! end do
243 ! end do
244
245 pemsv(:) = 0.0
246 pemsh(:) = 0.0
247
248 if (amsua) then
249 polar(1:4) = 0.0 ! polarizition: 0:vertical; 1:horizontal
250 polar(5:6) = 1.0
251 polar(7) = 0.0
252 polar(8:14)= 1.0
253 polar(15) = 0.0
254 end if
255 if (amsub) then
256 polar(1:5) = 0.0
257 end if
258
259 do n=n1,n2
260 knchan(n) = nchan
261 do kn=1,nchan
262 kprof((n-1)*nchan+kn) = n
263 kchan((n-1)*nchan+kn) = kn
264 kochan((n-1)*nchan+kn) = kn
265 indexn(kn,n) = (n-1)*nchan+kn
266 end do
267 end do
268
269 if (coefs(inst)%id_sensor == 2) then ! microwave
270 call da_trace("da_get_innov_vector_rttov",message="Calling gsi_emiss")
271 call gsi_emiss(inst,knchpf,kprof,kchan,knchan,indexn, &
272 iv%instid(inst)%satzen * deg2rad, &
273 iv%instid(inst)%satazi * deg2rad, polar, &
274 iv%instid(inst)%isflg,pemsv,pemsh,pems5, &
275 iv%instid(inst)%ts, &
276 iv%instid(inst)%soiltyp, &
277 iv%instid(inst)%tslb, &
278 iv%instid(inst)%smois, &
279 iv%instid(inst)%vegtyp, &
280 iv%instid(inst)%vegfra, &
281 iv%instid(inst)%snowh, &
282 sensor_type,n1n2,nchan, &
283 iv%instid(inst)%tb_inv(:,:),amsua,amsub,ssmi, &
284 iv%instid(inst)%u10, &
285 iv%instid(inst)%v10)
286 end if
287
288 do n=n1,n2
289 con_vars(n) % nlevels = nlevels
290 allocate (con_vars(n) % t(nlevels))
291 allocate (con_vars(n) % q(nlevels))
292
293 con_vars(n) % t(1:nlevels) = iv%instid(inst)%t(1:nlevels,n)
294 con_vars(n) % q(1:nlevels) = iv%instid(inst)%mr(1:nlevels,n)
295 con_vars(n) % ps = iv%instid(inst)%ps(n)
296
297 aux_vars(n) % t2m = iv%instid(inst)%t2m(n)
298 aux_vars(n) % q2m = iv%instid(inst)%mr2m(n)
299 aux_vars(n) % u10 = iv%instid(inst)%u10(n)
300 aux_vars(n) % v10 = iv%instid(inst)%v10(n)
301 aux_vars(n) % surftype = iv%instid(inst)%landsea_mask(n)
302 aux_vars(n) % surft = iv%instid(inst)%ts(n)
303 aux_vars(n) % satzen = iv%instid(inst)%satzen(n)
304 aux_vars(n) % satazi = iv%instid(inst)%satazi(n)
305
306 ! [1.3] Call RTM forward model
307
308 ! da_rttov_direct nominally an array version, but doesn't handle arrays
309 ! of surface flags properly
310 allocate(temp(nchan))
311 call da_rttov_direct (inst, iv%instid(inst)%isflg(n), nchan, 1, &
312 con_vars(n:n), aux_vars(n:n), &
313 ! iv%instid(inst)%tb_xb(:,n), pemsv, pemsh, &
314 temp(:), pemsv((n-1)*nchan+1:n*nchan), pemsh((n-1)*nchan+1:n*nchan), &
315 iv%instid(inst)%emiss(:,n))
316 iv%instid(inst)%tb_xb(:,n)=temp(:)
317 deallocate(temp)
318 end do
319
320 if (iv%instid(inst)%nchannels > nchan) then
321 iv%instid(inst)%emiss(nchan+1:iv%instid(inst)%nchannels,n1:n2) = 0.0
322 end if
323
324 if (coefs(inst)%id_sensor == 2) then ! microwave
325 do n=n1,n2
326 iv%instid(inst)%emiss(1:nchan,n) = pems5(1+(n-n1)*nchan:(n-n1+1)*nchan)
327 end do
328 end if
329
330 do n=n1,n2
331 deallocate (con_vars(n) % t)
332 deallocate (con_vars(n) % q)
333 end do
334
335 !----------------------------------------------------------------
336 ! [2.0] calculate components of innovation vector:
337 !----------------------------------------------------------------
338
339 do n=n1,n2
340 do k=1,nchan
341 if (iv%instid(inst)%tb_inv(k,n) > missing_r) then
342 iv%instid(inst)%tb_inv(k,n) = ob%instid(inst)%tb(k,n) - iv%instid(inst)%tb_xb(k,n)
343 else
344 iv%instid(inst)%tb_inv(k,n) = missing_r
345 end if
346 end do
347 end do
348
349 deallocate (v_p)
350 deallocate (clw)
351 deallocate (dpf)
352 deallocate (pf)
353 deallocate (pres)
354 deallocate (con_vars)
355 deallocate (aux_vars)
356 deallocate (zasat)
357 deallocate (zlsat)
358 deallocate (polar)
359 deallocate (knchan)
360 deallocate (kprof)
361 deallocate (kchan)
362 deallocate (kochan)
363 deallocate (indexn)
364 deallocate (pemsv)
365 deallocate (pemsh)
366 deallocate (pems5)
367
368 !------------------------------------------------------------------------
369 ! [3.0] Perform bias correction
370 !------------------------------------------------------------------------
371
372 if (biascorr) then
373 write(unit=stdout,fmt='(A,A)') 'Performing bias correction for ', &
374 trim(iv%instid(inst)%rttovid_string)
375 call da_biascorr(inst,ob,iv)
376 end if
377
378 call da_convert_zk (iv%instid(inst)%info)
379 end do ! end loop for sensor
380
381 !------------------------------------------------------------------------
382 ! [4.0] Perform QC check
383 !------------------------------------------------------------------------
384
385 if (qc_rad) then
386 call da_qc_rad(ob, iv)
387 end if
388
389 !------------------------------------------
390 ! 5.0 prepare bias statistics files
391 !------------------------------------------
392
393 if (biasprep) then
394 do inst = 1, iv%num_inst
395 call da_biasprep(inst,ob,iv)
396 end do
397 end if
398
399 if (trace_use) call da_trace_exit("da_get_innov_vector_rttov")
400 #else
401 call da_error(__FILE__,__LINE__, &
402 (/"Must compile with $RTTOV option for radiances"/))
403 #endif
404
405 end subroutine da_get_innov_vector_rttov
406
407