| blob | e627396308c3171a1f00b37037e45e794dbaa735 |
1 subroutine da_setup_obs_structures( grid, ob, iv, j_cost)
3 !---------------------------------------------------------------------------
4 ! Purpose: Allocate and read in components of observation structure.
5 !---------------------------------------------------------------------------
7 implicit none
9 type (domain), intent(inout) :: grid ! Model data
10 type (y_type), intent(out) :: ob ! Observation structure.
11 type (iv_type), intent(out) :: iv ! O-B structure.
12 type (j_type), intent(inout) :: j_cost ! Cost function.
14 integer :: i,j,icode
15 character (len=14) :: ob_name
17 ! thinning variables
18 integer :: istart,iend,jstart,jend
19 real :: rlonlat(4)
21 if (trace_use) call da_trace_entry("da_setup_obs_structures")
23 ! Default values for fmt_info, fmt_srfc, and fmt_each(YRG, 02/25/2009):
25 fmt_info = '(a12,1x,a19,1x,a40,1x,i6,3(f12.3,11x),6x,a40)'
26 fmt_srfc = '(7(:,f12.3,i4,f7.2))'
27 fmt_each = '(3(f12.3,i4,f7.2),11x,3(f12.3,i4,f7.2),11x,1(f12.3,i4,f7.2))'
29 call da_message((/'Set up observations (ob)'/))
31 ! Adjust obs switches
32 if (use_pseudo_rad .or. num_pseudo > 0) then
33 call da_message((/"Single OBS Test:: Turn off all the OBS switches ***"/))
34 use_synopobs = .false.
35 use_shipsobs = .false.
36 use_metarobs = .false.
37 use_soundobs = .false.
38 use_mtgirsobs = .false.
39 use_tamdarobs = .false.
40 use_bogusobs = .false.
41 use_pilotobs = .false.
42 use_airepobs = .false.
43 use_geoamvobs = .false.
44 use_polaramvobs = .false.
45 use_buoyobs = .false.
46 use_profilerobs = .false.
47 use_satemobs = .false.
48 use_gpspwobs = .false.
49 use_gpsztdobs = .false.
50 use_gpsrefobs = .false.
51 use_gpsephobs = .false.
52 use_ssmiretrievalobs = .false.
53 use_ssmitbobs = .false.
54 use_ssmt1obs = .false.
55 use_ssmt2obs = .false.
56 use_qscatobs = .false.
57 use_hirs2obs = .false.
58 use_hirs3obs = .false.
59 use_hirs4obs = .false.
60 use_mhsobs = .false.
61 use_msuobs = .false.
62 use_amsuaobs = .false.
63 use_amsubobs = .false.
64 use_mwtsobs = .false.
65 use_mwhsobs = .false.
66 use_atmsobs = .false.
67 use_airsobs = .false.
68 use_eos_amsuaobs = .false.
69 use_hsbobs = .false.
70 use_ahiobs = .false.
71 use_mwhs2obs = .false.
72 use_gmiobs = .false.
73 use_goesabiobs = .false.
74 use_obsgts = .false.
75 use_rad = .false.
76 use_airsretobs = .false.
77 use_rainobs = .false.
78 use_iasiobs = .false.
79 use_radarobs = .false.
80 use_radar_rv = .false.
81 use_radar_rf = .false.
82 use_lightningobs = .false.
83 use_lightning_w = .false.
84 use_lightning_div = .false.
85 use_lightning_qv = .false.
86 #if (WRF_CHEM == 1)
87 use_chemic_surfobs = .false.
88 #endif
89 ob_format = ob_format_ascii
90 end if
92 if (use_synopobs .OR. use_shipsobs .OR. use_metarobs .OR. use_pilotobs .OR. &
93 use_profilerobs .OR. use_buoyobs .OR. use_soundobs .OR. use_mtgirsobs .OR. use_bogusobs .OR. &
94 use_satemobs .OR. use_geoamvobs .OR. use_polaramvobs .OR. use_airepobs .OR. use_tamdarobs .OR. &
95 use_gpspwobs .OR. use_gpsztdobs .OR. use_gpsrefobs .OR. use_ssmiretrievalobs .OR. &
96 use_ssmitbobs .OR. use_ssmt1obs .OR. use_ssmt2obs .OR. use_qscatobs .OR. &
97 use_airsretobs .OR. use_lsac .OR. use_gpsephobs) then
99 use_obsgts = .true.
100 else
101 use_obsgts = .false.
102 end if
104 if (use_hirs2obs .OR. use_hirs3obs .OR. use_msuobs .OR. use_amsuaobs .OR. &
105 use_amsubobs .OR. use_airsobs .OR. use_eos_amsuaobs .OR. &
106 use_hsbobs .OR. use_kma1dvar .OR. use_filtered_rad .OR. &
107 use_ssmisobs .OR. use_hirs4obs .OR. use_mhsobs .OR. use_pseudo_rad .OR. &
108 use_mwtsobs .OR. use_mwhsobs .OR. use_atmsobs .OR. use_simulated_rad .OR. &
109 use_iasiobs .OR. use_seviriobs .OR. use_amsr2obs .OR. use_goesimgobs .OR. &
110 use_ahiobs .OR. use_mwhs2obs .OR. use_gmiobs .OR. use_goesabiobs) then
111 use_rad = .true.
112 else
113 use_rad = .false.
114 end if
116 ! test_dm_exact can be set to .true. to force DM_PARALLEL runs
117 ! to produce results that are bitwise-identical regardless of the number of
118 ! MPI tasks used. This is useful for validating that multi-processor runs
119 ! are not a source of bugs. Runtime will be much longer. This option is
120 ! automatically overridden to .false. for serial or 1-MPI-task runs.
122 if (test_dm_exact) then
123 if (num_procs == 1) then
124 test_dm_exact = .false.
125 write(unit=stdout,fmt='(A)') &
126 ' test_dm_exact overridden to .false. for serial or 1-MPI-task run'
127 end if
128 ! only implmenented for Sound and Synop, so switch other types off
129 use_shipsobs = .false.
130 use_metarobs = .false.
131 use_bogusobs = .false.
132 use_pilotobs = .false.
133 use_airepobs = .false.
134 use_geoamvobs = .false.
135 use_polaramvobs = .false.
136 use_buoyobs = .false.
137 use_profilerobs = .false.
138 use_satemobs = .false.
139 use_gpspwobs = .false.
140 use_gpsztdobs = .false.
141 use_gpsrefobs = .false.
142 use_gpsephobs = .false.
143 use_ssmiretrievalobs = .false.
144 use_ssmitbobs = .false.
145 use_ssmt1obs = .false.
146 use_ssmt2obs = .false.
147 use_qscatobs = .false.
148 use_hirs2obs = .false.
149 use_hirs3obs = .false.
150 use_hirs4obs = .false.
151 use_mhsobs = .false.
152 use_msuobs = .false.
153 use_amsuaobs = .false.
154 use_amsubobs = .false.
155 use_mwtsobs = .false.
156 use_mwhsobs = .false.
157 use_atmsobs = .false.
158 use_airsobs = .false.
159 use_eos_amsuaobs = .false.
160 use_hsbobs = .false.
161 use_ahiobs = .false.
162 use_mwhs2obs = .false.
163 use_gmiobs = .false.
164 use_goesabiobs = .false.
165 use_obsgts = .false.
166 use_rad = .false.
167 use_airsretobs = .false.
168 use_rainobs = .false.
169 use_iasiobs = .false.
170 use_radarobs = .false.
171 use_radar_rv = .false.
172 use_radar_rf = .false.
173 use_lightningobs = .false.
174 use_lightning_w = .false.
175 use_lightning_div = .false.
176 use_lightning_qv = .false.
177 #if (WRF_CHEM == 1)
178 use_chemic_surfobs = .false.
179 #endif
180 end if
182 if ( use_pseudo_rad ) then
183 use_rad = .true.
184 thinning = .false.
185 qc_rad = .false.
186 rtminit_nsensor = 1
187 rtminit_platform = pseudo_rad_platid
188 rtminit_satid = pseudo_rad_satid
189 rtminit_sensor = pseudo_rad_senid
190 end if
192 obs_use(1:num_ob_indexes) = .false.
193 if ( use_airepobs ) obs_use(airep) = .true.
194 if ( use_bogusobs ) obs_use(bogus) = .true.
195 if ( use_buoyobs ) obs_use(buoy) = .true.
196 if ( use_geoamvobs ) obs_use(geoamv) = .true.
197 if ( use_gpsephobs ) obs_use(gpseph) = .true.
198 if ( use_gpsrefobs ) obs_use(gpsref) = .true.
199 if ( use_gpsztdobs .or. use_gpspwobs ) obs_use(gpspw) = .true.
200 if ( use_metarobs ) obs_use(metar) = .true.
201 if ( use_mtgirsobs ) obs_use(mtgirs) = .true.
202 if ( use_pilotobs ) obs_use(pilot) = .true.
203 if ( use_polaramvobs ) obs_use(polaramv) = .true.
204 if ( use_profilerobs ) obs_use(profiler) = .true.
205 if ( use_qscatobs ) obs_use(qscat) = .true.
206 if ( use_radarobs ) obs_use(radar) = .true.
207 if ( use_lightningobs ) obs_use(lightning) = .true.
208 if ( use_rainobs ) obs_use(rain) = .true.
209 if ( use_satemobs ) obs_use(satem) = .true.
210 if ( use_shipsobs ) obs_use(ships) = .true.
211 if ( use_soundobs ) then
212 obs_use(sound) = .true.
213 obs_use(sonde_sfc) = .true.
214 end if
215 if ( use_ssmiretrievalobs ) obs_use(ssmi_rv) = .true.
216 if ( use_ssmitbobs ) obs_use(ssmi_tb) = .true.
217 if ( use_ssmt1obs ) obs_use(ssmt1) = .true.
218 if ( use_ssmt2obs ) obs_use(ssmt2) = .true.
219 if ( use_synopobs ) obs_use(synop) = .true.
220 if ( use_tamdarobs ) then
221 obs_use(tamdar) = .true.
222 obs_use(tamdar_sfc) = .true.
223 end if
225 if (sfc_assi_options < 1 .OR. sfc_assi_options > 2) then
226 write(unit=message(1),fmt='(A,I3)') &
227 'Invalid sfc_assi_option = ', sfc_assi_options
228 call da_error(__FILE__,__LINE__,message(1:1))
229 end if
231 if ( use_gpsrefobs ) then
232 if ( ob_format == ob_format_bufr ) then
233 ! when main conv input is from bufr, gpsro has to be from bufr too
234 ! overwrite the namelist setting
235 ob_format_gpsro = ob_format_bufr
236 else
237 if ( ob_format == ob_format_ascii .and. &
238 ob_format_gpsro == ob_format_bufr ) then
239 call da_message((/'ob_format=2 while ob_format_gpsro=1, will be reading GPSRO data from gpsro.bufr'/))
240 end if
241 end if
242 end if
244 ! convert pseudo_var to be in lowercase
245 do i = 1, len_trim(pseudo_var)
246 icode = ichar(pseudo_var(i:i))
247 if (icode>=65 .and. icode<=90) then
248 pseudo_var(i:i)=char(icode + 97 - 65)
249 end if
250 end do
252 if ( num_pseudo > 0 .and. &
253 (trim(adjustl(pseudo_var))=='qcw' .or. &
254 trim(adjustl(pseudo_var))=='qrn' .or. &
255 trim(adjustl(pseudo_var))=='qci' .or. &
256 trim(adjustl(pseudo_var))=='qsn' .or. &
257 trim(adjustl(pseudo_var))=='qgr') ) then
258 if ( cloud_cv_options < 2 ) then
259 write(unit=message(1),fmt='(A,A)') &
260 'cloud_cv_options should be 2 or 3 for pseudo_var = ', trim(pseudo_var)
261 call da_error(__FILE__,__LINE__,message(1:1))
262 end if
263 end if
265 pseudo_tpw = .false.
266 pseudo_ztd = .false.
267 pseudo_ref = .false.
268 if ( num_pseudo > 0 .and. trim(adjustl(pseudo_var))=='tpw' ) then
269 pseudo_tpw = .true.
270 use_gpspwobs = .true.
271 end if
272 if ( num_pseudo > 0 .and. trim(adjustl(pseudo_var))=='ztd' ) then
273 pseudo_ztd = .true.
274 use_gpsztdobs = .true.
275 end if
276 if ( num_pseudo > 0 .and. trim(adjustl(pseudo_var))=='ref' ) then
277 pseudo_ref = .true.
278 use_gpsrefobs = .true.
279 end if
280 ! if pseudo ob for conventional u, v, t, p, q
281 if ( num_pseudo > 0 ) then
282 pseudo_uvtpq = .true.
283 if ( pseudo_tpw .or. pseudo_ztd .or. pseudo_ref ) then
284 pseudo_uvtpq = .false.
285 end if
286 end if
287 if ( num_pseudo > 0 ) then
288 !to avoid errors when writing out filtered_obs for undefined variables
289 anal_type_qcobs = .false.
290 end if
292 if ( thin_conv .or. thin_conv_ascii ) then
293 thin_conv_opt(gpseph) = no_thin
294 thin_conv_opt(radar) = no_thin
295 thin_conv_opt(radiance) = no_thin
296 thin_conv_opt(rain) = no_thin
297 thin_conv_opt(lightning) = no_thin
298 if ( thin_conv .and. ob_format==ob_format_bufr ) then
299 ! gpsref horizontal thinning is not implemented for bufr input
300 thin_conv_opt(gpsref) = no_thin
301 end if
302 do i = 1, num_ob_indexes
303 ! resetting thin_conv_opt to no_thin for ob types that are not used
304 if ( thin_conv_opt(i) > no_thin .and. &
305 (.not. obs_use(i)) ) then
306 thin_conv_opt(i) = no_thin
307 end if
308 end do
309 if ( thin_conv .and. ob_format==ob_format_bufr .and. &
310 all(thin_conv_opt(:) <= no_thin) ) thin_conv = .false.
311 if ( thin_conv_ascii .and. ob_format==ob_format_ascii .and. &
312 all(thin_conv_opt(:) <= no_thin) ) thin_conv_ascii = .false.
313 end if
314 if ( ((.not. thin_conv) .and. ob_format==ob_format_bufr) .or. &
315 ((.not. thin_conv_ascii) .and. ob_format==ob_format_ascii) ) then
316 thin_conv_opt(:) = no_thin
317 end if
318 write(message(1),'(a,30i2)') 'thin_conv_opt = ', thin_conv_opt(:)
319 call da_message(message(1:1))
321 ! get lat/lon info (rlat_min,rlat_max,rlon_min,rlon_max,dlat_grid,dlon_grid)
322 ! to be used for making thinning grids for conv (ascii/bufr), radiance and rain ob types.
323 call init_constants_derived
324 rlat_min = r999
325 rlat_max = -r999
326 rlon_min = r999
327 rlon_max = -r999
328 istart=MINVAL( grid%i_start(1:grid%num_tiles) )
329 iend =MAXVAL( grid%i_end (1:grid%num_tiles) )
330 jstart=MINVAL( grid%j_start(1:grid%num_tiles) )
331 jend =MAXVAL( grid%j_end (1:grid%num_tiles) )
332 do i = istart, iend
333 do j = jstart, jend
334 rlat_min=min(rlat_min, grid%xb%lat(i,j))
335 rlat_max=max(rlat_max, grid%xb%lat(i,j))
336 if( grid%xb%lon(i,j) < 0.0) then
337 rlon_min=min(rlon_min, (r360+grid%xb%lon(i,j)))
338 rlon_max=max(rlon_max, (r360+grid%xb%lon(i,j)))
339 else
340 rlon_min=min(rlon_min, grid%xb%lon(i,j))
341 rlon_max=max(rlon_max, grid%xb%lon(i,j))
342 endif
343 enddo
344 enddo
345 #ifdef DM_PARALLEL
346 call mpi_reduce(rlat_min, rlonlat(1), 1, true_mpi_real, mpi_min, root, comm, ierr)
347 call mpi_reduce(rlon_min, rlonlat(2), 1, true_mpi_real, mpi_min, root, comm, ierr)
348 call mpi_reduce(rlat_max, rlonlat(3), 1, true_mpi_real, mpi_max, root, comm, ierr)
349 call mpi_reduce(rlon_max, rlonlat(4), 1, true_mpi_real, mpi_max, root, comm, ierr)
350 CALL mpi_bcast( rlonlat, 4 , true_mpi_real , root , comm, ierr )
351 rlat_min = rlonlat(1)
352 rlon_min = rlonlat(2)
353 rlat_max = rlonlat(3)
354 rlon_max = rlonlat(4)
355 #endif
356 dlat_grid = rlat_max - rlat_min
357 dlon_grid = rlon_max - rlon_min
359 !---------------------------------------------------------------------------
360 ! [1.0] Setup and read in fields from first guess:
361 !----------------------------------------------------------------------------
363 iv%missing = missing
364 ! iv%ptop = grid%xb%ptop
366 iv%total_rad_pixel = 0
367 iv%total_rad_channel = 0
369 iv%info(:)%nlocal = 0
370 iv%info(:)%ntotal = 0
371 iv%info(:)%thin_nlocal = 0
372 iv%info(:)%thin_ntotal = 0
373 do i=1,num_ob_indexes
374 iv%info(i)%plocal(:) = 0
375 iv%info(i)%ptotal(:) = 0
376 iv%info(i)%thin_plocal(:) = 0
377 iv%info(i)%thin_ptotal(:) = 0
378 end do
379 iv%num_inst = 0
381 ob%nlocal(:) = 0
382 ob%ntotal(:) = 0
383 ob%num_inst = 0
385 iv%info(:)%max_lev = 1
387 ! get FGAT time slots
389 allocate ( time_slots(0:num_fgat_time) )
390 call da_get_time_slots(num_fgat_time,time_window_min,analysis_date, &
391 time_window_max, time_slots, ifgat_ana)
393 if (use_obsgts) then
394 ! Conventional obs can be in BUFR or ascii format
395 if (ob_format == ob_format_bufr) then
396 call da_message((/'Using PREPBUFR format observation input'/))
397 call da_setup_obs_structures_bufr (grid, ob, iv)
398 else if (ob_format == ob_format_ascii) then
399 call da_message((/'Using ASCII format observation input'/))
400 call da_setup_obs_structures_ascii (ob, iv, grid)
401 else
402 write(unit=message(1),fmt='(A,I3)') 'Invalid ob_format = ', ob_format
403 call da_error(__FILE__,__LINE__,message(1:1))
404 end if
405 end if
407 if (use_radarobs) then
408 ! Radar obs are read from separate file(s)
409 call da_message((/'Using ASCII format radar observation input'/))
410 call da_setup_obs_structures_radar (grid, ob, iv)
411 end if
413 if (use_lightningobs) then
414 ! Lightning obs are read from separate file(s)
415 call da_message((/'Using ASCII format lightning observation input'/))
416 call da_setup_obs_structures_lightning (grid, ob, iv)
417 end if
419 if (use_rainobs .and. var4d) then
420 call da_message((/'Using ASCII format precipitation observation input'/))
421 call da_setup_obs_structures_rain (grid, ob, iv)
422 end if
424 #if (WRF_CHEM == 1)
425 if ( (use_chemic_surfobs) ) then
426 call da_message((/'Using ASCII format chem IC observation input'/))
427 call da_setup_obs_structures_chem_sfc (ob, iv, grid)
428 end if
429 #endif
431 if (use_rad) then
432 #if defined(CRTM) || defined(RTTOV)
433 ! Radiance files can only be in BUFR
434 call da_message((/'Using NCEP BUFR radiance 1b input'/))
435 if ( use_amsr2obs ) then
436 call da_message((/'Using AMSR2 radiance input in HDF5 format'/))
437 end if
438 if ( use_goesimgobs ) then
439 call da_message((/'Using GOES IMAGER radiance input in netcdf format'/))
440 end if
441 if ( use_goesabiobs ) then
442 call da_message((/'Using GOES ABI radiance input in netcdf format'/))
443 end if
444 if ( use_ahiobs ) then
445 call da_message((/'Using himawari AHI radiance input in netcdf format'/))
446 end if
447 if ( use_gmiobs ) then
448 call da_message((/'Using GMI radiance input in HDF5 format'/))
449 end if
450 call da_setup_radiance_structures(grid, ob, iv)
451 #else
452 call da_error(__FILE__,__LINE__,(/"options to use radiance data are turned on in wrfvar4 namelist, the code needs to be compiled with CRTM or RTTOV library"/))
453 #endif
454 end if
456 if ( num_pseudo > 0 ) then
457 call da_setup_pseudo_obs(grid, iv, ob)
458 end if
460 ! Summarize observations
462 write(unit=stdout, fmt='(a)') 'Observation summary'
463 do i=1,num_fgat_time
464 write(unit=stdout, fmt='(3x,2(a,i3))') 'ob time', i,' num_ob_indexes=',num_ob_indexes
465 do j=1,num_ob_indexes
466 if(j.eq.27) cycle ! skip tamdar_sfc (HA)
467 if (use_gpsztdobs .and.obs_names(j) == 'gpspw ' ) then
468 ob_name = 'gpsztd '
469 else
470 ob_name = obs_names(j)
471 endif
472 if(iv%info(j)%ptotal(i) - iv%info(j)%ptotal(i-1) > 0) then
473 if ( thin_conv_ascii .and. ob_format == ob_format_ascii ) then
474 if ( j == radiance ) cycle !temporary fix to not print incorrect counts for radiance
475 write(unit=stdout, fmt='(2x,2(a,i10),a,i8,a)') &
476 ob_name, iv%info(j)%ptotal(i) - iv%info(j)%ptotal(i-1), ' global, ', &
477 iv%info(j)%thin_ptotal(i) - iv%info(j)%thin_ptotal(i-1), ' global (post-thinning), ', &
478 iv%info(j)%thin_plocal(i) - iv%info(j)%thin_plocal(i-1), ' local (post-thinning)'
479 else
480 #if (WRF_CHEM == 1)
481 write(unit=stdout, fmt='(i3,3x,a,i10,a,i8,a,i8,a,i8,a)') &
482 j,ob_name, iv%info(j)%ptotal(i) - iv%info(j)%ptotal(i-1), ' global,', &
483 iv%info(j)%plocal(i) - iv%info(j)%plocal(i-1), ' local', &
484 iv%info(chemic_surf)%ntotal, ' chem total', iv%info(chemic_surf)%nlocal, ' chem local'
485 #else
486 write(unit=stdout, fmt='(6x,a,i10,a,i8,a)') &
487 ob_name, iv%info(j)%ptotal(i) - iv%info(j)%ptotal(i-1), ' global,', &
488 iv%info(j)%plocal(i) - iv%info(j)%plocal(i-1), ' local'
489 #endif
490 end if ! thin_conv_ascii
491 end if
492 end do
493 end do
494 write(unit=stdout, fmt='(a)') ' '
496 ! Get horizontal interpolation weights.
498 call da_setup_obs_interp_wts (iv)
500 if (trace_use) call da_trace_exit("da_setup_obs_structures")
502 end subroutine da_setup_obs_structures