ungrib/src/ungrib.F

   1 !*****************************************************************************!
   2 ! program ungrib                                                              !
   3 !                                                                             !
   4 ! Questions, comments, suggestions, even complaints should be directed to:    !
   5 !                        wrfhelp@ucar.edu                                     !
   6 ! Externals:                                                                  !
   7 !    Module TABLE                                                             !
   8 !    Module GRIDINFO                                                          !
   9 !    Module FILELIST                                                          !
  10 !    Subroutine READ_NAMELIST                                                 !
  11 !    Subroutine PARSE_TABLE                                                   !
  12 !    Subroutine CLEAR_STORAGE                                                 !
  13 !    Subroutine RD_GRIB                                                       !
  14 !    Subroutine RD_GRIB2                                                      !
  15 !    Subroutine PUT_STORAGE                                                   !
  16 !    Subroutine OUTPUT                                                        !
  17 !    Subroutine C_CLOSE                                                       !
  18 !    Subroutine RRPR                                                          !
  19 !    Subroutine DATINT                                                        !
  20 !    Subroutine FILE_DELETE                                                   !
  21 !                                                                             !
  22 ! Kevin W. Manning, NCAR/MMM  - original 'pregrid' code, 1998-2001            !
  23 ! Jim Bresch, Michael Duda, Dave Gill, NCAR/MMM - adapted for WPS, 2006       !
  24 !                                                                             !
  25 !*****************************************************************************!
  26 !                                                                             !
  27 !*****************************************************************************!
  28 !                                                                             !
  29 ! This program reads GRIB-formatted data and puts it into intermediate format !
  30 ! for passing data to a horizontal interpolation program. The intermediate    !
  31 ! format can be for WPS, SI, or MM5.                                          !
  32 !                                                                             !
  33 ! The program tries to read from files called "GRIBFILE.AAA", "GRIBFILE.AAB", !
  34 ! "GRIBFILE.AAC", ... "GRIBFILE.ABA", "GRIBFILE.ABB", ... "GRIBFILE.ZZZ" until!
  35 ! it cannot find a file.  This naming format allows for up to 17576 files,    !
  36 ! which should be enough for most applications.                               !
  37 !                                                                             !
  38 ! The program searches through those "GRIBFILE.???" files, and pulls out all  !
  39 ! the requested fields which fall between a starting date and an ending date. !
  40 ! It writes the fields from a given time period to a file named according to  !
  41 ! the date and hour, i.e., "FILE:YYYY-MM-DD_HH"                               !
  42 !                                                                             !
  43 !*****************************************************************************!
  44 program ungrib
  45
  46   use table
  47   use gridinfo
  48   use storage_module
  49   use filelist
  50   use datarray
  51   use module_debug
  52   use misc_definitions_module
  53   use stringutil
  54
  55   implicit none
  56
  57   interface
  58    subroutine read_namelist(hstart, hend, delta_time, ntimes,&
  59     ordered_by_date, debug_level, out_format, prefix,    &
  60     add_lvls, new_plvl, interp_type, ec_rec_len)
  61
  62     use misc_definitions_module
  63
  64     character(len=19) :: hstart, hend
  65     integer :: delta_time
  66     integer :: ntimes
  67     logical :: ordered_by_date
  68     integer :: debug_level
  69     character(len=3) :: out_format
  70     character(len=MAX_FILENAME_LEN) :: prefix
  71     logical :: add_lvls
  72     real, dimension(:) :: new_plvl
  73     integer :: interp_type, ec_rec_len
  74    end subroutine read_namelist
  75   end interface
  76
  77   integer :: nunit1 = 12
  78   character(LEN=132) :: gribflnm = 'GRIBFILE.AAA        '    ! won't work with len=12
  79
  80   integer :: debug_level
  81
  82   integer , parameter :: maxlvl = 250
  83
  84   real , dimension(maxlvl) :: plvl, new_plvl
  85   integer :: iplvl
  86   logical :: add_lvls
  87   integer :: interp_type, ec_rec_len
  88
  89   integer :: nlvl
  90
  91   real :: startlat, startlon, deltalat, deltalon
  92   real :: level
  93   character (LEN=9) ::  field
  94   character (LEN=3) ::  out_format
  95   character (LEN=MAX_FILENAME_LEN) ::  prefix
  96
  97   logical :: readit
  98
  99   integer, dimension(255) :: iuarr = 0
 100
 101   character (LEN=19) :: HSTART, HEND, HDATE
 102   character(LEN=19) :: hsave  = '0000-00-00_00:00:00'
 103   integer :: itime
 104   integer :: ntimes
 105   integer :: interval
 106   integer :: ierr
 107   logical :: ordered_by_date
 108   integer :: grib_version
 109   integer :: vtable_columns
 110
 111
 112   call mprintf(.true.,STDOUT,' *** Starting program ungrib.exe *** ')
 113   call mprintf(.true.,LOGFILE,' *** Starting program ungrib.exe *** ')
 114 ! -----------------
 115 ! Read the namelist, and return the information we want:
 116
 117   call read_namelist(hstart, hend, interval, ntimes, &
 118        ordered_by_date, debug_level, out_format, prefix, &
 119        add_lvls, new_plvl, interp_type, ec_rec_len)
 120
 121   call mprintf(.true.,INFORM,"Interval value: %i seconds or  %f hours", &
 122                i1=interval, f1=float(interval)/3600.)
 123
 124   call mprintf(.true.,STDOUT,'Path to intermediate files is %s',s1=get_path(prefix))
 125   call mprintf(.true.,LOGFILE,'Path to intermediate files is %s',s1=get_path(prefix))
 126
 127 ! -----------------
 128 ! Determine GRIB Edition number
 129   grib_version=0
 130 !   NCAR/RDA's EC 113 data are grib1 with a non-standard record length
 131   if ( ec_rec_len .eq. 0 ) then
 132     call edition_num(nunit1, gribflnm, grib_version, ierr)
 133     call mprintf((ierr.eq.3),ERROR,"GRIB file problem")
 134   else
 135     grib_version=1
 136     ierr=0
 137   endif
 138   if (grib_version.eq.2) then
 139      vtable_columns=12
 140 #if defined (USE_PNG) && (USE_JPEG2000)
 141      call mprintf(.true.,INFORM, &
 142         "Linked in png and jpeg libraries for Grib Edition 2")
 143 #else
 144      call mprintf(.true.,STDOUT,"WARNING - Grib Edition 2 data detected, and")
 145      call mprintf(.true.,STDOUT,"        - png and jpeg libs were NOT selected")
 146      call mprintf(.true.,STDOUT,"        - during the build.")
 147      call mprintf(.true.,STDOUT,"Stopping")
 148      call mprintf(.true.,LOGFILE,"WARNING - Grib Edition 2 data detected, and")
 149      call mprintf(.true.,LOGFILE,"        - png and jpeg libs were NOT selected")
 150      call mprintf(.true.,LOGFILE,"        - during the build.")
 151      call mprintf(.true.,LOGFILE,"Stopping")
 152      call mprintf(.true.,ERROR,"NEED_GRIB2_LIBS")
 153 #endif
 154   else
 155      vtable_columns=7
 156   endif
 157   call mprintf(.true.,INFORM,"Reading Grib Edition %i", i1=grib_version)
 158
 159 ! -----------------
 160 ! Read the "Vtable" file, and put the information contained therein into
 161 ! the module "table".
 162
 163   call parse_table(debug_level,vtable_columns)
 164
 165   call mprintf(.true.,DEBUG,"Parsed the vtable.")
 166
 167 ! -----------------
 168 ! Initialize the input filename to GRIBFILE.AA{character just before A}
 169 ! That way, when we update the filename below for the first time, it will
 170 ! have the correct value of GRIBFILE.AAA.
 171
 172   gribflnm(12:12)=char(ichar(gribflnm(12:12))-1)
 173
 174 ! -----------------
 175 ! LOOP2 cycles through the list of files named GRIBFILE.???, until it finds
 176 ! a non-existent file.  Then we exit
 177
 178   LOOP2 : DO
 179
 180      ! At the beginning of LOOP2 update the input filename.
 181      if (gribflnm(12:12).eq.'Z') then
 182         if (gribflnm(11:11).eq.'Z') then
 183           gribflnm(10:10) = char(ichar(gribflnm(10:10))+1)
 184           gribflnm(11:11) = 'A'
 185         else
 186           gribflnm(11:11) = char(ichar(gribflnm(11:11))+1)
 187         endif
 188         gribflnm(12:12) = 'A'
 189      else
 190         gribflnm(12:12) = char(ichar(gribflnm(12:12))+1)
 191      endif
 192
 193      ! Set READIT to .TRUE., meaning that we have not read any records yet
 194      ! from the file GRIBFLNM.
 195
 196      call mprintf(.true.,DEBUG,"Reading from gribflnm %s ",s1=gribflnm)
 197
 198      readit = .TRUE.  ! i.e., "Yes, we want to read a record."
 199
 200
 201 ! LOOP1 reads through the file GRIBFLNM, exiting under two conditions:
 202 !        1) We have hit the end-of-file
 203 !        2) We have read past the ending date HEND.
 204 !
 205 ! Condition 2 assumes that the data in file GRIBFLNM are ordered in time.
 206
 207      LOOP1 : DO
 208         ! At the beginning of LOOP1, we are at a new time period.
 209         ! Clear the storage arrays and associated level information.
 210         nlvl = 0
 211         plvl = -999.
 212         call clear_storage
 213
 214 ! LOOP0 reads through the file GRIBFLNM, looking for data of the current
 215 ! date.  It exits under the following conditions.
 216 !          1) We have hit the end-of-file
 217 !          2) The GRIBFLNM variable has been updated to a nonexistent file.
 218 !          3) We start reading a new date and the data are assumed to be
 219 !             ordered by date.
 220 !          4) We have a valid record and the data are not assumed to be
 221 !             ordered by date.
 222
 223         LOOP0 : DO
 224
 225            ! If we need to read a new grib record, then read one.
 226            if (READIT) then
 227
 228               if (grib_version.ne.2) then
 229                 call mprintf(.true.,DEBUG, &
 230                      "Calling rd_grib1 with iunit %i", i1=nunit1)
 231                 call mprintf(.true.,DEBUG, &
 232                      "flnm = %s",s1=gribflnm)
 233                  ! Read one record at a time from GRIB1 (and older Editions)
 234                  call rd_grib1(nunit1, gribflnm, level, field, &
 235                       hdate, ierr, iuarr, debug_level, ec_rec_len)
 236               else
 237
 238                  ! Read one file of records from GRIB2.
 239                  call mprintf(.true.,DEBUG,"Calling rd_grib2")
 240                  call rd_grib2(nunit1, gribflnm, hdate, &
 241                       grib_version, ierr, debug_level)
 242                  FIELD='NULL'
 243
 244               endif
 245
 246               call mprintf(.true.,DEBUG,"ierr = %i ",i1=ierr)
 247               if (ierr.eq.1) then
 248                  ! We have hit the end of a file.  Exit LOOP0 so we can
 249                  ! write output for date HDATE, and then exit LOOP1
 250                  ! to update the GRIBFLNM
 251                  hsave = hdate
 252                  exit LOOP0
 253               endif
 254
 255               if (ierr.eq.2) then
 256                  ! We have hit the end of all the files.  We can exit LOOP2
 257                  ! because there are no more files to read.
 258                  exit LOOP2
 259               endif
 260
 261               call mprintf(.true.,DEBUG, &
 262                "Read a record %s with date %s", s1=field,s2=hdate(1:13))
 263
 264            endif
 265
 266            call mprintf(.true.,DEBUG, &
 267             "hdate = %s , hsave = %s ",s1=hdate(1:13), s2=hsave(1:13) )
 268
 269 !           if (hdate < hstart) then
 270 !              ! The data read has a date HDATE earlier than the starting
 271 !              ! date HSTART.  So cycle LOOP0 to read the the next GRIB record.
 272 !              READIT = .TRUE.
 273 !              cycle LOOP0
 274 !           endif
 275
 276            if (FIELD.EQ.'NULL') then
 277               ! The data read does not match any fields requested
 278               ! in the Vtable.  So cycle LOOP0 to read the next GRIB record.
 279               READIT = .TRUE.
 280               cycle LOOP0
 281            endif
 282
 283            if (ordered_by_date .and. (hdate > hsave)) then
 284
 285               ! Exit LOOP0, because we started to read data from another
 286               ! date.
 287
 288               call mprintf(.true.,DEBUG, &
 289               "hdate %s > hsave %s so exit LOOP0",s1=hdate,s2=hsave)
 290
 291               ! We set READIT to FALSE because we have not yet processed
 292               ! the data from this record, and we will want to process this
 293               ! data on the next pass of LOOP1 (referring to the next time
 294               ! period of interest.
 295
 296               READIT = .FALSE.
 297
 298               exit LOOP0
 299
 300            endif
 301
 302 ! When we have reached this point, we have a data array ARRAY which has
 303 ! some data we want to save, with field name FIELD at pressure level
 304 ! LEVEL (Pa).  Dimensions of this data are map%nx and map%ny.  Put
 305 ! this data into storage.
 306
 307            if (((field == "SST").or.(field == "SKINTEMP")) .and. &
 308                 (level /= 200100.)) level = 200100.
 309            iplvl = int(level)
 310            call mprintf((.not.allocated(rdatarray)),ERROR, &
 311            "GRIB data slab not allocated in ungrib.F before call to put_storage.")
 312            call put_storage(iplvl,field, &
 313                 reshape(rdatarray(1:map%nx*map%ny),(/map%nx, map%ny/)),&
 314                 map%nx,map%ny)
 315            deallocate(rdatarray)
 316
 317            ! Since we processed the record we just read, we set
 318            ! READIT to .TRUE. so that LOOP0 will read the next record.
 319            READIT = .TRUE.
 320
 321            if (.not. ordered_by_date) then
 322               if (hdate >= hstart) then
 323                  hsave = hdate
 324               endif
 325               exit LOOP0
 326            endif
 327
 328         enddo LOOP0
 329
 330 ! When we have reached this point, we have either hit the end of file, or
 331 ! hit the end of the current date.  Either way, we want to output
 332 ! the data found for this date.  This current date is in HSAVE.
 333 ! However, if (HSAVE == 0000-00-00_00:00:00), no output is necessary,
 334 ! because that 0000 date is the flag for the very opening of a file.
 335
 336         if ((hsave(1:4).ne.'0000').and.(hsave.le.hend)) then
 337            if (debug_level .gt. 100) then
 338              print*, 'Calling output: '//hsave(1:13)
 339              call mprintf(.true.,DEBUG,"Calling output: %s ",s1=hsave(1:13))
 340            endif
 341            call output(hsave, nlvl, maxlvl, plvl, interval, 1, out_format, prefix, debug_level)
 342            hsave=hdate
 343
 344            ! If the next record we process has a date later than HEND,
 345            ! then time to exit LOOP1.
 346            if ((ordered_by_date) .and. (hdate.gt.hend)) exit LOOP1
 347
 348         else
 349            hsave = hdate
 350         endif
 351
 352         ! If we hit the end-of-file, its time to exit LOOP1 so we can
 353         ! increment the GRIBFLNM to read the next file.
 354         if (ierr.eq.1) exit LOOP1
 355
 356      enddo LOOP1
 357
 358 ! When we have reached this point, we read all the data we want to from
 359 ! file GRIBFLNM (either because we reached the end-of-file, or we
 360 ! read past the date HEND).  So we close the file and cycle LOOP2 to read
 361 ! the next file.
 362
 363      if (grib_version.ne.2) then
 364         call c_close(iuarr(nunit1), debug_level, ierr)
 365         iuarr(nunit1) = 0
 366      endif
 367      hsave = '0000-00-00_00:00:00'
 368
 369   enddo LOOP2
 370
 371 ! Now Reread, process, and reoutput.
 372
 373   call mprintf(.true.,INFORM,"First pass done, doing a reprocess")
 374
 375   call rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, &
 376             add_lvls, new_plvl, interp_type, &
 377             debug_level, out_format, prefix)
 378
 379 ! Make sure the filedates are in order, with an inefficient sort:
 380
 381   call sort_filedates
 382
 383 ! Interpolate temporally to fill in missing times:
 384
 385   call datint(filedates, nfiles, hstart, ntimes, interval, out_format, prefix)
 386
 387 ! Now delete the temporary files:
 388
 389   call file_delete(filedates, nfiles, trim(get_path(prefix))//'PFILE:', interval)
 390
 391 ! And Now we are done:
 392
 393    call mprintf(.true.,STDOUT,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
 394    call mprintf(.true.,STDOUT,'!  Successful completion of ungrib.   !')
 395 !  call mprintf(.true.,STDOUT,"!  We're hauling gear at Bandimere.   !")
 396    call mprintf(.true.,STDOUT,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
 397
 398    call mprintf(.true.,LOGFILE,' *** Successful completion of program ungrib.exe *** ')
 399
 400
 401
 402 contains
 403   subroutine sort_filedates
 404     implicit none
 405
 406     integer :: n
 407     logical :: done
 408     if (nfiles > 1) then
 409        done = .FALSE.
 410        do while ( .not. done)
 411           done = .TRUE.
 412           do n = 1, nfiles-1
 413              if (filedates(n) > filedates(n+1)) then
 414                 filedates(size(filedates)) = filedates(n)
 415                 filedates(n) = filedates(n+1)
 416                 filedates(n+1) = filedates(size(filedates))
 417                 filedates(size(filedates)) = '0000-00-00 00:00:00.0000'
 418                 done = .FALSE.
 419              endif
 420           enddo
 421        enddo
 422     endif
 423   end subroutine sort_filedates
 424
 425 end program ungrib