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