! 1/2018 Public Domain Wesley Ebisuzaki
!
!	Compute average rh from surface to 700 mb and write as grib file
!
!       input: grib file with RH on 1000, 925 850 and 700 mb
!       output: grib file with pressure weighted average RH from
!           surface to 700 mb
!
!       files:
!           in = input grib2 file that has 1000, 825, 850 and 700 mb RH
!                   (only one field at each level)
!           inv = inventory file
!           out = output grib2 with average RH
!
!       For RH below 1000 mb, use 1000 mb value
!
!       requires: wgrib2api built with wgrib2 v2.0.8 which supports
!          'surface - 700 mb'
!          grb_UNDEFINED
!
        use wgrib2api
        character (len=100) :: in, out, inv
        character (len=200) :: metadata
        real, allocatable :: sfc_prs(:,:), tmp(:,:), rh(:,:,:), ave_rh(:,:)
        integer, parameter :: nlevs = 4
	real, dimension (0:nlevs), parameter :: levels = (/ 2000.0, 1000.0, 925.0, 850.0, 700.0 /)
        character (len=9), dimension(0:nlevs), parameter :: &
            txt_levs = (/ ':2000 mb:', ':1000 mb:', ':925 mb: ', ':850 mb: ', ':700 mb: ' /)
        character (len=30), parameter :: out_level = 'surface - 700 mb'

	integer :: lev, nx, ny, i, j
	real:: factor, ave_val, lower_val

	in='/export/cpc-lw-webisuzak/wd51we/grib2/examples/gep19.aec'
	inv='gep19.aec.inv'
	out='ave_rh.grb'

!	make inventory
	i = grb2_mk_inv(in,inv)
	if (i /= 0) stop 1

!	read surface pressure
	i = grb2_inq(in, inv, ':PRES:', ':surface:', data2=sfc_prs, nx=nx, ny=ny)
	if (i /= 1) stop 2
	sfc_prs = 0.01 * sfc_prs		! convert from Pa to mb
	write(*,*) 'read sfc pressure (mb) done'

!	read rh 
	allocate(rh(nx,ny,0:nlevs))
	do lev = 1, nlevs
	   i = grb2_inq(in, inv, ':RH:', txt_levs(lev), data2=tmp, desc=metadata)
	   if (i /= 1) stop 3
	   rh(:,:,lev) = tmp
	enddo
	rh(:,:,0) = rh(:,:,1)			! rh(2000mb) = rh(1000mb)
	write(*,*) 'read rh finished'

!	now find the ave RH, surface-700 mb, pressure weighted
	allocate(ave_rh(nx,ny))
	ave_rh = 0.0
	do j = 1, ny
	    do i = 1, nx
		if (sfc_prs(i,j).gt.2000) stop 4
		if (sfc_prs(i,j).lt.levels(nlevs)) then
        	    ave_rh(i,j) = grb2_UNDEFINED
		    cycle
		endif
		do lev = 1, nlevs
		    if (sfc_prs(i,j) .gt. levels(lev-1)) then
!			include whole layer
			ave_rh(i,j) = ave_rh(i,j) + (rh(i,j,lev-1)+rh(i,j,lev))*0.5*(levels(lev-1)-levels(lev))
		    else if (sfc_prs(i,j) .gt. levels(lev)) then
!			include partial layer
			factor = (sfc_prs(i,j)-levels(lev)) / (levels(lev-1) - levels(lev))
			lower_val = factor*rh(i,j,lev-1) + (1.0-factor)*rh(i,j,lev)
		        ave_val = 0.5*(lower_val + rh(i,j,lev))
			ave_rh(i,j) = ave_rh(i,j) + ave_val * (sfc_prs(i,j) - levels(lev))
		    endif
		enddo
!	        normalize ave_rh 
		ave_rh(i,j) = ave_rh(i,j) / (sfc_prs(i,j) - levels(nlevs))
	    enddo
	enddo

!	write ave_rh
	i = grb2_wrt(out,in,1,data2=ave_rh,meta=metadata,level=out_level)
	if (i /= 0) stop 5

	stop
	end