diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F index 7f1ac1f921e8..7b39b3c29965 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F @@ -173,6 +173,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ edgeHaloComputeCounter, &! halo counters to reduce cellHaloComputeCounter ! halo updates during cycling + integer, dimension(:), pointer :: nEdgesArray + real (kind=RKIND) :: & normalThicknessFluxSum, &! sum of thickness flux in column thicknessSum, &! sum of thicknesses in column @@ -191,7 +193,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ uTemp real (kind=RKIND), dimension(:), allocatable :: & - btrvel_temp + btrvel_temp, & + bottomDepthEdge ! State Array Pointers real (kind=RKIND), dimension(:), pointer :: & @@ -382,6 +385,10 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracersNew, 2) + call mpas_pool_get_dimension(meshPool,'nEdgesArray', nEdgesArray) + + allocate(bottomDepthEdge(nEdgesAll+1)) + if (config_transport_tests_flow_id > 0) then ! This is a transport test. Write advection velocity from prescribed ! flow field. @@ -768,6 +775,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ enddo barotropicForcing(iEdge) = splitFact* & normalThicknessFluxSum/thicknessSum/dt + bottomDepthEdge(iEdge) = thicknessSum & + - 0.5_RKIND*(sshNew(cell1) + sshNew(cell2)) do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) ! These two steps are together here: @@ -786,6 +795,23 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ !$omp end do !$omp end parallel + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k, thicknessSum) + do iEdge = nEdgesOwned+1, nEdgesArray(4) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + thicknessSum = layerThickEdgeFlux(minLevelEdgeBot(iEdge),iEdge) + do k = minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge) + thicknessSum = thicknessSum + & + layerThickEdgeFlux(k,iEdge) + enddo + bottomDepthEdge(iEdge) = thicknessSum & + - 0.5_RKIND*(sshNew(cell1) + sshNew(cell2)) + enddo ! iEdge + !$omp end do + !$omp end parallel + deallocate(uTemp) call mpas_timer_start("se halo normalBaroclinicVelocity") @@ -1089,6 +1115,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ sshTend(iCell) = 0.0_RKIND do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) + if (maxLevelEdgeTop(iEdge).eq.0) cycle cell1 = cellsOnEdge(1, iEdge) cell2 = cellsOnEdge(2, iEdge) @@ -1096,20 +1123,9 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ sshEdge = 0.5_RKIND*(sshSubcycleCur(cell1) + & sshSubcycleCur(cell2)) - ! method 0: orig, works only without pbc: - !thicknessSum = sshEdge + & - ! refBottomDepthTopOfCell(maxLevelEdgeTop(iEdge)+1) - - ! method 1, matches method 0 without pbcs, works with pbcs. - thicknessSum = sshEdge + min(bottomDepth(cell1), & - bottomDepth(cell2)) - - ! method 2: may be better than method 1. - ! Take average of full thickness at two - ! neighboring cells. - !thicknessSum = sshEdge + 0.5* & - ! (bottomDepth(cell1) + bottomDepth(cell2)) - + ! Compute barotropic thickness at the edge. Note this + ! matches the sum of baroclinic thicknesses at this edge. + thicknessSum = sshEdge + bottomDepthEdge(iEdge) flux = ((1.0-config_btr_gam1_velWt1)* & normalBarotropicVelocitySubcycleCur(iEdge) & @@ -1149,16 +1165,14 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ !$omp do schedule(runtime) & !$omp private(cell1, cell2, sshEdge,thicknessSum,flux) do iEdge = 1, nEdges + if (maxLevelEdgeTop(iEdge).eq.0) cycle cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) sshEdge = 0.5_RKIND*(sshSubcycleCur(cell1) + & sshSubcycleCur(cell2)) - ! method 1, matches method 0 without pbcs, - ! works with pbcs. - thicknessSum = sshEdge + min(bottomDepth(cell1), & - bottomDepth(cell2)) + thicknessSum = sshEdge + bottomDepthEdge(iEdge) flux = ((1.0-config_btr_gam1_velWt1)* & normalBarotropicVelocitySubcycleCur(iEdge) & @@ -1370,6 +1384,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ sshTend(iCell) = 0.0_RKIND do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) + if (maxLevelEdgeTop(iEdge).eq.0) cycle cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) @@ -1387,21 +1402,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ sshEdge = 0.5_RKIND*(sshCell1 + sshCell2) - ! method 0: orig, works only without pbc: - !thicknessSum = sshEdge + & - !refBottomDepthTopOfCell(maxLevelEdgeTop(iEdge)+1) - - ! method 1, matches method 0 without pbcs, - ! but works with pbcs. - thicknessSum = sshEdge + & - min(bottomDepth(cell1), & - bottomDepth(cell2)) - - ! method 2: may be better than method 1. - ! take average of full thickness at two - ! neighboring cells - !thicknessSum = sshEdge + 0.5* & - ! (bottomDepth(cell1) + bottomDepth (cell2)) + thicknessSum = sshEdge + bottomDepthEdge(iEdge) flux = ((1.0-config_btr_gam3_velWt2)* & normalBarotropicVelocitySubcycleCur(iEdge) & @@ -1441,21 +1442,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ sshSubcycleNew(cell2) sshEdge = 0.5_RKIND * (sshCell1 + sshCell2) - ! method 0: orig, works only without pbc: - !thicknessSum = sshEdge + & - !refBottomDepthTopOfCell(maxLevelEdgeTop(iEdge)+1) - - ! method 1, matches method 0 without pbcs, - ! but works with pbcs. - thicknessSum = sshEdge + min(bottomDepth(cell1),& - bottomDepth(cell2)) - - ! method 2, better, I think. - ! take average of full thickness at two - ! neighboring cells - !thicknessSum = sshEdge + 0.5* & - ! (bottomDepth(cell1) + & - ! bottomDepth(cell2) ) + thicknessSum = sshEdge + bottomDepthEdge(iEdge) flux = ((1.0-config_btr_gam3_velWt2) * & normalBarotropicVelocitySubcycleCur(iEdge) & @@ -2695,6 +2682,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ call mpas_dmpar_exch_halo_field(effectiveDensityField) call mpas_timer_stop("se effective density halo") end if + deallocate(bottomDepthEdge) call mpas_timer_stop('se fini') call mpas_timer_stop("se timestep")