Lake information into the GEOS BCS workflow

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90

@@ -33,6 +33,7 @@ PROGRAM mkEASETilesParam
   use rmTinyCatchParaMod,    only : i_raster, j_raster 
   use rmTinyCatchParaMod,    only : RegridRasterReal     
   use rmTinyCatchParaMod,    only : Target_mean_land_elev
+  use rmTinyCatchParaMod,    only : LakeTopoCat_on_tiles_from_raster, AppendLakeVarsToTileNC4
   use process_hres_data,     only : histogram
   use LogRectRasterizeMod,   only : SRTM_maxcat, MAPL_UNDEF_R8   ! rasterize.F90
   use MAPL_SortMod
@@ -112,6 +113,11 @@ PROGRAM mkEASETilesParam
   character(len=512)     :: fname_mask
   character(len=400)     :: MAKE_BCS_INPUT_DIR
 
+  ! --- LakeTopoCat needed ---
+  real(REAL64), allocatable :: tile_lake_frac(:)
+  integer(1),   allocatable :: tile_is_lake_50(:)
+  integer :: n_tile
+
   ! --------------------------------------------------------------------------------------
 
   call get_environment_variable( "MAKE_BCS_INPUT_DIR", MAKE_BCS_INPUT_DIR )
@@ -976,7 +982,19 @@ PROGRAM mkEASETilesParam
 
   call MAPL_WriteTilingNC4('til/'//trim(gfile)//'.nc4', [EASELabel],[nc_ease],[nr_ease], &
        nc, nr, iTable, rTable)
