WPS/ungrib/src/ungrib.F90

   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   integer :: nunit1 = 12
  58   character(LEN=132) :: gribflnm = 'GRIBFILE.AAA        '    ! won't work with len=12
  59
  60   integer :: debug_level
  61
  62   integer , parameter :: maxlvl = 100
  63
  64   real , dimension(maxlvl) :: plvl
  65   integer :: iplvl
  66
  67   integer :: nlvl
  68
  69   real :: startlat, startlon, deltalat, deltalon
  70   real :: level
  71   character (LEN=9) ::  field
  72   character (LEN=3) ::  out_format
  73   character (LEN=MAX_FILENAME_LEN) ::  prefix
  74
  75   logical :: readit
  76
  77   integer, dimension(255) :: iuarr = 0
  78
  79   character (LEN=19) :: HSTART, HEND, HDATE
  80   character(LEN=19) :: hsave  = '0000-00-00_00:00:00'
  81   integer :: itime
  82   integer :: ntimes
  83   integer :: interval
  84   integer :: ierr
  85   logical :: ordered_by_date
  86   integer :: grib_version
  87   integer :: vtable_columns
  88
  89
  90   call mprintf(.true.,STDOUT,' *** Starting program ungrib.exe *** ')
  91   call mprintf(.true.,LOGFILE,' *** Starting program ungrib.exe *** ')
  92 ! -----------------
  93 ! Read the namelist, and return the information we want:
  94
  95   call read_namelist(hstart, hend, interval, ntimes, &
  96        ordered_by_date, debug_level, out_format, prefix)
  97
  98   call mprintf(.true.,INFORM,"Interval value: %i seconds or  %f hours", &
  99                i1=interval, f1=float(interval)/3600.)
 100
 101   call mprintf(.true.,STDOUT,'Path to intermediate files is %s',s1=get_path(prefix))
 102   call mprintf(.true.,LOGFILE,'Path to intermediate files is %s',s1=get_path(prefix))
 103
 104 ! -----------------
 105 ! Determine GRIB Edition number
 106   grib_version=0
 107   call edition_num(nunit1, gribflnm, grib_version, ierr)
 108   call mprintf((ierr.eq.3),ERROR,"GRIB file problem")
 109   if (grib_version.eq.2) then
 110      vtable_columns=11
 111 #if defined (USE_PNG) && (USE_JPEG2000)
 112      call mprintf(.true.,INFORM, &
 113         "Linked in png and jpeg libraries for Grib Edition 2")
 114 #else
 115      call mprintf(.true.,STDOUT,"WARNING - Grib Edition 2 data detected, and")
 116      call mprintf(.true.,STDOUT,"        - png and jpeg libs were NOT selected")
 117      call mprintf(.true.,STDOUT,"        - during the build.")
 118      call mprintf(.true.,STDOUT,"Stopping")
 119      call mprintf(.true.,LOGFILE,"WARNING - Grib Edition 2 data detected, and")
 120      call mprintf(.true.,LOGFILE,"        - png and jpeg libs were NOT selected")
 121      call mprintf(.true.,LOGFILE,"        - during the build.")
 122      call mprintf(.true.,LOGFILE,"Stopping")
 123      call mprintf(.true.,ERROR,"NEED_GRIB2_LIBS")
 124 #endif
 125   else
 126      vtable_columns=7
 127   endif
 128   call mprintf(.true.,INFORM,"Reading Grib Edition %i", i1=grib_version)
 129
 130 ! -----------------
 131 ! Read the "Vtable" file, and put the information contained therein into
 132 ! the module "table".
 133
 134   call parse_table(debug_level,vtable_columns)
 135
 136   call mprintf(.true.,DEBUG,"Parsed the vtable.")
 137
 138 ! -----------------
 139 ! Initialize the input filename to GRIBFILE.AA{character just before A}
 140 ! That way, when we update the filename below for the first time, it will
 141 ! have the correct value of GRIBFILE.AAA.
 142
 143   gribflnm(12:12)=char(ichar(gribflnm(12:12))-1)
 144
 145 ! -----------------
 146 ! LOOP2 cycles through the list of files named GRIBFILE.???, until it finds
 147 ! a non-existent file.  Then we exit
 148
 149   LOOP2 : DO
 150
 151      ! At the beginning of LOOP2 update the input filename.
 152      if (gribflnm(12:12).eq.'Z') then
 153         if (gribflnm(11:11).eq.'Z') then
 154           gribflnm(10:10) = char(ichar(gribflnm(10:10))+1)
 155           gribflnm(11:11) = 'A'
 156         else
 157           gribflnm(11:11) = char(ichar(gribflnm(11:11))+1)
 158         endif
 159         gribflnm(12:12) = 'A'
 160      else
 161         gribflnm(12:12) = char(ichar(gribflnm(12:12))+1)
 162      endif
 163
 164      ! Set READIT to .TRUE., meaning that we have not read any records yet
 165      ! from the file GRIBFLNM.
 166
 167      call mprintf(.true.,DEBUG,"Reading from gribflnm %s ",s1=gribflnm)
 168
 169      readit = .TRUE.  ! i.e., "Yes, we want to read a record."
 170
 171
 172 ! LOOP1 reads through the file GRIBFLNM, exiting under two conditions:
 173 !        1) We have hit the end-of-file
 174 !        2) We have read past the ending date HEND.
 175 !
 176 ! Condition 2 assumes that the data in file GRIBFLNM are ordered in time.
 177
 178      LOOP1 : DO
 179         ! At the beginning of LOOP1, we are at a new time period.
 180         ! Clear the storage arrays and associated level information.
 181         nlvl = 0
 182         plvl = -999.
 183         call clear_storage
 184
 185 ! LOOP0 reads through the file GRIBFLNM, looking for data of the current
 186 ! date.  It exits under the following conditions.
 187 !          1) We have hit the end-of-file
 188 !          2) The GRIBFLNM variable has been updated to a nonexistent file.
 189 !          3) We start reading a new date and the data are assumed to be
 190 !             ordered by date.
 191 !          4) We have a valid record and the data are not assumed to be
 192 !             ordered by date.
 193
 194         LOOP0 : DO
 195
 196            ! If we need to read a new grib record, then read one.
 197            if (READIT) then
 198
 199               if (grib_version.ne.2) then
 200                 call mprintf(.true.,DEBUG, &
 201                      "Calling rd_grib1 with iunit %i", i1=nunit1)
 202                 call mprintf(.true.,DEBUG, &
 203                      "flnm = %s",s1=gribflnm)
 204                  ! Read one record at a time from GRIB1 (and older Editions)
 205                  call rd_grib1(nunit1, gribflnm, level, field, &
 206                       hdate, ierr, iuarr, debug_level)
 207               else
 208
 209                  ! Read one file of records from GRIB2.
 210                  call mprintf(.true.,DEBUG,"Calling rd_grib2")
 211                  call rd_grib2(nunit1, gribflnm, hdate, &
 212                       grib_version, ierr, debug_level)
 213                  FIELD='NULL'
 214
 215               endif
 216
 217               call mprintf(.true.,DEBUG,"ierr = %i ",i1=ierr)
 218               if (ierr.eq.1) then
 219                  ! We have hit the end of a file.  Exit LOOP0 so we can
 220                  ! write output for date HDATE, and then exit LOOP1
 221                  ! to update the GRIBFLNM
 222                  hsave = hdate
 223                  exit LOOP0
 224               endif
 225
 226               if (ierr.eq.2) then
 227                  ! We have hit the end of all the files.  We can exit LOOP2
 228                  ! because there are no more files to read.
 229                  exit LOOP2
 230               endif
 231
 232               call mprintf(.true.,DEBUG, &
 233                "Read a record %s with date %s", s1=field,s2=hdate(1:13))
 234
 235            endif
 236
 237            call mprintf(.true.,DEBUG, &
 238             "hdate = %s , hsave = %s ",s1=hdate(1:13), s2=hsave(1:13) )
 239
 240 !           if (hdate < hstart) then
 241 !              ! The data read has a date HDATE earlier than the starting
 242 !              ! date HSTART.  So cycle LOOP0 to read the the next GRIB record.
 243 !              READIT = .TRUE.
 244 !              cycle LOOP0
 245 !           endif
 246
 247            if (FIELD.EQ.'NULL') then
 248               ! The data read does not match any fields requested
 249               ! in the Vtable.  So cycle LOOP0 to read the next GRIB record.
 250               READIT = .TRUE.
 251               cycle LOOP0
 252            endif
 253
 254            if (ordered_by_date .and. (hdate > hsave)) then
 255
 256               ! Exit LOOP0, because we started to read data from another
 257               ! date.
 258
 259               call mprintf(.true.,DEBUG, &
 260               "hdate %s > hsave %s so exit LOOP0",s1=hdate,s2=hsave)
 261
 262               ! We set READIT to FALSE because we have not yet processed
 263               ! the data from this record, and we will want to process this
 264               ! data on the next pass of LOOP1 (referring to the next time
 265               ! period of interest.
 266
 267               READIT = .FALSE.
 268
 269               exit LOOP0
 270
 271            endif
 272
 273 ! When we have reached this point, we have a data array ARRAY which has
 274 ! some data we want to save, with field name FIELD at pressure level
 275 ! LEVEL (Pa).  Dimensions of this data are map%nx and map%ny.  Put
 276 ! this data into storage.
 277
 278            if (((field == "SST").or.(field == "SKINTEMP")) .and. &
 279                 (level /= 200100.)) level = 200100.
 280            iplvl = int(level)
 281            call mprintf((.not.allocated(rdatarray)),ERROR, &
 282            "GRIB data slab not allocated in ungrib.F before call to put_storage.")
 283            call put_storage(iplvl,field, &
 284                 reshape(rdatarray(1:map%nx*map%ny),(/map%nx, map%ny/)),&
 285                 map%nx,map%ny)
 286            deallocate(rdatarray)
 287
 288            ! Since we processed the record we just read, we set
 289            ! READIT to .TRUE. so that LOOP0 will read the next record.
 290            READIT = .TRUE.
 291
 292            if (.not. ordered_by_date) then
 293               if (hdate >= hstart) then
 294                  hsave = hdate
 295               endif
 296               exit LOOP0
 297            endif
 298
 299         enddo LOOP0
 300
 301 ! When we have reached this point, we have either hit the end of file, or
 302 ! hit the end of the current date.  Either way, we want to output
 303 ! the data found for this date.  This current date is in HSAVE.
 304 ! However, if (HSAVE == 0000-00-00_00:00:00), no output is necessary,
 305 ! because that 0000 date is the flag for the very opening of a file.
 306
 307         if ((hsave(1:4).ne.'0000').and.(hsave.le.hend)) then
 308            if (debug_level .gt. 100) print*, 'Calling output: '//hsave(1:13)
 309            call output(hsave, nlvl, maxlvl, plvl, interval, 1, out_format, prefix, debug_level)
 310            hsave=hdate
 311
 312            ! If the next record we process has a date later than HEND,
 313            ! then time to exit LOOP1.
 314            if ((ordered_by_date) .and. (hdate.gt.hend)) exit LOOP1
 315
 316         else
 317            hsave = hdate
 318         endif
 319
 320         ! If we hit the end-of-file, its time to exit LOOP1 so we can
 321         ! increment the GRIBFLNM to read the next file.
 322         if (ierr.eq.1) exit LOOP1
 323
 324      enddo LOOP1
 325
 326 ! When we have reached this point, we read all the data we want to from
 327 ! file GRIBFLNM (either because we reached the end-of-file, or we
 328 ! read past the date HEND).  So we close the file and cycle LOOP2 to read
 329 ! the next file.
 330
 331      if (grib_version.ne.2) then
 332         call c_close(iuarr(nunit1), debug_level, ierr)
 333         iuarr(nunit1) = 0
 334      endif
 335      hsave = '0000-00-00_00:00:00'
 336
 337   enddo LOOP2
 338
 339 ! Now Reread, process, and reoutput.
 340
 341   call mprintf(.true.,INFORM,"First pass done, doing a reprocess")
 342
 343   call rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_format, prefix)
 344
 345 ! Make sure the filedates are in order, with an inefficient sort:
 346
 347   call sort_filedates
 348
 349 ! Interpolate temporally to fill in missing times:
 350
 351   call datint(filedates, nfiles, hstart, ntimes, interval, out_format, prefix)
 352
 353 ! Now delete the temporary files:
 354
 355   call file_delete(filedates, nfiles, trim(get_path(prefix))//'PFILE:', interval)
 356
 357 ! And Now we are done:
 358
 359    call mprintf(.true.,STDOUT,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
 360    call mprintf(.true.,STDOUT,'!  Successful completion of ungrib.   !')
 361 !  call mprintf(.true.,STDOUT,"!  We're hauling gear at Bandimere.   !")
 362    call mprintf(.true.,STDOUT,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
 363
 364    call mprintf(.true.,LOGFILE,' *** Successful completion of program ungrib.exe *** ')
 365
 366
 367
 368 contains
 369   subroutine sort_filedates
 370     implicit none
 371
 372     integer :: n
 373     logical :: done
 374     if (nfiles > 1) then
 375        done = .FALSE.
 376        do while ( .not. done)
 377           done = .TRUE.
 378           do n = 1, nfiles-1
 379              if (filedates(n) > filedates(n+1)) then
 380                 filedates(size(filedates)) = filedates(n)
 381                 filedates(n) = filedates(n+1)
 382                 filedates(n+1) = filedates(size(filedates))
 383                 filedates(size(filedates)) = '0000-00-00 00:00:00.0000'
 384                 done = .FALSE.
 385              endif
 386           enddo
 387        enddo
 388     endif
 389   end subroutine sort_filedates
 390
 391 end program ungrib