program africa_rd 
use netcdf
implicit none

! This is the name of the data file we will read.
character (len = *), parameter :: IFILE_NAME = "africa.nc" 
character (len = *), parameter :: OFILE_NAME = "hydraulic.nc"
integer :: incid, oncid

! We are reading 2D data, a 12 x 6 lon-lat grid. 
integer, parameter :: INDIMS = 4, INRECS = 2
integer, parameter :: INLATS = 97, INLONS = 97, INLEVS = 18 
character (len = *), parameter :: ITIME_NAME = "time"
character (len = *), parameter :: ILAT_NAME = "iy" 
character (len = *), parameter :: ILON_NAME = "jx" 
character (len = *), parameter :: ILEV_NAME = "kz"

! We are writing 3D data, lon lat lev dimension.

integer, parameter :: ONDIMS = 3
integer, parameter :: ONLAT = 97, ONLON = 97, ONLEV = 18
character (len = *), parameter :: ONLAT_NAME = "lat"
character (len = *), parameter :: ONLON_NAME = "lon"
character (len = *), parameter :: ONLEV_NAME = "nlevgrnd"

! pres_dimid : output parameter's dimension id

integer :: ilat_dimid, ilon_dimid, ilev_dimid, itime_dimid, pres_dimid, irec_dimid
integer :: istart(INDIMS), icount(INDIMS)

! Output data dimension ids.
integer :: onlat_dimid, onlon_dimid, onlev_dimid
integer :: ostart(ONDIMS), ocount(ONDIMS)

! For the lat lon coordinate netCDF variables. 
real :: ilats(INLATS), ilons(INLONS) 
integer :: ilat_varid, ilon_varid 


! We will read surface temperature and pressure fields. 
character (len = *), parameter :: ZWND_NAME = "ua" 
character (len = *), parameter :: MWND_NAME = "va" 
character (len = *), parameter :: ATMP_NAME = "ta"
character (len = *), parameter :: PRES_NAME = "ps"
integer :: zwnd_varid, mwnd_varid, atmp_varid , pres_varid
integer, dimension(:), allocatable :: idimids

! pres_varid is output variable id


! To check the units attributes.
character (len = *), parameter :: IUNITS = "units" 
character (len = *), parameter :: ZWND_UNITS = "m s-1" 
character (len = *), parameter :: MWND_UNITS = "m s-1"
character (len = *), parameter :: ATMP_UNITS = "K"
character (len = *), parameter :: PRES_UNITS = "HPa" ! Output's unit attribute
character (len = *), parameter :: ILAT_UNITS = "degrees_north" 
character (len = *), parameter :: ILON_UNITS = "degrees_east" 

integer, parameter :: MAX_ATT_LEN = 80
integer :: att_len
character*(MAX_ATT_LEN) :: zwnd_units_in, mwnd_units_in, atmp_units_in, pres_units_in  ! pres_units_in is output parameter
character*(MAX_ATT_LEN) :: ilat_units_in, ilon_units_in

! Read the data into these arrays.
real :: zwnd_in(INLONS, INLATS, INLEVS)
real :: mwnd_in(INLONS, INLATS, INLEVS)
real :: atmp_in(INLONS, INLATS, INLEVS)
real :: pres_out(INLONS, INLATS, INLEVS)
real :: ilon_in(INLONS)
real :: ilat_in(INLATS)



!integer :: ndims_in, nvars_in, ngatts_in, unlimdimid_in

! Loop indices 
integer :: ilat, ilon, ilev, i, irec

! Open the file.
call check( nf90_open(IFILE_NAME, nf90_nowrite, incid) )

! Get the varids of the latitude and longitude coordinate variables. 
call check( nf90_inq_varid(incid, ILAT_NAME, ilat_varid) )
call check( nf90_inq_varid(incid, ILON_NAME, ilon_varid) )

! Read the latitude and longitude data. 

call check( nf90_get_var(incid, ilat_varid, ilats) )
call check( nf90_get_var(incid, ilon_varid, ilons) )

! Get the varids of the pressure and temperature netCDF variables. 

call check( nf90_inq_varid(incid, ZWND_NAME, zwnd_varid) )
call check( nf90_inq_varid(incid, ATMP_NAME, atmp_varid) )
call check( nf90_inq_varid(incid, MWND_NAME, mwnd_varid) )

icount = (/ INLONS, INLATS, INLEVS, 1  /)
istart = (/ 1, 1, 1, 1 /)


! Read the zonal wind, meridional wind  and air 
!temperature data from the file.

do irec = 1, INRECS
   istart(4) = irec


call check( nf90_get_var(incid, zwnd_varid, zwnd_in, istart, icount) )
call check( nf90_get_var(incid, mwnd_varid, mwnd_in, istart, icount) )
call check( nf90_get_var(incid, atmp_varid, atmp_in, istart, icount) )
call check( nf90_get_var(incid, ilon_varid, ilon_in, istart, icount) )
call check( nf90_get_var(incid, ilat_varid, ilat_in, istart, icount) )

i = 0
do ilev = 1, INLEVS
   do ilat = 1, INLATS
      do ilon = 1, INLONS
!          print *, ilon, ilat, ilev, zwnd_in(ilon, ilat, ilev), mwnd_in(ilon,ilat,ilev), atmp_in(ilon,ilat,ilev)
	pres_out(ilon, ilat, ilev) = zwnd_in(ilon, ilat, ilev) * mwnd_in(ilon,ilat,ilev) * atmp_in(ilon,ilat,ilev)
	print *, ilon_in(ilon), ilat_in(ilat), pres_out(ilon, ilat, ilev)
      end do
   end do
end do

end do

call check(nf90_create(OFILE_NAME, nf90_clobber, incid))

call check(nf90_def_dim(incid, ILAT_NAME, INLATS, ilat_dimid))
call check(nf90_def_dim(incid, ILON_NAME, INLONS, ilon_dimid))

call check(nf90_def_var(incid, ILAT_NAME, NF90_REAL, ilat_dimid, ilat_varid))
call check(nf90_def_var(incid, ILON_NAME, NF90_REAL, ilon_dimid, ilon_varid))

call check(nf90_put_att(incid, ilat_varid, IUNITS, ILAT_UNITS) )
call check(nf90_put_att(incid, ilon_varid, IUNITS, ILON_UNITS) )

idimids = (/ ilon_dimid, ilat_dimid, ilev_dimid /)

call check(nf90_def_var(incid, PRES_NAME, NF90_REAL, idimids, pres_varid))

call check(nf90_put_att(incid, pres_varid, IUNITS, PRES_UNITS) )

call check( nf90_enddef(incid) )

call check(nf90_put_var(incid, ilat_varid, ilat))
call check(nf90_put_var(incid, ilon_varid, ilon))


icount = (/INLONS, INLATS, INLEVS, 1 /)
istart = (/ 1, 1, 1, 1 /)

do irec = 1, INRECS
   istart(4) = irec
      call check(nf90_put_var(incid, pres_varid, pres_out, istart, icount) )    
end do

!Close the file. This frees up any internal netCDF resources 
!associated with the file.

call check( nf90_close(incid) )

! If we got this far, everything worked as expected. Yipee!
print *,"*** SUCCESS creating", OFILE_NAME," ! "

contains

subroutine check(status)
integer, intent ( in) :: status

if(status /= nf90_noerr) then
print *, trim(nf90_strerror(status))
stop 2
end if
end subroutine check

end program africa_rd