program convert

!-------------------------------------------------------------------------------------
! Read in original gmted2010 30 arc second dataset, then convert to a 
! format expected by gridgen_sfc and the fv3gfs orography code.
!
! The original data does not have data north of 84N.  This program fills
! in these areas with an elevation of zero to create a fully global dataset.
!
! The original data has data voids in Lake Superior and the Caspisan Sea.
! These voids are fixed.
!
! The data is also converted to two-byte integers and big endian.
!
! This dataset is shifted south/west by 1/2 arc second compared to the gtopo30 data.
!
! The latitude of the pixel center of the south most row is: -90.0 + 1/240 - 1/7200
! The longitude of the pixel center of the west most column is: -180 + 1/240 -1/7200
!-------------------------------------------------------------------------------------

 implicit none

 integer(kind=4) :: i, j, istat, istart, iend, jstart, jend
 integer(kind=8) :: record
 integer(kind=8) :: idim_orig, jdim_orig, idim_global, jdim_global

 integer(kind=2), allocatable :: orog_global(:,:)

 real(kind=4), allocatable    :: orog_orig(:,:)

 real(kind=8) :: lat11_global, lon11_global

 idim_orig=43200_8
 jdim_orig=20880_8
 record=idim_orig*jdim_orig*4_8

 open (11, file="./gmted2010.30sec.flt", iostat=istat, access='direct', form="unformatted", recl=record)
 print*,'istat on open: ',istat

 allocate(orog_orig(idim_orig,jdim_orig))
 read (11, rec=1, iostat=istat) orog_orig
 print*,'istat on read: ',istat
 close (11)
 print*,maxval(orog_orig),minval(orog_orig)

 idim_global=43200_8
 jdim_global=21600_8
 allocate(orog_global(idim_global,jdim_global))
 orog_global = 0

 do j = 1, jdim_orig
 do i = 1, idim_orig
   orog_global(i,j+jdim_global-jdim_orig) = nint(orog_orig(i,j))
 enddo
 enddo

!----------------------------------------------------------------------------
! Fix void (with zero elevation) in Lake Superior between 47/48N and 86/87W.
! Reset to its neighboring points (179 m).
!----------------------------------------------------------------------------

 lat11_global = 90.0_8 - (1.0_8/240.0_8) - (1.0_8/7200.0_8)
 lon11_global = -180.0_8 + (1.0_8/240.0_8) - (1.0_8/7200.0_8)

 istart = nint ( ((-87.0_8-lon11_global) / (1.0_8/120.0_8) ) +  1.0 )
 istart = istart - 5
 iend   = nint ( ((-86.0_8-lon11_global) / (1.0_8/120.0_8) ) +  1.0 )
 iend   = iend + 5
 jstart = nint ( (( 48.0_8-lat11_global) / (-1.0_8/120.0_8) ) +  1.0 )
 jstart = jstart - 5
 jend   = nint ( (( 47.0_8-lat11_global) / (-1.0_8/120.0_8) ) +  1.0 )
 jend   = jend + 5

 do j = jstart, jend
 do i = istart, iend
   if (orog_global(i,j) == 0) orog_global(i,j) = 179
 enddo
 enddo

