MODULE liquidbody !----------------------------------------------------------------------- ! NAME ! liquidbody ! ! DESCRIPTION ! State and helpers for Generic PEM fields that are only needed when ! handing control back to the PCM. ! ! AUTHORS & DATE ! C. Metz, 07/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use numerics, only: dp, di, k4, eps ! DECLARATION ! ----------- implicit none ! VARIABLES ! --------- real(dp), protected :: slab_depth_PCM ! Slab ocean depth of the PCM wind-mixed layer [m] real(dp), protected :: soil_h2o_capacity_PCM ! Soil water bucket capacity of the PCM (soil moisture) [kg/m2] ! PCM restart fields read at PEM startup and rewritten at PEM handoff. real(dp), dimension(:,:), allocatable :: tslab_PCM real(dp), dimension(:), allocatable :: tsea_ice_PCM real(dp), dimension(:), allocatable :: tice_PCM real(dp), dimension(:), allocatable :: q_peren_h2o_PCM real(dp), dimension(:), allocatable :: qsurf_h2o_liq_PCM ! PARAMETERS ! ---------- real(dp), parameter :: ocean_tfreeze = 271.35_dp ! Slab-ocean t_freeze in phygeneric/ocean_slab_mod.F90 [K] real(dp), parameter, private :: slab2_thickness = 150._dp ! Slab-ocean second-layer thickness [m] contains !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !======================================================================= SUBROUTINE ini_liquidbody() !----------------------------------------------------------------------- ! NAME ! ini_liquidbody ! ! DESCRIPTION ! Allocate PCM handoff fields. ! ! AUTHORS & DATE ! C. Metz, 07/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: ngrid ! DECLARATION ! ----------- implicit none ! CODE ! ---- allocate(tslab_PCM(ngrid,2)) allocate(tsea_ice_PCM(ngrid)) allocate(tice_PCM(ngrid)) allocate(q_peren_h2o_PCM(ngrid)) allocate(qsurf_h2o_liq_PCM(ngrid)) tslab_PCM(:,:) = 0._dp tsea_ice_PCM(:) = 0._dp tice_PCM(:) = 0._dp q_peren_h2o_PCM(:) = 0._dp qsurf_h2o_liq_PCM(:) = 0._dp END SUBROUTINE ini_liquidbody !======================================================================= !======================================================================= SUBROUTINE end_liquidbody() !----------------------------------------------------------------------- ! NAME ! end_liquidbody ! ! DESCRIPTION ! Deallocate PCM handoff fields. ! ! AUTHORS & DATE ! C. Metz, 07/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none ! CODE ! ---- if (allocated(tslab_PCM)) deallocate(tslab_PCM) if (allocated(tsea_ice_PCM)) deallocate(tsea_ice_PCM) if (allocated(tice_PCM)) deallocate(tice_PCM) if (allocated(q_peren_h2o_PCM)) deallocate(q_peren_h2o_PCM) if (allocated(qsurf_h2o_liq_PCM)) deallocate(qsurf_h2o_liq_PCM) END SUBROUTINE end_liquidbody !======================================================================= !======================================================================= SUBROUTINE set_slab_depth_PCM(slab_depth_in) !----------------------------------------------------------------------- ! NAME ! set_slab_depth_PCM ! ! DESCRIPTION ! Setter for 'slab_depth_PCM'. ! ! AUTHORS & DATE ! C. Metz, 07/2026 ! ! NOTES ! Stores the PCM slab-ocean upper-layer depth read from callphys.def. !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use utility, only: real2str use display, only: print_msg, LVL_NFO use stoppage, only: stop_clean ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), intent(in) :: slab_depth_in ! CODE ! ---- slab_depth_PCM = slab_depth_in call print_msg('slab_depth = '//real2str(slab_depth_PCM),LVL_NFO) if (slab_depth_PCM <= 0._dp) call stop_clean(__FILE__,__LINE__,'''slab_depth'' must be positive!',1) END SUBROUTINE set_slab_depth_PCM !======================================================================= !======================================================================= SUBROUTINE set_soil_h2o_capacity_PCM(soil_h2o_capacity_in) !----------------------------------------------------------------------- ! NAME ! set_soil_h2o_capacity_PCM ! ! DESCRIPTION ! Setter for 'soil_h2o_capacity_PCM'. ! ! AUTHORS & DATE ! JB Clement, 08/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use utility, only: real2str use display, only: print_msg, LVL_NFO use stoppage, only: stop_clean ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), intent(in) :: soil_h2o_capacity_in ! CODE ! ---- soil_h2o_capacity_PCM = soil_h2o_capacity_in call print_msg('soil_h2o_capacity = '//real2str(soil_h2o_capacity_PCM),LVL_NFO) if (soil_h2o_capacity_PCM <= 0._dp) call stop_clean(__FILE__,__LINE__,'''soil_h2o_capacity'' must be positive!',1) END SUBROUTINE set_soil_h2o_capacity_PCM !======================================================================= !======================================================================= SUBROUTINE build4PCM_liquidbody(tslab4PCM,tsea_ice4PCM,tice4PCM,q_peren_h2o4PCM,qsurf_h2o_liq4PCM) !----------------------------------------------------------------------- ! NAME ! build4PCM_liquidbody ! ! DESCRIPTION ! Build the slab ocean and perennial water fields for the PCM. ! ! AUTHORS & DATE ! JB Clement, 08/2026 ! ! NOTES ! These fields are not evolved by the PEM: the incoming PCM state is ! handed back unchanged. ! In particular, ocean temperature is evolved entirely by the PCM. ! Indeed, sea ice balances rapidly, mostly at seasonal scale or in the ! span of a few years. It is also difficult to justify a PEM evolution ! model for ocean temperature given the complexity of the process. !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use hydrology, only: do_hydrology use display, only: print_msg, LVL_NFO ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), dimension(:,:), intent(out) :: tslab4PCM real(dp), dimension(:), intent(out) :: tsea_ice4PCM, tice4PCM real(dp), dimension(:), intent(out) :: q_peren_h2o4PCM, qsurf_h2o_liq4PCM ! CODE ! ---- ! Default initialization tslab4PCM(:,:) = 0._dp tsea_ice4PCM(:) = 0._dp tice4PCM(:) = 0._dp q_peren_h2o4PCM(:) = 0._dp qsurf_h2o_liq4PCM(:) = 0._dp if (.not. do_hydrology) return ! Build slab ovean and perennial lakes for the PCM call print_msg('> Building slab ocean and perennial lakes for the PCM',LVL_NFO) tslab4PCM(:,:) = tslab_PCM(:,:) tsea_ice4PCM(:) = tsea_ice_PCM(:) tice4PCM(:) = tice_PCM(:) q_peren_h2o4PCM(:) = q_peren_h2o_PCM(:) qsurf_h2o_liq4PCM(:) = qsurf_h2o_liq_PCM(:) END SUBROUTINE build4PCM_liquidbody !======================================================================= !======================================================================= SUBROUTINE tranform_lake_ocean4PCM(peren_h2o_liq_frac4PCM,is_ocean4PCM, & lakeice4PCM,seaice4PCM,seaice_frac4PCM, & tsurf4PCM,tsoil4PCM,tslab4PCM, & tsea_ice4PCM,tice4PCM, & q_peren_h2o4PCM,qsurf_h2o_liq4PCM) !----------------------------------------------------------------------- ! NAME ! tranform_lake_ocean4PCM ! ! DESCRIPTION ! Normalize all PCM handoff fields in cells whose water coverage ! crossed the lake/ocean threshold during the PEM leg. ! ! AUTHORS & DATE ! C. Metz, 07/2026 ! ! NOTES ! We follow the PCM's definitions of tslab, tsurf, etc. !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: ngrid, nslope, nsoil_PCM, cell_area use slopes, only: iflat use hydrology, only: is_ocean_PCM, lake_volume, lake_maxdepth, h2o_liq_density use soil, only: mlayer use display, only: print_msg, LVL_NFO, LVL_WRN use utility, only: int2str, real2str ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), dimension(:), intent(in) :: peren_h2o_liq_frac4PCM logical(k4), dimension(:), intent(in) :: is_ocean4PCM real(dp), dimension(:,:), intent(inout) :: tsurf4PCM real(dp), dimension(:), intent(inout) :: lakeice4PCM, seaice4PCM, seaice_frac4PCM real(dp), dimension(:,:,:), intent(inout) :: tsoil4PCM real(dp), dimension(:,:), intent(inout) :: tslab4PCM real(dp), dimension(:), intent(inout) :: tsea_ice4PCM, tice4PCM real(dp), dimension(:), intent(inout) :: q_peren_h2o4PCM, qsurf_h2o_liq4PCM ! LOCAL VARIABLES ! --------------- integer(di) :: ig, islope, isoil integer(di) :: n_lake2ocean, n_ocean2lake logical(k4) :: is_lake2ocean, is_ocean2lake real(dp) :: land_frac, land_area, actual_mass, available_mass, qsoil_fill real(dp) :: qsurf_old real(dp) :: slab1_middepth, slab2_middepth real(dp) :: lakeice_max ! CODE ! ---- n_lake2ocean = 0 n_ocean2lake = 0 slab1_middepth = 0.5_dp*slab_depth_PCM slab2_middepth = slab_depth_PCM + 0.5_dp*slab2_thickness do ig = 1,ngrid is_lake2ocean = (.not. is_ocean_PCM(ig)) .and. is_ocean4PCM(ig) is_ocean2lake = is_ocean_PCM(ig) .and. (.not. is_ocean4PCM(ig)) if (is_lake2ocean) then n_lake2ocean = n_lake2ocean + 1 ! Change lake ice to sea ice if there is any, and prepare slab temperatures. seaice_frac4PCM(ig) = 0._dp seaice4PCM(ig) = 0._dp tslab4PCM(ig,:) = 0._dp tsea_ice4PCM(ig) = 0._dp tice4PCM(ig) = 0._dp if (lakeice4PCM(ig) > eps) then seaice_frac4PCM(ig) = 1._dp seaice4PCM(ig) = max(0._dp,lakeice4PCM(ig)*peren_h2o_liq_frac4PCM(ig)) tslab4PCM(ig,:) = ocean_tfreeze tsea_ice4PCM(ig) = tsurf4PCM(ig,iflat) tice4PCM(ig) = tsurf4PCM(ig,iflat) ! tsurf stays as the ice surface temperature, since we assume its 100% covered else tslab4PCM(ig,1) = interp_tsoil_to_tslab(ig,iflat,tsoil4PCM,tsurf4PCM,slab1_middepth) tslab4PCM(ig,2) = interp_tsoil_to_tslab(ig,iflat,tsoil4PCM,tsurf4PCM,slab2_middepth) end if lakeice4PCM(ig) = 0._dp ! Discard dry-ground soil liquid and commute water debts qsurf_old = qsurf_h2o_liq4PCM(ig) land_frac = max(0._dp,1._dp - peren_h2o_liq_frac4PCM(ig)) if (qsurf_old > 0._dp .and. land_frac > eps) then ! cmetz later when we implement explicit tracking of soil moisture in PEM, we should saturate the PEM soil, taking water ! from the lake. right now we just forget about soil moisture in the PEM/PCM, since we do not track it end if qsurf_h2o_liq4PCM(ig) = q_peren_h2o4PCM(ig) q_peren_h2o4PCM(ig) = 0._dp else if (is_ocean2lake) then n_ocean2lake = n_ocean2lake + 1 ! Change sea ice to lake ice if there is any. lakeice4PCM(ig) = 0._dp if (peren_h2o_liq_frac4PCM(ig) > eps .and. seaice4PCM(ig) > eps .and. seaice_frac4PCM(ig) > eps) then ! Sea ice is spread over the lake area, capped at the maximum depth of the lake lakeice4PCM(ig) = min(seaice4PCM(ig)*seaice_frac4PCM(ig)/peren_h2o_liq_frac4PCM(ig),h2o_liq_density*lake_maxdepth(ig)) end if seaice4PCM(ig) = 0._dp seaice_frac4PCM(ig) = 0._dp ! Rebuild soil temperatures for the newly non-oceanic cell. do islope = 1,nslope if (lakeice4PCM(ig) > eps) then tsurf4PCM(ig,islope) = tsea_ice4PCM(ig) else tsurf4PCM(ig,islope) = tslab4PCM(ig,1) end if do isoil = 1,nsoil_PCM tsoil4PCM(ig,isoil,islope) = interp_tslab_to_tsoil(tslab4PCM(ig,:),tsurf4PCM(ig,islope),mlayer(isoil - 1)) end do end do tsea_ice4PCM(ig) = 0._dp tice4PCM(ig) = 0._dp tslab4PCM(ig,:) = 0._dp ! Fill the newly exposed dry-ground reservoir and commute water debts q_peren_h2o4PCM(ig) = qsurf_h2o_liq4PCM(ig) land_frac = max(0._dp,1._dp - peren_h2o_liq_frac4PCM(ig)) land_area = land_frac*cell_area(ig) qsoil_fill = 0._dp if (land_area > eps) then ! cmetz we need to fix this later when we implement explicit tracking of soil moisture in PEM ! right now we just saturate the soil and take the water from nowhere, which is OK since the PCM/PEM worlds are separate qsoil_fill = soil_h2o_capacity_PCM end if qsurf_h2o_liq4PCM(ig) = qsoil_fill end if end do call print_msg('Lake-to-ocean transitions: '//int2str(n_lake2ocean),LVL_NFO) call print_msg('Ocean-to-lake transitions: '//int2str(n_ocean2lake),LVL_NFO) END SUBROUTINE tranform_lake_ocean4PCM !======================================================================= !======================================================================= SUBROUTINE mask_lake_ocean4PCM(is_ocean4PCM,peren_h2o_liq_frac4PCM,lakeice4PCM, & q_peren_h2o4PCM,seaice4PCM,seaice_frac4PCM, & tslab4PCM,tsea_ice4PCM,tice4PCM) !----------------------------------------------------------------------- ! NAME ! mask_lake_ocean4PCM ! ! DESCRIPTION ! Enforce the mutual exclusivity of lakes and slab ocean variables ! in the fields handed over to the PCM. ! ! AUTHORS & DATE ! JB Clement, 08/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- logical(k4), dimension(:), intent(in) :: is_ocean4PCM real(dp), dimension(:), intent(inout) :: peren_h2o_liq_frac4PCM, lakeice4PCM, q_peren_h2o4PCM real(dp), dimension(:), intent(inout) :: seaice4PCM, seaice_frac4PCM, tsea_ice4PCM, tice4PCM real(dp), dimension(:,:), intent(inout) :: tslab4PCM ! CODE ! ---- where (is_ocean4PCM(:)) ! Lakes do not exist on the cells given to the slab ocean peren_h2o_liq_frac4PCM(:) = 0._dp lakeice4PCM(:) = 0._dp q_peren_h2o4PCM(:) = 0._dp else where ! The slab ocean does not exist on the continental cells seaice4PCM(:) = 0._dp seaice_frac4PCM(:) = 0._dp tsea_ice4PCM(:) = 0._dp tice4PCM(:) = 0._dp tslab4PCM(:,1) = 0._dp tslab4PCM(:,2) = 0._dp end where END SUBROUTINE mask_lake_ocean4PCM !======================================================================= !======================================================================= FUNCTION interp_tsoil_to_tslab(ig,islope,tsoil,tsurf,z) RESULT(temp) !----------------------------------------------------------------------- ! NAME ! interp_tsoil_to_tslab ! ! DESCRIPTION ! Interpolate the outgoing PCM soil profile to a slab depth. ! ! AUTHORS & DATE ! C. Metz, 07/2026 ! ! NOTES ! We know tsoil and tsurf, and we are mapping onto 2 slab depths. !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: nsoil_PCM use soil, only: mlayer ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- integer(di), intent(in) :: ig, islope real(dp), dimension(:,:,:), intent(in) :: tsoil real(dp), dimension(:,:), intent(in) :: tsurf real(dp), intent(in) :: z ! RESULT ! ------ real(dp) :: temp ! Interpolated temperature [K] ! LOCAL VARIABLES ! --------------- integer(di) :: isoil real(dp) :: upper_depth, lower_depth real(dp) :: upper_temp, lower_temp ! CODE ! ---- if (z <= 0._dp) then temp = tsurf(ig,islope) else if (z <= mlayer(0)) then temp = tsurf(ig,islope) + (tsoil(ig,1,islope) - tsurf(ig,islope))*z/max(mlayer(0),eps) else if (z >= mlayer(nsoil_PCM - 1)) then temp = tsoil(ig,nsoil_PCM,islope) ! This is nearly always going to be the case so we just set tslab = lowest soil temp else temp = tsoil(ig,nsoil_PCM,islope) do isoil = 2,nsoil_PCM if (z <= mlayer(isoil - 1)) then upper_depth = mlayer(isoil - 2) lower_depth = mlayer(isoil - 1) upper_temp = tsoil(ig,isoil - 1,islope) lower_temp = tsoil(ig,isoil,islope) temp = upper_temp + (lower_temp - upper_temp)*(z - upper_depth)/max(lower_depth - upper_depth,eps) exit end if end do end if ! Clamp it in case the soil was super cold or super hot temp = min(373._dp, max(ocean_tfreeze, temp)) END FUNCTION interp_tsoil_to_tslab !======================================================================= !======================================================================= FUNCTION interp_tslab_to_tsoil(tslab,tsurf,z) RESULT(temp) !----------------------------------------------------------------------- ! NAME ! interp_tslab_to_tsoil ! ! DESCRIPTION ! Interpolate the incoming slab profile onto a PCM soil depth. ! ! AUTHORS & DATE ! C. Metz, 07/2026 ! ! NOTES ! We know 3 points (tsurf, tslab1, and tslab2) and we are mapping onto soil depths. !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), dimension(:), intent(in) :: tslab real(dp), intent(in) :: tsurf, z ! RESULT ! ------ real(dp) :: temp ! Interpolated temperature [K] ! LOCAL VARIABLES ! --------------- real(dp) :: slab1_middepth, slab2_middepth ! CODE ! ---- slab1_middepth = 0.5_dp*slab_depth_PCM slab2_middepth = slab_depth_PCM + 0.5_dp*slab2_thickness if (z <= 0._dp) then temp = tsurf else if (z <= slab1_middepth) then temp = tsurf + (tslab(1) - tsurf)*z/max(slab1_middepth,eps) else if (z >= slab2_middepth) then temp = tslab(2) else temp = tslab(1) + (tslab(2) - tslab(1))*(z - slab1_middepth)/(slab2_middepth - slab1_middepth) end if END FUNCTION interp_tslab_to_tsoil !======================================================================= END MODULE liquidbody