From 402849a183c806d7548cf034e3ccf7815a222418 Mon Sep 17 00:00:00 2001
From: "Shan.Sun" <shan.sun@noaa.gov>
Date: Fri, 29 Sep 2023 04:20:56 +0000
Subject: [PATCH 01/15] -- Fix coastal lakes with revised code -- Link orog to
 /scratch1/NCEPDEV/global/glopara/fix/raw/orog/ to read MODISP lake fraction
 and GLDBV3 lake depth for now -- Add options to choose different lake
 datasets

Fixes #854
---
 driver_scripts/driver_grid.hera.sh           |   7 +-
 sorc/orog_mask_tools.fd/inland.fd/inland.F90 |  38 +-
 sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 | 501 ++++++++++++-------
 ush/fv3gfs_make_lake.sh                      |  36 +-
 4 files changed, 388 insertions(+), 194 deletions(-)

diff --git a/driver_scripts/driver_grid.hera.sh b/driver_scripts/driver_grid.hera.sh
index 63d1b872e..a054a3e24 100755
--- a/driver_scripts/driver_grid.hera.sh
+++ b/driver_scripts/driver_grid.hera.sh
@@ -106,10 +106,13 @@ export soil_type_src="bnu.v3.30s" #  Soil type data.
                                 # For Beijing Norm. Univ. data
                                 # 1) "bnu.v3.30s" for global 30s data.
 
+export reg_domain=na3km   # '3km', 'na3km', '13km', and 'na13km'
+export lake_data_srce=MODISP_GLDBV3 # 'GLDBV3', 'MODIS_GLOBATHY', 'MODISP_GLDBV3', and 'VIIRS_GLDBV3'
+
 if [ $gtype = uniform ]; then
   export res=96
-  export add_lake=false        # Add lake frac and depth to orography data.
-  export lake_cutoff=0.20      # lake frac < lake_cutoff ignored when add_lake=T
+  export add_lake=true         # Add lake frac and depth to orography data.
+  export lake_cutoff=0.50      # lake frac < lake_cutoff ignored when add_lake=T
 elif [ $gtype = stretch ]; then
   export res=96
   export stretch_fac=1.5       # Stretching factor for the grid
diff --git a/sorc/orog_mask_tools.fd/inland.fd/inland.F90 b/sorc/orog_mask_tools.fd/inland.fd/inland.F90
index 057c648d3..66398e2b2 100644
--- a/sorc/orog_mask_tools.fd/inland.fd/inland.F90
+++ b/sorc/orog_mask_tools.fd/inland.fd/inland.F90
@@ -15,7 +15,7 @@ PROGRAM inland_mask
 
   REAL, ALLOCATABLE :: inland(:,:,:)
   REAL, ALLOCATABLE :: land_frac(:,:,:)
-  INTEGER :: i_ctr, j_ctr, tile_beg, tile_end
+  INTEGER :: tile_beg, tile_end
   INTEGER :: cs_res, x_res, y_res
   CHARACTER(len=32) :: arg
   INTEGER :: stat
@@ -23,7 +23,7 @@ PROGRAM inland_mask
   REAL :: cutoff
   CHARACTER(len=1) :: reg
 
-  LOGICAL, ALLOCATABLE :: done(:,:,:)
+  LOGICAL, ALLOCATABLE :: done(:,:,:) 
 
   CALL getarg(0, arg) ! get the program name
   IF (iargc() /= 3 .AND. iargc() /= 4) THEN
@@ -78,14 +78,19 @@ PROGRAM inland_mask
 !! @author Ning Wang  
 SUBROUTINE mark_global_inland(cs_res)
   INTEGER, INTENT(IN) :: cs_res
+  INTEGER :: i_seed, j_seed
  
   ALLOCATE(done(cs_res,cs_res,6))
   ALLOCATE(inland(cs_res,cs_res,6))
   done = .false.
   inland = 1.0
-  i_ctr = cs_res/2; j_ctr = cs_res/2
 
-  CALL mark_global_inland_rec_d(i_ctr, j_ctr, 2, 0)
+  i_seed = cs_res/2; j_seed = cs_res/2
+  CALL mark_global_inland_rec_d(i_seed, j_seed, 2, 0)
+
+! to make sure black sea is excluded
+  i_seed = REAL(cs_res)/32.0*3; j_seed = i_seed 
+  CALL mark_global_inland_rec_d(i_seed, j_seed, 3, 0)
   
   DEALLOCATE(done)
 
@@ -99,6 +104,7 @@ END SUBROUTINE mark_global_inland
 SUBROUTINE mark_inland_reg(cs_res)
   INTEGER, INTENT(IN) :: cs_res
   INTEGER :: i_seed, j_seed
+  INTEGER :: i
 
   ALLOCATE(done(x_res,y_res,1))
   ALLOCATE(inland(x_res,y_res,1))
@@ -116,6 +122,30 @@ SUBROUTINE mark_inland_reg(cs_res)
   i_seed = x_res/3; j_seed = 1 
   CALL mark_regional_inland_rec_d(i_seed, j_seed, 1, 0)
 
+  j_seed = 1
+  DO i = 1, x_res
+    CALL mark_regional_inland_rec_d(i, j_seed, 1, 0)
+  ENDDO
+
+  j_seed = y_res
+  DO i = x_res/2, x_res
+    CALL mark_regional_inland_rec_d(i, j_seed, 1, 0)
+  ENDDO
+
+! set up additional 3 seeds for ESG CONUS grid 
+!  i_seed = 1600; j_seed = 1040 
+!  i_seed = x_res - 10; j_seed = y_res
+!  CALL mark_regional_inland_rec_d(i_seed, j_seed, 1, 0)
+
+!  i_seed = x_res - 60; j_seed = y_res
+!  CALL mark_regional_inland_rec_d(i_seed, j_seed, 1, 0)
+
+!  i_seed = x_res - 275; j_seed = y_res
+!  CALL mark_regional_inland_rec_d(i_seed, j_seed, 1, 0)
+
+!  i_seed = 500; j_seed = 1 
+!  CALL mark_regional_inland_rec_d(i_seed, j_seed, 1, 0)
+  
   DEALLOCATE(done)
 
 END SUBROUTINE mark_inland_reg
diff --git a/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 b/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
index 47f2fff60..f5587ea3d 100644
--- a/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
+++ b/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
@@ -39,19 +39,21 @@ PROGRAM lake_frac
     REAL, PARAMETER :: d2r = acos(-1.0) / 180.0
     REAL, PARAMETER :: r2d = 180.0 /acos(-1.0) 
     REAL, PARAMETER :: pi = acos(-1.0) 
-    REAL*8, PARAMETER :: gppdeg = 119.99444445
-    REAL*8, PARAMETER :: delta = 1.0 / 119.99444445
+!   REAL*8, PARAMETER :: gppdeg = 119.99444445
+!   REAL*8, PARAMETER :: delta = 1.0 / 119.99444445
+    REAL*8, PARAMETER :: gppdeg = 120.0
+    REAL*8, PARAMETER :: delta = 1.0 / 120.0
 
-    CHARACTER(len=32) :: arg
+    CHARACTER(len=32) :: arg, lakestatus_srce, lakedepth_srce
     CHARACTER(len=256) :: lakedata_path
     INTEGER :: stat
     
     CALL getarg(0, arg) ! get the program name
-    IF (iargc() /= 3 .AND. iargc() /= 4) THEN
+    IF (iargc() /= 5 .AND. iargc() /= 6) THEN
       PRINT*, 'Usage: ', trim(arg), & 
-       ' [tile_num (0:all tiles, 7:regional)] [resolution (48,96, ...)] [path to lake data file]'  
+       ' [tile_num (0:all tiles, 7:regional)] [resolution (48,96, ...)] [lake data path] [lake status source] [lake depth source]'  
       PRINT*, 'Or: ', trim(arg), & 
-       ' [tile_num (0:all tiles, 7:regional)] [resolution (48,96, ...)] [path to lake data file] [lake_cutoff]'  
+       ' [tile_num (0:all tiles, 7:regional)] [resolution (48,96, ...)] [lake data path] [lake status source] [lake depth source] [lake_cutoff]'
       STOP
     ENDIF
     CALL getarg(1, arg)
@@ -59,14 +61,20 @@ PROGRAM lake_frac
     CALL getarg(2, arg)
     READ(arg,*,iostat=stat) cs_res
     CALL getarg(3, lakedata_path)
+    CALL getarg(4, lakestatus_srce)
+    CALL getarg(5, lakedepth_srce)
 
-    IF (iargc() == 3) THEN
+    IF (iargc() == 5) THEN
       lake_cutoff = 0.20
     ELSE
-      CALL getarg(4, arg)
+      CALL getarg(6, arg)
       READ(arg,*,iostat=stat) lake_cutoff
     ENDIF
 
+    PRINT*, 'lake status source:', trim(lakestatus_srce) 
+    PRINT*, 'lake depth source:', trim(lakedepth_srce) 
+    PRINT*, 'lake cutoff:', lake_cutoff 
+
     IF (tile_req == 0) THEN
       tile_beg = 1; tile_end = 6
       PRINT*, 'Process tile 1 - 6 at resolution C',cs_res
@@ -91,17 +99,31 @@ PROGRAM lake_frac
     ELSE
       CALL  read_cubed_sphere_reg_grid(cs_res, cs_grid, 3, res_x, res_y) 
     ENDIF
-    ! compute source grid 
+
+    ! allocate and compute source grid 
     ALLOCATE(src_grid_lon(nlon), src_grid_lat(nlat))
-    DO i = 1, nlon
-      src_grid_lon(i) = -180.0 + delta * (i-1)
-    ENDDO
-    DO i = 1, nlat
-      src_grid_lat(i) = 90.0 - delta * (i-1) 
-    ENDDO
+
+    IF (lakestatus_srce == "GLDBV2" .OR. lakestatus_srce == "GLDBV3") THEN 
+      !  GLDB data points are at the lower right corners of the grid cells
+      DO i = 1, nlon
+        src_grid_lon(i) = -180.0 + delta * i
+      ENDDO
+      DO i = 1, nlat
+        src_grid_lat(i) = 90.0 - delta * i 
+      ENDDO
+    ENDIF 
+
+    IF (lakestatus_srce == "MODISP" .OR. lakestatus_srce == "VIIRS") THEN 
+      !  GLDB data points are at the upprt left corners of the grid cells
+      DO i = 1, nlon
+        src_grid_lon(i) = -180.0 + delta * (i-1)
+      ENDDO
+      DO i = 1, nlat
+        src_grid_lat(i) = 90.0 - delta * (i-1) 
+      ENDDO
+    ENDIF
 
     ! read in lake data file
-!    sfcdata_path = '/scratch1/NCEPDEV/global/glopara/fix/orog/'
     lakedata_path = trim(lakedata_path) // "/"
     ALLOCATE(lakestatus(nlon*nlat),lakedepth(nlon*nlat)) 
     PRINT*, 'Read in lake data file ...'
@@ -115,7 +137,7 @@ PROGRAM lake_frac
     PRINT*, 'Calculate lake fraction and average depth for cubed-sphere cells ...'
     CALL cal_lake_frac_depth(lakestatus,cs_lakestatus,lakedepth,cs_lakedepth)
 
-    ! write lake status (in terms of fraction) and lake depth(average)
+    ! write lake status (in terms of fraction) and lake depth(average, in meters)
     ! to a netcdf file
     IF (tile_req /= 7) THEN
       PRINT*, 'Write lake fraction/depth on cubed sphere grid cells to netCDF files ...'
@@ -157,7 +179,9 @@ SUBROUTINE cal_lake_frac_depth(lakestat,cs_lakestat,lakedpth,cs_lakedpth)
     INTEGER :: src_grid_lon_beg1,src_grid_lon_end1,src_grid_lon_beg2,src_grid_lon_end2
     INTEGER :: grid_ct, lake_ct, co_gc, tmp
 
-    REAL*8 :: lake_dpth_sum
+    INTEGER*1 :: lkst
+    INTEGER*2 :: lkdp
+    REAL*8 :: lake_dpth_sum, lake_avg_frac
     LOGICAL :: two_section, enclosure_cnvx 
 
     IF (tile_req /= 7) THEN
@@ -234,13 +258,11 @@ SUBROUTINE cal_lake_frac_depth(lakestat,cs_lakestat,lakedpth,cs_lakedpth)
         src_grid_lon_end = nint((180.0+lonmax)*gppdeg+0.5) 
        
         IF (src_grid_lat_beg > src_grid_lat_end) THEN
-          print*,'switch!!!'
           tmp = src_grid_lat_beg
           src_grid_lat_beg = src_grid_lat_end
           src_grid_lat_end = tmp
         ENDIF
         IF (src_grid_lon_beg > src_grid_lon_end) THEN
-          print*,'switch!!!'
           tmp = src_grid_lon_beg
           src_grid_lon_beg = src_grid_lon_end
           src_grid_lon_end = tmp
@@ -263,88 +285,65 @@ SUBROUTINE cal_lake_frac_depth(lakestat,cs_lakestat,lakedpth,cs_lakedpth)
 #endif
         lake_ct = 0; grid_ct = 0
         lake_dpth_sum = 0.0
+        lake_avg_frac = 0.0
         DO j = src_grid_lat_beg, src_grid_lat_end, stride_lat
           stride_lon = int(1.0/cos(src_grid_lat(j)*d2r)*REAL(stride_lat))
 #ifdef DIAG_N_VERBOSE
