! ****************************************************************************! ! 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 include "mpif.h" integer :: nunit1 = 12 character(LEN=132) :: gribflnm = 'GRIBFILE.AAA ' ! won't work with len=12 integer :: debug_level, KMAX, LUN1 integer , parameter :: maxlvl = 100 real, allocatable, save, dimension(:,:,:) :: TGRID, UGRID, & VGRID,QGRID, & CWMGRID, PRESGRID, & PINTGRID, ZGRID real , dimension(maxlvl) :: plvl integer :: iplvl, mype, npes, I,J integer :: nlvl, icount real :: startlat, startlon, deltalat, deltalon real :: level character (LEN=9) :: field character (LEN=3) :: out_format character (LEN=MAX_FILENAME_LEN) :: prefix logical :: readit, spectral 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, ierrmpi 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 *** ') call mpi_init(ierrmpi) call mpi_comm_rank(MPI_COMM_WORLD,mype,ierrmpi) call mpi_comm_size(MPI_COMM_WORLD,npes,ierrmpi) ! ----------------- ! Read the namelist, and return the information we want: call read_namelist(hstart, hend, interval, ntimes, & ordered_by_date, debug_level, out_format, prefix, spectral) 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 call edition_num(nunit1, gribflnm, grib_version, ierr) !!! hardwire grib version grib_version=2 !!! call mprintf((ierr.eq.3),ERROR,"GRIB file problem") if (grib_version.eq.2) then vtable_columns=11 #if defined (USE_PNG) && (USE_JPEG2000) call mprintf(.true.,INFORM, & "Linked in png and jpeg libraries for Grib Edition 2") #else call mprintf(.true.,STDOUT,"WARNING - Grib Edition 2 data detected, and") call mprintf(.true.,STDOUT," - png and jpeg libs were NOT selected") call mprintf(.true.,STDOUT," - during the build.") call mprintf(.true.,STDOUT,"Stopping") call mprintf(.true.,LOGFILE,"WARNING - Grib Edition 2 data detected, and") call mprintf(.true.,LOGFILE," - png and jpeg libs were NOT selected") call mprintf(.true.,LOGFILE," - during the build.") call mprintf(.true.,LOGFILE,"Stopping") call mprintf(.true.,ERROR,"NEED_GRIB2_LIBS") #endif 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 icount=-1 LOOP2 : DO print*, 'starting LOOP2' icount=icount+1 ! 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 print*, 'starting LOOP1' ! 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 ! write(0,*) 'start LOOP0 with MYPE: ', MYPE ! 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) ! write(0,*) 'back from rd_grib1 with field,hdate: ', field,hdate 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) 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 print*, 'exiting LOOP0...ierr=1 ', hsave 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. print*, 'exiting LOOP2...ierr=2' 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.") !mp call mprintf(.TRUE.,STDOUT,"(2)call put_storage...iplvl: %i field: %s ", i1=iplvl, s1=field) ! call mprintf(.TRUE.,STDOUT,"call put_storage...") !mp 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 ! print*, 'past 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 LUN1=27+icount print*, 'icount, lun1: ', icount, lun1 SPECT: if (spectral) then KMAX=64 if (.not. allocated(TGRID)) then ALLOCATE(TGRID(map%nx,map%ny,KMAX)) ALLOCATE(UGRID(map%nx,map%ny,KMAX)) ALLOCATE(VGRID(map%nx,map%ny,KMAX)) ALLOCATE(QGRID(map%nx,map%ny,KMAX)) ALLOCATE(CWMGRID(map%nx,map%ny,KMAX)) ALLOCATE(PRESGRID(map%nx,map%ny,KMAX)) ALLOCATE(PINTGRID(map%nx,map%ny,KMAX+1)) ALLOCATE(ZGRID(map%nx,map%ny,KMAX+1)) endif print*, 'call COF2GRD_NPS with LUN1', lun1, MYPE call COF2GRD_NPS(LUN1,map%nx, map%ny, KMAX,TGRID,ZGRID, & UGRID, VGRID, QGRID, CWMGRID, PRESGRID, PINTGRID) print*, 'return COF2GRD_NPS:', MYPE if (MYPE .eq. 0) then print*, 'to N_S flipping routine' call N_S_FLIP(TGRID,map%nx,map%ny,KMAX) call N_S_FLIP(UGRID,map%nx,map%ny,KMAX) call N_S_FLIP(VGRID,map%nx,map%ny,KMAX) call N_S_FLIP(QGRID,map%nx,map%ny,KMAX) call N_S_FLIP(PRESGRID,map%nx,map%ny,KMAX) do iplvl=1,KMAX field='TT' call put_storage(iplvl,field, & TGRID(:,:,iplvl),& map%nx,map%ny) field='UU' call put_storage(iplvl,field, & UGRID(:,:,iplvl),& map%nx,map%ny) field='VV' call put_storage(iplvl,field, & VGRID(:,:,iplvl),& map%nx,map%ny) field='SPECHUMD' call put_storage(iplvl,field, & QGRID(:,:,iplvl),& map%nx,map%ny) field='PRESSURE' call put_storage(iplvl,field, & PRESGRID(:,:,iplvl),& map%nx,map%ny) enddo call N_S_FLIP(PINTGRID,map%nx,map%ny,KMAX+1) call N_S_FLIP(ZGRID,map%nx,map%ny,KMAX+1) do iplvl=1,KMAX+1 field='PINT' call put_storage(iplvl,field, & PINTGRID(:,:,iplvl),& map%nx,map%ny) field='HGT' call put_storage(iplvl,field, & ZGRID(:,:,iplvl),& map%nx,map%ny) enddo do J=1,map%ny do I=1,map%nx if (ZGRID(I,J,KMAX) .lt. ZGRID(I,J,KMAX+1)) then write(0,*) 'I,J, ZGRID(I,J,1:5): ', I,J, ZGRID(I,J,KMAX-4:KMAX+1) write(0,*) 'I,J, PINTGRID(I,J,1:5): ', I,J, PINTGRID(I,J,KMAX-4:KMAX+1) endif enddo enddo endif ! MYPE=0 endif SPECT !what from? if (MYPE .eq. 0) then if (debug_level .gt. -1) print*, '--- Calling output: '//hsave(1:13) call output(hsave, nlvl, maxlvl, plvl, interval, 1, out_format, prefix, debug_level) hsave=hdate endif ! write(0,*) 'to barrier(1): ', MYPE call MPI_BARRIER(MPI_COMM_WORLD,ierrmpi) ! write(0,*) 'past barrier(1): ', MYPE ! 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 print*, 'hsave now: ', hsave 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 ! write(0,*) 'to barrier(2): ', MYPE call MPI_BARRIER(MPI_COMM_WORLD,ierrmpi) ! write(0,*) 'past barrier(2): ', MYPE ! 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' print*, 'end LOOP2' enddo LOOP2 ! Now Reread, process, and reoutput. call mprintf(.true.,INFORM,"First pass done, doing a reprocess") call MPI_BARRIER(MPI_COMM_WORLD,ierrmpi) if (MYPE .eq. 0) then call rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, 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) endif ! MYPE=0 ! 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 *** ') call close_logfile() call MPI_BARRIER(MPI_COMM_WORLD,ierrmpi) ! write(0,*) 'past barrier(5): ', MYPE call mpi_finalize(ierrmpi) 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 subroutine N_S_FLIP(A,nx,ny,nz) real A(nx,ny,nz),tmp(nx,ny,nz) integer:: I,J,JJ,K,nx,ny,nz do K=1,nz do J=1,ny JJ=ny-J+1 do I=1,nx tmp(I,JJ,K)=A(I,J,K) enddo enddo enddo A=tmp write(0,*) 'at K=1' K=1 do J=ny,1,-ny/20 write(0,617) (A(I,J,K),I=1,nx,nx/15) enddo K=nz write(0,*) 'at K=nz' do J=ny,1,-ny/20 write(0,617) (A(I,J,K),I=1,nx,nx/15) enddo 617 format(20(f8.1,1x)) end subroutine N_S_FLIP end program ungrib