var/obsproc/src/module_gpspw_caa.F90

   1 MODULE module_gpspw_caa
   2
   3 !------------------------------------------------------------------------------!
   4 ! Read in integrated precipitable water measurement from GPS stations coming
   5 ! from GST (gem file).
   6 !
   7 !  Y.-H. GUO, March 2001
   8 !------------------------------------------------------------------------------!
   9
  10 USE module_date
  11 USE module_type
  12 USE module_inside
  13 USE module_per_type
  14
  15 INCLUDE 'missing.inc'
  16
  17 !------------------------------------------------------------------------------!
  18
  19     TYPE station_name
  20          CHARACTER (LEN=10) :: full_name
  21          CHARACTER (LEN= 4) :: abbreviate_name
  22     END  TYPE station_name
  23
  24 !   Number of stations in thre file
  25
  26     INTEGER, PARAMETER      :: ns = 10
  27
  28 !   Name of station assigned as obs ID (5 digits)
  29
  30     TYPE (station_name), PARAMETER, DIMENSION (ns) :: station = &
  31                                   (/station_name ('BANCHAU   ','banc'), &
  32                                     station_name ('CHIAYI    ','chia'), &
  33                                     station_name ('CHENG-KUNG','chen'), &
  34                                     station_name ('HENGCHU   ','henc'), &
  35                                     station_name ('HSINCHU   ','hsin'), &
  36                                     station_name ('HUALIEN   ','hual'), &
  37                                     station_name ('ILAN      ','ilan'), &
  38                                     station_name ('PANGHU    ','pang'), &
  39                                     station_name ('SUNYI     ','sunm'), &
  40                                     station_name ('SUAO      ','suao')/)
  41
  42
  43 CONTAINS
  44 !------------------------------------------------------------------------------!
  45
  46 SUBROUTINE read_gpspw_caa (filename, filenum, obs, n_obs,                  &
  47                            total_number_of_obs, fatal_if_exceed_max_obs,   &
  48                            time_window_min, time_analysis, time_window_max,&
  49                            ins, jew, missing_flag, &
  50                            print_gpspw_read)
  51 !------------------------------------------------------------------------------!
  52
  53       IMPLICIT NONE
  54
  55       CHARACTER (LEN=*),            INTENT (in)  :: filename
  56       INTEGER,                      INTENT (in)  :: filenum
  57       INTEGER,                      INTENT (inout) :: n_obs
  58       INTEGER,                      INTENT (in)  :: total_number_of_obs
  59       LOGICAL,                      INTENT (in)  :: fatal_if_exceed_max_obs
  60       TYPE (report), DIMENSION (:), INTENT (out) :: obs
  61       INTEGER,                      INTENT (in)  :: ins, jew
  62       CHARACTER (LEN = 19),         INTENT (in)  :: time_window_min
  63       CHARACTER (LEN = 19),         INTENT (in)  :: time_analysis
  64       CHARACTER (LEN = 19),         INTENT (in)  :: time_window_max
  65       REAL,                         INTENT (out) :: missing_flag
  66       LOGICAL,                      INTENT (in)  :: print_gpspw_read
  67
  68       INTEGER                                 :: nyear
  69       INTEGER                                 :: io_error, n, ii
  70       INTEGER                                 :: length, sites_n
  71       INTEGER                                 :: time_n, juln_day, seconds
  72       INTEGER                                 :: obs_num, num_empty, num_outside
  73       CHARACTER (LEN=110)                     :: chr
  74       CHARACTER (LEN= 28)                     :: st_info
  75       REAL                                    :: lat_d,lat_m,lat_s, &
  76                                                  lon_d,lon_m,lon_s
  77       LOGICAL                                 :: outside_domain, outside_window
  78       LOGICAL                                 :: outside
  79
  80       REAL,                DIMENSION (ns)     :: height, lat, lon
  81       CHARACTER (LEN=10), DIMENSION (ns)      :: st_name
  82       INTEGER,             DIMENSION (300)    :: nyr, nmo, ndy, nhr, nmn, nsc
  83       REAL,                DIMENSION (300)    :: time
  84       REAL,                DIMENSION (ns,300) :: pw
  85
  86       CHARACTER (LEN = 80)                    :: file_name
  87       CHARACTER (LEN = 32), PARAMETER         :: proc_name = 'read_gpspw_caa '
  88       INTEGER                                 :: iunit
  89
  90
  91 !------------------------------------------------------------------------------!
  92
  93              WRITE (UNIT = 0, FMT = '(A)')  &
  94 '------------------------------------------------------------------------------'
  95
  96       WRITE (UNIT = 0, FMT = '(A,A,/)') &
  97      'READ GPS OBSERVATIONS IN FILE: ', TRIM (filename)
  98
  99 ! 0.  OPEN DIGANOSTIC FILE
 100 ! ========================
 101
 102       IF (print_gpspw_read) THEN
 103
 104           file_name = 'obs_gpspw_read.diag'
 105           iunit    = 999
 106
 107           OPEN (UNIT = iunit , FILE = file_name , FORM = 'FORMATTED'  , &
 108                 ACTION = 'WRITE' , STATUS = 'REPLACE', IOSTAT = io_error )
 109
 110           IF (io_error .NE. 0 ) THEN
 111               CALL error_handler (proc_name, &
 112              'Unable to open output diagnostic file. ' , file_name, .TRUE.)
 113           ELSE
 114               WRITE (UNIT = 0, FMT = '(A,A,/)') &
 115              "Diagnostics in file ", TRIM (file_name)
 116          ENDIF
 117
 118       ENDIF
 119
 120 ! 1.  OPEN INPUT FILE AND RESET
 121 ! =============================
 122
 123
 124       OPEN (UNIT = filenum, FILE = filename, FORM='formatted', ACTION = 'read',&
 125             IOSTAT = io_error )
 126
 127       REWIND (filenum)
 128
 129       !  Year of data
 130
 131       READ (time_analysis,'(I4)') nyear
 132
 133       pw       =  999.990
 134       time     =  999.990
 135       time_n   = -9999
 136       sites_n  = 0
 137       io_error = 0
 138
 139       obs_num = n_obs + 1
 140
 141 read: DO WHILE ( io_error >= 0 )
 142
 143       READ (filenum,'(A)', iostat = io_error) CHR
 144
 145 ! 1.  READ STATION INFORMATION:
 146 ! =============================
 147
 148       !  Get the station name:
 149
 150       IF (LEN_TRIM (chr) < 15 .AND. LEN_TRIM (chr) > 0) THEN
 151
 152           !  Search the index based on the station name:
 153
 154           ii = 0
 155
 156 search:   DO n = 1,ns
 157
 158              IF (chr(4:13) == station(n) % full_name) THEN
 159                  ii = n
 160                  EXIT search
 161              ENDIF
 162
 163           ENDDO search
 164
 165           !  If there is no data for this station, go to read next:
 166
 167           IF (ii > ns .OR. ii < 1) CYCLE read
 168
 169               !  Keep the station name, and counter + 1:
 170
 171               st_name(ii) = chr(4:13)
 172               sites_n     = sites_n + 1
 173
 174           ELSE IF ((ii>0 .AND. ii<=ns) .AND. time_n == -9999) THEN
 175
 176               !  Get the height, latitude, and longitude for the station:
 177
 178               st_info = chr(21:48)
 179
 180 stn_info:     SELECT CASE (st_info(1:2))
 181
 182               CASE ('HE')  stn_info
 183
 184                            READ (st_info,'(12X,F12.4)') height(ii)
 185                            if (height(ii) < 0.0) height(ii) = missing_r
 186
 187               CASE ('LA')  stn_info
 188
 189                            READ (st_info,'(12X,2F3.0,G10.6)') lat_d,lat_m,lat_s
 190                            lat (ii) = lat_d + lat_m/60. + lat_s/3600.
 191
 192               CASE ('LO')  stn_info
 193
 194                            read(st_info,'(12X,2F3.0,G10.6)') lon_d,lon_m,lon_s
 195                            lon(ii) = lon_d + lon_m/60. + lon_s/3600.
 196
 197               CASE DEFAULT  stn_info
 198
 199                       IF (print_gpspw_read) THEN
 200                       WRITE (UNIT = iunit, FMT = &
 201                       &    '(''II='',I2,2X,''site_n='',I2,2X,''STN_INFO:'',&
 202                       &    A,'' HEIGHT, LAT, and LON. ARE MISSING'',/)')  &
 203                       &    ii,sites_n,TRIM (st_info)
 204                       ENDIF
 205
 206                       ii = 0
 207
 208               END SELECT  stn_info
 209
 210           ENDIF
 211
 212
 213
 214 ! 2.  READ PW DATA:
 215 ! =================
 216
 217       IF (chr(1:4) == '====') THEN
 218
 219            time_n = 0
 220
 221       ELSE IF (time_n >= 0 .and. io_error >= 0) THEN
 222
 223            time_n = time_n + 1
 224
 225            READ (chr,'(11F9.3)') time(time_n),(pw(n,time_n),n=1,ns)
 226
 227            nyr (time_n) = nyear
 228            juln_day     = INT (time(time_n))
 229            seconds      = (time(time_n) - Juln_day) * 86400.
 230
 231            CALL julian_day    (nyear, nmo (time_n), ndy (time_n), juln_day, 2)
 232            CALL sec_to_hhmmss (nhr (time_n), nmn (time_n), nsc (time_n), &
 233                                seconds, 2)
 234
 235       ENDIF
 236
 237       ENDDO read
 238
 239       num_empty   = 0
 240       num_outside = 0
 241
 242 n_time:&
 243       DO ii = 1, time_n
 244
 245 stns:    DO n = 1,ns
 246
 247             !  Fill GPS PW data structure
 248
 249             obs(obs_num) % location % id        = adjustl (station(n)%abbreviate_name)
 250             obs(obs_num) % location % name      = adjustl (station(n)%full_name)
 251             obs(obs_num) % location % latitude  = lat(n)
 252             obs(obs_num) % location % longitude = lon(n)
 253
 254             obs(obs_num) % info % platform      = 'FM-111 GPSPW'
 255             obs(obs_num) % info % source        = 'TWP FROM TAIWAN'
 256             obs(obs_num) % info % elevation     =  height(n)
 257             obs(obs_num) % info % discard       = .false.
 258             obs(obs_num) % info % is_sound      = .false.
 259
 260             obs(obs_num) % valid_time % julian  = time(ii)
 261
 262             WRITE (obs(obs_num) % valid_time % date_char,'(I4,5I2.2)') &
 263                    nyr (ii), nmo (ii), ndy (ii), nhr (ii), nmn (ii), nsc (ii)
 264
 265             WRITE (obs(obs_num) % valid_time % date_mm5,'(I4,5('':'',I2.2))') &
 266                    nyr (ii), nmo (ii), ndy (ii), nhr (ii), nmn (ii), nsc (ii)
 267
 268             IF (pw(n,ii) == 0.000 .or. pw(n,ii) ==  999.990 ) then
 269                obs(obs_num) % ground % pw % qc   =  missing
 270                obs(obs_num) % ground % pw % data =  missing_r
 271                obs(obs_num) % ground % pw % error=  missing_r
 272             ELSE
 273                obs(obs_num) % ground % pw % data =  pw(n,ii)
 274                obs(obs_num) % ground % pw % qc   =  0
 275                obs(obs_num) % ground % pw % error=  0.2
 276             ENDIF
 277
 278             !  Reset slp data structure
 279
 280             obs(obs_num) % ground % slp % qc     =  missing
 281             obs(obs_num) % ground % slp % data   =  missing_r
 282             obs(obs_num) % ground % slp % error  =  missing_r
 283
 284             !  Check domain and time window
 285
 286             CALL inside_domain (obs(obs_num)%location%latitude  , &
 287                                 obs(obs_num)%location%longitude , &
 288                                 ins , jew , outside_domain , &
 289                                 obs(obs_num)%location%xjc, &
 290                                 obs(obs_num)%location%yic, &
 291                                 obs(obs_num)%location%xjd, &
 292                                 obs(obs_num)%location%yid)
 293
 294             CALL inside_window (obs(obs_num)%valid_time%date_char, &
 295                                 time_window_min, time_window_max, &
 296                                 outside_window, iunit)
 297
 298             outside = outside_domain .OR. outside_window
 299
 300
 301             IF (.NOT.outside) THEN
 302
 303                 !  Unlink upper levels
 304
 305                  NULLIFY (obs (obs_num) % surface)
 306
 307                  !  Print station information:
 308
 309                  IF (print_gpspw_read) THEN
 310
 311                      WRITE (UNIT = iunit, FMT = '(A,A5,1X,A23,2F9.3)', &
 312                             ADVANCE = 'no' )&
 313                            'Found Name and ID = ' , &
 314                             obs (obs_num) % location % id,          &
 315                             obs (obs_num) % location % name,        &
 316                             obs (obs_num) % location % latitude,    &
 317                             obs (obs_num) % location % longitude
 318                             WRITE (UNIT = iunit , FMT = '(A)' ) ' => GPS PW'
 319
 320                  ENDIF
 321
 322                  !  Increment the number of read station
 323
 324                  obs_num = obs_num + 1
 325
 326                  !  Check if there's enough memory to read a new data
 327
 328                  IF  ((obs_num .GT. total_number_of_obs) .AND.  &
 329                       (fatal_if_exceed_max_obs)) THEN
 330
 331                  CALL error_handler (proc_name, &
 332                 "Too many obs increase the max_number_of_obs_nml in namelist",&
 333                  "",.true.)
 334
 335                  ELSE IF ((obs_num .GT. total_number_of_obs) .AND. &
 336                           (.NOT. fatal_if_exceed_max_obs))    THEN
 337
 338                  CALL error_handler (proc_name, &
 339                 "Too many obs increase the max_number_of_obs_nml in namelist.",&
 340                  .false.)
 341
 342                  EXIT  n_time
 343
 344                  ENDIF
 345
 346             ELSE
 347                  num_outside = num_outside + 1
 348             ENDIF
 349
 350         ENDDO stns
 351
 352       ENDDO n_time
 353
 354       WRITE (UNIT = 0, FMT = '(A)') 'Have reached end of observations file.'
 355
 356       CLOSE (filenum)
 357
 358       obs_num  = obs_num - 1
 359
 360       !  Total number of observation accumulated per type
 361
 362       ngpspws (icor+0) = obs_num+num_empty+num_outside-n_obs
 363       ngpspws (icor+1) = num_empty
 364       ngpspws (icor+2) = num_outside
 365
 366       !  Print number of reports
 367
 368       WRITE (UNIT = 0, FMT = '(/,A)')  &
 369 '------------------------------------------------------------------------------'
 370
 371       WRITE (UNIT = 0, FMT = '(A)') 'GPS PW OBSERVATIONS READ:'
 372
 373       WRITE (UNIT = 0,FMT = '(/,A,I6)') &
 374      'Number of GPSPW reports: ',ngpspws (icor)
 375
 376       WRITE (UNIT = 0 , FMT = '(/,4(A,I8,/))' ) &
 377
 378           "Number of observations read:          ",obs_num+    &
 379                                                    num_empty+  &
 380                                                    num_outside-&
 381                                                    n_obs,      &
 382           "Number of empty observations:         ",num_empty,  &
 383           "Number of out of domain observations: ",num_outside,&
 384           "Number of observations for ingestion: ",obs_num-n_obs
 385
 386
 387       !  Total number of observation accumulated
 388
 389       n_obs = obs_num
 390
 391
 392       !  Missing data flag
 393
 394       missing_flag = missing_r
 395
 396       IF (print_gpspw_read) CLOSE (iunit)
 397
 398     END SUBROUTINE read_gpspw_caa
 399
 400 END MODULE module_gpspw_caa