-        IF (j == src_grid_lat_beg) THEN
-          PRINT*, 'lon index range and stride', &
-            src_grid_lon_beg,src_grid_lon_end,stride_lon
-          PRINT*, 'lon range ', &
-            src_grid_lon(src_grid_lon_beg),src_grid_lon(src_grid_lon_end)
-          IF (two_section == .true.) THEN
-          PRINT*, 'section1 index lon range and stride', &
-            src_grid_lon_beg1,src_grid_lon_end1,stride_lon
-          PRINT*, 'section1 lon range ', &
-            src_grid_lon(src_grid_lon_beg1),src_grid_lon(src_grid_lon_end1)
-          PRINT*, 'section2 index lon range and stride', &
-            src_grid_lon_beg2,src_grid_lon_end2,stride_lon
-          PRINT*, 'section2 lon range ', &
-            src_grid_lon(src_grid_lon_beg2),src_grid_lon(src_grid_lon_end2)
+          IF (j == src_grid_lat_beg) THEN
+            PRINT*, 'lon index range and stride', &
+              src_grid_lon_beg,src_grid_lon_end,stride_lon
+            PRINT*, 'lon range ', &
+              src_grid_lon(src_grid_lon_beg),src_grid_lon(src_grid_lon_end)
+            IF (two_section == .true.) THEN
+              PRINT*, 'section1 index lon range and stride', &
+                src_grid_lon_beg1,src_grid_lon_end1,stride_lon
+              PRINT*, 'section1 lon range ', &
+                src_grid_lon(src_grid_lon_beg1),src_grid_lon(src_grid_lon_end1)
+              PRINT*, 'section2 index lon range and stride', &
+                src_grid_lon_beg2,src_grid_lon_end2,stride_lon
+              PRINT*, 'section2 lon range ', &
+                src_grid_lon(src_grid_lon_beg2),src_grid_lon(src_grid_lon_end2)
+            ENDIF
           ENDIF
-        ENDIF
 #endif
-          IF (two_section .EQV. .false.) THEN
+          IF (two_section == .false.) THEN
             DO i =  src_grid_lon_beg, src_grid_lon_end, stride_lon
               p(1) = src_grid_lat(j); p(2) = src_grid_lon(i) 
               p(:) = p(:)*d2r
-              IF(enclosure_cnvx(v, 4, p, co_gc) .EQV. .true.) THEN
+              IF(enclosure_cnvx(v, 4, p, co_gc) == .true.) THEN
                 grid_ct = grid_ct+1
-                IF (lakestat((j-1)*nlon+i) /= 0) THEN
-                  lake_ct = lake_ct+1
-                  IF (lakedpth((j-1)*nlon+i) < 0) THEN
-                    IF (lakestat((j-1)*nlon+i) == 4) THEN
-                      lake_dpth_sum = lake_dpth_sum+30.0
-                    ELSE
-                      lake_dpth_sum = lake_dpth_sum+100.0
-                    ENDIF
-                  ELSE
-                    lake_dpth_sum = lake_dpth_sum+REAL(lakedpth((j-1)*nlon+i))
-                  ENDIF
-                ENDIF
+                lkst = lakestat((j-1)*nlon+i); lkdp = lakedpth((j-1)*nlon+i)
+                CALL lake_cell_comp(lkst, lkdp, lake_ct, lake_avg_frac, lake_dpth_sum)
               ENDIF   
             ENDDO
           ELSE 
             DO i =  src_grid_lon_beg1, src_grid_lon_end1, stride_lon
               p(1) = src_grid_lat(j); p(2) = src_grid_lon(i) 
               p(:) = p(:)*d2r
-              IF(enclosure_cnvx(v, 4, p, co_gc) .EQV. .true.) THEN
+              IF(enclosure_cnvx(v, 4, p, co_gc) == .true.) THEN
                 grid_ct = grid_ct+1
-                IF (lakestat((j-1)*nlon+i) /= 0) THEN
-                  lake_ct = lake_ct+1
-                  IF (lakedpth((j-1)*nlon+i) < 0) THEN
-                    IF (lakestat((j-1)*nlon+i) == 4) THEN
-                      lake_dpth_sum = lake_dpth_sum+30.0
-                    ELSE
-                      lake_dpth_sum = lake_dpth_sum+100.0
-                    ENDIF
-                  ELSE
-                    lake_dpth_sum = lake_dpth_sum+REAL(lakedpth((j-1)*nlon+i))
-                  ENDIF
-                ENDIF
+                lkst = lakestat((j-1)*nlon+i); lkdp = lakedpth((j-1)*nlon+i)
+                CALL lake_cell_comp(lkst, lkdp, lake_ct, lake_avg_frac, lake_dpth_sum)
               ENDIF   
             ENDDO
             DO i =  src_grid_lon_beg2, src_grid_lon_end2, stride_lon
               p(1) = src_grid_lat(j); p(2) = src_grid_lon(i) 
               p(:) = p(:)*d2r
-              IF(enclosure_cnvx(v, 4, p, co_gc) .EQV. .true.) THEN
+              IF(enclosure_cnvx(v, 4, p, co_gc) == .true.) THEN
                 grid_ct = grid_ct+1
-                IF (lakestat((j-1)*nlon+i) /= 0) THEN
-                  lake_ct = lake_ct+1
-                  IF (lakedpth((j-1)*nlon+i) < 0) THEN
-                    IF (lakestat((j-1)*nlon+i) == 4) THEN
-                      lake_dpth_sum = lake_dpth_sum+30.0
-                    ELSE
-                      lake_dpth_sum = lake_dpth_sum+100.0
-                    ENDIF
-                  ELSE
-                    lake_dpth_sum = lake_dpth_sum+REAL(lakedpth((j-1)*nlon+i))
-                  ENDIF
-                ENDIF
+                lkst = lakestat((j-1)*nlon+i); lkdp = lakedpth((j-1)*nlon+i)
+                CALL lake_cell_comp(lkst, lkdp, lake_ct, lake_avg_frac, lake_dpth_sum)
               ENDIF   
             ENDDO
           ENDIF
         ENDDO
-        cs_lakestat(cs_data_idx)=REAL(lake_ct)/REAL(grid_ct)
+        IF (lakestatus_srce == "GLDBV3" .OR. lakestatus_srce == "GLDBV2" .OR. &
+            lakestatus_srce == "VIIRS" ) THEN 
+          cs_lakestat(cs_data_idx)=REAL(lake_ct)/REAL(grid_ct)
+        ENDIF
+        IF (lakestatus_srce == "MODISP") THEN 
+          cs_lakestat(cs_data_idx)=lake_avg_frac/REAL(grid_ct)
+        ENDIF
         IF (lake_ct /= 0) THEN
           cs_lakedpth(cs_data_idx)=lake_dpth_sum/REAL(lake_ct)/10.0 !convert to meter
         ELSE
@@ -370,6 +369,36 @@ SUBROUTINE cal_lake_frac_depth(lakestat,cs_lakestat,lakedpth,cs_lakedpth)
 
 END SUBROUTINE cal_lake_frac_depth 
 
+SUBROUTINE lake_cell_comp(lkst, lkdp, lake_ct, lake_avg_frac, lake_dpth_sum)
+    INTEGER*1, INTENT(IN) :: lkst
+    INTEGER*2, INTENT(IN) :: lkdp
+    INTEGER, INTENT(OUT) :: lake_ct
+    REAL*8, INTENT(OUT) :: lake_avg_frac, lake_dpth_sum
+    
+    IF (lkst /= 0) THEN ! lake point 
+      lake_ct = lake_ct+1
+      IF (lakestatus_srce == "GLDBV3" .OR. lakestatus_srce == "GLDBV2") THEN 
+        IF (lkdp <= 0) THEN
+          IF (lkst == 4) THEN
+            lake_dpth_sum = lake_dpth_sum+30.0
+          ELSE
+            lake_dpth_sum = lake_dpth_sum+100.0
+          ENDIF
+        ELSE
+          lake_dpth_sum = lake_dpth_sum+REAL(lkdp)
+        ENDIF
+      ENDIF
+      IF (lakestatus_srce == "MODISP" .OR. lakestatus_srce == "VIIRS") THEN 
+        lake_avg_frac = lake_avg_frac + REAL(lkst) / 100.0 
+        IF (lkdp <= 0) THEN
+          lake_dpth_sum = lake_dpth_sum+100.0
+        ELSE
+          lake_dpth_sum = lake_dpth_sum+REAL(lkdp)
+        ENDIF
+      ENDIF 
+    ENDIF
+END SUBROUTINE lake_cell_comp
+
 !> Read the latitude and longitude for a cubed-sphere
 !! grid from the 'grid' files. For global grids, all
 !! six sides are returned.
@@ -400,7 +429,6 @@ SUBROUTINE read_cubed_sphere_grid(res, grid)
       tile_beg = tile_req; tile_end = tile_req
     ENDIF
     WRITE(res_ch,'(I4)') res
-!    gridfile_path = "/scratch1/NCEPDEV/global/glopara/fix/fix_fv3/C"//trim(adjustl(res_ch))//"/"
     gridfile_path = "./"
     DO i = tile_beg, tile_end
       WRITE(ich, '(I1)') i
@@ -521,14 +549,38 @@ SUBROUTINE read_lakedata(lakedata_path,lake_stat,lake_dpth,nlat,nlon)
 
     data_sz = nlon*nlat
 
-   ! read in 30 arc seconds lake data 
-    lakefile = trim(lakedata_path) // "GlobalLakeStatus.dat"
-    OPEN(10, file=lakefile, form='unformatted', access='direct', recl=data_sz*1)
+! read in 30 arc seconds lake data 
+   ! lake fraction data
+    lakefile = trim(lakedata_path) // "GlobalLakeStatus_GLDBv3release.dat"
+    IF (lakestatus_srce == "GLDBV2") THEN 
+      lakefile = trim(lakedata_path) // "GlobalLakeStatus.dat"
+    ENDIF
+    IF (lakestatus_srce == "GLDBV3") THEN 
+      lakefile = trim(lakedata_path) // "GlobalLakeStatus_GLDBv3release.dat"
+    ENDIF
+    IF (lakestatus_srce == "MODISP") THEN 
+!     lakefile = trim(lakedata_path) // "GlobalLakeStatus_MODIS15s.dat"
+      lakefile = trim(lakedata_path) // "GlobalLakeStatus_MODISp.dat"
+    ENDIF
+    IF (lakestatus_srce == "VIIRS") THEN 
+      lakefile = trim(lakedata_path) // "GlobalLakeStatus_VIIRS.dat"
+    ENDIF
+    OPEN(10, file=lakefile, form='unformatted', access='direct', status='old', recl=data_sz*1)
     READ(10,rec=1) lake_stat
     CLOSE(10)
 
-    lakefile = trim(lakedata_path) // "GlobalLakeDepth.dat"
-    OPEN(10, file=lakefile, form='unformatted', access='direct', recl=data_sz*2)
+   ! lake depth data
+    lakefile = trim(lakedata_path) // "GlobalLakeDepth_GLDBv3release.dat"
+    IF (lakedepth_srce == "GLDBV2") THEN 
+      lakefile = trim(lakedata_path) // "GlobalLakeDepth.dat"
+    ENDIF
+    IF (lakedepth_srce == "GLDBV3") THEN 
+      lakefile = trim(lakedata_path) // "GlobalLakeDepth_GLDBv3release.dat"
+    ENDIF
+    IF (lakedepth_srce == "GLOBATHY") THEN 
+      lakefile = trim(lakedata_path) // "GlobalLakeDepth_GLOBathy.dat"
+    ENDIF
+    OPEN(10, file=lakefile, form='unformatted', access='direct', status='old', recl=data_sz*2)
     READ(10,rec=1) lake_dpth
     CLOSE(10)
 
@@ -563,7 +615,6 @@ SUBROUTINE write_lakedata_to_orodata(cs_res, cs_lakestat, cs_lakedpth)
 
     INTEGER :: i, j
 
-!    include "netcdf.inc"
     tile_sz = cs_res*cs_res
 
     WRITE(res_ch,'(I4)') cs_res
@@ -598,28 +649,51 @@ SUBROUTINE write_lakedata_to_orodata(cs_res, cs_lakestat, cs_lakedpth)
       CALL nc_opchk(stat, "nf90_put_att: lake_frac:long_name") 
       stat = nf90_put_att(ncid, lake_frac_id,'unit','fraction')
       CALL nc_opchk(stat, "nf90_put_att: lake_frac:unit") 
-      write(string,'(a,f5.2)') 'based on GLDBv2 (Choulga et al. 2014); missing lakes &
-        added based on land_frac in this dataset; lake_frac cutoff is',lake_cutoff
+      IF (lakestatus_srce == "GLDBV3") THEN 
+        write(string,'(a,es8.1)') 'based on GLDBv3 (Choulga et al. 2019); missing lakes &
+        added based on land_frac in this dataset; lake_frac cutoff:',lake_cutoff
+      ELSE IF (lakestatus_srce == "GLDBV2") THEN 
+        write(string,'(a,es8.1)') 'based on GLDBv2 (Choulga et al. 2014); missing lakes &
+        added based on land_frac in this dataset; lake_frac cutoff:',lake_cutoff
+      ELSE IF (lakestatus_srce == "MODISP") THEN
+        write(string,'(a,es8.1)') 'based on MODIS (2011-2015) product updated with two &
+        Landsat products: the JRC water product (2016-2020) and the GLC-FCS30 (2020); &
+        the source data set was created by Chengquan Huang of UMD; &  
+        lake_frac cutoff:',lake_cutoff
+      ELSE IF (lakestatus_srce == "VIIRS") THEN
+        write(string,'(a,es8.1)') 'based on multi-year VIIRS global surface type & 
+        classification map (2012-2019); the source data set was created by &
+        Chengquan Huang of UMD and Michael Barlage of NOAA; &
+        lake_frac cutoff:',lake_cutoff
+      ENDIF
       stat = nf90_put_att(ncid, lake_frac_id,'description',trim(string))
       CALL nc_opchk(stat, "nf90_put_att: lake_frac:description") 
 #endif
       stat = nf90_def_var(ncid,"lake_depth",NF90_FLOAT,dimids,lake_depth_id)
       CALL nc_opchk(stat, "nf90_def_var: lake_depth")
 #ifdef ADD_ATT_FOR_NEW_VAR
-    stat = nf90_put_att(ncid, lake_depth_id,'coordinates','geolon geolat')
-    CALL nc_opchk(stat, "nf90_put_att: lake_depth:coordinates") 
+      stat = nf90_put_att(ncid, lake_depth_id,'coordinates','geolon geolat')
+      CALL nc_opchk(stat, "nf90_put_att: lake_depth:coordinates") 
       stat = nf90_put_att(ncid, lake_depth_id,'long_name','lake depth')
       CALL nc_opchk(stat, "nf90_put_att: lake_depth:long_name") 
       stat = nf90_put_att(ncid, lake_depth_id,'unit','meter')
       CALL nc_opchk(stat, "nf90_put_att: lake_depth:long_name") 
