!##############################################################################
Program griber

implicit none

character*10 c_date,c_hr,c_datatype
character*128 filein,fileout
character*2 cflag(10)
character*20 projection
character(len=20), allocatable :: fields(:),var3d(:), var2d(:), varsd(:)
integer, allocatable :: levs(:),irecs(:),idates(:),ifhrs(:)
real, allocatable :: a(:,:)

integer :: nx,ny,nrec,longdate,mode &
          ,ia,iargc,istr,iend,inc,i,j,lcg,ndate,nhr &
          ,nlev,slev,inproj,miss,startlev,endlev &
          ,iyyyy,imm,idd,ihh,fyyyy,fmm,fdd,fhh,leapyear &
          ,nvar3d,nvar2d,nvarsd,iplevs(1000),islevs(1000),datatype
real :: alat1,alon1,alat2,alon2,alov,aorient,dx,dy,reflat1,reflat2 &
       ,tinc,amiss

data amiss /-999./

!******************************************************************************
!Check for command line arguments
!******************************************************************************
IA = IARGC ( )
IF (IA .lt. 8) THEN
   WRITE(*,'("Usage:dgrib [-t #] [-d YYYYMMDDHH] [-h fcsthr] [-f gribfile]")')
   WRITE(*,'("-t option: # = 1-NcepReanalysis, 2-GDAS-FNL, 3-NARR")')
   WRITE(*,'("-d option: YYYYMMDDHH is the time")')
   WRITE(*,'("-h option: forecast hour (should be 0 for reanalysis)")')
   WRITE(*,'("-f option: gribfile is the name of the file to extract data")')
   STOP
ENDIF

!Set default date and forecast hour
c_date='99999999'
c_hr  ='99999999'

!Read command line arguments
do i=1,ia,2
   call ugetarg(i,cflag(i))
   if(cflag(i).eq.'-t')then
      call ugetarg(i+1,c_datatype)
   elseif(cflag(i).eq.'-d')then
      call ugetarg(i+1,c_date)
   elseif(cflag(i).eq.'-h')then
      call ugetarg(i+1,c_hr)
   elseif(cflag(i).eq.'-f')then
      call ugetarg(i+1,filein)
   endif
enddo

!wgrib only prints out the 2 digit year in short inventory
read(c_date(3:),'(i10)')ndate
read(c_hr,'(i10)')nhr
read(c_datatype,'(i10)')datatype
!Print out command line inputs
lcg=index(filein,' ')-1
print*,'gribfile=',filein(1:lcg)
print*,'c_date=',c_date,ndate
print*,'c_hr=',c_hr,nhr
print*,'c_datatype=',c_datatype,datatype

!******************************************************************************
! Set datatype:
! 1 = NCEPreanalysis 2.5deg
! 2 = GDAS-FNL 1.0deg
! 3 = NARR 32km
! 4 = ECMWF-ERA-Interim 1.5deg
! 5 = RAP Rapid Refresh Analysis 32km (awips grid 221)
! 6 = ERA-Interim (ECMWF) (Z,T,R,U,V naming convention)
! 7 = HRRR forecast grids (3km CONUS)
!     GFS forecast grids (multiple resolutions)
!******************************************************************************
!******************************************************************************
if(datatype==1)then  !NCEPreanalysis global 2.5deg lat-lon
 !3D atmospheric variables
 nvar3d=5
 allocate (var3d(nvar3d))
 var3d(1)='UGRD'
 var3d(2)='VGRD'
 var3d(3)='TMP'
 var3d(4)='HGT'
 var3d(5)='RH'
 !2D snow water and depth variables here if available
 nvar2d=2
 allocate (var2d(nvar2d))
 var2d(1)='NOTHING'          !Include accumulated snow mass (kg/m^2)
 var2d(2)='NOTHING'          !Include Snow depth [m]
 !SOIL moisture variables here if available
 nvarsd=2
 allocate (varsd(nvarsd))
 varsd(1)='NOTHING'          !Include Volumetric Soil Moisture (m3/m3)(fraction)
 varsd(2)='NOTHING'          !Include Soil temperature (K)

