!*****************************************************************************! ! program ungrib ! ! ! ! Questions, comments, suggestions, even complaints should be directed to: ! ! wrfhelp@ucar.edu ! ! Externals: ! ! Module TABLE ! ! Module GRIDINFO ! ! Module FILELIST ! ! Subroutine READ_NAMELIST ! ! Subroutine PARSE_TABLE ! ! Subroutine CLEAR_STORAGE ! ! Subroutine RD_GRIB ! ! Subroutine RD_GRIB2 ! ! Subroutine PUT_STORAGE ! ! Subroutine OUTPUT ! ! Subroutine C_CLOSE ! ! Subroutine RRPR ! ! Subroutine DATINT ! ! Subroutine FILE_DELETE ! ! ! ! Kevin W. Manning, NCAR/MMM - original 'pregrid' code, 1998-2001 ! ! Jim Bresch, Michael Duda, Dave Gill, NCAR/MMM - adapted for WPS, 2006 ! ! ! !*****************************************************************************! ! ! !*****************************************************************************! ! ! ! This program reads GRIB-formatted data and puts it into intermediate format ! ! for passing data to a horizontal interpolation program. The intermediate ! ! format can be for WPS, SI, or MM5. ! ! ! ! The program tries to read from files called "GRIBFILE.AAA", "GRIBFILE.AAB", ! ! "GRIBFILE.AAC", ... "GRIBFILE.ABA", "GRIBFILE.ABB", ... "GRIBFILE.ZZZ" until! ! it cannot find a file. This naming format allows for up to 17576 files, ! ! which should be enough for most applications. ! ! ! ! The program searches through those "GRIBFILE.???" files, and pulls out all ! ! the requested fields which fall between a starting date and an ending date. ! ! It writes the fields from a given time period to a file named according to ! ! the date and hour, i.e., "FILE:YYYY-MM-DD_HH" ! ! ! !*****************************************************************************! program ungrib use table use gridinfo use storage_module use filelist use datarray use module_debug use misc_definitions_module use stringutil implicit none interface subroutine read_namelist(hstart, hend, delta_time, ntimes,& ordered_by_date, debug_level, out_format, prefix, & add_lvls, new_plvl, interp_type, ec_rec_len, pmin) use misc_definitions_module character(len=19) :: hstart, hend integer :: delta_time integer :: ntimes logical :: ordered_by_date integer :: debug_level character(len=3) :: out_format character(len=MAX_FILENAME_LEN) :: prefix logical :: add_lvls real, dimension(:) :: new_plvl integer :: interp_type, ec_rec_len real :: pmin end subroutine read_namelist end interface integer :: nunit1 = 12 character(LEN=132) :: gribflnm = 'GRIBFILE.AAA ' ! won't work with len=12 integer :: debug_level integer , parameter :: maxlvl = 250 real , dimension(maxlvl) :: plvl, new_plvl integer :: iplvl logical :: add_lvls integer :: interp_type, ec_rec_len real :: pmin integer :: nlvl real :: startlat, startlon, deltalat, deltalon real :: level character (LEN=9) :: field character (LEN=3) :: out_format character (LEN=MAX_FILENAME_LEN) :: prefix logical :: readit integer, dimension(255) :: iuarr = 0 character (LEN=19) :: HSTART, HEND, HDATE character(LEN=19) :: hsave = '0000-00-00_00:00:00' integer :: itime integer :: ntimes integer :: interval integer :: ierr logical :: ordered_by_date integer :: grib_version integer :: vtable_columns call mprintf(.true.,STDOUT,' *** Starting program ungrib.exe *** ') call mprintf(.true.,LOGFILE,' *** Starting program ungrib.exe *** ') ! ----------------- ! Read the namelist, and return the information we want: call read_namelist(hstart, hend, interval, ntimes, & ordered_by_date, debug_level, out_format, prefix, & add_lvls, new_plvl, interp_type, ec_rec_len, pmin) call mprintf(.true.,INFORM,"Interval value: %i seconds or %f hours", & i1=interval, f1=float(interval)/3600.) call mprintf(.true.,STDOUT,'Path to intermediate files is %s',s1=get_path(prefix)) call mprintf(.true.,LOGFILE,'Path to intermediate files is %s',s1=get_path(prefix)) ! ----------------- ! Determine GRIB Edition number grib_version=0 ! NCAR/RDA's EC 113 data are grib1 with a non-standard record length if ( ec_rec_len .eq. 0 ) then call edition_num(nunit1, gribflnm, grib_version, ierr) call mprintf((ierr.eq.3),ERROR,"GRIB file problem") else grib_version=1 ierr=0 endif if (grib_version.eq.2) then vtable_columns=12 call mprintf(.true.,INFORM, & "Linked in png and jpeg libraries for Grib Edition 2") else vtable_columns=7 endif call mprintf(.true.,INFORM,"Reading Grib Edition %i", i1=grib_version) ! ----------------- ! Read the "Vtable" file, and put the information contained therein into ! the module "table". call parse_table(debug_level,vtable_columns) call mprintf(.true.,DEBUG,"Parsed the vtable.") ! ----------------- ! Initialize the input filename to GRIBFILE.AA{character just before A} ! That way, when we update the filename below for the first time, it will ! have the correct value of GRIBFILE.AAA. gribflnm(12:12)=char(ichar(gribflnm(12:12))-1) ! ----------------- ! LOOP2 cycles through the list of files named GRIBFILE.???, until it finds ! a non-existent file. Then we exit LOOP2 : DO ! At the beginning of LOOP2 update the input filename. if (gribflnm(12:12).eq.'Z') then if (gribflnm(11:11).eq.'Z') then gribflnm(10:10) = char(ichar(gribflnm(10:10))+1) gribflnm(11:11) = 'A' else gribflnm(11:11) = char(ichar(gribflnm(11:11))+1) endif gribflnm(12:12) = 'A' else gribflnm(12:12) = char(ichar(gribflnm(12:12))+1) endif ! Set READIT to .TRUE., meaning that we have not read any records yet ! from the file GRIBFLNM. call mprintf(.true.,DEBUG,"Reading from gribflnm %s ",s1=gribflnm) readit = .TRUE. ! i.e., "Yes, we want to read a record." ! LOOP1 reads through the file GRIBFLNM, exiting under two conditions: ! 1) We have hit the end-of-file ! 2) We have read past the ending date HEND. ! ! Condition 2 assumes that the data in file GRIBFLNM are ordered in time. LOOP1 : DO ! At the beginning of LOOP1, we are at a new time period. ! Clear the storage arrays and associated level information. nlvl = 0 plvl = -999. call clear_storage ! LOOP0 reads through the file GRIBFLNM, looking for data of the current ! date. It exits under the following conditions. ! 1) We have hit the end-of-file ! 2) The GRIBFLNM variable has been updated to a nonexistent file. ! 3) We start reading a new date and the data are assumed to be ! ordered by date. ! 4) We have a valid record and the data are not assumed to be ! ordered by date. LOOP0 : DO ! If we need to read a new grib record, then read one. if (READIT) then if (grib_version.ne.2) then call mprintf(.true.,DEBUG, & "Calling rd_grib1 with iunit %i", i1=nunit1) call mprintf(.true.,DEBUG, & "flnm = %s",s1=gribflnm) ! Read one record at a time from GRIB1 (and older Editions) call rd_grib1(nunit1, gribflnm, level, field, & hdate, ierr, iuarr, debug_level, ec_rec_len, pmin) else ! Read one file of records from GRIB2. call mprintf(.true.,DEBUG,"Calling rd_grib2") call rd_grib2(nunit1, gribflnm, hdate, & grib_version, ierr, debug_level, pmin) FIELD='NULL' endif call mprintf(.true.,DEBUG,"ierr = %i ",i1=ierr) if (ierr.eq.1) then ! We have hit the end of a file. Exit LOOP0 so we can ! write output for date HDATE, and then exit LOOP1 ! to update the GRIBFLNM hsave = hdate exit LOOP0 endif if (ierr.eq.2) then ! We have hit the end of all the files. We can exit LOOP2 ! because there are no more files to read. exit LOOP2 endif call mprintf(.true.,DEBUG, & "Read a record %s with date %s", s1=field,s2=hdate(1:13)) endif call mprintf(.true.,DEBUG, & "hdate = %s , hsave = %s ",s1=hdate(1:13), s2=hsave(1:13) ) ! if (hdate < hstart) then ! ! The data read has a date HDATE earlier than the starting ! ! date HSTART. So cycle LOOP0 to read the the next GRIB record. ! READIT = .TRUE. ! cycle LOOP0 ! endif if (FIELD.EQ.'NULL') then ! The data read does not match any fields requested ! in the Vtable. So cycle LOOP0 to read the next GRIB record. READIT = .TRUE. cycle LOOP0 endif if (ordered_by_date .and. (hdate > hsave)) then ! Exit LOOP0, because we started to read data from another ! date. call mprintf(.true.,DEBUG, & "hdate %s > hsave %s so exit LOOP0",s1=hdate,s2=hsave) ! We set READIT to FALSE because we have not yet processed ! the data from this record, and we will want to process this ! data on the next pass of LOOP1 (referring to the next time ! period of interest. READIT = .FALSE. exit LOOP0 endif ! When we have reached this point, we have a data array ARRAY which has ! some data we want to save, with field name FIELD at pressure level ! LEVEL (Pa). Dimensions of this data are map%nx and map%ny. Put ! this data into storage. if (((field == "SST").or.(field == "SKINTEMP")) .and. & (level /= 200100.)) level = 200100. iplvl = int(level) call mprintf((.not.allocated(rdatarray)),ERROR, & "GRIB data slab not allocated in ungrib.F before call to put_storage.") call put_storage(iplvl,field, & reshape(rdatarray(1:map%nx*map%ny),(/map%nx, map%ny/)),& map%nx,map%ny) deallocate(rdatarray) ! Since we processed the record we just read, we set ! READIT to .TRUE. so that LOOP0 will read the next record. READIT = .TRUE. if (.not. ordered_by_date) then if (hdate >= hstart) then hsave = hdate endif exit LOOP0 endif enddo LOOP0 ! When we have reached this point, we have either hit the end of file, or ! hit the end of the current date. Either way, we want to output ! the data found for this date. This current date is in HSAVE. ! However, if (HSAVE == 0000-00-00_00:00:00), no output is necessary, ! because that 0000 date is the flag for the very opening of a file. if ((hsave(1:4).ne.'0000').and.(hsave.le.hend)) then if (debug_level .gt. 100) then print*, 'Calling output: '//hsave(1:13) call mprintf(.true.,DEBUG,"Calling output: %s ",s1=hsave(1:13)) endif call output(hsave, nlvl, maxlvl, plvl, interval, 1, out_format, prefix, debug_level) hsave=hdate ! If the next record we process has a date later than HEND, ! then time to exit LOOP1. if ((ordered_by_date) .and. (hdate.gt.hend)) exit LOOP1 else hsave = hdate endif ! If we hit the end-of-file, its time to exit LOOP1 so we can ! increment the GRIBFLNM to read the next file. if (ierr.eq.1) exit LOOP1 enddo LOOP1 ! When we have reached this point, we read all the data we want to from ! file GRIBFLNM (either because we reached the end-of-file, or we ! read past the date HEND). So we close the file and cycle LOOP2 to read ! the next file. if (grib_version.ne.2) then call c_close(iuarr(nunit1), debug_level, ierr) iuarr(nunit1) = 0 endif hsave = '0000-00-00_00:00:00' enddo LOOP2 ! Now Reread, process, and reoutput. call mprintf(.true.,INFORM,"First pass done, doing a reprocess") call rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, & add_lvls, new_plvl, interp_type, & debug_level, out_format, prefix) ! Make sure the filedates are in order, with an inefficient sort: call sort_filedates ! Interpolate temporally to fill in missing times: call datint(filedates, nfiles, hstart, ntimes, interval, out_format, prefix) ! Now delete the temporary files: call file_delete(filedates, nfiles, trim(get_path(prefix))//'PFILE:', interval) ! And Now we are done: call mprintf(.true.,STDOUT,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') call mprintf(.true.,STDOUT,'! Successful completion of ungrib. !') ! call mprintf(.true.,STDOUT,"! We're hauling gear at Bandimere. !") call mprintf(.true.,STDOUT,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') call mprintf(.true.,LOGFILE,' *** Successful completion of program ungrib.exe *** ') contains subroutine sort_filedates implicit none integer :: n logical :: done if (nfiles > 1) then done = .FALSE. do while ( .not. done) done = .TRUE. do n = 1, nfiles-1 if (filedates(n) > filedates(n+1)) then filedates(size(filedates)) = filedates(n) filedates(n) = filedates(n+1) filedates(n+1) = filedates(size(filedates)) filedates(size(filedates)) = '0000-00-00 00:00:00.0000' done = .FALSE. endif enddo enddo endif end subroutine sort_filedates end program ungrib