| blob | 3cc5fd05d16d963942e456fd4dcdeaf1dd286fc0 |
2 MODULE module_decoded
4 !--------------------------------------------------------------------------
5 ! Read decoded observations (output format of the MM5 "gts_decoder" and "fetch"
6 ! facilities, which is also the input format
7 ! of MM5 "rawins" and "little_r" programs).
8 !
9 ! Keep only records with at least height or pressure in MM5 horizontal domain
10 ! as defined in namelist
11 !
12 ! Fill the observation data structure
13 !--------------------------------------------------------------------------
14 !
15 ! HISTORY:
16 !
17 ! D. GILL, April 1998
18 ! F. VANDENBERGHE, March 2001
19 !
20 ! 01/13/2003 - Updated for Profiler obs. S. R. H. Rizvi
21 !
22 ! 02/04/2003 - Updated for Buoy obs. S. R. H. Rizvi
23 !
24 ! 02/11/2003 - Reviewed and modified for Profiler
25 ! and Buoy obs. Y.-R. Guo
26 ! 03/30/2005 - Updated for MODIS obs. Syed RH Rizvi
27 ! 06/30/2006 - Updated for AIRS retrievals Syed RH Rizvi
28 ! 10/09/2006 - Updated for GPS Excess Phase Y.-R. Guo
29 !-----------------------------------------------------------------------------
30 ! DATA STRUCTURES
31 !--------------------------------------------------------------------------
33 USE module_type
34 use module_stntbl
36 !------------------------------------------------------------------------
37 ! FUNCTION
38 !------------------------------------------------------------------------
40 USE module_func
42 !------------------------------------------------------------------------
43 ! EXTERNAL
44 !------------------------------------------------------------------------
46 USE module_inside
47 use module_namelist, only: gts_from_mmm_archive, calc_psfc_from_QNH
49 !--------------------------------------------------------------------------
50 ! PARAMETERS
51 !---------------------------------------------------------------------------
53 INCLUDE 'missing.inc'
55 ! Define error return codes used by 'read_measurements' routine.
57 INTEGER , PARAMETER :: ok = 0 , &
58 eof_err = 1 , &
59 no_data = 2 , &
60 read_err = 3 , &
62 ! ssmi_qc_index is the index for good quality of SSMI Tb data from
63 ! the SSMI decoder program. (Y.-R. Guo, 09/03/2003)
65 ssmi_qc_index = 0
66 !
67 ! FORMAT STRINGS for input/output of data.
68 ! These format strings correspond to the data structures in this file
69 ! and are used for input and output of values in the 'report' structure
70 ! (first format string) and the 'measurement' structure (second format).
71 ! Note that report struct contains the first of a linked list of
72 ! measurements; this first meas is read using the 'meas_format'.
74 CHARACTER ( LEN = 120 ) , PARAMETER :: rpt_format = &
75 ' (2F20.5 , 2A40 , ' & ! format for location_type
76 // ' 2A40 , 1F20.5 , 5I10 , 3L10 , ' & ! format for source_info
77 // ' 2I10 , A20 , ' & ! fmt for valid_time
78 // ' 13(F13.5 , I7),' & ! fmt for 'terrestrial'
79 // '1(:,F13.5 , I7),' & ! fmt for PW other than GPS
80 // '7(:,F13.5 , I7))' ! fmt for Brightness Temp
82 CHARACTER ( LEN = 120 ) , PARAMETER :: meas_format = &
83 ' ( 10( F13.5 , I7 ) ) ' ! fmt for measurement rcd
85 CHARACTER ( LEN = 120 ) , PARAMETER :: end_format = &
86 ' ( 3 ( I7 ) ) ' ! fmt for end record
88 INTEGER :: N_air = 0, N_air_cut = 0
90 ! -------------------------------------------------------------------------
91 ! ROUTINES
92 ! -------------------------------------------------------------------------
94 CONTAINS
96 ! SUBROUTINE read_obs_gts ( file_name , file_num , obs , n_obs , &
97 ! SUBROUTINE read_measurements (file_num, surface, location, bad_data, &
98 ! SUBROUTINE dealloc_meas ( head )
99 ! SUBROUTINE surf_first ( surface , elevation )
100 ! SUBROUTINE diagnostics ( new , longitude )
101 ! SUBROUTINE insert_at ( surface , new , elevation )
103 !
104 !-------------------------------------------------------------------------
107 SUBROUTINE read_obs_gts ( file_name, obs , n_obs , &
108 total_number_of_obs , fatal_if_exceed_max_obs , print_gts_read , &
109 ins , jew , &
110 time_window_min, time_window_max, map_projection , missing_flag)
112 ! This routine opens file 'file_name' and reads all observations, and
113 ! measurements at all levels, from the file into the 'obs' observation array.
114 ! For each observation, calls read_measurements to read all measurements for
115 ! that one observation.
117 USE module_date
118 USE module_per_type
120 IMPLICIT NONE
122 CHARACTER ( LEN = * ) , INTENT ( IN ) :: file_name
123 INTEGER, INTENT (INOUT) :: n_obs
124 INTEGER, INTENT ( IN ) :: total_number_of_obs
125 LOGICAL, INTENT ( IN ) :: fatal_if_exceed_max_obs
126 LOGICAL, INTENT ( IN ) :: print_gts_read
127 TYPE (report), DIMENSION (:), INTENT (OUT) :: obs
128 INTEGER, INTENT ( IN ) :: ins, jew
130 CHARACTER (LEN = 19) , INTENT (IN) :: time_window_min
131 CHARACTER (LEN = 19) , INTENT (IN) :: time_window_max
132 INTEGER :: map_projection
133 REAL, INTENT (OUT) :: missing_flag
135 INTEGER :: file_num
136 CHARACTER ( LEN = 32 ) , PARAMETER :: proc_name = 'read_obs_gts '
137 INTEGER :: io_error, platform_error
138 INTEGER :: obs_num
139 INTEGER :: error_ret
141 INTEGER :: num_empty , num_outside , bad_count
142 LOGICAL :: outside_domain, outside_window
143 LOGICAL :: outside
145 TYPE(meas_data) :: dummy_middle
147 INTEGER :: icrs, k, iunit, levels
148 CHARACTER ( LEN = 80) :: filename
149 CHARACTER ( LEN = 160) :: error_message
150 CHARACTER ( LEN = 14) :: newstring
151 INTEGER :: nlevels, num_unknown, m_miss
152 TYPE ( measurement ) , POINTER :: current
153 real :: qnh, elev
154 real :: xla, xlo
155 integer :: ipos
156 !-----------------------------------------------------------------------------!
157 INCLUDE 'platform_interface.inc'
158 !-----------------------------------------------------------------------------!
160 ! INCLUDE 'error.inc'
161 ! INTERFACE
162 ! INCLUDE 'error.int'
163 ! END INTERFACE
165 !------------------------------------------------------------------------------!
167 WRITE (UNIT = 0, FMT = '(A)') &
168 '------------------------------------------------------------------------------'
170 WRITE (UNIT = 0, FMT = '(A,A,/)') &
171 'READ GTS OBSERVATIONS IN FILE ', TRIM (file_name)
173 ! empty or outside of the domain.
174 ! Initialize a couple of counters for how many observations are either
175 ! empty or outside of the domain.
177 num_unknown = 0
178 missing_flag = missing_r
179 num_empty = 0
180 num_outside = 0
181 m_miss = 0
183 ! Open file for writing diagnostics
184 IF (print_gts_read) THEN
186 filename = 'obs_gts_read.diag'
187 iunit = 999
189 OPEN ( UNIT = iunit , FILE = filename , FORM = 'FORMATTED' , &
190 ACTION = 'WRITE' , STATUS = 'REPLACE', IOSTAT = io_error )
192 IF ( io_error .NE. 0 ) THEN
193 CALL error_handler (proc_name, &
194 'Unable to open output diagnostic file:', filename, .true.)
195 ELSE
196 WRITE (UNIT = 0, FMT = '(A,A,/)') &
197 "Diagnostics in file ", TRIM (filename)
198 ENDIF
200 ENDIF
202 ! Open file for reading; handle any errors in open by quitting since
203 ! this is probably a simple-to-fix user mistake.
205 file_num = 99
206 OPEN ( UNIT = file_num , FILE = file_name , FORM = 'FORMATTED' , &
207 ACTION = 'READ' , IOSTAT = io_error )
209 IF ( io_error .NE. 0 ) THEN
210 CALL error_handler (proc_name, &
211 'Unable to open gts input observations file: ',file_name, .true.)
212 ENDIF
214 !
215 ! allow changing the obs date for testing
216 !
217 CALL GETENV('OBSPROC_TEST_DATE',newstring)
219 ! While we are not at the end of the observation file, keep reading data.
221 ! obs_num = 1
222 obs_num = n_obs + 1
224 read_obs : DO while ( io_error == 0 )
225 ! This is an array that we are filling. Are we beyond that limit yet?
227 IF (obs_num .GT. total_number_of_obs) THEN
229 error_message(1:60) = &
230 'Too many obs for the NAMELIST value of max_number_of_obs = '
232 WRITE (error_message(61:67),'(I7)') total_number_of_obs
234 CALL error_handler (proc_name, error_message (1:60), &
235 error_message(61:67),fatal_if_exceed_max_obs)
237 ! If fatal, code will stop above, otherwise the following lines exit
238 ! do loop and close file read
240 CLOSE ( file_num )
241 IF (print_gts_read) CLOSE ( iunit )
242 EXIT read_obs
244 END IF
246 ! initialize the vaiable that will be used later for checking
247 ! if the pw info is read or not
248 obs(obs_num)%ground%pw%data = missing_r
250 ! The first read is the "once only" type of information.
252 READ ( file_num , IOSTAT = io_error , FMT = rpt_format ) &
253 obs(obs_num)%location % latitude, obs(obs_num)%location % longitude, &
254 obs(obs_num)%location % id, obs(obs_num)%location % name, &
255 obs(obs_num)%info % platform, obs(obs_num)%info % source, &
256 obs(obs_num)%info % elevation, obs(obs_num)%info % num_vld_fld, &
257 obs(obs_num)%info % num_error, obs(obs_num)%info % num_warning, &
258 obs(obs_num)%info % seq_num, obs(obs_num)%info % num_dups, &
259 obs(obs_num)%info % is_sound, obs(obs_num)%info % bogus, &
260 obs(obs_num)%info % discard, &
261 obs(obs_num)%valid_time % sut, obs(obs_num)%valid_time % julian, &
262 obs(obs_num)%valid_time % date_char, &
263 obs(obs_num)%ground%slp%data, obs(obs_num)%ground%slp%qc, &
264 obs(obs_num)%ground%ref_pres%data, obs(obs_num)%ground%ref_pres%qc, &
265 obs(obs_num)%ground%ground_t%data, obs(obs_num)%ground%ground_t%qc, &
266 obs(obs_num)%ground%sst%data, obs(obs_num)%ground%sst%qc, &
267 obs(obs_num)%ground%psfc%data, obs(obs_num)%ground%psfc%qc, &
268 obs(obs_num)%ground%precip%data, obs(obs_num)%ground%precip%qc, &
269 obs(obs_num)%ground%t_max%data, obs(obs_num)%ground%t_max%qc, &
270 obs(obs_num)%ground%t_min%data, obs(obs_num)%ground%t_min%qc, &
271 obs(obs_num)%ground%t_min_night%data, obs(obs_num)%ground%t_min_night%qc,&
272 obs(obs_num)%ground%p_tend03%data, obs(obs_num)%ground%p_tend03%qc, &
273 obs(obs_num)%ground%p_tend24%data, obs(obs_num)%ground%p_tend24%qc, &
274 obs(obs_num)%ground%cloud_cvr%data, obs(obs_num)%ground%cloud_cvr%qc, &
275 obs(obs_num)%ground%ceiling%data, obs(obs_num)%ground%ceiling%qc, &
276 obs(obs_num)%ground%pw %data, obs(obs_num)%ground%pw%qc, &
277 obs(obs_num)%ground%tb19v %data, obs(obs_num)%ground%tb19v%qc, &
278 obs(obs_num)%ground%tb19h %data, obs(obs_num)%ground%tb19h%qc, &
279 obs(obs_num)%ground%tb22v %data, obs(obs_num)%ground%tb22v%qc, &
280 obs(obs_num)%ground%tb37v %data, obs(obs_num)%ground%tb37v%qc, &
281 obs(obs_num)%ground%tb37h %data, obs(obs_num)%ground%tb37h%qc, &
282 obs(obs_num)%ground%tb85v %data, obs(obs_num)%ground%tb85v%qc, &
283 obs(obs_num)%ground%tb85h %data, obs(obs_num)%ground%tb85h%qc
285 if ( io_error < 0 ) then ! end-of-file
286 write(unit=0, fmt='(A)') 'Have reached the end of observation file.'
287 close(file_num)
288 if ( print_gts_read ) close(iunit)
289 exit read_obs
290 end if
291 !
292 ! allow changing the obs date for testing (GETENV used above)
293 !
294 if (len_trim(newstring) .ne. 0 ) then
295 obs(obs_num)%valid_time%date_char = newstring
296 endif
298 !------------------------------------------------------------------------------
299 ! If a new data type without a WMO code (like the example of 'AWS SURFACE '
300 ! below) need to be processed, a assumed WMO code (like FM-16) must be
301 ! assigned, otherwise this new data type will be ignored.
302 !
303 ! For AFWA AWS SURFACE in platform, there is no WMO code given.
304 ! Here assume that it is FM-16 similar to metar data:
305 !
306 if (obs(obs_num)%info % platform(1:12) == 'AWS SURFACE ') then
307 obs(obs_num)%info % platform(1:12) = 'FM-16 AWSSFC'
308 endif
309 if (obs(obs_num)%info % platform(1:14) == 'FM-32 PROFILER') then
310 obs(obs_num)%info % platform(1:14) = 'FM-132 PROFILER'
311 endif
312 ! .............................................................................
313 ! Treatment of MODIS winds
314 if(index(obs(obs_num)%location%id, 'MODIS') > 0 .or. &
315 index(obs(obs_num)%location%id, 'modis') > 0 .or. &
316 index(obs(obs_num)%info%source, 'MODIS') > 0 .or. &
317 index(obs(obs_num)%info%source, 'modis') > 0 ) then
319 if( index(obs(obs_num)%info%platform(1:11), 'FM-88') > 0 ) then
320 obs(obs_num)%info%platform(1:11) = 'FM-88 MODIS'
321 !hcl-note: below might be specific to NCAR/MMM dataset
322 obs(obs_num)%location%name = obs(obs_num)%location%id
323 end if
325 end if
326 ! .............................................................................
328 if ( gts_from_mmm_archive ) then
329 ! distinguish BUOY and SHIP (both from BBXX reports assigned FM-13 SHIP)
330 ! by the name. A trick done in NCAR gts_decoder to leave a clue for
331 ! subsequent data processing
332 if ( obs(obs_num)%info%platform(1:10) == 'FM-13 SHIP' ) then
333 if ( index(obs(obs_num)%location%name, 'Platform Id >>>') > 0 ) then
334 obs(obs_num)%info%platform(1:10) = 'FM-18 BUOY'
335 obs(obs_num)%location%id = obs(obs_num)%location%name(17:21)
336 end if
337 end if
338 end if
340 call print_extra_obs
341 ! .............................................................................
342 ! To ignore the data type without WMO code:
344 if (IO_error > 0 ) then ! error reading header info
345 write(0,'("IO_ERROR=",i2,1x,i7,1x,a,2(f9.3,a),1x,a,1x,f11.3)') &
346 io_error, obs_num, obs(obs_num)%info % platform, &
347 obs(obs_num)%location%latitude, 'N',&
348 obs(obs_num)%location%longitude,'E ', &
349 obs(obs_num)%valid_time%date_char, &
350 obs(obs_num)%info % elevation
352 !hcl obs(obs_num)%info % elevation = missing_r
353 obs(obs_num)%info % discard = .true.
354 endif
356 READ (obs(obs_num) % info % platform (4:6), '(I3)', &
357 IOSTAT = platform_error) fm
358 if (platform_error /= 0) then
359 write(0,'(A)') &
360 "***WARNING: NO WMO CODE FOR THIS NEW TYPE OF OBS***"
361 write(0,'(A,I8,A,I4,A,A,A,A)') " ==> obs_num=",obs_num, &
362 " platform_error=",platform_error," platform=", &
363 obs(obs_num) % info % platform (1:12), " ID=", &
364 obs(obs_num) % location % id(1:6)
365 else if (fm <= 0) then
366 write(0,'(A,I8,A,A,A,A)') " ==> obs_num=",obs_num, &
367 " INVALID WMO CODE: platform=", &
368 obs(obs_num) % info % platform (1:12), " ID=", &
369 obs(obs_num) % location % id(1:6)
370 platform_error = -99
371 ! else if (fm == 0) then
372 ! obs(obs_num)%info % platform(1:12) = 'FM-42 AMDAR'
373 endif
375 !-----------------------------------------------------------------------------
377 ! Print read results
379 if ( print_gts_read ) then
380 WRITE (UNIT = iunit , FMT = '(A,1X,A,1X,A,1X,A,1X,2(F8.3,A),A,1X)',&
381 ADVANCE='no') 'Found' , &
382 obs(obs_num)%location%id (1: 5),&
383 obs(obs_num)%location%name (1:20),&
384 obs(obs_num)%info%platform (1: 12),&
385 obs(obs_num)%location%latitude, 'N',&
386 obs(obs_num)%location%longitude,'E ', &
387 obs(obs_num)%valid_time%date_char
388 end if
389 !
390 ! To avoid the floating error in latlon_to_ij calculation: (Guo 01/08/2005)
391 IF (IPROJ > 0) THEN
392 if (truelat1 > 0.0) then
393 if (obs(obs_num)%location%latitude == -90.0) then
394 write(0,'(/i7,1x,"modified the original lat =",f8.2," to -89.5"/)') &
395 obs_num, obs(obs_num)%location%latitude
396 obs(obs_num)%location%latitude = -89.5
397 endif
398 else if (truelat1 < 0.0) then
399 if (obs(obs_num)%location%latitude == 90.0) then
400 write(0,'(/i7,1x,"modified the original lat =",f8.2," to 89.5"/)') &
401 obs_num, obs(obs_num)%location%latitude
402 obs(obs_num)%location%latitude = 89.5
403 endif
404 endif
405 ENDIF
407 ! If PW obs weren't read, arrays are reset to 0, set to missing
409 IF (obs(obs_num)%ground%pw%data .LE. 0.) THEN
410 obs(obs_num)%ground%pw%data = missing_r
411 obs(obs_num)%ground%pw%qc = missing
412 obs(obs_num)%ground%pw%error = missing_r
413 ELSE
414 !------------------------------------------------------------
415 ! GPSPW is in cm and its QC-error in 0.1 mm
416 !
417 ! Note: the variable: pw used for either GPSPW or GPSZTD.
418 ! We NEVER DO the assimilation of both GPSPW or GPSZTD.
419 ! Y.-R. Guo 08/13/2003
420 !
421 ! From GPSPW_decoder, GPSZTD is also in cm and its QC-error
422 ! in the unit of 0.1mm
423 ! Y.-R. Guo 05/09/2008.
424 !
425 !------------------------------------------------------------
427 IF (fm == 111 .or. fm == 114) THEN
428 obs(obs_num)%ground%pw%error = &
429 REAL(obs(obs_num)%ground%pw%qc)/100. ! GPSPW, ZTD
430 ELSE
431 obs(obs_num)%ground%pw%error = missing_r
432 ENDIF
433 obs(obs_num)%ground%pw%qc = 0
434 ENDIF
436 ! PW QC for SSMI (AFWA data only)
438 ! IF (obs(obs_num)%ground%pw%qc .EQ. 5) obs(obs_num)%ground%pw%qc = 0
440 ! If Tb's obs weren't read, arrays are reset to 0, set to missing
441 ! If Tb were read, keep only data with QC = 5 (over water).
443 IF ((obs(obs_num)%ground%tb19v%data .LE. 0.) .OR. &
444 (obs(obs_num)%ground%tb19v%qc .NE. ssmi_qc_index)) THEN
445 obs(obs_num)%ground%tb19v%data = missing_r
446 obs(obs_num)%ground%tb19v%qc = missing
447 obs(obs_num)%ground%tb19v%error = missing_r
448 ELSE
449 obs(obs_num)%ground%tb19v%error = missing_r
450 obs(obs_num)%ground%tb19v%qc = 0
451 ENDIF
453 IF ((obs(obs_num)%ground%tb19h%data .LE. 0.) .OR. &
454 (obs(obs_num)%ground%tb19h%qc .NE. ssmi_qc_index)) THEN
455 obs(obs_num)%ground%tb19h%data = missing_r
456 obs(obs_num)%ground%tb19h%qc = missing
457 obs(obs_num)%ground%tb19h%error = missing_r
458 ELSE
459 obs(obs_num)%ground%tb19h%error = missing_r
460 obs(obs_num)%ground%tb19h%qc = 0
461 ENDIF
463 IF ((obs(obs_num)%ground%tb22v%data .LE. 0.) .OR. &
464 (obs(obs_num)%ground%tb22v%qc .NE. ssmi_qc_index)) THEN
465 obs(obs_num)%ground%tb22v%data = missing_r
466 obs(obs_num)%ground%tb22v%qc = missing
467 obs(obs_num)%ground%tb22v%error = missing_r
468 ELSE
469 obs(obs_num)%ground%tb22v%error = missing_r
470 obs(obs_num)%ground%tb22v%qc = 0
471 ENDIF
473 IF ((obs(obs_num)%ground%tb37v%data .LE. 0.) .OR. &
474 (obs(obs_num)%ground%tb37v%qc .NE. ssmi_qc_index)) THEN
475 obs(obs_num)%ground%tb37v%data = missing_r
476 obs(obs_num)%ground%tb37v%qc = missing
477 obs(obs_num)%ground%tb37v%error = missing_r
478 ELSE
479 obs(obs_num)%ground%tb37v%error = missing_r
480 obs(obs_num)%ground%tb37v%qc = 0
481 ENDIF
483 IF ((obs(obs_num)%ground%tb37h%data .LE. 0.) .OR. &
484 (obs(obs_num)%ground%tb37h%qc .NE. ssmi_qc_index)) THEN
485 obs(obs_num)%ground%tb37h%data = missing_r
486 obs(obs_num)%ground%tb37h%qc = missing
487 obs(obs_num)%ground%tb37h%error = missing_r
488 ELSE
489 obs(obs_num)%ground%tb37h%error = missing_r
490 obs(obs_num)%ground%tb37h%qc = 0
491 ENDIF
493 IF ((obs(obs_num)%ground%tb85v%data .LE. 0.) .OR. &
494 (obs(obs_num)%ground%tb85v%qc .NE. ssmi_qc_index)) THEN
495 obs(obs_num)%ground%tb85v%data = missing_r
496 obs(obs_num)%ground%tb85v%qc = missing
497 obs(obs_num)%ground%tb85v%error = missing_r
498 ELSE
499 obs(obs_num)%ground%tb85v%error = missing_r
500 obs(obs_num)%ground%tb85v%qc = 0
501 ENDIF
503 IF ((obs(obs_num)%ground%tb85h%data .LE. 0.) .OR. &
504 (obs(obs_num)%ground%tb85h%qc .NE. ssmi_qc_index)) THEN
505 obs(obs_num)%ground%tb85h%data = missing_r
506 obs(obs_num)%ground%tb85h%qc = missing
507 obs(obs_num)%ground%tb85h%error = missing_r
508 ELSE
509 obs(obs_num)%ground%tb85h%error = missing_r
510 obs(obs_num)%ground%tb85h%qc = 0
511 ENDIF
513 ! Initialize model coordinates yi and xj
515 obs(obs_num)%location%yic = missing_r
516 obs(obs_num)%location%yid = missing_r
517 obs(obs_num)%location%xjc = missing_r
518 obs(obs_num)%location%xjd = missing_r
520 ! The number of upper levels must be 0 before reading
522 obs (obs_num) % info % levels = 0
524 ! If there are troubles with the "once only" type of data, we keep trying
525 ! until we either come to the end of this report (and cycle) or we come
526 ! to the end of all of the data (exit).
528 IF ( io_error .GT. 0 .or. platform_error /= 0) THEN
530 WRITE (UNIT = 0, FMT = '(A,A,A,A)') 'Troubles with first line ', &
531 TRIM ( obs(obs_num)%location%id ) , &
532 ' ', TRIM ( obs(obs_num)%location%name )
534 ! Keep track of how many loops we are taking so this is not infinite.
536 bad_count = 0
538 DO WHILE ( io_error .GE. 0 )
540 bad_count = bad_count + 1
542 IF (bad_count .LT. 1000 ) THEN
543 ! READ (file_num, IOSTAT = io_error, FMT = meas_format) dummy_middle
544 READ (file_num, IOSTAT = io_error, FMT = meas_format) &
545 dummy_middle % pressure % data, dummy_middle % pressure % qc, &
546 dummy_middle % height % data, dummy_middle % height % qc, &
547 dummy_middle % temperature % data, dummy_middle % temperature % qc, &
548 dummy_middle % dew_point % data, dummy_middle % dew_point % qc, &
549 dummy_middle % speed % data, dummy_middle % speed % qc, &
550 dummy_middle % direction % data, dummy_middle % direction % qc, &
551 dummy_middle % u % data, dummy_middle % u % qc, &
552 dummy_middle % v % data, dummy_middle % v % qc, &
553 dummy_middle % rh % data, dummy_middle % rh % qc, &
554 dummy_middle % thickness % data, dummy_middle % thickness % qc
556 IF (eps_equal (dummy_middle%pressure%data, end_data_r , 1.)) THEN
557 READ (file_num , IOSTAT = io_error , FMT = end_format ) &
558 obs(obs_num)%info%num_vld_fld , &
559 obs(obs_num)%info%num_error , &
560 obs(obs_num)%info%num_warning
561 WRITE (UNIT = 0, FMT = '(A)') 'Starting to READ a new report.'
562 CYCLE read_obs
563 END IF
565 ELSE
567 WRITE (UNIT = 0, FMT = '(A)')&
568 'Too many attempts to read the data correctly. Exiting read loop.'
569 CLOSE ( file_num )
570 IF (print_gts_read) CLOSE ( iunit )
571 EXIT read_obs
573 END IF
575 END DO
577 ! While trying to find a good read, we came to the end of the file.
578 ! It happens to the best of us.
580 IF ( io_error .LT. 0 ) THEN
582 WRITE (UNIT = 0, FMT='(A)') 'Have reached end of observations file.'
583 CLOSE ( file_num )
584 IF (print_gts_read) CLOSE ( iunit )
585 EXIT read_obs
587 END IF
589 ! this scenario should be handled right after the first read
590 !ELSE IF ( io_error .LT. 0 ) THEN
592 ! ! No errors. This is the intended way to find the end of data mark.
594 ! WRITE (UNIT = 0, FMT = '(A)') 'Have reached end of observations file. '
596 ! CLOSE ( file_num )
597 ! IF (print_gts_read) CLOSE ( iunit )
598 ! EXIT read_obs
600 ELSE IF ( io_error .EQ. 0 ) THEN
602 IF ( domain_check_h ) then
604 CALL inside_domain (obs(obs_num)%location%latitude , &
605 obs(obs_num)%location%longitude , &
606 ins , jew , outside_domain , &
607 obs(obs_num)%location%xjc, &
608 obs(obs_num)%location%yic, &
609 obs(obs_num)%location%xjd, &
610 obs(obs_num)%location%yid)
611 else
612 outside_domain = .FALSE.
613 endif
615 ! The previous read ("once only" data) was valid.
616 ! If any of the data is suspicious, the easiest place to clean
617 ! it up is as soon as we read it in, so we do not have to track
618 ! through the array or accidently hit undefined values that are
619 ! not exactly "our" undefined values.
621 ! Sometimes the date and time are ingested in a silly way.
622 ! Set the valid time to a time guaranteed not to be within
623 ! the time window.
625 IF (INDEX (obs(obs_num)%valid_time%date_char , '*' ) .GT. 0 ) THEN
626 obs(obs_num)%valid_time%date_char = '19000101000000'
627 ELSE IF ((obs(obs_num)%valid_time%date_char( 1: 2) .NE. '19' ) .AND. &
628 (obs(obs_num)%valid_time%date_char( 1: 2) .NE. '20' ) ) THEN
629 obs(obs_num)%valid_time%date_char = '19000101000000'
630 ELSE IF ((obs(obs_num)%valid_time%date_char( 5: 6) .LT. '01' ) .OR. &
631 (obs(obs_num)%valid_time%date_char( 5: 6) .GT. '12' )) THEN
632 obs(obs_num)%valid_time%date_char = '19000101000000'
633 ELSE IF ((obs(obs_num)%valid_time%date_char( 7: 8) .LT. '01' ) .OR. &
634 (obs(obs_num)%valid_time%date_char( 7: 8) .GT. '31' )) THEN
635 obs(obs_num)%valid_time%date_char = '19000101000000'
636 ELSE IF ((obs(obs_num)%valid_time%date_char( 9:10) .LT. '00' ) .OR. &
637 (obs(obs_num)%valid_time%date_char( 9:10) .GT. '23' )) THEN
638 obs(obs_num)%valid_time%date_char = '19000101000000'
639 ELSE IF ((obs(obs_num)%valid_time%date_char(11:12) .LT. '00' ) .OR. &
640 (obs(obs_num)%valid_time%date_char(11:12) .GT. '59' )) THEN
641 obs(obs_num)%valid_time%date_char = '19000101000000'
642 ELSE IF ((obs(obs_num)%valid_time%date_char(13:14) .LT. '00' ) .OR. &
643 (obs(obs_num)%valid_time%date_char(13:14) .GT. '59' )) THEN
644 obs(obs_num)%valid_time%date_char = '19000101000000'
645 ELSE IF (((obs(obs_num)%valid_time%date_char( 5: 6) .EQ. '04' ) .OR. &
646 (obs(obs_num)%valid_time%date_char( 5: 6) .EQ. '06' ) .OR. &
647 (obs(obs_num)%valid_time%date_char( 5: 6) .EQ. '09' ) .OR. &
648 (obs(obs_num)%valid_time%date_char( 5: 6) .EQ. '11' )) .AND.&
649 (obs(obs_num)%valid_time%date_char( 7: 8) .GT. '30' ) ) THEN
650 obs(obs_num)%valid_time%date_char = '19000101000000'
651 ELSE IF ((obs(obs_num)%valid_time%date_char( 5: 6) .EQ. '02' ) .AND. &
652 (nfeb_ch ( obs(obs_num)%valid_time%date_char( 1: 4) ) .LT. &
653 obs(obs_num)%valid_time%date_char( 7: 8) ) ) THEN
654 obs(obs_num)%valid_time%date_char = '19000101000000'
655 END IF
657 CALL inside_window (obs(obs_num)%valid_time%date_char, &
658 time_window_min, time_window_max, &
659 outside_window, iunit)
661 outside = outside_domain .OR. outside_window
663 ! Date at MM5 V3 format
665 WRITE (obs(obs_num)%valid_time%date_mm5, &
666 FMT='(A4,"-",A2,"-",A2,"_",A2,":",A2,":",A2)') &
667 obs(obs_num)%valid_time%date_char ( 1: 4), &
668 obs(obs_num)%valid_time%date_char ( 5: 6), &
669 obs(obs_num)%valid_time%date_char ( 7: 8), &
670 obs(obs_num)%valid_time%date_char ( 9:10), &
671 obs(obs_num)%valid_time%date_char (11:12), &
672 obs(obs_num)%valid_time%date_char (13:14)
674 ! These tests are for the ground type data.
675 ! Missing data is OK, but sometimes it comes in as undefined,
676 ! which meant bad data. Set it all to missing.
678 IF ((obs(obs_num)%ground%slp%data .GT. (undefined1_r - 1.)) .OR. &
679 (obs(obs_num)%ground%slp%data .LT. (undefined2_r + 1.))) THEN
680 obs(obs_num)%ground%slp%data = missing_r
681 obs(obs_num)%ground%slp%qc = missing
682 END IF
683 IF ((obs(obs_num)%ground%ref_pres%data .GT. (undefined1_r - 1.)) .OR. &
684 (obs(obs_num)%ground%ref_pres%data .LT. (undefined2_r + 1.))) THEN
685 obs(obs_num)%ground%ref_pres%data = missing_r
686 obs(obs_num)%ground%ref_pres%qc = missing
687 END IF
688 IF ((obs(obs_num)%ground%ground_t%data .GT. (undefined1_r - 1.)) .OR. &
689 (obs(obs_num)%ground%ground_t%data .LT. (undefined2_r + 1.))) THEN
690 obs(obs_num)%ground%ground_t%data = missing_r
691 obs(obs_num)%ground%ground_t%qc = missing
692 END IF
693 IF ((obs(obs_num)%ground%sst%data .GT. (undefined1_r - 1.)) .OR. &
694 (obs(obs_num)%ground%sst%data .LT. (undefined2_r + 1.))) THEN
695 obs(obs_num)%ground%sst%data = missing_r
696 obs(obs_num)%ground%sst%qc = missing
697 END IF
698 IF ((obs(obs_num)%ground%psfc%data .GT. (undefined1_r - 1.)) .OR. &
699 (obs(obs_num)%ground%psfc%data .LT. (undefined2_r + 1.))) THEN
700 obs(obs_num)%ground%psfc%data = missing_r
701 obs(obs_num)%ground%psfc%qc = missing
702 END IF
703 IF ((obs(obs_num)%ground%precip%data .GT. (undefined1_r - 1.)) .OR. &
704 (obs(obs_num)%ground%precip%data .LT. (undefined2_r + 1.))) THEN
705 obs(obs_num)%ground%precip%data = missing_r
706 obs(obs_num)%ground%precip%qc = missing
707 END IF
708 IF ((obs(obs_num)%ground%t_max%data .GT. (undefined1_r - 1.)) .OR. &
709 (obs(obs_num)%ground%t_max%data .LT. (undefined2_r + 1.))) THEN
710 obs(obs_num)%ground%t_max%data = missing_r
711 obs(obs_num)%ground%t_max%qc = missing
712 END IF
713 IF ((obs(obs_num)%ground%t_min%data .GT. (undefined1_r - 1.)) .OR. &
714 (obs(obs_num)%ground%t_min%data .LT. (undefined2_r + 1.))) THEN
715 obs(obs_num)%ground%t_min%data = missing_r
716 obs(obs_num)%ground%t_min%qc = missing
717 END IF
718 IF ((obs(obs_num)%ground%t_min_night%data .GT. (undefined1_r - 1.)) &
719 .OR. (obs(obs_num)%ground%t_min_night%data .LT. (undefined2_r + 1.))) &
720 THEN
721 obs(obs_num)%ground%t_min_night%data = missing_r
722 obs(obs_num)%ground%t_min_night%qc = missing
723 END IF
724 IF ((obs(obs_num)%ground%p_tend03%data .GT. (undefined1_r - 1.)) .OR. &
725 (obs(obs_num)%ground%p_tend03%data .LT. (undefined2_r + 1.))) THEN
726 obs(obs_num)%ground%p_tend03%data = missing_r
727 obs(obs_num)%ground%p_tend03%qc = missing
728 END IF
729 IF ((obs(obs_num)%ground%p_tend24%data .GT. (undefined1_r - 1.)) .OR. &
730 (obs(obs_num)%ground%p_tend24%data .LT. (undefined2_r + 1.))) THEN
731 obs(obs_num)%ground%p_tend24%data = missing_r
732 obs(obs_num)%ground%p_tend24%qc = missing
733 END IF
734 IF ((obs(obs_num)%ground%cloud_cvr%data .GT. (undefined1_r - 1.)) .OR.&
735 (obs(obs_num)%ground%cloud_cvr%data .LT. (undefined2_r + 1.))) THEN
736 obs(obs_num)%ground%cloud_cvr%data = missing_r
737 obs(obs_num)%ground%cloud_cvr%qc = missing
738 END IF
739 IF ((obs(obs_num)%ground%ceiling%data .GT. (undefined1_r - 1.)) .OR. &
740 (obs(obs_num)%ground%ceiling%data .LT. (undefined2_r + 1.))) THEN
741 obs(obs_num)%ground%ceiling%data = missing_r
742 obs(obs_num)%ground%ceiling%qc = missing
743 END IF
745 ! Write obs id in diagnostic file
747 IF (print_gts_read) THEN
748 IF (outside_domain) THEN
749 WRITE (UNIT = iunit , FMT = '(A)' ) '=> OUT DOMAIN'
750 ELSE IF (outside_window) THEN
751 WRITE (UNIT = iunit , FMT = '(A)' ) '=> OUT WINDOW'
752 ELSE
753 WRITE (UNIT = iunit , FMT = '(A)' ) ''
754 ENDIF
755 ENDIF
757 ! Sort observation per type
759 READ (obs(obs_num) % info % platform (4:6), '(I3)')fm
760 CALL fm_decoder (fm, platform, &
761 synop=nsynops (icor), ship =nshipss (icor), &
762 metar=nmetars (icor), pilot=npilots (icor), &
763 sound=nsounds (icor), satem=nsatems (icor), &
764 satob=nsatobs (icor), airep=naireps (icor), &
765 gpspw=ngpspws (icor), gpszd=ngpsztd (icor), &
766 gpsrf=ngpsref (icor), gpsep=ngpseph (icor), &
767 ssmt1=nssmt1s (icor), &
768 ssmt2=nssmt2s (icor), ssmi =nssmis (icor), &
769 tovs =ntovss (icor), other=nothers (icor), &
770 amdar=namdars (icor), qscat=nqscats (icor), &
771 profl=nprofls (icor), buoy =nbuoyss (icor), &
772 bogus=nboguss (icor), airs = nairss(icor),tamdar=ntamdar(icor) )
774 ! Since no I/O errors, read 1 or more measurements.
775 ! Note that obs(obs_num)%surface is pointer to first node in linked
776 ! list, so it is initially not pointing to anything.
778 NULLIFY (obs(obs_num)%surface)
780 CALL read_measurements( file_num , obs(obs_num)%surface , &
781 obs(obs_num)%location , obs(obs_num)%info, outside , error_ret ,&
782 ins , jew , map_projection, &
783 obs (obs_num) % info % elevation, obs (obs_num) % info % levels, &
784 iunit, print_gts_read)
786 ! An error in the measurements read is handled in a couple of ways.
787 ! A flat out error in the read requires the process to start again
788 ! (cycle to read_obs). If there was no data, we need to clean up
789 ! a bit of stuff, and read the famous last three integers that
790 ! have some QC information.
792 ! write(0,'("N=",I6," info%levels=",3I8)') &
793 ! obs_num, obs (obs_num) % info % levels
795 IF (error_ret .EQ. read_err ) THEN
797 IF (ASSOCIATED (obs(obs_num)%surface ) ) THEN
798 ! dealloc entire linked list if it exists
799 CALL dealloc_meas (obs(obs_num)%surface)
800 END IF
802 WRITE (UNIT = 0, FMT = '(A,A,A,1X,A)') &
803 "Troubles with measurement lines ", &
804 TRIM (obs(obs_num)%location%id), &
805 TRIM (obs(obs_num)%location%name)
807 bad_count = 0
808 io_error = 0
810 DO WHILE (io_error .GE. 0)
812 bad_count = bad_count + 1
814 IF ( bad_count .LT. 1000 ) THEN
816 ! READ (file_num , IOSTAT = io_error , FMT = meas_format) dummy_middle
817 READ (file_num , IOSTAT = io_error , FMT = meas_format) &
818 dummy_middle % pressure % data, dummy_middle % pressure % qc,&
819 dummy_middle % height % data, dummy_middle % height % qc,&
820 dummy_middle % temperature % data, dummy_middle % temperature % qc,&
821 dummy_middle % dew_point % data, dummy_middle % dew_point % qc,&
822 dummy_middle % speed % data, dummy_middle % speed % qc,&
823 dummy_middle % direction % data, dummy_middle % direction % qc,&
824 dummy_middle % u % data, dummy_middle % u % qc,&
825 dummy_middle % v % data, dummy_middle % v % qc,&
826 dummy_middle % rh % data, dummy_middle % rh % qc,&
827 dummy_middle % thickness % data, dummy_middle % thickness % qc
829 IF (eps_equal (dummy_middle%pressure%data, end_data_r , 1.)) THEN
831 READ (file_num , IOSTAT = io_error , FMT = end_format ) &
832 obs(obs_num)%info%num_vld_fld , &
833 obs(obs_num)%info%num_error , &
834 obs(obs_num)%info%num_warning
836 WRITE (UNIT = 0, FMT = '(A)') 'Starting to READ a new report.'
837 CYCLE read_obs
839 END IF
841 ELSE
843 WRITE (UNIT = 0, FMT = '(A)') &
844 'Too many attempts to read the measurement data correctly.',&
845 'Exiting read loop.'
847 CLOSE ( file_num )
848 IF (print_gts_read) CLOSE ( iunit )
850 EXIT read_obs
851 END IF
853 END DO
855 IF (io_error .LT. 0) THEN
856 CLOSE ( file_num )
857 IF (print_gts_read) CLOSE ( iunit )
858 EXIT read_obs
859 END IF
861 ELSE IF (error_ret .EQ. no_data .and. &
862 eps_equal(obs(obs_num)%ground%pw %data,missing_r,1.0) .and.&
863 obs(obs_num)%ground%tb19v%qc .ne. 0 .and. &
864 obs(obs_num)%ground%tb19h%qc .ne. 0 .and. &
865 obs(obs_num)%ground%tb22v%qc .ne. 0 .and. &
866 obs(obs_num)%ground%tb37v%qc .ne. 0 .and. &
867 obs(obs_num)%ground%tb37h%qc .ne. 0 .and. &
868 obs(obs_num)%ground%tb85v%qc .ne. 0 .and. &
869 obs(obs_num)%ground%tb85h%qc .ne. 0 .and. &
870 eps_equal(obs(obs_num)%ground%slp%data,missing_r,1.0) ) THEN
872 ! IF (print_gts_read ) THEN
873 ! WRITE (UNIT = 0 , FMT = '(A)' ) ' => NO DATA'
874 ! END IF
876 READ (file_num , IOSTAT = io_error , FMT = end_format ) &
877 obs(obs_num)%info%num_vld_fld , &
878 obs(obs_num)%info%num_error , &
879 obs(obs_num)%info%num_warning
881 READ (obs(obs_num) % info % platform (4:6), '(I3)') fm
883 CALL fm_decoder (fm, platform, &
884 synop=nsynops (icor+1), ship =nshipss (icor+1), &
885 metar=nmetars (icor+1), pilot=npilots (icor+1), &
886 sound=nsounds (icor+1), satem=nsatems (icor+1), &
887 satob=nsatobs (icor+1), airep=naireps (icor+1), &
888 gpspw=ngpspws (icor+1), gpszd=ngpsztd (icor+1), &
889 gpsrf=ngpsref (icor+1), gpsep=ngpseph (icor+1), &
890 ssmt1=nssmt1s (icor+1), &
891 ssmt2=nssmt2s (icor+1), ssmi =nssmis (icor+1), &
892 tovs =ntovss (icor+1), other=nothers (icor+1), &
893 amdar=namdars (icor+1), qscat=nqscats (icor+1), &
894 profl=nprofls (icor+1), buoy =nbuoyss (icor+1), &
895 bogus=nboguss (icor+1), airs =nairss (icor+1), tamdar =ntamdar (icor+1) )
897 IF ( ASSOCIATED (obs(obs_num)%surface)) THEN
898 ! dealloc entire linked list if it exists
899 CALL dealloc_meas ( obs(obs_num)%surface)
900 END IF
902 num_empty = num_empty + 1
904 CYCLE read_obs
906 END IF
908 ! We can compare the observation location with the spacial domain
909 ! and time window that the analysis will require.
910 ! If we are significantly outside, we toss out the observation.
912 IF (outside) THEN
914 ! IF (print_gts_read) THEN
915 ! WRITE (UNIT = 0 , FMT = '(A)' ) ' => OUTSIDE'
916 ! END IF
918 READ (file_num , IOSTAT = io_error , FMT = end_format ) &
919 obs(obs_num)%info%num_vld_fld , &
920 obs(obs_num)%info%num_error , &
921 obs(obs_num)%info%num_warning
923 READ (obs(obs_num) % info % platform (4:6), '(I3)') fm
925 CALL fm_decoder (fm, platform, &
926 synop=nsynops (icor+2), ship =nshipss (icor+2), &
927 metar=nmetars (icor+2), pilot=npilots (icor+2), &
928 sound=nsounds (icor+2), satem=nsatems (icor+2), &
929 satob=nsatobs (icor+2), airep=naireps (icor+2), &
930 gpspw=ngpspws (icor+2), gpszd=ngpsztd (icor+2), &
931 gpsrf=ngpsref (icor+2), gpsep=ngpseph (icor+2), &
932 ssmt1=nssmt1s (icor+2), &
933 ssmt2=nssmt2s (icor+2), ssmi =nssmis (icor+2), &
934 tovs =ntovss (icor+2), other=nothers (icor+2), &
935 amdar=namdars (icor+2), qscat=nqscats (icor+2), &
936 profl=nprofls (icor+2), buoy =nbuoyss (icor+2), &
937 bogus=nboguss (icor+2), airs =nairss (icor+2),tamdar =ntamdar (icor+2) )
939 IF ( ASSOCIATED (obs(obs_num)%surface)) THEN
940 ! dealloc entire linked list if it exists
941 CALL dealloc_meas ( obs(obs_num)%surface)
942 END IF
944 num_outside = num_outside + 1
946 CYCLE read_obs
948 ELSE
950 ! Elevation can sometimes be undefined
952 IF ( (obs(obs_num)%info%elevation .GT. (undefined1_r - 1.)) .OR. &
953 (obs(obs_num)%info%elevation .LT. (undefined2_r + 1.)) ) THEN
955 obs(obs_num)%info % elevation = missing_r
957 ! set elevation to 0.0 for marine reports, excluding Great Lakes
958 if ( fm .eq. 13 .or. fm .eq. 18 .or. fm .eq. 19 .or. &
959 fm .eq. 33 .or. fm .eq. 36 ) then
960 if ( obs(obs_num)%location%latitude .lt. 41. .or. &
961 obs(obs_num)%location%latitude .gt. 50. .or. &
962 obs(obs_num)%location%longitude .lt. -95. .or. & ! ncep used -93
963 obs(obs_num)%location%longitude .gt. -75. ) then
964 obs(obs_num)%info % elevation = 0.
965 end if
966 else if (fm < 39) then
967 m_miss = m_miss + 1
968 write(0,'(I7,1X,A,1X,A,1X,A,1X,A,1X,2(F8.3,A),A,1X,f11.3)')&
969 m_miss,'Missing elevation(id,name,platform,lat,lon,date,elv:', &
970 obs(obs_num)%location%id (1: 5),&
971 obs(obs_num)%location%name (1:20),&
972 obs(obs_num)%info%platform (1: 12),&
973 obs(obs_num)%location%latitude, 'N',&
974 obs(obs_num)%location%longitude,'E ', &
975 obs(obs_num)%valid_time%date_char, &
976 obs(obs_num)%info % elevation
977 endif
979 ! assigning elevation info for ships and buoys located in the Great
980 ! Lakes.
981 !http://www.nco.ncep.noaa.gov/pmb/codes/nwprod/decoders/decod_dcmsfc/sorc/maelev.f
982 xla = obs(obs_num)%location%latitude
983 xlo = obs(obs_num)%location%longitude
984 if ( fm .eq. 13 .or. fm .eq. 33 .or. fm .eq. 36 ) then
985 ! for ships
986 if ( ( xlo .ge. -92.5 .and. xlo .le. -84.52 ) .and. &
987 ( xla .ge. 46.48 .and. xla .le. 49.0 ) ) then
988 !Ship located in Lake Superior
989 obs(obs_num)%info % elevation = 183.
990 else if ( ( xlo .ge. -88.1 .and. xlo .le. -84.8 ) .and. &
991 ( xla .ge. 41.2 .and. xla .le. 46.2 ) ) then
992 !Ship located in Lake Michigan
993 obs(obs_num)%info % elevation = 176.
994 else if ( ( xlo .ge. -84.8 .and. xlo .le. -79.79 ) .and. &
995 ( xla .ge. 43.0 .and. xla .le. 46.48 ) ) then
996 !Ship located in Lake Huron or Georgian Bay
997 obs(obs_num)%info % elevation = 176.
998 else if ( ( xlo .ge. -84.0 .and. xlo .le. -78.0 ) .and. &
999 ( xla .ge. 41.0 .and. xla .le. 42.9 ) ) then
1000 !Ship located in Lake Erie
1001 obs(obs_num)%info % elevation = 174.
1002 else if ( ( xlo .ge. -80.0 .and. xlo .le. -76.0 ) .and. &
1003 ( xla .ge. 43.1 .and. xla .le. 44.23 ) ) then
1004 !Ship located in Lake Ontario
1005 obs(obs_num)%info % elevation = 74.
1006 end if
1007 end if !end if ships
1008 if ( fm .eq. 18 .and. obs(obs_num)%location%id(1:2) .eq. '45' ) then
1009 ! for Great Lakes fixed buoys
1010 ! get station elevation from station table if available
1011 if ( use_msfc_tbl .and. num_stations_msfc > 0 ) then
1012 ipos = 0
1013 do while ( ipos < num_stations_msfc )
1014 ipos = ipos + 1
1015 if ( obs(obs_num)%location%id(1:5) == id_tbl(ipos) ) then
1016 obs(obs_num)%info % elevation = elev_tbl(ipos)
1017 exit
1018 end if
1019 end do
1020 end if !table info available
1021 end if ! end if buoys
1022 END IF ! missing elevation
1024 END IF ! not outside
1026 ! for gts_from_mmm_archive:
1027 ! METARs (FM-15) never report surface pressure. METARs usually
1028 ! report altimeter setting. The decoder places the altimeter
1029 ! setting into the little_r surface pressure data record.
1030 ! U.S. and Canadian METARs often contain SLP and temperatures
1031 ! in tenths of degree in the remarks section (RMK). These are
1032 ! decoded. JFB
1034 IF ((obs (obs_num)%info%platform(1:12) .EQ. 'FM-15 METAR ' .or. &
1035 obs (obs_num)%info%platform(1:12) .EQ. 'FM-16 SPECI ') .AND. &
1036 (ASSOCIATED (obs (obs_num)%surface ) ) ) THEN
1037 if ( calc_psfc_from_QNH .and. gts_from_mmm_archive ) then
1038 if ( obs(obs_num)%ground%psfc%data > 0.0 .and. &
1039 obs(obs_num)%info%elevation > 0.0 ) then
1040 QNH = obs(obs_num)%ground%psfc%data * 0.01 ! Pa to hPa
1041 elev = obs(obs_num)%info%elevation
1042 obs(obs_num)%ground%psfc%data = psfc_from_QNH(QNH,elev) &
1043 * 100.0 ! hPa to Pa
1044 obs(obs_num)%ground%psfc%qc = 0
1045 if ( associated(obs(obs_num)%surface) ) then
1046 obs(obs_num)%surface%meas%pressure%data = &
1047 obs(obs_num)%ground%psfc%data
1048 obs(obs_num)%surface%meas%pressure%qc = &
1049 obs(obs_num)%ground%psfc%qc
1050 end if ! associated data
1051 end if ! valid QNH and elev
1052 else
1053 obs(obs_num)%ground%psfc%data = missing_r
1054 obs(obs_num)%ground%psfc%qc = missing
1055 end if ! calc_psfc_from_QNH and gts_from_mmm_archive
1056 END IF ! metar
1058 ! for gts_from_mmm_archive:
1059 ! For ship data (FM-13 or FM-18), the elevation is missing and
1060 ! surface (or station) pressure is set to 1013.01 hPa. This
1061 ! value is a flag to indicate a suface report (the elevation of
1062 ! a ship or buoy does not have to be at sea level as in the case
1063 ! of Lake Michigan). JFB
1065 ! If this is a ship observation, we need to define the elevation
1066 ! as identical to the geopotential height, and set the height QC flag
1067 ! to ok. This is the only way to get SHIP data into the surface
1068 ! analysis. Since we are at sea level, we also set the pressure
1069 ! to equal to the sea level pressure.
1071 IF ( (obs(obs_num)%info%platform(1:10) == 'FM-13 SHIP') .or. &
1072 (obs(obs_num)%info%platform(1:10) == 'FM-18 BUOY') ) then
1073 if ( ASSOCIATED(obs(obs_num)%surface) ) then
1074 obs(obs_num)%surface%meas%height%data = &
1075 obs(obs_num)%info%elevation
1076 obs(obs_num)%surface%meas%height%qc = 0
1077 if ( (obs(obs_num)%info%elevation == 0.0) ) then
1078 !obs(obs_num)%surface%meas%height%data = &
1079 ! obs(obs_num)%info%elevation
1080 !obs(obs_num)%surface%meas%height%qc = 0
1081 obs(obs_num)%surface%meas%pressure%data = &
1082 obs(obs_num)%ground%slp%data
1083 obs(obs_num)%surface%meas%pressure%qc = 0
1084 else
1085 if ( eps_equal(obs(obs_num)%surface%meas%pressure%data, &
1086 101301.000, 1.) ) then
1087 ! replace 1013.01 with missing value
1088 obs(obs_num)%surface%meas%pressure%data = missing_r
1089 obs(obs_num)%surface%meas%pressure%qc = missing
1090 end if
1091 end if ! elev is 0
1092 end if ! has associated data
1093 end if ! end if ship or buoy
1095 ! FENG GAO 03/07/2014
1096 ! QuikSCAT nominal mission ended on November 23, 2009
1097 ! NOW FM-281 denote qscat only for research and
1098 ! ASCAT 10-m ocean surface wind for research and application
1100 IF ((obs (obs_num)%info%platform(1:6) .EQ. 'FM-281' ) .and. &
1101 (ASSOCIATED (obs (obs_num)%surface ) ) ) THEN
1102 if (obs(obs_num)%info%elevation .LT. 0.0) then
1103 obs(obs_num)%surface%meas%height%data = 10.0
1104 else
1105 obs(obs_num)%surface%meas%height%data = &
1106 obs(obs_num)%info%elevation
1107 end if
1108 obs(obs_num)%surface%meas%height%qc = 0
1109 END IF
1111 ! for gts_from_mmm_archive:
1112 ! SYNOP reports (FM-12) have a flag of 1013.01 in the surface
1113 ! pressure field to indicate surface data. However, the SYNOP
1114 ! code does include a station pressure group. The observed
1115 ! station pressure will be decoded and replace 1013.01 if it
1116 ! is available. JFB
1118 ! YRG 04/04/2009
1119 ! For SYNOP, if surface%meas%pressure%data = 101301.000
1120 ! (101301.000 is a fake value in NCAR archived LITTLE_R file)
1121 ! and the slp is missing (note if SLP is available, WRFVar
1122 ! will use the SLP to derive Psfc and ignore the original Psfc,
1123 ! see da_tools/da_obs_sfc_correction.inc),
1124 ! fill in surface%meas%pressure%data with ground%psfc%data:
1126 ! hcl-note: disagrees with the above slp to psfc procedure.
1127 ! if both slp and psfc are available for synop reports,
1128 ! psfc should be used in DA system
1129 ! hcl-note2: 101301 is simply a flag, real reports do not
1130 ! have 1/100 precision. Just replace 101301 with missing value
1132 IF ( (obs(obs_num)%info%platform(1:5).EQ.'FM-12') .and. &
1133 (ASSOCIATED(obs(obs_num)%surface)) ) THEN
1134 if ( eps_equal(obs(obs_num)%surface%meas%pressure%data, &
1135 101301.000, 1.) ) then
1136 obs(obs_num)%surface%meas%pressure%data = missing_r
1137 obs(obs_num)%surface%meas%pressure%qc = missing
1138 endif
1139 ENDIF !end if FM-12 synop
1141 ! This may be wasted print-out, but it is comforting to see.
1143 ! IF (print_gts_read ) THEN
1145 ! IF (obs(obs_num)%info%is_sound) THEN
1146 ! WRITE (UNIT = 0, FMT = '(A)') ' => SOUNDING'
1147 ! ELSE
1148 ! WRITE (UNIT = 0, FMT = '(A)') ' => SURFACE'
1149 ! ENDIF
1151 ! END IF
1153 END IF
1155 ! We have now read ground info and soundings, what follows in the
1156 ! standard format are three integers describing information gleaned
1157 ! from the program that generated the observational data.
1159 READ (file_num , IOSTAT = io_error , FMT = end_format ) &
1160 obs(obs_num)%info%num_vld_fld , &
1161 obs(obs_num)%info%num_error , &
1162 obs(obs_num)%info%num_warning
1164 ! Once again, after a read, was it successful. If not toss the whole
1165 ! thing out (this is handled through the dealloc_meas routine if any
1166 ! upper-air data was encountered). Discarding all of the ingested
1167 ! data may be a bit much, which is why the error print-out is provided.
1168 ! After the error is processed, the reading process is re-started.
1170 IF ( io_error .NE. 0 ) THEN
1172 error_message = &
1173 'Error trying to read last 3 integers in observation '&
1174 // TRIM ( obs(obs_num)%location%id ) &
1175 // TRIM ( obs(obs_num)%location%name ) // '.'
1177 CALL error_handler (proc_name, error_message, &
1178 ' Discarding entire and continuing.', .FALSE.)
1180 IF ( ASSOCIATED ( obs(obs_num)%surface ) ) THEN
1181 CALL dealloc_meas ( obs(obs_num)%surface)
1182 END IF
1184 CYCLE read_obs
1186 END IF
1188 ! Before we leave this loop, we make sure the surface level is the first
1189 ! level in the sounding.
1191 ! This test might results in removing levels
1193 IF (.NOT. eps_equal(obs(obs_num)%info%elevation, missing_r, 1.)) &
1194 CALL surf_first ( obs(obs_num)%surface , obs(obs_num)%info%elevation )
1196 ! Height and pressure when missing were replaced by ICAO values, restore
1198 CALL missing_hp ( obs(obs_num)%surface )
1200 ! Update the info % levels after surf_first test
1202 obs (obs_num) % info % levels = info_levels (obs(obs_num)%surface)
1204 ! Update the info % is_sound on the updated info % levels
1206 IF (obs (obs_num) % info % levels .GT. 1) THEN
1207 obs (obs_num) % info % is_sound = .TRUE.
1208 ELSE IF (obs (obs_num) % info % levels .EQ. 1) THEN
1209 obs (obs_num) % info % is_sound = .FALSE.
1211 ! Guo, 02/26/2003 -- When slp or pw is not available, then ...
1213 ELSE IF (eps_equal(obs(obs_num)%ground%pw %data,missing_r,1.0) .and.&
1214 obs(obs_num)%ground%tb19v%qc .ne. 0 .and. &
1215 obs(obs_num)%ground%tb19h%qc .ne. 0 .and. &
1216 obs(obs_num)%ground%tb22v%qc .ne. 0 .and. &
1217 obs(obs_num)%ground%tb37v%qc .ne. 0 .and. &
1218 obs(obs_num)%ground%tb37h%qc .ne. 0 .and. &
1219 obs(obs_num)%ground%tb85v%qc .ne. 0 .and. &
1220 obs(obs_num)%ground%tb85h%qc .ne. 0 .and. &
1221 eps_equal(obs(obs_num)%ground%slp%data,missing_r,1.0) ) THEN
1222 ! Station are expected with at least one level here
1223 obs(obs_num) % info % discard = .TRUE.
1224 CYCLE read_obs
1225 ENDIF
1227 ! Have read observation and its measurement(s) without error
1228 ! so continue to next observation.
1230 obs_num = obs_num + 1
1232 END DO read_obs
1235 ! The end of the observation file has been reached. Decrement the counter
1236 ! to get the total number of observations successfully read by the program.
1237 ! Output this information to the outside world. We can also provide
1238 ! the information on the observations that are NOT included in the analysis.
1240 obs_num = obs_num - 1
1242 ! Print out number of read observations per type
1244 WRITE (UNIT = 0, FMT = '(/,A)') &
1245 '------------------------------------------------------------------------------'
1246 WRITE (UNIT = 0, FMT = '(A)') 'GTS OBSERVATIONS READ:'
1248 WRITE (UNIT = 0, FMT = '(A)')
1249 WRITE (UNIT = 0, FMT = '(A,I7)') ' SYNOP reports:',nsynops (0)
1250 WRITE (UNIT = 0, FMT = '(A,I7)') ' SHIPS reports:',nshipss (0)
1251 WRITE (UNIT = 0, FMT = '(A,I7)') ' BUOYS reports:',nbuoyss (0)
1252 WRITE (UNIT = 0, FMT = '(A,I7)') ' BOGUS reports:',nboguss (0)
1253 WRITE (UNIT = 0, FMT = '(A,I7)') ' METAR reports:',nmetars (0)
1254 WRITE (UNIT = 0, FMT = '(A,I7)') ' PILOT reports:',npilots (0)
1255 WRITE (UNIT = 0, FMT = '(A,I7)') ' SOUND reports:',nsounds (0)
1256 WRITE (UNIT = 0, FMT = '(A,I7)') ' AMDAR reports:',namdars (0)
1257 WRITE (UNIT = 0, FMT = '(A,I7)') ' SATEM reports:',nsatems (0)
1258 WRITE (UNIT = 0, FMT = '(A,I7)') ' SATOB reports:',nsatobs (0)
1259 WRITE (UNIT = 0, FMT = '(A,I7)') ' GPSPW reports:',ngpspws (0)
1260 WRITE (UNIT = 0, FMT = '(A,I7)') ' GPSZD reports:',ngpsztd (0)
1261 WRITE (UNIT = 0, FMT = '(A,I7)') ' GPSRF reports:',ngpsref (0)
1262 WRITE (UNIT = 0, FMT = '(A,I7)') ' GPSEP reports:',ngpseph (0)
1263 WRITE (UNIT = 0, FMT = '(A,I7)') ' AIREP reports:',naireps (0)
1264 WRITE (UNIT = 0, FMT = '(A,I7)') 'TAMDAR reports:',ntamdar (0)
1265 WRITE (UNIT = 0, FMT = '(A,I7)') ' SSMT1 reports:',nssmt1s (0)
1266 WRITE (UNIT = 0, FMT = '(A,I7)') ' SSMT2 reports:',nssmt2s (0)
1267 WRITE (UNIT = 0, FMT = '(A,I7)') ' SSMI reports:',nssmis (0)
1268 WRITE (UNIT = 0, FMT = '(A,I7)') ' TOVS reports:',ntovss (0)
1269 WRITE (UNIT = 0, FMT = '(A,I7)') ' QSCAT reports:',nqscats (0)
1270 WRITE (UNIT = 0, FMT = '(A,I7)') ' PROFL reports:',nprofls (0)
1271 WRITE (UNIT = 0, FMT = '(A,I7)') ' AIRST reports:',nairss (0)
1272 WRITE (UNIT = 0, FMT = '(A,I7)') ' OTHER reports:',nothers (0)
1273 WRITE (UNIT = 0, FMT = '(A,I7)') ' Total reports:', &
1274 nsynops (0) + nshipss (0) + nmetars (0) + npilots (0) + nsounds (0)+&
1275 nsatems (0) + nsatobs (0) + naireps (0) + ntamdar (0)+ ngpspws (0) + ngpsztd (0)+&
1276 ngpsref (0) + ngpseph (0) + &
1277 nssmt1s (0) + nssmt2s (0) + nssmis (0) + ntovss (0) + nboguss (0)+&
1278 nothers (0) + namdars (0) + nqscats (0) + nprofls(0) + nbuoyss(0) +&
1279 nairss(0)
1281 ! Print number of observation ingested
1283 ! WRITE (0,'(/,A)') &
1284 !'------------------------------------------------------------------------------'
1285 ! WRITE (UNIT = 0, FMT = '(A)') 'INGESTED GTS OBSERVATIONS:'
1287 WRITE (UNIT = 0, FMT = '(/,4(A,i8,/))' ) &
1289 "Number of observations read: ",obs_num+ &
1290 num_empty+ &
1291 num_outside-&
1292 n_obs, &
1293 "Number of empty observations: ",num_empty, &
1294 "Number of observations out of domain: ",num_outside,&
1295 "Number of observations for ingestion: ",obs_num-n_obs
1297 ! Total number of observation accumulated
1299 n_obs = obs_num
1301 write(0,'(/"AIRCRAFT DATA: Total=",I7," Above cut_height=",I7)')&
1302 N_air, N_air_cut
1303 contains
1305 subroutine print_extra_obs
1307 READ (obs(obs_num) % info % platform (4:6), '(I3)', &
1308 IOSTAT = platform_error) fm
1309 ! synop 12,14 'SYNOP','SYNOP MOBIL'
1310 if (fm == 12 .or. fm ==14) return
1312 ! ship 13 'SHIP'
1313 if (fm == 13) return
1315 ! metar 15,16 'METAR','SPECI'
1316 if (fm == 15 .or. fm == 16) return
1318 ! buoy 18 'BUOY'
1319 if (fm == 18 .or. fm == 19) return
1321 ! pilot 32,33,34 'PILOT','PILOT SHIP','PILOT MOBIL'
1322 if (fm >= 32 .and. fm <= 34) return
1324 ! sound 35,36,37,38 'TEMP','TEMP SHIP, 'TEMP DROP','TEMP MOBIL'
1325 if (fm >= 35 .and. fm <= 38) return
1327 ! amdar 42 'AMDAR'
1328 if (fm == 42) return
1330 ! satem 86 'SATEM'
1331 if (fm == 86) return
1333 ! satob 88 'SATOB'
1334 if (fm == 88) return
1336 ! airep 96,97 'AIREP'
1337 if (fm == 96 .or. fm == 97) return
1339 ! tamdar 101 'TAMDAR'
1340 if (fm == 101) return
1342 ! gpspw 111 'GPSPW'
1343 if (fm == 111) return
1345 ! gpsztd 114 'GPSZD'
1346 if (fm == 114) return
1348 ! gpsref 116 'GPSRF'
1349 if (fm == 116) return
1351 ! gpseph 118 'GPSEP'
1352 if (fm == 118) return
1354 ! ssmt1 121 'SSMT1'
1355 if (fm == 121) return
1357 ! ssmt2 122 'SSMT2'
1358 if (fm == 122) return
1360 ! ssmi 125,126 'SSMI'
1361 if (fm == 125 .or. fm == 126) return
1363 ! tovs 131 'TOVS'
1364 if (fm == 131) return
1366 ! qscat 281 'Quikscat'
1367 if (fm == 281) return
1369 ! profl 132 'Profilers'
1370 if (fm == 132) return
1372 ! bogus 135 'Bogus'
1373 if (fm == 135) return
1374 ! AIRSRET 133 'Bogus'
1375 if (fm == 133) return
1378 ! other Any other code 'UNKNOWN'
1380 num_unknown = num_unknown + 1
1381 write(0,'(2I8," ID=",a," Name=",a," Platform=",a)') &
1382 num_unknown, obs_num, &
1383 obs(obs_num)%location % id(1:15), &
1384 obs(obs_num)%location % name, &
1385 obs(obs_num)%info % platform
1386 end subroutine print_extra_obs
1388 END SUBROUTINE read_obs_gts
1390 !\f
1391 !---------------------------------------------------------------------------
1393 SUBROUTINE read_measurements (file_num, surface, location, info, bad_data, &
1394 error, ins, jew, &
1395 map_projection, elevation, nlevels, iunit,&
1396 print_gts_read)
1398 ! This routine reads in 'measurements' at as many levels as there are in
1399 ! the report, then stops and returns when an end-of-measurements flag is
1400 ! found. If any reads produce error, return error code which causes entire
1401 ! observation to be discarded (ob is not discarded on eof error).
1403 USE module_icao
1405 IMPLICIT NONE
1407 INTEGER , INTENT ( IN ) :: file_num ! file to read
1408 TYPE ( measurement ) , POINTER :: surface ! ptr to 1st msmt
1409 TYPE ( location_type ) , INTENT ( IN ) :: location ! 5 digit ID, name
1410 TYPE ( source_info ) , INTENT ( INOUT ) :: info ! 5 digit ID, name
1411 LOGICAL , INTENT ( IN ) :: bad_data ! read, not store
1412 INTEGER , INTENT ( OUT ) :: error ! err and type
1414 INTEGER :: ins , jew, k
1415 INTEGER :: map_projection
1417 CHARACTER ( LEN = 32 ) , PARAMETER :: proc_name = 'read_measurements'
1418 INTEGER :: meas_count
1419 INTEGER :: io_error
1420 TYPE ( measurement ) , POINTER :: current
1422 CHARACTER ( LEN = 40 ) :: location_id , &
1423 location_name
1424 REAL , INTENT(IN) :: elevation
1425 REAL :: new_press, new_heightt, &
1426 ref_h
1427 INTEGER, INTENT (out) :: nlevels
1428 INTEGER, INTENT (in) :: iunit
1429 LOGICAL, INTENT (in) :: print_gts_read
1430 LOGICAL :: no_height, no_pressure
1431 LOGICAL :: no_temperature
1433 INTEGER :: icrs, fm
1435 !------------------------------------------------------------------------------!
1437 ! INCLUDE 'error.inc'
1438 ! INTERFACE
1439 ! INCLUDE 'error.int'
1440 ! END INTERFACE
1442 ! Initialize dummy pointers and counters and observation names, and such.
1444 ALLOCATE ( current )
1445 NULLIFY ( current%next )
1446 NULLIFY ( surface )
1447 error = ok
1448 meas_count = 0
1449 location_id = TRIM ( location%id )
1450 location_name = TRIM ( location%name )
1452 ! This loop continues until either an error occurs, or until the end of
1453 ! the measurement tag is found (the graceful exit).
1455 read_meas: DO
1457 ! Currently, this read puts in 12 pairs of data, a real observation
1458 ! value and the accompanying QC flag.
1460 !FV READ ( file_num , IOSTAT = io_error , FMT = meas_format ) &
1461 !FV current%meas
1463 READ ( file_num , IOSTAT = io_error , FMT = meas_format ) &
1464 current % meas % pressure % data, current % meas % pressure % qc, &
1465 current % meas % height % data, current % meas % height % qc, &
1466 current % meas % temperature % data, current % meas % temperature % qc, &
1467 current % meas % dew_point % data, current % meas % dew_point % qc, &
1468 current % meas % speed % data, current % meas % speed % qc, &
1469 current % meas % direction % data, current % meas % direction % qc, &
1470 current % meas % u % data, current % meas % u % qc, &
1471 current % meas % v % data, current % meas % v % qc, &
1472 current % meas % rh % data, current % meas % rh % qc, &
1473 current % meas % thickness % data, current % meas % thickness % qc
1475 ! An error < 0 means the end of the file (usually), and an error > 0
1476 ! is just a broken read. Describe the read error so that the calling
1477 ! routine knows what happened, then exit this loop (which is exiting
1478 ! this routine, basically).
1480 IF (io_error .GT. 0 ) THEN
1481 error = read_err
1482 ! CLOSE ( file_num )
1483 EXIT read_meas
1484 ELSE IF (io_error .LT. 0 ) THEN
1485 error = eof_err
1486 CLOSE ( file_num )
1487 IF (print_gts_read) CLOSE ( iunit )
1488 EXIT read_meas
1489 END IF
1491 ! If we know a priori that this data is bad, no tests are necessary on
1492 ! the various flags values.
1494 bad_loop_1 : IF (.NOT. bad_data) THEN
1496 ! A successful read, yahoo! As the data may not have the flags
1497 ! set up the way we want, go through directly after this read
1498 ! and make sure that any special values are all set to missing.
1500 IF ((current%meas%pressure%data .GE. ( undefined1_r - 10. ) ) .OR. &
1501 (current%meas%pressure%data .LE. ( undefined2_r + 10. ) ) .OR. &
1502 (current%meas%pressure%data .LE. 0.0) )THEN
1503 current%meas%pressure%data = missing_r
1504 current%meas%pressure%qc = missing
1505 END IF
1506 IF ((current%meas%height%data .GT. ( undefined1_r - 1. ) ) .OR. &
1507 (current%meas%height%data .LT. ( undefined2_r + 1. ) ) .OR. &
1508 (current%meas%height%data .GT. ( height_max_icao - 1.)) .OR. &
1509 (current%meas%height%data .GT. ( ABS (missing_r) - 1. ))) THEN
1510 current%meas%height%data = missing_r
1511 current%meas%height%qc = missing
1512 END IF
1513 IF ((current%meas%temperature%data .GT. ( undefined1_r - 1. ) ) .OR. &
1514 (current%meas%temperature%data .LT. ( undefined2_r + 1. ) ) ) THEN
1515 current%meas%temperature%data = missing_r
1516 current%meas%temperature%qc = missing
1517 END IF
1518 IF (current%meas%temperature%data .GT. ( 99999.0 - 1. ) ) THEN
1519 current%meas%temperature%data = missing_r
1520 current%meas%temperature%qc = missing
1521 END IF
1522 IF ((current%meas%dew_point%data .GT. ( undefined1_r - 1. ) ) .OR. &
1523 (current%meas%dew_point%data .LT. ( undefined2_r + 1. ) ) ) THEN
1524 current%meas%dew_point%data = missing_r
1525 current%meas%dew_point%qc = missing
1526 END IF
1527 IF ((current%meas%speed%data .GT. ( undefined1_r - 1. ) ) .OR. &
1528 (current%meas%speed%data .LT. ( undefined2_r + 1. ) ) ) THEN
1529 current%meas%speed%data = missing_r
1530 current%meas%speed%qc = missing
1531 END IF
1532 IF ((current%meas%direction%data .GT. ( undefined1_r - 1. ) ) .OR. &
1533 (current%meas%direction%data .LT. ( undefined2_r + 1. ) ) ) THEN
1534 current%meas%direction%data = missing_r
1535 current%meas%direction%qc = missing
1536 END IF
1537 IF ((current%meas%u%data .GT. ( undefined1_r - 1. ) ) .OR. &
1538 (current%meas%u%data .LT. ( undefined2_r + 1. ) ) ) THEN
1539 current%meas%u%data = missing_r
1540 current%meas%u%qc = missing
1541 END IF
1542 IF ((current%meas%v%data .GT. ( undefined1_r - 1. ) ) .OR. &
1543 (current%meas%v%data .LT. ( undefined2_r + 1. ) ) ) THEN
1544 current%meas%v%data = missing_r
1545 current%meas%v%qc = missing
1546 END IF
1547 IF ((current%meas%rh%data .GT. ( undefined1_r - 1. ) ) .OR. &
1548 (current%meas%rh%data .LT. ( undefined2_r + 1. ) ) ) THEN
1549 current%meas%rh%data = missing_r
1550 current%meas%rh%qc = missing
1551 END IF
1552 IF ((current%meas%thickness%data .GT. ( undefined1_r - 1. ) ) .OR. &
1553 (current%meas%thickness%data .LT. ( undefined2_r + 1. ) ) ) THEN
1554 current%meas%thickness%data = missing_r
1555 current%meas%thickness%qc = missing
1556 END IF
1558 current%meas%qv%data = missing_r
1559 current%meas%qv%qc = missing
1561 END IF bad_loop_1
1563 ! The data we just read in could have been the flag for the end of
1564 ! the measurement. This is the graceful way to exit this routine.
1565 ! If this is the end of the measurement section for this observation,
1566 ! set all of the data to the same end of measurement value,
1567 ! just in case there were some stray unset values in the generating
1568 ! program.
1570 IF (eps_equal (current%meas%pressure%data , end_data_r , 1.) .OR. &
1571 eps_equal (current%meas%height%data , end_data_r , 1.)) THEN
1572 current%meas%pressure%data = end_data_r
1573 current%meas%height%data = end_data_r
1574 current%meas%temperature%data = end_data_r
1575 current%meas%dew_point%data = end_data_r
1576 current%meas%speed%data = end_data_r
1577 current%meas%direction%data = end_data_r
1578 current%meas%u%data = end_data_r
1579 current%meas%v%data = end_data_r
1580 current%meas%rh%data = end_data_r
1581 current%meas%thickness%data = end_data_r
1582 current%meas%pressure%qc = end_data
1583 current%meas%height%qc = end_data
1584 current%meas%temperature%qc = end_data
1585 current%meas%dew_point%qc = end_data
1586 current%meas%speed%qc = end_data
1587 current%meas%direction%qc = end_data
1588 current%meas%u%qc = end_data
1589 current%meas%v%qc = end_data
1590 current%meas%rh%qc = end_data
1591 current%meas%thickness%qc = end_data
1592 current%meas%qv%data = end_data
1593 current%meas%qv%qc = end_data
1594 error = ok
1596 EXIT read_meas
1598 ! Don't copy record if either press. or height missing, and
1599 ! wind, temp, dew point and rh are all missing
1601 ELSEIF ((eps_equal(current%meas%pressure%data, missing_r , 1.) .OR. &
1602 eps_equal(current%meas%pressure%data, missing_r , 1.)) .AND. &
1603 eps_equal (current%meas%speed%data, missing_r , 1.) .AND. &
1604 eps_equal (current%meas%direction%data, missing_r , 1.) .AND. &
1605 eps_equal (current%meas%temperature%data, missing_r , 1.) .AND. &
1606 eps_equal (current%meas%dew_point%data, missing_r , 1.) .AND. &
1607 eps_equal (current%meas%rh%data, missing_r , 1.)) THEN
1609 CYCLE read_meas
1611 END IF
1613 ! If this is bad data, we needed to make sure that the ending measurement
1614 ! is the famous end_data flag so that we hit a correct exit from this
1615 ! loop and routine. We can just cycle the read loop again.
1617 IF (bad_data) THEN
1618 CYCLE read_meas
1619 END IF
1621 !
1622 ! If both pressure and height are missing, throw out data at this level
1623 !
1625 IF ((eps_equal ( current%meas%pressure%data , missing_r , 1.)) .AND. &
1626 (eps_equal ( current%meas%height %data , missing_r , 1.))) THEN
1627 CYCLE read_meas
1628 END IF
1630 if ( print_gts_read ) then
1631 IF (( eps_equal (current%meas%dew_point%data , missing_r , 1.)) .AND.&
1632 (.NOT.eps_equal (current%meas%rh %data , missing_r , 1.))) THEN
1633 WRITE (iunit,'(A,F10.2,/,A,F10.2)') &
1634 " Td = ",current%meas%dew_point%data,&
1635 " Rh = ",current%meas%rh%data
1636 ENDIF
1637 end if
1639 !
1640 ! Assign the SSMI error (AFWA only)
1641 !
1643 ! initialize the variable that might be used in module_err_ncep.F90
1644 ! for checking if the error is pre-assigned
1645 current%meas%speed%error = 0.0
1647 READ (info % platform (4:6), '(I3)') fm
1649 IF ((fm .EQ. 125) .AND. (current%meas%speed%qc .GT. missing)) THEN
1651 SELECT CASE (current%meas%speed%qc)
1653 CASE (0)
1654 current%meas%speed%error = 2. !m/s
1655 CASE (1)
1656 current%meas%speed%error = 5. !m/s
1657 CASE (2)
1658 current%meas%speed%error = 10. !m/s
1659 CASE (3)
1660 current%meas%speed%error = 20. !m/s
1661 CASE DEFAULT
1662 current%meas%speed%error = 20. !m/s
1663 END SELECT
1665 current%meas%speed%qc = 0
1667 ELSE IF ((fm == 97 .or. fm == 96 .or. fm == 42) .and. &
1668 (current%meas%height%qc == 0 ) ) then
1670 N_air = N_air + 1
1671 if (current%meas%height%data > aircraft_cut) then
1673 ! To convert the Aircraft observed height (> cutoff_height=3000m) to pressure:
1674 ! and discarded the observed height:
1675 N_air_cut = N_air_cut + 1
1676 call Aircraft_pressure(current%meas%height, current%meas%pressure)
1677 endif
1679 ! Y.-R. Guo, 03/20/2008: In RTOBS 2006091300 data:obs.2006091300.gz, there are
1680 ! two levels obs in FM-13 SHIP causing troubles in wrfvar.
1681 ! SHIP and BUOY, if pressure < 85000.0 Pa, discarded.
1682 ELSE IF ( fm == 13 .or. fm == 18 .or. fm == 19 ) THEN
1683 if (current%meas%pressure%data < 85000.0 .and. &
1684 current%meas%pressure%qc >= 0) then
1685 write(0,'(a,3x,a,2x,a,2x,2f13.5,2x,"Pressure=",f10.1,a,i8)') &
1686 'Discarded:', info%platform(1:12), trim(location%id), &
1687 location%latitude, location%longitude, &
1688 current%meas%pressure%data, " < 85000.0 Pa, qc=", &
1689 current%meas%pressure%qc
1690 CYCLE read_meas
1691 endif
1693 ENDIF
1695 ! Some pressure and height is needed for vertically inserting data
1697 IF (ASSOCIATED (current)) THEN
1699 ! Guo 01/26/2004: Very grossly h/p consistency check in case of both p and h
1700 ! are reported as good data (qc=0):
1701 ! hcl-note: disagree with checking qc in input data
1702 ! hcl-note: qc should be assigned by obsproc after
1703 ! hcl-note: applying quality check (if there is any)
1704 !-----------------------------------------------------------------------------
1705 ! Do no perform gross check on height for
1706 ! Sounde (FM-35) & AIRS (FM-133) retrievals profile data
1707 ! 07/07/2006 Syed RH Rizvi
1708 !
1709 ! Why ???? YRG modified again 11/09/2006
1710 !-----------------------------------------------------------------------------
1711 IF ( current%meas%height%qc == 0 .and. current%meas%pressure%qc == 0 &
1712 .and. .NOT.eps_equal(current%meas%height %data, missing_r, 1.) .and.&
1713 .NOT.eps_equal(current%meas%pressure%data, missing_r, 1.) )THEN
1715 ! ......if Pressure < 500 Pa, the reference height will be unrealistically
1716 ! decreased from the function of "Ref_height":
1718 if (current%meas%pressure%data >= 500.0) then
1719 Ref_h = Ref_height (current%meas%pressure%data)
1720 else
1721 Ref_h = Ref_height (500.0)
1722 endif
1723 ! ..........if the difference between the reported height and the reference
1724 ! height is greater than 12000m, discarded this level data:
1726 if (abs(Ref_h-current%meas%height%data) > 12000) then
1727 write(0,'("??? Pressure or Height reported inconsistent:")')
1728 write(0,'(3x,a,2x,a,2x,2f13.5)') &
1729 info%platform(1:12), trim(location%id), &
1730 location%latitude, location%longitude
1731 write(0,'("(height-Ref_h) > 12000m, p,h,ref_h:",3e15.5/)') &
1732 current%meas%pressure%data, current%meas%height%data, Ref_h
1733 CYCLE read_meas
1734 endif
1735 ENDIF
1737 IF (eps_equal (current%meas%pressure%data , missing_r , 1.)) THEN
1738 !hcl what is this?
1739 if (current%meas%height%data > (htop+100.)) CYCLE read_meas
1740 current%meas%pressure%data = Ref_pres (current%meas%height%data)
1741 current%meas%pressure%qc = missing
1742 ENDIF
1743 IF (eps_equal (current%meas%height%data , missing_r , 1.)) THEN
1744 !hcl what is this?
1745 if (current%meas%pressure%data < (ptop-500.)) CYCLE read_meas
1746 current%meas%height%data = Ref_height (current%meas%pressure%data)
1747 current%meas%height%qc = missing
1748 ENDIF
1749 ENDIF
1751 IF (ASSOCIATED (surface)) THEN
1752 IF (eps_equal (surface%meas%pressure%data , missing_r , 1.)) THEN
1753 surface%meas%pressure%data = Ref_pres (surface%meas%height%data)
1754 surface%meas%pressure%qc = missing
1755 ENDIF
1757 IF (eps_equal (surface%meas%height%data , missing_r , 1.)) THEN
1758 surface%meas%height%data = Ref_height (surface%meas%pressure%data)
1759 surface%meas%height%qc = missing
1760 ENDIF
1761 ENDIF
1763 ! Since it seems that everything went ok, insert this measurement ordered
1764 ! by pressure.
1766 meas_count = meas_count + 1
1768 ! Insert now
1770 CALL insert_at (surface , current , elevation)
1772 ! One level has been sucessfuly read
1774 nlevels = nlevels + 1
1777 ! Allocate space for another measurement, so we can go try and
1778 ! read another level in this routine.
1779 ! Initialize it to pointing to nothing.
1781 ALLOCATE ( current )
1782 NULLIFY ( current%next )
1784 END DO read_meas
1786 ! The last allocated measurement is not used (no matter how loop is exited)
1787 ! so deallocate space.
1789 DEALLOCATE ( current )
1791 ! If unable to read in at least one measurement, return error so that
1792 ! entire observation is discarded. If this was bad data, we forced it
1793 ! to skip over the observations without storing any data. That will be
1794 ! handled in the calling routine.
1796 IF ( ( meas_count .LT. 1 ) .AND. &
1797 ( error .EQ. ok ) .AND. &
1798 ( .NOT. bad_data ) ) THEN
1799 nlevels = 0
1800 error = no_data
1801 END IF
1803 ! This is some diagnostic print-out to state the problems encountered.
1804 ! Anything that is not expected issues a message. Any read errors mean
1805 ! to throw away the observation, though if there was no data,
1806 ! there was nothing to toss out anyways. If the error condition is
1807 ! not recognized, the program will stop from this routine.
1809 SELECT CASE ( error )
1811 CASE ( eof_err )
1812 CALL error_handler (proc_name, &
1813 ' Found EOF, expected measurement. Continuing. ', &
1814 TRIM(location_id) // ' ' // TRIM(location_name), &
1815 .FALSE.)
1817 CASE ( read_err )
1818 CALL error_handler (proc_name, &
1819 ' Error in measurement read. ' // &
1820 ' Discarding entire observation and continuing. ', &
1821 TRIM(location_id) // ' ' // TRIM(location_name), &
1822 .FALSE.)
1823 CALL dealloc_meas (surface)
1825 CASE (no_data , ok)
1827 CASE DEFAULT
1829 CALL error_handler (proc_name," Internal error: ","bad error number.",&
1830 .TRUE.)
1832 END SELECT
1834 END SUBROUTINE read_measurements
1836 !\f
1837 ! -----------------------------------------------------------------------
1839 SUBROUTINE dealloc_meas ( head )
1841 ! This deallocates all nodes in a linked list of measurements.
1843 IMPLICIT NONE
1845 TYPE ( measurement ) , POINTER :: head ! head of linked list
1847 TYPE ( measurement ) , POINTER :: previous &
1848 , temp
1849 INTEGER :: status
1851 ! Start at the head, kill everything that is pointed to. After the list is
1852 ! fully deallocated, disassociate the head pointer.
1854 IF ( ASSOCIATED ( head ) ) THEN
1856 previous => head
1857 list_loop : DO WHILE ( ASSOCIATED ( previous%next ) )
1858 temp => previous
1859 previous => previous%next
1860 DEALLOCATE ( temp , STAT = status)
1861 IF (status .NE. 0 ) THEN
1862 WRITE (UNIT = 0, FMT = '(A)') &
1863 'Error in DEALLOCATE, continuing by stopping DEALLOCATE on this list.'
1864 EXIT list_loop
1865 END IF
1866 END DO list_loop
1868 NULLIFY ( head )
1869 END IF
1871 END SUBROUTINE dealloc_meas
1873 !\f
1874 !---------------------------------------------------------------------------
1876 SUBROUTINE sub_info_levels ( surface, levels )
1879 ! This routine takes the sounding and makes sure that if a surface
1880 ! level exists, that it is the first level.
1882 IMPLICIT NONE
1884 TYPE ( measurement ) , POINTER :: surface
1885 INTEGER , INTENT(OUT) :: levels
1887 TYPE ( measurement ) , POINTER :: current
1889 ! Um, is there any data at all?
1891 levels = 0
1893 IF ( ASSOCIATED ( surface ) ) THEN
1895 levels = levels + 1
1897 current => surface%next
1899 DO WHILE ( ASSOCIATED ( current ) )
1901 levels = levels + 1
1902 current => current%next
1904 END DO
1906 END IF
1908 END SUBROUTINE sub_info_levels
1909 !---------------------------------------------------------------------------
1911 SUBROUTINE missing_hp ( surface )
1914 ! This routine takes the sounding and makes sure that if a surface
1915 ! level exists, that it is the first level.
1917 IMPLICIT NONE
1919 TYPE ( measurement ) , POINTER :: surface
1920 TYPE ( measurement ) , POINTER :: current
1922 ! Um, is there any data at all?
1924 IF ( ASSOCIATED ( surface ) ) THEN
1926 ! Alrighty, we have data, so loop through the sounding to see if their is
1927 ! a surface observation. We can't very well have the surface be the first
1928 ! level if we don't have one. Also, start looking at location #2
1929 ! (surface%next) so that we don't "fix" the ones that aren't broken.
1931 IF (surface%meas%height%qc == missing) &
1932 surface%meas%height%data = missing_r
1933 IF (surface%meas%pressure%qc == missing) &
1934 surface%meas%pressure%data = missing_r
1936 current => surface%next
1938 DO WHILE ( ASSOCIATED ( current ) )
1940 IF (current%meas%height%qc == missing) &
1941 current%meas%height%data = missing_r
1942 IF (current%meas%pressure%qc == missing) &
1943 current%meas%pressure%data = missing_r
1945 current => current%next
1947 END DO
1949 END IF
1951 END SUBROUTINE missing_hp
1953 !---------------------------------------------------------------------------
1955 SUBROUTINE surf_first ( surface , elevation )
1958 ! This routine takes the sounding and makes sure that if a surface
1959 ! level exists, that it is the first level.
1961 IMPLICIT NONE
1963 TYPE ( measurement ) , POINTER :: surface
1964 REAL , INTENT(IN) :: elevation
1966 TYPE ( measurement ) , POINTER :: current
1968 ! Um, is there any data at all?
1970 IF ( ASSOCIATED ( surface ) ) THEN
1972 ! Alrighty, we have data, so loop through the sounding to see if their is
1973 ! a surface observation. We can't very well have the surface be the first
1974 ! level if we don't have one. Also, start looking at location #2
1975 ! (surface%next) so that we don't "fix" the ones that aren't broken.
1977 current => surface%next
1979 find_sfc : DO WHILE ( ASSOCIATED ( current ) )
1981 IF ( eps_equal ( current%meas%height%data , elevation , 1. ) ) THEN
1982 surface => current
1983 EXIT find_sfc
1984 END IF
1986 current => current%next
1988 END DO find_sfc
1990 END IF
1992 END SUBROUTINE surf_first
1994 !---------------------------------------------------------------------------
1996 SUBROUTINE insert_at ( surface , new , elevation)
1998 ! This takes a new measurement (new) and inserts it in a linked list
1999 ! of measurements (surface points to first in list) in decreasing order of
2000 ! pressure value. If two levels' pressures are 'eps_equal', the levels
2001 ! are merged instead of being linked.
2003 USE module_obs_merge
2005 IMPLICIT NONE
2007 TYPE ( measurement ) , POINTER :: surface , new
2008 REAL , INTENT(IN) :: elevation
2010 TYPE ( measurement ) , POINTER :: current , previous , oldptr
2011 REAL :: new_pres , new_height
2012 CHARACTER ( LEN = 32 ) , PARAMETER :: name = 'insert_at'
2015 ! INCLUDE 'error.inc'
2016 ! INTERFACE
2017 ! INCLUDE 'error.int'
2018 ! END INTERFACE
2020 ! Initialize the variable to test the pressure and the place where the
2021 ! to-be-inserted measurement points.
2023 new_pres = new%meas%pressure%data
2024 new_height = new%meas%height%data
2026 NULLIFY ( new%next )
2028 ! The first check is to see if we are at the head of the linked list.
2029 ! This drops us through to exit the routine.
2031 IF ( .NOT. ASSOCIATED ( surface ) ) THEN
2033 surface => new
2035 ! We are either between a couple of values, after a last value,
2036 ! or we could need to be merged with a level.
2037 ! All those tests are handled in this else block.
2039 ELSE
2041 ! Initialize some dummy pointers to traverse to where we need to be.
2043 previous => surface
2044 current => surface
2046 ! Loop to find correct location to link in 'new'.
2047 ! The pressure is monotonically decreasing, so as soon as we find one
2048 ! where the current pressure is less than the new pressure,
2049 ! the new pressure goes just before it (or we run out of data looking!).
2050 ! Additionally, if both of the heights are equal AND the heights are
2051 ! the same as the input elevation of the station, then these need to be
2052 ! merged surface observations.
2054 still_some_data : DO WHILE ( ASSOCIATED ( current ) )
2056 IF ( current%meas%pressure%data .LT. new_pres ) EXIT still_some_data
2058 previous => current
2059 current => current%next
2061 END DO still_some_data
2063 ! There are several cases:
2064 ! 1) the new value has the same pressure as the previous value, or
2065 ! both heights are equal to the station elevation: merge them
2066 ! 2) ran out of data finding where to insert level: add it to the end
2067 ! 3) the new value has the same pressure as the current pressure value,
2068 ! or both heights are equal to the station elevation: merge them
2069 ! 4) new pressure is < the previous value: stick it at end of previous
2070 ! 5) new pressure > than previous: put at head of list
2071 ! ***** THE ORDER OF THE TESTS IS IMPORTANT *****
2073 IF ((eps_equal (previous%meas%pressure%data, new_pres , 1. )) .OR. &
2074 ((eps_equal (previous%meas%height%data , new_height , 1. )) .AND. &
2075 (eps_equal (previous%meas%height%data , elevation , 1. )))) THEN
2077 CALL merge_measurements (previous%meas , new%meas , 1)
2078 DEALLOCATE (new)
2080 ELSE IF (.NOT. ASSOCIATED (current)) THEN
2082 previous%next => new
2084 ELSE IF ((eps_equal (current%meas%pressure%data, new_pres , 1.)) .OR. &
2085 ((eps_equal (current%meas%height%data , new_height , 1.)) .AND. &
2086 (eps_equal (current%meas%height%data , elevation , 1.)))) THEN
2088 CALL merge_measurements (current%meas, new%meas , 1)
2090 DEALLOCATE (new)
2092 ELSE IF (previous%meas%pressure%data .GT. new_pres) THEN
2094 oldptr => previous%next
2095 previous%next => new
2096 new%next => oldptr
2098 ELSE IF (previous%meas%pressure%data .LT. new_pres) THEN
2100 ! If we aren't at head of list, have some internal (fatal) error.
2102 IF (.NOT. ASSOCIATED (previous, surface)) THEN
2103 CALL error_handler (name, 'Logic error in IF' ,"", .TRUE.)
2104 ELSE
2105 oldptr => surface
2106 surface => new
2107 new%next => oldptr
2108 END IF
2110 ELSE
2112 ! One of those "should never get here" logic errors, fatal.
2114 CALL error_handler (name, "Logic error in IF test: ",&
2115 "for where to put the new observation level.", .TRUE.)
2117 END IF
2119 END IF
2121 END SUBROUTINE insert_at
2123 !\f
2124 ! -------------------------------------------------------------------------
2126 SUBROUTINE output_obs ( obs , unit , file_name , num_obs , out_opt, forinput )
2128 ! Take the array of observations and write them including measurements
2129 ! at all levels. The two options (out_opt and forinput) are described
2130 ! below.
2132 ! If ( out_opt is 0 ) , write everything
2133 ! > 0 , write only non-discard data
2134 ! < 0 , write only discarded data
2136 ! If ( forinput is true ) output can be pipe back for input.
2138 IMPLICIT NONE
2140 TYPE ( report ) , INTENT ( IN ) , DIMENSION ( : ) :: obs
2141 INTEGER , INTENT ( IN ) :: num_obs
2142 INTEGER , INTENT ( IN ) :: out_opt
2143 INTEGER , INTENT ( IN ) :: unit
2144 CHARACTER ( LEN = * ) , INTENT ( IN ) :: file_name
2145 LOGICAL , INTENT ( IN ) :: forinput
2147 INTEGER :: i , iout
2148 TYPE ( measurement ) , POINTER :: next
2149 TYPE ( meas_data ) :: end_meas
2151 end_meas%pressure%data = end_data_r
2152 end_meas%height%data = end_data_r
2153 end_meas%temperature%data = end_data_r
2154 end_meas%dew_point%data = end_data_r
2155 end_meas%speed%data = end_data_r
2156 end_meas%direction%data = end_data_r
2157 end_meas%u%data = end_data_r
2158 end_meas%v%data = end_data_r
2159 end_meas%rh%data = end_data_r
2160 end_meas%thickness%data = end_data_r
2161 end_meas%pressure%qc = end_data
2162 end_meas%height%qc = end_data
2163 end_meas%temperature%qc = end_data
2164 end_meas%dew_point%qc = end_data
2165 end_meas%speed%qc = end_data
2166 end_meas%direction%qc = end_data
2167 end_meas%u%qc = end_data
2168 end_meas%v%qc = end_data
2169 end_meas%rh%qc = end_data
2170 end_meas%thickness%qc = end_data
2172 OPEN ( UNIT = unit , FILE = file_name , ACTION = 'write' , FORM = 'formatted' )
2174 iout = 0
2176 DO i = 1 , num_obs
2178 IF ( out_opt .EQ. 0 .OR. &
2179 ( out_opt .GT. 0 .AND. .NOT. obs(i)%info%discard ) .OR. &
2180 ( out_opt .LT. 0 .AND. obs(i)%info%discard ) ) THEN
2182 iout = iout + 1
2183 IF ( .NOT. forinput ) write(unit,*) '**************** Next Observation *******************'
2184 WRITE ( UNIT = unit , FMT = rpt_format ) &
2185 obs(i)%location % latitude, obs(i)%location % longitude, &
2186 obs(i)%location % id, obs(i)%location % name, &
2187 obs(i)%info % platform, obs(i)%info % source, &
2188 obs(i)%info % elevation, obs(i)%info % num_vld_fld, &
2189 obs(i)%info % num_error, obs(i)%info % num_warning, &
2190 obs(i)%info % seq_num, obs(i)%info % num_dups, &
2191 obs(i)%info % is_sound, obs(i)%info % bogus, &
2192 obs(i)%info % discard, &
2193 obs(i)%valid_time % sut, obs(i)%valid_time % julian, &
2194 obs(i)%valid_time % date_char, &
2195 obs(i)%ground%slp%data, obs(i)%ground%slp%qc,&
2196 obs(i)%ground%ref_pres%data, obs(i)%ground%ref_pres%qc,&
2197 obs(i)%ground%ground_t%data, obs(i)%ground%ground_t%qc,&
2198 obs(i)%ground%sst%data, obs(i)%ground%sst%qc,&
2199 obs(i)%ground%psfc%data, obs(i)%ground%psfc%qc,&
2200 obs(i)%ground%precip%data, obs(i)%ground%precip%qc,&
2201 obs(i)%ground%t_max%data, obs(i)%ground%t_max%qc,&
2202 obs(i)%ground%t_min%data, obs(i)%ground%t_min%qc,&
2203 obs(i)%ground%t_min_night%data, obs(i)%ground%t_min_night%qc,&
2204 obs(i)%ground%p_tend03%data, obs(i)%ground%p_tend03%qc,&
2205 obs(i)%ground%p_tend24%data, obs(i)%ground%p_tend24%qc, &
2206 obs(i)%ground%cloud_cvr%data, obs(i)%ground%cloud_cvr%qc, &
2207 obs(i)%ground%ceiling%data, obs(i)%ground%ceiling%qc
2208 ! obs(i)%location,
2209 ! obs(i)%valid_time,
2210 ! obs(i)%ground,
2212 next => obs(i)%surface
2213 DO WHILE ( ASSOCIATED ( next ) )
2214 if ( obs(i)%info%discard ) exit
2215 ! WRITE ( UNIT = unit , FMT = meas_format ) next%meas
2216 WRITE ( UNIT = unit , FMT = meas_format ) &
2217 next%meas % pressure % data, next%meas % pressure % qc, &
2218 next%meas % height % data, next%meas % height % qc, &
2219 next%meas % temperature % data, next%meas % temperature % qc, &
2220 next%meas % dew_point % data, next%meas % dew_point % qc, &
2221 next%meas % speed % data, next%meas % speed % qc, &
2222 next%meas % direction % data, next%meas % direction % qc, &
2223 next%meas % u % data, next%meas % u % qc, &
2224 next%meas % v % data, next%meas % v % qc, &
2225 next%meas % rh % data, next%meas % rh % qc, &
2226 next%meas % thickness % data, next%meas % thickness % qc
2228 next => next%next
2229 END DO
2230 ! WRITE ( UNIT = unit , FMT = meas_format ) end_meas
2231 WRITE ( UNIT = unit , FMT = meas_format ) &
2232 end_meas % pressure % data, end_meas % pressure % qc, &
2233 end_meas % height % data, end_meas % height % qc, &
2234 end_meas % temperature % data, end_meas % temperature % qc, &
2235 end_meas % dew_point % data, end_meas % dew_point % qc, &
2236 end_meas % speed % data, end_meas % speed % qc, &
2237 end_meas % direction % data, end_meas % direction % qc, &
2238 end_meas % u % data, end_meas % u % qc, &
2239 end_meas % v % data, end_meas % v % qc, &
2240 end_meas % rh % data, end_meas % rh % qc, &
2241 end_meas % thickness % data, end_meas % thickness % qc
2243 WRITE ( UNIT = unit , FMT = end_format ) obs(i)%info%num_vld_fld, &
2244 obs(i)%info%num_error, obs(i)%info%num_warning
2245 IF ( .NOT. forinput ) &
2246 write(unit,*) 'End of measurements for observation ' , i
2248 END IF
2250 END DO
2252 IF ( .NOT. forinput ) THEN
2253 write(unit,*) '======================================================='
2254 write(unit,*) 'Total Number of Measurements output ' , iout
2255 ENDIF
2257 ! This routine may be called again, with the same unit number, so CLOSE
2258 ! up the file so everything is handled cleanly.
2260 CLOSE ( unit )
2262 END SUBROUTINE output_obs
2264 Subroutine Aircraft_pressure(hh, pp)
2266 ! ---------------------------------------------------------------
2267 ! NOW COMPUTE THE PRESSURE OF THE FLIGHT LEVEL USING THE STANDARD
2268 ! ATMOSPHERE COMPUTATION, SET THE HEIGHT TO MISSING,
2269 ! THEN GROSS CHECK FOR METEOROLOGICAL LIMITS.
2270 ! ---------------------------------------------------------------
2272 implicit none
2274 Type (field), intent(inout) :: hh
2275 Type (field), intent(out) :: pp
2277 IF (HH%data > 11000.0) THEN
2279 PP%data = 226.3 * EXP(1.576106E-4 * (11000.00 - HH%data))
2281 ELSE IF (HH%data <= 11000.0) THEN
2283 PP%data = 1013.25 * (((288.15 - (0.0065 * HH%data))/288.15)**5.256)
2285 END IF
2286 PP%data = PP%data * 100.
2287 PP%qc = 0
2289 !hcl-note: disagree with throwing away the original flight level information
2290 !hcl HH%data = missing_r
2291 !hcl HH%qc = missing
2293 end Subroutine Aircraft_pressure
2295 function psfc_from_QNH(alt, elev) result(psfc)
2296 real, intent(in) :: alt ! altimeter setting/QNH
2297 real, intent(in) :: elev ! elevation
2298 real :: psfc
2299 psfc = (alt**0.190284-(((1013.25**0.190284)*0.0065/288.15)*elev))**5.2553026
2300 end function psfc_from_QNH
2302 END MODULE module_decoded