!******************************************************************************
elseif(datatype==2)then  !GFS GDAS-FNL global 1deg lat-lon
 !3D atmospheric variables
 nvar3d=5
 allocate (var3d(nvar3d))
 var3d(1)='UGRD'
 var3d(2)='VGRD'
 var3d(3)='TMP'
 var3d(4)='HGT'
 var3d(5)='RH'
 !2D snow water and depth variables here if available
 nvar2d=2
 allocate (var2d(nvar2d))
 var2d(1)='WEASD'          !Accum. snow [kg/m^2]
 var2d(2)='WEASD'          !Use WEASD to estimate Snow depth [m]
 !SOIL moisture variables here if available
 !The GFS/GDAS FNL 1deg analysis has 0-10cm, 10-40cm, 40-100cm, 100-200cm soil
 !depths in the dataset starting 2005-06-01. We write out top 2 soil levels.
 !expected by RAMS. Prior to and including 2005-05-31 there is only 1 soil level.
 !We accomodate for that below.
 nvarsd=2
 allocate (varsd(nvarsd))
 varsd(1)='SOILW'          !Volumetric Soil Moisture (m3/m3)(fraction)
 varsd(2)='TMP'            !Soil temperature (K)

!******************************************************************************
elseif(datatype==3)then !NARR data
 !3D atmospheric variables
 nvar3d=5
 allocate (var3d(nvar3d))
 var3d(1)='UGRD'
 var3d(2)='VGRD'
 var3d(3)='TMP'
 var3d(4)='HGT'
 var3d(5)='SPFH'
 !2D snow water and depth variables here if available
 nvar2d=2
 allocate (var2d(nvar2d))
 var2d(1)='WEASD'          !Accum. snow [kg/m^2]
 var2d(2)='SNOD'           !Snow depth [m]
 !SOIL moisture variables here if available
 !The NARR 32km data has 0-10cm, 10-40cm, 40-100cm, 100-200cm soil 
 !depths so write out top 2 levels expected by RAMS
 nvarsd=2
 allocate (varsd(nvarsd))
 varsd(1)='SOILW'          !Volumetric Soil Moisture (m3/m3)(fraction)
 varsd(2)='TSOIL'          !Soil temperature (K)

!******************************************************************************
elseif(datatype==4)then !ERA-Interim (ECMWF)
 !3D atmospheric variables
 nvar3d=5
 allocate (var3d(nvar3d))
 var3d(1)='UGRD'
 var3d(2)='VGRD'
 var3d(3)='TMP'
 var3d(4)='GP'
 var3d(5)='RH'
 !2D snow water and depth variables here if available
 nvar2d=2
 allocate (var2d(nvar2d))
 var2d(1)='SNOD'           !Use SNOD to estimate Accum. snow [kg/m^2]
 var2d(2)='SNOD'           !Snow depth [m]
 !SOIL moisture variables here if available
 !Top 2 ERA-Interim soil levels are 7-28cm and 0-7cm down
 nvarsd=2
 allocate (varsd(nvarsd))
 varsd(1)='VSOILM'         !Volumetric Soil Moisture (m3/m3)(fraction)
 varsd(2)='TSOIL'          !Soil temperature (K)

!******************************************************************************
elseif(datatype==5)then !RAP Rapid Refresh Analysis 32km (awips grid 221)
 !3D atmospheric variables
 nvar3d=5
 allocate (var3d(nvar3d))
 var3d(1)='UGRD'
 var3d(2)='VGRD'
 var3d(3)='TMP'
 var3d(4)='HGT'
 var3d(5)='RH'
 !2D snow water and depth variables here if available
 nvar2d=2
 allocate (var2d(nvar2d))
 var2d(1)='WEASD'          !Accum. snow [kg/m^2]
 var2d(2)='SNOD'           !Snow depth [m]
 !SOIL moisture variables here if available
 nvarsd=2
 allocate (varsd(nvarsd))
 varsd(1)='SOILW'          !Volumetric Soil Moisture (m3/m3)(fraction)
 varsd(2)='TSOIL'          !Soil temperature (K)

