PROGRAM main_obsproc !-----------------------------------------------------------------------------! ! Read decoded observations (output format of the MM5 "gts_decoder" and "fetch" ! facilities, which is also the input format ! of MM5 "rawins" and "little_r" programs). ! Keep only records with at least height or pressure ! in MM5 horizontal domain as defined in namelist ! Fill with pressure when missing, ! Sort observations by location and time, ! Merge space duplicate stations (same type, same location, same time), ! Remove time duplicate stations in time (same type & location, different ! time), ! Fill with height when missing, ! Check soundings (vertical consistency and super adiabatic), ! Remove observations above MM5 lid as defined by ptop in namelist, ! Estimate observational error, ! Write out for inputting into MM5 3D-VAR ! !-----------------------------------------------------------------------------! ! ! HISTORY: ! ! F. VANDENBERGHE, March 2001 ! ! 01/13/2003 - Updated for Profiler obs. S. R. H. Rizvi ! ! 02/04/2003 - Updated for Buoy obs. S. R. H. Rizvi ! ! 02/11/2003 - Reviewed and modified for Profiler ! and Buoy obs. Y.-R. Guo ! ! 05/23/2003 - GPS ZTD code added L. Cucurull ! 08/13/2003 - Reviewed Y.-R. Guo ! 09/XX/2005 - Modified to output in PREPBUFR format J. Drake ! 06/30/2006 - Updated for AIRS retrievals Syed RH Rizvi ! 10/09/2006 - Updated for GPS Excess Phase S.-Y. Chen !------------------------------------------------------------------------------ USE module_mm5 USE module_map USE map_utils USE module_namelist USE module_decoded USE module_type USE module_func USE module_icao USE module_sort USE module_duplicate USE module_per_type USE module_recoverp USE module_recoverh USE module_diagnostics USE module_qc USE module_err_ncep USE module_err_afwa USE module_complete USE module_thin_ob USE module_write !------------------------------------------------------------------------------! IMPLICIT NONE CHARACTER (LEN = 80) :: nml_filename CHARACTER (LEN = 80) :: title, caption LOGICAL :: exist INTEGER :: iunit, ii, fm_code, ntime INTEGER :: ins, jew INTEGER :: loop_index, number_of_obs INTEGER :: total_dups_loc, total_dups_time INTEGER :: total_dups INTEGER :: map_projection REAL :: missing_flag TYPE (report), DIMENSION (:), ALLOCATABLE :: obs, obs_copy TYPE (measurement), POINTER :: next, current INTEGER, DIMENSION (:), ALLOCATABLE :: index real :: lat, lon, xjd, yid, xxi, yyj, xxi1, yyj1 CHARACTER (LEN = 19) :: time_min CHARACTER (LEN = 19) :: time_fg CHARACTER (LEN = 19) :: time_max CHARACTER (LEN = 31) :: SLOT_TITLE LOGICAL :: dis0, dis1, dis2 !------------------------------------------------------------------------------! ! 1. READ THE NAMELIST ! ==================== iunit = 99 nml_filename = 'namelist.obsproc' WRITE (UNIT = 0, FMT = '(/,A,A,/,A)') & ' READ NAMELIST FILE: ',TRIM (nml_filename), & ' ------------------' CALL get_namelist (nml_filename, iunit) ! 2. SET THE MM5 MAP PROJECTION COEFICIENTS ! ========================================== ! 2.1 Height at model lid ! ------------------- htop = 0.; h_tropo = 0. ! Tropopause height: h_tropo = Ref_height (pis0) htop = Ref_height (ptop) height_max_icao = h_from_p_icao (0.0) ! write(0,'(/5(A,f12.1)/)') & ! " htop=", htop, ", h_tropo=", h_tropo, & ! ", height_max_icao=", height_max_icao, & ! ", ptop=", ptop, ", pis0=", pis0 ! 2.2 Map coeficients ! --------------- CALL setup (domain_check_h, iproj, phic, xlonc, truelat1, truelat2, & #ifdef BKG dis, xcntr, ycntr, xn, pole, psi1, c2) #else maxnes, nestix, nestjx, dis, numc, nesti, nestj, & ixc, jxc, xcntr, ycntr, xn, pole, psi1, c2, & xim11, xjm11) #endif if (fg_format == 'WRF') then lat = map_info%lat1 lon = map_info%lon1 call LLXY(lat, lon, xxi, yyj) write(0,'(/" LLXY:lat, lon, (dot) xxi, yyj:",4f12.5 )') & lat, lon, xxi, yyj call latlon_to_ij(map_info, lat, lon, xxi, yyj) write(0,'( "latlon_to_ij:lat, lon, (cross)xxi, yyj:",4f12.5/)') & lat, lon, xxi, yyj endif ! 2.3 Grid dimensions for the current domain ! -------------------------------------- ins = nestix (idd) jew = nestjx (idd) ! 3. LOAD GTS OBSERVATIONS ! ========================= icor = 0; caption (7*(icor+0)+1:7*(icor+1)) = " READ" caption (7*(icor+1)+1:7*(icor+2)) = " EMPTY" caption (7*(icor+2)+1:7*(icor+3)) = "OUTSIDE" ! 3.1 Reset the number of ingested obs per type ! ----------------------------------------- nsynops = 0; nshipss = 0; nmetars = 0; npilots = 0; nsounds = 0; nsatems = 0; nsatobs = 0; naireps = 0; ngpspws = 0; namdars = 0; nssmt1s = 0; nssmt2s = 0; nssmis = 0; ntovss = 0; nqscats = 0; nothers = 0; nprofls = 0; nbuoyss = 0; ngpsztd = 0; ngpsref = 0; nboguss = 0; nairss = 0; ngpseph = 0; ntamdar = 0; ! 3.2 Reset the total number of ingested obs ! -------------------------------------- number_of_obs = 0 total_dups_time = 0 total_dups_loc = 0 ! 3.3 Allocate memory ! -------------- ALLOCATE (obs (max_number_of_obs)) ALLOCATE (index (max_number_of_obs)) ! 3.4 Read gts observations ! --------------------- INQUIRE (FILE = obs_gts_filename, EXIST = exist ) IF (exist .and. LEN(TRIM(obs_gts_filename))>0) THEN ! Read data from input file iunit = 99 CALL read_obs_gts (obs_gts_filename, iunit, obs, number_of_obs, & max_number_of_obs, fatal_if_exceed_max_obs, print_gts_read, & ins, jew, time_window_min, time_window_max, & map_projection, missing_flag) ! Reset unused memory for subsequent reading DO loop_index = number_of_obs+1, max_number_of_obs NULLIFY (obs (loop_index) % surface) ENDDO ELSE WRITE (0,'(/,A,/)') "No decoded observation file to read." ENDIF ! 3.6 Check if any data have been loaded ! ---------------------------------- IF (number_of_obs .GT. 0) THEN icor = icor + 3 ! 4. OBSERVATION SORTING ! ======================= ! Sort the observations so that they can be merged easily. Really ! what happens is that the index array is uniquely sorted by location ! (except for observations that are from the same "place"). This ! puts duplicate location observations next to each other. ! Then, merge the observations to (try to) remove all duplicates and ! build composite data. ! Then, sort obervations upon time, type and location ! 4.1 Recover the missing pressure based on the observed heights under ! hydrostatic assumption or the model reference state. ! --------------------------------------------------- CALL recover_pressure_from_height (max_number_of_obs , & obs, number_of_obs, print_recoverp) CALL check_obs (max_number_of_obs, obs, number_of_obs, 'pressure') ! 4.2 Sort station per location ! ------------------------- CALL sort_obs (obs ,number_of_obs , compare_loc, index ) ! 4.3 Merge duplicate stations (same location and same time) ! ------------------------------------------------------ ! Because data merging is based on the pressure, the missing ! pressure have been recovered (Sec. 4.1) caption (7*icor+1:7*(icor+1)) = "LOCDUPL" CALL check_duplicate_loc (obs ,index ,number_of_obs, total_dups_loc, & time_analysis, print_duplicate_loc) icor = icor + 1 ! 4.4 Remove duplicate stations (same location and different time) ! ------------------------------------------------------------ caption (7*icor+1:7*(icor+1)) = "TIMDUPL" if (use_for /= '4DVAR' ) & CALL check_duplicate_time (obs ,index ,number_of_obs, total_dups_time, & time_analysis, print_duplicate_time) icor = icor + 1 ! 4.5 Total Number of duplicate stations removed ! ------------------------------------------ total_dups = total_dups_loc + total_dups_time ! 4.7 Sort obs chronologically ! ------------------------ CALL sort_obs (obs ,number_of_obs , compare_tim, index ) ! Temporarily set back index to increasing order ! index = (/ ( loop_index , loop_index = 1 , max_number_of_obs ) /) ! 5. DIAGNOSTICS, ESTIMATE HEIGHT WHEN MISSING, OBSERVATIONAL ERROR ! ================================================================== ! 5.1 Wind components, dew point relative humidity and mixing ratio ! ------------------------------------------------------------- ! Gross error check on wind module (<0), direction (>360) and ! Temperature (>0), dew point (>0&>" CALL obs_err_ncep (max_number_of_obs, obs, number_of_obs) ENDIF ! 6. CHECK COMPLETNESS AND QUALITY CONTROL PURELY BASED ON THE OBS DATA ! ====================================================================== ! 6.1 Vertical consistency check and dry convective adjustment QC ! ----------------------------------------------------------- CALL proc_qc1 (max_number_of_obs , obs , number_of_obs, & qc_test_vert_consistency , qc_test_convective_adj , & print_qc_vert, print_qc_conv) ! 6.2 Vertical domain check ! --------------------- CALL proc_qc2 (max_number_of_obs , obs , number_of_obs, & qc_test_above_lid, print_qc_lid) ! 6.2 Check completness (Remove levels without data and too low or too high) ! --------------------------------------------------------------------- caption (7*icor+1:7*(icor+1)) = "UNCOMPL" CALL check_completness (max_number_of_obs , & obs, number_of_obs, remove_above_lid, & print_uncomplete) icor = icor + 1 ! 6.3 Print statistics on obs ! ----------------------- caption (7*icor+1:7*(icor+1)+1) = "INGESTD" title = "INGESTED OBSERVATION AFTER CHECKS:" CALL print_per_type (TRIM (title), TRIM (caption), icor) icor = icor+1 ! 7. Reduce the QC flags and thin the satellite data ! ================================================== ! 7.1 Reduce QC from 8 to 3 digits ! ---------------------------- CALL qc_reduction (max_number_of_obs, obs, number_of_obs) ! 7.2 Thin the SATOB data before writing out IF (Thining_SATOB) & CALL THIN_OB (max_number_of_obs, obs, number_of_obs, & nestix(idd), nestjx(idd), missing_flag, & 'SATOB', 88, 1000.0) IF (Thining_SSMI) THEN CALL THIN_OB (max_number_of_obs, obs, number_of_obs, & nestix(idd), nestjx(idd), missing_flag, & 'SSMI_Rtvl', 125) CALL THIN_OB (max_number_of_obs, obs, number_of_obs, & nestix(idd), nestjx(idd), missing_flag, & 'SSMI_Tb', 126) ENDIF IF (Thining_QSCAT) & CALL THIN_OB (max_number_of_obs, obs, number_of_obs, & nestix(idd), nestjx(idd), missing_flag, & 'Qscatcat', 281) ! 8. OUTPUT DIRECT OBSERVATIONS (U,V,T,QV) IN ASCII and PREBUFR FILE ! =================================================================== ! 8.1 Keep a copy of the obs ! ------------------------ allocate (obs_copy(1:max_number_of_obs)) obs_copy = obs ! 8.2 Loop over the time slots ! ------------------------ do ntime = 1, num_time_slots ! 8.3 Set the time_min, time_fg, and time_max for the time slot: ntime ! ---------------------------------------------------------------- if ( num_time_slots > 1 ) then ! 8.3.1 the number of time slots greater than 1 for FGAT and 4DVAR ! .......................................................... if ( ntime == 1 .or. ntime == num_time_slots ) then idt = slot_len / 2 - 1 else idt = slot_len - 1 endif if (ntime == 1) then time_min = time_window_min time_fg = time_window_min else call geth_newdate (time_min, time_max, 1) if (ntime == num_time_slots ) then time_fg = time_window_max else call geth_newdate (time_fg, time_min, slot_len/2) endif endif call geth_newdate (time_max, time_min, idt) else ! 8.3.2 the number of time slots equal to 1 for 3DVAR, FGAT and 4DVAR ! .......................................................... time_min = time_window_min call geth_newdate (time_max, time_window_max, -1) time_fg = time_analysis endif ! 8.3.3 Print the time slot information ! ............................... write(0,'(//a,i2,4(2x,a))') "slot=",ntime, & "time_min, time_fg, time_max:", time_min, time_fg, time_max ! 8.4 Get the original total obs back ! ------------------------------- deallocate (obs) allocate(obs(1:max_number_of_obs)) obs = obs_copy ! 8.5 Loop over all the obs ! --------------------- do ii = 1,number_of_obs ! 8.6 determine the obs "discard" flag ! -------------------------------- ! 8.7.1 Initialize the "discard" switch" for a obs ! ........................................... dis0 = obs(ii)%info%discard dis1 = .false. dis2 = .false. ! 8.7.2 Kick out the specific type of obs based on namelist settings ! ............................................................ read(obs(ii)%info%platform(4:6),'(i3)') fm_code if ( (.not.write_synop .and. (fm_code==12 .or. fm_code==14)) .or. & (.not.write_ship .and. fm_code==13) .or. & (.not.write_metar .and. (fm_code==15 .or. fm_code==16)) .or. & (.not.write_buoy .and. (fm_code==18 .or. fm_code==19)) .or. & (.not.write_pilot .and. (fm_code>=32 .and. fm_code<=34)) .or. & (.not.write_sound .and. (fm_code>=35 .and. fm_code<=38)) .or. & (.not.write_amdar .and. fm_code==42) .or. & (.not.write_satem .and. fm_code==86) .or. & (.not.write_satob .and. fm_code==88) .or. & (.not.write_airep .and. (fm_code==96 .or. fm_code==97)) .or. & (.not.write_tamdar.and. fm_code==101) .or. & (.not.write_gpspw .and. fm_code==111) .or. & (.not.write_gpsztd.and. fm_code==114) .or. & (.not.write_gpsref.and. fm_code==116) .or. & (.not.write_gpseph.and. fm_code==118) .or. & (.not.write_ssmt1 .and. fm_code==121) .or. & (.not.write_ssmt2 .and. fm_code==122) .or. & (.not.write_ssmi .and. (fm_code==125 .or. fm_code==126)).or. & (.not.write_tovs .and. fm_code==131) .or. & (.not.write_qscat .and. fm_code==281) .or. & (.not.write_profl .and. fm_code==132) .or. & (.not.write_bogus .and. fm_code==135) .or. & (.not.write_airs .and. fm_code==133) ) & dis1 = .true. ! 8.7.3 Kick out the obs outside the time slot ! ...................................... CALL inside_window (obs(ii)%valid_time%date_char, & time_min, time_max, & dis2) ! 8.7.4 Reset the obs "discard" flag ! ............................ obs(ii)%info%discard = dis0 .or. dis1 .or.dis2 ! 8.8 end of the loop over obs ! ------------------------ enddo ! 8.9 Time duplicate check within a slot for 4DVAR ! -------------------------------------------- if (use_for == '4DVAR' ) then CALL sort_obs (obs ,number_of_obs , compare_loc, index ) CALL check_duplicate_time (obs ,index ,number_of_obs, total_dups_time, & time_fg, print_duplicate_time) endif ! 8.10 Print report per platform type and count the # levels per stations ! ------------------------------------------------------------------ write(SLOT_TITLE,'("OBSERVATIONS FOR OUTPUT SLOT ",I2.2)') ntime CALL sort_platform (max_number_of_obs, obs, number_of_obs, & nsynops (icor), nshipss (icor), nmetars (icor), & npilots (icor), nsounds (icor), nsatems (icor), & nsatobs (icor), naireps (icor), ngpspws (icor), & ngpsztd (icor), ngpsref (icor), ngpseph (icor), & nssmt1s (icor), nssmt2s (icor), nssmis (icor), & ntovss (icor), nothers (icor), namdars (icor), & nqscats (icor), nprofls (icor), nbuoyss (icor), & nboguss (icor), nairss (icor), ntamdar (icor), & SLOT_TITLE) ! 8.11 Determine output type ! ---------------------------------------------- IF (output_ob_format .eq. 1 .or. output_ob_format .eq. 3) THEN ! 8.11.1 Output observations in PREPBUFR ! ------------------------------- CALL output_prep (max_number_of_obs, obs, number_of_obs, index, & prepbufr_table_filename, & prepbufr_output_filename, & nsynops (icor), nshipss (icor), nmetars (icor), & npilots (icor), nsounds (icor), nsatems (icor), & nsatobs (icor), naireps (icor), ngpspws (icor), & ngpsztd (icor), ngpsref (icor), ngpseph (icor), & nssmt1s (icor), nssmt2s (icor), nssmis (icor), & ntovss (icor), nothers (icor), namdars (icor), & nqscats (icor), nprofls (icor), nbuoyss (icor), & nboguss (icor), missing_flag, time_analysis) ENDIF IF (output_ob_format .eq. 2 .or. output_ob_format .eq. 3) THEN ! 8.11.2 Output observations in 3D-VAR Version 3.1 gts format ! ---------------------------------------------------- CALL output_gts_31 (max_number_of_obs, obs, number_of_obs, index, & nsynops (icor), nshipss (icor), nmetars (icor), & npilots (icor), nsounds (icor), nsatems (icor), & nsatobs (icor), naireps (icor), ngpspws (icor), & ngpsztd (icor), ngpsref (icor), ngpseph (icor), & nssmt1s (icor), nssmt2s (icor), nssmis (icor), & ntovss (icor), nothers (icor), namdars (icor), & nqscats (icor), nprofls (icor), nbuoyss (icor), & nboguss (icor), nairss (icor), ntamdar(icor), missing_flag, time_fg) CALL output_ssmi_31 (max_number_of_obs, obs, number_of_obs, index, & nssmis (icor), & missing_flag, time_fg) ENDIF enddo ! 9. END ! ======= ENDIF STOP "99999" END PROGRAM main_obsproc