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