da_crtm_init.inc
References to this file elsewhere.
1 subroutine da_crtm_init(iv,ob)
2 !------------------------------------------------------------------------------
3 ! PURPOSE: interface to the initialization subroutine of RTTOV.
4 !
5 ! METHOD: 1. read CRTM coefs files; 2. allocate radiance structure
6 !
7 ! HISTORY: 10/15/2006 added crtm initialization Tomislava Vukicevic, ATOC, University of Colorado
8 ! 11/09/2006 Updated Zhiquan Liu
9 !------------------------------------------------------------------------------
10
11 implicit none
12
13 type (iv_type), intent (inout) :: iv
14 type (y_type) , intent (inout) :: ob
15
16 #ifdef CRTM
17
18 !
19 ! local arguments
20 !-------------------
21 integer :: n, j, ichan
22
23 !
24 ! CRTM local ---------------------------------------------------
25 !
26 integer :: Error_Status
27 ! character( 256 ) :: SpcCoeff_File
28 ! character( 256 ) :: TauCoeff_File
29 character( 256 ) :: AerosolCoeff_File
30 character( 256 ) :: CloudCoeff_File
31 character( 256 ) :: EmisCoeff_File
32 character( 256 ) :: File_Path
33 ! character( 80 ), pointer :: Sensor_Descriptor(:)
34 !
35 ! end of CRTM local
36 !----------------------------------------------------------------
37 ! input parameters of RTTOV_SETUP
38 !----------------------------------
39 ! integer :: err_unit ! Logical error unit (<0 for default)
40 ! integer :: verbosity_level ! (<0 for default)
41 integer :: nsensor, unit_factor_rad
42 ! integer, allocatable :: sensor(:,:) ! instrument id
43
44 ! output parameters of RTTOV_SETUP
45 !-------------------------------------
46 ! integer, allocatable :: errorstatus(:) ! return code
47 ! integer :: coef_errorstatus
48 ! integer :: i_amsua, i_amsub
49 integer :: error
50 integer, allocatable :: nscan(:), nchanl(:)
51
52 ! local variables
53 !----------------
54 ! integer :: nprofiles, nfrequencies, nchannels, nbtout
55 ! integer , pointer :: lprofiles(:)
56 ! real , pointer :: surfem(:)
57 ! real , pointer :: emissivity(:)
58 ! integer , pointer :: channels (:)
59 ! integer , pointer :: polarisations(:,:)
60 ! integer :: mxchn
61 integer :: idum, wmo_sensor_id, sensor_type, iost
62 integer :: iunit
63 character(len=filename_len) :: filename
64
65 ! local variables for tuning error factor
66 character(len=20) :: rttovid_string
67 integer :: num_tot
68 real :: joa, jo, trace, factor
69
70 !--------------------------------------------------------------
71 ! 1.0 setup RTTOV instrument triplets from namelist parameter
72 !--------------------------------------------------------------
73
74 call da_trace_entry("da_crtm_init")
75
76 !--------------------------------------------------------------
77 ! 1.0 setup from namelist parameter
78 !--------------------------------------------------------------
79 nsensor = rtminit_nsensor
80 allocate (nscan(nsensor))
81 allocate (nchanl(nsensor))
82
83 !----------------------------------------------------------------
84 ! 2.0 set up some common varibles for innovation/observation structure
85 !----------------------------------------------------------------
86 iv % num_inst = nsensor
87 ob % num_inst = nsensor
88
89 allocate (iv%instid(1:nsensor))
90 allocate (ob%instid(1:nsensor))
91 allocate (satinfo(1:nsensor))
92
93 iv%instid(1:nsensor)%num_rad = 0
94 ob%instid(1:nsensor)%num_rad = 0
95
96 loop_sensor: do n = 1, nsensor
97
98 iv%instid(n)%platform_id = rtminit_platform(n)
99 iv%instid(n)%satellite_id = rtminit_satid(n)
100 iv%instid(n)%sensor_id = rtminit_sensor(n)
101 if ( rtminit_satid(n) < 10 ) then
102 write(iv%instid(n)%rttovid_string, '(a,i1,a)') &
103 trim( rttov_platform_name(rtminit_platform(n)) )//'-', &
104 rtminit_satid(n), &
105 '-'//trim( rttov_inst_name(rtminit_sensor(n)) )
106 else
107 write(iv%instid(n)%rttovid_string, '(a,i2.2,a)') &
108 trim( rttov_platform_name(rtminit_platform(n)) )//'-', &
109 rtminit_satid(n), &
110 '-'//trim( rttov_inst_name(rtminit_sensor(n)) )
111 end if
112
113 if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'msu' ) then
114 nchanl(n) = 4
115 nscan(n) = 11
116 else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'hirs' ) then
117 nchanl(n) = 19
118 nscan(n) = 56
119 else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'amsua' ) then
120 nchanl(n) = 15
121 nscan(n) = 30
122 else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'amsub' ) then
123 nchanl(n) = 5
124 nscan(n) = 90
125 else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'airs' ) then
126 nchanl(n) = 281
127 nscan(n) = 90
128 else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'hsb' ) then
129 nchanl(n) = 4
130 nscan(n) = 90
131 else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'mhs' ) then
132 nchanl(n) = 5
133 nscan(n) = 90
134 else if ( trim( crtm_sensor_name(rtminit_sensor(n))) == 'ssmis' ) then
135 nchanl(n) = 24
136 nscan(n) = 180
137 else
138 call da_error(__FILE__,__LINE__, &
139 (/"Unrecorgnized instrument"/))
140 end if
141
142 iv%instid(n)%nchan = nchanl(n)
143 ob%instid(n)%nchan = nchanl(n)
144
145 allocate ( iv%instid(n)%ichan(1:nchanl(n)), stat = error )
146 if( error /= 0 ) then
147 call da_error(__FILE__,__LINE__, &
148 (/"Memory allocation error to iv%instid(n)%ichan"/))
149 end if
150
151 allocate ( ob%instid(n)%ichan(1:nchanl(n)), stat = error )
152 if( error /= 0 ) then
153 call da_error(__FILE__,__LINE__, &
154 (/"Memory allocation error to ob%instid(n)%ichan"/))
155 end if
156
157 call da_get_unit(iunit)
158 filename='radiance_info/'//adjustl(trim(iv%instid(n)%rttovid_string))//'.info'
159 open(unit=iunit,file=filename, form='formatted',iostat = iost, status='old')
160 if (iost /= 0) then
161 message(1)="Cannot open radiance info file "//adjustl(filename)
162 call da_error(__FILE__,__LINE__,message(1:1))
163 end if
164
165 allocate ( satinfo(n) % ichan(nchanl(n)) )
166 allocate ( satinfo(n) % iuse (nchanl(n)) )
167 allocate ( satinfo(n) % error(nchanl(n)) )
168 allocate ( satinfo(n) % polar(nchanl(n)) )
169
170 read(iunit,*)
171 do j = 1, nchanl(n)
172 read(iunit,'(1x,5i5,2e18.10)') &
173 wmo_sensor_id, &
174 satinfo(n)%ichan(j), &
175 sensor_type, &
176 satinfo(n)%iuse(j) , &
177 idum, &
178 satinfo(n)%error(j), &
179 satinfo(n)%polar(j)
180 iv%instid(n)%ichan(j) = satinfo(n)%ichan(j)
181 ob%instid(n)%ichan(j) = satinfo(n)%ichan(j)
182 end do
183 call da_free_unit(iunit)
184
185 end do loop_sensor
186
187 !---------------------------------------------------------------------
188 ! 3.0 get CRTM sensor descriptor to use in CRTM_Set_ChannelInfo later
189 !---------------------------------------------------------------------
190 allocate(Sensor_Descriptor(nsensor))
191 call da_crtm_sensor_descriptor(nsensor,Sensor_Descriptor)
192 allocate(ChannelInfo(nsensor))
193
194 ! CRTM load coefficients
195 !-----------------------------------------------------------
196 ! 3.1 call CRTM_Init to load coefficients and fill ChannelInfo structure
197 !-----------------------------------------------------------
198 ! input:
199 !SpcCoeff_File = 'amsua_n15-n16.SpcCoeff.bin'
200 !TauCoeff_File = 'amsua_n15-n16.TauCoeff.bin'
201 AerosolCoeff_File = 'AerosolCoeff.bin'
202 CloudCoeff_File = 'CloudCoeff.bin'
203 EmisCoeff_File = 'EmisCoeff.bin'
204 File_Path = 'crtm_coeffs/'
205 !----------------------------------------------------------------
206 ! ChannelInfo structure contains on output:
207 !
208 ! n_channels - integer, total number of channels
209 ! Sensor_Index - integer
210 ! Channel_Index - integer pointer, index of the channels loaded during initialization
211 ! Sensor_Channel - integer pointer, the sensor channel #
212 ! SensorID - character pointer, character string containing satellite and sensor descr
213 ! example: amsre_aqua (Appendix B in User guide)
214 ! WMO_Satellite_ID - integer pointer
215 ! WMO_Sensor_ID - integer pointer
216 !-----------------------------------------------------------------
217
218 Error_Status = CRTM_Init( ChannelInfo, &
219 SensorID = Sensor_Descriptor, &
220 AerosolCoeff_File = AerosolCoeff_File, &
221 EmisCoeff_File = EmisCoeff_File , &
222 CloudCoeff_File = CloudCoeff_File, &
223 File_Path = File_Path)
224
225 ! Error_Status = CRTM_Set_ChannelInfo(Sensor_Descriptor(1),ChannelInfo)
226
227 if ( Error_Status /= 0 ) then
228 call da_error(__FILE__,__LINE__, &
229 (/"Error in initializing CRTM"/))
230 END IF
231
232 ! do n = 1, nsensor
233 ! write (6,*) 'in da_crtm_init: ChannelInfo content'
234 ! write (6,*) 'Sensor_Index ',ChannelInfo(n)%Sensor_Index
235 ! write (6,*) 'n_channels ',ChannelInfo(n)%n_channels
236 ! write (6,*) 'Channel_Index ',ChannelInfo(n)%Channel_Index(:)
237 ! write (6,*) 'Sensor_Channel ',ChannelInfo(n)%Sensor_Channel(:)
238 ! write (6,*) 'SensorID ',ChannelInfo(n)%SensorID
239 ! end do
240
241 !-------------------------------------------------------
242 ! 4.0 read bias correction coefs files
243 !-------------------------------------------------------
244
245 if (read_biascoef) then
246 loop_sensor2: do n = 1, nsensor
247
248 allocate ( satinfo(n) % scanbias (nchanl(n),nscan(n)) )
249 allocate ( satinfo(n) % scanbias_b(nchanl(n),nscan(n),18) )
250 allocate ( satinfo(n) % bcoef (nchanl(n),4) )
251 allocate ( satinfo(n) % bcoef0 (nchanl(n)) )
252 allocate ( satinfo(n) % error_std (nchanl(n)) )
253
254 satinfo(n) % error_std(:) = 500.0
255 satinfo(n) % scanbias(:,:) = 0.0
256 satinfo(n) % scanbias_b(:,:,:) = 0.0
257 satinfo(n) % bcoef(:,:) = 0.0
258 satinfo(n) % bcoef0(:) = 0.0
259
260 if ( index(iv%instid(n)%rttovid_string,'eos') > 0 ) cycle ! not implemented
261 if ( index(iv%instid(n)%rttovid_string,'hirs') > 0 ) cycle ! not implemented
262 if ( index(iv%instid(n)%rttovid_string,'ssmis') > 0) cycle ! not implemented
263 ! if ( index(iv%instid(n)%rttovid_string,'mhs') > 0) cycle ! not implemented
264
265 ! new bias coefs files
266 !----------------------------------
267 call da_read_biascoef(iv%instid(n)%rttovid_string, &
268 nchanl(n),nscan(n),18,4,global, &
269 satinfo(n)%scanbias, &
270 satinfo(n)%scanbias_b, &
271 satinfo(n)%bcoef, &
272 satinfo(n)%bcoef0, &
273 satinfo(n)%error_std)
274
275 end do loop_sensor2
276 end if
277
278 !-------------------------------------------------------
279 ! 5.0 read error factor file
280 !-------------------------------------------------------
281 if (use_error_factor_rad) then
282
283 do n=1,rtminit_nsensor
284 allocate ( satinfo(n)%error_factor(1:nchanl(n)) )
285 satinfo(n)%error_factor(:) = 1.0
286 end do
287
288 call da_get_unit(unit_factor_rad)
289 open(unit_factor_rad, file='radiance_error.factor', &
290 form='formatted',iostat = iost, status='old')
291
292 if (iost /= 0) then
293 call da_error(__FILE__,__LINE__, &
294 (/"Cannot open radiance error factor file: radiance_error.factor"/))
295 end if
296
297 read(unit_factor_rad, *)
298 do
299 read(unit_factor_rad,fmt='(a15,i3,i8,3f15.5,f8.3)',iostat=iost) &
300 rttovid_string, ichan, num_tot, joa,jo,trace,factor
301 if ( iost == 0 ) then
302 do n=1,rtminit_nsensor
303 if ( index(rttovid_string,trim(iv%instid(n)%rttovid_string))>0 ) then
304 satinfo(n)%error_factor(ichan) = factor
305 write(6,'(a,i5,a,f10.3)') trim(rttovid_string)//' Channel ', ichan, ' Error Factor = ', factor
306 exit
307 end if
308 end do
309 else
310 exit
311 end if
312 end do
313 close(unit_factor_rad)
314 call da_free_unit(unit_factor_rad)
315
316 end if
317
318 !---------------------------------------------
319 ! 6.0 get FGAT time slots
320 !---------------------------------------------
321 allocate ( time_slots(0:num_fgat_time) )
322 call da_get_time_slots(num_fgat_time,time_window_min,analysis_date, &
323 time_window_max, time_slots)
324
325
326 deallocate(nscan)
327 deallocate(nchanl)
328
329 call da_trace_exit("da_crtm_init")
330 #endif
331 end subroutine da_crtm_init