-      stat = nf90_put_att(ncid, lake_depth_id,'description', &
+      IF (lakedepth_srce == "GLDBV3") THEN 
+        stat = nf90_put_att(ncid, lake_depth_id,'description', &
+        'based on GLDBv3 (Choulga et al. 2019); missing depth set to 10m &
+        (except to 211m in Caspian Sea); spurious large pos. depths are left unchanged.')
+      ELSE IF (lakedepth_srce == "GLDBV2") THEN 
+        stat = nf90_put_att(ncid, lake_depth_id,'description', &
         'based on GLDBv2 (Choulga et al. 2014); missing depth set to 10m &
         (except to 211m in Caspian Sea); spurious large pos. depths are left unchanged.')
+      ELSE IF (lakedepth_srce == "GLOBATHY") THEN
+        stat = nf90_put_att(ncid, lake_depth_id,'description', &
+        'based on GLOBathy data resampled and projected to the MODIS domain.')
+      ENDIF
       CALL nc_opchk(stat, "nf90_put_att: lake_depth:description") 
 #endif
 
-      write(string,'(a,es8.1)') 'land_frac and lake_frac are adjusted such that their sum '// &
-                                'is 1 at points where inland=1; land_frac cutoff is',land_cutoff
+      write(string,'(a,es8.1)') 'land_frac and lake_frac are adjusted such that their sum is 1 at points where inland=1; land_frac cutoff is',land_cutoff
       stat = nf90_put_att(ncid, land_frac_id,'description',trim(string))
       CALL nc_opchk(stat, "nf90_put_att: land_frac:description") 
 
@@ -660,46 +734,37 @@ SUBROUTINE write_lakedata_to_orodata(cs_res, cs_lakestat, cs_lakedpth)
       lake_frac (:) = cs_lakestat ((tile_num-1)*tile_sz+1:tile_num*tile_sz)
       lake_depth(:) = cs_lakedepth((tile_num-1)*tile_sz+1:tile_num*tile_sz)
 
-! add Caspian Sea and Aral Sea to lake_frac and lake_depth
-      IF (tile_num == 2 .or. tile_num == 3) THEN
-        DO i = 1, tile_sz
-          IF (land_frac(i) < 0.9 .AND. lake_frac(i) < 0.1) THEN
-            IF (geolat(i) > 35.0 .AND. geolat(i) <= 50.0 .AND. &
-                geolon(i) > 45.0 .AND. geolon(i) <= 55.0) THEN
-              lake_frac(i) = 1.-land_frac(i)
-              lake_depth(i) = 211.0
-            ENDIF
-            IF (geolat(i) > 35.0 .AND. geolat(i) <= 50.0 .AND. &
-                geolon(i) > 57.0 .AND. geolon(i) <= 63.0) THEN
-              lake_frac(i) = 1.-land_frac(i)
-              lake_depth(i) = 10.0
-            ENDIF
-          ENDIF
-        ENDDO
-      ENDIF
-
-! adjust land_frac and lake_frac, and make sure land_frac+lake_frac=1 at inland points
-      DO i = 1, tile_sz
-        if (lake_frac(i) >= lake_cutoff) then ! non-zero lake_frac dominates over land_frac
-          land_frac(i) = max(0., min(1., 1.-lake_frac(i)))
-        elseif (inland(i) == 1.) then ! land_frac dominates over lake_frac at inland points
-          lake_frac(i) = max(0., min(1., 1.-land_frac(i)))
-        end if
+! include Caspian Sea and Aral Sea if GLDB data set is used, and 
+! exclude lakes in the coastal areas of Antarctica if MODIS data set is used  
+      CALL include_exclude_lakes(lake_frac,land_frac,lake_depth,geolat,geolon,tile_num)
 
 ! epsil is "numerical" nonzero min for lake_frac/land_frac
 ! lake_cutoff/land_cutoff is practical min for lake_frac/land_frac
-        if (min(lake_cutoff,land_cutoff) < epsil) then
-          print *,'lake_cutoff/land_cutoff cannot be smaller than epsil, reset...'
-          lake_cutoff=max(epsil,lake_cutoff)
-          land_cutoff=max(epsil,land_cutoff)
-        end if
+      IF (min(lake_cutoff,land_cutoff) < epsil) then
+        PRINT *,'lake_cutoff/land_cutoff cannot be smaller than epsil, reset...'
+        lake_cutoff=max(epsil,lake_cutoff)
+        land_cutoff=max(epsil,land_cutoff)
+      ENDIF
 
-        if (lake_frac(i)<   lake_cutoff) lake_frac(i)=0.
-        if (lake_frac(i)>1.-lake_cutoff) lake_frac(i)=1.
-        if (inland(i) == 1.) land_frac(i) = 1.-lake_frac(i)
+! adjust land_frac and lake_frac, and make sure land_frac+lake_frac=1 at inland points
+      DO i = 1, tile_sz
         if (land_frac(i)<   land_cutoff) land_frac(i)=0.
         if (land_frac(i)>1.-land_cutoff) land_frac(i)=1.
 
+        if (inland(i) /= 1.) then 
+          lake_frac(i) = 0.
+        endif
+
+        if (lake_frac(i) < lake_cutoff) lake_frac(i)=0.
+        if (lake_frac(i) > 1.-epsil) lake_frac(i)=1.
+      ENDDO
+
+! finalize land_frac/slmsk based on modified lake_frac
+      DO i = 1, tile_sz
+        if (inland(i) == 1.) then
+          land_frac(i) = 1. - lake_frac(i)
+        end if
+
         if (lake_frac(i) < lake_cutoff) then
           lake_depth(i)=0.
         elseif (lake_frac(i) >= lake_cutoff .and. lake_depth(i)==0.) then
@@ -707,6 +772,7 @@ SUBROUTINE write_lakedata_to_orodata(cs_res, cs_lakestat, cs_lakedpth)
         end if
         slmsk(i) = nint(land_frac(i))
       ENDDO
+
 ! write 2 new variables      
       stat = nf90_put_var(ncid, lake_frac_id, lake_frac, &
              start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
@@ -794,43 +860,64 @@ SUBROUTINE write_reg_lakedata_to_orodata(cs_res, tile_x_dim, tile_y_dim, cs_lake
     CALL nc_opchk(stat, "nf90_redef")
 
     IF (nf90_inq_varid(ncid, "lake_frac",lake_frac_id) /= 0) THEN
-    stat = nf90_def_var(ncid,"lake_frac",NF90_FLOAT,dimids,lake_frac_id)
-    CALL nc_opchk(stat, "nf90_def_var: lake_frac")
+      stat = nf90_def_var(ncid,"lake_frac",NF90_FLOAT,dimids,lake_frac_id)
+      CALL nc_opchk(stat, "nf90_def_var: lake_frac")
 #ifdef ADD_ATT_FOR_NEW_VAR
-    stat = nf90_put_att(ncid, lake_frac_id,'coordinates','geolon geolat')
-    CALL nc_opchk(stat, "nf90_put_att: lake_frac:coordinates") 
-    stat = nf90_put_att(ncid, lake_frac_id,'long_name','lake fraction')
-    CALL nc_opchk(stat, "nf90_put_att: lake_frac:long_name") 
-    stat = nf90_put_att(ncid, lake_frac_id,'unit','fraction')
-    CALL nc_opchk(stat, "nf90_put_att: lake_frac:unit") 
-    stat = nf90_put_att(ncid, lake_frac_id,'description', &
-      'based on GLDBv2 (Choulga et al. 2014); missing lakes &
-       added based on land_frac in this dataset.')
-    CALL nc_opchk(stat, "nf90_put_att: lake_frac:description") 
+      stat = nf90_put_att(ncid, lake_frac_id,'coordinates','geolon geolat')
+      CALL nc_opchk(stat, "nf90_put_att: lake_frac:coordinates") 
+      stat = nf90_put_att(ncid, lake_frac_id,'long_name','lake fraction')
+      CALL nc_opchk(stat, "nf90_put_att: lake_frac:long_name") 
+      stat = nf90_put_att(ncid, lake_frac_id,'unit','fraction')
+      CALL nc_opchk(stat, "nf90_put_att: lake_frac:unit") 
+      IF (lakestatus_srce == "GLDBV3") THEN 
+        write(string,'(a,es8.1)') 'based on GLDBv3 (Choulga et al. 2019); missing lakes &
+        added based on land_frac in this dataset; lake_frac cutoff:',lake_cutoff
+      ELSE IF (lakestatus_srce == "GLDBV2") THEN 
+        write(string,'(a,es8.1)') 'based on GLDBv2 (Choulga et al. 2014); missing lakes &
+        added based on land_frac in this dataset; lake_frac cutoff:',lake_cutoff
+      ELSE IF (lakestatus_srce == "MODISP") THEN
+        write(string,'(a,es8.1)') 'based on MODIS (2011-2015) product updated with two &
+        Landsat products: the JRC water product (2016-2020) and the GLC-FCS30 (2020); &
+        the source data set was created by Chengquan Huang of UMD; &  
+        lake_frac cutoff:',lake_cutoff
+      ELSE IF (lakestatus_srce == "VIIRS") THEN
+        write(string,'(a,es8.1)') 'based on multi-year VIIRS global surface type & 
+        classification map (2012-2019); the source data set was created by &
+        Chengquan Huang of UMD and Michael Barlage of NOAA; &
+        lake_frac cutoff:',lake_cutoff
+      ENDIF
+      stat = nf90_put_att(ncid, lake_frac_id,'description',trim(string))
+      CALL nc_opchk(stat, "nf90_put_att: lake_frac:description") 
 #endif
     ENDIF
     IF (nf90_inq_varid(ncid, "lake_depth",lake_depth_id) /= 0) THEN
-    stat = nf90_def_var(ncid,"lake_depth",NF90_FLOAT,dimids,lake_depth_id)
-    CALL nc_opchk(stat, "nf90_def_var: lake_depth")
+      stat = nf90_def_var(ncid,"lake_depth",NF90_FLOAT,dimids,lake_depth_id)
+      CALL nc_opchk(stat, "nf90_def_var: lake_depth")
 #ifdef ADD_ATT_FOR_NEW_VAR
-    stat = nf90_put_att(ncid, lake_depth_id,'coordinates','geolon geolat')
-    CALL nc_opchk(stat, "nf90_put_att: lake_depth:coordinates") 
-    stat = nf90_put_att(ncid, lake_depth_id,'long_name','lake depth')
-    CALL nc_opchk(stat, "nf90_put_att: lake_depth:long_name") 
-    stat = nf90_put_att(ncid, lake_depth_id,'unit','meter')
-    CALL nc_opchk(stat, "nf90_put_att: lake_depth:long_name") 
-    stat = nf90_put_att(ncid, lake_depth_id,'description', &
-      'based on GLDBv2 (Choulga et al. 2014); missing depth set to 10m &
-      (except to 211m in Caspian Sea); spurious large pos. depths are left unchanged.')
-    CALL nc_opchk(stat, "nf90_put_att: lake_depth:description") 
-#endif
+      stat = nf90_put_att(ncid, lake_depth_id,'coordinates','geolon geolat')
+      CALL nc_opchk(stat, "nf90_put_att: lake_depth:coordinates") 
+      stat = nf90_put_att(ncid, lake_depth_id,'long_name','lake depth')
+      CALL nc_opchk(stat, "nf90_put_att: lake_depth:long_name") 
+      stat = nf90_put_att(ncid, lake_depth_id,'unit','meter')
+      CALL nc_opchk(stat, "nf90_put_att: lake_depth:long_name") 
+      IF (lakedepth_srce == "GLDBV3") THEN 
+        stat = nf90_put_att(ncid, lake_depth_id,'description', &
+        'based on GLDBv3 (Choulga et al. 2019); missing depth set to 10m &
+        (except to 211m in Caspian Sea); spurious large pos. depths are left unchanged.')
+      ELSE IF (lakedepth_srce == "GLDBV2") THEN 
+        stat = nf90_put_att(ncid, lake_depth_id,'description', &
+        'based on GLDBv2 (Choulga et al. 2014); missing depth set to 10m &
+        (except to 211m in Caspian Sea); spurious large pos. depths are left unchanged.')
+      ELSE IF (lakedepth_srce == "GLOBATHY") THEN
+        stat = nf90_put_att(ncid, lake_depth_id,'description', &
+        'based on GLOBathy data resampled and projected to the MODIS domain.')
+      ENDIF
+      CALL nc_opchk(stat, "nf90_put_att: lake_depth:description") 
     ENDIF
-    write(string,'(a,es8.1)') 'land_frac and lake_frac are adjusted '// &
-     'such that their sum is 1 at points where inland=1; land_frac '// &
-     'cutoff is ',land_cutoff
+#endif
+    write(string,'(a,es8.1)') 'land_frac and lake_frac are adjusted such that their sum is 1 at points where inland=1; land_frac cutoff is',land_cutoff
     stat = nf90_put_att(ncid, land_frac_id,'description',trim(string))
     CALL nc_opchk(stat, "nf90_put_att: land_frac:description") 
-
     write(string,'(a)') 'slmsk = nint(land_frac)'
     stat = nf90_put_att(ncid, slmsk_id,'description',trim(string))
     CALL nc_opchk(stat, "nf90_put_att: slmsk:description") 
@@ -870,44 +957,36 @@ SUBROUTINE write_reg_lakedata_to_orodata(cs_res, tile_x_dim, tile_y_dim, cs_lake
     lake_frac(:)  = cs_lakestat((tile_num-1)*tile_sz+1:tile_num*tile_sz)
     lake_depth(:) = cs_lakedepth((tile_num-1)*tile_sz+1:tile_num*tile_sz)
 
-! add Caspian Sea and Aral Sea to lake_frac and lake_depth
-    DO i = 1, tile_sz
-      IF (land_frac(i) < 0.9 .AND. lake_frac(i) < 0.1) THEN
-        IF (geolat(i) > 35.0 .AND. geolat(i) <= 50.0 .AND. &
-            geolon(i) > 45.0 .AND. geolon(i) <= 55.0) THEN
-          lake_frac(i) = 1.-land_frac(i)
-          lake_depth(i) = 211.0
-        ENDIF
-        IF (geolat(i) > 35.0 .AND. geolat(i) <= 50.0 .AND. &
-            geolon(i) > 57.0 .AND. geolon(i) <= 63.0) THEN
-          lake_frac(i) = 1.-land_frac(i)
-          lake_depth(i) = 10.0
-        ENDIF
-      ENDIF
-    ENDDO
-
-! adjust land_frac and lake_frac, and make sure land_frac+lake_frac=1 at inland points
-    DO i = 1, tile_sz
-      if (lake_frac(i) >= lake_cutoff) then ! non-zero lake_frac dominates over land_frac
-        land_frac(i) = max(0., min(1., 1.-lake_frac(i)))
-      elseif (inland(i) == 1.) then ! land_frac dominates over lake_frac at inland points
-        lake_frac(i) = max(0., min(1., 1.-land_frac(i)))
-      end if
+! include Caspian Sea and Aral Sea if GLDB data set is used, and 
+! exclude lakes in the coastal areas of Antarctica if MODIS data set is used  
+    CALL include_exclude_lakes(lake_frac,land_frac,lake_depth,geolat,geolon,tile_num)
 
 ! epsil is "numerical" nonzero min for lake_frac/land_frac
 ! lake_cutoff/land_cutoff is practical min for lake_frac/land_frac
-      if (min(lake_cutoff,land_cutoff) < epsil) then
+    IF (min(lake_cutoff,land_cutoff) < epsil) then
         print *,'lake_cutoff/land_cutoff cannot be smaller than epsil, reset...'
         lake_cutoff=max(epsil,lake_cutoff)
         land_cutoff=max(epsil,land_cutoff)
-      end if
+    ENDIF 
 
-      if (lake_frac(i)<   lake_cutoff) lake_frac(i)=0.
-      if (lake_frac(i)>1.-lake_cutoff) lake_frac(i)=1.
-      if (inland(i) == 1.) land_frac(i) = 1.-lake_frac(i)
+! adjust land_frac and lake_frac, and make sure land_frac+lake_frac=1 at inland points
+    DO i = 1, tile_sz
       if (land_frac(i)<   land_cutoff) land_frac(i)=0.
       if (land_frac(i)>1.-land_cutoff) land_frac(i)=1.
 
+      if (inland(i) /= 1.) then
+        lake_frac(i) = 0.
+      endif
+
+      if (lake_frac(i) < lake_cutoff) lake_frac(i)=0.
+      if (lake_frac(i) > 1.-epsil) lake_frac(i)=1.
+    ENDDO
+
+    DO i = 1, tile_sz
+      if (inland(i) == 1.) then
+        land_frac(i) = 1. - lake_frac(i)
+      end if
+
       if (lake_frac(i) < lake_cutoff) then
         lake_depth(i)=0.
       elseif (lake_frac(i) >= lake_cutoff .and. lake_depth(i)==0.) then
@@ -937,12 +1016,66 @@ SUBROUTINE write_reg_lakedata_to_orodata(cs_res, tile_x_dim, tile_y_dim, cs_lake
     stat = nf90_close(ncid)
     CALL nc_opchk(stat, "nf90_close") 
   
-    DEALLOCATE(lake_frac, lake_depth)
-    DEALLOCATE(geolon, geolat)
-    DEALLOCATE(land_frac, slmsk)
-
 END SUBROUTINE write_reg_lakedata_to_orodata
 
+SUBROUTINE include_exclude_lakes(lake_frac,land_frac,lake_depth,geolat,geolon,tile_num)
+    REAL, INTENT(INOUT) :: lake_frac(cs_res*cs_res), lake_depth(cs_res*cs_res)
+    REAL, INTENT(IN) :: land_frac(cs_res*cs_res)
+    REAL, INTENT(IN) :: geolat(cs_res*cs_res), geolon(cs_res*cs_res)
+    INTEGER, INTENT(IN) :: tile_num
+    
+    INTEGER :: tile_sz
+
+    tile_sz = cs_res*cs_res
+! add Caspian Sea and Aral Sea 
+    IF (tile_num == 2 .OR. tile_num == 3 .OR. tile_num == 7) THEN
+      IF (lakedepth_srce == "GLDBV3" .OR. lakedepth_srce == "GLDBV2") THEN 
+        IF (lakestatus_srce == "GLDBV3" .OR. lakestatus_srce == "GLDBV2") THEN
+          DO i = 1, tile_sz
+            IF (land_frac(i) < 0.9 .AND. lake_frac(i) < 0.1) THEN
+              IF (geolat(i) > 35.0 .AND. geolat(i) <= 50.0 .AND. &
+                  geolon(i) > 45.0 .AND. geolon(i) <= 55.0) THEN
+                lake_frac(i) = 1.-land_frac(i)
+                lake_depth(i) = 211.0
+              ENDIF
+              IF (geolat(i) > 35.0 .AND. geolat(i) <= 50.0 .AND. &
+                  geolon(i) > 57.0 .AND. geolon(i) <= 63.0) THEN
+                lake_frac(i) = 1.-land_frac(i)
+                lake_depth(i) = 10.0
+              ENDIF
+            ENDIF
+          ENDDO
+        ENDIF
+        IF (lakestatus_srce == "MODISP" .OR. lakestatus_srce == "VIIRS") THEN
+          DO i = 1, tile_sz
+            IF (land_frac(i) < 0.9) THEN
+              IF (geolat(i) > 35.0 .AND. geolat(i) <= 50.0 .AND. &
+                  geolon(i) > 45.0 .AND. geolon(i) <= 55.0) THEN
+                lake_depth(i) = 211.0
+              ENDIF
+              IF (geolat(i) > 35.0 .AND. geolat(i) <= 50.0 .AND. &
+                  geolon(i) > 57.0 .AND. geolon(i) <= 63.0) THEN
+                lake_depth(i) = 10.0
+              ENDIF
+            ENDIF
+          ENDDO
+        ENDIF
+      ENDIF
+    ENDIF
+! remove lakes below 60 deg south 
+    IF (tile_num == 6) THEN
+      IF (lakestatus_srce == "MODISP" .OR. lakestatus_srce == "VIIRS") THEN 
+        DO i = 1, tile_sz
+          IF (geolat(i) < -60.0) THEN
+            lake_frac(i) = 0.0
+            lake_depth(i) = 0.0
+          ENDIF
+        ENDDO
+      ENDIF
+    ENDIF
+     
+END SUBROUTINE include_exclude_lakes
+
 !> Check NetCDF error code
 !!
 !! @param[in] stat Error code.
diff --git a/ush/fv3gfs_make_lake.sh b/ush/fv3gfs_make_lake.sh
index cfb12aba6..a3a59a3d6 100755
--- a/ush/fv3gfs_make_lake.sh
+++ b/ush/fv3gfs_make_lake.sh
@@ -13,6 +13,31 @@ if [ $gtype != uniform ] && [ $gtype != regional_gfdl ]; then
   echo "lakefrac has only been implemented for 'uniform' and 'regional_gfdl'."
   exit 0
 fi
+echo "lake_data_srce = $lake_data_srce"
+if [ $lake_data_srce == GLDBV3 ]; then
+  lakestatusrc="GLDBV3"
+  lakedepthsrc="GLDBV3"
+fi
+if [ $lake_data_srce == GLDBV2 ]; then
+  lakestatusrc="GLDBV2"
+  lakedepthsrc="GLDBV2"
+fi
+if [ $lake_data_srce == MODISP_GLOBATHY ]; then
+  lakestatusrc="MODISP"
+  lakedepthsrc="GLOBATHY"
+fi
+if [ $lake_data_srce == VIIRS_GLOBATHY ]; then
+  lakestatusrc="VIIRS"
+  lakedepthsrc="GLOBATHY"
+fi
+if [ $lake_data_srce == MODISP_GLDBV3 ]; then
+  lakestatusrc="MODISP"
+  lakedepthsrc="GLDBV3"
+fi
+if [ $lake_data_srce == VIIRS_GLDBV3 ]; then
+  lakestatusrc="VIIRS"
+  lakedepthsrc="GLDBV3"
+fi
 
 exe_add_lake=$exec_dir/lakefrac
 if [ ! -s $exe_add_lake ]; then
@@ -54,7 +79,7 @@ if [ $gtype == uniform ]; then
   done
 fi
 
-if [ $gtype == regional_gfdl ]; then
+if [ $gtype == regional_gfdl ] || [ $gtype == regional_esg ]; then
   tile_beg=7
   tile_end=7
   tile=7
@@ -66,7 +91,7 @@ fi
 
 # create inland mask and save it to the orography files
 
-cutoff=0.99
+cutoff=0.75
 rd=7
 if [ $gtype == uniform ]; then
   $APRUN $exe_inland $res $cutoff $rd g
@@ -74,6 +99,9 @@ fi
 if [ $gtype == regional_gfdl ]; then
   $APRUN $exe_inland $res $cutoff $rd r
 fi
+if [ $gtype == regional_esg ]; then
+  $APRUN $exe_inland $res $cutoff $rd r
+fi
 err=$?
 if [ $err != 0 ]; then
   set +x
@@ -81,12 +109,12 @@ if [ $err != 0 ]; then
   exit $err
 fi
 
-# create lake data for FV3 grid and save it to the orography files
+# create fractional lake data for FV3 grid and save it to the orography files
 
 tile=$tile_beg
 while [ $tile -le $tile_end ]; do
   outfile=oro.C${res}.tile${tile}.nc
-  $APRUN $exe_add_lake ${tile} ${res} ${indir} ${lake_cutoff}
+  $APRUN $exe_add_lake ${tile} ${res} ${indir} ${lakestatusrc} ${lakedepthsrc} ${lake_cutoff}
   err=$?
   if [ $err != 0 ]; then
     set +x

From 7dee97bc3a91289d81022f936594eb9711cb6340 Mon Sep 17 00:00:00 2001
From: "Shan.Sun" <shan.sun@noaa.gov>
Date: Thu, 5 Oct 2023 21:56:28 +0000
Subject: [PATCH 02/15] Modified lakefrac.F90 to ensure compatibility with the
 GNU compiler.

Co-authored-by: Ning Wang ning.wang@noaa.gov
---
 sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 | 22 ++++++++++++--------
 1 file changed, 13 insertions(+), 9 deletions(-)

diff --git a/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 b/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
index f5587ea3d..38f792793 100644
--- a/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
+++ b/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
@@ -51,9 +51,11 @@ PROGRAM lake_frac
     CALL getarg(0, arg) ! get the program name
     IF (iargc() /= 5 .AND. iargc() /= 6) THEN
       PRINT*, 'Usage: ', trim(arg), & 
-       ' [tile_num (0:all tiles, 7:regional)] [resolution (48,96, ...)] [lake data path] [lake status source] [lake depth source]'  
+       ' [tile_num (0:all tiles, 7:regional)] [resolution (48,96, ...)] &
+         [lake data path] [lake status source] [lake depth source]'  
       PRINT*, 'Or: ', trim(arg), & 
-       ' [tile_num (0:all tiles, 7:regional)] [resolution (48,96, ...)] [lake data path] [lake status source] [lake depth source] [lake_cutoff]'
+       ' [tile_num (0:all tiles, 7:regional)] [resolution (48,96, ...)] &
+         [lake data path] [lake status source] [lake depth source] [lake_cutoff]'
       STOP
     ENDIF
     CALL getarg(1, arg)
@@ -294,7 +296,7 @@ SUBROUTINE cal_lake_frac_depth(lakestat,cs_lakestat,lakedpth,cs_lakedpth)
               src_grid_lon_beg,src_grid_lon_end,stride_lon
             PRINT*, 'lon range ', &
               src_grid_lon(src_grid_lon_beg),src_grid_lon(src_grid_lon_end)
-            IF (two_section == .true.) THEN
+            IF (two_section .eqv. .true.) THEN
               PRINT*, 'section1 index lon range and stride', &
                 src_grid_lon_beg1,src_grid_lon_end1,stride_lon
               PRINT*, 'section1 lon range ', &
@@ -306,11 +308,11 @@ SUBROUTINE cal_lake_frac_depth(lakestat,cs_lakestat,lakedpth,cs_lakedpth)
             ENDIF
           ENDIF
 #endif
-          IF (two_section == .false.) THEN
+          IF (two_section .eqv. .false.) THEN
             DO i =  src_grid_lon_beg, src_grid_lon_end, stride_lon
               p(1) = src_grid_lat(j); p(2) = src_grid_lon(i) 
               p(:) = p(:)*d2r
-              IF(enclosure_cnvx(v, 4, p, co_gc) == .true.) THEN
+              IF(enclosure_cnvx(v, 4, p, co_gc) .eqv. .true.) THEN
                 grid_ct = grid_ct+1
                 lkst = lakestat((j-1)*nlon+i); lkdp = lakedpth((j-1)*nlon+i)
                 CALL lake_cell_comp(lkst, lkdp, lake_ct, lake_avg_frac, lake_dpth_sum)
@@ -320,7 +322,7 @@ SUBROUTINE cal_lake_frac_depth(lakestat,cs_lakestat,lakedpth,cs_lakedpth)
             DO i =  src_grid_lon_beg1, src_grid_lon_end1, stride_lon
               p(1) = src_grid_lat(j); p(2) = src_grid_lon(i) 
               p(:) = p(:)*d2r
-              IF(enclosure_cnvx(v, 4, p, co_gc) == .true.) THEN
+              IF(enclosure_cnvx(v, 4, p, co_gc) .eqv. .true.) THEN
                 grid_ct = grid_ct+1
                 lkst = lakestat((j-1)*nlon+i); lkdp = lakedpth((j-1)*nlon+i)
                 CALL lake_cell_comp(lkst, lkdp, lake_ct, lake_avg_frac, lake_dpth_sum)
@@ -329,7 +331,7 @@ SUBROUTINE cal_lake_frac_depth(lakestat,cs_lakestat,lakedpth,cs_lakedpth)
             DO i =  src_grid_lon_beg2, src_grid_lon_end2, stride_lon
               p(1) = src_grid_lat(j); p(2) = src_grid_lon(i) 
               p(:) = p(:)*d2r
-              IF(enclosure_cnvx(v, 4, p, co_gc) == .true.) THEN
+              IF(enclosure_cnvx(v, 4, p, co_gc) .eqv. .true.) THEN
                 grid_ct = grid_ct+1
                 lkst = lakestat((j-1)*nlon+i); lkdp = lakedpth((j-1)*nlon+i)
                 CALL lake_cell_comp(lkst, lkdp, lake_ct, lake_avg_frac, lake_dpth_sum)
@@ -693,7 +695,8 @@ SUBROUTINE write_lakedata_to_orodata(cs_res, cs_lakestat, cs_lakedpth)
       CALL nc_opchk(stat, "nf90_put_att: lake_depth:description") 
 #endif
 
-      write(string,'(a,es8.1)') 'land_frac and lake_frac are adjusted such that their sum is 1 at points where inland=1; land_frac cutoff is',land_cutoff
+      write(string,'(a,es8.1)') 'land_frac and lake_frac are adjusted such that &
+         their sum is 1 at points where inland=1; land_frac cutoff is',land_cutoff
       stat = nf90_put_att(ncid, land_frac_id,'description',trim(string))
       CALL nc_opchk(stat, "nf90_put_att: land_frac:description") 
 
@@ -915,7 +918,8 @@ SUBROUTINE write_reg_lakedata_to_orodata(cs_res, tile_x_dim, tile_y_dim, cs_lake
       CALL nc_opchk(stat, "nf90_put_att: lake_depth:description") 
     ENDIF
 #endif
-    write(string,'(a,es8.1)') 'land_frac and lake_frac are adjusted such that their sum is 1 at points where inland=1; land_frac cutoff is',land_cutoff
+    write(string,'(a,es8.1)') 'land_frac and lake_frac are adjusted such that &
+      their sum is 1 at points where inland=1; land_frac cutoff is',land_cutoff
     stat = nf90_put_att(ncid, land_frac_id,'description',trim(string))
     CALL nc_opchk(stat, "nf90_put_att: land_frac:description") 
     write(string,'(a)') 'slmsk = nint(land_frac)'

From 93f05ce62ad23d72e87d02d5279e4cda5c0e9b47 Mon Sep 17 00:00:00 2001
From: "Shan.Sun" <shan.sun@noaa.gov>
Date: Fri, 6 Oct 2023 21:04:43 +0000
Subject: [PATCH 03/15] Adding doxygen and revising driver script.

---
 driver_scripts/driver_grid.hera.sh           |  3 +--
 sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 | 20 +++++++++++++++++---
 2 files changed, 18 insertions(+), 5 deletions(-)

diff --git a/driver_scripts/driver_grid.hera.sh b/driver_scripts/driver_grid.hera.sh
index a054a3e24..3e38a6521 100755
--- a/driver_scripts/driver_grid.hera.sh
+++ b/driver_scripts/driver_grid.hera.sh
@@ -106,8 +106,7 @@ export soil_type_src="bnu.v3.30s" #  Soil type data.
                                 # For Beijing Norm. Univ. data
                                 # 1) "bnu.v3.30s" for global 30s data.
 
-export reg_domain=na3km   # '3km', 'na3km', '13km', and 'na13km'
-export lake_data_srce=MODISP_GLDBV3 # 'GLDBV3', 'MODIS_GLOBATHY', 'MODISP_GLDBV3', and 'VIIRS_GLDBV3'
+export lake_data_srce=MODISP_GLDBV3 # 'GLDBV3', 'MODISP_GLOBATHY', 'MODISP_GLDBV3', and 'VIIRS_GLDBV3'
 
 if [ $gtype = uniform ]; then
   export res=96
diff --git a/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 b/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
index 38f792793..5953b128c 100644
--- a/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
+++ b/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
@@ -15,7 +15,6 @@
 !!  - Ning Wang, Apr. 2019: Extended the program to process the same lake data 
 !!                         for FV3 stand-alone regional (SAR) model.
 !!     
-!! @return 0 for successful completion and for error.
 !#define DIAG_N_VERBOSE
 #define ADD_ATT_FOR_NEW_VAR
 PROGRAM lake_frac
@@ -39,8 +38,6 @@ PROGRAM lake_frac
     REAL, PARAMETER :: d2r = acos(-1.0) / 180.0
     REAL, PARAMETER :: r2d = 180.0 /acos(-1.0) 
     REAL, PARAMETER :: pi = acos(-1.0) 
-!   REAL*8, PARAMETER :: gppdeg = 119.99444445
-!   REAL*8, PARAMETER :: delta = 1.0 / 119.99444445
     REAL*8, PARAMETER :: gppdeg = 120.0
     REAL*8, PARAMETER :: delta = 1.0 / 120.0
 
@@ -371,6 +368,14 @@ SUBROUTINE cal_lake_frac_depth(lakestat,cs_lakestat,lakedpth,cs_lakedpth)
 
 END SUBROUTINE cal_lake_frac_depth 
 
+!> Compute cumulatively the lake fraction and lake depth for a cell 
+!!
+!! @param[in] lkst and lkdp are lake status and lake depth values from 
+!! a grid point in the source file.
+!! @param[out] lake_ct and lake_avg_frac are lake fraction info accumulated 
+!! for a cell. 
+!! @param[out] lake_depth_sum is the lake depth accumulation for the cell. 
+!! @author Ning Wang
 SUBROUTINE lake_cell_comp(lkst, lkdp, lake_ct, lake_avg_frac, lake_dpth_sum)
     INTEGER*1, INTENT(IN) :: lkst
     INTEGER*2, INTENT(IN) :: lkdp
@@ -1022,6 +1027,15 @@ SUBROUTINE write_reg_lakedata_to_orodata(cs_res, tile_x_dim, tile_y_dim, cs_lake
   
 END SUBROUTINE write_reg_lakedata_to_orodata
 
+!> Include Caspian Sea and Aral Sea if GLDB dataset is used, and 
+!! exclude lakes in the coastal areas of Antarctica if MODIS dataset is used.
+!!
+!! @param[inout] lake_frac, lake_depth, lake fraction and lake depth arrays
+!! of the given tile. 
+!! @param[in] land_frac, land fraction array of the given tile.
+!! @param[in] geolat, geolon, latitude and longitude arrays of the given tile.
+!! @param[in] tile_num, tile number of the given tile.
+!! @author Ning Wang
 SUBROUTINE include_exclude_lakes(lake_frac,land_frac,lake_depth,geolat,geolon,tile_num)
     REAL, INTENT(INOUT) :: lake_frac(cs_res*cs_res), lake_depth(cs_res*cs_res)
     REAL, INTENT(IN) :: land_frac(cs_res*cs_res)

From 2d19ebe0c1c81f12a48e918c25c27b7c3bddb5a7 Mon Sep 17 00:00:00 2001
From: "Shan.Sun" <shan.sun@noaa.gov>
Date: Fri, 13 Oct 2023 20:51:43 +0000
Subject: [PATCH 04/15] Revise doxygen info

---
 sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 | 21 ++++++++++----------
 1 file changed, 11 insertions(+), 10 deletions(-)

diff --git a/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 b/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
index 5953b128c..0b70d50b5 100644
--- a/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
+++ b/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
@@ -370,11 +370,11 @@ END SUBROUTINE cal_lake_frac_depth
 
 !> Compute cumulatively the lake fraction and lake depth for a cell 
 !!
-!! @param[in] lkst and lkdp are lake status and lake depth values from 
-!! a grid point in the source file.
-!! @param[out] lake_ct and lake_avg_frac are lake fraction info accumulated 
-!! for a cell. 
-!! @param[out] lake_depth_sum is the lake depth accumulation for the cell. 
+!! @param[in] lkst lake status value from a grid point in the source data.
+!! @param[in] lkdp lake depth value from a grid point in the source data.
+!! @param[out] lake_ct lake points number accumulated for the cell 
+!! @param[out] lake_avg_frac lake fraction value accumulated for the cell. 
+!! @param[out] lake_dpth_sum is the lake depth value accumulated for the cell. 
 !! @author Ning Wang
 SUBROUTINE lake_cell_comp(lkst, lkdp, lake_ct, lake_avg_frac, lake_dpth_sum)
     INTEGER*1, INTENT(IN) :: lkst
@@ -1030,11 +1030,12 @@ END SUBROUTINE write_reg_lakedata_to_orodata
 !> Include Caspian Sea and Aral Sea if GLDB dataset is used, and 
 !! exclude lakes in the coastal areas of Antarctica if MODIS dataset is used.
 !!
-!! @param[inout] lake_frac, lake_depth, lake fraction and lake depth arrays
-!! of the given tile. 
-!! @param[in] land_frac, land fraction array of the given tile.
-!! @param[in] geolat, geolon, latitude and longitude arrays of the given tile.
-!! @param[in] tile_num, tile number of the given tile.
+!! @param[inout] lake_frac lake fraction array of the given tile. 
+!! @param[inout] lake_depth lake depth array of the given tile. 
+!! @param[in] land_frac land fraction array of the given tile.
+!! @param[in] geolat latitude array of the given tile.
+!! @param[in] geolon longitude array of the given tile.
+!! @param[in] tile_num tile number of the given tile.
 !! @author Ning Wang
 SUBROUTINE include_exclude_lakes(lake_frac,land_frac,lake_depth,geolat,geolon,tile_num)
     REAL, INTENT(INOUT) :: lake_frac(cs_res*cs_res), lake_depth(cs_res*cs_res)

From adc319d608d118ae5e9a731b5b075853d6a712d9 Mon Sep 17 00:00:00 2001
From: "Shan.Sun" <shan.sun@noaa.gov>
Date: Mon, 16 Oct 2023 22:28:37 +0000
Subject: [PATCH 05/15] Addressing doxygen warning.

Co-authored-by: Ning Wang ning.wang@noaa.gov
---
 sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 b/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
index 0b70d50b5..9d9608381 100644
--- a/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
+++ b/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
@@ -15,6 +15,7 @@
 !!  - Ning Wang, Apr. 2019: Extended the program to process the same lake data 
 !!                         for FV3 stand-alone regional (SAR) model.
 !!     
+!! @return 0 for success.
 !#define DIAG_N_VERBOSE
 #define ADD_ATT_FOR_NEW_VAR
 PROGRAM lake_frac
@@ -53,7 +54,7 @@ PROGRAM lake_frac
       PRINT*, 'Or: ', trim(arg), & 
        ' [tile_num (0:all tiles, 7:regional)] [resolution (48,96, ...)] &
          [lake data path] [lake status source] [lake depth source] [lake_cutoff]'
-      STOP
+      STOP -1
     ENDIF
     CALL getarg(1, arg)
     READ(arg,*,iostat=stat) tile_req
@@ -151,7 +152,7 @@ PROGRAM lake_frac
     DEALLOCATE(lakestatus,lakedepth)
     DEALLOCATE(src_grid_lat, src_grid_lon)
 
-    STOP
+    STOP 0
 CONTAINS
 
 !> Calculate lake fraction and depth on the model grid from
@@ -1110,7 +1111,7 @@ SUBROUTINE nc_opchk(stat,opname)
    IF (stat .NE.0)  THEN
      msg = trim(opname) // ' Error, status code and message:'
      PRINT*,trim(msg), stat, nf90_strerror(stat)
-     STOP
+     STOP -1
    END IF
 
 END SUBROUTINE nc_opchk

From 45f59f3c34d6cdaeb7b44fb5a63d5217593d5d8c Mon Sep 17 00:00:00 2001
From: "Shan.Sun" <shan.sun@noaa.gov>
Date: Wed, 18 Oct 2023 22:44:34 +0000
Subject: [PATCH 06/15] Update scripts in driver_scripts for the new datasets.

---
 driver_scripts/driver_grid.jet.sh    | 5 +++--
 driver_scripts/driver_grid.orion.sh  | 5 +++--
 driver_scripts/driver_grid.wcoss2.sh | 5 +++--
 3 files changed, 9 insertions(+), 6 deletions(-)

diff --git a/driver_scripts/driver_grid.jet.sh b/driver_scripts/driver_grid.jet.sh
index fe2d28dc5..aa7d31fbb 100755
--- a/driver_scripts/driver_grid.jet.sh
+++ b/driver_scripts/driver_grid.jet.sh
@@ -103,11 +103,12 @@ export soil_type_src="bnu.v3.30s" # Soil type data.
                                 # 3) "statsgo.conus.30s" for CONUS 30s data
                                 # 4) "statsgo.nh.30s" for NH 30s data
                                 # 5) "statsgo.30s" for global 30s data
+export lake_data_srce=MODISP_GLDBV3 # 'GLDBV3', 'MODISP_GLOBATHY', 'MODISP_GLDBV3', and 'VIIRS_GLDBV3'
 
 if [ $gtype = uniform ]; then
   export res=96
-  export add_lake=false        # Add lake frac and depth to orography data.
-  export lake_cutoff=0.20      # lake frac < lake_cutoff ignored when add_lake=T
+  export add_lake=true         # Add lake frac and depth to orography data.
+  export lake_cutoff=0.50      # lake frac < lake_cutoff ignored when add_lake=T
 elif [ $gtype = stretch ]; then
   export res=96
   export stretch_fac=1.5       # Stretching factor for the grid
diff --git a/driver_scripts/driver_grid.orion.sh b/driver_scripts/driver_grid.orion.sh
index c0190f4d9..1fe9a4bca 100755
--- a/driver_scripts/driver_grid.orion.sh
+++ b/driver_scripts/driver_grid.orion.sh
@@ -102,11 +102,12 @@ export soil_type_src="bnu.v3.30s" #  Soil type data.
                                 # 3) "statsgo.conus.30s" for CONUS 30s data
                                 # 4) "statsgo.nh.30s" for NH 30s data
                                 # 5) "statsgo.30s" for global 30s data
+export lake_data_srce=MODISP_GLDBV3 # 'GLDBV3', 'MODISP_GLOBATHY', 'MODISP_GLDBV3', and 'VIIRS_GLDBV3'
 
 if [ $gtype = uniform ]; then
   export res=96
-  export add_lake=false        # Add lake frac and depth to orography data.
-  export lake_cutoff=0.20      # lake frac < lake_cutoff ignored when add_lake=T
+  export add_lake=true         # Add lake frac and depth to orography data.
+  export lake_cutoff=0.50      # lake frac < lake_cutoff ignored when add_lake=T
 elif [ $gtype = stretch ]; then
   export res=96
   export stretch_fac=1.5       # Stretching factor for the grid
diff --git a/driver_scripts/driver_grid.wcoss2.sh b/driver_scripts/driver_grid.wcoss2.sh
index cbb5ef188..7cb09eb00 100755
--- a/driver_scripts/driver_grid.wcoss2.sh
+++ b/driver_scripts/driver_grid.wcoss2.sh
@@ -101,11 +101,12 @@ export soil_type_src="bnu.v3.30s"  # Soil type data
                                 # 3) "statsgo.conus.30s" for CONUS 30s data
                                 # 4) "statsgo.nh.30s" for NH 30s data
                                 # 5) "statsgo.30s" for global 30s data
+export lake_data_srce=MODISP_GLDBV3 # 'GLDBV3', 'MODISP_GLOBATHY', 'MODISP_GLDBV3', and 'VIIRS_GLDBV3'
 
 if [ $gtype = uniform ]; then
   export res=96
-  export add_lake=false        # Add lake frac and depth to orography data.
-  export lake_cutoff=0.20      # lake frac < lake_cutoff ignored when add_lake=T
+  export add_lake=true         # Add lake frac and depth to orography data.
+  export lake_cutoff=0.50      # lake frac < lake_cutoff ignored when add_lake=T
 elif [ $gtype = stretch ]; then
   export res=96
   export stretch_fac=1.5       # Stretching factor for the grid

From 7a680cbf0f4e02e201788698de34ad326024881f Mon Sep 17 00:00:00 2001
From: "Shan.Sun" <shan.sun@noaa.gov>
Date: Thu, 19 Oct 2023 16:53:49 +0000
Subject: [PATCH 07/15] Adding comments

---
 driver_scripts/driver_grid.hera.sh   | 4 +++-
 driver_scripts/driver_grid.jet.sh    | 5 ++++-
 driver_scripts/driver_grid.orion.sh  | 5 ++++-
 driver_scripts/driver_grid.wcoss2.sh | 5 ++++-
 4 files changed, 15 insertions(+), 4 deletions(-)

diff --git a/driver_scripts/driver_grid.hera.sh b/driver_scripts/driver_grid.hera.sh
index 3e38a6521..03b3bd31b 100755
--- a/driver_scripts/driver_grid.hera.sh
+++ b/driver_scripts/driver_grid.hera.sh
@@ -106,7 +106,9 @@ export soil_type_src="bnu.v3.30s" #  Soil type data.
                                 # For Beijing Norm. Univ. data
                                 # 1) "bnu.v3.30s" for global 30s data.
 
-export lake_data_srce=MODISP_GLDBV3 # 'GLDBV3', 'MODISP_GLOBATHY', 'MODISP_GLDBV3', and 'VIIRS_GLDBV3'
+# choose dataset sources for lakefrac & lakedepth so that lake_data_srce=LakeFrac_LakeDepth; 
+# available options are 'MODISP_GLDBV3', 'MODISP_GLOBATHY' & 'VIIRS_GLDBV3';
+export lake_data_srce=MODISP_GLDBV3 
 
 if [ $gtype = uniform ]; then
   export res=96
diff --git a/driver_scripts/driver_grid.jet.sh b/driver_scripts/driver_grid.jet.sh
index aa7d31fbb..1dafc5d00 100755
--- a/driver_scripts/driver_grid.jet.sh
+++ b/driver_scripts/driver_grid.jet.sh
@@ -103,7 +103,10 @@ export soil_type_src="bnu.v3.30s" # Soil type data.
                                 # 3) "statsgo.conus.30s" for CONUS 30s data
                                 # 4) "statsgo.nh.30s" for NH 30s data
                                 # 5) "statsgo.30s" for global 30s data
-export lake_data_srce=MODISP_GLDBV3 # 'GLDBV3', 'MODISP_GLOBATHY', 'MODISP_GLDBV3', and 'VIIRS_GLDBV3'
+
+# choose dataset sources for lakefrac & lakedepth so that lake_data_srce=LakeFrac_LakeDepth; 
+# available options are 'MODISP_GLDBV3', 'MODISP_GLOBATHY' & 'VIIRS_GLDBV3';
+export lake_data_srce=MODISP_GLDBV3 
 
 if [ $gtype = uniform ]; then
   export res=96
diff --git a/driver_scripts/driver_grid.orion.sh b/driver_scripts/driver_grid.orion.sh
index 1fe9a4bca..605c89aab 100755
--- a/driver_scripts/driver_grid.orion.sh
+++ b/driver_scripts/driver_grid.orion.sh
@@ -102,7 +102,10 @@ export soil_type_src="bnu.v3.30s" #  Soil type data.
                                 # 3) "statsgo.conus.30s" for CONUS 30s data
                                 # 4) "statsgo.nh.30s" for NH 30s data
                                 # 5) "statsgo.30s" for global 30s data
-export lake_data_srce=MODISP_GLDBV3 # 'GLDBV3', 'MODISP_GLOBATHY', 'MODISP_GLDBV3', and 'VIIRS_GLDBV3'
+
+# choose dataset sources for lakefrac & lakedepth so that lake_data_srce=LakeFrac_LakeDepth; 
+# available options are 'MODISP_GLDBV3', 'MODISP_GLOBATHY' & 'VIIRS_GLDBV3';
+export lake_data_srce=MODISP_GLDBV3 
 
 if [ $gtype = uniform ]; then
   export res=96
diff --git a/driver_scripts/driver_grid.wcoss2.sh b/driver_scripts/driver_grid.wcoss2.sh
index 7cb09eb00..58a689b03 100755
--- a/driver_scripts/driver_grid.wcoss2.sh
+++ b/driver_scripts/driver_grid.wcoss2.sh
@@ -101,7 +101,10 @@ export soil_type_src="bnu.v3.30s"  # Soil type data
                                 # 3) "statsgo.conus.30s" for CONUS 30s data
                                 # 4) "statsgo.nh.30s" for NH 30s data
                                 # 5) "statsgo.30s" for global 30s data
-export lake_data_srce=MODISP_GLDBV3 # 'GLDBV3', 'MODISP_GLOBATHY', 'MODISP_GLDBV3', and 'VIIRS_GLDBV3'
+
+# choose dataset sources for lakefrac & lakedepth so that lake_data_srce=LakeFrac_LakeDepth; 
+# available options are 'MODISP_GLDBV3', 'MODISP_GLOBATHY' & 'VIIRS_GLDBV3';
+export lake_data_srce=MODISP_GLDBV3 
 
 if [ $gtype = uniform ]; then
   export res=96

From 73b193a470d4abd88c40cb12ad9608538754d4aa Mon Sep 17 00:00:00 2001
From: "Shan.Sun" <shan.sun@noaa.gov>
Date: Fri, 20 Oct 2023 04:07:09 +0000
Subject: [PATCH 08/15] Clean up scripts

---
 driver_scripts/driver_grid.hera.sh   |  2 +-
 driver_scripts/driver_grid.jet.sh    |  2 +-
 driver_scripts/driver_grid.orion.sh  |  2 +-
 driver_scripts/driver_grid.wcoss2.sh |  2 +-
 ush/fv3gfs_make_lake.sh              | 20 ++++++++------------
 5 files changed, 12 insertions(+), 16 deletions(-)

diff --git a/driver_scripts/driver_grid.hera.sh b/driver_scripts/driver_grid.hera.sh
index 03b3bd31b..85dcb6360 100755
--- a/driver_scripts/driver_grid.hera.sh
+++ b/driver_scripts/driver_grid.hera.sh
@@ -107,7 +107,7 @@ export soil_type_src="bnu.v3.30s" #  Soil type data.
                                 # 1) "bnu.v3.30s" for global 30s data.
 
 # choose dataset sources for lakefrac & lakedepth so that lake_data_srce=LakeFrac_LakeDepth; 
-# available options are 'MODISP_GLDBV3', 'MODISP_GLOBATHY' & 'VIIRS_GLDBV3';
+# available options are 'MODISP_GLDBV3', 'MODISP_GLOBATHY', 'VIIRS_GLDBV3', 'VIIRS_GLOBATHY' & 'GLDBV3'
 export lake_data_srce=MODISP_GLDBV3 
 
 if [ $gtype = uniform ]; then
diff --git a/driver_scripts/driver_grid.jet.sh b/driver_scripts/driver_grid.jet.sh
index 1dafc5d00..0fe8b113e 100755
--- a/driver_scripts/driver_grid.jet.sh
+++ b/driver_scripts/driver_grid.jet.sh
@@ -105,7 +105,7 @@ export soil_type_src="bnu.v3.30s" # Soil type data.
                                 # 5) "statsgo.30s" for global 30s data
 
 # choose dataset sources for lakefrac & lakedepth so that lake_data_srce=LakeFrac_LakeDepth; 
-# available options are 'MODISP_GLDBV3', 'MODISP_GLOBATHY' & 'VIIRS_GLDBV3';
+# available options are 'MODISP_GLDBV3', 'MODISP_GLOBATHY', 'VIIRS_GLDBV3', 'VIIRS_GLOBATHY' & 'GLDBV3'
 export lake_data_srce=MODISP_GLDBV3 
 
 if [ $gtype = uniform ]; then
diff --git a/driver_scripts/driver_grid.orion.sh b/driver_scripts/driver_grid.orion.sh
index 605c89aab..0ef5240a7 100755
--- a/driver_scripts/driver_grid.orion.sh
+++ b/driver_scripts/driver_grid.orion.sh
@@ -104,7 +104,7 @@ export soil_type_src="bnu.v3.30s" #  Soil type data.
                                 # 5) "statsgo.30s" for global 30s data
 
 # choose dataset sources for lakefrac & lakedepth so that lake_data_srce=LakeFrac_LakeDepth; 
-# available options are 'MODISP_GLDBV3', 'MODISP_GLOBATHY' & 'VIIRS_GLDBV3';
+# available options are 'MODISP_GLDBV3', 'MODISP_GLOBATHY', 'VIIRS_GLDBV3', 'VIIRS_GLOBATHY' & 'GLDBV3'
 export lake_data_srce=MODISP_GLDBV3 
 
 if [ $gtype = uniform ]; then
diff --git a/driver_scripts/driver_grid.wcoss2.sh b/driver_scripts/driver_grid.wcoss2.sh
index 58a689b03..74e829130 100755
--- a/driver_scripts/driver_grid.wcoss2.sh
+++ b/driver_scripts/driver_grid.wcoss2.sh
@@ -103,7 +103,7 @@ export soil_type_src="bnu.v3.30s"  # Soil type data
                                 # 5) "statsgo.30s" for global 30s data
 
 # choose dataset sources for lakefrac & lakedepth so that lake_data_srce=LakeFrac_LakeDepth; 
-# available options are 'MODISP_GLDBV3', 'MODISP_GLOBATHY' & 'VIIRS_GLDBV3';
+# available options are 'MODISP_GLDBV3', 'MODISP_GLOBATHY', 'VIIRS_GLDBV3', 'VIIRS_GLOBATHY' & 'GLDBV3'
 export lake_data_srce=MODISP_GLDBV3 
 
 if [ $gtype = uniform ]; then
diff --git a/ush/fv3gfs_make_lake.sh b/ush/fv3gfs_make_lake.sh
index a3a59a3d6..9d8d784d7 100755
--- a/ush/fv3gfs_make_lake.sh
+++ b/ush/fv3gfs_make_lake.sh
@@ -14,28 +14,24 @@ if [ $gtype != uniform ] && [ $gtype != regional_gfdl ]; then
   exit 0
 fi
 echo "lake_data_srce = $lake_data_srce"
-if [ $lake_data_srce == GLDBV3 ]; then
-  lakestatusrc="GLDBV3"
+if [ $lake_data_srce == MODISP_GLDBV3 ]; then
+  lakestatusrc="MODISP"
   lakedepthsrc="GLDBV3"
 fi
-if [ $lake_data_srce == GLDBV2 ]; then
-  lakestatusrc="GLDBV2"
-  lakedepthsrc="GLDBV2"
-fi
 if [ $lake_data_srce == MODISP_GLOBATHY ]; then
   lakestatusrc="MODISP"
   lakedepthsrc="GLOBATHY"
 fi
-if [ $lake_data_srce == VIIRS_GLOBATHY ]; then
+if [ $lake_data_srce == VIIRS_GLDBV3 ]; then
   lakestatusrc="VIIRS"
-  lakedepthsrc="GLOBATHY"
-fi
-if [ $lake_data_srce == MODISP_GLDBV3 ]; then
-  lakestatusrc="MODISP"
   lakedepthsrc="GLDBV3"
 fi
-if [ $lake_data_srce == VIIRS_GLDBV3 ]; then
+if [ $lake_data_srce == VIIRS_GLOBATHY ]; then
   lakestatusrc="VIIRS"
+  lakedepthsrc="GLOBATHY"
+fi
+if [ $lake_data_srce == GLDBV3 ]; then
+  lakestatusrc="GLDBV3"
   lakedepthsrc="GLDBV3"
 fi
 

From 3e2808d4f19769a9c233ddb31d3e44fb9f3e935b Mon Sep 17 00:00:00 2001
From: "Shan.Sun" <shan.sun@noaa.gov>
Date: Fri, 20 Oct 2023 20:41:39 +0000
Subject: [PATCH 09/15] Adding add_lake=true to c96.uniform.sh

---
 reg_tests/grid_gen/c96.uniform.sh | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/reg_tests/grid_gen/c96.uniform.sh b/reg_tests/grid_gen/c96.uniform.sh
index 6811c71af..d4f569d51 100755
--- a/reg_tests/grid_gen/c96.uniform.sh
+++ b/reg_tests/grid_gen/c96.uniform.sh
@@ -1,9 +1,9 @@
 #!/bin/bash
 
 #-----------------------------------------------------------------------
-# Create a C96 global uniform grid.  Compare output to a set
-# of baseline files using the 'nccmp' utility.  This script is
-# run by the machine specific driver script.
+# Create a C96 global uniform grid including inland lakes. Compare 
+# output to a set of baseline files using the 'nccmp' utility.  This
+# script is run by the machine specific driver script.
 #-----------------------------------------------------------------------
 
 set -x
@@ -13,6 +13,8 @@ export out_dir=${WORK_DIR}/c96.uniform
 
 export res=96
 export gtype=uniform
+export add_lake=true
+export lake_data_srce=MODISP_GLDBV3
 
 NCCMP=${NCCMP:-$(which nccmp)}
 

From 3754527be29acdfc446d38f462cf0014782e6e7f Mon Sep 17 00:00:00 2001
From: "Shan.Sun" <shan.sun@noaa.gov>
Date: Tue, 24 Oct 2023 15:04:42 +0000
Subject: [PATCH 10/15] Increase wall clock time

---
 reg_tests/grid_gen/driver.hera.sh   | 2 +-
 reg_tests/grid_gen/driver.jet.sh    | 2 +-
 reg_tests/grid_gen/driver.wcoss2.sh | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/reg_tests/grid_gen/driver.hera.sh b/reg_tests/grid_gen/driver.hera.sh
index c0211cc73..66673ab34 100755
--- a/reg_tests/grid_gen/driver.hera.sh
+++ b/reg_tests/grid_gen/driver.hera.sh
@@ -66,7 +66,7 @@ export OMP_NUM_THREADS=24
 #-----------------------------------------------------------------------------
 
 LOG_FILE1=${LOG_FILE}01
-TEST1=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:10:00 -A $PROJECT_CODE -q $QUEUE -J c96.uniform \
+TEST1=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:15:00 -A $PROJECT_CODE -q $QUEUE -J c96.uniform \
       -o $LOG_FILE1 -e $LOG_FILE1 ./c96.uniform.sh)
 
 #-----------------------------------------------------------------------------
diff --git a/reg_tests/grid_gen/driver.jet.sh b/reg_tests/grid_gen/driver.jet.sh
index 6641098bc..db387c9d4 100755
--- a/reg_tests/grid_gen/driver.jet.sh
+++ b/reg_tests/grid_gen/driver.jet.sh
@@ -64,7 +64,7 @@ export OMP_NUM_THREADS=24
 #-----------------------------------------------------------------------------
 
 LOG_FILE1=${LOG_FILE}01
-TEST1=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:10:00 -A $PROJECT_CODE -q $QUEUE -J c96.uniform \
+TEST1=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:20:00 -A $PROJECT_CODE -q $QUEUE -J c96.uniform \
       --partition=xjet -o $LOG_FILE1 -e $LOG_FILE1 ./c96.uniform.sh)
 
 #-----------------------------------------------------------------------------
diff --git a/reg_tests/grid_gen/driver.wcoss2.sh b/reg_tests/grid_gen/driver.wcoss2.sh
index 398765e1c..35c077f34 100755
--- a/reg_tests/grid_gen/driver.wcoss2.sh
+++ b/reg_tests/grid_gen/driver.wcoss2.sh
@@ -67,7 +67,7 @@ rm -fr $WORK_DIR
 #-----------------------------------------------------------------------------
 
 LOG_FILE1=${LOG_FILE}01
-TEST1=$(qsub -V -o $LOG_FILE1 -e $LOG_FILE1 -q $QUEUE -A $PROJECT_CODE -l walltime=00:10:00 \
+TEST1=$(qsub -V -o $LOG_FILE1 -e $LOG_FILE1 -q $QUEUE -A $PROJECT_CODE -l walltime=00:15:00 \
         -N c96.uniform -l select=1:ncpus=30:mem=40GB $PWD/c96.uniform.sh)
 
 #-----------------------------------------------------------------------------

From abc106d381fcd2680ba02361e06d5bf54799e310 Mon Sep 17 00:00:00 2001
From: George Gayno <George.Gayno@noaa.gov>
Date: Wed, 1 Nov 2023 18:44:35 +0000
Subject: [PATCH 11/15] Update 'readthedocs' for new input lake data.

Fixes #854.
---
 docs/source/ufs_utils.rst | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/docs/source/ufs_utils.rst b/docs/source/ufs_utils.rst
index 314c179b9..7810f3817 100644
--- a/docs/source/ufs_utils.rst
+++ b/docs/source/ufs_utils.rst
@@ -375,8 +375,13 @@ Program inputs and outputs
 
       * grid file - the "grid" file from the make_hgrid or regional_esg programs  - CRES_grid.tile#.nc - (NetCDF)
       * orography file - the orography file including the 'inland' flag record from the inland program - oro.CRES.tile#.nc (NetCDF)
-      * lake status code file - GlobalLakeStatus.dat (located in `./fix/fix_orog <https://noaa-ufs-srw-pds.s3.amazonaws.com/index.html#fix/fix_orog/>`_). See GlobalLakeStatus.txt for the defintion of each code.
-      * lake depth file - GlobalLakeDepth.dat (located in `./fix/fix_orog <https://noaa-ufs-srw-pds.s3.amazonaws.com/index.html#fix/fix_orog/>`_). See GlobalLakeDepth.txt for a description of this file.
+      * lake status code file - One of the following files. (located in `./fix/fix_orog <https://noaa-ufs-srw-pds.s3.amazonaws.com/index.html#fix/fix_orog/>`_). See GlobalLakeStatus.txt for a description of the file format.
+          * GlobalLakeStatus_MOSISp.dat
+          * GlobalLakeStatus_GLDBv3release.dat
+          * GlobalLakeStatus_VIIRS.dat
+      * lake depth file - One of the following files. (located in `./fix/fix_orog <https://noaa-ufs-srw-pds.s3.amazonaws.com/index.html#fix/fix_orog/>`_). See GlobalLakeDepth.txt for a description of this file.
+          * GlobalLakeDepth_GLDBv3release.dat
+          * GlobalLakeDepth_GLOBathy.dat
 
 **Output data:**
 

From 4925caec5ce61f9b04d52d182327b823a639f309 Mon Sep 17 00:00:00 2001
From: "Shan.Sun" <shan.sun@noaa.gov>
Date: Thu, 2 Nov 2023 03:59:56 +0000
Subject: [PATCH 12/15] Introducing 'binary_lake': when binary_lake=1 or 0,
 grid points with '1 > lake_frac >= lake_cutoff' to have lake_frac= 1 or 0;

---
 driver_scripts/driver_grid.hera.sh           |  3 +-
 driver_scripts/driver_grid.jet.sh            |  3 +-
 driver_scripts/driver_grid.orion.sh          |  3 +-
 driver_scripts/driver_grid.wcoss2.sh         |  3 +-
 sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 | 35 +++++++++++---------
 ush/fv3gfs_driver_grid.sh                    |  6 ++--
 ush/fv3gfs_make_lake.sh                      |  2 +-
 7 files changed, 32 insertions(+), 23 deletions(-)

diff --git a/driver_scripts/driver_grid.hera.sh b/driver_scripts/driver_grid.hera.sh
index 85dcb6360..e78db425b 100755
--- a/driver_scripts/driver_grid.hera.sh
+++ b/driver_scripts/driver_grid.hera.sh
@@ -113,7 +113,8 @@ export lake_data_srce=MODISP_GLDBV3
 if [ $gtype = uniform ]; then
   export res=96
   export add_lake=true         # Add lake frac and depth to orography data.
-  export lake_cutoff=0.50      # lake frac < lake_cutoff ignored when add_lake=T
+  export lake_cutoff=0.50      # ignore lake_frac < lake_cutoff when add_lake=T
+  export binary_lake=1         # return '1 > lake_frac >= lake_cutoff' as 1/0 when binary_lake=1/0
 elif [ $gtype = stretch ]; then
   export res=96
   export stretch_fac=1.5       # Stretching factor for the grid
diff --git a/driver_scripts/driver_grid.jet.sh b/driver_scripts/driver_grid.jet.sh
index 0fe8b113e..c8ebc51f4 100755
--- a/driver_scripts/driver_grid.jet.sh
+++ b/driver_scripts/driver_grid.jet.sh
@@ -111,7 +111,8 @@ export lake_data_srce=MODISP_GLDBV3
 if [ $gtype = uniform ]; then
   export res=96
   export add_lake=true         # Add lake frac and depth to orography data.
-  export lake_cutoff=0.50      # lake frac < lake_cutoff ignored when add_lake=T
+  export lake_cutoff=0.50      # ignore lake_frac < lake_cutoff when add_lake=T
+  export binary_lake=1         # return '1 > lake_frac >= lake_cutoff' as 1/0 when binary_lake=1/0
 elif [ $gtype = stretch ]; then
   export res=96
   export stretch_fac=1.5       # Stretching factor for the grid
diff --git a/driver_scripts/driver_grid.orion.sh b/driver_scripts/driver_grid.orion.sh
index 0ef5240a7..bd12e944e 100755
--- a/driver_scripts/driver_grid.orion.sh
+++ b/driver_scripts/driver_grid.orion.sh
@@ -110,7 +110,8 @@ export lake_data_srce=MODISP_GLDBV3
 if [ $gtype = uniform ]; then
   export res=96
   export add_lake=true         # Add lake frac and depth to orography data.
-  export lake_cutoff=0.50      # lake frac < lake_cutoff ignored when add_lake=T
+  export lake_cutoff=0.50      # ignore lake_frac < lake_cutoff when add_lake=T
+  export binary_lake=1         # return '1 > lake_frac >= lake_cutoff' as 1/0 when binary_lake=1/0
 elif [ $gtype = stretch ]; then
   export res=96
   export stretch_fac=1.5       # Stretching factor for the grid
diff --git a/driver_scripts/driver_grid.wcoss2.sh b/driver_scripts/driver_grid.wcoss2.sh
index 74e829130..4d55efef5 100755
--- a/driver_scripts/driver_grid.wcoss2.sh
+++ b/driver_scripts/driver_grid.wcoss2.sh
@@ -109,7 +109,8 @@ export lake_data_srce=MODISP_GLDBV3
 if [ $gtype = uniform ]; then
   export res=96
   export add_lake=true         # Add lake frac and depth to orography data.
-  export lake_cutoff=0.50      # lake frac < lake_cutoff ignored when add_lake=T
+  export lake_cutoff=0.50      # ignore lake_frac < lake_cutoff when add_lake=T
+  export binary_lake=1         # return '1 > lake_frac >= lake_cutoff' as 1/0 when binary_lake=1/0
 elif [ $gtype = stretch ]; then
   export res=96
   export stretch_fac=1.5       # Stretching factor for the grid
diff --git a/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 b/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
index 9d9608381..161213297 100644
--- a/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
+++ b/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
@@ -32,7 +32,7 @@ PROGRAM lake_frac
     REAL, ALLOCATABLE :: cs_lakestatus(:), cs_lakedepth(:)
     REAL, ALLOCATABLE :: src_grid_lon(:), src_grid_lat(:)
 
-    INTEGER :: tile_req, tile_beg, tile_end
+    INTEGER :: tile_req, tile_beg, tile_end, binary_lake
     REAL :: lake_cutoff
 
     INTEGER, PARAMETER :: nlat = 21600, nlon = 43200
@@ -47,13 +47,13 @@ PROGRAM lake_frac
     INTEGER :: stat
     
     CALL getarg(0, arg) ! get the program name
-    IF (iargc() /= 5 .AND. iargc() /= 6) THEN
+    IF (iargc() /= 5 .AND. iargc() /= 7) THEN
       PRINT*, 'Usage: ', trim(arg), & 
        ' [tile_num (0:all tiles, 7:regional)] [resolution (48,96, ...)] &
          [lake data path] [lake status source] [lake depth source]'  
       PRINT*, 'Or: ', trim(arg), & 
        ' [tile_num (0:all tiles, 7:regional)] [resolution (48,96, ...)] &
-         [lake data path] [lake status source] [lake depth source] [lake_cutoff]'
+         [lake data path] [lake status source] [lake depth source] [lake_cutoff] [binary_lake]'
       STOP -1
     ENDIF
     CALL getarg(1, arg)
@@ -66,14 +66,18 @@ PROGRAM lake_frac
 
     IF (iargc() == 5) THEN
       lake_cutoff = 0.20
+      binary_lake = 0
     ELSE
       CALL getarg(6, arg)
       READ(arg,*,iostat=stat) lake_cutoff
+      CALL getarg(7, arg)
+      READ(arg,*,iostat=stat) binary_lake
     ENDIF
 
     PRINT*, 'lake status source:', trim(lakestatus_srce) 
     PRINT*, 'lake depth source:', trim(lakedepth_srce) 
     PRINT*, 'lake cutoff:', lake_cutoff 
+    PRINT*, 'binary lake:', binary_lake
 
     IF (tile_req == 0) THEN
       tile_beg = 1; tile_end = 6
@@ -612,7 +616,7 @@ SUBROUTINE write_lakedata_to_orodata(cs_res, cs_lakestat, cs_lakedpth)
     INTEGER :: stat, ncid, x_dimid, y_dimid, varid, dimids(2)
     INTEGER :: lake_frac_id, lake_depth_id
     INTEGER :: land_frac_id, slmsk_id, inland_id, geolon_id, geolat_id
-    CHARACTER(len=256) :: filename,string
+    CHARACTER(len=256) :: filename,string,lakeinfo
     CHARACTER(len=1) :: ich
     CHARACTER(len=4) res_ch
     REAL :: lake_frac(cs_res*cs_res),lake_depth(cs_res*cs_res)
@@ -657,22 +661,19 @@ SUBROUTINE write_lakedata_to_orodata(cs_res, cs_lakestat, cs_lakedpth)
       CALL nc_opchk(stat, "nf90_put_att: lake_frac:long_name") 
       stat = nf90_put_att(ncid, lake_frac_id,'unit','fraction')
       CALL nc_opchk(stat, "nf90_put_att: lake_frac:unit") 
+      write(lakeinfo,'(a,f4.2,a,i1)') ' lake_frac cutoff=',lake_cutoff,'; binary_lake=',binary_lake
       IF (lakestatus_srce == "GLDBV3") THEN 
-        write(string,'(a,es8.1)') 'based on GLDBv3 (Choulga et al. 2019); missing lakes &
-        added based on land_frac in this dataset; lake_frac cutoff:',lake_cutoff
+        write(string,'(2a)') 'based on GLDBv3 (Choulga et al. 2019); missing lakes & added based on land_frac in this dataset;',trim(lakeinfo)
       ELSE IF (lakestatus_srce == "GLDBV2") THEN 
-        write(string,'(a,es8.1)') 'based on GLDBv2 (Choulga et al. 2014); missing lakes &
-        added based on land_frac in this dataset; lake_frac cutoff:',lake_cutoff
+        write(string,'(2a)') 'based on GLDBv2 (Choulga et al. 2014); missing lakes & added based on land_frac in this dataset;',trim(lakeinfo)
       ELSE IF (lakestatus_srce == "MODISP") THEN
-        write(string,'(a,es8.1)') 'based on MODIS (2011-2015) product updated with two &
+        write(string,'(2a)') 'based on MODIS (2011-2015) product updated with two &
         Landsat products: the JRC water product (2016-2020) and the GLC-FCS30 (2020); &
-        the source data set was created by Chengquan Huang of UMD; &  
-        lake_frac cutoff:',lake_cutoff
+        the source data set was created by Chengquan Huang of UMD;',trim(lakeinfo)
       ELSE IF (lakestatus_srce == "VIIRS") THEN
-        write(string,'(a,es8.1)') 'based on multi-year VIIRS global surface type & 
+        write(string,'(a,f4.2,a,i1)') 'based on multi-year VIIRS global surface type & 
         classification map (2012-2019); the source data set was created by &
-        Chengquan Huang of UMD and Michael Barlage of NOAA; &
-        lake_frac cutoff:',lake_cutoff
+        Chengquan Huang of UMD and Michael Barlage of NOAA;',trim(lakeinfo)
       ENDIF
       stat = nf90_put_att(ncid, lake_frac_id,'description',trim(string))
       CALL nc_opchk(stat, "nf90_put_att: lake_frac:description") 
@@ -764,7 +765,11 @@ SUBROUTINE write_lakedata_to_orodata(cs_res, cs_lakestat, cs_lakedpth)
           lake_frac(i) = 0.
         endif
 
-        if (lake_frac(i) < lake_cutoff) lake_frac(i)=0.
+        if (lake_frac(i) < lake_cutoff) then
+           lake_frac(i)=0.
+        elseif (binary_lake == 1) then
+           lake_frac(i)=1.
+        end if
         if (lake_frac(i) > 1.-epsil) lake_frac(i)=1.
       ENDDO
 
diff --git a/ush/fv3gfs_driver_grid.sh b/ush/fv3gfs_driver_grid.sh
index 9347749f0..c0d456ac7 100755
--- a/ush/fv3gfs_driver_grid.sh
+++ b/ush/fv3gfs_driver_grid.sh
@@ -47,9 +47,9 @@ export res=${res:-96}           # resolution of tile: 48, 96, 128, 192, 384, 768
 export gtype=${gtype:-uniform}  # grid type: uniform, stretch, nest, regional_gfdl,
                                 # or regional_esg
 
-export add_lake=${add_lake:-false}      # add lake fraction and depth.  uniform only.
-export lake_cutoff=${lake_cutoff:-0.20} # lake fractions less than lake_cutoff 
-                                        # are ignored.
+export add_lake=${add_lake:-false}      # add lake fraction and depth. uniform only.
+export lake_cutoff=${lake_cutoff:-0.20} # ignore lake_frac < lake_cutoff when add_lake=T
+export binary_lake=${binary_lake:0}     # return '1 > lake_frac >= lake_cutoff' as 1/0 when binary_lake=1/0
 
 export make_gsl_orog=${make_gsl_orog:-false} # when true, create GSL drag suite orog files.
 
diff --git a/ush/fv3gfs_make_lake.sh b/ush/fv3gfs_make_lake.sh
index 9d8d784d7..f8f293587 100755
--- a/ush/fv3gfs_make_lake.sh
+++ b/ush/fv3gfs_make_lake.sh
@@ -110,7 +110,7 @@ fi
 tile=$tile_beg
 while [ $tile -le $tile_end ]; do
   outfile=oro.C${res}.tile${tile}.nc
-  $APRUN $exe_add_lake ${tile} ${res} ${indir} ${lakestatusrc} ${lakedepthsrc} ${lake_cutoff}
+  $APRUN $exe_add_lake ${tile} ${res} ${indir} ${lakestatusrc} ${lakedepthsrc} ${lake_cutoff} ${binary_lake}
   err=$?
   if [ $err != 0 ]; then
     set +x

From 4a7eeb31b2b4c2655aa1bd8b0fbd115c3a3e9c54 Mon Sep 17 00:00:00 2001
From: "Shan.Sun" <shan.sun@noaa.gov>
Date: Thu, 2 Nov 2023 04:41:25 +0000
Subject: [PATCH 13/15] Fixed a typo

---
 sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 b/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
index 161213297..7816e21ad 100644
--- a/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
+++ b/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
@@ -671,7 +671,7 @@ SUBROUTINE write_lakedata_to_orodata(cs_res, cs_lakestat, cs_lakedpth)
         Landsat products: the JRC water product (2016-2020) and the GLC-FCS30 (2020); &
         the source data set was created by Chengquan Huang of UMD;',trim(lakeinfo)
       ELSE IF (lakestatus_srce == "VIIRS") THEN
-        write(string,'(a,f4.2,a,i1)') 'based on multi-year VIIRS global surface type & 
+        write(string,'(2a)') 'based on multi-year VIIRS global surface type & 
         classification map (2012-2019); the source data set was created by &
         Chengquan Huang of UMD and Michael Barlage of NOAA;',trim(lakeinfo)
       ENDIF

From adc79a390942e0694acc71f747e11383dee5a218 Mon Sep 17 00:00:00 2001
From: "Shan.Sun" <shan.sun@noaa.gov>
Date: Thu, 2 Nov 2023 15:04:54 +0000
Subject: [PATCH 14/15] Improving the previous commit: change default values
 and comments

---
 driver_scripts/driver_grid.hera.sh           | 4 ++--
 driver_scripts/driver_grid.jet.sh            | 4 ++--
 driver_scripts/driver_grid.orion.sh          | 4 ++--
 driver_scripts/driver_grid.wcoss2.sh         | 4 ++--
 reg_tests/grid_gen/c96.uniform.sh            | 2 ++
 sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 | 4 ++--
 ush/fv3gfs_driver_grid.sh                    | 4 ++--
 7 files changed, 14 insertions(+), 12 deletions(-)

diff --git a/driver_scripts/driver_grid.hera.sh b/driver_scripts/driver_grid.hera.sh
index e78db425b..6975857a3 100755
--- a/driver_scripts/driver_grid.hera.sh
+++ b/driver_scripts/driver_grid.hera.sh
@@ -113,8 +113,8 @@ export lake_data_srce=MODISP_GLDBV3
 if [ $gtype = uniform ]; then
   export res=96
   export add_lake=true         # Add lake frac and depth to orography data.
-  export lake_cutoff=0.50      # ignore lake_frac < lake_cutoff when add_lake=T
-  export binary_lake=1         # return '1 > lake_frac >= lake_cutoff' as 1/0 when binary_lake=1/0
+  export lake_cutoff=0.50      # return 0 if lake_frac <  lake_cutoff & add_lake=T
+  export binary_lake=1         # return 1 if lake_frac >= lake_cutoff & add_lake=T
 elif [ $gtype = stretch ]; then
   export res=96
   export stretch_fac=1.5       # Stretching factor for the grid
diff --git a/driver_scripts/driver_grid.jet.sh b/driver_scripts/driver_grid.jet.sh
index c8ebc51f4..8cfce4d89 100755
--- a/driver_scripts/driver_grid.jet.sh
+++ b/driver_scripts/driver_grid.jet.sh
@@ -111,8 +111,8 @@ export lake_data_srce=MODISP_GLDBV3
 if [ $gtype = uniform ]; then
   export res=96
   export add_lake=true         # Add lake frac and depth to orography data.
-  export lake_cutoff=0.50      # ignore lake_frac < lake_cutoff when add_lake=T
-  export binary_lake=1         # return '1 > lake_frac >= lake_cutoff' as 1/0 when binary_lake=1/0
+  export lake_cutoff=0.50      # return 0 if lake_frac <  lake_cutoff & add_lake=T
+  export binary_lake=1         # return 1 if lake_frac >= lake_cutoff & add_lake=T
 elif [ $gtype = stretch ]; then
   export res=96
   export stretch_fac=1.5       # Stretching factor for the grid
diff --git a/driver_scripts/driver_grid.orion.sh b/driver_scripts/driver_grid.orion.sh
index bd12e944e..458abc318 100755
--- a/driver_scripts/driver_grid.orion.sh
+++ b/driver_scripts/driver_grid.orion.sh
@@ -110,8 +110,8 @@ export lake_data_srce=MODISP_GLDBV3
 if [ $gtype = uniform ]; then
   export res=96
   export add_lake=true         # Add lake frac and depth to orography data.
-  export lake_cutoff=0.50      # ignore lake_frac < lake_cutoff when add_lake=T
-  export binary_lake=1         # return '1 > lake_frac >= lake_cutoff' as 1/0 when binary_lake=1/0
+  export lake_cutoff=0.50      # return 0 if lake_frac <  lake_cutoff & add_lake=T
+  export binary_lake=1         # return 1 if lake_frac >= lake_cutoff & add_lake=T
 elif [ $gtype = stretch ]; then
   export res=96
   export stretch_fac=1.5       # Stretching factor for the grid
diff --git a/driver_scripts/driver_grid.wcoss2.sh b/driver_scripts/driver_grid.wcoss2.sh
index 4d55efef5..c5bb48919 100755
--- a/driver_scripts/driver_grid.wcoss2.sh
+++ b/driver_scripts/driver_grid.wcoss2.sh
@@ -109,8 +109,8 @@ export lake_data_srce=MODISP_GLDBV3
 if [ $gtype = uniform ]; then
   export res=96
   export add_lake=true         # Add lake frac and depth to orography data.
-  export lake_cutoff=0.50      # ignore lake_frac < lake_cutoff when add_lake=T
-  export binary_lake=1         # return '1 > lake_frac >= lake_cutoff' as 1/0 when binary_lake=1/0
+  export lake_cutoff=0.50      # return 0 if lake_frac <  lake_cutoff & add_lake=T
+  export binary_lake=1         # return 1 if lake_frac >= lake_cutoff & add_lake=T
 elif [ $gtype = stretch ]; then
   export res=96
   export stretch_fac=1.5       # Stretching factor for the grid
diff --git a/reg_tests/grid_gen/c96.uniform.sh b/reg_tests/grid_gen/c96.uniform.sh
index d4f569d51..09e18c843 100755
--- a/reg_tests/grid_gen/c96.uniform.sh
+++ b/reg_tests/grid_gen/c96.uniform.sh
@@ -15,6 +15,8 @@ export res=96
 export gtype=uniform
 export add_lake=true
 export lake_data_srce=MODISP_GLDBV3
+export lake_cutoff=0.50
+export binary_lake=1
 
 NCCMP=${NCCMP:-$(which nccmp)}
 
diff --git a/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 b/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
index 7816e21ad..92ae4269c 100644
--- a/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
+++ b/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
@@ -65,8 +65,8 @@ PROGRAM lake_frac
     CALL getarg(5, lakedepth_srce)
 
     IF (iargc() == 5) THEN
-      lake_cutoff = 0.20
-      binary_lake = 0
+      lake_cutoff = 0.50
+      binary_lake = 1
     ELSE
       CALL getarg(6, arg)
       READ(arg,*,iostat=stat) lake_cutoff
diff --git a/ush/fv3gfs_driver_grid.sh b/ush/fv3gfs_driver_grid.sh
index c0d456ac7..81ac33551 100755
--- a/ush/fv3gfs_driver_grid.sh
+++ b/ush/fv3gfs_driver_grid.sh
@@ -48,8 +48,8 @@ export gtype=${gtype:-uniform}  # grid type: uniform, stretch, nest, regional_gf
                                 # or regional_esg
 
 export add_lake=${add_lake:-false}      # add lake fraction and depth. uniform only.
-export lake_cutoff=${lake_cutoff:-0.20} # ignore lake_frac < lake_cutoff when add_lake=T
-export binary_lake=${binary_lake:0}     # return '1 > lake_frac >= lake_cutoff' as 1/0 when binary_lake=1/0
+export lake_cutoff=${lake_cutoff:-0.50} # return 0 if lake_frac <  lake_cutoff & add_lake=T
+export binary_lake=${binary_lake:-1}    # return 1 if lake_frac >= lake_cutoff & add_lake=T
 
 export make_gsl_orog=${make_gsl_orog:-false} # when true, create GSL drag suite orog files.
 

From 58e25cd2d20f7ecd7e2e5ce37630e1daf5339ae2 Mon Sep 17 00:00:00 2001
From: George Gayno <George.Gayno@noaa.gov>
Date: Fri, 3 Nov 2023 12:38:12 +0000
Subject: [PATCH 15/15] Minor update so code compiles with GNU.

Fixes #854.
---
 sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 b/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
index 92ae4269c..ad0f6e184 100644
--- a/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
+++ b/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90
@@ -663,9 +663,11 @@ SUBROUTINE write_lakedata_to_orodata(cs_res, cs_lakestat, cs_lakedpth)
       CALL nc_opchk(stat, "nf90_put_att: lake_frac:unit") 
       write(lakeinfo,'(a,f4.2,a,i1)') ' lake_frac cutoff=',lake_cutoff,'; binary_lake=',binary_lake
       IF (lakestatus_srce == "GLDBV3") THEN 
-        write(string,'(2a)') 'based on GLDBv3 (Choulga et al. 2019); missing lakes & added based on land_frac in this dataset;',trim(lakeinfo)
+        write(string,'(2a)') 'based on GLDBv3 (Choulga et al. 2019); missing lakes & added based on land_frac in this dataset;', &
+                              trim(lakeinfo)
       ELSE IF (lakestatus_srce == "GLDBV2") THEN 
-        write(string,'(2a)') 'based on GLDBv2 (Choulga et al. 2014); missing lakes & added based on land_frac in this dataset;',trim(lakeinfo)
+        write(string,'(2a)') 'based on GLDBv2 (Choulga et al. 2014); missing lakes & added based on land_frac in this dataset;', &
+                              trim(lakeinfo)
       ELSE IF (lakestatus_srce == "MODISP") THEN
         write(string,'(2a)') 'based on MODIS (2011-2015) product updated with two &
         Landsat products: the JRC water product (2016-2020) and the GLC-FCS30 (2020); &