!******************************************************************************
elseif(datatype==6)then !ERA-Interim (ECMWF) (Z,T,R,U,V naming convention)
 !3D atmospheric variables
 nvar3d=5
 allocate (var3d(nvar3d))
 var3d(1)='U'
 var3d(2)='V'
 var3d(3)='T'
 var3d(4)='Z'
 var3d(5)='R'
 !2D snow water and depth variables here if available
 nvar2d=2
 allocate (var2d(nvar2d))
 var2d(1)='NOTHING'           !Use SNOD to estimate Accum. snow [kg/m^2]
 var2d(2)='NOTHING'           !Snow depth [m]
 !SOIL moisture variables here if available
 !Top 2 ERA-Interim soil levels are 7-28cm and 0-7cm down
 nvarsd=2
 allocate (varsd(nvarsd))
 varsd(1)='NOTHING'         !Volumetric Soil Moisture (m3/m3)(fraction)
 varsd(2)='NOTHING'          !Soil temperature (K)

!******************************************************************************
elseif(datatype==7)then !HRRR and GFS forecast grids
 !3D atmospheric variables
 nvar3d=5
 allocate (var3d(nvar3d))
 var3d(1)='UGRD'
 var3d(2)='VGRD'
 var3d(3)='TMP'
 var3d(4)='HGT'
 var3d(5)='RH'
 !2D snow water and depth variables here if available
 nvar2d=2
 allocate (var2d(nvar2d))
 var2d(1)='WEASD'          !Accum. snow [kg/m^2]
 var2d(2)='SNOD'           !Snow depth [m]
 !SOIL moisture variables here if available
 nvarsd=2
 allocate (varsd(nvarsd))
 varsd(1)='SOILW'          !Volumetric Soil Moisture (m3/m3)(fraction)
 varsd(2)='TSOIL'          !Soil temperature (K)
endif

!******************************************************************************
!Read some of the grib metadata
call grib_query (filein(1:lcg))            !Open the grib file
call grib_queryc('projection',projection)  !get Grid projection
call grib_queryf('lat1',alat1)             !get SW lat
call grib_queryf('lon1',alon1)             !get SW lon
call grib_queryf('lat2',alat2)             !get NE lat
call grib_queryf('lon2',alon2)             !get NE lon
call grib_queryf('dx',dx)                  !get Delta-X (km or deg)
call grib_queryf('dy',dy)                  !get Delta-Y (km or deg)
call grib_queryi('nx',nx)                  !get # of X grid points
call grib_queryi('ny',ny)                  !get # of Y grid points
allocate (a(nx,ny))
call grib_queryi('nrec',nrec)              !get # of records
call grib_queryi('longdate',longdate)      !get the date/time
call grib_queryi('mode',mode)              !get grid mode info (see grid tables)
call grib_queryf('orient',aorient)         
call grib_queryf('lov',alov)               !get grid orientation for LC projection
call grib_queryf('latin1',reflat1)         !get reference lat 1 for LC projection
call grib_queryf('latin2',reflat2)         !get reference lat 2 for LC projection
allocate (irecs(nrec))
allocate (levs(nrec))
allocate (idates(nrec))
allocate (ifhrs(nrec))
allocate (fields(nrec))
call grib_queryi('irecs',irecs)
call grib_queryi('levs',levs)
call grib_queryi('idates',idates)
call grib_queryi('ifhrs',ifhrs)
call grib_queryc('fields',fields)

print*,'DX,DY: ',dx,dy
print*,'NX,NY: ',nx,ny
print*,'LAT1,LON1: ',alat1,alon1
print*,'LAT2,LON2: ',alat2,alon2
print*,'PROJECTION: ',projection
print*,'LOV: ',alov                   !??? grid rotation
print*,'Orient: ',aorient             !??? grid rotation
print*,'Mode: ',mode
print*,'Nrec: ',nrec
print*,'Irecs: ',irecs(nrec)

!Do the dprep stuff

!First some grid navigation bookeeping
if (projection(1:6)=="latlon") then
   inproj=1
   !Will have to flip the latlon projections,
   !Will modify start lat/lon to be -90;-180
   !RAMS wants coordinates -180 to 180 longitude (not 0 to 360)
   if(mode==128)then
     alat1=-alat1
     alon1=alon1-180.
     alat2=-alat2
     if(alon2 <  180.) alon2=alon2+180.
     if(alon2 >= 180.) alon2=alon2-180.
   endif
