Add unit testing to ocean_merge program (ufs-community#997)

Refactor the program into subroutines.

Add unit testing.

Update the c96.uniform grid_gen regression test to exercise the ocean_merge
program.

Fixes ufs-community#944.

cmake/mpiexec.hercules

@@ -6,7 +6,7 @@
 #   $2+ - Executable and its arguments
 #
 
-ACCOUNT=
+ACCOUNT=fv3-cpu
 QOS=debug
 
 NP=$1

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

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.

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

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