MODULE lake_ocean_transformation !----------------------------------------------------------------------- ! NAME ! lake_ocean_transformation ! ! DESCRIPTION ! State and helpers for Generic PEM fields that are only needed when ! handing control back to the PCM. ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use numerics, only: dp, di, k4, eps ! DECLARATION ! ----------- implicit none ! PARAMETERS ! ---------- real(dp), parameter, private :: mx_eau_sol_PCM = 150._dp ! From phygeneric/watercommon_h.F90 [kg/m2] real(dp), parameter, private :: 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] real(dp), protected :: slab_depth_PCM ! Slab ocean depth of the PCM wind-mixed layer [m] ! 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_liq_PCM real(dp), dimension(:), allocatable :: qsurf_h2o_liq_PCM contains !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !======================================================================= SUBROUTINE ini_lake_ocean_transformation() !----------------------------------------------------------------------- ! Allocate PCM handoff fields. !----------------------------------------------------------------------- use geometry, only: ngrid implicit none allocate(tslab_PCM(ngrid,2)) allocate(tsea_ice_PCM(ngrid)) allocate(tice_PCM(ngrid)) allocate(q_peren_h2o_liq_PCM(ngrid)) allocate(qsurf_h2o_liq_PCM(ngrid)) tslab_PCM(:,:) = 0._dp tsea_ice_PCM(:) = 0._dp tice_PCM(:) = 0._dp q_peren_h2o_liq_PCM(:) = 0._dp qsurf_h2o_liq_PCM(:) = 0._dp END SUBROUTINE ini_lake_ocean_transformation !======================================================================= !======================================================================= SUBROUTINE end_lake_ocean_transformation() !----------------------------------------------------------------------- ! Deallocate PCM handoff fields. !----------------------------------------------------------------------- implicit none 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_liq_PCM)) deallocate(q_peren_h2o_liq_PCM) if (allocated(qsurf_h2o_liq_PCM)) deallocate(qsurf_h2o_liq_PCM) END SUBROUTINE end_lake_ocean_transformation !======================================================================= !======================================================================= SUBROUTINE set_slab_depth_PCM(slab_depth_in) !----------------------------------------------------------------------- ! Store the PCM slab-ocean upper-layer depth read from callphys.def. !----------------------------------------------------------------------- implicit none real(dp), intent(in) :: slab_depth_in slab_depth_PCM = slab_depth_in END SUBROUTINE set_slab_depth_PCM !======================================================================= !======================================================================= SUBROUTINE tranform_lake_ocean4PCM(peren_h2o_liq_frac4PCM,is_ocean4PCM, & lakeice4PCM,seaice4PCM,seaice_frac4PCM, & tsurf4PCM,tsoil4PCM,tslab4PCM, & tsea_ice4PCM,tice4PCM, & q_peren_h2o_liq4PCM,qsurf_h2o_liq4PCM) !----------------------------------------------------------------------- ! Normalize all PCM handoff fields in cells whose water coverage crossed ! the lake/ocean threshold during the PEM leg. ! We follow the PCM's definitions of tslab, tsurf, etc. !----------------------------------------------------------------------- use geometry, only: ngrid, nslope, nsoil_PCM, cell_area use slopes, only: iflat use hydrology, only: is_ocean_PCM, lake_volume, h2o_liq_density use soil, only: mlayer use display, only: print_msg, LVL_NFO use utility, only: int2str implicit none 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_h2o_liq4PCM, qsurf_h2o_liq4PCM 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 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_h2o_liq4PCM(ig) q_peren_h2o_liq4PCM(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 lakeice4PCM(ig) = seaice4PCM(ig)*seaice_frac4PCM(ig)/peren_h2o_liq_frac4PCM(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_h2o_liq4PCM(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 = mx_eau_sol_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 !======================================================================= !======================================================================= FUNCTION interp_tsoil_to_tslab(ig,islope,tsoil,tsurf,z) RESULT(temp) !----------------------------------------------------------------------- ! Interpolate the outgoing PCM soil profile to a slab depth. ! We know tsoil and tsurf, and we are mapping onto 2 slab depths !----------------------------------------------------------------------- use geometry, only: nsoil_PCM use soil, only: mlayer implicit none integer(di), intent(in) :: ig, islope real(dp), dimension(:,:,:), intent(in) :: tsoil real(dp), dimension(:,:), intent(in) :: tsurf real(dp), intent(in) :: z real(dp) :: temp integer(di) :: isoil real(dp) :: upper_depth, lower_depth real(dp) :: upper_temp, lower_temp 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) !----------------------------------------------------------------------- ! Interpolate the incoming slab profile onto a PCM soil depth. ! We know 3 points (tsurf, tslab1, and tslab2) and we are mapping onto soil depths !----------------------------------------------------------------------- implicit none real(dp), dimension(:), intent(in) :: tslab real(dp), intent(in) :: tsurf, z real(dp) :: temp real(dp) :: slab1_middepth, slab2_middepth 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 lake_ocean_transformation