From e3212ecb99584a25237390fdf09d34671a738790 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Tue, 28 Nov 2023 12:23:24 -0800 Subject: [PATCH] Rebased w/d LTS branch --- .../mpas_ocn_time_integration_lts.F | 376 ++++++++++++------ .../mpas_ocn_time_integration_rk4.F | 2 +- .../src/shared/mpas_ocn_tidal_forcing.F | 4 +- .../src/shared/mpas_ocn_wetting_drying.F | 123 ++++-- 4 files changed, 337 insertions(+), 168 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_lts.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_lts.F index 073984852368..d123c69aecb5 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_lts.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_lts.F @@ -40,6 +40,7 @@ module ocn_time_integration_lts use ocn_equation_of_state use ocn_time_average_coupled use ocn_time_varying_forcing + use ocn_wetting_drying implicit none private @@ -153,7 +154,7 @@ subroutine ocn_time_integrator_lts(domain,dt)!{{{ layerThicknessFirstStage, & ! layer thickness at first stage of LTS layerThicknessSecondStage ! layer thickness at second stage of LTS - ! LTS objects + ! LTS Objects real (kind=RKIND) :: & dtFine, & ! fine dt, defined as dt / M alpha, alphaHat, beta, betaHat, gam, gamHat ! interpolation coefficients for the predictor step of LTS @@ -388,6 +389,10 @@ subroutine ocn_time_integrator_lts(domain,dt)!{{{ ! BEGIN LTS SCHEME !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + if (config_use_wetting_drying) then + wettingVelocityFactor(:, :) = 0.0_RKIND + end if + ! -------------------------------- COMPUTE SLOW TEND FOR ALL LTS REGIONS ------------------------------- call mpas_timer_start("lts slow tend computation") @@ -400,19 +405,10 @@ subroutine ocn_time_integrator_lts(domain,dt)!{{{ ! -------------------------------------- STEP ------------------------------------------ ! --- compute the first stage of LTS for interface layer 1, interface layer 2 and coarse - ! --- update halos for diagnostic variables. - call mpas_timer_start("lts diagnostic halo update") - call mpas_dmpar_field_halo_exch(domain, 'normalizedRelativeVorticityEdge') - if (config_mom_del4 > 0.0_RKIND) then - call mpas_dmpar_field_halo_exch(domain, 'divergence') - call mpas_dmpar_field_halo_exch(domain, 'relativeVorticity') - end if - call mpas_timer_stop("lts diagnostic halo update") - ! --- compute tendencies for thickness and fast tendencies only for velocity call mpas_timer_start("lts compute tendencies") - call ocn_lts_tends(statePool, LTSPool, tendPool, 1, 0, 1, 1, 1) + call ocn_lts_tends(statePool, LTSPool, tendPool, dt, 1, 0, 1, 1, 1) call mpas_timer_stop("lts compute tendencies") @@ -428,18 +424,21 @@ subroutine ocn_time_integrator_lts(domain,dt)!{{{ do iRegion =1,nRegions do ie = 1, nEdgesInLTSRegion(iRegion,2) iEdge = edgesInLTSRegion(iRegion,2,ie) - normalVelocityFirstStage(:,iEdge) = normalVelocityCur(:,iEdge) + dt * normalVelocityTend(:,iEdge) + normalVelocityFirstStage(:,iEdge) = ( normalVelocityCur(:,iEdge) + dt * normalVelocityTend(:,iEdge) ) & + * (1.0_RKIND - wettingVelocityFactor(:, iEdge)) end do end do ! --- coarse do ie = 1, nEdgesInLTSRegion(2,1) iEdge = edgesInLTSRegion(2,1,ie) - normalVelocityFirstStage(:,iEdge) = normalVelocityCur(:,iEdge) + dt * normalVelocityTend(:,iEdge) + normalVelocityFirstStage(:,iEdge) = ( normalVelocityCur(:,iEdge) + dt * normalVelocityTend(:,iEdge) ) & + * (1.0_RKIND - wettingVelocityFactor(:, iEdge)) end do ! --- fine layers close to interface layers do ie = 1, nEdgesInLTSRegion(1,3) iEdge = edgesInLTSRegion(1,3,ie) - normalVelocityFirstStage(:,iEdge) = normalVelocityCur(:,iEdge) + dt * normalVelocityTend(:,iEdge) + normalVelocityFirstStage(:,iEdge) = ( normalVelocityCur(:,iEdge) + dt * normalVelocityTend(:,iEdge) ) & + * (1.0_RKIND - wettingVelocityFactor(:, iEdge)) end do @@ -478,19 +477,11 @@ subroutine ocn_time_integrator_lts(domain,dt)!{{{ ! -------------------------------------- STEP ------------------------------------------ ! --- compute the second stage of LTS for interface layer 1, interface layer 2 and coarse - ! --- update halos for diagnostic variables. - call mpas_timer_start("lts diagnostic halo update") - call mpas_dmpar_field_halo_exch(domain, 'normalizedRelativeVorticityEdge') - if (config_mom_del4 > 0.0_RKIND) then - call mpas_dmpar_field_halo_exch(domain, 'divergence') - call mpas_dmpar_field_halo_exch(domain, 'relativeVorticity') - end if - call mpas_timer_stop("lts diagnostic halo update") call mpas_timer_start("lts compute tendencies") ! TENDENCIES COMPUTATION --- - call ocn_lts_tends(statePool, LTSPool, tendPool, 3, 0, 0, 1, 1) + call ocn_lts_tends(statePool, LTSPool, tendPool, dt, 3, 0, 0, 1, 1) call mpas_timer_stop("lts compute tendencies") @@ -505,18 +496,20 @@ subroutine ocn_time_integrator_lts(domain,dt)!{{{ ! --- coarse do ie = 1, nEdgesInLTSRegion(2,1) iEdge = edgesInLTSRegion(2,1,ie) - normalVelocitySecondStage(:,iEdge) = weightOld * normalVelocityCur(:,iEdge) & - + weightNew * normalVelocityFirstStage(:,iEdge) & - + weightTend * dt * normalVelocityTend(:,iEdge) + normalVelocitySecondStage(:,iEdge) = ( weightOld * normalVelocityCur(:,iEdge) & + + weightNew * normalVelocityFirstStage(:,iEdge) & + + weightTend * dt * normalVelocityTend(:,iEdge) ) & + * (1.0_RKIND - wettingVelocityFactor(:, iEdge)) end do ! --- interface layers do iRegion = 1, nRegions do ie = 1, nEdgesInLTSRegion(iRegion,2) iEdge = edgesInLTSRegion(iRegion,2,ie) - normalVelocitySecondStage(:,iEdge) = weightOld * normalVelocityCur(:,iEdge) & - + weightNew * normalVelocityFirstStage(:,iEdge) & - + weightTend * dt * normalVelocityTend(:,iEdge) + normalVelocitySecondStage(:,iEdge) = ( weightOld * normalVelocityCur(:,iEdge) & + + weightNew * normalVelocityFirstStage(:,iEdge) & + + weightTend * dt * normalVelocityTend(:,iEdge) ) & + * (1.0_RKIND - wettingVelocityFactor(:, iEdge)) end do end do @@ -590,20 +583,10 @@ subroutine ocn_time_integrator_lts(domain,dt)!{{{ call mpas_timer_stop("lts prognostic halo update") - ! --- update halos for diagnostic variables - call mpas_timer_start("lts diagnostic halo update") - call mpas_dmpar_field_halo_exch(domain, 'normalizedRelativeVorticityEdge') - if (config_mom_del4 > 0.0_RKIND) then - call mpas_dmpar_field_halo_exch(domain, 'divergence') - call mpas_dmpar_field_halo_exch(domain, 'relativeVorticity') - end if - call mpas_timer_stop("lts diagnostic halo update") - - call mpas_timer_start("lts compute tendencies") ! --- compute tendencies for thickness and fast tendencies only for velocity - call ocn_lts_tends(statePool, LTSPool, tendPool, 1, 1, 1, 0, 1) + call ocn_lts_tends(statePool, LTSPool, tendPool, dtFine, 1, 1, 1, 0, 1) call mpas_timer_stop("lts compute tendencies") @@ -620,22 +603,28 @@ subroutine ocn_time_integrator_lts(domain,dt)!{{{ do ie = 1, nEdgesInLTSRegion(1,2) iEdge = edgesInLTSRegion(1,2,ie) normalVelocityCur(:,iEdge) = normalVelocityNew(:,iEdge) - normalVelocityTendSum2nd(:,iEdge) = normalVelocityTendSum2nd(:,iEdge) + normalVelocityTend(:,iEdge) + normalVelocityTendSum2nd(:,iEdge) = ( 1.0_RKIND / (3.0_RKIND * dt * weightTendSum2nd) * normalVelocityCur(:,iEdge) & + + normalVelocityTendSum2nd(:,iEdge) + normalVelocityTend(:,iEdge) ) & + * (1.0_RKIND - wettingVelocityFactor(:, iEdge)) end do ! --- interface layer 2 (correction) do ie = 1, nEdgesInLTSRegion(2,2) iEdge = edgesInLTSRegion(2,2,ie) - normalVelocityTendSum2nd(:,iEdge) = normalVelocityTendSum2nd(:,iEdge) + normalVelocityTend(:,iEdge) + normalVelocityTendSum2nd(:,iEdge) = ( 1.0_RKIND / (3.0_RKIND * dt * weightTendSum2nd) * normalVelocityCur(:,iEdge) & + + normalVelocityTendSum2nd(:,iEdge) + normalVelocityTend(:,iEdge) ) & + * (1.0_RKIND - wettingVelocityFactor(:, iEdge)) end do ! --- fine (soln advancement) do ie = 1, nEdgesInLTSRegion(1,1) iEdge = edgesInLTSRegion(1,1,ie) - normalVelocityFirstStage(:,iEdge) = normalVelocityCur(:,iEdge) + dtFine * normalVelocityTend(:,iEdge) !soln update for the fine + normalVelocityFirstStage(:,iEdge) = ( normalVelocityCur(:,iEdge) + dtFine * normalVelocityTend(:,iEdge) ) & !soln update for the fine + * (1.0_RKIND - wettingVelocityFactor(:, iEdge)) end do do ie = 1, nEdgesInLTSRegion(1,3) iEdge = edgesInLTSRegion(1,3,ie) - normalVelocityFirstStage(:,iEdge) = normalVelocityCur(:,iEdge) + dtFine * normalVelocityTend(:,iEdge) !soln update for the fine + normalVelocityFirstStage(:,iEdge) = ( normalVelocityCur(:,iEdge) + dtFine * normalVelocityTend(:,iEdge) ) &!soln update for the fine + * (1.0_RKIND - wettingVelocityFactor(:, iEdge)) end do ! --- LAYER THICKNESS @@ -702,21 +691,10 @@ subroutine ocn_time_integrator_lts(domain,dt)!{{{ call mpas_dmpar_field_halo_exch(domain, 'layerThickness', timeLevel=3) call mpas_timer_stop("lts prognostic halo update") - - ! --- update halos for diagnostic variables - call mpas_timer_start("lts diagnostic halo update") - call mpas_dmpar_field_halo_exch(domain, 'normalizedRelativeVorticityEdge') - if (config_mom_del4 > 0.0_RKIND) then - call mpas_dmpar_field_halo_exch(domain, 'divergence') - call mpas_dmpar_field_halo_exch(domain, 'relativeVorticity') - end if - call mpas_timer_stop("lts diagnostic halo update") - - call mpas_timer_start("lts compute tendencies") ! TENDENCIES COMPUTATION --- - call ocn_lts_tends(statePool, LTSPool, tendPool, 3, 1, 1, 0, 1) + call ocn_lts_tends(statePool, LTSPool, tendPool, dtFine, 3, 1, 1, 0, 1) call mpas_timer_stop("lts compute tendencies") @@ -733,24 +711,30 @@ subroutine ocn_time_integrator_lts(domain,dt)!{{{ do ie = 1, nEdgesInLTSRegion(1,2) iEdge = edgesInLTSRegion(1,2,ie) normalVelocityFirstStage(:,iEdge) = normalVelocityNew(:,iEdge) - normalVelocityTendSum1st(:,iEdge) = normalVelocityTendSum1st(:,iEdge) + normalVelocityTend(:,iEdge) + normalVelocityTendSum1st(:,iEdge) = ( 1.0_RKIND / (3.0_RKIND * dt * weightTendSum1st) * normalVelocityCur(:,iEdge) & + + normalVelocityTendSum1st(:,iEdge) + normalVelocityTend(:,iEdge) ) & + * (1.0_RKIND - wettingVelocityFactor(:, iEdge)) end do ! --- interface layer 2 (correction) do ie = 1, nEdgesInLTSRegion(2,2) iEdge = edgesInLTSRegion(2,2,ie) - normalVelocityTendSum1st(:,iEdge) = normalVelocityTendSum1st(:,iEdge) + normalVelocityTend(:,iEdge) + normalVelocityTendSum1st(:,iEdge) = ( 1.0_RKIND / (3.0_RKIND * dt * weightTendSum1st) * normalVelocityCur(:,iEdge) & + + normalVelocityTendSum1st(:,iEdge) + normalVelocityTend(:,iEdge) ) & + * (1.0_RKIND - wettingVelocityFactor(:, iEdge)) end do ! --- fine (soln advancement) do ie = 1, nEdgesInLTSRegion(1,1) iEdge = edgesInLTSRegion(1,1,ie) - normalVelocitySecondStage(:,iEdge) = weightOld * normalVelocityCur(:,iEdge) + weightNew * normalVelocityFirstStage(:,iEdge) & - + weightTend * dtFine * normalVelocityTend(:,iEdge) !soln update for the fine + normalVelocitySecondStage(:,iEdge) = ( weightOld * normalVelocityCur(:,iEdge) + weightNew * normalVelocityFirstStage(:,iEdge) & + + weightTend * dtFine * normalVelocityTend(:,iEdge) ) & !soln update for the fine + * (1.0_RKIND - wettingVelocityFactor(:, iEdge)) end do do ie = 1, nEdgesInLTSRegion(1,3) iEdge = edgesInLTSRegion(1,3,ie) - normalVelocitySecondStage(:,iEdge) = weightOld * normalVelocityCur(:,iEdge) + weightNew * normalVelocityFirstStage(:,iEdge) & - + weightTend * dtFine * normalVelocityTend(:,iEdge) !soln update for the fine + normalVelocitySecondStage(:,iEdge) = ( weightOld * normalVelocityCur(:,iEdge) + weightNew * normalVelocityFirstStage(:,iEdge) & + + weightTend * dtFine * normalVelocityTend(:,iEdge) ) & !soln update for the fine + * (1.0_RKIND - wettingVelocityFactor(:, iEdge)) end do ! --- LAYER THICKNESS @@ -820,20 +804,9 @@ subroutine ocn_time_integrator_lts(domain,dt)!{{{ call mpas_timer_stop("lts prognostic halo update") - - ! --- update halos for diagnostic variables - call mpas_timer_start("lts diagnostic halo update") - call mpas_dmpar_field_halo_exch(domain, 'normalizedRelativeVorticityEdge') - if (config_mom_del4 > 0.0_RKIND) then - call mpas_dmpar_field_halo_exch(domain, 'divergence') - call mpas_dmpar_field_halo_exch(domain, 'relativeVorticity') - end if - call mpas_timer_stop("lts diagnostic halo update") - - call mpas_timer_start("lts compute tendencies") - call ocn_lts_tends(statePool, LTSPool, tendPool, 4, 1, 1, 0, 1) + call ocn_lts_tends(statePool, LTSPool, tendPool, dtFine, 4, 1, 1, 0, 1) call mpas_timer_stop("lts compute tendencies") @@ -850,24 +823,30 @@ subroutine ocn_time_integrator_lts(domain,dt)!{{{ do ie = 1, nEdgesInLTSRegion(1,2) iEdge = edgesInLTSRegion(1,2,ie) normalVelocitySecondStage(:,iEdge) = normalVelocityNew(:,iEdge) - normalVelocityTendSum3rd(:,iEdge) = normalVelocityTendSum3rd(:,iEdge) + normalVelocityTend(:,iEdge) + normalVelocityTendSum3rd(:,iEdge) = ( 1.0_RKIND / (3.0_RKIND * dt * weightTendSum3rd) * normalVelocityCur(:,iEdge) & + + normalVelocityTendSum3rd(:,iEdge) + normalVelocityTend(:,iEdge) ) & + * (1.0_RKIND - wettingVelocityFactor(:, iEdge)) end do ! --- interface layer 2 (correction) do ie = 1, nEdgesInLTSRegion(2,2) iEdge = edgesInLTSRegion(2,2,ie) - normalVelocityTendSum3rd(:,iEdge) = normalVelocityTendSum3rd(:,iEdge) + normalVelocityTend(:,iEdge) + normalVelocityTendSum3rd(:,iEdge) = ( 1.0_RKIND / (3.0_RKIND * dt * weightTendSum3rd) * normalVelocityCur(:,iEdge) & + + normalVelocityTendSum3rd(:,iEdge) + normalVelocityTend(:,iEdge) ) & + * (1.0_RKIND - wettingVelocityFactor(:, iEdge)) end do ! --- fine (soln advancement) do ie = 1, nEdgesInLTSRegion(1,1) iEdge = edgesInLTSRegion(1,1,ie) - normalVelocityCur(:,iEdge) = (1.0_RKIND / 3.0_RKIND) * normalVelocityCur(:,iEdge) + (2.0_RKIND / 3.0_RKIND) * normalVelocitySecondStage(:,iEdge) & - + (2.0_RKIND / 3.0_RKIND) * dtFine * normalVelocityTend(:,iEdge) !soln update for the fine + normalVelocityCur(:,iEdge) = ( (1.0_RKIND / 3.0_RKIND) * normalVelocityCur(:,iEdge) + (2.0_RKIND / 3.0_RKIND) * normalVelocitySecondStage(:,iEdge) & + + (2.0_RKIND / 3.0_RKIND) * dtFine * normalVelocityTend(:,iEdge) ) & !soln update for the fine + * (1.0_RKIND - wettingVelocityFactor(:, iEdge)) end do do ie = 1, nEdgesInLTSRegion(1,3) iEdge = edgesInLTSRegion(1,3,ie) - normalVelocityCur(:,iEdge) = (1.0_RKIND / 3.0_RKIND) * normalVelocityCur(:,iEdge) + (2.0_RKIND / 3.0_RKIND) * normalVelocitySecondStage(:,iEdge) & - + (2.0_RKIND / 3.0_RKIND) * dtFine * normalVelocityTend(:,iEdge) !soln update for the fine + normalVelocityCur(:,iEdge) = ( (1.0_RKIND / 3.0_RKIND) * normalVelocityCur(:,iEdge) + (2.0_RKIND / 3.0_RKIND) * normalVelocitySecondStage(:,iEdge) & + + (2.0_RKIND / 3.0_RKIND) * dtFine * normalVelocityTend(:,iEdge) ) & !soln update for the fine + * (1.0_RKIND - wettingVelocityFactor(:, iEdge)) end do ! --- LAYER THICKNESS @@ -943,18 +922,9 @@ subroutine ocn_time_integrator_lts(domain,dt)!{{{ ! --- compute the third stage of LTS for coarse ! - ! --- update halos for diagnostic variables - call mpas_timer_start("lts diagnostic halo update") - call mpas_dmpar_field_halo_exch(domain, 'normalizedRelativeVorticityEdge') - if (config_mom_del4 > 0.0_RKIND) then - call mpas_dmpar_field_halo_exch(domain, 'divergence') - call mpas_dmpar_field_halo_exch(domain, 'relativeVorticity') - end if - call mpas_timer_stop("lts diagnostic halo update") - call mpas_timer_start("lts compute tendencies") - call ocn_lts_tends(statePool, LTSPool, tendPool, 4, 0, 0, 1, 0) + call ocn_lts_tends(statePool, LTSPool, tendPool, dt, 4, 0, 0, 1, 0) call mpas_timer_stop("lts compute tendencies") @@ -968,8 +938,9 @@ subroutine ocn_time_integrator_lts(domain,dt)!{{{ ! --- coarse do ie = 1, nEdgesInLTSRegion(2,1) iEdge = edgesInLTSRegion(2,1,ie) - normalVelocityNew(:,iEdge) = (1.0_RKIND / 3.0_RKIND) * normalVelocityCur(:,iEdge) + (2.0_RKIND / 3.0_RKIND) * normalVelocitySecondStage(:,iEdge) & - + (2.0_RKIND / 3.0_RKIND) * dt * normalVelocityTend(:,iEdge) + normalVelocityNew(:,iEdge) = ( (1.0_RKIND / 3.0_RKIND) * normalVelocityCur(:,iEdge) + (2.0_RKIND / 3.0_RKIND) * normalVelocitySecondStage(:,iEdge) & + + (2.0_RKIND / 3.0_RKIND) * dt * normalVelocityTend(:,iEdge) ) & + * (1.0_RKIND - wettingVelocityFactor(:, iEdge)) end do ! --- LAYER THICKNESS @@ -999,8 +970,12 @@ subroutine ocn_time_integrator_lts(domain,dt)!{{{ do ie = 1, nEdgesInLTSRegion(iRegion,2) iEdge = edgesInLTSRegion(iRegion,2,ie) - normalVelocityNew(:,iEdge) = normalVelocityCur(:,iEdge) + weightTendSum2nd * dtFine * normalVelocityTendSum2nd(:,iEdge) & - + weightTendSum1st * dtFine * normalVelocityTendSum1st(:,iEdge) + weightTendSum3rd * dtFine * normalVelocityTendSum3rd(:,iEdge) + ! The corrections are computed like this (taking normalVelocityTendSum2nd as an example): + ! normalVelocityTendSum2nd(:,iEdge) = ( 1 / (3*dt*weightTendSum2nd) * normalVelocityCur(:,iEdge) & + ! + normalVelocityTendSum2nd(:,iEdge) + normalVelocityTend(:,iEdge) ) + normalVelocityNew(:,iEdge) = weightTendSum2nd * dtFine * normalVelocityTendSum2nd(:,iEdge) & + + weightTendSum1st * dtFine * normalVelocityTendSum1st(:,iEdge) & + + weightTendSum3rd * dtFine * normalVelocityTendSum3rd(:,iEdge) end do ! --- LAYER THICKNESS @@ -1045,6 +1020,23 @@ subroutine ocn_time_integrator_lts(domain,dt)!{{{ call mpas_timer_start("lts cleanup phase") + ! verify that cells are not dry at conclusion of time step + if (config_use_wetting_drying) then + call mpas_timer_start("lts check wet cells") + + ! ensure existing layerThickness is valid + if (config_verify_not_dry) then + call ocn_wetting_drying_verify(block, config_drying_min_cell_height, err) + end if + + call mpas_timer_stop("lts check wet cells") + end if + + ! direct application of tidal boundary condition + if (config_use_tidal_forcing .and. trim(config_tidal_forcing_type) == 'direct') then + call ocn_lts_apply_direct_tidal_bc(statePool, forcingPool, verticalMeshPool, 2) + end if + if (config_prescribe_velocity) then do iEdge = 1, nEdgesAll normalVelocityNew(:, iEdge) = normalVelocityCur(:, iEdge) @@ -1089,7 +1081,6 @@ subroutine ocn_time_integrator_lts(domain,dt)!{{{ SSHGradient(indexSSHGradientMeridional, iCell) = gradSSHMeridional(iCell) end do - ! DIAGNOSTICS UPDATE --- call ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, verticalMeshPool, scratchPool, tracersPool, 2) @@ -1131,7 +1122,7 @@ subroutine ocn_time_integration_lts_init(domain)!{{{ type (domain_type), intent(inout) :: domain type (block_type), pointer :: block - type (mpas_pool_type), pointer :: LTSPool, meshPool + type (mpas_pool_type), pointer :: LTSPool integer, dimension(:), allocatable :: isLTSRegionEdgeAssigned integer :: i, iCell, iEdge, iRegion, coarseRegions, fineRegions, fineRegionsM1 integer, dimension(:), pointer :: LTSRegion @@ -1144,7 +1135,6 @@ subroutine ocn_time_integration_lts_init(domain)!{{{ block => domain % blocklist call mpas_pool_get_subpool(block % structs, 'LTS', LTSPool) - call mpas_pool_get_subpool(block % structs, ',mesh', meshPool) call mpas_pool_get_array(LTSPool, 'LTSRegion', LTSRegion) call mpas_pool_get_array(LTSPool, 'cellsInLTSRegion', cellsInLTSRegion) @@ -1191,7 +1181,7 @@ subroutine ocn_time_integration_lts_init(domain)!{{{ ! by the cell in the LTS region closest to the fine region, ! see Figure 3 in "Conservative explicit local time-stepping schemes for ! the shallow water equations" by Hoang et al. (halo edges however - ! are owned by whatever processor they are initially assigned to) + ! are owned by whatever process they are initially assigned to) allocate(isLTSRegionEdgeAssigned(nEdgesOwned)) isLTSRegionEdgeAssigned(:) = 0 @@ -1248,7 +1238,7 @@ end subroutine ocn_time_integration_lts_init!}}} !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine ocn_lts_tends(statePool, LTSPool, tendPool, timeLevelIn, computeOnFineBig, computeOnFineSmall, computeOnCoarse, computeOnInterface)!{{{ + subroutine ocn_lts_tends(statePool, LTSPool, tendPool, dt, timeLevelIn, computeOnFineBig, computeOnFineSmall, computeOnCoarse, computeOnInterface)!{{{ integer, intent(in) :: timeLevelIn, computeOnFineBig, computeOnFineSmall, computeOnCoarse, computeOnInterface @@ -1258,6 +1248,8 @@ subroutine ocn_lts_tends(statePool, LTSPool, tendPool, timeLevelIn, computeOnFin type (mpas_pool_type), intent(in) :: & LTSPool !< Input: LTS data + real (kind=RKIND), intent(in) :: dt + type (mpas_pool_type), intent(inout) :: & tendPool !< Input: tendency variables @@ -1266,8 +1258,8 @@ subroutine ocn_lts_tends(statePool, LTSPool, tendPool, timeLevelIn, computeOnFin real (kind=RKIND), dimension(:), pointer :: ssh real (kind=RKIND), dimension(:,:), pointer :: normalVelocityTend, layerThicknessTend, normalVelocity, layerThickness - integer :: iEdge, cell1, cell2, k, ie, iRegion, nRegions, ic, i, iCell, kmin, kmax - real (kind=RKIND) :: invdcEdge, betaSelfAttrLoad, flux, ssh_sal_on + integer :: iEdge, cell1, cell2, k, ie, iRegion, nRegions, ic, i, iCell, nCells, kmin, kmax + real (kind=RKIND) :: invdcEdge, betaSelfAttrLoad, flux, ssh_sal_on, tidal_pot_for_on betaSelfAttrLoad = config_self_attraction_and_loading_beta nRegions = 2 @@ -1289,27 +1281,64 @@ subroutine ocn_lts_tends(statePool, LTSPool, tendPool, timeLevelIn, computeOnFin else ssh_sal_on = 0.0_RKIND endif + if (config_use_tidal_potential_forcing) then + tidal_pot_for_on = 1.0_RKIND + else + tidal_pot_for_on = 0.0_RKIND + end if ! inline computation of the diagnostics ! layerThickEdgeFlux - do iEdge = 1, nEdgesAll - kmin = minLevelEdgeBot(iEdge) - kmax = maxLevelEdgeTop(iEdge) - cell1 = cellsOnEdge(1,iEdge) - cell2 = cellsOnEdge(2,iEdge) - do k = 1,nVertLevels - ! initialize layerThicknessEdgeMean to avoid divide by - ! zero and NaN problems. - layerThickEdgeFlux(k,iEdge) = -1.0e34_RKIND + if(config_thickness_flux_type == 'centered') then + + do iEdge = 1, nEdgesAll + kmin = minLevelEdgeBot(iEdge) + kmax = maxLevelEdgeTop(iEdge) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + do k = 1,nVertLevels + ! initialize layerThicknessEdgeMean to avoid divide by + ! zero and NaN problems. + layerThickEdgeFlux(k,iEdge) = -1.0e34_RKIND + end do + do k = kmin,kmax + ! central differenced + layerThickEdgeFlux(k,iEdge) = 0.5_RKIND * & + (layerThickness(k,cell1) + & + layerThickness(k,cell2)) + end do end do - do k = kmin,kmax - ! central differenced - layerThickEdgeFlux(k,iEdge) = 0.5_RKIND * & - (layerThickness(k,cell1) + & - layerThickness(k,cell2)) + + else if(config_thickness_flux_type == 'upwind') then + + do iEdge = 1, nEdgesAll + kmin = minLevelEdgeBot(iEdge) + kmax = maxLevelEdgeTop(iEdge) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + do k=1,nVertLevels + ! initialize layerThicknessEdgeFlux to avoid divide by + ! zero and NaN problems. + layerThickEdgeFlux(k,iEdge) = -1.0e34_RKIND + end do + do k = kmin,kmax + if (normalVelocity(k,iEdge) > 0.0_RKIND) then + layerThickEdgeFlux(k,iEdge) = layerThickness(k,cell1) + elseif (normalVelocity(k,iEdge) < 0.0_RKIND) then + layerThickEdgeFlux(k,iEdge) = layerThickness(k,cell2) + else + layerThickEdgeFlux(k,iEdge) = max(layerThickness(k,cell1), layerThickness(k,cell2)) + end if + end do end do - end do + + else + + call mpas_log_write('ERROR: config_thickness_flux_type selected not implemented for LTS', messageType=MPAS_LOG_CRIT) + call abort + + end if ! ssh do iCell = 1, nCellsAll @@ -1322,10 +1351,13 @@ subroutine ocn_lts_tends(statePool, LTSPool, tendPool, timeLevelIn, computeOnFin ssh(iCell) = zTop(minLevelCell(iCell),iCell) end do - ! inline computation of the tendencies + if (config_use_wetting_drying) then + call ocn_lts_wd(statePool, dt, timeLevelIn) + end if + ! inline computation of the tendencies normalVelocityTend(:,:) = 0.0_RKIND - layerThicknessTend(:, :) = 0.0_RKIND + layerThicknessTend(:,:) = 0.0_RKIND ! interface if (computeOnInterface == 1) then @@ -1339,9 +1371,9 @@ subroutine ocn_lts_tends(statePool, LTSPool, tendPool, timeLevelIn, computeOnFin kMin = minLevelEdgeBot(iEdge) kMax = maxLevelEdgeTop(iEdge) do k=kMin,kMax - normalVelocityTend(k,iEdge) = normalVelocityTend(k,iEdge) & - - edgeMask(k,iEdge) * invdcEdge * ( gravity * ( (ssh(cell2) - ssh(cell1)) & - - (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad * (ssh(cell2) - ssh(cell1)) ) ) + normalVelocityTend(k,iEdge) = normalVelocityTend(k,iEdge) & + - edgeMask(k,iEdge) * invdcEdge * ( gravity * ( (ssh(cell2) - ssh(cell1)) & + - tidal_pot_for_on * (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad * (ssh(cell2) - ssh(cell1)) ) ) end do end do ! thickness tendency @@ -1350,7 +1382,7 @@ subroutine ocn_lts_tends(statePool, LTSPool, tendPool, timeLevelIn, computeOnFin do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) - flux = normalVelocity(k, iEdge) * dvEdge(iEdge) * layerThickEdgeFlux(k, iEdge) + flux = normalVelocity(k, iEdge) * dvEdge(iEdge) * layerThickEdgeFlux(k, iEdge) * (1.0_RKIND - wettingVelocityFactor(k, iEdge)) layerThicknessTend(k, iCell) = layerThicknessTend(k, iCell) + edgeSignOnCell(i, iCell) * flux * invAreaCell(iCell) end do end do @@ -1370,7 +1402,7 @@ subroutine ocn_lts_tends(statePool, LTSPool, tendPool, timeLevelIn, computeOnFin do k=kMin,kMax normalVelocityTend(k,iEdge) = normalVelocityTend(k,iEdge) & - edgeMask(k,iEdge) * invdcEdge * ( gravity * ( (ssh(cell2) - ssh(cell1)) & - - (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad * (ssh(cell2) - ssh(cell1)) ) ) + - tidal_pot_for_on * (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad * (ssh(cell2) - ssh(cell1)) ) ) end do end do ! thickness tendency @@ -1379,7 +1411,7 @@ subroutine ocn_lts_tends(statePool, LTSPool, tendPool, timeLevelIn, computeOnFin do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) - flux = normalVelocity(k, iEdge) * dvEdge(iEdge) * layerThickEdgeFlux(k, iEdge) + flux = normalVelocity(k, iEdge) * dvEdge(iEdge) * layerThickEdgeFlux(k, iEdge) * (1.0_RKIND - wettingVelocityFactor(k, iEdge)) layerThicknessTend(k, iCell) = layerThicknessTend(k, iCell) + edgeSignOnCell(i, iCell) * flux * invAreaCell(iCell) end do end do @@ -1398,7 +1430,7 @@ subroutine ocn_lts_tends(statePool, LTSPool, tendPool, timeLevelIn, computeOnFin do k=kMin,kMax normalVelocityTend(k,iEdge) = normalVelocityTend(k,iEdge) & - edgeMask(k,iEdge) * invdcEdge * ( gravity * ( (ssh(cell2) - ssh(cell1)) & - - (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad * (ssh(cell2) - ssh(cell1)) ) ) + - tidal_pot_for_on * (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad * (ssh(cell2) - ssh(cell1)) ) ) end do end do ! thickness tendency @@ -1407,7 +1439,7 @@ subroutine ocn_lts_tends(statePool, LTSPool, tendPool, timeLevelIn, computeOnFin do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) - flux = normalVelocity(k, iEdge) * dvEdge(iEdge) * layerThickEdgeFlux(k, iEdge) + flux = normalVelocity(k, iEdge) * dvEdge(iEdge) * layerThickEdgeFlux(k, iEdge) * (1.0_RKIND - wettingVelocityFactor(k, iEdge)) layerThicknessTend(k, iCell) = layerThicknessTend(k, iCell) + edgeSignOnCell(i, iCell) * flux * invAreaCell(iCell) end do end do @@ -1424,9 +1456,9 @@ subroutine ocn_lts_tends(statePool, LTSPool, tendPool, timeLevelIn, computeOnFin kMin = minLevelEdgeBot(iEdge) kMax = maxLevelEdgeTop(iEdge) do k=kMin,kMax - normalVelocityTend(k,iEdge) = normalVelocityTend(k,iEdge) & - - edgeMask(k,iEdge) * invdcEdge * ( gravity * ( (ssh(cell2) - ssh(cell1)) & - - (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad * (ssh(cell2) - ssh(cell1)) ) ) + normalVelocityTend(k,iEdge) = normalVelocityTend(k,iEdge) & + - edgeMask(k,iEdge) * invdcEdge * ( gravity * ( (ssh(cell2) - ssh(cell1)) & + - tidal_pot_for_on * (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad * (ssh(cell2) - ssh(cell1)) ) ) end do end do ! thickness tendency @@ -1435,7 +1467,7 @@ subroutine ocn_lts_tends(statePool, LTSPool, tendPool, timeLevelIn, computeOnFin do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) - flux = normalVelocity(k, iEdge) * dvEdge(iEdge) * layerThickEdgeFlux(k, iEdge) + flux = normalVelocity(k, iEdge) * dvEdge(iEdge) * layerThickEdgeFlux(k, iEdge) * (1.0_RKIND - wettingVelocityFactor(k, iEdge)) layerThicknessTend(k, iCell) = layerThicknessTend(k, iCell) + edgeSignOnCell(i, iCell) * flux * invAreaCell(iCell) end do end do @@ -1446,6 +1478,90 @@ end subroutine ocn_lts_tends!}}} !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine ocn_lts_wd(statePool, dt, timeLevelIn) + + type (mpas_pool_type), intent(in) :: & + statePool !< Input: state variables + + real (kind=RKIND), intent(in) :: dt + + integer, intent(in) :: timeLevelIn + + real (kind=RKIND), dimension(:, :), pointer :: layerThickness + real (kind=RKIND), dimension(:, :), pointer :: normalVelocity + real (kind=RKIND) :: divOutFlux, layerThicknessTmp + + integer :: iEdge, k, iCell, i + + call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, timeLevelIn) + call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, timeLevelIn) + + do iCell = 1, nCellsAll + do k = minLevelCell(iCell), maxLevelCell(iCell) + divOutFlux = 0.0_RKIND + layerThicknessTmp = layerThickness(k, iCell) + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + if (k <= maxLevelEdgeTop(iEdge) .and. k >= minLevelEdgeBot(iEdge)) then + ! only consider divergence flux leaving the cell + if ( normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) < 0.0_RKIND ) then + divOutFlux = divOutFlux + & + normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) * & + layerThickEdgeFlux(k, iEdge) * dvEdge(iEdge) * & + invAreaCell(iCell) + end if + end if + end do + layerThicknessTmp = layerThicknessTmp + dt * divOutFlux + + call ocn_wetting_velocity_factor_on_cell_edges(wettingVelocityFactor, layerThicknessTmp, normalVelocity, iCell, k) + + end do + end do + + call mpas_dmpar_exch_halo_field(wettingVelocityField) + + end subroutine ocn_lts_wd!}}} + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + subroutine ocn_lts_apply_direct_tidal_bc(statePool, forcingPool, verticalMeshPool, timeLevelIn)!{{{ + + type (mpas_pool_type), intent(inout) :: statePool + + type (mpas_pool_type), intent(in) :: forcingPool, verticalMeshPool + + integer, intent(in) :: timeLevelIn + + integer :: iCell, k + + real (kind=RKIND) :: totalDepth + real (kind=RKIND), dimension(:,:), pointer :: layerThicknessNew, restingThickness + real (kind=RKIND), dimension(:), pointer :: tidalInputMask, tidalBCValue + + call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessNew, timeLevelIn) + call mpas_pool_get_array(forcingPool, 'tidalInputMask', tidalInputMask) + call mpas_pool_get_array(forcingPool, 'tidalBCValue', tidalBCValue) + call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness) + + do iCell=1, nCellsOwned + if (tidalInputMask(iCell) == 1.0_RKIND) then + ! compute total depth for relative thickness contribution + totalDepth = 0.0_RKIND + do k = minLevelCell(iCell), maxLevelCell(iCell) + totalDepth = totalDepth + restingThickness(k,iCell) + end do + ! only modify layer thicknesses on tidal boundary + do k = minLevelCell(iCell), maxLevelCell(iCell) + layerThicknessNew(k, iCell) = tidalInputMask(iCell)*(tidalBCValue(iCell) + bottomDepth(iCell))*(restingThickness(k,iCell)/totalDepth) + end do + end if + end do + + end subroutine ocn_lts_apply_direct_tidal_bc!}}} + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + real function sphere_distance(lat1, lon1, lat2, lon2, radius) ! Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F index 8a3c633536c0..af583448157c 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F @@ -744,7 +744,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ #endif call ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, verticalMeshPool, scratchPool, tracersPool, 2) - ! Update the effective desnity in land ice if we're coupling to land ice + ! Update the effective density in land ice if we're coupling to land ice call ocn_effective_density_in_land_ice_update(forcingPool, & statePool, err) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tidal_forcing.F b/components/mpas-ocean/src/shared/mpas_ocn_tidal_forcing.F index 98a6f99e77bd..56e2e03f26e7 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tidal_forcing.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tidal_forcing.F @@ -109,7 +109,7 @@ subroutine ocn_tidal_forcing_layer_thickness(forcingPool, layerThicknessTend, er err = 0 - if ( .not. tidalfluxOn ) return + if ( .not. tidalFluxOn ) return call mpas_timer_start("tidal thickness tendency") @@ -219,7 +219,7 @@ subroutine ocn_tidal_forcing_build_array(domain, meshPool, forcingPool, statePoo character (len=StrKIND), pointer :: xtime - if ( .not. tidalfluxOn ) return + if ( .not. tidalFluxOn ) return call mpas_pool_get_array(statePool, 'ssh', ssh, 1) call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 1) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F b/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F index 0ef4c7851154..92cb65b2bbaf 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F @@ -51,6 +51,7 @@ module ocn_wetting_drying !-------------------------------------------------------------------- public :: ocn_wetting_drying_verify, ocn_prevent_drying_rk4 + public :: ocn_wetting_velocity_factor_on_cell_edges !-------------------------------------------------------------------- ! @@ -353,16 +354,13 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness integer :: cell1, cell2, iEdge, iCell, k, i - real (kind=RKIND) :: divFlux, divOutFlux + real (kind=RKIND) :: divOutFlux real (kind=RKIND) :: layerThickness - real (kind=RKIND) :: hCrit, hRampMin, hEdgeTotal character (len=100) :: log_string err = 0 - hCrit = config_drying_min_cell_height - if (.not. config_zero_drying_velocity) return ! need predicted transport velocity to limit drying flux @@ -386,37 +384,7 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness end do layerThickness = layerThickness + dt * divOutFlux - ! if layer thickness is too small, limit divergence flux outwards with - ! opposite velocity - if (layerThickness <= & - hCrit + config_drying_safety_height) then - do i = 1, nEdgesOnCell(iCell) - iEdge = edgesOnCell(i, iCell) - if (k <= maxLevelEdgeBot(iEdge) .and. k >= minLevelEdgeTop(iEdge)) then - if ( normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) <= 0.0_RKIND ) then - wettingVelocityFactor(k, iEdge) = 1.0_RKIND - end if - end if - end do - elseif (config_zero_drying_velocity_ramp .and. & - (layerThickness > & - hCrit + config_drying_safety_height) .and. & - (layerThickness <= config_zero_drying_velocity_ramp_hmax)) then - - hRampMin = config_zero_drying_velocity_ramp_hmin - ! Following O'Dea et al. (2021), if total upwinded wct is less than - ! 2*critical thickness, apply damping at each edge - do i = 1, nEdgesOnCell(iCell) - iEdge = edgesOnCell(i, iCell) - if (k <= maxLevelEdgeBot(iEdge) .and. k >= minLevelEdgeTop(iEdge)) then - if ( normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) <= 0.0_RKIND ) then - wettingVelocityFactor(k, iEdge) = 1.0_RKIND - & - tanh(50.0_RKIND * (layerThickness - hRampMin)/hRampMin) - end if - end if - end do - - end if + call ocn_wetting_velocity_factor_on_cell_edges(wettingVelocityFactor, layerThickness, normalVelocity, iCell, k) end do end do @@ -425,6 +393,91 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness end subroutine ocn_wetting_drying_wettingVelocity !}}} +!*********************************************************************** +! +! routine ocn_wetting_velocity_factor_on_cell_edges +! +!> \brief Computes wettingVelocityFactor at edges of iCell +!> \author Giacomo Capodaglio +!> \date 09/06/2023 +!> \details +!> This routine computes wettingVelocityFactor at the edges of iCell +! +!----------------------------------------------------------------------- + subroutine ocn_wetting_velocity_factor_on_cell_edges(wettingVelocityFactor, layerThicknessTmp, normalVelocity, iCell, k)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), intent(in) :: & + layerThicknessTmp !< Input: layer thickness to use for limiting + + real (kind=RKIND), dimension(:,:), intent(in) :: & + normalVelocity !< Input: transport + + integer, intent(in) :: & + iCell !< Input: wettingVelocityFactor is computed at the edges of iCell + + integer, intent(in) :: & + k !< Input: wettingVelocityFactor is computed at layer k + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + wettingVelocityFactor !< Input/Output: velocity wettingVelocityFactor + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: i, iEdge + + real (kind=RKIND) :: hCrit, hRampMin + + hCrit = config_drying_min_cell_height + + ! if layer thickness is too small, limit divergence flux outwards with + ! opposite velocity + if (layerThicknessTmp <= & + hCrit + config_drying_safety_height) then + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + if (k <= maxLevelEdgeBot(iEdge) .and. k >= minLevelEdgeTop(iEdge)) then + if ( normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) <= 0.0_RKIND ) then + wettingVelocityFactor(k, iEdge) = 1.0_RKIND + end if + end if + end do + elseif (config_zero_drying_velocity_ramp .and. & + (layerThicknessTmp > & + hCrit + config_drying_safety_height) .and. & + (layerThicknessTmp <= config_zero_drying_velocity_ramp_hmax)) then + + hRampMin = config_zero_drying_velocity_ramp_hmin + ! Following O'Dea et al. (2021), if total upwinded wct is less than + ! 2*critical thickness, apply damping at each edge + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + if (k <= maxLevelEdgeBot(iEdge) .and. k >= minLevelEdgeTop(iEdge)) then + if ( normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) <= 0.0_RKIND ) then + wettingVelocityFactor(k, iEdge) = 1.0_RKIND - & + tanh(50.0_RKIND * (layerThicknessTmp - hRampMin)/hRampMin) + end if + end if + end do + + end if + + end subroutine ocn_wetting_velocity_factor_on_cell_edges!}}} end module ocn_wetting_drying