elseif (projection(1:7)=="Lambert") then
   inproj=2
   alat2=0.0 !NE lat
   alon2=0.0 !NE lon
   dx=dx*1000.
   dy=dy*1000.
   aorient=alov
elseif (projection(1:5)=="polar") then
   inproj=3
   !Assume the polar stereo is true polar stereo
   reflat1=90.
   reflat2=amiss
endif

!Quick check to see if there are any fields matching time wanted
j=0
do i=1,nrec
   if(('HGT'==fields(i).or.'GP'==fields(i).or.'Z'==fields(i)) &
       .and.ndate==idates(i) &
       .and.nhr==ifhrs(i))then
      j=1
      exit
   endif
enddo
if(j==0)stop "NOTHING FOUND"

!**************************************************************************
!Determine the number of Pressure levels we are going to use
nlev=0
do i=1,nrec
   if(levs(i)==0)cycle  !don't count the surface fields here
   if((fields(i)=='HGT'.or.fields(i)=='GP'.or.fields(i)=='Z') &
       .and.ndate==idates(i) &
       .and.nhr==ifhrs(i))then
      nlev=nlev+1
      iplevs(nlev)=levs(i)
   endif
enddo

if(iplevs(1) < iplevs(nlev))then
   nlev=0
   do i=nrec,1,-1
      if(levs(i)==0)cycle  !don't count the surface fields here
      if((fields(i)=='HGT'.or.fields(i)=='GP'.or.fields(i)=='Z') &
          .and.ndate==idates(i) &
          .and.nhr==ifhrs(i))then
         nlev=nlev+1
         iplevs(nlev)=levs(i)
      endif
   enddo
endif
print*,'Number of Pressure Levels= ',nlev

!**************************************************************************
!Determine the number of SOIL levels we are going to use.
!Only add soil (moisture and/or temperature) if the soil moisture is there.
slev=0
do i=1,nrec
   if((fields(i)=='SOILW'.or.fields(i)=='VSOILM').and.ndate==idates(i) &
       .and.nhr==ifhrs(i))then
      slev=slev+1
      islevs(slev)=levs(i)
   endif
enddo
!Invert soil levels so level=1 is the bottom soil level (similar to RAMS)
if(abs(islevs(1)) < abs(islevs(slev)))then
   slev=0
   do i=nrec,1,-1
      if((fields(i)=='SOILW'.or.fields(i)=='VSOILM').and.ndate==idates(i) &
          .and.nhr==ifhrs(i))then
         slev=slev+1
         islevs(slev)=levs(i)
      endif
   enddo
endif
print*,'Number of Soil Levels= ',slev

!**************************************************************************
!Set times for writing to file and dp header info
read(c_date(1:4) ,'(i4)')iyyyy
read(c_date(5:6) ,'(i2)')imm
read(c_date(7:8) ,'(i2)')idd
read(c_date(9:10),'(i2)')ihh

leapyear=0
!leap year check. Every 4th year is a leap year except years divisible by 100
!unless they are also divisible by 400.
if(  mod(iyyyy, 4)==0) then
  !we could be a leap year. Keep checking.
  leapyear=1
  if(mod(iyyyy, 100)==0) then
    !maybe not
    leapyear=0
    if(mod(iyyyy,400)==0) then
      !yes leap year
      leapyear=1
    endif
  endif
endif

!adjust date if using forecast grids
fyyyy = iyyyy
fmm = imm
fhh = ihh + mod(nhr,24)
fdd = idd + (nhr/24)
if(fhh>=24) then
  fdd = fdd+fhh/24
  fhh = mod(fhh,24)
endif
if(imm==1 .and. fdd>31)then
 fmm=2
 fdd=fdd-31
elseif(imm==2 .and. fdd>28 .and. leapyear==0)then
 fmm=3
 fdd=fdd-28
elseif(imm==2 .and. fdd>29 .and. leapyear==1)then
 fmm=3
 fdd=fdd-29
