MODULE module_gpspw_caa !------------------------------------------------------------------------------! ! Read in integrated precipitable water measurement from GPS stations coming ! from GST (gem file). ! ! Y.-H. GUO, March 2001 !------------------------------------------------------------------------------! USE module_date USE module_type USE module_inside USE module_per_type INCLUDE 'missing.inc' !------------------------------------------------------------------------------! TYPE station_name CHARACTER (LEN=10) :: full_name CHARACTER (LEN= 4) :: abbreviate_name END TYPE station_name ! Number of stations in thre file INTEGER, PARAMETER :: ns = 10 ! Name of station assigned as obs ID (5 digits) TYPE (station_name), PARAMETER, DIMENSION (ns) :: station = & (/station_name ('BANCHAU ','banc'), & station_name ('CHIAYI ','chia'), & station_name ('CHENG-KUNG','chen'), & station_name ('HENGCHU ','henc'), & station_name ('HSINCHU ','hsin'), & station_name ('HUALIEN ','hual'), & station_name ('ILAN ','ilan'), & station_name ('PANGHU ','pang'), & station_name ('SUNYI ','sunm'), & station_name ('SUAO ','suao')/) CONTAINS !------------------------------------------------------------------------------! SUBROUTINE read_gpspw_caa (filename, filenum, obs, n_obs, & total_number_of_obs, fatal_if_exceed_max_obs, & time_window_min, time_analysis, time_window_max,& ins, jew, missing_flag, & print_gpspw_read) !------------------------------------------------------------------------------! IMPLICIT NONE CHARACTER (LEN=*), INTENT (in) :: filename INTEGER, INTENT (in) :: filenum INTEGER, INTENT (inout) :: n_obs INTEGER, INTENT (in) :: total_number_of_obs LOGICAL, INTENT (in) :: fatal_if_exceed_max_obs TYPE (report), DIMENSION (:), INTENT (out) :: obs INTEGER, INTENT (in) :: ins, jew CHARACTER (LEN = 19), INTENT (in) :: time_window_min CHARACTER (LEN = 19), INTENT (in) :: time_analysis CHARACTER (LEN = 19), INTENT (in) :: time_window_max REAL, INTENT (out) :: missing_flag LOGICAL, INTENT (in) :: print_gpspw_read INTEGER :: nyear INTEGER :: io_error, n, ii INTEGER :: length, sites_n INTEGER :: time_n, juln_day, seconds INTEGER :: obs_num, num_empty, num_outside CHARACTER (LEN=110) :: chr CHARACTER (LEN= 28) :: st_info REAL :: lat_d,lat_m,lat_s, & lon_d,lon_m,lon_s LOGICAL :: outside_domain, outside_window LOGICAL :: outside REAL, DIMENSION (ns) :: height, lat, lon CHARACTER (LEN=10), DIMENSION (ns) :: st_name INTEGER, DIMENSION (300) :: nyr, nmo, ndy, nhr, nmn, nsc REAL, DIMENSION (300) :: time REAL, DIMENSION (ns,300) :: pw CHARACTER (LEN = 80) :: file_name CHARACTER (LEN = 32), PARAMETER :: proc_name = 'read_gpspw_caa ' INTEGER :: iunit !------------------------------------------------------------------------------! WRITE (UNIT = 0, FMT = '(A)') & '------------------------------------------------------------------------------' WRITE (UNIT = 0, FMT = '(A,A,/)') & 'READ GPS OBSERVATIONS IN FILE: ', TRIM (filename) ! 0. OPEN DIGANOSTIC FILE ! ======================== IF (print_gpspw_read) THEN file_name = 'obs_gpspw_read.diag' iunit = 999 OPEN (UNIT = iunit , FILE = file_name , FORM = 'FORMATTED' , & ACTION = 'WRITE' , STATUS = 'REPLACE', IOSTAT = io_error ) IF (io_error .NE. 0 ) THEN CALL error_handler (proc_name, & 'Unable to open output diagnostic file. ' , file_name, .TRUE.) ELSE WRITE (UNIT = 0, FMT = '(A,A,/)') & "Diagnostics in file ", TRIM (file_name) ENDIF ENDIF ! 1. OPEN INPUT FILE AND RESET ! ============================= OPEN (UNIT = filenum, FILE = filename, FORM='formatted', ACTION = 'read',& IOSTAT = io_error ) REWIND (filenum) ! Year of data READ (time_analysis,'(I4)') nyear pw = 999.990 time = 999.990 time_n = -9999 sites_n = 0 io_error = 0 obs_num = n_obs + 1 read: DO WHILE ( io_error >= 0 ) READ (filenum,'(A)', iostat = io_error) CHR ! 1. READ STATION INFORMATION: ! ============================= ! Get the station name: IF (LEN_TRIM (chr) < 15 .AND. LEN_TRIM (chr) > 0) THEN ! Search the index based on the station name: ii = 0 search: DO n = 1,ns IF (chr(4:13) == station(n) % full_name) THEN ii = n EXIT search ENDIF ENDDO search ! If there is no data for this station, go to read next: IF (ii > ns .OR. ii < 1) CYCLE read ! Keep the station name, and counter + 1: st_name(ii) = chr(4:13) sites_n = sites_n + 1 ELSE IF ((ii>0 .AND. ii<=ns) .AND. time_n == -9999) THEN ! Get the height, latitude, and longitude for the station: st_info = chr(21:48) stn_info: SELECT CASE (st_info(1:2)) CASE ('HE') stn_info READ (st_info,'(12X,F12.4)') height(ii) if (height(ii) < 0.0) height(ii) = missing_r CASE ('LA') stn_info READ (st_info,'(12X,2F3.0,G10.6)') lat_d,lat_m,lat_s lat (ii) = lat_d + lat_m/60. + lat_s/3600. CASE ('LO') stn_info read(st_info,'(12X,2F3.0,G10.6)') lon_d,lon_m,lon_s lon(ii) = lon_d + lon_m/60. + lon_s/3600. CASE DEFAULT stn_info IF (print_gpspw_read) THEN WRITE (UNIT = iunit, FMT = & & '(''II='',I2,2X,''site_n='',I2,2X,''STN_INFO:'',& & A,'' HEIGHT, LAT, and LON. ARE MISSING'',/)') & & ii,sites_n,TRIM (st_info) ENDIF ii = 0 END SELECT stn_info ENDIF ! 2. READ PW DATA: ! ================= IF (chr(1:4) == '====') THEN time_n = 0 ELSE IF (time_n >= 0 .and. io_error >= 0) THEN time_n = time_n + 1 READ (chr,'(11F9.3)') time(time_n),(pw(n,time_n),n=1,ns) nyr (time_n) = nyear juln_day = INT (time(time_n)) seconds = (time(time_n) - Juln_day) * 86400. CALL julian_day (nyear, nmo (time_n), ndy (time_n), juln_day, 2) CALL sec_to_hhmmss (nhr (time_n), nmn (time_n), nsc (time_n), & seconds, 2) ENDIF ENDDO read num_empty = 0 num_outside = 0 n_time:& DO ii = 1, time_n stns: DO n = 1,ns ! Fill GPS PW data structure obs(obs_num) % location % id = adjustl (station(n)%abbreviate_name) obs(obs_num) % location % name = adjustl (station(n)%full_name) obs(obs_num) % location % latitude = lat(n) obs(obs_num) % location % longitude = lon(n) obs(obs_num) % info % platform = 'FM-111 GPSPW' obs(obs_num) % info % source = 'TWP FROM TAIWAN' obs(obs_num) % info % elevation = height(n) obs(obs_num) % info % discard = .false. obs(obs_num) % info % is_sound = .false. obs(obs_num) % valid_time % julian = time(ii) WRITE (obs(obs_num) % valid_time % date_char,'(I4,5I2.2)') & nyr (ii), nmo (ii), ndy (ii), nhr (ii), nmn (ii), nsc (ii) WRITE (obs(obs_num) % valid_time % date_mm5,'(I4,5('':'',I2.2))') & nyr (ii), nmo (ii), ndy (ii), nhr (ii), nmn (ii), nsc (ii) IF (pw(n,ii) == 0.000 .or. pw(n,ii) == 999.990 ) then obs(obs_num) % ground % pw % qc = missing obs(obs_num) % ground % pw % data = missing_r obs(obs_num) % ground % pw % error= missing_r ELSE obs(obs_num) % ground % pw % data = pw(n,ii) obs(obs_num) % ground % pw % qc = 0 obs(obs_num) % ground % pw % error= 0.2 ENDIF ! Reset slp data structure obs(obs_num) % ground % slp % qc = missing obs(obs_num) % ground % slp % data = missing_r obs(obs_num) % ground % slp % error = missing_r ! Check domain and time window CALL inside_domain (obs(obs_num)%location%latitude , & obs(obs_num)%location%longitude , & ins , jew , outside_domain , & obs(obs_num)%location%xjc, & obs(obs_num)%location%yic, & obs(obs_num)%location%xjd, & obs(obs_num)%location%yid) CALL inside_window (obs(obs_num)%valid_time%date_char, & time_window_min, time_window_max, & outside_window, iunit) outside = outside_domain .OR. outside_window IF (.NOT.outside) THEN ! Unlink upper levels NULLIFY (obs (obs_num) % surface) ! Print station information: IF (print_gpspw_read) THEN WRITE (UNIT = iunit, FMT = '(A,A5,1X,A23,2F9.3)', & ADVANCE = 'no' )& 'Found Name and ID = ' , & obs (obs_num) % location % id, & obs (obs_num) % location % name, & obs (obs_num) % location % latitude, & obs (obs_num) % location % longitude WRITE (UNIT = iunit , FMT = '(A)' ) ' => GPS PW' ENDIF ! Increment the number of read station obs_num = obs_num + 1 ! Check if there's enough memory to read a new data IF ((obs_num .GT. total_number_of_obs) .AND. & (fatal_if_exceed_max_obs)) THEN CALL error_handler (proc_name, & "Too many obs increase the max_number_of_obs_nml in namelist",& "",.true.) ELSE IF ((obs_num .GT. total_number_of_obs) .AND. & (.NOT. fatal_if_exceed_max_obs)) THEN CALL error_handler (proc_name, & "Too many obs increase the max_number_of_obs_nml in namelist.",& .false.) EXIT n_time ENDIF ELSE num_outside = num_outside + 1 ENDIF ENDDO stns ENDDO n_time WRITE (UNIT = 0, FMT = '(A)') 'Have reached end of observations file.' CLOSE (filenum) obs_num = obs_num - 1 ! Total number of observation accumulated per type ngpspws (icor+0) = obs_num+num_empty+num_outside-n_obs ngpspws (icor+1) = num_empty ngpspws (icor+2) = num_outside ! Print number of reports WRITE (UNIT = 0, FMT = '(/,A)') & '------------------------------------------------------------------------------' WRITE (UNIT = 0, FMT = '(A)') 'GPS PW OBSERVATIONS READ:' WRITE (UNIT = 0,FMT = '(/,A,I6)') & 'Number of GPSPW reports: ',ngpspws (icor) WRITE (UNIT = 0 , FMT = '(/,4(A,I8,/))' ) & "Number of observations read: ",obs_num+ & num_empty+ & num_outside-& n_obs, & "Number of empty observations: ",num_empty, & "Number of out of domain observations: ",num_outside,& "Number of observations for ingestion: ",obs_num-n_obs ! Total number of observation accumulated n_obs = obs_num ! Missing data flag missing_flag = missing_r IF (print_gpspw_read) CLOSE (iunit) END SUBROUTINE read_gpspw_caa END MODULE module_gpspw_caa