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