diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F index a719bd24e135..46283c4e4b29 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F @@ -208,8 +208,9 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ real (kind=RKIND), dimension(:), allocatable :: & btrvel_temp, & sshTemp, & + deltaSSH, & wctEdge, &! flux water column thickness at edge - bottomDepthEdge + baroclinicThickness ! State Array Pointers real (kind=RKIND), dimension(:), pointer :: & @@ -276,6 +277,10 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ ! Forcing pool real (kind=RKIND), dimension(:), pointer :: tidalPotentialEta + real (kind=RKIND) :: totalDepth + real (kind=RKIND), dimension(:), pointer :: tidalInputMask, tidalBCValue + real (kind=RKIND), dimension(:,:), pointer :: restingThickness + real (kind=RKIND), dimension(:), pointer :: & seaIcePressure, atmosphericPressure @@ -406,7 +411,7 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracersNew, 2) - allocate(bottomDepthEdge(nEdgesAll+1)) + allocate(baroclinicThickness(nEdgesAll+1)) if (config_transport_tests_flow_id > 0) then ! This is a transport test. Write advection velocity from prescribed @@ -908,20 +913,18 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ !$omp end do !$omp end parallel - bottomDepthEdge(:) = 0.0_RKIND + baroclinicThickness(:) = 0.0_RKIND !$omp parallel !$omp do schedule(runtime) & - !$omp private(cell1, cell2, k, thicknessSum) + !$omp private(cell1, cell2, k) do iEdge = 1, nEdgesHalo(config_num_halos+1) cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) - thicknessSum = layerThickEdgeFlux(minLevelEdgeBot(iEdge),iEdge) + baroclinicThickness(iEdge) = layerThickEdgeFlux(minLevelEdgeBot(iEdge),iEdge) do k = minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge) - thicknessSum = thicknessSum + & + baroclinicThickness(iEdge) = baroclinicThickness(iEdge) + & layerThickEdgeFlux(k,iEdge) enddo - bottomDepthEdge(iEdge) = thicknessSum & - - 0.5_RKIND*(sshNew(cell1) + sshNew(cell2)) enddo ! iEdge !$omp end do !$omp end parallel @@ -1015,6 +1018,8 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ btrvel_temp(:) = 0.0_RKIND allocate(wctEdge(nEdgesAll+1)) wctEdge(:) = 0.0_RKIND + allocate(deltaSSH(nCellsAll)) + deltaSSH(:) = 0.0_RKIND cellHaloComputeCounter = 0 edgeHaloComputeCounter = 0 @@ -1089,6 +1094,30 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ end if ! tidal potential forcing + ! direct application of tidal boundary condition + if (config_use_tidal_forcing .and. trim(config_tidal_forcing_type) == 'direct') then + call mpas_pool_get_array(forcingPool, & + 'sshSubcycleCurWithTides', & + sshSubcycleCurWithTides) + call mpas_pool_get_array(forcingPool, 'tidalInputMask', tidalInputMask) + call mpas_pool_get_array(forcingPool, 'tidalBCValue', tidalBCValue) + !$omp parallel + !$omp do schedule(runtime) + do iCell=1, nCellsAll + ! boolean mask for now, could generalize to tappered sponge layer + ! only modify layer thicknesses on tidal boundary + do k = minLevelCell(iCell), maxLevelCell(iCell) + sshSubcycleCurWithTides(iCell) = sshSubcycleCur(iCell) + & + (tidalInputMask(iCell)*(tidalBCValue(iCell) - sshSubcycleCur(iCell))) + end do + end do + !$omp end do + !$omp end parallel + call mpas_pool_get_array(forcingPool, & + 'sshSubcycleCurWithTides', & + sshSubcycleCur) + end if + if (edgeHaloComputeCounter <= 1) then nEdges = nEdgesOwned else @@ -1187,8 +1216,16 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ ! Compute barotropic thickness at the edge. Note this ! matches the sum of baroclinic thicknesses at this edge. - call ocn_diagnostic_solve_wctEdge(btrvel_temp, sshSubcycleCur, & - bottomDepthEdge, wctEdge) + !$omp parallel + !$omp do schedule(runtime) + do iCell = 1, nCellsAll + deltaSSH(iCell) = sshSubcycleCur(iCell) - sshNew(iCell) + end do + !$omp end do + !$omp end parallel + + call ocn_diagnostic_solve_wctEdge(btrvel_temp, deltaSSH, & + baroclinicThickness, wctEdge) !$omp parallel !$omp do schedule(runtime) & @@ -1316,6 +1353,38 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ sshSubcycleNew) end if ! tidal potential forcing + ! direct application of tidal boundary condition + if (config_use_tidal_forcing .and. trim(config_tidal_forcing_type) == 'direct') then + call mpas_pool_get_array(forcingPool, & + 'sshSubcycleCurWithTides', & + sshSubcycleCurWithTides) + call mpas_pool_get_array(forcingPool, & + 'sshSubcycleNewWithTides', & + sshSubcycleNewWithTides) + call mpas_pool_get_array(forcingPool, 'tidalInputMask', tidalInputMask) + call mpas_pool_get_array(forcingPool, 'tidalBCValue', tidalBCValue) + !$omp parallel + !$omp do schedule(runtime) + do iCell=1, nCellsAll + ! boolean mask for now, could generalize to tappered sponge layer + ! only modify layer thicknesses on tidal boundary + do k = minLevelCell(iCell), maxLevelCell(iCell) + sshSubcycleCurWithTides(iCell) = sshSubcycleCur(iCell) + & + (tidalInputMask(iCell)*(tidalBCValue(iCell) - sshSubcycleCur(iCell))) + sshSubcycleNewWithTides(iCell) = sshSubcycleNew(iCell) + & + (tidalInputMask(iCell)*(tidalBCValue(iCell) - sshSubcycleNew(iCell))) + end do + end do + !$omp end do + !$omp end parallel + call mpas_pool_get_array(forcingPool, & + 'sshSubcycleCurWithTides', & + sshSubcycleCur) + call mpas_pool_get_array(forcingPool, & + 'sshSubcycleNewWithTides', & + sshSubcycleNew) + end if + ! Need to initialize btr_vel_temp over one more halo ! than we are computing over nEdges = nEdgesHalo(min(edgeHaloComputeCounter, & @@ -1437,12 +1506,13 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ !$omp parallel !$omp do schedule(runtime) - do iCell = 1, nCells + do iCell = 1, nCellsAll ! SSH is linear combination of SSHold and SSHnew sshTemp(iCell) = (1-config_btr_gam2_SSHWt1)* & sshSubcycleCur(iCell) + & config_btr_gam2_SSHWt1 * & sshSubcycleNew(iCell) + deltaSSH(iCell) = sshTemp(iCell) - sshNew(iCell) end do ! cell loop for ssh !$omp end do !$omp end parallel @@ -1458,8 +1528,8 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ !$omp end do !$omp end parallel - call ocn_diagnostic_solve_wctEdge(btrvel_temp, sshTemp, & - bottomDepthEdge, wctEdge) + call ocn_diagnostic_solve_wctEdge(btrvel_temp, deltaSSH, & + baroclinicThickness, wctEdge) !$omp parallel !$omp do schedule(runtime) & @@ -1489,18 +1559,19 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ !$omp parallel !$omp do schedule(runtime) - do iCell = 1, nCells + do iCell = 1, nCellsAll ! SSH is linear combination of SSHold and SSHnew sshTemp(iCell) = (1-config_btr_gam2_SSHWt1)* & sshSubcycleCur(iCell) + & config_btr_gam2_SSHWt1 * & sshSubcycleNew(iCell) + deltaSSH(iCell) = sshTemp(iCell) - sshNew(iCell) end do ! cell loop for ssh !$omp end do !$omp end parallel - call ocn_diagnostic_solve_wctEdge(btrvel_temp, sshTemp, & - bottomDepthEdge, wctEdge) + call ocn_diagnostic_solve_wctEdge(btrvel_temp, deltaSSH, & + baroclinicThickness, wctEdge) ! compute barotropic thickness flux on edges !$omp parallel @@ -1542,7 +1613,7 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ !$omp parallel !$omp do schedule(runtime) - do iEdge = 1, nEdgesAll + do iEdge = 1, nEdgesOwned normalBarotropicVelocityNew(iEdge) = & normalBarotropicVelocityNew(iEdge) + & normalBarotropicVelocitySubcycleNew(iEdge) @@ -1559,6 +1630,7 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ deallocate(btrvel_temp) deallocate(wctEdge) + deallocate(deltaSSH) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! END Barotropic subcycle loop @@ -1570,7 +1642,7 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ !$omp parallel !$omp do schedule(runtime) - do iEdge = 1, nEdgesAll + do iEdge = 1, nEdgesOwned barotropicThicknessFlux(iEdge) = & barotropicThicknessFlux(iEdge) & /(nBtrSubcycles*config_btr_subcycle_loop_factor) @@ -2831,6 +2903,31 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ #endif end if + ! direct application of tidal boundary condition + if (config_use_tidal_forcing .and. trim(config_tidal_forcing_type) == 'direct') then + call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness) + call mpas_pool_get_array(forcingPool, 'tidalInputMask', tidalInputMask) + call mpas_pool_get_array(forcingPool, 'tidalBCValue', tidalBCValue) + do iCell=1, nCells + ! artificially assumes boolean mask for now, could generalize to tappered sponge layer + 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) + !(1.0_RKIND - tidalInputMask(iCell))*layerThicknessNew(k, iCell) ! generalized tappered assumption code + end do + end if + end do + end if + call ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, verticalMeshPool, & scratchPool, tracersPool, 2) @@ -2953,7 +3050,7 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ call mpas_dmpar_exch_halo_field(effectiveDensityField) call mpas_timer_stop("se effective density halo") end if - deallocate(bottomDepthEdge) + deallocate(baroclinicThickness) ! get current time & setup num TS iteration nowTime = mpas_get_clock_time(domain%clock, MPAS_NOW, err)