! Public Domain 4/2017 Wesley Ebisuzaki
!
!	The grib file contains many small grids which need to be merged onto a 
!       larger grid.  I have already prepared a template of larger grid which 
!       uses the same grid points as the smaller grids
!
!	The template was created by
!
! wgrib2 Z__C_RJTD_20170420043500_NOWC_GPV_Ggis0p25km_Prr05lv_Aper5min_FH0000-0030_grib2.bin \
!   -new_grid_winds earth -new_grid latlon 120.506250:2280:0.012500 21.004166:3081:0.0083333334 junk
! wgrib2 junk -rpn 0 -grib japan_big_grid.grb2
!
!	The big grid is saved in merged.grb2
!	GrADS cannot handle merged.grb2 because it uses a JMA local template. So you
!	make merged_grb2 changing the pdt to 0.  (This loses metadata.)
!
!	wgrib2 merged.grb2 -set_pdt +0 -grib mreged_.grb2

	use wgrib2api

	real, allocatable :: big_grid(:,:), big_lat(:,:), big_lon(:,:)
	real, allocatable :: grid(:,:), lat(:,:), lon(:,:)
	integer :: big_nx, big_ny, nx, ny
	integer :: iret, n_fields, i0, j0, i, ix, iy, narg
	real (kind=8) :: dlat, dlon, lat0, lon0, latn, lonn
	character (len=200) :: metadata, time, nxny
	character (len=400) :: file, template, options(10), varname, outFile

!	read command line
	narg = command_argument_count()
	if (narg /= 6) then
	    write(*,*) 'file template variablename ftime small_grid_size outputfile'
	    stop
	endif

	do i = 1, narg
	   call get_command_argument(i, options(i))
	enddo

 	file = options(1)
	template = options(2)
	varname = options(3)
	time = options(4)
	nxny =	':npts=' // trim(options(5)) // ':'
	outFile = options(6)

	write(*,*)
	write(*,*) 'input data file=',trim(file)
	write(*,*) 'template file=',trim(template)
	write(*,*) 'output data file=',trim(outFile)

!	make index for template of big grid
	iret = grb2_mk_inv(template, '@mem:0')
	if (iret.ne.0) stop 1

!	read template for big grid
	iret = grb2_inq(template,'@mem:0',':APCP:',data2=big_grid,lat=big_lat,lon=big_lon,nx=big_nx,ny=big_ny)
	if (iret.ne.1) stop 2

!	figure out grid parameters
	lon0 = big_lon(1,1)
	lonn = big_lon(big_nx,1)
	lat0 = big_lat(1,1)
	latn = big_lat(1,big_ny)

	dlon = (lonn-lon0) / (big_nx - 1)
	dlat = (latn-lat0) / (big_ny - 1)

	write(*,*) 'lower-left lat/lon=',big_lat(1,1), big_lon(1,1)
        big_grid = 9.999e20						! set to undefined values

!	make index file for data file
	iret = grb2_mk_inv(file, '@mem:1')
	if (iret.ne.0) stop 5

!	find number of small fields that match variable, time and number of grid points
	n_fields = grb2_inq(file,'@mem:1', varname, time, nxny)

	write(*,*) 'found ',n_fields, ' that match ', trim(varname), ' ', trim(time), ' ', trim(nxny)

	do i = 1, n_fields
	   iret = grb2_inq(file,'@mem:1', varname, time, nxny , sequential=i-1, data2=grid, &
		lat=lat, lon=lon, nx=nx, ny=ny, desc=metadata)
	   if (iret.ne.1) stop 6

!	   find position of small grid within large grid
	   i0 = int( (lon(1,1)-lon0)/dlon + 0.5d0)
	   j0 = int( (lat(1,1)-lat0)/dlat + 0.5d0)

	   write(*,*) 'grid ',i, ' i0/j0=',i0,j0
!	   if small grid is outside of big grid, cycle
	   if (i0 < 0 .or. j0 < 0) cycle
	   if (i0 + nx > big_nx) cycle
	   if (j0 + ny > big_ny) cycle

!	   test to see if grid matches
!	   precision is limited by single precision lat/lon returned by wgrib2api

	   write(*,*) 'grid ',i, lon(1,1),big_lon(i0+1,j0+1), lat(1,1),big_lat(i0+1,j0+1)
	   if ( abs(lon(1,1) - big_lon(i0+1,j0+1)) .gt. 0.01) cycle
	   if ( abs(lat(1,1) - big_lat(i0+1,j0+1)) .gt. 0.01) cycle

	   if ( abs(lon(nx,ny) - big_lon(i0+nx,j0+ny)) .gt. 0.01) cycle
	   if ( abs(lat(nx,ny) - big_lat(i0+nx,j0+ny)) .gt. 0.01) cycle


	   write(*,*) 'processing small grid ', i

	   do iy = 1, ny
	       do ix = 1, nx
	          big_grid(ix + i0, iy + j0) = grid(ix, iy)
	       enddo
	   enddo

	enddo

!	write big_grid with 20 bit precision

	metadata=trim(metadata) // 'grib_max_bits=24:encode 20 bits:'
	iret = grb2_wrt(outFile,template,1,data2=big_grid,meta=metadata)
	write(*,*) 'metadata=',trim(metadata)
	write(*,*) 'ouput is in ',trim(outFile)
	if (iret.ne.0) stop 10

	stop
	end