MODULE lake_mod !----------------------------------------------------------------------- ! NAME ! lake_mod ! ! DESCRIPTION ! Physics of the perennial lakes, i.e. the liquid water fraction ! 'peren_h2o_liq_frac' of the continental cells. ! ! AUTHORS & DATE ! JB Clement, 10/2026 ! ! NOTES ! The lakes are infinite reservoirs whose exchanges with the rest of ! the model are recorded by 'q_peren_h2o_lake'. They have no surface ! temperature of their own: the energy balance is the cell-averaged one. !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none ! PARAMETERS ! ---------- real, parameter :: newlake_frac = 1.e-3 ! Fraction of a new lake formed by the PEM ('lake_min_frac') [0-1] private :: compute_lake_energy_balance contains !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ pure elemental real function dry_ground_fraction(liq_frac,ice_frac) real, intent(in) :: liq_frac, ice_frac ! Lakes and glaciers are disjoint unless their summed areas exceed the cell. dry_ground_fraction = 1. - min(liq_frac + ice_frac,1.) end function dry_ground_fraction pure elemental real function exposed_lake_fraction(liq_frac,ice_frac) real, intent(in) :: liq_frac, ice_frac exposed_lake_fraction = max(0.,min(liq_frac + ice_frac,1.) - ice_frac) end function exposed_lake_fraction !======================================================================= SUBROUTINE compute_lake_weighted_TI(ngrid,nsoil,is_ocean) !----------------------------------------------------------------------- ! NAME ! compute_lake_weighted_TI ! ! DESCRIPTION ! Compute the thermal inertia of the cell weighted by the lake fraction. ! ! AUTHORS & DATE ! C. Metz, 06/2026 ! ! NOTES ! Oceanic grid cells are handled by the oceanic branch of the physics, ! so their thermal inertia must stay the oceanic one even if they carry ! a lake fraction. That can happen when the PEM turns an oceanic cell ! into a lake. !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use comsoil_h, only: inertiedat, lake_weighted_TI use surfdat_h, only: peren_h2o_liq_frac use callkeys_mod, only: lake_TI ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- integer, intent(in) :: ngrid integer, intent(in) :: nsoil logical, dimension(ngrid), intent(in) :: is_ocean ! True = ocean grid cell ! LOCAL VARIABLES ! --------------- integer :: isoil real, dimension(ngrid) :: lake_frac ! Lake fraction, set to 0 on oceanic cells ! CODE ! ---- lake_frac(:) = merge(0.,peren_h2o_liq_frac(:),is_ocean(:)) do isoil = 1,nsoil lake_weighted_TI(:,isoil) = lake_frac(:)*lake_TI + (1. - lake_frac(:))*inertiedat(:,isoil) end do END SUBROUTINE compute_lake_weighted_TI !======================================================================= !======================================================================= SUBROUTINE compute_lake_ice(ptimestep,pcapcal,peren_h2o_liq_frac,ztsurf,psnow,dq_peren_h2o_lake, & q_peren_h2o_lake,pdtsurf_hyd,lake_ice,dlake_ice,icecover_sublim4pem) !----------------------------------------------------------------------- ! NAME ! compute_lake_ice ! ! DESCRIPTION ! Treats ice for partial (perennial) liquid water area coverage. ! We solve the lake ice growth/melting with the implicit scheme, ! assuming a uniform surface temperature and growth can only occur ! due to conduction flux. The equations are: ! 1. rho_ice*Lf*dhice/dt = (T_freezing - T_surface_final)/R, where ! R = k_ice/h_ice + k_snow/h_snow; ! 2. T_surface_final = T_surface_provisional ! + rho_ice*Lf*dhice/dt*A_liq*delta_t/capcal. ! For freezing, solve (1) for T_surface_final and plug into (2), then ! solve for dhice/dt. ! For melting, set T_surface_final = T_freeze and solve (2) for dhice/dt. ! ! AUTHORS & DATE ! C. Metz, 06/2026 ! ! NOTES ! The lake ice state is kept separate from slab ocean sea_ice. ! We assume the lakes are not salty (no freezing point depression), ! to keep consistent with snowmelt above. ! Lakes cover 100% of 'peren_h2o_liq_frac'. 'lake_ice' [kg/m2] is only ! over the icy part, not grid-cell averaged, whereas 'q_peren_h2o_lake' ! and 'dq_peren_h2o_lake' are per m2 of cell. !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use watercommon_h, only: T_h2o_ice_liq, RLFTT, rhowaterice use ocean_slab_h, only: sno_cond, sno_den, ice_cond, ice_den ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real, intent(in) :: ptimestep ! Physical timestep [s] real, intent(in) :: pcapcal ! Surface heat capacity [J/(m^2 * K)] real, intent(in) :: peren_h2o_liq_frac ! Liquid water surface area coverage fraction [unitless] real, intent(in) :: ztsurf ! Tsurf you would get in the absence of hydrology (from other forcings) [K] real, intent(in) :: psnow ! Snow on top of the lake [kg/m^2] real, intent(in) :: dq_peren_h2o_lake ! Variation rate of the perennial lake reservoir [kg/m^2/s] real, intent(inout) :: q_peren_h2o_lake ! Tracker for the perennial liquid [kg/m^2] real, intent(inout) :: pdtsurf_hyd ! d_tsurf_dt [K/s] from hydrology real, intent(inout) :: lake_ice ! Lake ice thickness, per m2 of lake [kg/m^2] real, intent(inout) :: dlake_ice ! Tendency of lake ice, per m2 of lake [kg/m2/s] real, intent(out) :: icecover_sublim4pem ! Sublimation rate of the bare lake ice sent to the PEM, per m2 of lake [kg/m2/s] ! PARAMETERS ! ---------- real, parameter :: ice_thermal_conductivity = ice_cond/ice_den ! [W/(m*K)] from slab ocean model ! LOCAL VARIABLES ! --------------- real :: snow_thermal_conductivity ! [W/(m*K)] value taken from slab ocean model real :: hice ! Lake ice thickness [m] real :: tsurf_prov ! Tsurf you would get in the absence of lake ice real :: freezing_temp ! Lake freezing point real :: latent_energy_for_phase_change real :: d_hice ! Fractional change in ice thickness during the timestep [m] real :: hice_prov ! Provisional ice thickness after top growth/melt [m] real :: hsnow ! Snow thickness [m] real :: r_total ! Combined thermal resistance of ice + snow real :: leftover_energy real :: lake_ice_init real :: evap_amount ! Net water leaving the lake [kg/m2] real :: evap_icecover ! Sublimation of the ice cover, per m2 of lake [kg/m2] ! CODE ! ---- freezing_temp = T_h2o_ice_liq ! Assume fresh water snow_thermal_conductivity = sno_cond/sno_den lake_ice_init = lake_ice icecover_sublim4pem = 0. ! Turbdiff has already partitioned its single atmospheric flux between the ! snowpack and the lake. Exchange the lake share with the perennial inventory ! before applying lake-ice thermodynamics. Frost is taken care of on hydrol's ! existing snow pathway, while dew joins the lake. ! Whatever leaves the lake for the atmosphere is debited from the lake inventory, be it evaporated from the ! open water or sublimated from the lake ice: 'q_peren_h2o_lake' records every exchange between the lake and ! the rest of the model, which is what the water budget needs to see. ! 'lake_ice' is only the state of the frozen lid, thinned here by the part that sublimates from it. It is ! not an inventory of its own (freezing and melting just move water inside the lake), so it stays out of ! the water budget. q_peren_h2o_lake = q_peren_h2o_lake + dq_peren_h2o_lake*ptimestep evap_amount = -dq_peren_h2o_lake*ptimestep if (evap_amount > 0. .and. lake_ice > 0.) then ! The lid sublimates, entirely if there is not enough of it evap_icecover = evap_amount/peren_h2o_liq_frac icecover_sublim4pem = min(evap_icecover,lake_ice)/ptimestep lake_ice = max(lake_ice - evap_icecover,0.) end if ! Now handle thermodynamic ice melt/freeze based on tsurf tendencies tsurf_prov = ztsurf + pdtsurf_hyd*ptimestep hsnow = psnow/sno_den ! We assume the same amount of snow that is on the dry ground would also sit atop the lake ice hice = lake_ice/rhowaterice ! Convert from kg/m2 to m if (hice <= 0.) then ! Case A: no ice exists yet (open water cooling and initial freezing) if (tsurf_prov < freezing_temp) then latent_energy_for_phase_change = (freezing_temp - tsurf_prov)*pcapcal hice = latent_energy_for_phase_change/(rhowaterice*RLFTT*peren_h2o_liq_frac) pdtsurf_hyd = pdtsurf_hyd + (freezing_temp - tsurf_prov)/ptimestep end if else ! Case B: ice already exists (slab ice thermodynamics and complete melt check) if (tsurf_prov > freezing_temp) then ! Melt ice from the top if surface temperature exceeds freezing d_hice = pcapcal*(freezing_temp - tsurf_prov)/(rhowaterice*RLFTT*peren_h2o_liq_frac) hice_prov = hice + d_hice if (hice_prov < 0.) then ! Case b1: we melt ice, but we ran out of ice leftover_energy = (tsurf_prov - freezing_temp)*pcapcal - (hice*rhowaterice*RLFTT*peren_h2o_liq_frac) hice = 0. pdtsurf_hyd = pdtsurf_hyd + (freezing_temp + (leftover_energy/pcapcal) - tsurf_prov)/ptimestep else ! Case b2: we melt ice pdtsurf_hyd = pdtsurf_hyd + (freezing_temp - tsurf_prov)/ptimestep hice = hice_prov end if else ! Case b3: freeze ice, from the bottom r_total = (hice/ice_thermal_conductivity) + (hsnow/snow_thermal_conductivity) d_hice = ptimestep*(freezing_temp - tsurf_prov)/((RLFTT*rhowaterice)*((peren_h2o_liq_frac*ptimestep/pcapcal) + r_total)) pdtsurf_hyd = pdtsurf_hyd + (freezing_temp - tsurf_prov)/(ptimestep + (pcapcal*r_total)/peren_h2o_liq_frac) hice = hice + d_hice end if end if lake_ice = hice*rhowaterice ! Convert from m to kg/m2 dlake_ice = (lake_ice - lake_ice_init)/ptimestep END SUBROUTINE compute_lake_ice !======================================================================= !======================================================================= SUBROUTINE correct_lake_evap(rhoch,temp_a,qa,qsat_s,ps,sw_dn,lw_dn,emis,frac,beta_surround,lake_ice,snow,z1,z0,area,beta_lake,temp_lake,f_ibl) !----------------------------------------------------------------------- ! NAME ! correct_lake_evap ! ! DESCRIPTION ! Compute the correction of the lake evaporation based on: ! 1. an equilibrium lake surface temperature obtained from the lake ! energy balance (non-linearised Penman approach); ! 2. a reduction coefficient of the lake exchange by the internal ! boundary layer growing over the lake (Elliott 1958), if the ! surroundings are drier than the lake. ! ! AUTHORS & DATE ! JB Clement, 10/2026 ! ! NOTES ! Without this correction, the lake evaporation is overestimated in ! cells with a small lake fraction because (i) the lakes evaporate at ! the cell surface temperature, with no evaporative cooling of their own ! temperature, and (ii) under the cell-mean air humidity. ! ! The ratio of the lake evaporation to the one at the cell surface ! temperature is blended with the lake fraction so that the correction ! vanishes when a lake covers the cell. ! ! An ice-covered lake is treated the same way, with the albedo of the ! ice cover and the heat conducted from the lake water through the lake ! ice and the snow. !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use watercommon_h, only: Psat_water, T_h2o_ice_liq use callkeys_mod, only: alb_ocean, lake_ibl_max, albedosnow use ocean_slab_h, only: alb_ice_min, h_alb_ice, snow_min, ice_cond, sno_cond ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real, intent(in) :: rhoch ! rho*Ch*|U| [kg/m2/s] real, intent(in) :: temp_a ! Temperature of the first layer brought to the surface pressure [K] real, intent(in) :: qa ! Water vapour of the first layer [kg/kg] real, intent(in) :: qsat_s ! Saturation water vapour at the cell surface temperature [kg/kg] real, intent(in) :: ps ! Surface pressure [Pa] real, intent(in) :: sw_dn ! Incident SW flux at the surface [W/m2] real, intent(in) :: lw_dn ! Incident LW flux at the surface [W/m2] real, intent(in) :: emis ! Surface emissivity real, intent(in) :: frac ! Lake fraction of the cell [0-1] real, intent(in) :: beta_surround ! Evaporation efficiency of the lake surroundings [0-1] real, intent(in) :: lake_ice ! Lake ice, per m2 of lake [kg/m2] real, intent(in) :: snow ! Snow on top of the lake ice [kg/m2] real, intent(in) :: z1 ! Height of the first layer [m] real, intent(in) :: z0 ! Roughness length [m] real, intent(in) :: area ! Cell area [m2] real, intent(out) :: beta_lake ! Evaporation efficiency of the lakes [0-1] real, intent(out) :: temp_lake ! Equilibrium surface temperature of the lakes [K] real, intent(out) :: f_ibl ! Reduction coefficient of the lake exchange by the internal boundary layer [0-1] ! PARAMETERS ! ---------- integer, parameter :: nitermax = 50 ! Maximum number of iterations real, parameter :: dtemp_conv = 1.e-3 ! Convergence of the iterations [K] real, parameter :: temp_lake_min = 10. ! Lower bound of the lake temperature [K] real, parameter :: dtemp_bracket = 10. ! First step of the bracket search [K] ! LOCAL VARIABLES ! --------------- integer :: iter real :: fetch, ibl, rhoch_lake, rad_in real :: alb, alb_ice, frac_snow real :: cond_icecover ! Conductance of the ice cover [W/m2/K] real :: psat, qsat real :: phi, dphi, dtemp real :: temp_1, temp_2, step real :: evap_lake, evap_bulk, corr ! CODE ! ---- beta_lake = 1. temp_lake = temp_a f_ibl = 1. if (rhoch <= 0.) return ! No turbulent exchange with the atmosphere ! Internal boundary layer over the lake fetch = 0.5*sqrt(frac*area) ! We assume one square lake per cell and take the mean fetch ibl = min(lake_ibl_max,0.75*z0**0.2*fetch**0.8) ! Elliott 1958; Garratt 1990 ! The bulk formula assumes that the dry cell-mean air 'qa' lies at 'z1'. The air moistened by the lake ! fills the IBL, so this dry air is only found at the IBL top. ! If the IBL is thinner than the first layer, the air at 'z1' is still dry and the bulk formula holds. ! Otherwise, with logarithmic profiles of wind and humidity, the flux across the IBL relative to the bulk ! flux at 'z1' is weakened because the same humidity difference spans a thicker layer. if (ibl > z1) f_ibl = log(1. + z1/z0)/log(1. + ibl/z0) ! No IBL forms over a uniformly wet surface if the surroundings (snow/glaciers or saturated soil) evaporate as well f_ibl = 1. - (1. - beta_surround)*(1. - f_ibl) ! Heat and vapour are exchanged with the reduced coefficient rhoch_lake = f_ibl*rhoch ! Surface of the lake (open water or ice cover) frac_snow = max(0.,min(1.,snow/snow_min)) alb = alb_ocean + frac_snow*(albedosnow - alb_ocean) cond_icecover = 0. if (lake_ice > 0.) then ! Albedo law of the lake ice in hydrol alb_ice = albedosnow - (albedosnow - alb_ice_min)*exp(-lake_ice/h_alb_ice) alb = albedosnow*frac_snow + alb_ice*(1. - frac_snow) ! Heat conducted from the lake water through the lake ice and the snow cond_icecover = 1./(lake_ice/ice_cond + max(snow,0.)/sno_cond) end if rad_in = (1. - alb)*sw_dn + emis*lw_dn ! Radiation absorbed by the lake ! Bracket of the root of the lake energy balance, searched from 'temp_a' call compute_lake_energy_balance(temp_a,temp_a,qa,ps,rad_in,emis,rhoch_lake,cond_icecover,phi,dphi) temp_1 = temp_a temp_2 = temp_a step = dtemp_bracket if (phi > 0.) then ! The root is above 'temp_a': energy surplus, the lake is warmer than the air do iter = 1,nitermax temp_1 = temp_2 temp_2 = temp_2 + step step = 2.*step call compute_lake_energy_balance(temp_2,temp_a,qa,ps,rad_in,emis,rhoch_lake,cond_icecover,phi,dphi) if (phi <= 0.) exit end do else ! The root is below 'temp_a': energy deficit, the lake is colder than the air (e.g. evaporative cooling) do iter = 1,nitermax if (phi >= 0. .or. temp_1 <= temp_lake_min) exit temp_2 = temp_1 temp_1 = max(temp_1 - step,temp_lake_min) step = 2.*step call compute_lake_energy_balance(temp_1,temp_a,qa,ps,rad_in,emis,rhoch_lake,cond_icecover,phi,dphi) end do end if ! Equilibrium temperature of the lake surface ! Solved by Newton iterations from 'temp_a', with bisection when leaving the bracket ! The first iteration is the Penman equation (Penman 1948; Monteith & Unsworth 2013), the next ones remove its linearisation of 'qsat' do iter = 1,nitermax call compute_lake_energy_balance(temp_lake,temp_a,qa,ps,rad_in,emis,rhoch_lake,cond_icecover,phi,dphi) if (phi > 0.) then temp_1 = max(temp_1,temp_lake) else temp_2 = min(temp_2,temp_lake) end if dtemp = -phi/dphi if (temp_lake + dtemp < temp_1 .or. temp_lake + dtemp > temp_2) dtemp = 0.5*(temp_1 + temp_2) - temp_lake temp_lake = temp_lake + dtemp if (abs(dtemp) < dtemp_conv) exit end do ! The surface of open water cannot be colder than freezing temperature and the lake ice cannot be warmer if (lake_ice > 0.) then temp_lake = min(temp_lake,T_h2o_ice_liq) else temp_lake = max(temp_lake,T_h2o_ice_liq) end if ! Lake evaporation relative to the one at the cell surface temperature call Psat_water(temp_lake,ps,psat,qsat) evap_lake = rhoch_lake*(qsat - qa) ! Lake at its own temperature, with the reduced exchange evap_bulk = rhoch*(qsat_s - qa) ! Bulk formula of the model: lake at the cell surface temperature corr = 1. if (evap_bulk > 0.) corr = max(0.,min(evap_lake/evap_bulk,1.)) ! A lake covering the cell already sets the cell temperature and humidity, so the correction vanishes beta_lake = frac + (1. - frac)*corr END SUBROUTINE correct_lake_evap !======================================================================= !======================================================================= SUBROUTINE compute_lake_energy_balance(temp,temp_a,qa,ps,rad_in,emis,rhoch_lake,cond_icecover,phi,dphi) !----------------------------------------------------------------------- ! NAME ! compute_lake_energy_balance ! ! DESCRIPTION ! Energy balance of the lake surface without storage at temperature ! 'temp' and its derivative with respect to 'temp'. ! ! AUTHORS & DATE ! JB Clement, 10/2026 ! ! NOTES ! The heat storage of the lake is neglected since the lakes have no ! temperature of their own. The evaporation is nearly linear in the ! available energy and the stored heat is given back over the diurnal ! and annual cycles, so that it mainly shifts the evaporation in time ! and the time-integrated evaporation is preserved to first order. ! ! The energy balance is strictly decreasing with the temperature, so ! its root is unique. However, it is not concave when the saturation ! humidity is capped to 1 (air above the boiling point), hence the ! Newton iterations are kept within a bracket of the root by bisection ! to avoid them to overshoot to negative temperatures. !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use watercommon_h, only: RLVTT, RCPD, T_h2o_ice_liq, Psat_water, Lcpdqsat_water use radcommon_h, only: sigma use comcstfi_mod, only: cppd_ref ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real, intent(in) :: temp ! Lake surface temperature [K] real, intent(in) :: temp_a ! Temperature of the first layer brought to the surface pressure [K] real, intent(in) :: qa ! Water vapour of the first layer [kg/kg] real, intent(in) :: ps ! Surface pressure [Pa] real, intent(in) :: rad_in ! Radiative flux absorbed by the lake [W/m2] real, intent(in) :: emis ! Surface emissivity real, intent(in) :: rhoch_lake ! rho*Ch*|U| over the lake [kg/m2/s] real, intent(in) :: cond_icecover ! Conductance of the ice cover (0 for open water) [W/m2/K] real, intent(out) :: phi ! Energy balance [W/m2] real, intent(out) :: dphi ! Derivative of the energy balance [W/m2/K] ! LOCAL VARIABLES ! --------------- real :: psat, qsat, dqsat, dlnpsat ! CODE ! ---- call Psat_water(temp,ps,psat,qsat) call Lcpdqsat_water(temp,ps,psat,qsat,dqsat,dlnpsat) dqsat = dqsat*RCPD/RLVTT ! From L/cp*dqsat/dT to dqsat/dT ! Absorbed radiation - thermal emission - sensible heat - latent heat + conduction through the ice cover, without heat storage phi = rad_in - emis*sigma*temp**4 - cppd_ref*rhoch_lake*(temp - temp_a) - RLVTT*rhoch_lake*(qsat - qa) + cond_icecover*(T_h2o_ice_liq - temp) ! The latent heat term follows the saturation curve, whose slope vanishes above the boiling point ('qsat' capped to 1) dphi = -4.*emis*sigma*temp**3 - cppd_ref*rhoch_lake - RLVTT*rhoch_lake*dqsat - cond_icecover END SUBROUTINE compute_lake_energy_balance !======================================================================= END MODULE lake_mod