program dtu2nc

! Rewrite the ASCII DTU MSS and ERR grids to NetCDF.
! In addition:
! - Remove the duplication at the longitude edges
! - Add both MSS and ERR grids in a single NetCDF file
!
! Example syntax:
! zcat DTU18MSS_1min.mss.gz DTU18ERR_1min.err.gz | dtu2nc DTU18MSS_1min.nc
!-
! 13-Feb-2019 - Created by Remko Scharroo (c) EUMETSAT
!-----------------------------------------------------------------------
use typesizes
use netcdf
use rads_netcdf
use rads_misc
use rads_geo
use rads_time
integer :: nx, ny, i
character(80) :: filename
integer(fourbyteint) :: ncid,varid(4),dimid(2)
integer(fourbyteint), parameter :: deflate_level=4
real(eightbytereal) :: x0, x1, dx, y0, y1, dy, y
real(eightbytereal), parameter :: scale_mss=1d-4, scale_err=1d-3
real(eightbytereal), allocatable :: mss(:,:), dh(:)
logical :: wgs

! Get the output filename

call getarg (1,filename)
wgs = index(filename,'WGS') > 0

! Read the MSS header

read (*,*) y0, y1, x0, x1, dx, dy
nx = nint((x1 - x0)/dx + 1)
ny = nint((y1 - y0)/dy + 1)

! Adjust boundaries to its proper values

dy = 180d0 / ny
dx = dy
x0 = 0d0 + dx / 2
x1 = 360d0 - dx / 2
y0 = -90d0 + dy / 2
y1 = 90d0 - dy / 2
allocate (mss(nx,ny),dh(ny))

! Compute the ellipsoid adjustment

do i = 1, ny
	y = y1 - (i-1) * dy
	dh(i) = dhellips (1, y)
enddo

! Write out the meta data for the new NetCDF grid

call getarg(1,filename)
call nfs(nf90_create(filename,nf90_write+nf90_nofill+nf90_hdf5,ncid))
call nfs(nf90_put_att(ncid,nf90_global,'Conventions','CF-1.7'))
call nfs(nf90_put_att(ncid,nf90_global,'title',filename(:5) // ' mean sea surface'))
call nfs(nf90_put_att(ncid,nf90_global,'Model',filename(:5)))
call nfs(nf90_put_att(ncid,nf90_global,'source',filename(:13)//'.mss.gz '//filename(:5)//'ERR'//filename(9:13)//'.err.gz'))
call nfs(nf90_put_att(ncid,nf90_global,'institution','Danish Technical University'))

if (wgs) then
	call nfs(nf90_put_att(ncid,nf90_global,'semi_major_axis',ae_wgs84))
	call nfs(nf90_put_att(ncid,nf90_global,'inverse_flattening',1d0/f_wgs84))
	call nfs(nf90_put_att(ncid,nf90_global,'history', &
		'dtu2nc: converted to WGS84 ellipsoid, joined MSS and ERR grids, and used internal compression'))
else
	call nfs(nf90_put_att(ncid,nf90_global,'semi_major_axis',ae_topex))
	call nfs(nf90_put_att(ncid,nf90_global,'inverse_flattening',1d0/f_topex))
	call nfs(nf90_put_att(ncid,nf90_global,'history', &
		'dtu2nc: joined MSS and ERR grids and used internal compression'))
endif
call nfs(nf90_put_att(ncid,nf90_global,'creation_time',timestamp('T')//'Z'))

call nf90_def_axis(ncid,'lon','longitude','degrees_east',nx-2,x0,x1,dimid(1),varid(1),'X','longitude')
call nf90_def_axis(ncid,'lat','latitude','degrees_north',ny,y0,y1,dimid(2),varid(2),'Y','latitude')

call nfs(nf90_def_var(ncid,'mean_sea_surf_sol2',nf90_int,dimid,varid(3)))
call nfs(nf90_put_att(ncid,varid(3),'long_name','mean sea surface height'))
call nfs(nf90_put_att(ncid,varid(3),'units','m'))
call nfs(nf90_put_att(ncid,varid(3),'scale_factor',scale_mss))

call nfs(nf90_def_var(ncid,'mean_sea_surf_sol2_acc',nf90_int2,dimid,varid(4)))
call nfs(nf90_put_att(ncid,varid(4),'long_name','error of mean sea surface height'))
call nfs(nf90_put_att(ncid,varid(4),'units','m'))
call nfs(nf90_put_att(ncid,varid(4),'scale_factor',scale_err))

! Switch on compression

do i = 1,4
	call nfs(nf90_def_var_deflate(ncid,varid(i),1,1,deflate_level))
enddo
call nfs(nf90_enddef(ncid))

! Write coordinate axes

call nf90_put_axis(ncid,varid(1))
call nf90_put_axis(ncid,varid(2))

! Read the MSS grid

read (*,*) mss

! Do some stats on the suspicious cells

!write (*,'(a,2f10.4)') "minimum and maximum difference cell 1 and nx-1 = ",minval(mss(1,:)-mss(nx-1,:)),maxval(mss(1,:)-mss(nx-1,:))
!write (*,'(a,2f10.4)') "minimum and maximum difference cell 2 and nx   = ",minval(mss(2,:)-mss(nx  ,:)),maxval(mss(2,:)-mss(nx  ,:))
!write (*,'(a,2f10.4)') "minimum and maximum difference cell 1 and surr = ", &
!	minval(mss(1,:)-mss(nx-2,:)/2-mss(nx,:)/2),maxval(mss(1,:)-mss(nx-2,:)/2-mss(nx,:)/2)
!do i = 1,ny
!	write (*,'(i6,8f10.4)') ny-i+1,y1-(i-1)*dy,mss(nx-3:nx,i),mss(3,i),mss(nx-2,i)-mss(nx,i),mss(nx-1,i)-mss(nx-2,i)/2-mss(nx,i)/2
!enddo

! If needed, change ellipsoid

if (wgs) then
	forall (i = 1:ny)
		mss(:,i) = mss(:,i) - dh(i)
	end forall
endif

! Write MSS grid, upside down

call nfs(nf90_put_var(ncid,varid(3),nint4(mss(2:nx-1,ny:1:-1)/scale_mss)))

! Skip the ERR header

read (*,*)

! Read and store the ERR grid

read (*,*) mss
call nfs(nf90_put_var(ncid,varid(4),nint2(mss(2:nx-1,ny:1:-1)/scale_err)))

! Update time stamp

call nfs(nf90_redef(ncid))
call nfs(nf90_put_att(ncid,nf90_global,'creation_time',timestamp('T')//'Z'))

! Close output files

call nfs(nf90_close(ncid))

deallocate (mss, dh)

end program dtu2nc