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