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