| blob | 17a8d4c63581709631162111a7ec5088358e12f0 |
1 subroutine da_get_innov_vector_crtm ( it, grid, ob, iv )
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 ! HISTORY: 12/15/2008 effective radius unit is micron Zhiquan Liu
11 !---------------------------------------------------------------------------
13 implicit none
15 integer, intent(in) :: it ! External iteration.
16 type (domain), intent(in) :: grid ! first guess state.
17 type (y_type), intent(inout) :: ob ! Observation structure.
18 type (iv_type), intent(inout) :: iv ! O-B structure.
20 integer, parameter :: AIRS_Max_Channels = 281
22 integer :: n, icld ! Loop counter.
23 integer :: i, j, k ! Index dimension.
24 integer :: l ! Index dimension.
25 integer :: num_levs ! Number of obs levels.
26 real :: dx, dxm ! Interpolation weights.
27 real :: dy, dym ! Interpolation weights.
28 integer :: alloc_status(40)
30 real, allocatable :: model_u10(:)
31 real, allocatable :: model_v10(:)
32 real, allocatable :: model_psfc(:)
33 real :: model_ptop
34 real, allocatable :: model_ts(:)
35 real, allocatable :: model_elv(:)
36 real, allocatable :: model_smois(:)
37 real, allocatable :: model_tslb(:)
38 real, allocatable :: model_snowh(:)
39 real, allocatable :: model_vegfra(:)
40 real :: model_isltyp, model_ivgtyp
41 integer :: model_isflg
43 real :: model_qcw(kms:kme)
44 real :: model_qci(kms:kme)
45 real :: model_rho(kms:kme)
46 real, allocatable :: model_snow(:) ! snow water equivalent, different from model_snowh,
47 ! used in calculating reff_water
48 real, allocatable :: em_mspps(:) ! emissivity caluclated using MSPPS algorithm
49 real :: ts_mspps ! surface temperature calcualted using MSPPS algorithm
50 character(len=19) :: date_char
52 integer :: month, crtm_climat
54 ! real :: model_tm(kms:kme)
56 integer :: inst, nchanl, n1,n2
57 integer :: ipred, npred, gammapred
59 ! variables for computing clwp (clw/clwp) and cloud ice (ci/cip)
60 real :: clw(kms:kme), dpf(kms:kme), ci(kms:kme)
61 real :: clwp, cip
63 real :: total_od
65 ! CRTM local varaibles and types
66 integer :: Allocate_Status
67 integer :: n_layers, n_absorbers, n_clouds, n_aerosols
68 !! for ecmwf cloud
69 real :: t_tropp(ims:ime,jms:jme),temp1
70 integer :: k_tropp(ims:ime,jms:jme)
71 type (CRTM_RTSolution_type), allocatable :: RTSolution(:,:)
72 type (CRTM_Atmosphere_type) :: Atmosphere(1)
73 type (CRTM_Surface_type) :: Surface(1)
74 type (CRTM_Geometry_type) :: GeometryInfo(1)
75 type (CRTM_Options_type) :: Options(1)
76 type (CRTM_RTSolution_type), allocatable :: RTSolution_K(:,:)
77 type (CRTM_Atmosphere_type), allocatable :: Atmosphere_K(:,:)
78 type (CRTM_Surface_type), allocatable :: Surface_K(:,:)
80 real :: t(1), a(1), p(1)
82 integer, allocatable :: wrf_to_crtm_mw(:)
83 integer :: n_vegtype
85 !! for crtm cloud
86 real :: qcw(1),qrn(1), qci(1),qsn(1),qgr(1)
87 logical :: calc_tb_clr
88 type (CRTM_RTSolution_type), allocatable :: RTSolution_clr(:,:)
90 ! Initializations for Cloud Detection
91 integer :: ilev, jlev, nclouds
92 real, allocatable :: hessian(:,:)
93 real*8, allocatable :: eignvec(:,:), eignval(:)
94 real :: rad_clr, rad_ovc_ilev, rad_ovc_jlev
96 integer :: Band_Size(5), Bands(AIRS_Max_Channels,5)
97 !For Zhuge and Zou cloud detection
98 real, allocatable :: geoht_full(:,:,:)
99 real :: geoht_pixel(kts:min(kte,kme-1))
100 real :: tt_pixel(kts:min(kte,kme-1))
101 real :: pp_pixel(kts:min(kte,kme-1))
102 Band_Size(1:5) = (/86, 0, 0, 16, 0 /)
103 Bands(:,:) = 0
104 Bands(1:Band_Size(1),1) = &
105 & (/ & !& 1, 6, 7, 10, 11, 15, 16, 17, 20, 21, &
106 & & !& 22, 24, 27, 28, 30, 36, 39, 40, 42, 51, &
107 & & !& 52, 54, 55, 56, 59, 62, 63, 68, 69, 71, &
108 & & !& 72, 73, 74, 75, 76, 77, 78, 79, 80, 82, &
109 & 92, 93, 98, 99, 101, 104, 105, & !& 83, 84, 86, 92, 93, 98, 99, 101, 104, 105, &
110 & 108, 110, 111, 113, 116, 117, 123, 124, 128, 129, &
111 & 138, 139, 144, 145, 150, 151, 156, 157, 159, 162, &
112 & 165, 168, 169, 170, 172, 173, 174, 175, 177, 179, &
113 & 180, 182, 185, 186, 190, 192, 198, 201, 204, & !& 180, 182, 185, 186, 190, 192, 193, 198, 201, 204, &
114 & 207, 210, 215, 216, 221, 226, 227, & !& 207, 210, 213, 215, 216, 218, 221, 224, 226, 227, &
115 & 232, 252, 253, 256, 257, 261, & !& 232, 239, 248, 250, 251, 252, 253, 256, 257, 261, &
116 & 262, 267, 272, 295, 299, 305, 310, & !& 262, 267, 272, 295, 299, 300, 305, 308, 309, 310, &
117 & 321, 325, 333, 338, 355, 362, 375, 453, 475, & !& 318, 321, 325, 333, 338, 355, 362, 375, 453, 475, &
118 & 484, 497, 528, 587, 672, 787, 791, 843, 870, 914, &
119 & 950 /)
121 ! Bands(1:Band_Size(2),2) = &
122 !& (/ 1003, 1012, 1019, 1024, 1030, 1038, 1048, 1069, 1079, 1082, &
123 !& 1083, 1088, 1090, 1092, 1095, 1104, 1111, 1115, 1116, 1119, &
124 !& 1120, 1123, 1130, 1138, 1142, 1178, 1199, 1206, 1221, 1237, &
125 !& 1252, 1260, 1263, 1266, 1278, 1285 /)
127 ! Bands(1:Band_Size(3),3) = &
128 !& (/ 1301, 1304, 1329, 1371, 1382, 1415, 1424, 1449, 1455, & !& 1290, 1301, 1304, 1329, 1371, 1382, 1415, 1424, 1449, 1455, &
129 !& 1466, 1477, 1500, 1519, 1538, 1545, & !& 1466, 1471, 1477, 1479, 1488, 1500, 1519, 1520, 1538, 1545, &
130 !& 1565, 1574, 1583, 1593, 1627, 1636, 1652, 1669, & !& 1565, 1574, 1583, 1593, 1614, 1627, 1636, 1644, 1652, 1669, &
131 !& 1694, 1708, 1723, 1740, 1748, 1756, & !& 1674, 1681, 1694, 1708, 1717, 1723, 1740, 1748, 1751, 1756, &
132 !& 1766, 1771, 1777, 1783, 1794, 1800, 1806, & !& 1763, 1766, 1771, 1777, 1780, 1783, 1794, 1800, 1803, 1806, &
133 !& 1826, 1843 /) !& 1812, 1826, 1843 /)
135 Bands(1:Band_Size(4),4) = &
136 & (/ 1852, 1865, 1866, 1868, 1869, 1872, 1873, 1876, & !& 1852, 1865, 1866, 1867, 1868, 1869, 1872, 1873, 1875, 1876,
137 & 1881, 1882, 1883, 1911, 1917, 1918, & !& 1877, 1881, 1882, 1883, 1884, 1897, 1901, 1911, 1917, 1918, &
138 & 1924, 1928 /) !& 1921, 1923, 1924, 1928, 1937 /)
140 ! Bands(1:Band_Size(5),5) = &
141 !& (/ 1938, 1939, 1941, 1946, 1947, 1948, 1958, 1971, 1973, 1988, &
142 !& 1995, 2084, 2085, 2097, 2098, 2099, 2100, 2101, 2103, 2104, &
143 !& 2106, 2107, 2108, 2109, 2110, 2111, 2112, 2113, 2114, 2115, &
144 !& 2116, 2117, 2118, 2119, 2120, 2121, 2122, 2123, 2128, 2134, &
145 !& 2141, 2145, 2149, 2153, 2164, 2189, 2197, 2209, 2226, 2234, &
146 !& 2280, 2318, 2321, 2325, 2328, 2333, 2339, 2348, 2353, 2355, &
147 !& 2363, 2370, 2371, 2377 /)
149 ! WHY? use argument it
150 if (it==0) then; write(unit=stdout,fmt='(A)') "WHY? have argument it to"//__FILE__; end if
152 alloc_status (:) = 0
154 if (trace_use) call da_trace_entry("da_get_innov_vector_crtm")
156 ! CRTM allocation
158 ! CRTM Atmosphere
159 n_layers = (kte-kts)+1 ! number of vertical levels
160 n_absorbers = 2
161 n_aerosols = 0
162 n_clouds = 0
163 if ( crtm_cloud ) n_clouds = 6
165 call CRTM_Atmosphere_Create ( Atmosphere(1), &
166 n_layers, &
167 n_absorbers, &
168 n_clouds, &
169 n_aerosols )
170 IF ( .NOT. CRTM_Atmosphere_Associated( Atmosphere(1) ) ) THEN
171 call da_error(__FILE__,__LINE__, &
172 (/"Error in allocating CRTM Atmosphere Structure"/))
173 END IF
175 Atmosphere(1)%Absorber_ID(1)=H2O_ID
176 Atmosphere(1)%Absorber_ID(2)=O3_ID
177 Atmosphere(1)%Absorber_Units(1) = MASS_MIXING_RATIO_UNITS
178 Atmosphere(1)%Absorber_Units(2) = VOLUME_MIXING_RATIO_UNITS
180 if (crtm_cloud) then
181 Atmosphere(1)%Cloud(1)%Type=WATER_CLOUD
182 Atmosphere(1)%Cloud(2)%Type=ICE_CLOUD
183 Atmosphere(1)%Cloud(3)%Type=RAIN_CLOUD
184 Atmosphere(1)%Cloud(4)%Type=SNOW_CLOUD
185 Atmosphere(1)%Cloud(5)%Type=GRAUPEL_CLOUD
186 Atmosphere(1)%Cloud(6)%Type=HAIL_CLOUD
187 end if
189 select case (trim(grid%mminlu))
190 case ('USGS')
191 n_vegtype = USGS_n_type
192 case ('MODIFIED_IGBP_MODIS_NOAH')
193 n_vegtype = IGBP_n_type
194 case default
195 call da_error(__FILE__,__LINE__,(/"Unknown land use type"/))
196 end select
197 allocate(wrf_to_crtm_mw(n_vegtype))
198 if ( n_vegtype == USGS_n_type ) then
199 wrf_to_crtm_mw = usgs_to_crtm_mw
200 else if ( n_vegtype == IGBP_n_type ) then
201 wrf_to_crtm_mw = igbp_to_crtm_mw
202 end if
203 if ( use_clddet_zz ) then
204 allocate ( geoht_full(ims:ime,jms:jme,kms:kme-1) )
205 do k = kms, kme-1
206 do j = jms, jme
207 do i = ims, ime
208 geoht_full(i,j,k) = 0.5 * ( grid%ph_2(i,j,k) + grid%phb(i,j,k) + &
209 grid%ph_2(i,j,k+1) + grid%phb(i,j,k+1) ) / gravity
210 end do
211 end do
212 end do
213 end if
214 !------------------------------------------------------
215 ! [1.0] calculate the background bright temperature
216 !-------------------------------------------------------
217 do inst = 1, iv%num_inst ! loop for sensor
218 ! if ( iv%instid(inst)%num_rad < 1 ) cycle
219 if (iv%instid(inst)%info%n2 < iv%instid(inst)%info%n1) cycle
220 num_levs = kte-kts+1
221 nchanl = ChannelInfo(inst)%n_channels
223 if ( ChannelInfo(inst)%Sensor_Type == 2 ) then ! INFRARED_SENSOR=2
224 if ( index(grid%mminlu, trim(CRTM_IRlandCoeff_Classification())) <= 0 ) then
225 call da_error(__FILE__,__LINE__, &
226 (/"modify crtm_irland_coef to use the coef file that is consistent with your model land use type"/))
227 end if
228 end if
230 ! Allocate forward model solution RTSolution array to number of channels
231 allocate (RTSolution(ChannelInfo(inst)%n_Channels,1), STAT = Allocate_Status )
232 if ( Allocate_Status /= 0 ) then
233 call da_error(__FILE__,__LINE__, (/"Error in allocating RTSolution"/))
234 end if
236 call CRTM_RTSolution_Create( RTSolution , & !Output
237 Atmosphere(1)%n_layers ) !Input
238 if ( .NOT.(ANY(crtm_rtsolution_associated(RTSolution))) ) then
239 call da_error(__FILE__,__LINE__, &
240 (/"Error in allocating CRTM RTSolution Stucture"/))
241 end if
243 calc_tb_clr = .false.
244 if ( crtm_cloud .and. &
245 ( trim( crtm_sensor_name(rtminit_sensor(inst))) == 'amsr2' .or. &
246 trim( crtm_sensor_name(rtminit_sensor(inst))) == 'abi' .or. &
247 trim( crtm_sensor_name(rtminit_sensor(inst))) == 'ahi') ) then
248 !Tb_clear_sky is only needed for symmetric obs error model
249 !symmetric obs error model only implemented for amsr2 & abi/ahi for now
250 calc_tb_clr = .true.
251 end if
253 if ( calc_tb_clr ) then
254 allocate (RTSolution_clr(ChannelInfo(inst)%n_Channels,1), STAT = Allocate_Status )
255 if ( Allocate_Status /= 0 ) then
256 call da_error(__FILE__,__LINE__, (/"Error in allocating RTSolution_clr"/))
257 end if
258 call CRTM_RTSolution_Create( RTSolution_clr , & !Output
259 Atmosphere(1)%n_layers ) !Input
260 if ( .NOT.(ANY(crtm_rtsolution_associated(RTSolution_clr))) ) then
261 call da_error(__FILE__,__LINE__, &
262 (/"Error in allocating CRTM RTSolution_clr Stucture"/))
263 end if
264 end if !calc_tb_clr
266 if (use_crtm_kmatrix) then
267 allocate (RTSolution_K(ChannelInfo(inst)%n_Channels,1), &
268 Atmosphere_K(ChannelInfo(inst)%n_Channels,1), &
269 Surface_K(ChannelInfo(inst)%n_Channels,1), &
270 STAT = Allocate_Status )
271 if ( Allocate_Status /= 0 ) then
272 call da_error(__FILE__,__LINE__, (/"Error in allocating RTSolution_K"/))
273 end if
274 call CRTM_RTSolution_Create( RTSolution_K , & !Output
275 Atmosphere(1)%n_layers ) !Input
276 if ( .NOT.(ANY(crtm_rtsolution_associated(RTSolution_K))) ) then
277 call da_error(__FILE__,__LINE__, &
278 (/"Error in allocating CRTM RTSolution_K Stucture"/))
279 end if
280 end if
282 call CRTM_Surface_Create ( Surface(1), & ! Output
283 nchanl ) ! Input
284 IF ( .NOT. CRTM_Surface_Associated( Surface(1) ) ) THEN
285 call da_error(__FILE__,__LINE__, &
286 (/"Error in allocating CRTM Surface Structure"/))
287 end if
289 ! CRTM Options structure
290 call CRTM_Options_Create( Options(1), & ! Output
291 nchanl ) ! Input
292 IF ( .NOT. CRTM_Options_Associated( Options(1) ) ) THEN
293 call da_error(__FILE__,__LINE__, &
294 (/"Error in allocatting CRTM Options Structure"/))
295 endif
296 if ( use_antcorr(inst) ) Options(1)%Use_Antenna_Correction = .true.
298 ! do n= 1, iv%instid(inst)%num_rad ! loop for pixel
300 n1 = iv%instid(inst)%info%n1
301 n2 = iv%instid(inst)%info%n2
303 allocate (model_u10(n1:n2))
304 allocate (model_v10(n1:n2))
305 allocate (model_psfc(n1:n2))
306 allocate (model_ts(n1:n2))
307 allocate (model_elv(n1:n2))
308 allocate (model_smois(n1:n2))
309 allocate (model_tslb(n1:n2))
310 allocate (model_snowh(n1:n2))
311 allocate (model_snow(n1:n2))
312 allocate (model_vegfra(n1:n2))
314 model_u10(:) = 0.0
315 model_v10(:) = 0.0
316 model_psfc(:) = 0.0
317 model_ts(:) = 0.0
318 model_elv(:) = 0.0
319 model_smois(:) = 0.0
320 model_tslb(:) = 0.0
321 model_snowh(:) = 0.0
322 model_snow(:) = 0.0
323 model_vegfra(:) = 0.0
325 ! Gamma correction from VarBC
326 !----------------------------
327 !#ifdef CRTM_MODIF
328 if ( use_varbc .or. freeze_varbc ) then
329 gammapred = iv%instid(inst)%varbc_info%gammapred
330 do k = 1, nchanl
331 npred = iv%instid(inst)%varbc(k)%npred
332 if (npred <= 0) cycle ! VarBC channels only
333 if (iv%instid(inst)%varbc_info%npredmax < gammapred) cycle ! Gamma channels only
334 if (iv%instid(inst)%varbc(k)%pred_use(gammapred) < 0) cycle ! Gamma channels only
335 do ipred = 1, npred
336 if (iv%instid(inst)%varbc(k)%ipred(ipred) /= gammapred) cycle
337 RTSolution(k,1)%Gamma = iv%instid(inst)%varbc(k)%param(ipred)
338 end do
339 end do
340 end if
341 RTSolution_K(:,1)%Gamma = 0.0
342 !#endif
344 if ( use_mspps_emis(inst) ) allocate(em_mspps(nchanl))
346 pixel_loop: do n=n1,n2
347 ! if ( n > iv%instid(inst)%num_rad ) exit
349 ! [1.1] Get horizontal interpolation weights:
351 i = iv%instid(inst)%info%i(1,n)
352 j = iv%instid(inst)%info%j(1,n)
353 dx = iv%instid(inst)%info%dx(1,n)
354 dy = iv%instid(inst)%info%dy(1,n)
355 dxm = iv%instid(inst)%info%dxm(1,n)
356 dym = iv%instid(inst)%info%dym(1,n)
358 ! determine CRTM climatology based on latitude and month
359 date_char=iv%instid(inst)%info%date_char(n)
360 read(unit=date_char(6:7), fmt='(i2.2)') month
361 call da_det_crtm_climat(iv%instid(inst)%info%lat(1,n), month, crtm_climat)
362 iv%instid(inst)%crtm_climat(n) = crtm_climat
363 Atmosphere(1)%Climatology = iv%instid(inst)%crtm_climat(n)
365 ! determine surface type of obs location
366 !-----------------------------------------
367 call da_detsurtyp ( grid%xb%snow, grid%xb%xice, grid%xb%landmask, &
368 grid%xb%ivgtyp, grid%xb%isltyp, &
369 ims, ime, jms, jme, &
370 i, j, dx, dy, dxm, dym, &
371 model_isflg,model_ivgtyp, model_isltyp, &
372 Surface(1)%Water_Coverage, Surface(1)%Ice_Coverage, &
373 Surface(1)%Land_Coverage, Surface(1)%Snow_Coverage )
375 call da_interp_2d_partial (grid%xb%snow, iv%instid(inst)%info, 1, n, n, model_snow(n:n))
377 ! [1.2] Interpolate horizontally to ob:
378 do k=kts,kte ! from bottom to top
379 call da_interp_2d_partial (grid%xb%p(:,:,k), iv%instid(inst)%info,k,n,n,p)
380 call da_interp_2d_partial (grid%xb%t(:,:,k), iv%instid(inst)%info,k,n,n,t)
381 call da_interp_2d_partial (grid%xb%q(:,:,k), iv%instid(inst)%info,k,n,n,a)
383 Atmosphere(1)%Pressure(kte-k+1) = 0.01 * p(1) ! convert Pa to hPa
384 Atmosphere(1)%Temperature(kte-k+1) = t(1)
385 Atmosphere(1)%Absorber(kte-k+1,1) = 1000.0 * a(1)/(1.0-a(1)) ! in g/kg mixing ratio
387 ! NOTE: WRF high-level q values seems too big, replaced by constants
388 if (p(1)*0.01 < 75.0) Atmosphere(1)%Absorber(kte-k+1,1) = 0.001
390 call da_interp_2d_partial (grid%xb%qcw(:,:,k), iv%instid(inst)%info,k,n,n, model_qcw(kte-k+1:kte-k+1))
392 if (crtm_cloud) then
394 call da_interp_2d_partial (grid%xb%qci(:,:,k), iv%instid(inst)%info,k,n,n,qci)
395 call da_interp_2d_partial (grid%xb%qci(:,:,k), iv%instid(inst)%info,k,n,n,model_qci(kte-k+1:kte-k+1))
397 call da_interp_2d_partial (grid%xb%qrn(:,:,k), iv%instid(inst)%info,k,n,n,qrn)
399 call da_interp_2d_partial (grid%xb%qsn(:,:,k), iv%instid(inst)%info, k,n,n,qsn)
401 call da_interp_2d_partial (grid%xb%qgr(:,:,k), iv%instid(inst)%info, k,n,n,qgr)
403 Atmosphere(1)%Cloud(1)%Water_Content(kte-k+1)=model_qcw(kte-k+1)
404 Atmosphere(1)%Cloud(2)%Water_Content(kte-k+1)=qci(1)
405 Atmosphere(1)%Cloud(3)%Water_Content(kte-k+1)=qrn(1)
406 Atmosphere(1)%Cloud(4)%Water_Content(kte-k+1)=qsn(1)
407 Atmosphere(1)%Cloud(5)%Water_Content(kte-k+1)=qgr(1)
408 Atmosphere(1)%Cloud(6)%Water_Content(kte-k+1)=0.0
410 call da_interp_2d_partial (grid%xb%rho(:,:,k), iv%instid(inst)%info, k,n,n, &
411 model_rho(k:k) )
413 call da_cld_eff_radius(Atmosphere(1)%Temperature(kte-k+1),model_rho(k),&
414 Atmosphere(1)%Cloud(2)%Water_Content(kte-k+1), & !qci
415 Atmosphere(1)%Cloud(3)%Water_Content(kte-k+1), & !qrn
416 Atmosphere(1)%Cloud(4)%Water_Content(kte-k+1), & !qsn
417 Atmosphere(1)%Cloud(5)%Water_Content(kte-k+1), & !qgr
418 model_snow(n), &
419 Surface(1)%Ice_Coverage, Surface(1)%Land_Coverage, 1, &
420 Atmosphere(1)%Cloud(1)%Effective_Radius(kte-k+1), &
421 Atmosphere(1)%Cloud(2)%Effective_Radius(kte-k+1), &
422 Atmosphere(1)%Cloud(3)%Effective_Radius(kte-k+1), &
423 Atmosphere(1)%Cloud(4)%Effective_Radius(kte-k+1), &
424 Atmosphere(1)%Cloud(5)%Effective_Radius(kte-k+1) )
426 ! reset the da_cld_eff_radius calcualted effective radius to constants if desired
428 Atmosphere(1)%Cloud(1)%Effective_Radius(kte-k+1)=10.0 ! in micron
429 Atmosphere(1)%Cloud(2)%Effective_Radius(kte-k+1)=30.0 ! in micron
430 !Atmosphere(1)%Cloud(3)%Effective_Radius(kte-k+1)=300
431 !Atmosphere(1)%Cloud(4)%Effective_Radius(kte-k+1)=600
432 !Atmosphere(1)%Cloud(5)%Effective_Radius(kte-k+1)=600
433 Atmosphere(1)%Cloud(6)%Effective_Radius(kte-k+1)=600
435 end if
436 end do
437 if (use_clddet_zz) then
438 ! Find tropopause temperature for Zhuge and Zou Cloud Detection
439 do k = kts, min(kte,kme-1)
440 tt_pixel(k) = Atmosphere(1)%Temperature(kte-k+1)
441 pp_pixel(k) = Atmosphere(1)%Pressure(kte-k+1)
442 call da_interp_2d_partial ( geoht_full(:,:,k), iv%instid(inst)%info, k, n, n, geoht_pixel(k) )
443 end do
444 call da_trop_wmo ( tt_pixel, geoht_pixel, pp_pixel, (min(kte,kme-1)-kts+1), tropt = iv%instid(inst)%tropt(n) )
445 end if
447 call da_interp_2d_partial (grid%xb%u10, iv%instid(inst)%info, 1, n, n, model_u10(n:n))
448 call da_interp_2d_partial (grid%xb%v10, iv%instid(inst)%info, 1, n, n, model_v10(n:n))
449 call da_interp_2d_partial (grid%xb%psfc, iv%instid(inst)%info, 1, n, n, model_psfc(n:n))
451 model_psfc(n) = 0.01*model_psfc(n) ! convert to hPa
452 model_ptop = 0.01*grid%xb%ptop
454 ! get CRTM levels (0.005hPa at top) /model full level
455 Atmosphere(1)%Level_Pressure(0) = model_ptop ! to sigma level 51-->sigmaf=0
456 Atmosphere(1)%Level_Pressure(Atmosphere(1)%n_Layers) = model_psfc(n) ! to sigma level 1->sigmaf=1
458 ! calculating full-level P from half-level P using
459 ! linear interpolation: (pf-pf1)/(sigmaf-sigmaf1)=(ph2-ph1)/(sigmah2-sigmah1)
460 ! note: sigma is from bottem to top; pressure is from top to bottem
461 do k=kts,kte-1 ! from top to bottem
462 Atmosphere(1)%Level_Pressure(k)= (Atmosphere(1)%Pressure(k+1) - Atmosphere(1)%Pressure(k))/ &
463 (grid%xb%sigmah(kte-k) - grid%xb%sigmah(kte-k+1)) * &
464 (grid%xb%sigmaf(kte-k+1) - grid%xb%sigmah(kte-k+1)) + &
465 Atmosphere(1)%Pressure(k)
466 end do
468 ! check for pressure monotonicity
469 do k = kte, kts+1, -1 ! from bottom to top
470 if ( Atmosphere(1)%Level_Pressure(k) <= Atmosphere(1)%Level_Pressure(k-1) ) then
471 write(message(1),'(a,2i5.0,a)') ' Skipping the pixel at loc ', i, j, &
472 ' where the pressure is not monotonic'
473 call da_warning(__FILE__,__LINE__,message(1:1))
474 iv%instid(inst)%tb_inv(:,n) = missing_r
475 iv%instid(inst)%info%proc_domain(:,n) = .false.
476 cycle pixel_loop
477 end if
478 end do
479 !if ( all(ob%instid(inst)%tb(1:nchanl,n) < 0.) ) then
480 ! write(message(1),'(a,2i5.0,a)') ' Skipping the pixel at loc ', i, j, &
481 ! ' where all observed BTs are < 0'
482 ! call da_warning(__FILE__,__LINE__,message(1:1))
483 ! iv%instid(inst)%tb_inv(:,n) = missing_r
484 ! iv%instid(inst)%info%proc_domain(:,n) = .false.
485 ! cycle pixel_loop
486 !end if
488 ! convert cloud content unit from kg/kg to kg/m^2
489 if (crtm_cloud) then
490 do k=kts,kte
491 do icld=1,Atmosphere(1)%n_Clouds
492 Atmosphere(1)%Cloud(icld)%Water_Content(k)= Atmosphere(1)%Cloud(icld)%Water_Content(k)* &
493 (Atmosphere(1)%Level_Pressure(k)- Atmosphere(1)%Level_Pressure(k-1))*100.0/gravity
494 end do
495 end do
496 end if
498 if ( model_isflg == 0 ) then ! over sea using SST
499 call da_interp_2d_partial (grid%xb % tgrn, iv%instid(inst)%info, 1, n, n, model_ts(n:n))
500 else
501 call da_interp_2d_partial (grid%xb % tsk, iv%instid(inst)%info, 1, n, n, model_ts(n:n))
502 end if
504 call da_interp_2d_partial (grid%xb % terr, iv%instid(inst)%info, 1, n, n, model_elv(n:n))
506 ! variables for emissivity calculations
507 !----------------------------------------
508 call da_interp_2d_partial (grid%xb % smois, iv%instid(inst)%info, 1, n, n, model_smois(n:n) )
509 call da_interp_2d_partial (grid%xb % tslb, iv%instid(inst)%info, 1, n, n, model_tslb(n:n) )
510 call da_interp_2d_partial (grid%xb % snowh, iv%instid(inst)%info, 1, n, n, model_snowh(n:n) )
511 call da_interp_2d_partial (grid%xb % vegfra, iv%instid(inst)%info, 1, n, n, model_vegfra(n:n) )
513 ! model_snowh(n) = model_snowh(n)*100.0 ! convert from m to mm
514 model_vegfra(n) = 0.01*model_vegfra(n) ! convert range to 0~1
516 ! ADD for computing cloud liquid water path (mm) from guess
517 clwp = 0.0
518 do k = kts,kte ! from top to bottom
519 dpf(k) = 100.0*(Atmosphere(1)%level_pressure(k) - Atmosphere(1)%level_pressure(k-1))
520 clw (k) = model_qcw(k)*dpf(k)/gravity ! kg/m2 or mm
521 if (Atmosphere(1)%pressure(k)<100.0) clw(k) = 0.0
522 clwp = clwp + clw(k)
523 end do
525 ! computing ice path from guess
526 if (crtm_cloud) then
527 cip = 0.0
528 do k = kts,kte ! from top to bottom
529 dpf(k) = 100.0*(Atmosphere(1)%level_pressure(k) - Atmosphere(1)%level_pressure(k-1))
530 ci(k) = model_qci(k)*dpf(k)/gravity ! kg/m2 or mm
531 cip = cip + ci(k)
532 end do
533 endif
535 ! CRTM GeometryInfo Structure
536 GeometryInfo(1)%Sensor_Zenith_Angle=iv%instid(inst)%satzen(n)
537 GeometryInfo(1)%Source_Zenith_Angle=iv%instid(inst)%solzen(n)
538 GeometryInfo(1)%iFOV=iv%instid(inst)%scanpos(n)
539 ! GeometryInfo(1)%Satellite_Height=830.0
540 ! GeometryInfo(1)%Sensor_Scan_Angle=
541 ! GeometryInfo(1)%Sensor_Zenith_Angle=
542 ! GeometryInfo(1)%Sensor_Scan_Angle=
543 ! GeometryInfo(1)%Source_Zenith_Angle=
545 ! CRTM Surface parameter data
547 if ( use_mspps_ts(inst) ) then
548 ! only for AMSU-A over land
549 if ( trim(crtm_sensor_name(rtminit_sensor(inst))) == 'amsua' &
550 .and. model_isflg == 2 ) then
551 call da_mspps_ts(ob%instid(inst)%tb(1:nchanl,n), nchanl, &
552 iv%instid(inst)%satzen(n), ts_mspps)
553 ! ts_mspps is initilaized as negative values in the
554 ! da_mspps_ts subroutine. Apply only valid values here.
555 if ( ts_mspps > 0.0 ) then
556 model_ts(n) = ts_mspps
557 end if
558 end if
559 end if
561 if (Surface(1)%Land_Coverage > 0.0) then
562 Surface(1)%Land_Type = min(11,nint(model_ivgtyp))
563 Surface(1)%Vegetation_Type = max(1,wrf_to_crtm_mw(nint(model_ivgtyp)))
564 Surface(1)%Soil_Type = max(1,wrf_to_crtm_soil(nint(model_isltyp)))
565 ! glacial land ice soil type and glacial vegetation type are not valid land types
566 if ( surface(1)%Soil_Type == 9 .or. surface(1)%Vegetation_Type == 13 ) then
567 surface(1)%ice_coverage = min(surface(1)%ice_coverage + surface(1)%land_coverage, 1.0)
568 surface(1)%land_coverage = 0
569 end if
570 ! reset water-dominated land type
571 if ( Surface(1)%Land_Type == grid%iswater ) Surface(1)%Land_Type = 1
572 Surface(1)%Land_Temperature=model_ts(n) ! K
573 Surface(1)%Soil_Moisture_Content= model_smois(n) !0.05 ! volumetric water content (g/cm**3)
574 ! Surface(1)%Canopy_Water_Content=0.05 ! gravimetric water content
575 Surface(1)%Vegetation_Fraction=model_vegfra(n)
576 Surface(1)%Soil_Temperature=model_tslb(n)
577 end if
578 if (Surface(1)%Water_Coverage > 0.0) then
579 ! Surface%Water_Type=SEA_WATER ! (Currently NOT used)
580 Surface(1)%Water_Temperature=model_ts(n) ! K
581 Surface(1)%Wind_Speed=sqrt(model_u10(n)**2+model_v10(n)**2) ! m/sec
582 ! surface(1)%Wind_Direction=0.0 ! NOT used
583 Surface(1)%Salinity=33.0 ! ppmv
584 end if
586 if (Surface(1)%Snow_Coverage > 0.0) then
587 Surface(1)%Snow_Type=2 !NEW_SNOW
588 Surface(1)%Snow_Temperature=model_ts(n) ! K
589 Surface(1)%Snow_Depth=model_snowh(n) ! mm
590 ! Surface(1)%Snow_Density=0.2 ! g/cm**3
591 ! Surface(1)%Snow_Grain_Size=2.0 ! mm
592 end if
593 if (Surface(1)%Ice_Coverage > 0.0) then
594 ! Surface(1)%Ice_Type=FRESH_ICE ! NO Table offered, single example is FRESH_ICE
595 Surface(1)%Ice_Temperature=model_ts(n) ! K
596 Surface(1)%Ice_Thickness=10.0 ! mm
597 ! Surface(1)%Ice_Density=0.9 ! g/cm**3
598 ! Surface(1)%Ice_Roughness=0.0 ! NO Table offered, single example is ZERO
600 end if
601 if (nchanl > 0) then
602 Surface(1)%SensorData%n_channels = nchanl
603 Surface(1)%SensorData%Sensor_Channel = ChannelInfo(inst)%Sensor_Channel
604 Surface(1)%SensorData%Sensor_Id = ChannelInfo(inst)%Sensor_Id
605 Surface(1)%SensorData%WMO_Satellite_Id = ChannelInfo(inst)%WMO_Satellite_Id
606 Surface(1)%SensorData%WMO_Sensor_Id = ChannelInfo(inst)%WMO_Sensor_Id
607 Surface(1)%SensorData%Tb(1:nchanl) = ob%instid(inst)%tb(1:nchanl,n)
608 end if
610 ! user emissivity options
611 Options(1)%use_emissivity = .false.
612 if ( use_mspps_emis(inst) ) then
613 ! Only for AMSU-A over land
614 if ( trim(crtm_sensor_name(rtminit_sensor(inst))) == 'amsua' &
615 .and. model_isflg == 2 ) then
616 call da_mspps_emis(ob%instid(inst)%tb(1:nchanl,n), nchanl, em_mspps)
617 Options(1)%use_emissivity = .true.
618 Options(1)%emissivity(1:nchanl) = em_mspps(1:nchanl)
619 end if
620 end if
622 if (use_crtm_kmatrix) then
624 ! CRTM surface/atmosphere K initialization
625 do l = 1, ChannelInfo(inst)%n_Channels
626 ! -- Copy the adjoint atmosphere structure
627 Atmosphere_K(l,1) = Atmosphere(1)
629 ! -- Copy the adjoint surface structure
630 Surface_K(l,1) = Surface(1)
631 end do
633 ! -- Zero the Adjoint outputs
634 ! Important: adjoint variables must be initialized
635 call CRTM_Atmosphere_Zero( Atmosphere_K )
636 call CRTM_Surface_Zero( Surface_K )
638 ! Assign tb = R^-1 Re :
639 RTSolution_K(:,1)%brightness_temperature = 1.
640 RTSolution_K(:,1)%radiance = 0.
642 ! [1.3] Call RTM K-Matrix model
643 call da_crtm_k(1, nchanl, 1, Atmosphere, &
644 Surface, &
645 RTSolution_K,&
646 GeometryInfo, &
647 ChannelInfo(inst), &
648 Atmosphere_K,&
649 Surface_K, &
650 RTSolution, &
651 Options)
652 else
654 ! [1.3] Call RTM forward model
655 call da_crtm_direct(1, nchanl, 1, Atmosphere, &
656 Surface, &
657 GeometryInfo, &
658 ChannelInfo(inst:inst), &
659 RTSolution, &
660 Options)
662 end if
664 ! Compute Overcast Radiances for MMR in Auligné (2014).or. PF in Xu et al., (2016).or. MW03 ecmwf method in McNally AP and Watts PD (2003)
665 !----------------------------------------------------------------
666 if (use_clddet==1 .or. use_clddet==2.or. use_clddet==3 ) then
667 do i = 1, nchanl
668 !crtm_2.2.3 contains Upwelling_Overcast_Radiance
669 iv%instid(inst)%rad_ovc(i,kts:kte,n)=RTSolution(i,1)%Upwelling_Overcast_Radiance(:)
670 end do
671 end if
673 !---------------------------------------------------------------------------
674 ! Precondition satellite control variable (Cloud Cover(s)):
675 !---------------------------------------------------------------------------
676 if (.false.) then
677 ! if (use_satcv(2)) then
678 nclouds = iv%instid(inst)%cv_index(n)%nclouds
679 allocate (hessian(nclouds,nclouds))
680 allocate (eignvec(nclouds,nclouds))
681 allocate (eignval(nclouds))
682 hessian(:,:)= 0.0
683 do ilev=1, nclouds
684 do jlev=ilev, nclouds
685 do k = 1, nchanl
686 if (ALL(iv%instid(inst)%ichan(k) /= Bands(:,1))) cycle ! Only Channels in Band 1
687 rad_clr = iv%instid(inst)%rad_xb(k,n)
688 rad_ovc_ilev = iv%instid(inst)%rad_ovc(k,kte-nclouds+ilev,n)
689 rad_ovc_jlev = iv%instid(inst)%rad_ovc(k,kte-nclouds+jlev,n)
690 hessian(ilev,jlev) = hessian(ilev,jlev) + (rad_ovc_ilev-rad_clr) * (rad_ovc_jlev-rad_clr)
691 end do
692 hessian(jlev,ilev) = hessian(ilev,jlev)
693 end do
694 end do
696 call da_eof_decomposition(nclouds, hessian, eignvec, eignval)
698 ! if (ANY(eignval <= 0)) write(unit=stdout,fmt='(3A,I4,A,100F18.5)') &
699 ! 'SATCV: non-positive Hessian for ', trim(iv%instid(inst)%rttovid_string), ' ,pixel ',n,'--> Eigenvalues =',eignval
701 do ilev = 1, nclouds
702 do jlev = ilev, nclouds
703 iv%instid(inst)%cv_index(n)%vtox(ilev,jlev) = &
704 sum( eignvec(ilev,:) * sqrt(1.0/eignval(:)) * eignvec(jlev,:), mask = eignval >0 )
706 iv%instid(inst)%cv_index(n)%vtox(jlev,ilev) = iv%instid(inst)%cv_index(n)%vtox(ilev,jlev)
707 end do
708 end do
710 deallocate(hessian,eignvec,eignval)
711 end if
713 !----------------------------------------------------------------
714 ! [2.0] calculate components of innovation vector:
715 !----------------------------------------------------------------
716 do k = 1, nchanl
717 iv%instid(inst)%tb_xb(k,n) = RTSolution(k,1)%Brightness_Temperature
718 iv%instid(inst)%rad_xb(k,n) = RTSolution(k,1)%Radiance
719 iv%instid(inst)%emiss(k,n) = RTSolution(k,1)%surface_emissivity
721 if (use_pseudo_rad .or. use_simulated_rad) then ! input is innovation
722 ob%instid(inst)%tb(k,n) = missing_r
723 if ( iv%instid(inst)%tb_inv(k,n) > missing_r ) &
724 ob%instid(inst)%tb(k,n) = RTSolution(k,1)%Brightness_Temperature + iv%instid(inst)%tb_inv(k,n)
725 else
726 if ( iv%instid(inst)%tb_inv(k,n) > missing_r ) &
727 iv%instid(inst)%tb_inv(k,n) = ob%instid(inst)%tb(k,n) - iv%instid(inst)%tb_xb(k,n)
728 end if
730 if (use_crtm_kmatrix) then
731 ! surface Jacobian
732 iv%instid(inst)%ts_jacobian(k,n) = Surface_k(k,1)%water_temperature
733 iv%instid(inst)%windspeed_jacobian(k,n) = Surface_k(k,1)%wind_speed
734 iv%instid(inst)%emiss_jacobian(k,n) = RTSolution_k(k,1)%surface_emissivity
735 !#ifdef CRTM_MODIF
736 iv%instid(inst)%gamma_jacobian(k,n) = RTSolution_k(k,1)%Gamma
737 !#endif
738 end if
739 end do
741 !----------------------------------------------------------------
742 ! [3.0] store base state (and Jacobian) to innovation structure
743 !----------------------------------------------------------------
744 ! full level pressures
745 iv%instid(inst)%pf(0,n) = Atmosphere(1)%level_pressure(0)
746 if (use_crtm_kmatrix) then
747 ! PS Jacobian
748 do l=1,nchanl
749 iv%instid(inst)%ps_jacobian(l,n) = Atmosphere_k(l,1)%level_pressure(Atmosphere(1)%n_layers)
750 end do
751 end if
752 do k=1,Atmosphere(1)%n_layers
753 iv%instid(inst)%pm(k,n) = Atmosphere(1)%pressure(k)
754 iv%instid(inst)%pf(k,n) = Atmosphere(1)%level_pressure(k)
755 iv%instid(inst)%tm(k,n) = Atmosphere(1)%temperature(k)
756 iv%instid(inst)%qm(k,n) = Atmosphere(1)%absorber(k,1)
758 if (use_crtm_kmatrix) then
759 ! T, Q Jacobian
760 do l=1,nchanl
761 iv%instid(inst)%t_jacobian(l,k,n) = Atmosphere_k(l,1)%temperature(k)
762 iv%instid(inst)%q_jacobian(l,k,n) = Atmosphere_k(l,1)%absorber(k,1)
763 end do
764 end if
765 end do
767 if (crtm_cloud) then
768 do k=1,Atmosphere(1)%n_layers
769 iv%instid(inst)%qcw(k,n) = Atmosphere(1)%cloud(1)%water_content(k)
770 iv%instid(inst)%qci(k,n) = Atmosphere(1)%cloud(2)%water_content(k)
771 iv%instid(inst)%qrn(k,n) = Atmosphere(1)%cloud(3)%water_content(k)
772 iv%instid(inst)%qsn(k,n) = Atmosphere(1)%cloud(4)%water_content(k)
773 iv%instid(inst)%qgr(k,n) = Atmosphere(1)%cloud(5)%water_content(k)
774 iv%instid(inst)%qhl(k,n) = Atmosphere(1)%cloud(6)%water_content(k)
775 iv%instid(inst)%rcw(k,n) = Atmosphere(1)%cloud(1)%effective_radius(k)
776 iv%instid(inst)%rci(k,n) = Atmosphere(1)%cloud(2)%effective_radius(k)
777 iv%instid(inst)%rrn(k,n) = Atmosphere(1)%cloud(3)%effective_radius(k)
778 iv%instid(inst)%rsn(k,n) = Atmosphere(1)%cloud(4)%effective_radius(k)
779 iv%instid(inst)%rgr(k,n) = Atmosphere(1)%cloud(5)%effective_radius(k)
780 iv%instid(inst)%rhl(k,n) = Atmosphere(1)%cloud(6)%effective_radius(k)
782 if (use_crtm_kmatrix) then
783 ! Cloud Jacobian
784 do l=1,nchanl
785 iv%instid(inst)%water_jacobian(l,k,n) = Atmosphere_k(l,1)%cloud(1)%water_content(k)
786 iv%instid(inst)%ice_jacobian(l,k,n) = Atmosphere_k(l,1)%cloud(2)%water_content(k)
787 iv%instid(inst)%rain_jacobian(l,k,n) = Atmosphere_k(l,1)%cloud(3)%water_content(k)
788 iv%instid(inst)%snow_jacobian(l,k,n) = Atmosphere_k(l,1)%cloud(4)%water_content(k)
789 iv%instid(inst)%graupel_jacobian(l,k,n) = Atmosphere_k(l,1)%cloud(5)%water_content(k)
790 iv%instid(inst)%hail_jacobian(l,k,n) = Atmosphere_k(l,1)%cloud(6)%water_content(k)
792 iv%instid(inst)%water_r_jacobian(l,k,n) = Atmosphere_k(l,1)%cloud(1)%effective_radius(k)
793 iv%instid(inst)%ice_r_jacobian(l,k,n) = Atmosphere_k(l,1)%cloud(2)%effective_radius(k)
794 iv%instid(inst)%rain_r_jacobian(l,k,n) = Atmosphere_k(l,1)%cloud(3)%effective_radius(k)
795 iv%instid(inst)%snow_r_jacobian(l,k,n) = Atmosphere_k(l,1)%cloud(4)%effective_radius(k)
796 iv%instid(inst)%graupel_r_jacobian(l,k,n)= Atmosphere_k(l,1)%cloud(5)%effective_radius(k)
797 iv%instid(inst)%hail_r_jacobian(l,k,n) = Atmosphere_k(l,1)%cloud(6)%effective_radius(k)
798 end do
799 end if
800 end do
801 end if
803 if ( calc_weightfunc .and. use_crtm_kmatrix ) then
804 do l = 1, nchanl
805 total_od = 0.0
806 do k = 1, Atmosphere(1)%n_layers
807 iv%instid(inst)%lod(l,k,n) = RTSolution(l,1)%layer_optical_depth(k)
808 iv%instid(inst)%lod_jacobian(l,k,n) = RTSolution_K(l,1)%layer_optical_depth(k)
809 total_od = total_od + RTSolution(l,1)%layer_optical_depth(k)
810 iv%instid(inst)%trans(l,k,n) = exp(-total_od)
811 iv%instid(inst)%trans_jacobian(l,k,n) = -iv%instid(inst)%lod_jacobian(l,k,n)*exp(iv%instid(inst)%lod(l,k,n))
812 end do
813 iv%instid(inst)%der_trans(l,1,n) = 0.0
814 do k= 2, Atmosphere(1)%n_layers
815 iv%instid(inst)%der_trans(l,k,n) = abs( (iv%instid(inst)%trans(l,k-1,n)-iv%instid(inst)%trans(l,k,n))/ &
816 (log(iv%instid(inst)%pm(k-1,n))-log(iv%instid(inst)%pm(k,n))) )
817 end do
818 end do
819 end if
821 if (use_clddet==3) then !begin ecmwf cloud detection
822 !kmin_t,kmax_p corresponding tropopause and KBPL
823 ! find the k index of tropopause
824 t_tropp = 999.0 ! initialize
825 k_tropp = kte ! initialize the k index of tropopause
826 do k = kte, kts, -1
827 do j = jts, jte
828 do i = its, ite
829 if ( grid%xb%t(i,j,k) < t_tropp(i,j) .and. &
830 grid%xb%p(i,j,k) > 5000.0 ) then ! tropopause should not be higher than 50 hPa
831 t_tropp(i,j) = grid%xb%t(i,j,k)
832 k_tropp(i,j) = k
833 iv%instid(inst)%kmin_t(n)=Atmosphere(1)%n_layers-k_tropp(i,j)+1
834 end if
835 end do
836 end do
837 end do
838 iv%instid(inst)%kmax_p(n)=Atmosphere(1)%n_layers-grid%kpbl(i,j)+1
839 !write(unit=stdout,fmt='(A)') &
840 ! 'preparing ECMWF cloud detection input using cloud sensitivity to rank channels'
841 do l=1,nchanl
842 temp1=0
843 do k=1,Atmosphere(1)%n_layers !rad_ovc(l,1,n) top
844 iv%instid(inst)%sensitivity_ratio(l,k,n) = abs(iv%instid(inst)%rad_xb(l,n)-iv%instid(inst)%rad_ovc(l,k,n))
845 iv%instid(inst)%sensitivity_ratio(l,k,n) = iv%instid(inst)%sensitivity_ratio(l,k,n)/iv%instid(inst)%rad_xb(l,n)
846 end do
847 end do
848 do l=1,nchanl
849 temp1=0
850 do k=Atmosphere(1)%n_layers,1,-1 ! sensitivity_ratio(l,1,n) top
851 if( iv%instid(inst)%sensitivity_ratio(l,k,n) > 0.01 ) then
852 iv%instid(inst)%p_chan_level(l,n) = real(k) ! the peak level
853 exit
854 end if
855 end do
856 end do
857 end if !end ecmwf cloud detection
859 !----------------------------------------------
860 ! [4.0] store surface information to innovation structure
861 !----------------------------------------------
862 iv%instid(inst)%u10(n) = model_u10(n)
863 iv%instid(inst)%v10(n) = model_v10(n)
864 iv%instid(inst)%t2m(n) = 0.01*missing_r !model_t2m
865 iv%instid(inst)%mr2m(n) = 0.01*missing_r !model_mr2m
866 iv%instid(inst)%ps(n) = model_psfc(n)
867 iv%instid(inst)%ts(n) = model_ts(n)
868 iv%instid(inst)%smois(n) = model_smois(n)
869 iv%instid(inst)%tslb(n) = model_tslb(n)
870 iv%instid(inst)%snowh(n) = model_snowh(n)
871 iv%instid(inst)%isflg(n) = model_isflg
872 iv%instid(inst)%elevation(n) = model_elv(n)
873 iv%instid(inst)%soiltyp(n) = model_isltyp
874 iv%instid(inst)%vegtyp(n) = model_ivgtyp
875 iv%instid(inst)%vegfra(n) = model_vegfra(n)
876 iv%instid(inst)%clwp(n) = clwp
877 if (crtm_cloud) iv%instid(inst)%cip(n) = cip
878 iv%instid(inst)%water_coverage(n) = Surface(1)%water_coverage
879 iv%instid(inst)%land_coverage(n) = Surface(1)%land_coverage
880 iv%instid(inst)%ice_coverage(n) = Surface(1)%ice_coverage
881 iv%instid(inst)%snow_coverage(n) = Surface(1)%snow_coverage
883 if ( calc_tb_clr ) then
884 atmosphere(1)%n_clouds = 0
885 call da_crtm_direct(1, nchanl, 1, Atmosphere, &
886 Surface, &
887 GeometryInfo, &
888 ChannelInfo(inst:inst), &
889 RTSolution_clr, &
890 Options)
891 do k = 1, nchanl
892 iv%instid(inst)%tb_xb_clr(k,n) = RTSolution_clr(k,1)%Brightness_Temperature
893 end do
894 ! resetting n_clouds
895 atmosphere(1)%n_clouds = n_clouds
896 end if !calc_tb_clr
898 end do pixel_loop ! end loop for pixels
900 if ( use_mspps_emis(inst) ) deallocate(em_mspps)
902 deallocate (model_u10)
903 deallocate (model_v10)
904 deallocate (model_psfc)
905 deallocate (model_ts)
906 deallocate (model_tslb)
907 deallocate (model_snowh)
908 deallocate (model_snow)
909 deallocate (model_elv)
910 deallocate (model_vegfra)
911 deallocate (model_smois)
913 call CRTM_RTSolution_Destroy(RTSolution)
914 deallocate( RTSolution, STAT = Allocate_Status )
915 if (Allocate_Status /= 0) &
916 call da_error(__FILE__,__LINE__, &
917 (/"Error in deallocating RTSolution"/))
919 if ( calc_tb_clr ) then
920 call CRTM_RTSolution_Destroy(RTSolution_clr)
921 deallocate( RTSolution_clr, STAT = Allocate_Status )
922 if (Allocate_Status /= 0) &
923 call da_error(__FILE__,__LINE__, &
924 (/"Error in deallocating RTSolution_clr"/))
925 end if !calc_tb_clr
927 if (use_crtm_kmatrix) then
928 call CRTM_RTSolution_Destroy(RTSolution_K)
929 deallocate( RTSolution_K, STAT = Allocate_Status )
930 if (Allocate_Status /= 0) &
931 call da_error(__FILE__,__LINE__, &
932 (/"Error in deallocating RTSolution_K"/))
933 end if
935 call CRTM_Options_Destroy(Options)
936 IF ( CRTM_Options_Associated( Options(1) ) ) &
937 call da_error(__FILE__,__LINE__, &
938 (/"Error in deallocating CRTM Options Structure"/))
940 call CRTM_Surface_Destroy(Surface)
941 IF ( CRTM_Surface_Associated( Surface(1) ) ) &
942 call da_error(__FILE__,__LINE__, &
943 (/"Error in deallocating CRTM Surface Structure"/))
945 if (use_crtm_kmatrix) then
946 call CRTM_Surface_Destroy(Surface_K)
947 IF ( ANY( CRTM_Surface_Associated( Surface_K(1,:)) ) ) &
948 call da_error(__FILE__,__LINE__, &
949 (/"Error in deallocatting CRTM Surface_K Structure"/))
950 endif
952 if (use_crtm_kmatrix) then
953 call CRTM_Atmosphere_Destroy( Atmosphere_K )
954 IF ( ANY( CRTM_Atmosphere_Associated( Atmosphere_K(1,:)) ) ) &
955 call da_error(__FILE__,__LINE__, &
956 (/"Error in deallocatting CRTM Atmosphere_K Structure"/))
957 endif
959 if (use_crtm_kmatrix) then
960 deallocate( Atmosphere_K, Surface_K, STAT = Allocate_Status )
961 if ( Allocate_Status /= 0 ) &
962 call da_error(__FILE__,__LINE__, &
963 (/"Error in deallocatting CRTM Surface_K Structure"/))
964 endif
966 end do ! end loop for sensor
968 deallocate (wrf_to_crtm_mw)
969 if ( use_clddet_zz ) deallocate ( geoht_full )
970 call CRTM_Atmosphere_Destroy (Atmosphere)
971 IF ( CRTM_Atmosphere_Associated( Atmosphere(1) ) ) &
972 call da_error(__FILE__,__LINE__, &
973 (/"Error in deallocating CRTM Atmosphere Structure"/))
975 if (trace_use) call da_trace_exit("da_get_innov_vector_crtm")
977 end subroutine da_get_innov_vector_crtm