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