da_get_innov_vector_crtm.inc
References to this file elsewhere.
1 subroutine da_get_innov_vector_crtm ( 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 foreward RTM to get simulated bright temperature
9 ! 3. obs BT - simulated BT
10 !
11 ! HISTORY: 11/14/2005 - Creation Zhiquan Liu
12 !---------------------------------------------------------------------------
13
14 IMPLICIT NONE
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, icld ! Loop counter.
23 INTEGER :: i, j, k ! Index dimension.
24 INTEGER :: num_levs ! Number of obs levels.
25 REAL :: dx, dxm ! Interpolation weights.
26 REAL :: dy, dym ! Interpolation weights.
27 INTEGER :: alloc_status(40)
28
29 REAL :: model_u10, model_v10
30 REAL :: model_psfc, model_ptop
31 REAL :: model_ts, model_elv
32 REAL :: model_smois, model_tslb, model_snowh
33 REAL :: model_isltyp, model_ivgtyp, model_vegfra
34 INTEGER :: model_isflg
35
36 REAL, DIMENSION(xp%kms:xp%kme) :: v_p, & ! Model value p at ob hor. location.
37 model_qcw
38
39 ! REAL, DIMENSION(xp%kms:xp%kme) :: model_tm
40
41 INTEGER :: inst, nchanl
42 INTEGER :: nchan_emis
43
44 ! variables for computing clwp
45 REAL, dimension(xp%kms:xp%kme) :: vtm, clw, dlnpf
46 REAL :: clwp
47
48 #if defined(CRTM)
49 ! CRTM local varaibles and types
50 INTEGER :: wmo_sensor_id,Error_Status, Allocate_Status
51 TYPE( CRTM_RTSolution_type ), ALLOCATABLE :: RTSolution(:)
52 TYPE (CRTM_ChannelInfo_type) :: ChannelInfo
53 TYPE( CRTM_Atmosphere_type ) :: Atmosphere
54 TYPE( CRTM_Surface_type ) :: Surface
55 TYPE( CRTM_GeometryInfo_type ) :: GeometryInfo
56 #endif
57 !---------------------------------------------------------
58
59 !JRB use argument it
60 if (it==0) then; print *,"why have argument it to"//__FILE__; end if
61
62 alloc_status (:) = 0
63
64 IF ( iv%num_inst < 1 ) return
65
66 #if !defined(CRTM)
67 call da_error(__FILE__,__LINE__, &
68 (/"Must compile with $CRTM option for radiances"/))
69 #else
70
71 if (trace_use) call da_trace_entry("da_get_innov_vector_crtm")
72
73 !----------------------------------------------------------------------------
74 ! CRTM allocation
75 !
76 ! Atmosphere structure
77 Atmosphere%n_Layers=(xp%kte-xp%kts)+1 ! number of vertical levels
78 Atmosphere%n_Absorbers=2
79 Atmosphere%n_Clouds=0
80 Atmosphere%n_Aerosols=0
81 if (crtm_cloud) Atmosphere%n_Clouds=6
82
83 Error_Status = CRTM_Allocate_Atmosphere( Atmosphere%n_Layers, &
84 Atmosphere%n_Absorbers, &
85 Atmosphere%n_Clouds, &
86 Atmosphere%n_Aerosols, &
87 Atmosphere)
88 if ( Error_Status /= 0 ) THEN
89 call da_error(__FILE__,__LINE__, &
90 (/"Error in allocatting CRTM Atmosphere Structure"/))
91 endif
92
93 Atmosphere%Absorber_ID(1)=H2O_ID
94 Atmosphere%Absorber_ID(2)=O3_ID
95
96 if (crtm_cloud) then
97 Atmosphere%Cloud(1)%Type=WATER_CLOUD
98 Atmosphere%Cloud(2)%Type=ICE_CLOUD
99 Atmosphere%Cloud(3)%Type=RAIN_CLOUD
100 Atmosphere%Cloud(4)%Type=SNOW_CLOUD
101 Atmosphere%Cloud(5)%Type=GRAUPEL_CLOUD
102 Atmosphere%Cloud(6)%Type=HAIL_CLOUD
103 end if
104
105 !------------------------------------------------------
106 ! [1.0] calculate the background bright temperature
107 !-------------------------------------------------------
108 do inst = 1, iv%num_inst ! loop for sensor
109 ! if ( iv%instid(inst)%num_rad < 1 ) cycle
110 if (iv%ob_numb(iv%current_ob_time)%radiance(inst) <= &
111 iv%ob_numb(iv%current_ob_time-1)%radiance(inst) ) cycle
112 write(UNIT=stdout,FMT='(A,A)') ' Computing Innovation : ',iv%instid(inst)%rttovid_string
113 num_levs = xp%kte-xp%kts+1
114 ! CRTM channel information structure
115 Error_Status = CRTM_Set_ChannelInfo(Sensor_Descriptor(inst),ChannelInfo)
116 if ( Error_Status /= 0) then
117 call da_error(__FILE__,__LINE__, &
118 (/"Error in calling CRTM_Set_ChannelInfo"/))
119 endif
120 nchanl = ChannelInfo%n_channels
121
122 ! Allocate forward model solution RTSolution array to number of channels
123 ALLOCATE( RTSolution( ChannelInfo%n_Channels ), &
124 STAT = Allocate_Status )
125 IF ( Allocate_Status /= 0 ) THEN
126 call da_error(__FILE__,__LINE__, &
127 (/"Error in allocatting RTSolution"/))
128 END IF
129
130 ! CRTM Surface Structure
131 if (trim(crtm_sensor_name(rtminit_sensor(inst))) =='amsua') then
132 nchan_emis=4
133 wmo_sensor_id=WMO_AMSUA
134 elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='amsub') then
135 nchan_emis=2
136 wmo_sensor_id=WMO_AMSUB
137 elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='amsre') then
138 nchan_emis=12
139 wmo_sensor_id=WMO_AMSRE
140 elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='ssmi') then
141 nchan_emis=7
142 wmo_sensor_id=WMO_SSMI
143 else
144 nchan_emis=0
145 wmo_sensor_id=INVALID_WMO_SENSOR_ID
146 endif
147
148 Error_Status = CRTM_Allocate_Surface( nchan_emis, & ! Input
149 Surface) ! Output
150 if ( Error_Status /= 0 ) THEN
151 call da_error(__FILE__,__LINE__, &
152 (/"Error in allocatting CRTM Surface Structure Structure"/))
153 endif
154
155 if (iv%ob_numb(iv%current_ob_time)%radiance(inst) <= &
156 iv%ob_numb(iv%current_ob_time-1)%radiance(inst) ) cycle
157
158 ! do n= 1, iv%instid(inst)%num_rad ! loop for pixel
159 do n= iv%ob_numb(iv%current_ob_time-1)%radiance(inst)+1, &
160 iv%ob_numb(iv%current_ob_time)%radiance(inst)
161 ! if ( n > iv%instid(inst)%num_rad ) exit
162 model_ts = 0.0
163 model_u10 = 0.0
164 model_v10 = 0.0
165 model_psfc = 0.0
166
167 ! [1.1] Get horizontal interpolation weights:
168
169 i = iv%instid(inst)%loc_i(n)
170 j = iv%instid(inst)%loc_j(n)
171 dx = iv%instid(inst)%loc_dx(n)
172 dy = iv%instid(inst)%loc_dy(n)
173 dxm = iv%instid(inst)%loc_dxm(n)
174 dym = iv%instid(inst)%loc_dym(n)
175 ! horizontal interpolate xb pressure to ob position for every xb layer
176 ! get CRTM pressure Layers
177 do k=xp%kts,xp%kte ! from bottem to top
178 v_p(k) = dym*(dxm*xb%p(i,j ,k) + dx*xb%p(i+1,j ,k)) &
179 + dy *(dxm*xb%p(i,j+1,k) + dx*xb%p(i+1,j+1,k))
180 v_p(k) = 0.01*v_p(k) ! convert Pa to hPa
181 Atmosphere%Pressure(xp%kte-k+1)=v_p(k) ! from top to bottom
182 enddo
183
184 ! [1.2] Interpolate horizontally to ob:
185 do k=xp%kts,xp%kte ! from bottem to top
186 call da_interp_lin_2d( xb%t(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
187 i, j, dx, dy, dxm, dym, & ! temperature (in K)
188 Atmosphere%Temperature(xp%kte-k+1) )
189 call da_interp_lin_2d( xb%q(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
190 i, j, dx, dy, dxm, dym, & ! mixture ratio (in kg/kg)
191 Atmosphere%Absorber(xp%kte-k+1,1) )
192 Atmosphere%Absorber(xp%kte-k+1,1) = 1000.*Atmosphere%Absorber(xp%kte-k+1,1) ! in g/kg
193 ! NOTE: WRF high-level q values seems too big, replaced by constants
194 if ( v_p(k) < 75. ) Atmosphere%Absorber(xp%kte-k+1,1) = 0.001
195
196 call da_interp_lin_2d( xb%qcw(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
197 i, j, dx, dy, dxm, dym, & ! cloud water
198 model_qcw(xp%kte-k+1) )
199
200 if (crtm_cloud) then
201 Atmosphere%Cloud(1)%Effective_Radius(xp%kte-k+1)=10.
202 Atmosphere%Cloud(1)%Water_Content(xp%kte-k+1) = &
203 model_qcw(xp%kte-k+1)
204
205 Atmosphere%Cloud(2)%Effective_Radius(xp%kte-k+1)=200.
206 call da_interp_lin_2d( xb%qci(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
207 i, j, dx, dy, dxm, dym, & ! cloud ice
208 Atmosphere%Cloud(2)%Water_Content(xp%kte-k+1) )
209
210 Atmosphere%Cloud(3)%Effective_Radius(xp%kte-k+1)=200.
211 call da_interp_lin_2d( xb%qrn(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
212 i, j, dx, dy, dxm, dym, & ! rain
213 Atmosphere%Cloud(3)%Water_Content(xp%kte-k+1) )
214
215 Atmosphere%Cloud(4)%Effective_Radius(xp%kte-k+1)=200.
216 call da_interp_lin_2d( xb%qsn(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
217 i, j, dx, dy, dxm, dym, & ! snow
218 Atmosphere%Cloud(4)%Water_Content(xp%kte-k+1) )
219
220 Atmosphere%Cloud(5)%Effective_Radius(xp%kte-k+1)=200.
221 call da_interp_lin_2d( xb%qgr(:,:,k), xp%ims, xp%ime, xp%jms, xp%jme, &
222 i, j, dx, dy, dxm, dym, & ! grauple
223 Atmosphere%Cloud(5)%Water_Content(xp%kte-k+1) )
224
225 Atmosphere%Cloud(6)%Effective_Radius(xp%kte-k+1)=200.
226 Atmosphere%Cloud(6)%Water_Content(xp%kte-k+1) = 0.
227
228 end if
229 enddo
230
231 ! determine surface type of obs location
232 !-----------------------------------------
233 call da_detsurtyp ( xb%snow, xb%xice, xb%landmask, &
234 xb%ivgtyp, xb%isltyp, &
235 xp%ims, xp%ime, xp%jms, xp%jme, &
236 i, j, dx, dy, dxm, dym, &
237 model_isflg,model_ivgtyp, model_isltyp, &
238 Surface%Water_Coverage, Surface%Ice_Coverage, &
239 Surface%Land_Coverage, Surface%Snow_Coverage )
240
241 call da_interp_lin_2d( xb % u10, xp%ims, xp%ime, xp%jms, xp%jme, &
242 i, j, dx, dy, dxm, dym, &
243 model_u10 )
244 call da_interp_lin_2d( xb % v10, xp%ims, xp%ime, xp%jms, xp%jme, &
245 i, j, dx, dy, dxm, dym, &
246 model_v10 )
247 call da_interp_lin_2d( xb % psfc, xp%ims, xp%ime, xp%jms, xp%jme, &
248 i, j, dx, dy, dxm, dym, &
249 model_psfc )
250
251 model_psfc = 0.01*model_psfc ! convert to hPa
252 model_ptop = 0.01*xb%ptop
253
254 ! get CRTM levels (0.005hPa at top) /model full level
255 Atmosphere%Level_Pressure(0)=model_ptop ! to sigma level 51-->sigmaf=0
256 Atmosphere%Level_Pressure(Atmosphere%n_Layers)=model_psfc ! to sigma level 1->sigmaf=1
257 do k=xp%kts+1,xp%kte
258 Atmosphere%Level_Pressure(xp%kte-k+1)= &
259 xb%sigmaf(k)*(model_psfc-model_ptop)+model_ptop
260 enddo
261
262 ! convert cloud content unit from kg/kg to kg/m^2
263 if (crtm_cloud) then
264 do k=xp%kts,xp%kte
265 do icld=1,Atmosphere%n_Clouds
266 Atmosphere%Cloud(icld)%Water_Content(k)= Atmosphere%Cloud(icld)%Water_Content(k)* &
267 (Atmosphere%Level_Pressure(k)- Atmosphere%Level_Pressure(k-1))*100./gravity
268 enddo
269 enddo
270 end if
271
272 if ( model_isflg == 0 ) then ! over sea using SST
273 call da_interp_lin_2d( xb % tgrn, xp%ims, xp%ime, xp%jms, xp%jme, &
274 i, j, dx, dy, dxm, dym, &
275 model_ts )
276 else
277 call da_interp_lin_2d( xb % tsk, xp%ims, xp%ime, xp%jms, xp%jme, &
278 i, j, dx, dy, dxm, dym, &
279 model_ts )
280 end if
281
282 call da_interp_lin_2d( xb % terr, xp%ims, xp%ime, xp%jms, xp%jme, &
283 i, j, dx, dy, dxm, dym, &
284 model_elv )
285
286 ! variables for emissivity calculations
287 !----------------------------------------
288 call da_interp_lin_2d( xb % smois, xp%ims, xp%ime, xp%jms, xp%jme, &
289 i, j, dx, dy, dxm, dym, &
290 model_smois )
291 call da_interp_lin_2d( xb % tslb, xp%ims, xp%ime, xp%jms, xp%jme, &
292 i, j, dx, dy, dxm, dym, &
293 model_tslb )
294 call da_interp_lin_2d( xb % snowh, xp%ims, xp%ime, xp%jms, xp%jme, &
295 i, j, dx, dy, dxm, dym, &
296 model_snowh )
297 call da_interp_lin_2d( xb % vegfra, xp%ims, xp%ime, xp%jms, xp%jme, &
298 i, j, dx, dy, dxm, dym, &
299 model_vegfra )
300
301 ! model_snowh = model_snowh*100.0 ! convert from m to mm
302 model_vegfra = 0.01*model_vegfra ! convert range to 0~1
303
304 do k=xp%kts,xp%kte
305 ! ADD for computing cloud liquid water path
306 vtm (k) = (1.+0.608*0.001*atmosphere%absorber(k,1))*atmosphere%temperature(k) ! virtual T
307 ! convert kg/kg to g/m3
308 ! clw (k) = model_qcw(k)*v_p(k)*100.*0.0289644/8.31451/model_tm(k)
309 clw (k) = 0.34836*model_qcw(k)*atmosphere%pressure(k)/atmosphere%temperature(k)
310 if ( atmosphere%pressure(k)<100. ) clw (k) = 0.
311 enddo
312
313 ! ADD for computing cloud liquid water path (mm) from guess
314 clwp = 0.
315 do k = xp%kts,xp%kte ! from top to bottom
316 dlnpf(k) = log(atmosphere%level_pressure(k)) - log(atmosphere%level_pressure(k-1))
317 clw (k) = 29.27095*clw(k)*vtm(k)*dlnpf(k)
318 clwp = clwp + clw(k)
319 end do
320 clwp = 0.001*clwp ! kg/m2 = mm
321
322 ! CRTM GeometryInfo Structure
323 GeometryInfo%Sensor_Zenith_Angle=iv%instid(inst)%satzen(n)
324 GeometryInfo%Source_Zenith_Angle=iv%instid(inst)%solzen(n)
325 ! GeometryInfo%Satellite_Height=830.
326 ! GeometryInfo%Sensor_Scan_Angle=
327 ! GeometryInfo%Sensor_Zenith_Angle=
328 ! GeometryInfo%Sensor_Scan_Angle=
329 ! GeometryInfo%Source_Zenith_Angle=
330
331 ! CRTM Surface parameter data
332
333 if (Surface%Land_Coverage > 0.) then
334 Surface%Land_Type=GRASS_SOIL ! land type (User guide appendix 3)
335 Surface%Land_Temperature=model_ts ! K
336 Surface%Soil_Moisture_Content= model_smois !0.05 ! volumetric water content (g/cm**3)
337 !Surface%Canopy_Water_Content=0.05 ! gravimetric water content
338 Surface%Vegetation_Fraction=model_vegfra
339 Surface%Soil_Temperature=model_tslb
340 endif
341 if (Surface%Water_Coverage > 0.) then
342 !Surface%Water_Type=SEA_WATER ! (Currently NOT used)
343 Surface%Water_Temperature=model_ts ! K
344 Surface%Wind_Speed=sqrt(model_u10**2+model_v10**2) ! m/sec
345 !surface%Wind_Direction=0.0 ! NOT used
346 Surface%Salinity=33. ! ppmv
347 endif
348 if (Surface%Snow_Coverage > 0.) then
349 Surface%Snow_Type=NEW_SNOW ! User guide appendix 3
350 Surface%Snow_Temperature=model_ts ! K
351 Surface%Snow_Depth=model_snowh ! mm
352 !Surface%Snow_Density=0.2 ! g/cm**3
353 !Surface%Snow_Grain_Size=2.0 ! mm
354 endif
355 if (Surface%Ice_Coverage > 0.) then
356 !Surface%Ice_Type=FRESH_ICE ! NO Table offered, single example is FRESH_ICE
357 Surface%Ice_Temperature=model_ts ! K
358 Surface%Ice_Thickness=10. ! mm
359 !Surface%Ice_Density=0.9 ! g/cm**3
360 !Surface%Ice_Roughness=0. ! NO Table offered, single example is ZERO
361 endif
362 if (nchan_emis > 0) then
363 Surface%SensorData%n_channels = nchan_emis
364 Surface%SensorData%Sensor_ID = wmo_sensor_id
365 Surface%SensorData%Tb(1:nchan_emis) = ob%instid(inst)%tb(1:nchan_emis,n)
366 end if
367
368 ! [1.3] Call RTM foreward model
369 call da_crtm_direct(nchanl, Atmosphere, &
370 Surface, &
371 GeometryInfo, &
372 ChannelInfo, &
373 RTSolution)
374
375 !Error_Status = CRTM_Forward (Atmosphere, &
376 ! Surface, &
377 ! GeometryInfo, &
378 ! ChannelInfo, &
379 ! RTSolution)
380 !if (n<=10) then
381 !write(6,'(f10.3)') Atmosphere%level_pressure(0)
382 !do k=1,Atmosphere%n_layers
383 ! write(6,'(4f10.3)') Atmosphere%Level_pressure(k),Atmosphere%pressure(k), &
384 ! Atmosphere%temperature(k),Atmosphere%absorber(k,1)
385 !enddo
386 ! write(6,'(15f8.3)') RTSolution(1:nchanl)%Brightness_Temperature
387 ! write(6,'(15f8.3)') RTSolution(1:nchanl)%surface_emissivity
388 !end if
389
390 !if ( Error_Status /= 0 ) THEN
391 ! call da_error(__FILE__,__LINE__, &
392 ! (/"Error in calling CRTM_Forward"/))
393 !endif
394
395 !----------------------------------------------------------------
396 ! [2.0] calculate components of innovation vector:
397 !----------------------------------------------------------------
398 do k = 1, nchanl
399 if ( iv%instid(inst)%tb_inv(k,n) > missing_r ) then
400 iv%instid(inst)%tb_xb(k,n) = rtsolution(k)%Brightness_Temperature
401 iv%instid(inst)%tb_inv(k,n) = &
402 ob%instid(inst)%tb(k,n) - rtsolution(k)%Brightness_Temperature
403 else
404 iv%instid(inst)%tb_xb(k,n) = rtsolution(k)%Brightness_Temperature
405 iv%instid(inst)%tb_inv(k,n) = missing_r
406 end if
407 iv%instid(inst)%emiss(k,n) = rtsolution(k)%surface_emissivity
408 end do
409
410 !----------------------------------------------------------------
411 ! [3.0] store base state to innovation structure
412 !----------------------------------------------------------------
413 ! full level pressures
414 iv%instid(inst)%pf(0,n) = atmosphere%level_pressure(0)
415 do k=1,atmosphere%n_layers
416 iv%instid(inst)%pm(k,n) = atmosphere%pressure(k)
417 iv%instid(inst)%pf(k,n) = atmosphere%level_pressure(k)
418 iv%instid(inst)%tm(k,n) = atmosphere%temperature(k)
419 iv%instid(inst)%qm(k,n) = atmosphere%absorber(k,1)
420 if (crtm_cloud) then
421 iv%instid(inst)%qcw(k,n) = atmosphere%cloud(1)%water_content(k)
422 iv%instid(inst)%qci(k,n) = atmosphere%cloud(2)%water_content(k)
423 iv%instid(inst)%qrn(k,n) = atmosphere%cloud(3)%water_content(k)
424 iv%instid(inst)%qsn(k,n) = atmosphere%cloud(4)%water_content(k)
425 iv%instid(inst)%qgr(k,n) = atmosphere%cloud(5)%water_content(k)
426 end if
427 enddo
428 iv%instid(inst)%u10(n) = model_u10
429 iv%instid(inst)%v10(n) = model_v10
430 iv%instid(inst)%t2m(n) = 0.01*missing_r !model_t2m
431 iv%instid(inst)%mr2m(n) = 0.01*missing_r !model_mr2m
432 iv%instid(inst)%ps(n) = model_psfc
433 iv%instid(inst)%ts(n) = model_ts
434 iv%instid(inst)%smois(n) = model_smois
435 iv%instid(inst)%tslb(n) = model_tslb
436 iv%instid(inst)%snowh(n) = model_snowh
437 iv%instid(inst)%isflg(n) = model_isflg
438 iv%instid(inst)%elevation(n) = model_elv
439 iv%instid(inst)%soiltyp(n) = model_isltyp
440 iv%instid(inst)%vegtyp(n) = model_ivgtyp
441 iv%instid(inst)%vegfra(n) = model_vegfra
442 iv%instid(inst)%clwp(n) = clwp
443 iv%instid(inst)%water_coverage(n)=surface%water_coverage
444 iv%instid(inst)%land_coverage(n)=surface%land_coverage
445 iv%instid(inst)%ice_coverage(n)=surface%ice_coverage
446 iv%instid(inst)%snow_coverage(n)=surface%snow_coverage
447
448 end do ! end loop for pixels
449
450 deallocate( rtsolution, STAT = Allocate_Status )
451 IF ( Allocate_Status /= 0 ) THEN
452 call da_error(__FILE__,__LINE__, &
453 (/"Error in deallocatting RTSolution"/))
454 END IF
455
456 Error_Status = CRTM_Destroy_Surface(Surface)
457 if ( Error_Status /= 0 ) THEN
458 call da_error(__FILE__,__LINE__, &
459 (/"Error in deallocatting CRTM Surface Structure"/))
460 endif
461 end do ! end loop for sensor
462
463 !------------------------------------------
464 ! 4.0 perfoming bias correction files
465 !------------------------------------------
466
467 if (biascorr) then
468 do inst = 1, iv%num_inst
469 WRITE(UNIT=stdout,FMT='(A,A)') 'Performing bias correction for ', &
470 trim(iv%instid(inst)%rttovid_string)
471 call da_biascorr(inst,ob,iv)
472 end do
473 end if
474
475 !------------------------------------------------------------------------
476 ! [5.0] Perform QC check
477 !------------------------------------------------------------------------
478 if (qc_rad) then
479 call da_qc_crtm(ob, iv)
480 end if
481 !------------------------------------------
482 ! 6.0 prepare bias statistics files
483 !------------------------------------------
484
485 if (biasprep) then
486 do inst = 1, iv%num_inst
487 WRITE(UNIT=stdout,FMT='(A,A)') 'Preparing bias statistics files for ', &
488 trim(iv%instid(inst)%rttovid_string)
489 call da_biasprep(inst,ob,iv)
490 end do
491 end if
492
493 Error_Status = CRTM_Destroy_Atmosphere( Atmosphere )
494 if ( Error_Status /= 0 ) THEN
495 call da_error(__FILE__,__LINE__, &
496 (/"Error in deallocatting CRTM Atmosphere Structure"/))
497 endif
498
499 if (trace_use) call da_trace_exit("da_get_innov_vector_crtm")
500 #endif
501
502 end subroutine da_get_innov_vector_crtm
503