-
+
+  ! ---  LakeTopoCat: compute lake fraction in tile space and append to tile nc4 ---   
+  n_tile = n_landlakelandice
+  allocate(tile_lake_frac(n_tile))
+  allocate(tile_is_lake_50(n_tile))
+  tile_lake_frac  = 0.0_REAL64
+  tile_is_lake_50 = 0_1
+
+  call LakeTopoCat_on_tiles_from_raster(n_tile, nc, nr, tileid_index, tile_lake_frac, tile_is_lake_50)
+  call AppendLakeVarsToTileNC4('til/'//trim(gfile)//'.nc4', tile_lake_frac, tile_is_lake_50, n_tile)
+
+  deallocate( tile_lake_frac, tile_is_lake_50)  
+  ! --- End LakeTopoCat ---   
   deallocate( tileid_index, catid_index,veg )
   deallocate( tile_area, ease_grid_area, tile_elev, my_land, all_id )
 

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rmTinyCatchParaMod.F90

@@ -55,6 +55,7 @@ module rmTinyCatchParaMod
   public REFORMAT_VEGFILES
   public Get_MidTime, Time_Interp_Fac
   public ascat_r0, jpl_canoph, NC_VarID, init_bcs_config  
+  public LakeTopoCat_on_tiles_from_raster, AppendLakeVarsToTileNC4
 
   ! The following variables define the details of the BCS version (data sources).
   ! Initialize to dummy values here and set to desired values in init_bcs_config().
@@ -902,6 +903,238 @@ END SUBROUTINE modis_alb_on_tiles
   ! subroutine pick_cat(sam,clr)
   !
   !----------------------------------------------------------------------
+  !----------------------------------------------------------------------  
+  !----------------------------------------------------------------------  
+  SUBROUTINE LakeTopoCat_on_tiles_from_raster(n_tile, nc_rst, nr_rst, tile_id, tile_lake_frac, tile_is_lake_50, rc)
+
+    ! Map HydroLAKES-TopoCat lake fraction (30 arc-sec, 10x10 deg tiles) to tile space
+    ! using the raster tile-id map directly.
+    !
+    ! Why NOT use rmap/create_mapping here?
+    ! ------------------------------------
+    ! 1) In mkCatchParam, the "tile creation" step (Step 01) happens BEFORE any
+    !    precomputed rmap exists (those are built later for other products).
+    ! 2) create_mapping() is explicitly designed for *land* tiles only (1..n_land),
+    !    which is perfect for MODIS snow albedo etc., but lakes need to be computed
+    !    for *all* tiles in the tile file (land + lake + ice + ocean).
+    ! 3) Lake inputs are already on the same global 30-arcsec grid as the raster
+    !    tile-id map (43200x21600), so we can map each 30" pixel directly to a tile
+    !    via tile_id(iG,jG) with unit weights (or count weights if desired).
+    !
+    ! Inputs:
+    !   n_tile   : number of tiles in the tile file (all types)
+    !   nc_rst,nr_rst : dims of tile_id raster grid (must match global 30" for this routine)
+    !   tile_id  : raster->tile mapping on 30-arcsec global grid; tile_id(iG,jG) in [1..n_tile]
+    !
+    ! Outputs:
+    !   tile_lake_frac  : mean lake fraction per tile (0..1)
+    !   tile_is_lake_50 : 1 if tile_lake_frac >= 0.5, else 0
+    !
+    ! Biljana Orescanin Feb 2026, SSAI@NASA
+
+    implicit none
+
+    integer,           intent(in)  :: n_tile
+    integer,           intent(in)  :: nc_rst, nr_rst
+    integer,           intent(in)  :: tile_id(1:nc_rst, 1:nr_rst)
+    real(REAL64),      intent(out) :: tile_lake_frac(1:n_tile)
+    integer(1),        intent(out) :: tile_is_lake_50(1:n_tile)
+    integer,           intent(out), optional :: rc
+
+    integer(kind=4), parameter :: nc_10 = 1200
+    integer(kind=4), parameter :: nr_10 = 1200
+
+    integer :: status, ncid, varid
+    integer :: ix, jx, ii, jj
+    integer :: iLL, jLL, iG, jG, tid
+    character*512 :: fname
+    character*2 :: HH, VV
+    character*512      :: MAKE_BCS_INPUT_DIR
+
+    real, allocatable :: lake_frac(:,:)                 ! 10x10 deg chunk, 30"
+    real(REAL64), allocatable :: sum_w(:), sum_fw(:)    ! counts + weighted fraction
+
+    call get_environment_variable("MAKE_BCS_INPUT_DIR", MAKE_BCS_INPUT_DIR)
+
+    ! ------------------------------------------------------------------
+    ! IMPORTANT:
+    ! LakeTopoCat inputs are defined on the native global 30-arcsec grid
+    ! (43200 x 21600). This routine assumes that the tile_id raster
+    ! (nc_rst x nr_rst) is also on that same 30" grid.
+    !
+    ! For workflows using GEOS5_10arcsec_mask, nx=43200 and ny=21600,
+    ! so lake fractions can be mapped directly via tile_id(iG,jG).
+    !
+    ! For coarser or alternative masks (e.g., 8640x4320), there is no
+    ! consistent 30" raster-to-tile alignment here. In those cases we
+    ! skip LakeTopoCat generation rather than attempting an implicit
+    ! remap that could introduce silent errors.
+    !
+    ! Lake variables are therefore OPTIONAL and only added when the
+    ! raster resolution matches the native 30" grid.
+    ! ------------------------------------------------------------------
+    if (nc_rst /= 43200 .or. nr_rst /= 21600) then
+       print *, 'NOTE: Skipping LakeTopoCat (requires 43200x21600). Got ', nc_rst, nr_rst
+       tile_lake_frac  = 0.0_REAL64
+       tile_is_lake_50 = 0_1
+       if (present(rc)) rc = 1
+       return
+    endif
+    allocate(lake_frac(1:nc_10, 1:nr_10))
+    allocate(sum_w (1:n_tile))
+    allocate(sum_fw(1:n_tile))
+
+    sum_w  = 0.0_8
+    sum_fw = 0.0_8
+
+    ! Loop through 36x18 tiles (HxxVyy)
+    do jx = 1, 18
+      do ix = 1, 36
+
+        write(HH,'(i2.2)') ix
+        write(VV,'(i2.2)') jx
+
+        fname = trim(MAKE_BCS_INPUT_DIR)//'/lake/lake_mask/v1/LakeTopoCat_30arcsec_H'//HH//'V'//VV//'.nc4'
+
+        status = NF_OPEN(trim(fname), NF_NOWRITE, ncid)
+        if (status /= NF_NOERR) then
+           print *, 'ERROR: Lake mask tile file not found: ', trim(fname)
+           print *, 'STOPPING.'
+           stop
+        endif
+
+        status = NF_INQ_VARID(ncid, 'lake_presence_frac', varid) ; VERIFY_(status)
+        status = NF_GET_VARA_REAL(ncid, varid, (/1,1/), (/nc_10,nr_10/), lake_frac) ; VERIFY_(status)
+        status = NF_CLOSE(ncid) ; VERIFY_(status)
+
+        ! Lower-left global indices of this 10x10deg chunk (1-based)
+        iLL = (ix-1)*nc_10 + 1
+        jLL = (jx-1)*nr_10 + 1
+
+        ! Loop through chunk pixels
+        do jj = 1, nr_10
+          do ii = 1, nc_10
+
+            iG = ii + iLL - 1
+            jG = jj + jLL - 1
+
+            tid = tile_id(iG, jG)
+            if (tid < 1 .or. tid > n_tile) cycle
+
+            ! unit weight: each 30" pixel counts once for its tile
+            sum_fw(tid) = sum_fw(tid) + real(lake_frac(ii,jj), REAL64)
+            sum_w (tid) = sum_w (tid) + 1.0_8
+
+          end do
+        end do
+
+      end do
+    end do
+
+    ! Finalize fraction and binary flag
+    tile_lake_frac = 0.0_8
+    where (sum_w > 0.0_8)
+      tile_lake_frac = sum_fw / sum_w
+    endwhere
+
+    tile_is_lake_50 = 0_1
+    do tid = 1, n_tile
+      if (tile_lake_frac(tid) >= 0.5_8) tile_is_lake_50(tid) = 1_1
+    end do
+
+    if (present(rc)) rc = 0
+
+    if (allocated(lake_frac)) deallocate(lake_frac)
+    if (allocated(sum_w))     deallocate(sum_w)
+    if (allocated(sum_fw))    deallocate(sum_fw)    
+
+  END SUBROUTINE LakeTopoCat_on_tiles_from_raster
+  !----------------------------------------------------------------------  
+  !----------------------------------------------------------------------
+  SUBROUTINE AppendLakeVarsToTileNC4(tilefile, tile_lake_frac, tile_is_lake_50, n_tile)
+
+    implicit none
+    character(*), intent(in) :: tilefile
+    integer,      intent(in) :: n_tile
+    real(REAL64), intent(in) :: tile_lake_frac(1:n_tile)
+    integer(1),   intent(in) :: tile_is_lake_50(1:n_tile)
+
+    integer :: status, ncid, ndims, nvars, ngatts, unlimdimid
+    integer :: dimid_tile, dimlen, d
+    integer :: varid_frac, varid_l50
+    character(len=NF_MAX_NAME) :: dname
+    character(len=64), parameter :: msg_l50 = 'Flag: tile_lake_frac >= 0.5'
+    character(len=64), parameter :: msg_frac = 'Lake fraction in tile'
+    integer, allocatable :: tmp_l50(:)
+    logical :: need_redef
+
+    status = NF_OPEN(trim(tilefile), NF_WRITE, ncid) ; VERIFY_(status)
+
+    ! ---- find tile dimension by matching length ----
+    status = NF_INQ(ncid, ndims, nvars, ngatts, unlimdimid) ; VERIFY_(status)
+
+    dimid_tile = -1
+    do d = 1, ndims
+      status = NF_INQ_DIM(ncid, d, dname, dimlen) ; VERIFY_(status)
+      if (dimlen == n_tile) then
+        dimid_tile = d
+        exit
+      endif
+    enddo
+
+    if (dimid_tile < 0) then
+      print *, 'ERROR: could not find tile dimension of length ', n_tile, ' in ', trim(tilefile)
+      stop
+    endif
+
+    ! ---- ensure tile_lake_frac exists ----
+    need_redef = .false.
+    status = NF_INQ_VARID(ncid, 'tile_lake_frac', varid_frac)
+    if (status /= NF_NOERR) need_redef = .true.
+
+    ! ---- ensure tile_is_lake_50pct exists ----
+    status = NF_INQ_VARID(ncid, 'tile_is_lake_50pct', varid_l50)
+    if (status /= NF_NOERR) need_redef = .true.
+
+    if (need_redef) then
+      status = NF_REDEF(ncid) ; VERIFY_(status)
+
+      status = NF_INQ_VARID(ncid, 'tile_lake_frac', varid_frac)
+      if (status /= NF_NOERR) then
+        status = NF_DEF_VAR(ncid, 'tile_lake_frac', NF_DOUBLE, 1, (/dimid_tile/), varid_frac)
+        VERIFY_(status)
+        ! optional metadata:
+        status = NF_PUT_ATT_TEXT(ncid, varid_frac, 'long_name', len_trim(msg_frac), msg_frac) ; VERIFY_(status)
+        status = NF_PUT_ATT_TEXT(ncid, varid_frac, 'units',      1, '1')                      ; VERIFY_(status)
+      endif
+
+      status = NF_INQ_VARID(ncid, 'tile_is_lake_50pct', varid_l50)
+      if (status /= NF_NOERR) then
+        ! Use NF_INT so we can write with NF_PUT_VARA_INT safely.
+        status = NF_DEF_VAR(ncid, 'tile_is_lake_50pct', NF_INT, 1, (/dimid_tile/), varid_l50)
+        VERIFY_(status)
+        ! optional metadata:
+        status = NF_PUT_ATT_TEXT(ncid, varid_l50, 'long_name', len_trim(msg_l50), msg_l50) ; VERIFY_(status)
+        status = NF_PUT_ATT_TEXT(ncid, varid_l50, 'units',      1, '1')                            ; VERIFY_(status)
+      endif
+
+      status = NF_ENDDEF(ncid) ; VERIFY_(status)
+    endif
+
+    ! ---- write data ----
+    status = NF_PUT_VARA_DOUBLE(ncid, varid_frac, (/1/), (/n_tile/), tile_lake_frac)
+    VERIFY_(status)
+
+    allocate(tmp_l50(n_tile), stat=status); VERIFY_(status)
+    tmp_l50 = int(tile_is_lake_50)   ! 0/1
+    status = NF_PUT_VARA_INT(ncid, varid_l50, (/1/), (/n_tile/), tmp_l50)
+    VERIFY_(status)
+    deallocate(tmp_l50)
+
+    status = NF_CLOSE(ncid) ; VERIFY_(status)
+
+  END SUBROUTINE AppendLakeVarsToTileNC4
+  !----------------------------------------------------------------------  
 
   SUBROUTINE supplemental_tile_attributes(nx,ny,regrid,dateline,fnameTil, Rst_id)
 
@@ -949,6 +1182,12 @@ SUBROUTINE supplemental_tile_attributes(nx,ny,regrid,dateline,fnameTil, Rst_id)
     character(len=128)                     :: gName(2)
     logical,           allocatable         :: IsOcean(:)
 
+    ! --- LakeTopoCat needed ---
+    real(REAL64), allocatable :: tile_lake_frac(:)
+    integer(1),   allocatable :: tile_is_lake_50(:)
+    integer :: n_tile
+    integer :: rc_lake
+
     ! -----------------------------------------------------
     !
     ! get elevation (q0) from "gtopo30" raster file ("srtm30_withKMS_2.5x2.5min.data")
@@ -993,6 +1232,13 @@ SUBROUTINE supplemental_tile_attributes(nx,ny,regrid,dateline,fnameTil, Rst_id)
     open (10,file=fname,status='old',action='read',form='formatted')
     read (10,*)ip
 
+    ! --- LakeTopoCat needed ---
+    n_tile = ip
+    allocate(tile_lake_frac(n_tile))
+    allocate(tile_is_lake_50(n_tile))
+    tile_lake_frac  = 0.0_REAL64
+    tile_is_lake_50 = 0    
+
     allocate(id(       ip))
     allocate(i_index(  ip))
     allocate(j_index(  ip))
@@ -1210,11 +1456,20 @@ SUBROUTINE supplemental_tile_attributes(nx,ny,regrid,dateline,fnameTil, Rst_id)
 
     fname=trim(fnameTil)//'.nc4'
     call MAPL_WriteTilingNC4(fname,  gName(1:n_grid), im(1:n_grid), jm(1:n_grid), nx, ny, iTable, rTable, N_PfafCat=SRTM_maxcat, rc=status)
-
+
+    ! --- LakeTopoCat: compute lake fraction in tile space and append to tile nc4 ---
+    ! LakeTopoCat tiles are on the global 30-arcsec grid (43200x21600).
+    ! Step 01 runs with nx=43200, ny=21600 for GEOS5_10arcsec_mask workflows.
+    call LakeTopoCat_on_tiles_from_raster(n_tile, nx, ny, Rst_id, tile_lake_frac, tile_is_lake_50)
+    if (rc_lake == 0) then
+      call AppendLakeVarsToTileNC4(fname, tile_lake_frac, tile_is_lake_50, n_tile)
+    endif    
+
     deallocate (rTable, iTable)
     deallocate (limits)
     deallocate (catid)
     deallocate (q0)
+    deallocate(tile_lake_frac, tile_is_lake_50)
     if(regrid) then
        deallocate(raster)
     endif