diff --git a/cmake/mpiexec.hercules b/cmake/mpiexec.hercules
index 332b33e29..23bec1047 100755
--- a/cmake/mpiexec.hercules
+++ b/cmake/mpiexec.hercules
@@ -6,7 +6,7 @@
 #   $2+ - Executable and its arguments
 #
 
-ACCOUNT=
+ACCOUNT=fv3-cpu
 QOS=debug
 
 NP=$1
diff --git a/reg_tests/grid_gen/c96.uniform.sh b/reg_tests/grid_gen/c96.uniform.sh
index 2be0ba3ca..f5e29aeee 100755
--- a/reg_tests/grid_gen/c96.uniform.sh
+++ b/reg_tests/grid_gen/c96.uniform.sh
@@ -17,6 +17,7 @@ export add_lake=true
 export lake_data_srce=MODISP_GLDBV3
 export lake_cutoff=0.50
 export binary_lake=1
+export ocn=500
 
 NCCMP=${NCCMP:-$(which nccmp)}
 
@@ -41,7 +42,7 @@ echo "Ending at: " `date`
 # Compare output to baseline set of data.
 #-----------------------------------------------------------------------------
 
-cd $out_dir/C96
+cd $out_dir/C96.mx500
 
 test_failed=0
 for files in *tile*.nc ./sfc/*tile*.nc
diff --git a/sorc/ocean_merge.fd/CMakeLists.txt b/sorc/ocean_merge.fd/CMakeLists.txt
index 661060a19..21b26b280 100644
--- a/sorc/ocean_merge.fd/CMakeLists.txt
+++ b/sorc/ocean_merge.fd/CMakeLists.txt
@@ -1,6 +1,5 @@
-list(APPEND fortran_src
-  merge_lake_ocnmsk.f90
-)
+set(lib_src read_write.F90 merge.F90 utils.F90 namelist.F90)
+set(exe_src merge_lake_ocnmsk.F90)
 
 if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel|IntelLLVM)$")
   set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -r8 -i4 -convert big_endian")
@@ -12,12 +11,17 @@ elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU)$")
 endif()
 
 set(exe_name ocean_merge)
-add_executable(${exe_name} ${fortran_src})
+
+add_library(om_lib STATIC ${lib_src})
+add_executable(${exe_name} ${exe_src})
+
 target_link_libraries(
-  ${exe_name}
-  
+  om_lib
+  PUBLIC
   NetCDF::NetCDF_Fortran)
 
+target_link_libraries(${exe_name} PRIVATE om_lib)
+
 install(TARGETS ${exe_name})
 
 # If doxygen documentation we enabled, build it.
diff --git a/sorc/ocean_merge.fd/merge.F90 b/sorc/ocean_merge.fd/merge.F90
new file mode 100644
index 000000000..36c007436
--- /dev/null
+++ b/sorc/ocean_merge.fd/merge.F90
@@ -0,0 +1,69 @@
+!> Determine the water mask by merging the lake mask with the mapped ocean
+!! mask from MOM6. Both are on the FV3 grid. During merge, the ocean mask
+!! dominates the lake mask if there is a conflict. After the merge, the remaining
+!! non-water fraction is the land fraction.
+!! 
+!! @param[in] lon The "east/west" dimension of the model grid.
+!! @param[in] lat The "north/south" dimension of the model grid.
+!! @param[in] binary_lake When '1', lake fraction is either 0 or 1. Otherwise, it is a fraction.
+!! @param[in] lat2d Latitude of the model grid points.
+!! @param[in] ocn_frac Fraction of the grid point that is ocean.
+!! @param[inout] lake_frac Fraction of the grid point that is lake.
+!! @param[inout] lake_depth Lake depth in meters.
+!! @param[out] land_frac Fraction of the grid point that is land.
+!! @param[out] slmsk Land/sea mask. '0' if less than 50% land. Otherwise, '1'.
+!!
+!! @author Shan Sun
+!! @author Rahul Mahajan
+!! @author Sanath Kumar
+ subroutine merge(lon, lat, binary_lake, lat2d, ocn_frac, &
+                  lake_frac, lake_depth, land_frac, slmsk)
+
+ implicit none
+
+ integer, intent(in)        :: lon, lat, binary_lake
+
+ real, intent(in)           :: lat2d(lon,lat)
+ real, intent(in)           :: ocn_frac(lon,lat)
+ real, intent(inout)        :: lake_frac(lon,lat)
+ real, intent(inout)        :: lake_depth(lon,lat)
+ real, intent(out)          :: land_frac(lon,lat)
+ real, intent(out)          :: slmsk(lon,lat)
+
+ real, parameter            :: min_land=1.e-4, def_lakedp=10.
+
+ integer                    :: i, j, nodp_pt, lake_pt
+
+ nodp_pt=0
+ lake_pt=0
+
+ do i=1,lon
+ do j=1,lat
+   if (binary_lake.eq.1) lake_frac(i,j)=nint(lake_frac(i,j))     ! using integer lake_frac
+   if (lat2d(i,j).le.-60.) lake_frac(i,j)=0.             ! ignore lakes on Antarctica
+   land_frac(i,j)=1.-ocn_frac(i,j)
+   if (land_frac(i,j) <    min_land) land_frac(i,j)=0.   ! ignore land  < min_land
+   if (land_frac(i,j) > 1.-min_land) land_frac(i,j)=1.   ! ignore water < min_land
+   if (1.-land_frac(i,j) > 0.) lake_frac(i,j)=0.         ! ocn dominates
+
+   if (lake_frac(i,j) > 0.) then
+     lake_pt=lake_pt+1            ! calculating total lake points
+     if (binary_lake.eq.1) then
+       land_frac(i,j)=0.
+     else
+       land_frac(i,j)=1.-lake_frac(i,j)
+     end if
+     if (lake_depth(i,j) <= 0.) then
+       lake_depth(i,j)=def_lakedp ! set missing lake depth to default value
+       nodp_pt=nodp_pt+1          ! calculating total lake points without depth
+     end if
+   else
+     lake_depth(i,j)=0.
+   end if
+   slmsk(i,j) = nint(land_frac(i,j)) ! nint got the land pts correct
+ end do
+ end do
+
+ write(*,'(a,i8,a,i8,a)') 'Total lake point ',lake_pt,' where ',nodp_pt,' has no depth'
+
+ end subroutine merge
diff --git a/sorc/ocean_merge.fd/merge_lake_ocnmsk.F90 b/sorc/ocean_merge.fd/merge_lake_ocnmsk.F90
new file mode 100644
index 000000000..ee2bccf6e
--- /dev/null
+++ b/sorc/ocean_merge.fd/merge_lake_ocnmsk.F90
@@ -0,0 +1,60 @@
+!> @file
+!! @brief Determines the water mask by merging the lake mask with
+!! the mapped ocean mask from MOM6.
+!! @author Shan Sun
+!! @author Rahul Mahajan
+!! @author Sanath Kumar
+
+!> Determine the water mask by merging the lake mask with the mapped ocean
+!! mask from MOM6, both are on the FV3 grid. During merge, the ocean mask
+!! dominates the lake mask if there is a conflict. After the merge, the remaining
+!! non-water fraction is the land fraction.
+!! 
+!! @return 0 for success, error code otherwise.
+!! @author Shan Sun
+!! @author Rahul Mahajan
+!! @author Sanath Kumar
+program merge_lake_ocnmsk
+
+  implicit none
+
+  character(len=120) :: pth1
+  character(len=120) :: pth2,pth3
+  character(len=10)  :: atmres,ocnres
+
+  integer :: binary_lake
+  integer :: lat,lon,tile
+
+  real, allocatable :: lake_frac(:,:),lake_depth(:,:),land_frac(:,:),ocn_frac(:,:),slmsk(:,:),lat2d(:,:)
+
+  print*,"- BEGIN OCEAN MERGE PROGRAM."
+
+  call read_nml(pth1, pth2, atmres, ocnres, pth3,binary_lake)
+
+  do tile=1,6
+
+    call read_grid_dims(pth1, atmres, ocnres, tile, lon, lat)
+
+    if (tile==1) then
+      write(6,*) 'lat=',lat,'lon=',lon
+      allocate (lake_frac(lon,lat),lake_depth(lon,lat),land_frac(lon,lat), &
+                ocn_frac(lon,lat),slmsk(lon,lat),lat2d(lon,lat))
+    endif
+
+    call read_ocean_frac(pth1,atmres,ocnres,tile,lon,lat,ocn_frac)
+
+    call read_lake_mask(pth2,atmres,tile,lon,lat,lake_frac, &
+                        lake_depth,lat2d)
+
+    call merge(lon, lat, binary_lake, lat2d, ocn_frac, lake_frac, lake_depth, land_frac, slmsk)
+
+    call write_data(atmres,ocnres,pth3,tile,lon,lat,land_frac, &
+                    lake_frac,lake_depth,slmsk)
+
+  end do ! tile
+
+  deallocate (lake_frac,lake_depth,land_frac,ocn_frac,slmsk,lat2d)
+
+  print*,"- NORMAL TERMINATION."
+
+end program merge_lake_ocnmsk
diff --git a/sorc/ocean_merge.fd/merge_lake_ocnmsk.f90 b/sorc/ocean_merge.fd/merge_lake_ocnmsk.f90
deleted file mode 100644
index 9cba4088d..000000000
--- a/sorc/ocean_merge.fd/merge_lake_ocnmsk.f90
+++ /dev/null
@@ -1,184 +0,0 @@
-!> @file
-!! @brief Determines the water mask by merging the lake mask with
-!! the mapped ocean mask from MOM6.
-!! @author Shan Sun
-!! @author Rahul Mahajan
-!! @author Sanath Kumar
-
-!> Determine the water mask by merging the lake mask with the mapped ocean
-!! mask from MOM6, both are on the FV3 grid. During merge, the ocean mask
-!! dominates the lake mask if there is a conflict. After the merge, the remaining
-!! non-water fraction is the land fraction.
-!! 
-!! @return 0 for success, error code otherwise.
-!! @author Shan Sun
-!! @author Rahul Mahajan
-!! @author Sanath Kumar
-program merge_lake_ocnmsk
-  use netcdf
-
-  implicit none
-
-  character(len=120) :: pth1
-  character(len=120) :: pth2,pth3
-  character(len=10)  :: atmres,ocnres
-  real, parameter    :: min_land=1.e-4, def_lakedp=10.
- ! this variable is now renamed as binary_lake and is passed in from the name
- ! list
- ! logical, parameter :: int_lake=.true.  
- ! all instances of int_lake was changed to binary_lake  
-  integer :: binary_lake
-
-  character(len=250) :: flnm
-  integer :: ncid,ndims,nvars,natts,lat,lon,v1id,v2id,v3id,v4id,start(2),count(2),i,j,latid,lonid,ncid4, dims(2),tile,nodp_pt
-  integer :: lake_pt,vlat
-  real, allocatable :: lake_frac(:,:),lake_depth(:,:),land_frac(:,:),ocn_frac(:,:),slmsk(:,:),lat2d(:,:)
-
-
-  call read_nml(pth1, pth2, atmres, ocnres, pth3,binary_lake)
-  
-  print *, pth1
-  nodp_pt=0
-  lake_pt=0  
-  do tile=1,6
-    write(flnm,'(5a,i1,a)') trim(pth1),trim(atmres),'.',trim(ocnres),'.tile',tile,'.nc'
-    call handle_err (nf90_open (flnm, NF90_NOWRITE, ncid))
-    call handle_err (nf90_inquire (ncid, ndimensions=ndims, nvariables=nvars, nattributes=natts))
-    write(6,*) 'flnm_ocn=',flnm,' ncid=',ncid, ' ndims=',ndims, 'nvars=',nvars,' natts=',natts
-    call handle_err (nf90_inq_dimid (ncid, 'grid_xt', latid))  ! RM: lon is no longer in this file; try grid_xt
-    call handle_err (nf90_inq_dimid (ncid, 'grid_yt', lonid))  ! RM: lat is no longer in this file; try grid_tyt
-    call handle_err (nf90_inquire_dimension (ncid, latid, len=lat))
-    call handle_err (nf90_inquire_dimension (ncid, lonid, len=lon))
-    if (tile==1) then
-      write(6,*) 'lat=',lat,'lon=',lon
-      allocate (lake_frac(lon,lat),lake_depth(lon,lat),land_frac(lon,lat),ocn_frac(lon,lat),slmsk(lon,lat),lat2d(lon,lat))
-      start(1:2) = (/1,1/)
-      count(1:2) = (/lon,lat/)
-    end if
-    call handle_err (nf90_inq_varid(ncid, 'land_frac', v1id))
-    call handle_err (nf90_get_var (ncid, v1id, ocn_frac, start=start, count=count))
-
-   !write(flnm,'(3a,i1,a)') trim(pth2),trim(atmres),'_orog_data.tile',tile,'.nc'
-    write(flnm,'(4a,i1,a)') trim(pth2),'oro.',trim(atmres),'.tile',tile,'.nc'
-    print *,' flnm2=',trim(flnm)
-    call handle_err (nf90_open (flnm, NF90_NOWRITE, ncid))
-    call handle_err (nf90_inquire (ncid, ndimensions=ndims, nvariables=nvars, nattributes=natts))
-    write(6,*) 'flnm_lake=',flnm,' ncid=',ncid, ' ndims=',ndims, 'nvars=',nvars,' natts=',natts
-   !if (tile==1) then
-   !  call handle_err (nf90_inq_dimid (ncid, 'lat', latid))
-   !  call handle_err (nf90_inq_dimid (ncid, 'lon', lonid))
-   !  call handle_err (nf90_inquire_dimension (ncid, latid, len=lat))
-   !  call handle_err (nf90_inquire_dimension (ncid, lonid, len=lon))
-   !  write(6,*) 'lat=',lat,'lon=',lon
-   !  start(1:2) = (/1,1/)
-   !  count(1:2) = (/lon,lat/)
-   !end if
-    call handle_err (nf90_inq_varid(ncid, 'lake_frac', v2id))
-    call handle_err (nf90_inq_varid(ncid, 'lake_depth',v3id))
-    call handle_err (nf90_inq_varid(ncid, 'geolat'    ,vlat))
-    call handle_err (nf90_get_var (ncid, v2id, lake_frac, start=start, count=count))
-    call handle_err (nf90_get_var (ncid, v3id, lake_depth,start=start, count=count))
-    call handle_err (nf90_get_var (ncid, vlat, lat2d,     start=start, count=count))
-
-    do i=1,lon
-    do j=1,lat
-      if (binary_lake.eq.1) lake_frac(i,j)=nint(lake_frac(i,j))     ! using integer lake_frac
-      if (lat2d(i,j).le.-60.) lake_frac(i,j)=0.             ! ignore lakes on Antarctica
-      land_frac(i,j)=1.-ocn_frac(i,j)
-      if (land_frac(i,j) <    min_land) land_frac(i,j)=0.   ! ignore land  < min_land
-      if (land_frac(i,j) > 1.-min_land) land_frac(i,j)=1.   ! ignore water < min_land
-      if (1.-land_frac(i,j) > 0.) lake_frac(i,j)=0.         ! ocn dominates
-
-      if (lake_frac(i,j) > 0.) then
-        lake_pt=lake_pt+1            ! calculating total lake points
-        if (binary_lake.eq.1) then
-          land_frac(i,j)=0.
-        else
-          land_frac(i,j)=1.-lake_frac(i,j)
-        end if
-        if (lake_depth(i,j) <= 0.) then
-          lake_depth(i,j)=def_lakedp ! set missing lake depth to default value
-          nodp_pt=nodp_pt+1          ! calculating total lake points without depth
-        end if
-      else
-        lake_depth(i,j)=0.
-      end if
-!      slmsk(i,j) = ceiling(land_frac(i,j))! "ceiling is used for orog smoothing"
-       slmsk(i,j) = nint(land_frac(i,j)) ! nint got the land pts correct
-    end do
-    end do
-
-    write(flnm,'(4a,i1,a)') trim(atmres),'.',trim(ocnres),'.tile',tile,'.nc'
-    print *,'output=',trim(flnm)
-    call handle_err (nf90_create (path=trim(pth3)//trim(flnm), &
-    cmode=or(NF90_CLOBBER, NF90_64BIT_OFFSET), ncid=ncid4))   ! netcdf3
-
-    call handle_err (nf90_def_dim (ncid4,'lon', lon, dims(1)))
-    call handle_err (nf90_def_dim (ncid4,'lat', lat, dims(2)))
-    call handle_err (nf90_def_var (ncid4,'land_frac', nf90_float, dims(1:2), v1id))
-    call handle_err (nf90_def_var (ncid4,'lake_frac', nf90_float, dims(1:2), v2id))
-    call handle_err (nf90_def_var (ncid4,'lake_depth',nf90_float, dims(1:2), v3id))
-    call handle_err (nf90_def_var (ncid4,'slmsk',     nf90_float, dims(1:2), v4id))
-
-    call handle_err (nf90_enddef  (ncid4))
-
-    call handle_err (nf90_put_var (ncid4, v1id,land_frac))
-    call handle_err (nf90_put_var (ncid4, v2id,lake_frac))
-    call handle_err (nf90_put_var (ncid4, v3id,lake_depth))
-    call handle_err (nf90_put_var (ncid4, v4id,slmsk))
-    call handle_err (nf90_close(ncid4))
-
-!   do i=1,48
-!   do j=1,48
-!     write(*,'(2i4,4f6.1)') i,j,land_frac(i,j),lake_frac(i,j),lake_depth(i,j),slmsk(i,j)
-!   end do
-!   end do
-
-  end do ! tile
-  write(*,'(a,i8,a,i8,a)') 'total lake point ',lake_pt,' where ',nodp_pt,' has no depth'
-
-end program merge_lake_ocnmsk
-
-!> Handle netCDF errors.
-!!
-!! @param[in] ret NetCDF return code.
-!! @author Shan Sun
-subroutine handle_err (ret)
-  use netcdf
-  integer, intent(in) :: ret
-
-  if (ret /= NF90_NOERR) then
-    write(6,*) nf90_strerror (ret)
-    stop 999
-  end if
-end subroutine handle_err
-
-!> Read program namelist.
-!!
-!! @param[out] ocean_mask_dir Directory containing MOM6 ocean mask file.
-!! @param[out] lake_mask_dir Directory containing the lake mask file.
-!! @param[out] out_dir Directory where output file will be written.
-!! @param[out] atmres Atmosphere grid resolution.
-!! @param[out] ocnres Ocean grid resolution.
-!! @param[out] binary_lake or fractional lake 
-!! @author Rahul Mahajan
-!! @author Sanath Kumar
-subroutine read_nml(ocean_mask_dir, lake_mask_dir, atmres,ocnres,out_dir,binary_lake)
-
-  integer :: unit=7, io_status
-
-  character(len=120), intent(out) :: ocean_mask_dir
-  character(len=120), intent(out) :: lake_mask_dir
-  character(len=120), intent(out) :: out_dir
-  character(len=10),  intent(out) :: atmres,ocnres
-  integer, intent(out):: binary_lake
-
-  namelist/mask_nml/ocean_mask_dir, lake_mask_dir, atmres, ocnres,out_dir,binary_lake
-  open(unit=unit, file='input.nml', iostat=io_status )
-  read(unit,mask_nml, iostat=io_status )
-  close(unit)
-  if (io_status > 0) then
-        print *,'Error reading input.nml' 
-        call handle_err(-1)
-  end if      
-end subroutine read_nml
diff --git a/sorc/ocean_merge.fd/namelist.F90 b/sorc/ocean_merge.fd/namelist.F90
new file mode 100644
index 000000000..f66896eae
--- /dev/null
+++ b/sorc/ocean_merge.fd/namelist.F90
@@ -0,0 +1,33 @@
+!> Read program namelist.
+!!
+!! @param[out] ocean_mask_dir Directory containing MOM6 ocean mask file.
+!! @param[out] lake_mask_dir Directory containing the lake mask file.
+!! @param[out] out_dir Directory where output file will be written.
+!! @param[out] atmres Atmosphere grid resolution.
+!! @param[out] ocnres Ocean grid resolution.
+!! @param[out] binary_lake When '1', treat lake fraction as either 0 or 1. Otherwise, 
+!! it is a fraction.
+!!
+!! @author Rahul Mahajan
+!! @author Sanath Kumar
+subroutine read_nml(ocean_mask_dir, lake_mask_dir, atmres,ocnres,out_dir,binary_lake)
+
+  implicit none
+
+  integer :: unit=7, io_status
+
+  character(len=120), intent(out) :: ocean_mask_dir
+  character(len=120), intent(out) :: lake_mask_dir
+  character(len=120), intent(out) :: out_dir
+  character(len=10),  intent(out) :: atmres,ocnres
+  integer, intent(out):: binary_lake
+
+  namelist/mask_nml/ocean_mask_dir, lake_mask_dir, atmres, ocnres,out_dir,binary_lake
+  open(unit=unit, file='input.nml', iostat=io_status )
+  read(unit,mask_nml, iostat=io_status )
+  close(unit)
+  if (io_status > 0) then
+        print *,'FATAL ERROR reading input.nml' 
+        call handle_err(-1)
+  end if      
+end subroutine read_nml
diff --git a/sorc/ocean_merge.fd/read_write.F90 b/sorc/ocean_merge.fd/read_write.F90
new file mode 100644
index 000000000..b445c525b
--- /dev/null
+++ b/sorc/ocean_merge.fd/read_write.F90
@@ -0,0 +1,200 @@
+!> Read the grid dimensions from the MOM6 ocean mask NetCDF file.
+!!
+!! @param[in] pth1 Directory path to file.
+!! @param[in] atmres Atmospheric resolution.
+!! @param[in] ocnres Ocean resolution in decimal percent.
+!! @param[in] tile Tile number.
+!! @param[out] lon E/W dimension of tile.
+!! @param[out] lat N/S dimension of tile.
+!!
+!! @author Shan Sun
+!! @author Rahul Mahajan
+ subroutine read_grid_dims(pth1, atmres, ocnres, tile, lon, lat)
+
+ use netcdf
+
+ implicit none
+
+ character(len=*), intent(in)   :: pth1
+ character(len=*), intent(in)   :: atmres
+ character(len=*), intent(in)   :: ocnres
+
+ integer, intent(in)            :: tile
+ integer, intent(out)           :: lon, lat
+
+ character(len=250)             :: flnm
+
+ integer                        :: ncid, ndims, nvars, natts
+ integer                        :: latid, lonid
+
+ write(flnm,'(5a,i1,a)') trim(pth1),trim(atmres),'.',trim(ocnres),'.tile',tile,'.nc'
+
+ print*,'- READ GRID DIMESIONS FROM: ',trim(flnm)
+
+ call handle_err (nf90_open (flnm, NF90_NOWRITE, ncid))
+ call handle_err (nf90_inquire (ncid, ndimensions=ndims, nvariables=nvars, nattributes=natts))
+ call handle_err (nf90_inquire (ncid, ndimensions=ndims, nvariables=nvars, nattributes=natts))
+ call handle_err (nf90_inq_dimid (ncid, 'grid_xt', latid))  ! RM: lon is no longer in this file; try grid_xt
+ call handle_err (nf90_inq_dimid (ncid, 'grid_yt', lonid))  ! RM: lat is no longer in this file; try grid_yt
+ call handle_err (nf90_inquire_dimension (ncid, latid, len=lat))
+ call handle_err (nf90_inquire_dimension (ncid, lonid, len=lon))
+ call handle_err (nf90_close (ncid))
+
+ print*,'- DIMENSIONS ARE: ',lon, lat
+
+ end subroutine read_grid_dims
+
+!> Read the ocean fraction from the MOM6 ocean NetCDF file.
+!!
+!! @param[in] pth1 Directory path to file.
+!! @param[in] atmres Atmospheric resolution.
+!! @param[in] ocnres Ocean resolution.
+!! @param[in] tile Tile number.
+!! @param[in] lon E/W dimension of tile.
+!! @param[in] lat N/S dimension of tile.
+!! @param[out] ocn_frac ocean fraction in decimal percent.
+!!
+!! @author Shan Sun
+!! @author Rahul Mahajan
+ subroutine read_ocean_frac(pth1,atmres,ocnres,tile,lon,lat,ocn_frac)
+
+ use netcdf
+
+ implicit none
+
+ character(len=*), intent(in)   :: pth1
+ character(len=*), intent(in)   :: atmres
+ character(len=*), intent(in)   :: ocnres
+ 
+ integer,          intent(in)   :: lat
+ integer,          intent(in)   :: lon
+ integer,          intent(in)   :: tile
+
+ real, intent(out)              :: ocn_frac(lon,lat)
+
+ character(len=300)             :: flnm
+
+ integer                        :: ncid, v1id, start(2), count(2)
+
+ write(flnm,'(5a,i1,a)') trim(pth1),trim(atmres),'.',trim(ocnres),'.tile',tile,'.nc'
+
+ print*,'-READ OCEAN FRACTION FROM: ',trim(flnm)
+
+ start(1:2) = (/1,1/)
+ count(1:2) = (/lon,lat/)
+
+ call handle_err (nf90_open (flnm, NF90_NOWRITE, ncid))
+
+! The file record is named 'land_frac', but the data is ocean fraction.
+
+ call handle_err (nf90_inq_varid(ncid, 'land_frac', v1id))
+ call handle_err (nf90_get_var (ncid, v1id, ocn_frac, start=start, count=count))
+ call handle_err (nf90_close (ncid))
+
+ end subroutine read_ocean_frac
+
+!> Read lake fraction, lake depth and latitude from a NetCDF file.
+!!
+!! @param[in] pth2 Directory path to the file.
+!! @param[in] atmres Atmospheric resolution.
+!! @param[in] tile Tile number.
+!! @param[in] lon E/W dimension of tile.
+!! @param[in] lat N/S dimension of tile.
+!! @param[out] lake_frac Lake fraction in decimal percent.
+!! @param[out] lake_depth Lake depth in meters.
+!! @param[out] lat2d Latitude in degrees.
+!!
+!! @author Shan Sun
+!! @author Rahul Mahajan
+ subroutine read_lake_mask(pth2,atmres,tile,lon,lat,lake_frac, &
+                           lake_depth,lat2d)
+
+ use netcdf
+
+ implicit none
+
+ character(len=*), intent(in) :: pth2, atmres
+
+ integer, intent(in)          :: tile, lon, lat
+
+ real, intent(out)            :: lake_frac(lon,lat)
+ real, intent(out)            :: lake_depth(lon,lat)
+ real, intent(out)            :: lat2d(lon,lat)
+
+ character(len=250)           :: flnm
+
+ integer                      :: ncid, ndims, nvars, natts
+ integer                      :: v2id, v3id, vlat
+ integer                      :: start(2), count(2)
+
+ write(flnm,'(4a,i1,a)') trim(pth2),'oro.',trim(atmres),'.tile',tile,'.nc'
+ print *,'- READ LAKE DEPTH, FRACTION AND LATITUDE FROM: ',trim(flnm)
+ call handle_err (nf90_open (flnm, NF90_NOWRITE, ncid))
+ call handle_err (nf90_inquire (ncid, ndimensions=ndims, nvariables=nvars, nattributes=natts))
+ call handle_err (nf90_inq_varid(ncid, 'lake_frac', v2id))
+ call handle_err (nf90_inq_varid(ncid, 'lake_depth',v3id))
+ call handle_err (nf90_inq_varid(ncid, 'geolat'    ,vlat))
+ start(1:2) = (/1,1/)
+ count(1:2) = (/lon,lat/)
+ call handle_err (nf90_get_var (ncid, v2id, lake_frac, start=start, count=count))
+ call handle_err (nf90_get_var (ncid, v3id, lake_depth,start=start, count=count))
+ call handle_err (nf90_get_var (ncid, vlat, lat2d,     start=start, count=count))
+ call handle_err (nf90_close (ncid))
+
+ end subroutine read_lake_mask
+!> Write the merged data to a NetCDF file. Each tile is in
+!! its own file.
+!!
+!! @param[in] atmres Atmospheric resolution.
+!! @param[in] ocnres Ocean resolution.
+!! @param[in] pth3 Directory path to output file.
+!! @param[in] tile Tile number.
+!! @param[in] lon E/W dimension of tile.
+!! @param[in] lat N/S dimension of tile.
+!! @param[in] land_frac Land fraction in decimal percent.
+!! @param[in] lake_frac Lake fraction.
+!! @param[in] lake_depth Lake depth in meters.
+!! @param[in] slmsk Land/sea mask - 0-non-land; 1-land.
+!!
+!! @author Shan Sun
+!! @author Rahul Mahajan
+!! @author Sanath Kumar
+ subroutine write_data(atmres,ocnres,pth3,tile,lon,lat,land_frac, &
+                       lake_frac,lake_depth,slmsk)
+
+ use netcdf
+
+ implicit none
+
+ character(len=*), intent(in)       :: atmres, ocnres, pth3
+ 
+ integer, intent(in)                :: tile, lon, lat
+
+ real, intent(in)                   :: land_frac(lon,lat), lake_frac(lon,lat)
+ real, intent(in)                   :: lake_depth(lon,lat), slmsk(lon,lat)
+
+ character(len=250) :: flnm
+
+ integer            :: ncid4, dims(2), v1id, v2id, v3id, v4id
+ 
+ write(flnm,'(4a,i1,a)') trim(atmres),'.',trim(ocnres),'.tile',tile,'.nc'
+ print *,'- OUTPUT DATA TO FILE: ',trim(flnm)
+ call handle_err (nf90_create (path=trim(pth3)//trim(flnm), &
+    cmode=or(NF90_CLOBBER, NF90_64BIT_OFFSET), ncid=ncid4))   ! netcdf3
+
+ call handle_err (nf90_def_dim (ncid4,'lon', lon, dims(1)))
+ call handle_err (nf90_def_dim (ncid4,'lat', lat, dims(2)))
+ call handle_err (nf90_def_var (ncid4,'land_frac', nf90_float, dims(1:2), v1id))
+ call handle_err (nf90_def_var (ncid4,'lake_frac', nf90_float, dims(1:2), v2id))
+ call handle_err (nf90_def_var (ncid4,'lake_depth',nf90_float, dims(1:2), v3id))
+ call handle_err (nf90_def_var (ncid4,'slmsk',     nf90_float, dims(1:2), v4id))
+
+ call handle_err (nf90_enddef  (ncid4))
+
+ call handle_err (nf90_put_var (ncid4, v1id,land_frac))
+ call handle_err (nf90_put_var (ncid4, v2id,lake_frac))
+ call handle_err (nf90_put_var (ncid4, v3id,lake_depth))
+ call handle_err (nf90_put_var (ncid4, v4id,slmsk))
+ call handle_err (nf90_close(ncid4))
+
+ end subroutine write_data
diff --git a/sorc/ocean_merge.fd/utils.F90 b/sorc/ocean_merge.fd/utils.F90
new file mode 100644
index 000000000..22e2fadbd
--- /dev/null
+++ b/sorc/ocean_merge.fd/utils.F90
@@ -0,0 +1,16 @@
+!> Check NetCDF return code. If an error is indicated,
+!! stop program.
+!!
+!! @param[in] ret NetCDF return code.
+!! @author Shan Sun
+subroutine handle_err (ret)
+  use netcdf
+  implicit none
+  integer, intent(in) :: ret
+
+  if (ret /= NF90_NOERR) then
+    write(6,*) '- FATAL ERROR.'
+    write(6,*) nf90_strerror (ret)
+    stop 999
+  end if
+end subroutine handle_err
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 9a5850084..4299cb7cb 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -42,4 +42,5 @@ add_subdirectory(sfc_climo_gen)
 add_subdirectory(cpld_gridgen)
 add_subdirectory(emcsfc_snow2mdl)
 add_subdirectory(ocnice_prep)
+add_subdirectory(ocean_merge)
 add_subdirectory(orog)
diff --git a/tests/ocean_merge/CMakeLists.txt b/tests/ocean_merge/CMakeLists.txt
new file mode 100644
index 000000000..2aaa785b9
--- /dev/null
+++ b/tests/ocean_merge/CMakeLists.txt
@@ -0,0 +1,39 @@
+# This is the cmake build file for the tests directory of the
+# UFS_UTILS project.
+#
+# George Gayno
+
+if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel|IntelLLVM)$")
+  set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -r8 -assume byterecl")
+elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU)$")
+  set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -ffree-line-length-0 -fdefault-real-8")
+endif()
+
+# Copy necessary test files from the source data directory to the
+# build data directory.
+
+execute_process( COMMAND ${CMAKE_COMMAND} -E copy
+    ${CMAKE_CURRENT_SOURCE_DIR}/data/input.nml ${CMAKE_CURRENT_BINARY_DIR}/input.nml)
+
+execute_process( COMMAND ${CMAKE_COMMAND} -E copy
+    ${CMAKE_CURRENT_SOURCE_DIR}/data/oro.C48.tile1.nc ${CMAKE_CURRENT_BINARY_DIR}/oro.C48.tile1.nc)
+
+execute_process( COMMAND ${CMAKE_COMMAND} -E copy
+    ${CMAKE_CURRENT_SOURCE_DIR}/data/C48.mx500.tile1.nc ${CMAKE_CURRENT_BINARY_DIR}/C48.mx500.tile1.nc)
+
+add_executable(ftst_read_nml ftst_read_nml.F90)
+target_link_libraries(ftst_read_nml om_lib)
+
+add_executable(ftst_merge ftst_merge.F90)
+target_link_libraries(ftst_merge om_lib)
+
+add_executable(ftst_read_lake_mask ftst_read_lake_mask.F90)
+target_link_libraries(ftst_read_lake_mask om_lib)
+
+add_executable(ftst_read_ocean_frac ftst_read_ocean_frac.F90)
+target_link_libraries(ftst_read_ocean_frac om_lib)
+
+add_test(NAME ocean_merge-ftst_read_nml COMMAND ftst_read_nml)
+add_test(NAME ocean_merge-ftst_merge COMMAND ftst_merge)
+add_test(NAME ocean_merge-ftst_read_lake_mask COMMAND ftst_read_lake_mask)
+add_test(NAME ocean_merge-ftst_read_ocean_frac COMMAND ftst_read_ocean_frac)
diff --git a/tests/ocean_merge/data/C48.mx500.tile1.nc b/tests/ocean_merge/data/C48.mx500.tile1.nc
new file mode 100644
index 000000000..d3999075f
Binary files /dev/null and b/tests/ocean_merge/data/C48.mx500.tile1.nc differ
diff --git a/tests/ocean_merge/data/input.nml b/tests/ocean_merge/data/input.nml
new file mode 100644
index 000000000..ab86929b0
--- /dev/null
+++ b/tests/ocean_merge/data/input.nml
@@ -0,0 +1,8 @@
+     &mask_nml
+     ocean_mask_dir="/ocean/mask/dir"
+     ocnres="mx500"
+     lake_mask_dir="/lake/mask/dir"
+     atmres="C96"
+     out_dir="/out/dir"
+     binary_lake=1
+    /
diff --git a/tests/ocean_merge/data/oro.C48.tile1.nc b/tests/ocean_merge/data/oro.C48.tile1.nc
new file mode 100644
index 000000000..64366c355
Binary files /dev/null and b/tests/ocean_merge/data/oro.C48.tile1.nc differ
diff --git a/tests/ocean_merge/ftst_merge.F90 b/tests/ocean_merge/ftst_merge.F90
new file mode 100644
index 000000000..6d9c40728
--- /dev/null
+++ b/tests/ocean_merge/ftst_merge.F90
@@ -0,0 +1,207 @@
+! Unit test for the merge routine.
+! 
+! Test several combinations of lake and ocean attributes
+! and check the 'merged' value for correctness.
+ program ftst_merge
+
+ implicit none
+
+ integer, parameter   :: lon = 1
+ integer, parameter   :: lat = 1
+
+ integer              :: binary_lake
+
+ real, parameter      :: epsilon=0.00001
+ real                 :: lat2d(lon,lat)
+ real                 :: ocn_frac(lon,lat)
+ real                 :: lake_frac(lon,lat)
+ real                 :: lake_depth(lon,lat)
+ real                 :: land_frac(lon,lat)
+ real                 :: slmsk(lon,lat)
+
+ print*,"Begin test of merge routine"
+
+! Test point 1
+! Some ocean. No lake. Ocean info retained.
+
+ print*,'Test point 1.'
+
+ binary_lake = 0         ! Keep fractional lake.
+ lat2d(1,1) = 30.0       ! Point at 30N.
+ ocn_frac(1,1) = .75     ! .75 ocean/.25 land
+ lake_frac(1,1) = 0.0    ! No lake.
+ lake_depth(1,1) = 0.0   ! No lake.
+ land_frac(1,1) = -99.   ! Based on ocn_frac, should be .25.
+ slmsk(1,1) = -99.       ! Should be zero (nint of land_frac)
+
+ call merge(lon, lat, binary_lake, lat2d, ocn_frac, &
+            lake_frac, lake_depth, land_frac, slmsk)
+
+ if ( abs(lake_frac(1,1) - 0.) > epsilon) stop 2
+ if ( abs(lake_depth(1,1) - 0.) > epsilon) stop 2
+ if ( abs(land_frac(1,1) - .25) > epsilon) stop 2
+ if ( abs(slmsk(1,1) - 0.) > epsilon) stop 2
+
+! Test point 2
+! No ocean. Some lake. Lake info retained.
+
+ print*,'Test point 2.'
+
+ binary_lake = 0          ! Keep fractional lake.
+ lat2d(1,1) = 30.0        ! Point at 30N.
+ ocn_frac(1,1) = 0.0      ! 0 ocean/1 land.
+ lake_frac(1,1) = .75     ! Some lake.
+ lake_depth(1,1) = 100.0  ! Lake depth.
+ land_frac(1,1) = -99.    ! Should be .25
+ slmsk(1,1) = -99.        ! Should be zero (nint of land_frac).
+
+ call merge(lon, lat, binary_lake, lat2d, ocn_frac, &
+            lake_frac, lake_depth, land_frac, slmsk)
+
+ if ( abs(lake_frac(1,1) - .75) > epsilon) stop 4
+ if ( abs(lake_depth(1,1) - 100.) > epsilon) stop 4
+ if ( abs(land_frac(1,1) - .25) > epsilon) stop 4
+ if ( abs(slmsk(1,1) - 0.) > epsilon) stop 4
+
+! Test point 3
+! Some lake and some ocean. Ocean should dominate. Lake
+! should be removed.
+
+ print*,'Test point 3.'
+
+ binary_lake = 0           ! Keep fractional lake.
+ lat2d(1,1) = 30.0         ! Point at 30N.
+ ocn_frac(1,1) = .45       ! Some ocean.
+ lake_frac(1,1) = .75      ! Some lake.
+ lake_depth(1,1) = 100.0   ! Lake depth
+ land_frac(1,1) = -99.     ! Should be 1 minus ocn_frac
+ slmsk(1,1) = -99.         ! Should be 1 (nint of land_frac).
+
+ call merge(lon, lat, binary_lake, lat2d, ocn_frac, &
+            lake_frac, lake_depth, land_frac, slmsk)
+
+ if ( abs(lake_frac(1,1) - 0.) > epsilon) stop 6
+ if ( abs(lake_depth(1,1) - 0.) > epsilon) stop 6
+ if ( abs(land_frac(1,1) - .55) > epsilon) stop 6
+ if ( abs(slmsk(1,1) - 1.) > epsilon) stop 6
+
+! Test point 4
+! No ocean. Some lake. Lake will be retainted. Lake has 
+! a missing depth that should be given a default value.
+
+ print*,'Test point 4.'
+
+ binary_lake = 0           ! Keep fractional lake.
+ lat2d(1,1) = 30.0         ! Point at 30N.
+ ocn_frac(1,1) = 0.0       ! No ocean.
+ lake_frac(1,1) = .75      ! Some lake.
+ lake_depth(1,1) = -9.     ! Lake depth is missing.
+ land_frac(1,1) = -99.     ! Should be 1 minus lake_frac
+ slmsk(1,1) = -99.         ! Should be 0 (nint of land_frac).
+
+ call merge(lon, lat, binary_lake, lat2d, ocn_frac, &
+            lake_frac, lake_depth, land_frac, slmsk)
+
+ if ( abs(lake_frac(1,1) - .75) > epsilon) stop 8
+ if ( abs(lake_depth(1,1) - 10.) > epsilon) stop 8
+ if ( abs(land_frac(1,1) - .25) > epsilon) stop 8
+ if ( abs(slmsk(1,1) - 0.) > epsilon) stop 8
+
+! Test point 5
+! Some ocean (but very small percentage). No lake. 
+! The ocean % is below the min threshold, so it will be
+! removed and the point becomes all land.
+
+ print*,'Test point 5.'
+
+ binary_lake = 0          ! Keep fractional lake.
+ lat2d(1,1) = 30.0        ! Point at 30N.
+ ocn_frac(1,1) = 1.e-6    ! Very small % of ocean.
+ lake_frac(1,1) = 0.0     ! No lake.
+ lake_depth(1,1) = 0.0    ! No lake.
+ land_frac(1,1) = -99.    ! Should be 1.
+ slmsk(1,1) = -99.        ! Should be 1.
+
+ call merge(lon, lat, binary_lake, lat2d, ocn_frac, &
+            lake_frac, lake_depth, land_frac, slmsk)
+
+ if ( abs(lake_frac(1,1) - 0.) > epsilon) stop 10
+ if ( abs(lake_depth(1,1) - 0.) > epsilon) stop 10
+ if ( abs(land_frac(1,1) - 1.) > epsilon) stop 10
+ if ( abs(slmsk(1,1) - 1.) > epsilon) stop 10
+
+! Test point 6
+! Almost all ocean. Some lake. Since the ocean %
+! is very close to 1, it will be rounded up to
+! exactly 1. Since ocean dominates, the lake
+! will be removed.
+
+ print*,'Test point 6.'
+
+ binary_lake = 0           ! Keep fractional lake.
+ lat2d(1,1) = 30.0         ! Point at 30N
+ ocn_frac(1,1) = 0.99999   ! Very high % of ocean.
+ lake_frac(1,1) = 0.24     ! Some lake.
+ lake_depth(1,1) = 50.0    ! Lake depth.
+ land_frac(1,1) = -99.     ! Should be 0.
+ slmsk(1,1) = -99.         ! Should be 0.
+
+ call merge(lon, lat, binary_lake, lat2d, ocn_frac, &
+            lake_frac, lake_depth, land_frac, slmsk)
+
+ if ( abs(lake_frac(1,1) - 0.) > epsilon) stop 12
+ if ( abs(lake_depth(1,1) - 0.) > epsilon) stop 12
+ if ( abs(land_frac(1,1) - 0.) > epsilon) stop 12
+ if ( abs(slmsk(1,1) - 0.) > epsilon) stop 12
+
+! Test point 7
+! No ocean. Some lake, but it is near Antarctica.
+! Lakes near Antarctica are to be removed. Point
+! becomes all land.
+
+ print*,'Test point 7.'
+
+ binary_lake = 0           ! Keep fractional lake.
+ lat2d(1,1) = -70.0        ! Point at 70S.
+ ocn_frac(1,1) = 0.0       ! No ocean.
+ lake_frac(1,1) = .75      ! Some lake.
+ lake_depth(1,1) = 100.0   ! Lake depth.
+ land_frac(1,1) = -99.     ! Should be 1.
+ slmsk(1,1) = -99.         ! Should be 1.
+
+ call merge(lon, lat, binary_lake, lat2d, ocn_frac, &
+            lake_frac, lake_depth, land_frac, slmsk)
+
+ if ( abs(lake_frac(1,1) - 0.) > epsilon) stop 14
+ if ( abs(lake_depth(1,1) - 0.) > epsilon) stop 14
+ if ( abs(land_frac(1,1) - 1.) > epsilon) stop 14
+ if ( abs(slmsk(1,1) - 1.) > epsilon) stop 14
+
+! Test point 8
+! No ocean. Some lake. Assume binary yes/no lake.
+! In this case, the lake fraction will become zero and
+! the land fraction will become 1.
+
+ print*,'Test point 8.'
+
+ binary_lake = 1           ! Keep fractional lake.
+ lat2d(1,1) = 30.0         ! Point at 30N.
+ ocn_frac(1,1) = 0.0       ! No ocean.
+ lake_frac(1,1) = .15      ! Less than 50% lake.
+ lake_depth(1,1) = 100.0   ! Lake depth.
+ land_frac(1,1) = -99.     ! Should be 1.
+ slmsk(1,1) = -99.         ! Should be 1.
+
+ call merge(lon, lat, binary_lake, lat2d, ocn_frac, &
+            lake_frac, lake_depth, land_frac, slmsk)
+
+ if ( abs(lake_frac(1,1) - 0.) > epsilon) stop 16
+ if ( abs(lake_depth(1,1) - 0.) > epsilon) stop 16
+ if ( abs(land_frac(1,1) - 1.) > epsilon) stop 16
+ if ( abs(slmsk(1,1) - 1.) > epsilon) stop 16
+
+ print*, "OK"
+
+ print*, "SUCCESS!"
+
+ end program ftst_merge
diff --git a/tests/ocean_merge/ftst_read_lake_mask.F90 b/tests/ocean_merge/ftst_read_lake_mask.F90
new file mode 100644
index 000000000..3754455cc
--- /dev/null
+++ b/tests/ocean_merge/ftst_read_lake_mask.F90
@@ -0,0 +1,42 @@
+! Unit test for the read_lake_mask routine.
+!
+! Reads a 6x4 version of the lake mask NetCDF file and
+! checks values from the lake fraction, lake depth
+! and latitude records. If differences exceed a
+! threshold, then the test fails.
+!
+ program read_lake_info
+
+ implicit none
+
+ integer, parameter     :: lon = 6
+ integer, parameter     :: lat = 4
+ integer, parameter     :: tile = 1
+
+ character(len=3)       :: pth2="./"
+ character(len=3)       :: atmres="C48"
+
+ real, parameter        :: thresh = 0.001
+ real                   :: lake_frac(lon,lat)
+ real                   :: lake_depth(lon,lat)
+ real                   :: lat2d(lon,lat)
+
+ print*,"Call routine read_lake_mask."
+
+ call read_lake_mask(pth2,atmres,tile,lon,lat,lake_frac, &
+                           lake_depth,lat2d)
+
+ print*,"Check records."
+
+ if (abs(lake_depth(1,1)-0.0) > thresh) stop 2
+ if (abs(lake_depth(5,2)-68.3817) > thresh) stop 4
+ if (abs(lake_frac(4,3)-1.0) > thresh) stop 6
+ if (abs(lake_frac(4,4)-0.0) > thresh) stop 8
+ if (abs(lat2d(1,4)-0.399218) > thresh) stop 10
+ if (abs(lat2d(6,1)-(-1.87402)) > thresh) stop 12
+
+ print*,"OK"
+
+ print*,"SUCCESS"
+
+ end program read_lake_info
diff --git a/tests/ocean_merge/ftst_read_nml.F90 b/tests/ocean_merge/ftst_read_nml.F90
new file mode 100644
index 000000000..9ccb3156d
--- /dev/null
+++ b/tests/ocean_merge/ftst_read_nml.F90
@@ -0,0 +1,32 @@
+! Unit test for the read_nml routine.
+!
+! Read a sample namelist and check the data in each
+! variable against expected values.
+!
+ program read_namelist
+
+ implicit none
+
+ character(len=120) :: ocean_mask_dir
+ character(len=120) :: lake_mask_dir
+ character(len=120) :: out_dir
+ character(len=10)  :: atmres,ocnres
+ 
+ integer            :: binary_lake
+
+ print*,"Call routine read_nml."
+
+ call read_nml(ocean_mask_dir, lake_mask_dir, atmres,ocnres,out_dir,binary_lake)
+
+ if (trim(ocean_mask_dir) /= "/ocean/mask/dir") stop 2
+ if (trim(lake_mask_dir) /= "/lake/mask/dir") stop 4
+ if (trim(atmres) /= "C96") stop 6
+ if (trim(ocnres) /= "mx500") stop 8
+ if (trim(out_dir) /= "/out/dir") stop 10
+ if (binary_lake /= 1) stop 12
+
+ print*, "OK"
+
+ print*, "SUCCESS!"
+
+ end program read_namelist
diff --git a/tests/ocean_merge/ftst_read_ocean_frac.F90 b/tests/ocean_merge/ftst_read_ocean_frac.F90
new file mode 100644
index 000000000..3e83377ed
--- /dev/null
+++ b/tests/ocean_merge/ftst_read_ocean_frac.F90
@@ -0,0 +1,37 @@
+! Unit test for the read_ocean_frac routine.
+!
+! Reads a 6x5 version of the MOM6 ocean mask file and
+! checks values from the ocean fraction record.
+! If differences exceed a threshold, then the test fails.
+!
+ program read_ocean_info
+
+ implicit none
+
+ integer, parameter     :: lon = 6
+ integer, parameter     :: lat = 5
+ integer, parameter     :: tile = 1
+
+ character(len=3)       :: pth1="./"
+ character(len=3)       :: atmres="C48"
+ character(len=5)       :: ocnres="mx500"
+
+ real, parameter        :: thresh = 0.001
+ real                   :: ocn_frac(lon,lat)
+
+ print*,"Call routine read_ocean_frac."
+
+ call read_ocean_frac(pth1,atmres,ocnres,tile,lon,lat,ocn_frac)
+
+ print*,"Check records."
+
+ if (abs(ocn_frac(1,1)-1.0) > thresh) stop 2
+ if (abs(ocn_frac(6,5)-0.0) > thresh) stop 4
+ if (abs(ocn_frac(4,5)-0.0917) > thresh) stop 6
+ if (abs(ocn_frac(3,1)-0.162347) > thresh) stop 8
+
+ print*,"OK"
+
+ print*,"SUCCESS"
+
+ end program read_ocean_info