Introduction

The wgrib2 utility can do calculations but some calculations are more easily done in Fortran. This example came from a user question on how to calculate the average RH from the surface to 700 mb. I am using a sample file with the surface pressure, RH for the levels, 1000, 925, 850 and 700 mb. The proper tool for the job is wgrib2api.

     1	! 1/2018 Public Domain Wesley Ebisuzaki
     2	!
     3	!	Compute average rh from surface to 700 mb and write as grib file
     4	!
     5	!       input: grib file with RH on 1000, 925 850 and 700 mb
     6	!       output: grib file with pressure weighted average RH from
     7	!           surface to 700 mb
     8	!
     9	!       files:
    10	!           in = input grib2 file that has 1000, 825, 850 and 700 mb RH
    11	!                   (only one field at each level)
    12	!           inv = inventory file
    13	!           out = output grib2 with average RH
    14	!
    15	!       For RH below 1000 mb, use 1000 mb value
    16	!
    17	!       requires: wgrib2api built with wgrib2 v2.0.8 which supports
    18	!          'surface - 700 mb'
    19	!          grb_UNDEFINED
    20	!
    21	        use wgrib2api
    22	        character (len=100) :: in, out, inv
    23	        character (len=200) :: metadata
    24	        real, allocatable :: sfc_prs(:,:), tmp(:,:), rh(:,:,:), ave_rh(:,:)
    25	        integer, parameter :: nlevs = 4
    26		real, dimension (0:nlevs), parameter :: levels = (/ 2000.0, 1000.0, 925.0, 850.0, 700.0 /)
    27	        character (len=9), dimension(0:nlevs), parameter :: &
    28	            txt_levs = (/ ':2000 mb:', ':1000 mb:', ':925 mb: ', ':850 mb: ', ':700 mb: ' /)
    29	        character (len=30), parameter :: out_level = 'surface - 700 mb'
    30	
    31		integer :: lev, nx, ny, i, j
    32		real:: factor, ave_val, lower_val
    33	
    34		in='/export/cpc-lw-webisuzak/wd51we/grib2/examples/gep19.aec'
    35		inv='gep19.aec.inv'
    36		out='ave_rh.grb'
    37	
    38	!	make inventory
    39		i = grb2_mk_inv(in,inv)
    40		if (i /= 0) stop 1
    41	
    42	!	read surface pressure
    43		i = grb2_inq(in, inv, ':PRES:', ':surface:', data2=sfc_prs, nx=nx, ny=ny)
    44		if (i /= 1) stop 2
    45		sfc_prs = 0.01 * sfc_prs		! convert from Pa to mb
    46		write(*,*) 'read sfc pressure (mb) done'
    47	
    48	!	read rh 
    49		allocate(rh(nx,ny,0:nlevs))
    50		do lev = 1, nlevs
    51		   i = grb2_inq(in, inv, ':RH:', txt_levs(lev), data2=tmp, desc=metadata)
    52		   if (i /= 1) stop 3
    53		   rh(:,:,lev) = tmp
    54		enddo
    55		rh(:,:,0) = rh(:,:,1)			! rh(2000mb) = rh(1000mb)
    56		write(*,*) 'read rh finished'
    57	
    58	!	now find the ave RH, surface-700 mb, pressure weighted
    59		allocate(ave_rh(nx,ny))
    60		ave_rh = 0.0
    61		do j = 1, ny
    62		    do i = 1, nx
    63			if (sfc_prs(i,j).gt.2000) stop 4
    64			if (sfc_prs(i,j).lt.levels(nlevs)) then
    65	        	    ave_rh(i,j) = grb2_UNDEFINED
    66			    cycle
    67			endif
    68			do lev = 1, nlevs
    69			    if (sfc_prs(i,j) .gt. levels(lev-1)) then
    70	!			include whole layer
    71				ave_rh(i,j) = ave_rh(i,j) + (rh(i,j,lev-1)+rh(i,j,lev))*0.5*(levels(lev-1)-levels(lev))
    72			    else if (sfc_prs(i,j) .gt. levels(lev)) then
    73	!			include partial layer
    74				factor = (sfc_prs(i,j)-levels(lev)) / (levels(lev-1) - levels(lev))
    75				lower_val = factor*rh(i,j,lev-1) + (1.0-factor)*rh(i,j,lev)
    76			        ave_val = 0.5*(lower_val + rh(i,j,lev))
    77				ave_rh(i,j) = ave_rh(i,j) + ave_val * (sfc_prs(i,j) - levels(lev))
    78			    endif
    79			enddo
    80	!	        normalize ave_rh 
    81			ave_rh(i,j) = ave_rh(i,j) / (sfc_prs(i,j) - levels(nlevs))
    82		    enddo
    83		enddo
    84	
    85	!	write ave_rh
    86		i = grb2_wrt(out,in,1,data2=ave_rh,meta=metadata,level=out_level)
    87		if (i /= 0) stop 5
    88	
    89		stop
    90		end

line 21:     needed use wgrib2api
lines 38-40: make inventory file
lines 42-46: read surface pressure, get nx, ny of the grid
lines 48-56: read RH for 1000, 925, 850 and 700 mb
             save in RH(:,:,level)
             for computations make RH for 2000 mb for the case
              where surface pressure is greater than 1000 mb.
lines 58-83: compute surface-700 mb average RH.
line 63:     if surface pressure is 2000 mb, something is really wrong
lines 64-67: if surface pressure is less than 700, ave RH is undefined
              note: grib2_UNDEFINED requires wgrib2api v2.0.7
lines 68-79: contribution of each layer to ave_rh
lines 85-87: write out ave_rh
              note: level_out requires wgrib2api v2.0.7

Comments

I needed to add couple of fixes to wgrib2 to make this example work. First I added grb2_UNDEFINED and grb2_DEFINED_VAL(x) to wgrib2api. Next I had to enhance wgrib2 so that -set_lev would handle "surface - 700 mb".

In the fortran code, reading and writing grib2 is done is a simple manner.

Code