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