fixing broken lake tiles in offline #335 and addresses #347 (#348)

# CABLE

## Description

This PR aims to achieve two tasks
- removing the hard-wired lake index from the offline code (including GW
module)
- maintaining lakes at saturation in offline model (soil-snow only)

Fixes #341
Addresses #347

## Type of change

- [x] Bug fix
- [x] Science advancement 

## Checklist

- [x] The new content is accessible and located in the appropriate
section.
- [ ] I have checked that links are valid and point to the intended
content.
- [ ] I have checked my code/text and corrected any misspellings

Please add a reviewer when ready for review.


<!-- readthedocs-preview cable start -->
----
📚 Documentation preview 📚:
https://cable--348.org.readthedocs.build/en/348/

<!-- readthedocs-preview cable end -->

src/offline/cbl_model_driver_offline.F90

@@ -53,6 +53,7 @@ SUBROUTINE cbm( ktau,dels, air, bgc, canopy, met,
 USE cable_phys_constants_mod, ONLY : CEMLEAF => EMLEAF
 USE cable_phys_constants_mod, ONLY : CEMSOIL => EMSOIL
 USE cable_phys_constants_mod, ONLY : CSBOLTZ => SBOLTZ
+USE cable_phys_constants_mod, ONLY : density_liq
     !mrd561
     USE cable_gw_hydro_module, ONLY : sli_hydrology,&
          soil_snow_gw
@@ -105,6 +106,20 @@ SUBROUTINE cbm( ktau,dels, air, bgc, canopy, met,
 !iFor testing
 cable_user%soil_struc="default"
 
+!At start of each time step ensure that lakes surface soil layer is at/above field capacity.
+!Diagnose any water needed to maintain this - this will be removed from 
+!runoff, drainage and/or deepest soil layer in surfbv
+!For offline case retain the water imbalance between timesteps - permits
+!balance to be maintained in the longer term. This differs to the coupled model
+!where %wb_lake is zero'd each time step (and river outflow is rescaled)
+WHERE( veg%iveg == lakes_cable .AND. ssnow%wb(:,1) < soil%sfc ) 
+  ssnow%wbtot1(:)  = REAL( ssnow%wb(:,1) ) * density_liq * soil%zse(1)
+  ssnow%wb(:,1) = soil%sfc
+  ssnow%wbtot2  = REAL( ssnow%wb(:,1) ) * density_liq * soil%zse(1)
+ENDWHERE
+ssnow%wb_lake = ssnow%wb_lake + MAX( ssnow%wbtot2 - ssnow%wbtot1, 0.)
+
+
 CALL ruff_resist( veg, rough, ssnow, canopy, veg%vlai, veg%hc, canopy%vlaiw )
 
 !jhan: this call to define air may be redundant

src/science/canopy/cbl_SurfaceWetness.F90

@@ -11,6 +11,7 @@ SUBROUTINE Surf_wetness_fact( cansat, canopy, ssnow,veg, met, soil, dels )
 
     USE cable_common_module
     USE cable_def_types_mod
+    USE grid_constants_mod_cbl, ONLY : lakes_cable
 ! physical constants
 USE cable_phys_constants_mod, ONLY : CTFRZ   => TFRZ
     !H!USE cable_gw_hydro_module, ONLY : calc_srf_wet_fraction
@@ -75,11 +76,11 @@ SUBROUTINE Surf_wetness_fact( cansat, canopy, ssnow,veg, met, soil, dels )
 
           IF( ssnow%snowd(i) > 0.1) ssnow%wetfac(i) = 0.9
 
-          IF ( veg%iveg(i) == 16 .and. met%tk(i) >= Ctfrz + 5. )   &
-               ssnow%wetfac(i) = 1.0 ! lakes: hard-wired number to be removed
+          IF ( veg%iveg(i) == lakes_cable .and. met%tk(i) >= Ctfrz + 5. )   &
+               ssnow%wetfac(i) = 1.0
 
-          IF( veg%iveg(i) == 16 .and. met%tk(i) < Ctfrz + 5. )   &
-               ssnow%wetfac(i) = 0.7 ! lakes: hard-wired number to be removed
+          IF( veg%iveg(i) == lakes_cable .and. met%tk(i) < Ctfrz + 5. )   &
+               ssnow%wetfac(i) = 0.7
 
        ENDDO
        ! owetfac introduced to reduce sharp changes in dry regions,

src/science/gw_hydro/cable_gw_hydro.F90

@@ -326,6 +326,7 @@ END SUBROUTINE remove_transGW
   !-------------------------------------------------------------------------
   SUBROUTINE ovrlndflx (dels, ssnow, soil,veg, canopy,sli_call )
     USE cable_common_module, ONLY : gw_params,cable_user
+    USE grid_constants_mod_cbl, ONLY : lakes_cable
 
     IMPLICIT NONE
     REAL, INTENT(IN)                         :: dels ! integration time step (s)
@@ -394,7 +395,7 @@ SUBROUTINE ovrlndflx (dels, ssnow, soil,veg, canopy,sli_call )
 
     !add back to the lakes to keep saturated instead of drying
     DO i=1,mp
-       IF (veg%iveg(i) .EQ. 16) THEN
+       IF (veg%iveg(i) .EQ. lakes_cable) THEN
           ssnow%fwtop(i) = ssnow%fwtop(i) + ssnow%rnof1(i)
           ssnow%rnof1(i) = 0._r_2
        END IF
@@ -1273,6 +1274,8 @@ END SUBROUTINE calc_equilibrium_water_content
 
   SUBROUTINE calc_srf_wet_fraction(ssnow,soil,met,veg)
 
+    USE grid_constants_mod_cbl, ONLY : lakes_cable
+
     IMPLICIT NONE
     TYPE(soil_snow_type), INTENT(INOUT)      :: ssnow  ! soil+snow variables
     TYPE(soil_parameter_type), INTENT(IN)    :: soil ! soil parameters
@@ -1297,11 +1300,11 @@ SUBROUTINE calc_srf_wet_fraction(ssnow,soil,met,veg)
        DO i=1,mp
           IF( ssnow%snowd(i) > 0.1) ssnow%wetfac(i) = 0.9
 
-          IF ( veg%iveg(i) == 16 .AND. met%tk(i) >= CTFRZ + 5. )   &
-               ssnow%wetfac(i) = 1.0 ! lakes: hard-wired number to be removed
+          IF ( veg%iveg(i) == lakes_cable .AND. met%tk(i) >= CTFRZ + 5. )   &
+               ssnow%wetfac(i) = 1.0 
 
-          IF( veg%iveg(i) == 16 .AND. met%tk(i) < CTFRZ + 5. )   &
-               ssnow%wetfac(i) = 0.7 ! lakes: hard-wired number to be removed
+          IF( veg%iveg(i) == lakes_cable .AND. met%tk(i) < CTFRZ + 5. )   &
+               ssnow%wetfac(i) = 0.7
        END DO
 
     ELSEIF (cable_user%gw_model) THEN
@@ -1366,11 +1369,11 @@ SUBROUTINE calc_srf_wet_fraction(ssnow,soil,met,veg)
 
           IF( ssnow%snowd(i) > 0.1) ssnow%wetfac(i) = 0.9
 
-          IF ( veg%iveg(i) == 16 .AND. met%tk(i) >= Ctfrz + 5. )   &
-               ssnow%wetfac(i) = 1.0 ! lakes: hard-wired number to be removed
+          IF ( veg%iveg(i) == lakes_cable .AND. met%tk(i) >= Ctfrz + 5. )   &
+               ssnow%wetfac(i) = 1.0 
 
-          IF( veg%iveg(i) == 16 .AND. met%tk(i) < Ctfrz + 5. )   &
-               ssnow%wetfac(i) = 0.7 ! lakes: hard-wired number to be removed
+          IF( veg%iveg(i) == lakes_cable .AND. met%tk(i) < Ctfrz + 5. )   &
+               ssnow%wetfac(i) = 0.7
 
        ENDDO
        ! owetfac introduced to reduce sharp changes in dry regions,
@@ -1713,6 +1716,7 @@ END SUBROUTINE aquifer_recharge
 
   SUBROUTINE subsurface_drainage(ssnow,soil,veg,dzmm)
     USE cable_common_module
+    USE grid_constants_mod_cbl, ONLY : lakes_cable
 
     IMPLICIT NONE
 
@@ -1841,7 +1845,7 @@ SUBROUTINE subsurface_drainage(ssnow,soil,veg,dzmm)
 
        !Keep "lakes" saturated forcing qhz = 0.  runoff only from lakes
        !overflowing
-       IF (soil%isoilm(i) .EQ. 9 .OR. veg%iveg(i) .GE. 16) THEN
+       IF (soil%isoilm(i) .EQ. 9 .OR. veg%iveg(i) .GE. lakes_cable) THEN
           ssnow%qhz(i) = 0._r_2
           ssnow%qhlev(i,:) = 0._r_2
        END IF

src/science/soilsnow/cbl_surfbv.F90

@@ -10,6 +10,8 @@ SUBROUTINE surfbv (dels, met, ssnow, soil, veg, canopy )
 
 USE smoisturev_mod,               ONLY: smoisturev
     USE cable_common_module
+    USE grid_constants_mod_cbl, ONLY : lakes_cable
+
 IMPLICIT NONE
 
     REAL, INTENT(IN) :: dels ! integration time step (s)
@@ -102,7 +104,7 @@ SUBROUTINE surfbv (dels, met, ssnow, soil, veg, canopy )
 
     !  Rescale drainage to remove water added to lakes (wb_lake)
     ssnow%sinfil = 0.0
-    WHERE( veg%iveg == 16 )
+    WHERE( veg%iveg == lakes_cable )
        ssnow%sinfil  = MIN( ssnow%rnof1, ssnow%wb_lake ) ! water that can be extracted from the rnof1
        ssnow%rnof1   = MAX( 0.0, ssnow%rnof1 - ssnow%sinfil )
        ssnow%wb_lake = MAX( 0.0, ssnow%wb_lake - ssnow%sinfil)