PROGRAM PlayWithPrecip

! Written by Kate T-C
! March 9, 2009  For the Fortran short course
! This program opens a netCDF file from NCEP reanalysis with 5 years of 
!     precipitation rates, reads in the precip rates, converts to mm/day,
!     and then calculates the zonal average and prints it to the terminal.

use netcdf
implicit none

character(len=80) :: infile = 'NCEP.Precip.0100-1204.nc'
integer :: ncid, err, PvarID=0, i, j
integer, parameter :: lon = 192, lat=92, time=1827
real, dimension(lon,lat,time) :: precip = 0.0
real, dimension(lon,lat) :: precipTA = 0.0
real, dimension(lat) :: precipZA = 0.0

    err = nf90_open(infile, NF90_NOWRITE, ncid) 
    print *, 'File Opened'
    
    err = nf90_inq_varid(ncid, 'prate', PvarID)
    print *, 'Got Var ID: ', PvarID
    
    err = nf90_get_var(ncid,PvarID,precip)
    print *, 'Precip loaded'

	! Convert to Real using scale and offset
	precip = precip * 1.0000000116861e-07 + 0.00327649991959333
	print *, 'precip max:', maxVal(precip), '  precip min:', minVal(precip)
	
	! Convert from Kg/m**2/s to mm/day
	precip = precip * 3600 * 24
	
	print *, 'precip max:', maxVal(precip), '  precip min:', minVal(precip)
	
	err = nf90_close(ncid)
	print *, 'File closed'
	
	! Do a zonal average
	do i=1,lat   ! average longitudes for each latitude
	  do j=1,lon    ! average time at each longitude
	    precipTA(j,i) = sum(precip(j,i,:))/time
	  enddo
	  	precipZA(i) = sum(precipTA(:,i))/lon
	enddo
	
	print *, 'Precip Zonal Average:'
	print *, precipZA

END PROGRAM PlayWithPrecip