da_rttov_init.inc
References to this file elsewhere.
1 subroutine da_rttov_init(iv,ob)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: interface to the initialization subroutine of RTTOV8_5.
5 !
6 ! METHOD: 1. read RTTOV coefs files; 2. allocate radiance structure
7 !
8 ! HISTORY: 07/28/2005 - Creation Zhiquan Liu
9 ! 03/22/2006 add error tuning factor read in Zhiquan Liu
10 !
11 !---------------------------------------------------------------------------
12
13 implicit none
14
15 type (ob_type), intent (inout) :: iv
16 type (y_type) , intent (inout) :: ob
17
18 #ifdef RTTOV
19 #include "rttov_setup.interface"
20 #include "rttov_readscattcoeffs.interface"
21
22 ! local arguments
23 !-------------------
24 integer :: i, k,j,bunit, iscan,ichan,npos,unit_factor_rad
25 real :: ave
26
27 !
28 ! input parameters of RTTOV_SETUP
29 !----------------------------------
30 integer :: err_unit ! Logical error unit (<0 for default)
31 integer :: verbosity_level ! (<0 for default)
32 integer :: nsensor
33 integer, allocatable :: sensor(:,:) ! instrument id
34 ! integer, allocatable :: channels(:,:) ! channel list per instrument
35
36 ! output parameters of RTTOV_SETUP
37 !-------------------------------------
38 integer, allocatable :: errorstatus(:) ! return code
39 integer :: coef_errorstatus, i_amsua, i_amsub
40 ! type( rttov_coef ), allocatable :: coefs(:) ! coefficients
41 integer , allocatable :: nscan(:), nchan(:)
42
43 ! local variables
44 !----------------
45 integer :: nprofiles, nfrequencies, nchannels, nbtout
46 integer , allocatable :: lprofiles(:)
47 real , allocatable :: surfem(:), emissivity(:)
48 integer , allocatable :: channels (:), polarisations(:,:)
49 integer :: mxchn
50 integer , allocatable :: coefs_channels (:,:)
51
52 integer :: idum, wmo_sensor_id, sensor_type, iost, iunit
53
54 ! local variables for tuning error factor
55 character(len=20) :: rttovid_string
56 integer :: num_tot
57 real :: joa, jo, trace, factor
58 character(len=filename_len) :: filename
59
60 if (trace_use) call da_trace_entry("da_rttov_init")
61
62 !--------------------------------------------------------------
63 ! 1.0 setup RTTOV instrument triplets from namelist parameter
64 !--------------------------------------------------------------
65
66 mxchn = 300
67 err_unit = stderr
68 verbosity_level = rtminit_print
69 nsensor = rtminit_nsensor
70
71 allocate (errorstatus(nsensor))
72 allocate (coefs(nsensor))
73 allocate (sensor(3,nsensor))
74 allocate (coefs_channels(mxchn,nsensor))
75 allocate (nscan(nsensor))
76 allocate (nchan(nsensor))
77
78 sensor (1,1:nsensor) = rtminit_platform (1:nsensor)
79 sensor (2,1:nsensor) = rtminit_satid (1:nsensor)
80 sensor (3,1:nsensor) = rtminit_sensor (1:nsensor)
81 coefs_channels(:,:) = 0
82
83 if (print_detail_radiance) then
84 write(unit=message(1),fmt='(A,I5)') 'err_unit = ', err_unit
85 write(unit=message(2),fmt='(A,I5)') 'verbosity_level = ', verbosity_level
86 write(unit=message(3),fmt='(A,I5)') 'nsensor = ', nsensor
87 write(unit=message(4),fmt='(A,10I5)') 'sensor (1,1:nsensor) = ', sensor (1,1:nsensor)
88 write(unit=message(5),fmt='(A,10I5)') 'sensor (2,1:nsensor) = ', sensor (2,1:nsensor)
89 write(unit=message(6),fmt='(A,10I5)') 'sensor (3,1:nsensor) = ', sensor (3,1:nsensor)
90 call da_message(message(1:6))
91 end if
92
93 !----------------------------------------------------------------
94 ! 2.0 set up some common variables for innovation/observation structure
95 !----------------------------------------------------------------
96
97 iv % num_inst = nsensor
98 ob % num_inst = nsensor
99
100 allocate (iv%instid(1:nsensor))
101 allocate (ob%instid(1:nsensor))
102 allocate (satinfo(1:nsensor))
103
104 iv%instid(1:nsensor)%num_rad = 0
105 ob%instid(1:nsensor)%num_rad = 0
106
107 loop_sensor: do i = 1, nsensor
108
109 iv%instid(i)%platform_id = rtminit_platform(i)
110 iv%instid(i)%satellite_id = rtminit_satid(i)
111 iv%instid(i)%sensor_id = rtminit_sensor(i)
112 if ( rtminit_satid(i) < 10 ) then
113 write(iv%instid(i)%rttovid_string, '(a,i1,a)') &
114 trim( platform_name(rtminit_platform(i)) )//'-', &
115 rtminit_satid(i), &
116 '-'//trim( inst_name(rtminit_sensor(i)) )
117 else
118 write(iv%instid(i)%rttovid_string, '(a,i2.2,a)') &
119 trim( platform_name(rtminit_platform(i)) )//'-', &
120 rtminit_satid(i), &
121 '-'//trim( inst_name(rtminit_sensor(i)) )
122 end if
123
124 if ( trim( inst_name(rtminit_sensor(i))) == 'msu' ) then
125 nchan(i) = 4
126 nscan(i) = 11
127 else if ( trim( inst_name(rtminit_sensor(i))) == 'hirs' ) then
128 nchan(i) = 19
129 nscan(i) = 56
130 else if ( trim( inst_name(rtminit_sensor(i))) == 'amsua' ) then
131 nchan(i) = 15
132 nscan(i) = 30
133 else if ( trim( inst_name(rtminit_sensor(i))) == 'amsub' ) then
134 nchan(i) = 5
135 nscan(i) = 90
136 else if ( trim( inst_name(rtminit_sensor(i))) == 'airs' ) then
137 nchan(i) = 281
138 nscan(i) = 90
139 else if ( trim( inst_name(rtminit_sensor(i))) == 'hsb' ) then
140 nchan(i) = 4
141 nscan(i) = 90
142 else
143 call da_error(__FILE__,__LINE__,(/"Unrecognized instrument"/))
144 end if
145
146 iv%instid(i)%nchan = nchan(i)
147 ob%instid(i)%nchan = nchan(i)
148
149 allocate (iv%instid(i)%ichan(1:nchan(i)))
150 allocate (ob%instid(i)%ichan(1:nchan(i)))
151
152 call da_get_unit(iunit)
153 filename='wrfvar_run/'//adjustl(trim(iv%instid(i)%rttovid_string))//'.info'
154 open(unit=iunit,file=filename, &
155 form='formatted',iostat = iost, status='old')
156
157 if (iost /= 0) then
158 message(1)="Cannot open radiance info file "//adjustl(filename)
159 call da_error(__FILE__,__LINE__,message(1:1))
160 end if
161
162 allocate (satinfo(i) % ichan(nchan(i)))
163 allocate (satinfo(i) % iuse (nchan(i)))
164 allocate (satinfo(i) % error(nchan(i)))
165 allocate (satinfo(i) % polar(nchan(i)))
166
167 read(unit=iunit,fmt=*)
168 do j = 1, nchan(i)
169 read(unit=iunit,fmt='(1x,5i5,2e18.10)') &
170 wmo_sensor_id, &
171 satinfo(i)%ichan(j), &
172 sensor_type, &
173 satinfo(i)%iuse(j) , &
174 idum, &
175 satinfo(i)%error(j), &
176 satinfo(i)%polar(j)
177 iv%instid(i)%ichan(j) = satinfo(i)%ichan(j)
178 ob%instid(i)%ichan(j) = satinfo(i)%ichan(j)
179 ! only load coefs of selected channels for AIRS
180 coefs_channels(j,i) = satinfo(i)%ichan(j)
181 end do
182
183 close(iunit)
184 call da_free_unit(iunit)
185
186 ! else
187 ! do j = 1, nchan(i)
188 ! iv%instid(i)%ichan(j) = j
189 ! ob%instid(i)%ichan(j) = j
190 ! coefs_channels(j,n) = j
191 ! end do
192 ! end if
193
194 end do loop_sensor
195
196 !-----------------------------------------------------------
197 ! 3.0 call rttov_setup for reading clear sky coefficients
198 !-----------------------------------------------------------
199
200 call rttov_setup (&
201 errorstatus, &! out
202 err_unit, &! in
203 verbosity_level, &! in
204 nsensor, &! in
205 coefs, &! out
206 sensor ) !, &! in
207 ! coefs_channels ) ! in Optional not work with amsua+b
208
209 if (any(errorstatus(:) /= errorstatus_success)) then
210 call da_error(__FILE__,__LINE__,(/"rttov_setup fatal error"/))
211 end if
212
213 !-------------------------------------------------------------
214 ! 4.0 read coef file for cloud/rain absorption/scattering
215 !-------------------------------------------------------------
216
217 if (rttov_scatt) then
218 i_amsua = 0
219 i_amsub = 0
220 do i=1,nsensor
221 if (trim(inst_name(rtminit_sensor(i))) == 'amsua') &
222 i_amsua = i
223 if ( trim(inst_name(rtminit_sensor(i))) == 'amsub') &
224 i_amsub = i
225 end do
226
227 if (i_amsua /= 0 .and. i_amsub == 0) then
228 n_scatt_coef = 1
229 allocate (coefs_scatt(n_scatt_coef))
230 allocate (coefs_scatt_instname(n_scatt_coef))
231 coefs_scatt_instname(1) = 'amsua'
232 call rttov_readscattcoeffs(coef_errorstatus, coefs(i_amsua), &
233 coefs_scatt(1))
234 if (coef_errorstatus /= errorstatus_success ) then
235 call da_error(__FILE__,__LINE__,(/"rttov_readscattcoeffs fatal error"/))
236 end if
237 end if
238 if (i_amsua == 0 .and. i_amsub /= 0) then
239 n_scatt_coef = 1
240 allocate (coefs_scatt(n_scatt_coef))
241 allocate (coefs_scatt_instname(n_scatt_coef))
242 coefs_scatt_instname(1) = 'amsub'
243 call rttov_readscattcoeffs(coef_errorstatus, coefs(i_amsub), &
244 coefs_scatt(1))
245 if (coef_errorstatus /= errorstatus_success ) then
246 call da_error(__FILE__,__LINE__,(/"rttov_readscattcoeffs fatal error"/))
247 end if
248 end if
249
250 if (i_amsua /= 0 .and. i_amsub /= 0) then
251 n_scatt_coef = 2
252 allocate (coefs_scatt(n_scatt_coef))
253 allocate (coefs_scatt_instname(n_scatt_coef))
254 coefs_scatt_instname(1) = 'amsua'
255 call rttov_readscattcoeffs(coef_errorstatus, coefs(i_amsua), &
256 coefs_scatt(1))
257 if (coef_errorstatus /= errorstatus_success ) then
258 call da_error(__FILE__,__LINE__,(/"rttov_readscattcoeffs fatal error"/))
259 end if
260 coefs_scatt_instname(2) = 'amsub'
261 call rttov_readscattcoeffs(coef_errorstatus, coefs(i_amsub), &
262 coefs_scatt(2))
263 if (coef_errorstatus /= errorstatus_success ) then
264 call da_error(__FILE__,__LINE__,(/"rttov_readscattcoeffs fatal error"/))
265 end if
266 end if
267
268 if ( i_amsua == 0 .and. i_amsub == 0 ) n_scatt_coef = 0
269 end if
270
271 !----------------------------------------------------------------
272 ! 5.0 set up some common varibles for innovation/observation structure
273 !----------------------------------------------------------------
274
275 iv%instid(1:nsensor)%nlevels = coefs(1:nsensor)%nlevels
276
277 call da_message((/"Read in the RTTOV coef files for the following sensors"/))
278
279 loop_sensor2: do i = 1, nsensor
280
281 !---------------------------------------------------
282 ! 5.1 get more information about sensor
283 !---------------------------------------------------
284
285 nprofiles = 1
286 write(unit=message(1),fmt='(3x,2a,2x,a,i5,2x,a,i5)') " ", &
287 iv%instid(i)%rttovid_string, 'nchan=',nchan(i), 'nscan=',nscan(i)
288
289 call rttov_setupchan(nprofiles, nchan(i), coefs(i), & ! in
290 nfrequencies, nchannels, nbtout ) ! out
291
292 allocate (lprofiles(nfrequencies))
293 allocate (channels (nfrequencies))
294 allocate (polarisations(nchannels, 3))
295 allocate (emissivity( nchannels))
296 allocate (surfem ( nchannels))
297 surfem(:) = 0.
298
299 call rttov_setupindex(nchan(i), nprofiles, nfrequencies, & ! in
300 nchannels, nbtout, coefs(i), surfem, & ! in
301 lprofiles, channels, polarisations, & ! out
302 emissivity )
303
304 write (unit=message(2),fmt=*) &
305 ' nprofiles = ', nprofiles, &
306 'nchan = ', nchan(i), &
307 'nfrequencies = ', nfrequencies, &
308 'nchannels = ', nchannels, &
309 'nbtout = ', nbtout
310 ! JRB problems with PGI
311 ! write (unit=message(3),fmt=*) ' lprofiles = ', lprofiles
312 ! write (unit=message(4),fmt=*) ' channels = ', channels
313
314 call da_message(message(1:2))
315
316 deallocate (lprofiles)
317 deallocate (channels)
318 deallocate (polarisations)
319 deallocate (emissivity)
320 deallocate (surfem)
321
322 !-------------------------------------------------------
323 ! 6.0 read bias correction coefs files
324 !-------------------------------------------------------
325
326 if ( read_biascoef .and. &
327 ( trim(iv%instid(i)%rttovid_string) /= 'eos-2-airs' ) &
328 .and. index(iv%instid(i)%rttovid_string,'hirs') <= 0 ) then
329
330 allocate ( satinfo(i) % scanbias (nchan(i),nscan(i)) )
331 allocate ( satinfo(i) % scanbias_b(nchan(i),nscan(i),18) )
332 allocate ( satinfo(i) % bcoef (nchan(i),4) )
333 allocate ( satinfo(i) % bcoef0 (nchan(i)) )
334 allocate ( satinfo(i) % error_std (nchan(i)) )
335
336 satinfo(i) % error_std(:) = 500.
337 satinfo(i) % scanbias(:,:) = 0.
338 satinfo(i) % scanbias_b(:,:,:) = 0.
339 satinfo(i) % bcoef(:,:) = 0.
340 satinfo(i) % bcoef0(:) = 0.
341
342 ! new bias coefs files
343 !----------------------------------
344 call da_read_biascoef(iv%instid(i)%rttovid_string, &
345 nchan(i),nscan(i),18,4,global, &
346 satinfo(i)%scanbias, &
347 satinfo(i)%scanbias_b, &
348 satinfo(i)%bcoef, &
349 satinfo(i)%bcoef0, &
350 satinfo(i)%error_std)
351
352 end if
353
354 end do loop_sensor2
355
356 if (use_error_factor_rad) then
357 do i=1,rtminit_nsensor
358 allocate ( satinfo(i)%error_factor(1:nchan(i)) )
359 satinfo(i)%error_factor(:) = 1.0
360 end do
361
362 call da_get_unit(unit_factor_rad)
363 filename='wrfvar_run/radiance_error.factor'
364 open(unit_factor_rad, file=filename, &
365 form='formatted',iostat = iost, status='old')
366
367
368 if (iost /= 0) then
369 call da_error(__FILE__,__LINE__, &
370 (/"Cannot open radiance error factor file "//adjustl(filename)/))
371 end if
372
373 read(unit_factor_rad, *)
374 do
375 read(unit_factor_rad,fmt='(a15,i3,i8,3f15.5,f8.3)',iostat=iost) &
376 rttovid_string, ichan, num_tot, joa,jo,trace,factor
377 if ( iost == 0 ) then
378 do i=1,rtminit_nsensor
379 if ( trim(rttovid_string)== trim(iv%instid(i)%rttovid_string) ) then
380 satinfo(i)%error_factor(ichan) = factor
381 if (print_detail_radiance) &
382 write(unit=stdout,fmt='(a,i5,a,f10.3)') &
383 trim(rttovid_string)//' Channel ', ichan, &
384 ' Error Factor = ', factor
385 exit
386 end if
387 end do
388 else
389 exit
390 end if
391 end do
392 close(unit_factor_rad)
393 call da_free_unit(unit_factor_rad)
394 end if
395
396 deallocate (errorstatus)
397 deallocate (sensor)
398 deallocate (coefs_channels)
399 deallocate (nscan)
400 deallocate (nchan)
401
402 allocate (time_slots(0:num_fgat_time))
403 call da_get_time_slots(num_fgat_time,time_window_min,analysis_date, &
404 time_window_max, time_slots)
405
406 if (trace_use) call da_trace_exit("da_rttov_init")
407
408 #endif
409
410 end subroutine da_rttov_init