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 (ob_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))) == 'ssmis' ) then
132 nchanl(n) = 24
133 nscan(n) = 180
134 else
135 call da_error(__FILE__,__LINE__, &
136 (/"Unrecorgnized instrument"/))
137 end if
138
139 iv%instid(n)%nchan = nchanl(n)
140 ob%instid(n)%nchan = nchanl(n)
141
142 allocate ( iv%instid(n)%ichan(1:nchanl(n)), stat = error )
143 if( error /= 0 ) then
144 call da_error(__FILE__,__LINE__, &
145 (/"Memory allocation error to iv%instid(n)%ichan"/))
146 end if
147
148 allocate ( ob%instid(n)%ichan(1:nchanl(n)), stat = error )
149 if( error /= 0 ) then
150 call da_error(__FILE__,__LINE__, &
151 (/"Memory allocation error to ob%instid(n)%ichan"/))
152 end if
153
154 call da_get_unit(iunit)
155 filename='wrfvar_run/'//adjustl(trim(iv%instid(n)%rttovid_string))//'.info'
156 open(unit=iunit,file=filename, form='formatted',iostat = iost, status='old')
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(n) % ichan(nchanl(n)) )
163 allocate ( satinfo(n) % iuse (nchanl(n)) )
164 allocate ( satinfo(n) % error(nchanl(n)) )
165 allocate ( satinfo(n) % polar(nchanl(n)) )
166
167 read(iunit,*)
168 do j = 1, nchanl(n)
169 read(iunit,'(1x,5i5,2e18.10)') &
170 wmo_sensor_id, &
171 satinfo(n)%ichan(j), &
172 sensor_type, &
173 satinfo(n)%iuse(j) , &
174 idum, &
175 satinfo(n)%error(j), &
176 satinfo(n)%polar(j)
177 iv%instid(n)%ichan(j) = satinfo(n)%ichan(j)
178 ob%instid(n)%ichan(j) = satinfo(n)%ichan(j)
179 end do
180 call da_free_unit(iunit)
181
182 end do loop_sensor
183
184 !---------------------------------------------------------------------
185 ! 3.0 get CRTM sensor descriptor to use in CRTM_Set_ChannelInfo later
186 !---------------------------------------------------------------------
187 allocate(Sensor_Descriptor(nsensor))
188 call da_crtm_sensor_descriptor(nsensor,Sensor_Descriptor)
189 allocate(ChannelInfo(nsensor))
190
191 ! CRTM load coefficients
192 !-----------------------------------------------------------
193 ! 3.1 call CRTM_Init to load coefficients and fill ChannelInfo structure
194 !-----------------------------------------------------------
195 ! INPUT:
196 !SpcCoeff_File = 'amsua_n15-n16.SpcCoeff.bin'
197 !TauCoeff_File = 'amsua_n15-n16.TauCoeff.bin'
198 AerosolCoeff_File = 'AerosolCoeff.bin'
199 CloudCoeff_File = 'CloudCoeff.bin'
200 EmisCoeff_File = 'EmisCoeff.bin'
201 File_Path = 'crtm_coeffs/' ! path to coefficient files (if not specified current directory is used)
202 !----------------------------------------------------------------
203 ! ChannelInfo structure contains on output:
204 !
205 ! n_channels - integer, total number of channels
206 ! Sensor_Index - integer
207 ! Channel_Index - integer pointer, index of the channels loaded during initialization
208 ! Sensor_Channel - integer pointer, the sensor channel #
209 ! SensorID - character pointer, character string containing satellite and sensor descr
210 ! example: amsre_aqua (Appendix B in User guide)
211 ! WMO_Satellite_ID - integer pointer
212 ! WMO_Sensor_ID - integer pointer
213 !-----------------------------------------------------------------
214
215 Error_Status = CRTM_Init( ChannelInfo, &
216 SensorID = Sensor_Descriptor, &
217 AerosolCoeff_File = AerosolCoeff_File, &
218 EmisCoeff_File = EmisCoeff_File , &
219 CloudCoeff_File = CloudCoeff_File, &
220 File_Path = File_Path)
221
222 ! Error_Status = CRTM_Set_ChannelInfo(Sensor_Descriptor(1),ChannelInfo)
223
224 IF ( Error_Status /= 0 ) THEN
225 call da_error(__FILE__,__LINE__, &
226 (/"Error in initializing CRTM"/))
227 END IF
228
229 ! do n = 1, nsensor
230 ! write (6,*) 'in da_crtm_init: ChannelInfo content'
231 ! write (6,*) 'Sensor_Index ',ChannelInfo(n)%Sensor_Index
232 ! write (6,*) 'n_channels ',ChannelInfo(n)%n_channels
233 ! write (6,*) 'Channel_Index ',ChannelInfo(n)%Channel_Index(:)
234 ! write (6,*) 'Sensor_Channel ',ChannelInfo(n)%Sensor_Channel(:)
235 ! write (6,*) 'SensorID ',ChannelInfo(n)%SensorID
236 ! end do
237
238 !-------------------------------------------------------
239 ! 4.0 read bias correction coefs files
240 !-------------------------------------------------------
241
242 if (read_biascoef) then
243 loop_sensor2: do n = 1, nsensor
244
245 allocate ( satinfo(n) % scanbias (nchanl(n),nscan(n)) )
246 allocate ( satinfo(n) % scanbias_b(nchanl(n),nscan(n),18) )
247 allocate ( satinfo(n) % bcoef (nchanl(n),4) )
248 allocate ( satinfo(n) % bcoef0 (nchanl(n)) )
249 allocate ( satinfo(n) % error_std (nchanl(n)) )
250
251 satinfo(n) % error_std(:) = 500.
252 satinfo(n) % scanbias(:,:) = 0.
253 satinfo(n) % scanbias_b(:,:,:) = 0.
254 satinfo(n) % bcoef(:,:) = 0.
255 satinfo(n) % bcoef0(:) = 0.
256
257 if ( index(iv%instid(n)%rttovid_string,'eos') > 0 ) cycle ! not implemented
258 if ( index(iv%instid(n)%rttovid_string,'hirs') > 0 ) cycle ! not implemented
259 if ( index(iv%instid(n)%rttovid_string,'ssmis') > 0) cycle ! not implemented
260
261 ! new bias coefs files
262 !----------------------------------
263 call da_read_biascoef(iv%instid(n)%rttovid_string, &
264 nchanl(n),nscan(n),18,4,global, &
265 satinfo(n)%scanbias, &
266 satinfo(n)%scanbias_b, &
267 satinfo(n)%bcoef, &
268 satinfo(n)%bcoef0, &
269 satinfo(n)%error_std)
270
271 end do loop_sensor2
272 end if
273
274 !-------------------------------------------------------
275 ! 5.0 read error factor file
276 !-------------------------------------------------------
277 if (use_error_factor_rad) then
278
279 do n=1,rtminit_nsensor
280 allocate ( satinfo(n)%error_factor(1:nchanl(n)) )
281 satinfo(n)%error_factor(:) = 1.0
282 end do
283
284 call da_get_unit(unit_factor_rad)
285 open(unit_factor_rad, file='wrfvar_run/radiance_error.factor', &
286 form='formatted',iostat = iost, status='old')
287
288 if (iost /= 0) then
289 call da_error(__FILE__,__LINE__, &
290 (/"Cannot open radiance error factor file: radiance_error.factor"/))
291 end if
292
293 read(unit_factor_rad, *)
294 do
295 read(unit_factor_rad,fmt='(a15,i3,i8,3f15.5,f8.3)',iostat=iost) &
296 rttovid_string, ichan, num_tot, joa,jo,trace,factor
297 if ( iost == 0 ) then
298 do n=1,rtminit_nsensor
299 if ( index(rttovid_string,trim(iv%instid(n)%rttovid_string))>0 ) then
300 satinfo(n)%error_factor(ichan) = factor
301 write(6,'(a,i5,a,f10.3)') trim(rttovid_string)//' Channel ', ichan, ' Error Factor = ', factor
302 exit
303 end if
304 end do
305 else
306 exit
307 end if
308 end do
309 close(unit_factor_rad)
310 call da_free_unit(unit_factor_rad)
311
312 end if
313
314 !---------------------------------------------
315 ! 6.0 get FGAT time slots
316 !---------------------------------------------
317 allocate ( time_slots(0:num_fgat_time) )
318 call da_get_time_slots(num_fgat_time,time_window_min,analysis_date, &
319 time_window_max, time_slots)
320
321
322 deallocate(nscan)
323 deallocate(nchanl)
324
325 call da_trace_exit("da_crtm_init")
326 #endif
327 end subroutine da_crtm_init