!----------------------------------------------------------------------------
! Fix voids in Caspian Sea with zero elevation.  Reset to minus 27 meters.
!
! - 43/45N  48/50E
! - 42/43N  49/51E
! - 41/42N  50/52E
! - 40/41N  51/52E
! - 38/40N  50/52E
! - 38/39N  52/53E
! - 37/38N  51/53E
!----------------------------------------------------------------------------

 istart = nint ( ((48.0_8-lon11_global) / (1.0_8/120.0_8) ) +  1.0 )
 istart = istart - 5
 iend   = nint ( ((50.0_8-lon11_global) / (1.0_8/120.0_8) ) +  1.0 )
 iend   = iend + 5
 jstart = nint ( ((45.0_8-lat11_global) / (-1.0_8/120.0_8) ) +  1.0 )
 jstart = jstart - 5
 jend   = nint ( ((43.0_8-lat11_global) / (-1.0_8/120.0_8) ) +  1.0 )
 jend   = jend + 5

 do j = jstart, jend
 do i = istart, iend
   if (orog_global(i,j) == 0) orog_global(i,j) = -27
 enddo
 enddo

 istart = nint ( ((49.0_8-lon11_global) / (1.0_8/120.0_8) ) +  1.0 )
 istart = istart - 5
 iend   = nint ( ((51.0_8-lon11_global) / (1.0_8/120.0_8) ) +  1.0 )
 iend   = iend + 5
 jstart = nint ( ((43.0_8-lat11_global) / (-1.0_8/120.0_8) ) +  1.0 )
 jstart = jstart - 5
 jend   = nint ( ((42.0_8-lat11_global) / (-1.0_8/120.0_8) ) +  1.0 )
 jend   = jend + 5

 do j = jstart, jend
 do i = istart, iend
   if (orog_global(i,j) == 0) orog_global(i,j) = -27
 enddo
 enddo

 istart = nint ( ((50.0_8-lon11_global) / (1.0_8/120.0_8) ) +  1.0 )
 istart = istart - 5
 iend   = nint ( ((52.0_8-lon11_global) / (1.0_8/120.0_8) ) +  1.0 )
 iend   = iend + 5
 jstart = nint ( ((42.0_8-lat11_global) / (-1.0_8/120.0_8) ) +  1.0 )
 jstart = jstart - 5
 jend   = nint ( ((41.0_8-lat11_global) / (-1.0_8/120.0_8) ) +  1.0 )
 jend   = jend + 5

 do j = jstart, jend
 do i = istart, iend
   if (orog_global(i,j) == 0) orog_global(i,j) = -27
 enddo
 enddo

 istart = nint ( ((51.0_8-lon11_global) / (1.0_8/120.0_8) ) +  1.0 )
 istart = istart - 5
 iend   = nint ( ((52.0_8-lon11_global) / (1.0_8/120.0_8) ) +  1.0 )
 iend   = iend + 5
 jstart = nint ( ((41.0_8-lat11_global) / (-1.0_8/120.0_8) ) +  1.0 )
 jstart = jstart - 5
 jend   = nint ( ((40.0_8-lat11_global) / (-1.0_8/120.0_8) ) +  1.0 )
 jend   = jend + 5

 do j = jstart, jend
 do i = istart, iend
   if (orog_global(i,j) == 0) orog_global(i,j) = -27
 enddo
 enddo

 istart = nint ( ((50.0_8-lon11_global) / (1.0_8/120.0_8) ) +  1.0 )
 istart = istart - 5
 iend   = nint ( ((52.0_8-lon11_global) / (1.0_8/120.0_8) ) +  1.0 )
 iend   = iend + 5
 jstart = nint ( ((40.0_8-lat11_global) / (-1.0_8/120.0_8) ) +  1.0 )
 jstart = jstart - 5
 jend   = nint ( ((38.0_8-lat11_global) / (-1.0_8/120.0_8) ) +  1.0 )
 jend   = jend + 5

 do j = jstart, jend
 do i = istart, iend
   if (orog_global(i,j) == 0) orog_global(i,j) = -27
 enddo
 enddo

 istart = nint ( ((52.0_8-lon11_global) / (1.0_8/120.0_8) ) +  1.0 )
 istart = istart - 5
 iend   = nint ( ((53.0_8-lon11_global) / (1.0_8/120.0_8) ) +  1.0 )
 iend   = iend + 5
 jstart = nint ( ((39.0_8-lat11_global) / (-1.0_8/120.0_8) ) +  1.0 )
 jstart = jstart - 5
 jend   = nint ( ((38.0_8-lat11_global) / (-1.0_8/120.0_8) ) +  1.0 )
 jend   = jend + 5

 do j = jstart, jend
 do i = istart, iend
   if (orog_global(i,j) == 0) orog_global(i,j) = -27
 enddo
 enddo

 istart = nint ( ((51.0_8-lon11_global) / (1.0_8/120.0_8) ) +  1.0 )
 istart = istart - 5
 iend   = nint ( ((53.0_8-lon11_global) / (1.0_8/120.0_8) ) +  1.0 )
 iend   = iend + 5
 jstart = nint ( ((38.0_8-lat11_global) / (-1.0_8/120.0_8) ) +  1.0 )
 jstart = jstart - 5
 jend   = nint ( ((37.0_8-lat11_global) / (-1.0_8/120.0_8) ) +  1.0 )
 jend   = jend + 5

 do j = jstart, jend
 do i = istart, iend
   if (orog_global(i,j) == 0) orog_global(i,j) = -27
 enddo
 enddo

!----------------------------------------------------------------------------
! Write out final data.
!----------------------------------------------------------------------------

 record=idim_global*jdim_global*2_8
 open (12, file="./gmted2010.30sec.int", iostat=istat, access='direct', form="unformatted", recl=record, convert='big_endian')
 print*,'istat on open: ',istat
 write(12, rec=1) orog_global
 print*,'istat on write: ',istat
 close(12)

 stop
 end program convert