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 (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_rad) 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 end do loop_sensor
187
188 !-----------------------------------------------------------
189 ! 3.0 call rttov_setup for reading clear sky coefficients
190 !-----------------------------------------------------------
191
192 call rttov_setup (&
193 errorstatus, &! out
194 err_unit, &! in
195 verbosity_level, &! in
196 nsensor, &! in
197 coefs, &! out
198 sensor ) !, &! in
199 ! coefs_channels ) ! in Optional not work with amsua+b
200
201 if (any(errorstatus(:) /= errorstatus_success)) then
202 call da_error(__FILE__,__LINE__,(/"rttov_setup fatal error"/))
203 end if
204
205 !-------------------------------------------------------------
206 ! 4.0 read coef file for cloud/rain absorption/scattering
207 !-------------------------------------------------------------
208
209 if (rttov_scatt) then
210 i_amsua = 0
211 i_amsub = 0
212 do i=1,nsensor
213 if (trim(inst_name(rtminit_sensor(i))) == 'amsua') &
214 i_amsua = i
215 if ( trim(inst_name(rtminit_sensor(i))) == 'amsub') &
216 i_amsub = i
217 end do
218
219 if (i_amsua /= 0 .and. i_amsub == 0) then
220 n_scatt_coef = 1
221 allocate (coefs_scatt(n_scatt_coef))
222 allocate (coefs_scatt_instname(n_scatt_coef))
223 coefs_scatt_instname(1) = 'amsua'
224 call rttov_readscattcoeffs(coef_errorstatus, coefs(i_amsua), &
225 coefs_scatt(1))
226 if (coef_errorstatus /= errorstatus_success ) then
227 call da_error(__FILE__,__LINE__,(/"rttov_readscattcoeffs fatal error"/))
228 end if
229 end if
230 if (i_amsua == 0 .and. i_amsub /= 0) then
231 n_scatt_coef = 1
232 allocate (coefs_scatt(n_scatt_coef))
233 allocate (coefs_scatt_instname(n_scatt_coef))
234 coefs_scatt_instname(1) = 'amsub'
235 call rttov_readscattcoeffs(coef_errorstatus, coefs(i_amsub), &
236 coefs_scatt(1))
237 if (coef_errorstatus /= errorstatus_success ) then
238 call da_error(__FILE__,__LINE__,(/"rttov_readscattcoeffs fatal error"/))
239 end if
240 end if
241
242 if (i_amsua /= 0 .and. i_amsub /= 0) then
243 n_scatt_coef = 2
244 allocate (coefs_scatt(n_scatt_coef))
245 allocate (coefs_scatt_instname(n_scatt_coef))
246 coefs_scatt_instname(1) = 'amsua'
247 call rttov_readscattcoeffs(coef_errorstatus, coefs(i_amsua), &
248 coefs_scatt(1))
249 if (coef_errorstatus /= errorstatus_success ) then
250 call da_error(__FILE__,__LINE__,(/"rttov_readscattcoeffs fatal error"/))
251 end if
252 coefs_scatt_instname(2) = 'amsub'
253 call rttov_readscattcoeffs(coef_errorstatus, coefs(i_amsub), &
254 coefs_scatt(2))
255 if (coef_errorstatus /= errorstatus_success ) then
256 call da_error(__FILE__,__LINE__,(/"rttov_readscattcoeffs fatal error"/))
257 end if
258 end if
259
260 if ( i_amsua == 0 .and. i_amsub == 0 ) n_scatt_coef = 0
261 end if
262
263 !----------------------------------------------------------------
264 ! 5.0 set up some common varibles for innovation/observation structure
265 !----------------------------------------------------------------
266
267 iv%instid(1:nsensor)%nlevels = coefs(1:nsensor)%nlevels
268
269 call da_message((/"Read in the RTTOV coef files for the following sensors"/))
270
271 loop_sensor2: do i = 1, nsensor
272
273 !---------------------------------------------------
274 ! 5.1 get more information about sensor
275 !---------------------------------------------------
276
277 nprofiles = 1
278 write(unit=message(1),fmt='(3x,2a,2x,a,i5,2x,a,i5)') " ", &
279 iv%instid(i)%rttovid_string, 'nchan=',nchan(i), 'nscan=',nscan(i)
280
281 call rttov_setupchan(nprofiles, nchan(i), coefs(i), & ! in
282 nfrequencies, nchannels, nbtout ) ! out
283
284 allocate (lprofiles(nfrequencies))
285 allocate (channels (nfrequencies))
286 allocate (polarisations(nchannels, 3))
287 allocate (emissivity( nchannels))
288 allocate (surfem ( nchannels))
289 surfem(:) = 0.
290
291 call rttov_setupindex(nchan(i), nprofiles, nfrequencies, & ! in
292 nchannels, nbtout, coefs(i), surfem, & ! in
293 lprofiles, channels, polarisations, & ! out
294 emissivity )
295
296 write (unit=message(2),fmt=*) &
297 ' nprofiles = ', nprofiles, &
298 'nchan = ', nchan(i), &
299 'nfrequencies = ', nfrequencies, &
300 'nchannels = ', nchannels, &
301 'nbtout = ', nbtout
302 ! JRB problems with PGI
303 ! write (unit=message(3),fmt=*) ' lprofiles = ', lprofiles
304 ! write (unit=message(4),fmt=*) ' channels = ', channels
305
306 call da_message(message(1:2))
307
308 deallocate (lprofiles)
309 deallocate (channels)
310 deallocate (polarisations)
311 deallocate (emissivity)
312 deallocate (surfem)
313
314 !-------------------------------------------------------
315 ! 6.0 read bias correction coefs files
316 !-------------------------------------------------------
317
318 if ( read_biascoef ) then
319
320 allocate ( satinfo(i) % scanbias (nchan(i),nscan(i)) )
321 allocate ( satinfo(i) % scanbias_b(nchan(i),nscan(i),18) )
322 allocate ( satinfo(i) % bcoef (nchan(i),4) )
323 allocate ( satinfo(i) % bcoef0 (nchan(i)) )
324 allocate ( satinfo(i) % error_std (nchan(i)) )
325
326 satinfo(i) % error_std(:) = 500.
327 satinfo(i) % scanbias(:,:) = 0.
328 satinfo(i) % scanbias_b(:,:,:) = 0.
329 satinfo(i) % bcoef(:,:) = 0.
330 satinfo(i) % bcoef0(:) = 0.
331
332 if ( index(iv%instid(i)%rttovid_string,'eos') > 0 ) cycle ! not implemented
333 if ( index(iv%instid(i)%rttovid_string,'hirs') > 0 ) cycle ! not implemented
334 if ( index(iv%instid(i)%rttovid_string,'ssmis') > 0) cycle ! not implemented
335
336 ! new bias coefs files
337 !----------------------------------
338 call da_read_biascoef(iv%instid(i)%rttovid_string, &
339 nchan(i),nscan(i),18,4,global, &
340 satinfo(i)%scanbias, &
341 satinfo(i)%scanbias_b, &
342 satinfo(i)%bcoef, &
343 satinfo(i)%bcoef0, &
344 satinfo(i)%error_std)
345
346 end if
347
348 end do loop_sensor2
349
350 if (use_error_factor_rad) then
351 do i=1,rtminit_nsensor
352 allocate ( satinfo(i)%error_factor(1:nchan(i)) )
353 satinfo(i)%error_factor(:) = 1.0
354 end do
355
356 call da_get_unit(unit_factor_rad)
357 filename='wrfvar_run/radiance_error.factor'
358 open(unit_factor_rad, file=filename, &
359 form='formatted',iostat = iost, status='old')
360
361
362 if (iost /= 0) then
363 call da_error(__FILE__,__LINE__, &
364 (/"Cannot open radiance error factor file "//adjustl(filename)/))
365 end if
366
367 read(unit_factor_rad, *)
368 do
369 read(unit_factor_rad,fmt='(a15,i3,i8,3f15.5,f8.3)',iostat=iost) &
370 rttovid_string, ichan, num_tot, joa,jo,trace,factor
371 if ( iost == 0 ) then
372 do i=1,rtminit_nsensor
373 if ( index(rttovid_string,trim(iv%instid(i)%rttovid_string))>0 ) then
374 satinfo(i)%error_factor(ichan) = factor
375 write(unit=stdout,fmt='(a,i5,a,f10.3)') &
376 trim(rttovid_string)//' Channel ', ichan, &
377 ' Error Factor = ', factor
378 exit
379 end if
380 end do
381 else
382 exit
383 end if
384 end do
385 close(unit_factor_rad)
386 call da_free_unit(unit_factor_rad)
387 end if
388
389 deallocate (errorstatus)
390 deallocate (sensor)
391 deallocate (coefs_channels)
392 deallocate (nscan)
393 deallocate (nchan)
394
395 allocate (time_slots(0:num_fgat_time))
396 call da_get_time_slots(num_fgat_time,time_window_min,analysis_date, &
397 time_window_max, time_slots)
398
399 if (trace_use) call da_trace_exit("da_rttov_init")
400
401 #endif
402
403 end subroutine da_rttov_init