! Public Domain 4/2017 Wesley Ebisuzaki
!
!	The grib file contains 2000 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 merged_.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

	character (len=200) :: metadata, time
	character (len=400) :: file, template

	file = 'Z__C_RJTD_20170420043500_NOWC_GPV_Ggis0p25km_Prr05lv_Aper5min_FH0000-0030_grib2.bin'
	template = 'japan_big_grid.grb2'
!        time=':anl:'
	time=':55 min fcst:'

!	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
	if (big_nx /= 2280) stop 3
	if (big_ny /= 3081) stop 4
	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 grid size
	n_fields = grb2_inq(file,'@mem:1', ':APCP:', time, ':npts=1600:')
	write(*,*) 'found ',n_fields, ' that match :APCP:, :npts=1600: and ',trim(time)

	do i = 1, n_fields
	   iret = grb2_inq(file,'@mem:1', ':APCP:', time, ':npts=1600:', sequential=i-1, data2=grid, &
		lat=lat, lon=lon, nx=nx, ny=ny, desc=metadata)
	   if (iret.ne.1) stop 6
	   if (nx /= 40) stop 7
	   if (ny /= 40) stop 8

!	   find position of small grid within large grid
	   i0 = int( (lon(1,1) - 120.506250d0) / 0.012500d0 + 0.5d0)
	   j0 = int( (lat(1,1) - 21.004166d0) / 0.00833333333 + 0.5d0)

!	   write(*,*) 'lat/lon=',lat(1,1), lon(1,1), ' i0/j0=',i0,j0

	   if (i0.ge.0 .and. i0+40 .le. big_nx .and. j0.ge.0 .and. j0+40 .le.big_ny) then
	       do iy = 1, 40
	           do ix = 1, 40
		       big_grid(ix + i0, iy + j0) = grid(ix, iy)
	           enddo
	       enddo
	    else
		write(*,*) 'outsided of big domain i0/j0=',i0,j0
	    endif

	enddo

!	write big_grid with 20 bit precision

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

	stop
	end