elseif(imm==3 .and. fdd>31)then
 fmm=4
 fdd=fdd-31
elseif(imm==4 .and. fdd>30)then
 fmm=5
 fdd=fdd-30
elseif(imm==5 .and. fdd>31)then
 fmm=6
 fdd=fdd-31
elseif(imm==6 .and. fdd>30)then
 fmm=7
 fdd=fdd-30
elseif(imm==7 .and. fdd>31)then
 fmm=8
 fdd=fdd-31
elseif(imm==8 .and. fdd>31)then
 fmm=9
 fdd=fdd-31
elseif(imm==9 .and. fdd>30)then
 fmm=10
 fdd=fdd-30
elseif(imm==10 .and. fdd>31)then
 fmm=11
 fdd=fdd-31
elseif(imm==11 .and. fdd>30)then
 fmm=12
 fdd=fdd-30
elseif(imm==12 .and. fdd>31)then
 fmm=1
 fdd=fdd-31
 fyyyy=iyyyy+1
endif

write(fileout,45 ) 'dp-p',fyyyy,'-',fmm,'-',fdd,'-',fhh,'00'

print*,'Writing file:',fileout
45   format(a4,i4.4,a1,i2.2,a1,i2.2,a1,i2.2,a2)

open(10,file=fileout,status='unknown',form='formatted')

!Do the following 2 lines to make values between -180 to +180 for RAMS
if(alon1>180.) alon1=alon1-360.
if(aorient>180.) aorient=aorient-360.

WRITE(10,900)iyyyy,imm,idd,ihh*100,nhr*100,nlev,nx,ny &
     ,inproj,dx,dx,alat1,alon1,alat2,alon2,reflat1,aorient &
     ,(iplevs(i),i=1,nlev)
900  FORMAT(8I6,/,i1,2x,8(F10.4,2x),/,(12I6))

!*********************************************************************
!*********************************************************************
!     Go through each variable at each level for writing.

!Uncomment if you are pulling 3D variables from grib (Saleeby 4-18-06)
do j=1,nlev
 do i=1,nvar3d
   a=amiss
   call getfield(a,var3d(i),nrec,irecs,fields,levs,idates,ifhrs &
                ,iplevs(j),ndate,nhr)
   call prepfield(a,nx,ny,projection,var3d(i),mode,amiss,'3d',datatype,i)
   write(10,'(8f12.5)')a
 enddo
enddo

!Uncomment if you are pulling SOIL variables from grib (Saleeby 4-18-06)
!Datasets vary in the number of soil levels they have. Write out the top
!2 levels if we have 2+. If we only have 1 soil level, write it 2 times.
!RAMS wants to ingest 2 soil levels. If data fields are empty RAMS uses
!a default.
do i=1,nvarsd
  startlev=max(1,slev-1)
  endlev=max(1,slev)
  do j=startlev,endlev
     a=amiss
     call getfield(a,varsd(i),nrec,irecs,fields,levs,idates,ifhrs &
                  ,islevs(j),ndate,nhr)
     call prepfield(a,nx,ny,projection,varsd(i),mode,amiss,'sd',datatype,i)
     write(10,'(8f12.5)')a
  enddo
  !Write out at least 2 levels. If only 1 soil level exists, repeat it.
  if(startlev==endlev)then
    do j=startlev,endlev
      a=amiss
      call getfield(a,varsd(i),nrec,irecs,fields,levs,idates,ifhrs &
                   ,islevs(j),ndate,nhr)
      call prepfield(a,nx,ny,projection,varsd(i),mode,amiss,'sd',datatype,i)
      write(10,'(8f12.5)')a
    enddo
  endif
enddo

!Uncomment if you are pulling 2D variables from grib (Saleeby 4-18-06)
do i=1,nvar2d
   a=amiss
   j=0 !specify j=0 for the surface
   call getfield(a,var2d(i),nrec,irecs,fields,levs,idates,ifhrs &
                ,j,ndate,nhr)
   call prepfield(a,nx,ny,projection,var2d(i),mode,amiss,'2d',datatype,i)
   write(10,'(8f12.5)')a
enddo

!*********************************************************************
!*********************************************************************

close(10)

