da_get_innov_vector_crtmk.inc
References to this file elsewhere.
1 subroutine da_get_innov_vector_crtmk ( 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 :: l ! 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, allocatable :: model_u10(:)
30 real, allocatable :: model_v10(:)
31 real, allocatable :: model_psfc(:)
32 real :: model_ptop
33 real, allocatable :: model_ts(:)
34 real, allocatable :: model_elv(:)
35 real, allocatable :: model_smois(:)
36 real, allocatable :: model_tslb(:)
37 real, allocatable :: model_snowh(:)
38 real, allocatable :: model_vegfra(:)
39 real :: model_isltyp, model_ivgtyp
40 integer :: model_isflg
41
42 real :: v_p(kms:kme) ! Model value p at ob hor. location.
43 real :: model_qcw(kms:kme)
44 real :: model_rho(kms:kme)
45 real, allocatable :: model_snow(:) ! snow water equivalent, different from model_snowh,
46 ! used in calculating reff_water
47
48 ! real :: model_tm(kms:kme)
49
50 integer :: inst, nchanl, n1,n2
51
52 ! variables for computing clwp
53 real :: clw(kms:kme), dpf(kms:kme)
54 real :: clwp
55
56 ! CRTM local varaibles and types
57 integer :: wmo_sensor_id,Error_Status, Allocate_Status
58
59 type (CRTM_RTSolution_type), allocatable :: RTSolution(:,:), RTSolution_K(:,:)
60 type (CRTM_Atmosphere_type) :: Atmosphere(1)
61 type (CRTM_Surface_type) :: Surface(1)
62 type (CRTM_Atmosphere_type), allocatable :: Atmosphere_K(:,:)
63 type (CRTM_Surface_type), allocatable :: Surface_K(:,:)
64 type (CRTM_GeometryInfo_type) :: GeometryInfo(1)
65
66 real :: t(1), a(1)
67
68 ! WHY? use argument it
69 if (it==0) then; write(unit=stdout,fmt='(A)') "WHY? have argument it to"//__FILE__; end if
70
71 alloc_status (:) = 0
72
73 if (trace_use) call da_trace_entry("da_get_innov_vector_crtmk")
74
75 ! CRTM allocation
76
77 ! Atmosphere structure
78 Atmosphere(1)%n_Layers=(kte-kts)+1 ! number of vertical levels
79 Atmosphere(1)%n_Absorbers=2
80 Atmosphere(1)%n_Clouds=0
81 Atmosphere(1)%n_Aerosols=0
82 if (crtm_cloud) Atmosphere(1)%n_Clouds=6
83
84 Error_Status = CRTM_Allocate_Atmosphere( Atmosphere(1)%n_Layers, &
85 Atmosphere(1)%n_Absorbers, &
86 Atmosphere(1)%n_Clouds, &
87 Atmosphere(1)%n_Aerosols, &
88 Atmosphere)
89 if (Error_Status /= 0) then
90 call da_error(__FILE__,__LINE__, &
91 (/"Error in allocating CRTM Atmosphere Structure"/))
92 end if
93
94 Atmosphere(1)%Absorber_ID(1)=H2O_ID
95 Atmosphere(1)%Absorber_ID(2)=O3_ID
96
97 if (crtm_cloud) then
98 Atmosphere(1)%Cloud(1)%Type=WATER_CLOUD
99 Atmosphere(1)%Cloud(2)%Type=ICE_CLOUD
100 Atmosphere(1)%Cloud(3)%Type=RAIN_CLOUD
101 Atmosphere(1)%Cloud(4)%Type=SNOW_CLOUD
102 Atmosphere(1)%Cloud(5)%Type=GRAUPEL_CLOUD
103 Atmosphere(1)%Cloud(6)%Type=HAIL_CLOUD
104 end if
105
106 !------------------------------------------------------
107 ! [1.0] calculate the background bright temperature
108 !-------------------------------------------------------
109 do inst = 1, iv%num_inst ! loop for sensor
110 ! if ( iv%instid(inst)%num_rad < 1 ) cycle
111 if (iv%instid(inst)%info%n2 < iv%instid(inst)%info%n1) cycle
112 num_levs = kte-kts+1
113 ! CRTM channel information structure
114 ! Error_Status = CRTM_Set_ChannelInfo(Sensor_Descriptor(inst),ChannelInfo)
115 ! if (Error_Status /= 0) then
116 ! call da_error(__FILE__,__LINE__, (/"Error in calling CRTM_Set_ChannelInfo"/))
117 ! end if
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 RTSolution_K(ChannelInfo(inst)%n_Channels,1), &
123 Atmosphere_K(ChannelInfo(inst)%n_Channels,1), &
124 Surface_K(ChannelInfo(inst)%n_Channels,1), &
125 STAT = Allocate_Status )
126 if ( Allocate_Status /= 0 ) then
127 call da_error(__FILE__,__LINE__, (/"Error in allocating RTSolution"/))
128 end if
129
130 ! CRTM Surface Structure
131 if (trim(crtm_sensor_name(rtminit_sensor(inst))) =='amsua') then
132 wmo_sensor_id=WMO_AMSUA
133 elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='amsub') then
134 wmo_sensor_id=WMO_AMSUB
135 elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='amsre') then
136 wmo_sensor_id=WMO_AMSRE
137 elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='ssmi') then
138 wmo_sensor_id=WMO_SSMI
139 else
140 wmo_sensor_id=INVALID_WMO_SENSOR_ID
141 end if
142
143 Error_Status = CRTM_Allocate_Surface( nchanl, & ! Input
144 Surface) ! Output
145 if (Error_Status /= 0) then
146 call da_error(__FILE__,__LINE__, &
147 (/"Error in allocating CRTM Surface Structure Structure"/))
148 end if
149 ! do n= 1, iv%instid(inst)%num_rad ! loop for pixel
150
151 n1 = iv%instid(inst)%info%n1
152 n2 = iv%instid(inst)%info%n2
153
154 allocate (model_u10(n1:n2))
155 allocate (model_v10(n1:n2))
156 allocate (model_psfc(n1:n2))
157 allocate (model_ts(n1:n2))
158 allocate (model_elv(n1:n2))
159 allocate (model_smois(n1:n2))
160 allocate (model_tslb(n1:n2))
161 allocate (model_snowh(n1:n2))
162 allocate (model_snow(n1:n2))
163 allocate (model_vegfra(n1:n2))
164
165 model_u10(:) = 0.0
166 model_v10(:) = 0.0
167 model_psfc(:) = 0.0
168 model_ts(:) = 0.0
169 model_elv(:) = 0.0
170 model_smois(:) = 0.0
171 model_tslb(:) = 0.0
172 model_snowh(:) = 0.0
173 model_snow(:) = 0.0
174 model_vegfra(:) = 0.0
175
176 do n=n1,n2
177 ! if ( n > iv%instid(inst)%num_rad ) exit
178
179 ! [1.1] Get horizontal interpolation weights:
180
181 i = iv%instid(inst)%info%i(1,n)
182 j = iv%instid(inst)%info%j(1,n)
183 dx = iv%instid(inst)%info%dx(1,n)
184 dy = iv%instid(inst)%info%dy(1,n)
185 dxm = iv%instid(inst)%info%dxm(1,n)
186 dym = iv%instid(inst)%info%dym(1,n)
187
188 ! horizontal interpolate grid%xb pressure to ob position for every grid%xb layer
189 ! get CRTM pressure Layers
190 do k=kts,kte ! from bottem to top
191 v_p(k) = dym*(dxm*grid%xb%p(i,j ,k) + dx*grid%xb%p(i+1,j ,k)) &
192 + dy *(dxm*grid%xb%p(i,j+1,k) + dx*grid%xb%p(i+1,j+1,k))
193 v_p(k) = 0.01*v_p(k) ! convert Pa to hPa
194 Atmosphere(1)%Pressure(kte-k+1)=v_p(k) ! from top to bottom
195 end do
196
197 ! determine surface type of obs location
198 !-----------------------------------------
199 call da_detsurtyp ( grid%xb%snow, grid%xb%xice, grid%xb%landmask, &
200 grid%xb%ivgtyp, grid%xb%isltyp, &
201 ims, ime, jms, jme, &
202 i, j, dx, dy, dxm, dym, &
203 model_isflg,model_ivgtyp, model_isltyp, &
204 Surface(1)%Water_Coverage, Surface(1)%Ice_Coverage, &
205 Surface(1)%Land_Coverage, Surface(1)%Snow_Coverage )
206
207 call da_interp_lin_2d_partial (grid%xb%snow, iv%instid(inst)%info, 1, n, n, model_snow(n:n) )
208
209 ! [1.2] Interpolate horizontally to ob:
210 do k=kts,kte ! from bottom to top
211 call da_interp_lin_2d_partial (grid%xb%t(:,:,k), iv%instid(inst)%info,k,n,n,t)
212 call da_interp_lin_2d_partial (grid%xb%q(:,:,k), iv%instid(inst)%info,k,n,n,a)
213 Atmosphere(1)%Temperature(kte-k+1) = t(1)
214 Atmosphere(1)%Absorber(kte-k+1,1) = a(1)
215 Atmosphere(1)%Absorber(kte-k+1,1) = 1000.0*Atmosphere(1)%Absorber(kte-k+1,1) ! in g/kg
216
217 ! NOTE: WRF high-level q values seems too big, replaced by constants
218 if (v_p(k) < 75.0) Atmosphere(1)%Absorber(kte-k+1,1) = 0.001
219
220 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))
221
222 if (crtm_cloud) then
223 Atmosphere(1)%Cloud(1)%Water_Content(kte-k+1) = model_qcw(kte-k+1)
224
225 call da_interp_lin_2d_partial (grid%xb%qci(:,:,k), iv%instid(inst)%info,k,n,n, &
226 Atmosphere(1)%Cloud(2)%Water_Content(kte-k+1:kte-k+1) )
227
228 call da_interp_lin_2d_partial (grid%xb%qrn(:,:,k), iv%instid(inst)%info,k,n,n, &
229 Atmosphere(1)%Cloud(3)%Water_Content(kte-k+1:kte-k+1) )
230
231 call da_interp_lin_2d_partial (grid%xb%qsn(:,:,k), iv%instid(inst)%info, k,n,n, &
232 Atmosphere(1)%Cloud(4)%Water_Content(kte-k+1:kte-k+1) )
233
234 call da_interp_lin_2d_partial (grid%xb%qgr(:,:,k), iv%instid(inst)%info, k,n,n, &
235 Atmosphere(1)%Cloud(5)%Water_Content(kte-k+1:kte-k+1) )
236
237 Atmosphere(1)%Cloud(6)%Water_Content(kte-k+1) = 0.0
238
239 call da_interp_lin_2d_partial (grid%xb%rho(:,:,k), iv%instid(inst)%info, k,n,n, &
240 model_rho(k:k) )
241
242 call da_cld_eff_radius(Atmosphere(1)%Temperature(kte-k+1),model_rho(k),&
243 Atmosphere(1)%Cloud(2)%Water_Content(kte-k+1), & !qci
244 Atmosphere(1)%Cloud(3)%Water_Content(kte-k+1), & !qrn
245 Atmosphere(1)%Cloud(4)%Water_Content(kte-k+1), & !qsn
246 Atmosphere(1)%Cloud(5)%Water_Content(kte-k+1), & !qgr
247 model_snow(n), &
248 Surface(1)%Ice_Coverage, Surface(1)%Land_Coverage, 1, &
249 Atmosphere(1)%Cloud(1)%Effective_Radius(kte-k+1), &
250 Atmosphere(1)%Cloud(2)%Effective_Radius(kte-k+1), &
251 Atmosphere(1)%Cloud(3)%Effective_Radius(kte-k+1), &
252 Atmosphere(1)%Cloud(4)%Effective_Radius(kte-k+1), &
253 Atmosphere(1)%Cloud(5)%Effective_Radius(kte-k+1) )
254
255 ! reset the da_cld_eff_radius calcualted effective radius to constants if desired
256
257 Atmosphere(1)%Cloud(1)%Effective_Radius(kte-k+1)=0.01 ! in mm
258 Atmosphere(1)%Cloud(2)%Effective_Radius(kte-k+1)=0.03 ! in mm
259 !Atmosphere(1)%Cloud(3)%Effective_Radius(kte-k+1)=0.3
260 !Atmosphere(1)%Cloud(4)%Effective_Radius(kte-k+1)=0.6
261 !Atmosphere(1)%Cloud(5)%Effective_Radius(kte-k+1)=0.6
262 Atmosphere(1)%Cloud(6)%Effective_Radius(kte-k+1)=0.6
263
264 end if
265 end do
266
267 call da_interp_lin_2d_partial (grid%xb%u10, iv%instid(inst)%info, 1, n, n, model_u10(n:n))
268 call da_interp_lin_2d_partial (grid%xb%v10, iv%instid(inst)%info, 1, n, n, model_v10(n:n))
269 call da_interp_lin_2d_partial (grid%xb%psfc, iv%instid(inst)%info, 1, n, n, model_psfc(n:n))
270
271 model_psfc(n) = 0.01*model_psfc(n) ! convert to hPa
272 model_ptop = 0.01*grid%xb%ptop
273
274 ! get CRTM levels (0.005hPa at top) /model full level
275 Atmosphere(1)%Level_Pressure(0)=model_ptop ! to sigma level 51-->sigmaf=0
276 Atmosphere(1)%Level_Pressure(Atmosphere(1)%n_Layers)=model_psfc(n) ! to sigma level 1->sigmaf=1
277
278 do k=kts+1,kte
279 Atmosphere(1)%Level_Pressure(kte-k+1)= grid%xb%sigmaf(k)*(model_psfc(n)-model_ptop)+model_ptop
280 end do
281
282 ! convert cloud content unit from kg/kg to kg/m^2
283 if (crtm_cloud) then
284 do k=kts,kte
285 do icld=1,Atmosphere(1)%n_Clouds
286 Atmosphere(1)%Cloud(icld)%Water_Content(k)= Atmosphere(1)%Cloud(icld)%Water_Content(k)* &
287 (Atmosphere(1)%Level_Pressure(k)- Atmosphere(1)%Level_Pressure(k-1))*100.0/gravity
288 end do
289 end do
290 end if
291
292 if ( model_isflg == 0 ) then ! over sea using SST
293 call da_interp_lin_2d_partial (grid%xb % tgrn, iv%instid(inst)%info, 1, n, n, model_ts(n:n))
294 else
295 call da_interp_lin_2d_partial (grid%xb % tsk, iv%instid(inst)%info, 1, n, n, model_ts(n:n))
296 end if
297
298 call da_interp_lin_2d_partial (grid%xb % terr, iv%instid(inst)%info, 1, n, n, model_elv(n:n))
299
300 ! variables for emissivity calculations
301 !----------------------------------------
302 call da_interp_lin_2d_partial (grid%xb % smois, iv%instid(inst)%info, 1, n, n, model_smois(n:n) )
303 call da_interp_lin_2d_partial (grid%xb % tslb, iv%instid(inst)%info, 1, n, n, model_tslb(n:n) )
304 call da_interp_lin_2d_partial (grid%xb % snowh, iv%instid(inst)%info, 1, n, n, model_snowh(n:n) )
305 call da_interp_lin_2d_partial (grid%xb % vegfra, iv%instid(inst)%info, 1, n, n, model_vegfra(n:n) )
306
307 ! model_snowh(n) = model_snowh(n)*100.0 ! convert from m to mm
308 model_vegfra(n) = 0.01*model_vegfra(n) ! convert range to 0~1
309
310 ! ADD for computing cloud liquid water path (mm) from guess
311 clwp = 0.0
312 do k = kts,kte ! from top to bottom
313 dpf(k) = 100.0*(Atmosphere(1)%level_pressure(k) - Atmosphere(1)%level_pressure(k-1))
314 clw (k) = model_qcw(k)*dpf(k)/gravity ! kg/m2 or mm
315 if (Atmosphere(1)%pressure(k)<100.0) clw(k) = 0.0
316 clwp = clwp + clw(k)
317 end do
318
319 ! CRTM GeometryInfo Structure
320 GeometryInfo(1)%Sensor_Zenith_Angle=iv%instid(inst)%satzen(n)
321 GeometryInfo(1)%Source_Zenith_Angle=iv%instid(inst)%solzen(n)
322 ! GeometryInfo(1)%Satellite_Height=830.0
323 ! GeometryInfo(1)%Sensor_Scan_Angle=
324 ! GeometryInfo(1)%Sensor_Zenith_Angle=
325 ! GeometryInfo(1)%Sensor_Scan_Angle=
326 ! GeometryInfo(1)%Source_Zenith_Angle=
327
328 ! CRTM Surface parameter data
329
330 if (Surface(1)%Land_Coverage > 0.0) then
331 Surface(1)%Land_Type=GRASS_SOIL ! land type (User guide appendix 3)
332 Surface(1)%Land_Temperature=model_ts(n) ! K
333 Surface(1)%Soil_Moisture_Content= model_smois(n) !0.05 ! volumetric water content (g/cm**3)
334 ! Surface(1)%Canopy_Water_Content=0.05 ! gravimetric water content
335 Surface(1)%Vegetation_Fraction=model_vegfra(n)
336 Surface(1)%Soil_Temperature=model_tslb(n)
337 end if
338 if (Surface(1)%Water_Coverage > 0.0) then
339 ! Surface%Water_Type=SEA_WATER ! (Currently NOT used)
340 Surface(1)%Water_Temperature=model_ts(n) ! K
341 Surface(1)%Wind_Speed=sqrt(model_u10(n)**2+model_v10(n)**2) ! m/sec
342 ! surface(1)%Wind_Direction=0.0 ! NOT used
343 Surface(1)%Salinity=33.0 ! ppmv
344 end if
345
346 if (Surface(1)%Snow_Coverage > 0.0) then
347 Surface(1)%Snow_Type=NEW_SNOW ! User guide appendix 3
348 Surface(1)%Snow_Temperature=model_ts(n) ! K
349 Surface(1)%Snow_Depth=model_snowh(n) ! mm
350 ! Surface(1)%Snow_Density=0.2 ! g/cm**3
351 ! Surface(1)%Snow_Grain_Size=2.0 ! mm
352 end if
353 if (Surface(1)%Ice_Coverage > 0.0) then
354 ! Surface(1)%Ice_Type=FRESH_ICE ! NO Table offered, single example is FRESH_ICE
355 Surface(1)%Ice_Temperature=model_ts(n) ! K
356 Surface(1)%Ice_Thickness=10.0 ! mm
357 ! Surface(1)%Ice_Density=0.9 ! g/cm**3
358 ! Surface(1)%Ice_Roughness=0.0 ! NO Table offered, single example is ZERO
359
360 end if
361 if (nchanl > 0) then
362 Surface(1)%SensorData%n_channels = nchanl
363 Surface(1)%SensorData%Sensor_ID = wmo_sensor_id
364 Surface(1)%SensorData%Tb(1:nchanl) = ob%instid(inst)%tb(1:nchanl,n)
365 end if
366
367 ! CRTM surface/atmosphere K initialization
368 do l = 1, ChannelInfo(inst)%n_Channels
369 ! -- Copy the adjoint atmosphere structure
370 Error_Status = CRTM_Assign_Atmosphere( Atmosphere, Atmosphere_K(l,:) )
371
372 if ( Error_Status /= 0 ) then
373 call da_error(__FILE__,__LINE__, &
374 (/"Error copying Atmosphere_K structure"/))
375 end if
376
377 ! -- Copy the adjoint surface structure
378 Error_Status = CRTM_Assign_Surface( Surface, Surface_K(l,:) )
379
380 if ( Error_Status /= 0 ) then
381 call da_error(__FILE__,__LINE__, &
382 (/"Error copying Surface_K structure"/))
383 end if
384 end do
385
386 ! -- Zero the Adjoint outputs
387 ! Important: adjoint variables must be initialized
388 call CRTM_Zero_Atmosphere( Atmosphere_K )
389 call CRTM_Zero_Surface( Surface_K )
390
391 ! [1.3] assign tb = R^-1 Re :
392
393 RTSolution_K(:,1)%brightness_temperature = 1.
394 RTSolution_K(:,1)%radiance = 0.
395
396 ! [1.3] Call RTM K-Matrix model
397
398 call da_crtm_k(1, nchanl, 1, Atmosphere, &
399 Surface, &
400 RTSolution_K,&
401 GeometryInfo, &
402 ChannelInfo(inst), &
403 Atmosphere_K,&
404 Surface_K, &
405 RTSolution) !, &
406 !Options = Options)
407
408 !----------------------------------------------------------------
409 ! [2.0] calculate components of innovation vector:
410 !----------------------------------------------------------------
411 do k = 1, nchanl
412 if ( iv%instid(inst)%tb_inv(k,n) > missing_r ) then
413 iv%instid(inst)%tb_xb(k,n) = RTSolution(k,1)%Brightness_Temperature
414 iv%instid(inst)%tb_inv(k,n) = &
415 ob%instid(inst)%tb(k,n) - RTSolution(k,1)%Brightness_Temperature
416 else
417 iv%instid(inst)%tb_xb(k,n) = RTSolution(k,1)%Brightness_Temperature
418 iv%instid(inst)%tb_inv(k,n) = missing_r
419 end if
420 iv%instid(inst)%emiss(k,n) = RTSolution(k,1)%surface_emissivity
421 ! surface Jacobian
422 iv%instid(inst)%ts_jacobian(k,n) = Surface_k(k,1)%water_temperature
423 iv%instid(inst)%windspeed_jacobian(k,n) = Surface_k(k,1)%wind_speed
424 iv%instid(inst)%emiss_jacobian(k,n) = RTSolution_k(k,1)%surface_emissivity
425 end do
426
427 !----------------------------------------------------------------
428 ! [3.0] store base state and Jacobian to innovation structure
429 !----------------------------------------------------------------
430 ! full level pressures
431 iv%instid(inst)%pf(0,n) = Atmosphere(1)%level_pressure(0)
432 do k=1,Atmosphere(1)%n_layers
433 iv%instid(inst)%pm(k,n) = Atmosphere(1)%pressure(k)
434 iv%instid(inst)%pf(k,n) = Atmosphere(1)%level_pressure(k)
435 iv%instid(inst)%tm(k,n) = Atmosphere(1)%temperature(k)
436 iv%instid(inst)%qm(k,n) = Atmosphere(1)%absorber(k,1)
437 ! T, Q Jacobian
438 do l=1,nchanl
439 iv%instid(inst)%t_jacobian(l,k,n) = Atmosphere_k(l,1)%temperature(k)
440 iv%instid(inst)%q_jacobian(l,k,n) = Atmosphere_k(l,1)%absorber(k,1)
441 end do
442 if (crtm_cloud) then
443 iv%instid(inst)%qcw(k,n) = Atmosphere(1)%cloud(1)%water_content(k)
444 iv%instid(inst)%qci(k,n) = Atmosphere(1)%cloud(2)%water_content(k)
445 iv%instid(inst)%qrn(k,n) = Atmosphere(1)%cloud(3)%water_content(k)
446 iv%instid(inst)%qsn(k,n) = Atmosphere(1)%cloud(4)%water_content(k)
447 iv%instid(inst)%qgr(k,n) = Atmosphere(1)%cloud(5)%water_content(k)
448 iv%instid(inst)%qhl(k,n) = Atmosphere(1)%cloud(6)%water_content(k)
449 iv%instid(inst)%rcw(k,n) = Atmosphere(1)%cloud(1)%effective_radius(k)
450 iv%instid(inst)%rci(k,n) = Atmosphere(1)%cloud(2)%effective_radius(k)
451 iv%instid(inst)%rrn(k,n) = Atmosphere(1)%cloud(3)%effective_radius(k)
452 iv%instid(inst)%rsn(k,n) = Atmosphere(1)%cloud(4)%effective_radius(k)
453 iv%instid(inst)%rgr(k,n) = Atmosphere(1)%cloud(5)%effective_radius(k)
454 iv%instid(inst)%rhl(k,n) = Atmosphere(1)%cloud(6)%effective_radius(k)
455 ! Cloud Jacobian
456 do l=1,nchanl
457 iv%instid(inst)%water_jacobian(l,k,n) = &
458 Atmosphere_k(l,1)%cloud(1)%water_content(k)
459 iv%instid(inst)%ice_jacobian(l,k,n) = &
460 Atmosphere_k(l,1)%cloud(2)%water_content(k)
461 iv%instid(inst)%rain_jacobian(l,k,n) = &
462 Atmosphere_k(l,1)%cloud(3)%water_content(k)
463 iv%instid(inst)%snow_jacobian(l,k,n) = &
464 Atmosphere_k(l,1)%cloud(4)%water_content(k)
465 iv%instid(inst)%graupel_jacobian(l,k,n) = &
466 Atmosphere_k(l,1)%cloud(5)%water_content(k)
467 iv%instid(inst)%hail_jacobian(l,k,n) = &
468 Atmosphere_k(l,1)%cloud(6)%water_content(k)
469
470 iv%instid(inst)%water_r_jacobian(l,k,n) = &
471 Atmosphere_k(l,1)%cloud(1)%effective_radius(k)
472 iv%instid(inst)%ice_r_jacobian(l,k,n) = &
473 Atmosphere_k(l,1)%cloud(2)%effective_radius(k)
474 iv%instid(inst)%rain_r_jacobian(l,k,n) = &
475 Atmosphere_k(l,1)%cloud(3)%effective_radius(k)
476 iv%instid(inst)%snow_r_jacobian(l,k,n) = &
477 Atmosphere_k(l,1)%cloud(4)%effective_radius(k)
478 iv%instid(inst)%graupel_r_jacobian(l,k,n) = &
479 Atmosphere_k(l,1)%cloud(5)%effective_radius(k)
480 iv%instid(inst)%hail_r_jacobian(l,k,n) = &
481 Atmosphere_k(l,1)%cloud(6)%effective_radius(k)
482 end do
483 end if
484 end do
485
486 !----------------------------------------------
487 ! [5.0] store surface information to innovation structure
488 !----------------------------------------------
489 iv%instid(inst)%u10(n) = model_u10(n)
490 iv%instid(inst)%v10(n) = model_v10(n)
491 iv%instid(inst)%t2m(n) = 0.01*missing_r !model_t2m
492 iv%instid(inst)%mr2m(n) = 0.01*missing_r !model_mr2m
493 iv%instid(inst)%ps(n) = model_psfc(n)
494 iv%instid(inst)%ts(n) = model_ts(n)
495 iv%instid(inst)%smois(n) = model_smois(n)
496 iv%instid(inst)%tslb(n) = model_tslb(n)
497 iv%instid(inst)%snowh(n) = model_snowh(n)
498 iv%instid(inst)%isflg(n) = model_isflg
499 iv%instid(inst)%elevation(n) = model_elv(n)
500 iv%instid(inst)%soiltyp(n) = model_isltyp
501 iv%instid(inst)%vegtyp(n) = model_ivgtyp
502 iv%instid(inst)%vegfra(n) = model_vegfra(n)
503 iv%instid(inst)%clwp(n) = clwp
504 iv%instid(inst)%water_coverage(n) = Surface(1)%water_coverage
505 iv%instid(inst)%land_coverage(n) = Surface(1)%land_coverage
506 iv%instid(inst)%ice_coverage(n) = Surface(1)%ice_coverage
507 iv%instid(inst)%snow_coverage(n) = Surface(1)%snow_coverage
508 end do ! end loop for pixels
509
510 deallocate (model_u10)
511 deallocate (model_v10)
512 deallocate (model_psfc)
513 deallocate (model_ts)
514 deallocate (model_tslb)
515 deallocate (model_snowh)
516 deallocate (model_snow)
517 deallocate (model_elv)
518 deallocate (model_vegfra)
519 deallocate (model_smois)
520
521 deallocate( RTSolution, RTSolution_K, STAT = Allocate_Status )
522
523 if (Allocate_Status /= 0) then
524 call da_error(__FILE__,__LINE__, &
525 (/"Error in deallocating RTSolution"/))
526 end if
527
528 Error_Status = CRTM_Destroy_Surface(Surface)
529 if (Error_Status /= 0) then
530 call da_error(__FILE__,__LINE__, &
531 (/"Error in deallocating CRTM Surface Structure"/))
532 end if
533
534 Error_Status = CRTM_Destroy_Surface(Surface_K)
535 if ( Error_Status /= 0 ) THEN
536 call da_error(__FILE__,__LINE__, &
537 (/"Error in deallocatting CRTM Surface_K Structure"/))
538 endif
539
540 Error_Status = CRTM_Destroy_Atmosphere( Atmosphere_K )
541 if ( Error_Status /= 0 ) THEN
542 call da_error(__FILE__,__LINE__, &
543 (/"Error in deallocatting CRTM Atmosphere_K Structure"/))
544 endif
545
546 deallocate( Atmosphere_K, Surface_K, STAT = Allocate_Status )
547 if ( Allocate_Status /= 0 ) THEN
548 call da_error(__FILE__,__LINE__, &
549 (/"Error in deallocatting CRTM Surface Structure"/))
550 endif
551
552 end do ! end loop for sensor
553
554 !------------------------------------------
555 ! 4.0 perfoming bias correction files
556 !------------------------------------------
557
558 if (biascorr) then
559 do inst = 1, iv%num_inst
560 write(unit=stdout,fmt='(A,A)') 'Performing bias correction for ', &
561 trim(iv%instid(inst)%rttovid_string)
562 call da_biascorr(inst,ob,iv)
563 end do
564 write(unit=stdout,fmt='(A)') " "
565 end if
566
567 !------------------------------------------------------------------------
568 ! [5.0] Perform QC check
569 !------------------------------------------------------------------------
570 if (qc_rad) then
571 call da_qc_crtm(ob, iv)
572 end if
573 !------------------------------------------
574 ! 6.0 prepare bias statistics files
575 !------------------------------------------
576
577 if (biasprep) then
578 do inst = 1, iv%num_inst
579 write(unit=stdout,fmt='(a,a)') 'Preparing bias statistics files for ', &
580 trim(iv%instid(inst)%rttovid_string)
581 call da_biasprep(inst,ob,iv)
582 end do
583 write(unit=stdout,fmt='(A)') " "
584 end if
585
586 Error_Status = CRTM_Destroy_Atmosphere (Atmosphere)
587 if (Error_Status /= 0) then
588 call da_error(__FILE__,__LINE__, &
589 (/"Error in deallocating CRTM Atmosphere Structure"/))
590 end if
591
592 if (trace_use) call da_trace_exit("da_get_innov_vector_crtmk")
593 #else
594 call da_error(__FILE__,__LINE__, &
595 (/"Must compile with $CRTM option for radiances"/))
596 #endif
597
598 end subroutine da_get_innov_vector_crtmk
599