!Write the tag file to indicate file is complete
open(10,file=fileout(1:len_trim(fileout))//'.tag',status='unknown' &
    ,form='formatted')
write(10,'(a)')'done'
close(10)

deallocate (a)
deallocate (fields)
deallocate (ifhrs)
deallocate (idates)
deallocate (levs)
deallocate (irecs)
deallocate (var3d)
if(allocated(var2d)) deallocate (var2d)
if(allocated(varsd)) deallocate (varsd)

END

!************************************************************************
Subroutine getfield (a,type,nrec,irecs,fields,levs,idates,ifhrs &
                     ,iplev,ndate,nhr)

implicit none

real, dimension(*) :: a
character*(*) type, fields(*)
integer, dimension(*) :: irecs,levs,idates,ifhrs
integer :: iplev,ndate,nhr,nrec,i

do i=1,nrec,1

   if(type==fields(i).and.iplev==levs(i).and.ndate==idates(i).and. &
      nhr==ifhrs(i))then
      call grib_get(a,irecs(i))
      print*,'matched ',fields(i),levs(i),idates(i),ifhrs(i)
      exit
   endif

enddo

return
END SUBROUTINE getfield

!************************************************************************
Subroutine prepfield (a,nx,ny,projection,type,mode,amiss,dimtype &
                     ,datatype,vars)

implicit none

integer nx,ny,datatype,vars
real, dimension(nx*ny) :: a
character*(*) type,projection,dimtype

real, allocatable :: b(:)
real amiss
integer mode

integer i

if(projection=='latlon'.and.mode==128)then
   print*,'   flipping the projection'
   allocate (b(nx*ny))
   call flip(nx,ny,a,b)
   a=b
   deallocate (b)
endif

do i=1,nx*ny
   if(a(i) >= 1.e10)a(i)=amiss
enddo

!Turn from relative humidity in percent to a fraction
if(type=='RH' .or. type=='R')then
   print*,'   factoring RH'
   do i=1,nx*ny
      if(a(i) > amiss)a(i)=a(i)/100.
   enddo
!ECMWF uses GEOPOTENTIAL as the geopotential height multiplied by the WMO
!defined gravity constant of 9.80665 m/s**2 which is constant for all 
!latitudes and all heights. We must adjust for geopotential height.
elseif(type=='GP' .or. type=='Z')then
   print*,'   convert geopotential to geopotential-height'
   do i=1,nx*ny
      if(a(i) > amiss)a(i)=a(i)/9.80665
   enddo
!Turn from Pa to hPa or millibars
elseif(type=='PRMSL'.or.type=='PRES')then
   print*,'   factoring to mb'
   do i=1,nx*ny
      if(a(i) > amiss)a(i)=a(i)/100.
   enddo
!For GDAS-FNL snowdepth, Convert the 2nd WEASD to a snow-depth(m) 
!assuming density of 10:1 (Convert mm-liquid-equiv to m-snowdepth)
elseif(datatype==2.and.dimtype=='2d'.and.vars==2.and.type=='WEASD')then
   do i=1,nx*ny
      if(a(i) > amiss)a(i)=a(i)/100.
   enddo
!For ERA-Interim snowmass, Convert the SNOD snow-depth(m) 
!assuming density of 10:1 (Convert m-snowdepth to mm-liquid-equiv)
elseif(datatype==4.and.dimtype=='2d'.and.vars==1.and.type=='SNOD')then
   do i=1,nx*ny
      if(a(i) > amiss)a(i)=a(i)*100.
   enddo
endif

return
END SUBROUTINE prepfield

!************************************************************************
Subroutine flip (nx,ny,vin,vout)

implicit none

integer nx,ny,i,j,ii,jj
real, dimension(nx,ny) :: vin, vout

!Adjust starting longitude. Grib starts at 0.0deg lon
!while RAMS starts at -180.0deg
do i=1,nx
  ii=nx/2+i
  if(ii.gt.nx)ii=i-nx/2
  !Flip N/S for grib1 format
  do j=1,ny
    jj=ny-j+1
    vout(i,j)=vin(ii,jj)
  enddo
enddo

return
END SUBROUTINE flip
