From 6bab166a9ef0fa527a258bfe3238c0cf9f504653 Mon Sep 17 00:00:00 2001 From: Damrongsak Wirasaet Date: Tue, 8 Aug 2023 13:54:46 -0700 Subject: [PATCH 01/16] Modify diagontic_solve. Compute ssh before layer thickness edge --- .../src/shared/mpas_ocn_diagnostics.F | 24 +++++++++++++------ 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F index b0401b420e00..dca43e1e4991 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F @@ -235,6 +235,16 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, scratchPoo ! endif #endif + ! DW: Move here + ! + ! Z-coordinates + ! This section must be placed in the code before computing the density. + ! + ! inputs : layerThickness + ! outputs : zMid, zTop, ssh) + call ocn_diagnostic_solve_z_coordinates(layerThickness, zMid, zTop, ssh) + + ! inputs: layerThickness, normalVelocity ! output: layerThickEdgeMean, layerThickEdgeDrag, layerThickEdgeFlux call ocn_diagnostic_solve_layerThicknessEdge(normalVelocity, & @@ -334,13 +344,13 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, scratchPoo call ocn_diagnostic_solve_surface_pressure(forcingPool, atmosphericPressure, & seaIcePressure, surfacePressure) - ! - ! Z-coordinates - ! This section must be placed in the code before computing the density. - ! - ! inputs : layerThickness - ! outputs : zMid, zTop, ssh) - call ocn_diagnostic_solve_z_coordinates(layerThickness, zMid, zTop, ssh) +! ! +! ! Z-coordinates +! ! This section must be placed in the code before computing the density. +! ! +! ! inputs : layerThickness +! ! outputs : zMid, zTop, ssh) +! call ocn_diagnostic_solve_z_coordinates(layerThickness, zMid, zTop, ssh) ! ! equation of state From 315e3f443d7b28f3f967bba5aa148f1639bb6dba Mon Sep 17 00:00:00 2001 From: Damrongsak Wirasaet Date: Tue, 15 Aug 2023 14:26:38 -0700 Subject: [PATCH 02/16] 1. compute ssh before edge thickness 2. use the min bathymetry pixels for the subgrid edge bathymetry --- .../mode_init/mpas_ocn_init_parabolic_bowl.F | 281 ++++++++++++++---- .../src/shared/mpas_ocn_diagnostics.F | 27 +- 2 files changed, 242 insertions(+), 66 deletions(-) diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F index 94f6c9b31c90..09b859f26196 100644 --- a/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F +++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F @@ -177,6 +177,14 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ integer :: slice, nSlice integer :: i,j real (kind=RKIND) :: deltaZ +! DW + integer:: nsubgridCellEdge, iEdgeSegment + real (kind=RKIND), dimension(:,:), allocatable :: cellEdgeBathymetryValues + real (kind=RKIND), dimension(:), allocatable:: dsEdge + integer:: jj + real (kind=RKIND), dimension(:), allocatable:: xSubgridCell, ySubgridCell + + iErr = 0 @@ -374,7 +382,7 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ do iCell = 1, nCellsSolve ! Set up vertical grid maxLevelCell(iCell) = nVertLevels ; ! sigma coordinates - end do + end do ! Set velocity @@ -475,10 +483,12 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ subgridSshCellTableRange(1,iCell) = -maxval(subgridBathymetryValues(1:nSubgridCell)) + config_drying_min_cell_height + eps subgridSshCellTableRange(2,iCell) = -minval(subgridBathymetryValues(1:nSubgridCell)) - + ! Evaluate subgrid bathymetry !--------------------------------------------------------------- - bottomDepth(iCell) = sum(subgridBathymetryValues(1:nSubgridCell))/real(nSubgridCell,RKIND) + ! DW + bottomDepth(iCell) = sum(subgridBathymetryValues(1:nSubgridCell)*subgridAreas(1:nSubgridCell))/sum(subgridAreas(1:nSubgridCell)) +!! real(nSubgridCell,RKIND) subgridCellBathymetryMin(iCell) = maxval(subgridBathymetryValues(1:nSubgridCell)) !subgridCellBathymetryMin(iCell) = bottomDepth(iCell) @@ -486,22 +496,25 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ !--------------------------------------------------------------- call ocn_init_vertical_integration(iCell,subgridSshCellTableRange, nSubgridCell, subgridBathymetryValues, subgridAreas, subgridWetVolumeCellTable) - !print*, iCell, areaCell(iCell)-sum(subgridAreas(1:nSubgridCell)) - - ! Evaluate wet layerThickness average - !--------------------------------------------------------------- - call ocn_init_wet_average(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, layerThickness(1,iCell)) - - if (iCell == 3030) then - print*, subgridSshCellTableRange(1,iCell) - print*, subgridSshCellTableRange(2,iCell) - print*, nSubgridTableLevels - deltaZ = (subgridSshCellTableRange(2,iCell) - subgridSshCellTableRange(1,iCell))/(nSubgridTableLevels-1) - print*, deltaZ - do i = 1,nSubgridTableLevels - print*, (real(i,RKIND)-1.0_RKIND)*deltaZ + subgridSshCellTableRange(1,iCell), subgridWetVolumeCellTable(i,iCell) - enddo - endif + ! print*, iCell, areaCell(iCell)-sum(subgridAreas(1:nSubgridCell)) + + call ocn_init_wet_average_ssh(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, ssh(iCell))!{{{ + +!!$ ! Evaluate wet layerThickness average +!!$ !--------------------------------------------------------------- +!!$ ! call ocn_init_wet_average(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, layerThickness(1,iCell)) +!!$ call ocn_init_grid_average(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, layerThickness(1,iCell)) +!!$ +!!$ if (iCell == 3030) then +!!$ print*, subgridSshCellTableRange(1,iCell) +!!$ print*, subgridSshCellTableRange(2,iCell) +!!$ print*, nSubgridTableLevels +!!$ deltaZ = (subgridSshCellTableRange(2,iCell) - subgridSshCellTableRange(1,iCell))/(nSubgridTableLevels-1) +!!$ print*, deltaZ +!!$ do i = 1,nSubgridTableLevels +!!$ print*, (real(i,RKIND)-1.0_RKIND)*deltaZ + subgridSshCellTableRange(1,iCell), subgridWetVolumeCellTable(i,iCell) +!!$ enddo +!!$ endif enddo @@ -512,6 +525,16 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! print*, "Begin edges", nEdgesSolve + nSubgridCellEdge=config_parabolic_bowl_subgrid_level + allocate( cellEdgeBathymetryValues(3,nSubgridCellEdge) ) + allocate( dsEdge(nSubgridCellEdge) ) + cellEdgeBathymetryValues = -99999 ; + + allocate(xSubgridCell(maxEdges*nSubgridTriPerSlice)) + allocate(ySubgridCell(maxEdges*nSubgridTriPErSlice)) + +! open( unit = 22, file = 'edgediag.dat', action = 'write' ) ; + do iEdge = 1,nEdgesSolve ! Evaluate subgrid bathymetry at centers of sub-triangles for edge slices ! (all subdivided triangles for each edge slice are gathered into @@ -540,29 +563,60 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ x(3) = xCell(cellsOnEdge(slice,iEdge)) y(3) = yCell(cellsOnEdge(slice,iEdge)) - call ocn_init_evaluate_subgrid_data(x, y, nSubgridTriPerSlice, nSubgridEdge, rSubgridPoints, sSubgridPoints, subgridBathymetryValues, subgridAreas) - + + call ocn_init_evaluate_subgrid_data(x, y, nSubgridTriPerSlice, nSubgridEdge, rSubgridPoints, sSubgridPoints, & + subgridBathymetryValues, subgridAreas, xSubgridCell=xSubgridCell, ySubgridCell=ySubgridCell ) ; enddo + ! DW: Use the higher values of the pair of subcells along the + ! cell edge + cellEdgeBathymetryValues(1,:) = subgridBathymetryValues(1:2*nSubgridCellEdge - 1:2) + if ( nslice > 1 ) then + cellEdgeBathymetryValues(2,:) = & + subgridBathymetryValues(nSubgridTriPerSlice+2*nsubgridCellEdge - 1:nSubgridTriPerSlice+1:-2) + else + cellEdgeBathymetryValues(2,:) = cellEdgeBathymetryValues(1,:) + endif + + do iEdgeSegment = 1, nSubgridCellEdge + cellEdgeBathymetryValues(3,iEdgeSegment) = minval(cellEdgeBathymetryValues(1:2,iEdgeSegment) ) + end do + dsEdge(:) = sqrt( (x(2) - x(1))*(x(2) - x(1)) + (y(2) - y(1))*(y(2)- y(1)) )/nSubgridCellEdge + +! write(22,'(4I8,4F17.7)') iEdge, CellsOnEdge(1:2,iEdge), nSubgridTriPerSlice, x(1:2), y(1:2) +! do jj = 1, nSubgridEdge +! write(22,'(3F17.7)') xSubgridCell(jj), ySubgridCell(jj), subgridBathymetryValues(jj) +! end do +! write(22,*) ' ' ; +! +! do iEdgeSegment = 1, nSubgridCellEdge +! write(22,'(E17.7)') cellEdgeBathymetryValues(3,iEdgeSegment) ; +! end do +! write(22,*) '' ; + ! Evaluate bounds of look-up table range !--------------------------------------------------------------- - subgridSshEdgeTableRange(1,iEdge) = -maxval(subgridBathymetryValues(1:nSubgridEdge)) + config_drying_min_cell_height + eps - subgridSshEdgeTableRange(2,iEdge) = -minval(subgridBathymetryValues(1:nSubgridEdge)) +! subgridSshEdgeTableRange(1,iEdge) = -maxval(subgridBathymetryValues(1:nSubgridEdge)) + config_drying_min_cell_height + eps +! subgridSshEdgeTableRange(2,iEdge) = -minval(subgridBathymetryValues(1:nSubgridEdge)) + subgridSshEdgeTableRange(1,iEdge) = -maxval(CellEdgeBathymetryValues(3,1:nSubgridCellEdge)) + config_drying_min_cell_height + eps + subgridSshEdgeTableRange(2,iEdge) = -minval(CellEdgeBathymetryValues(3,1:nSubgridCellEdge)) + ! Evaluate subgrid bathymetry !--------------------------------------------------------------- - subgridEdgeBathymetryMean(iEdge) = sum(subgridBathymetryValues(1:nSubgridEdge))/real(nSubgridEdge,RKIND) - subgridEdgeBathymetryMin(iEdge) = minval(subgridBathymetryValues(1:nSubgridEdge)) +! subgridEdgeBathymetryMean(iEdge) = sum(subgridBathymetryValues(1:nSubgridEdge))/real(nSubgridEdge,RKIND) +! subgridEdgeBathymetryMin(iEdge) = minval(subgridBathymetryValues(1:nSubgridEdge)) + subgridEdgeBathymetryMean(iEdge) = sum(cellEdgeBathymetryValues(3,1:nSubgridCellEdge))/real(nSubgridCellEdge,RKIND) + subgridEdgeBathymetryMin(iEdge) = maxval(cellEdgeBathymetryValues(3,1:nSubgridCellEdge)) ! Vertical integration of wet fraction !--------------------------------------------------------------- - - call ocn_init_vertical_integration(iEdge,subgridSshEdgeTableRange, nSubgridEdge, subgridBathymetryValues, subgridAreas, subgridWetVolumeEdgeTable) - - enddo - - +! call ocn_init_vertical_integration(iEdge,subgridSshEdgeTableRange, nSubgridEdge, subgridBathymetryValues, subgridAreas, subgridWetVolumeEdgeTable) + call ocn_init_vertical_integration( iEdge, subgridSshEdgeTableRange, & + nSubgridCellEdge, cellEdgeBathymetryValues(3,:), dsEdge, subgridWetVolumeEdgeTable ) + end do +! close(22) ; !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Vertex @@ -617,8 +671,11 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ ! Evaluate subgrid bathymetry !--------------------------------------------------------------- - subgridVertexBathymetryMean(iVertex) = sum(subgridBathymetryValues(1:nSubgridVertex))/real(nSubgridVertex,RKIND) - subgridVertexBathymetryMin(iVertex) = minval(subgridBathymetryValues(1:nSubgridVertex)) + subgridVertexBathymetryMean(iVertex) = sum(subgridBathymetryValues(1:nSubgridVertex)*subgridAreas(1:nSubgridVertex))/sum(subgridAreas(1:nSubgridVertex)) + ! real(nSubgridVertex,RKIND) +! subgridVertexBathymetryMin(iVertex) = minval(subgridBathymetryValues(1:nSubgridVertex)) + subgridVertexBathymetryMin(iVertex) = maxval(subgridBathymetryValues(1:nSubgridVertex)) + ! Vertical integration of wet fraction !--------------------------------------------------------------- @@ -656,12 +713,12 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ ! if (config_use_subgrid_wetting_drying) then - call ocn_subgrid_ssh_lookup(layerThickness(1,iCell),& - subgridWetVolumeCellTable(:,iCell),& - subgridSshCellTableRange(:,iCell),& - bottomDepth(iCell),& - subgridCellBathymetryMin(iCell),& - ssh(iCell)) +!! call ocn_subgrid_ssh_lookup(layerThickness(1,iCell),& +!! subgridWetVolumeCellTable(:,iCell),& +!! subgridSshCellTableRange(:,iCell),& +!! bottomDepth(iCell),& +!! subgridCellBathymetryMin(iCell),& +!! ssh(iCell)) !call ocn_init_parabolic_bowl_ssh(xCell(iCell),yCell(iCell),bottomDepth(iCell),ssh(iCell)) !ssh(iCell) = - bottomDepth(iCell) + & @@ -701,13 +758,24 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ block_ptr => block_ptr % next end do - do iCell = 1,nCellsSolve - call ocn_subgrid_layer_thickness_lookup(ssh(iCell), & + ! do iCell = 1,nCellsSolve + ! call ocn_subgrid_layer_thickness_lookup(ssh(iCell), & + ! subgridWetVolumeCellTable(:,iCell), & + ! subgridSshCellTableRange(:,iCell),& + ! bottomDepth(iCell),& + ! subgridLayerThicknessDebug(iCell)) + ! enddo + + ! DW + do iCell = 1,nCellsSolve + call ocn_subgrid_layer_thickness_lookup(ssh(iCell), & subgridWetVolumeCellTable(:,iCell), & subgridSshCellTableRange(:,iCell),& bottomDepth(iCell),& - subgridLayerThicknessDebug(iCell)) - enddo + LayerThickness(1,iCell)) + enddo + ! + !do iCell = 1,nCellsSolve ! call ocn_subgrid_ssh_lookup(layerThickness(1,iCell),& ! subgridWetVolumeCellTable(:,iCell),& @@ -834,7 +902,7 @@ end subroutine ocn_init_define_subgrid_points!}}} ! !----------------------------------------------------------------------- - subroutine ocn_init_evaluate_subgrid_data(xTri, yTri, nSubgridTri, nSubgridCV, rSubgridPoints, sSubgridPoints, subgridBathymetryValues, subgridAreas, subgridSshValues)!{{{ + subroutine ocn_init_evaluate_subgrid_data(xTri, yTri, nSubgridTri, nSubgridCV, rSubgridPoints, sSubgridPoints, subgridBathymetryValues, subgridAreas, subgridSshValues, xSubgridCell, ySubgridCell )!{{{ implicit none real (kind=RKIND), intent(inout) :: xTri(3), yTri(3) @@ -844,12 +912,16 @@ subroutine ocn_init_evaluate_subgrid_data(xTri, yTri, nSubgridTri, nSubgridCV, r real (kind=RKIND), dimension(:), intent(inout) :: subgridBathymetryValues real (kind=RKIND), dimension(:), intent(inout) :: subgridAreas real (kind=RKIND), dimension(:), intent(inout),optional :: subgridSshValues + real (kind=RKIND), dimension(:), intent(inout), optional :: xSubgridCell, ySubgridCell + + real (kind=RKIND) :: rCenter, sCenter real (kind=RKIND) :: xCenter, yCenter real (kind=RKIND) :: xSubgridPoints(3), ySubgridPoints(3) real (kind=RKIND) :: x, y real (kind=RKIND) :: area + integer :: iPt,i !-------------------------------------------------------------------- @@ -858,14 +930,21 @@ subroutine ocn_init_evaluate_subgrid_data(xTri, yTri, nSubgridTri, nSubgridCV, r ! (ensure counter-clockwise numering) call ocn_init_tri_area(xTri(:), yTri(:), area) if (area < 0.0_RKIND) then - x = xTri(2) - y = yTri(2) - xTri(2) = xTri(3) - yTri(2) = yTri(3) - xTri(3) = x - yTri(3) = y - !print*, area - area = abs(area) +! x = xTri(2) +! y = yTri(2) +! xTri(2) = xTri(3) +! yTri(2) = yTri(3) +! xTri(3) = x +! yTri(3) = y +! !print*, area +! area = abs(area) + x = xTri(1) + y = yTri(1) + xTri(1) = xTri(2) + yTri(1) = yTri(2) + xTri(2) = x ; + yTri(2) = y ; + area = abs(area) ; endif !print("(3(f10.3,1x))"), (xTri(i)/1000, i=1,3) @@ -900,7 +979,13 @@ subroutine ocn_init_evaluate_subgrid_data(xTri, yTri, nSubgridTri, nSubgridCV, r endif !print("(I5,I5,1x,5(f10.3,1x))"), nSubgridCV, iPt, x/1000, y/1000, subgridBathymetryValues(nSubgridCV), subgridAreas(nSubgridCV) + if ( present(xSubgridCell) ) then + xSubgridCell(nSubgridCV) = xCenter ; + endif + if ( present(ySubgridCell) ) then + ySubgridCell(nSubgridCV) = yCenter ; + end if enddo !-------------------------------------------------------------------- @@ -1097,7 +1182,7 @@ end subroutine ocn_init_vertical_integration!}}} !> ! !----------------------------------------------------------------------- - + ! subroutine ocn_init_wet_average(nSubgridCV, subgridBathymetryValues, subgridSshValues, subgridAreas, subgridWetAverage)!{{{ implicit none @@ -1113,9 +1198,58 @@ subroutine ocn_init_wet_average(nSubgridCV, subgridBathymetryValues, subgridSshV !-------------------------------------------------------------------- - averageDepth = sum(subgridBathymetryValues(1:nSubgridCV))/real(nSubgridCV,RKIND) + subgridWetAverage = 0.0_RKIND + wetArea = 0.0_RKIND + do tri = 1,nSubgridCV + + layerThicknessValue = subgridBathymetryValues(tri) + subgridSshValues(tri) + !if (layerThicknessValue < config_drying_min_cell_height + eps) then + ! layerThicknessValue = config_drying_min_cell_height + eps + !endif + !subgridWetAverage = subgridWetAverage + layerThicknessValue*subgridAreas(tri) + !wetArea = wetArea + subgridAreas(tri) + if ( layerThicknessValue > config_drying_min_cell_height ) then + subgridWetAverage = subgridWetAverage + layerThicknessValue*subgridAreas(tri) + WetArea = WetArea + subgridAreas(tri) ; + endif + enddo + + if ( WetArea > 0.D0 ) then + subgridWetAverage = subgridWetAverage/wetArea ; + else + subgridWetAverage = config_drying_min_cell_height + eps ; + endif + + if ( subgridWetAverage < config_drying_min_cell_height + eps ) then + subgridWetAverage = config_drying_min_cell_height + eps ; + endif + + !-------------------------------------------------------------------- + end subroutine ocn_init_wet_average!}}} + + + ! DW + subroutine ocn_init_grid_average(nSubgridCV, subgridBathymetryValues, subgridSshValues, subgridAreas, subgridWetAverage)!{{{ + + implicit none + + integer, intent(in) :: nSubgridCV + real (kind=RKIND), dimension(:), intent(in) :: subgridBathymetryValues, subgridAreas + real (kind=RKIND), dimension(:), intent(in) :: subgridSshValues + real (kind=RKIND), intent(inout) :: subgridWetAverage + + real (kind=RKIND) :: deltaZ, ssh, pVal + real (kind=RKIND) :: averageDepth, wetArea, layerThicknessValue + integer :: lev, tri + + + !------------------------------------------------------------------- + ! averageDepth = sum(subgridBathymetryValues(1:nSubgridCV))/real(nSubgridCV,RKIND) + averageDepth = sum(subgridBathymetryValues(1:nSubgridCV)*subgridAreas(1:nSubgridCV))/sum(subgridAreas(1:nSubgridCV)) + ! real(nSubgridCV,RKIND) + subgridWetAverage = 0.0_RKIND wetArea = 0.0_RKIND do tri = 1,nSubgridCV @@ -1130,10 +1264,45 @@ subroutine ocn_init_wet_average(nSubgridCV, subgridBathymetryValues, subgridSshV subgridWetAverage = subgridWetAverage/wetArea - !-------------------------------------------------------------------- + end subroutine ocn_init_grid_average!}}} - end subroutine ocn_init_wet_average!}}} + + ! DW + subroutine ocn_init_wet_average_ssh( nSubgridCV, subgridBathymetryValues, subgridSshValues, subgridAreas, sshWetAverage)!{{{ + implicit none + + integer, intent(in) :: nSubgridCV + real (kind=RKIND), dimension(:), intent(in) :: subgridBathymetryValues, subgridAreas + real (kind=RKIND), dimension(:), intent(in) :: subgridSshValues + real (kind=RKIND), intent(inout) :: sshWetAverage + + real (kind=RKIND) :: deltaZ, ssh, pVal + real (kind=RKIND) :: averageDepth, wetArea, layerThicknessValue + integer :: lev, tri + + + sshWetAverage = 0.0_RKIND + wetArea = 0.0_RKIND + + do tri = 1, nsubgridCV + layerThicknessValue = subgridBathymetryValues(tri) + subgridSshValues(tri) + + if (layerThicknessValue > config_drying_min_cell_height + eps ) then + sshWetAverage = sshWetAverage + subgridSshValues(tri)*subgridAreas(tri) + wetArea = wetArea + subgridAreas(tri) + endif + end do + + if ( WetArea > 0.0_RKIND ) then + sshWetAverage = sshWetAverage/WetArea ; + else + sshWetAverage = -maxval(subgridBathymetryValues(1:nsubgridCV)) & + + config_drying_min_cell_height ; + endif + + return ; + end subroutine ocn_init_wet_average_ssh !*********************************************************************** ! diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F index dca43e1e4991..9e93f7c56cfb 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F @@ -845,19 +845,26 @@ subroutine ocn_diagnostic_solve_layerThicknessEdge(normalVelocity, & subgridEdgeBathymetryMean(iEdge), & layerThickEdgeFlux(1,iEdge)) else - call ocn_subgrid_layer_thickness_lookup(ssh(cell1), & - subgridWetVolumeEdgeTable(:,iEdge), & - subgridSshEdgeTableRange(:,iEdge), & - subgridEdgeBathymetryMean(iEdge), & - layerThicknessCell1) - call ocn_subgrid_layer_thickness_lookup(ssh(cell2), & +!!$ call ocn_subgrid_layer_thickness_lookup(ssh(cell1), & +!!$ subgridWetVolumeEdgeTable(:,iEdge), & +!!$ subgridSshEdgeTableRange(:,iEdge), & +!!$ subgridEdgeBathymetryMean(iEdge), & +!!$ layerThicknessCell1) +!!$ call ocn_subgrid_layer_thickness_lookup(ssh(cell2), & +!!$ subgridWetVolumeEdgeTable(:,iEdge), & +!!$ subgridSshEdgeTableRange(:,iEdge), & +!!$ subgridEdgeBathymetryMean(iEdge), & +!!$ layerThicknessCell2) +!!$ layerThickEdgeFlux(1,iEdge) = & +!!$ max(layerThicknessCell1, & +!!$ layerThicknessCell2) + + sshMean = (ssh(cell1) + ssh(cell2))/2.0_RKIND ; + call ocn_subgrid_layer_thickness_lookup(sshMean, & subgridWetVolumeEdgeTable(:,iEdge), & subgridSshEdgeTableRange(:,iEdge), & subgridEdgeBathymetryMean(iEdge), & - layerThicknessCell2) - layerThickEdgeFlux(1,iEdge) = & - max(layerThicknessCell1, & - layerThicknessCell2) + layerThickEdgeFlux(1,iedge)) end if end do #ifndef MPAS_OPENACC From 1d453e1ff8472ee32c94cf00c8197a2c2d8e7f73 Mon Sep 17 00:00:00 2001 From: Damrongsak Wirasaet Date: Fri, 18 Aug 2023 09:46:15 -0700 Subject: [PATCH 03/16] Set subgrid edge pixels to maximum bottom elevation pixels --- .../src/mode_init/Registry_parabolic_bowl.xml | 6 ++- .../mode_init/mpas_ocn_init_parabolic_bowl.F | 51 ++++++++++--------- 2 files changed, 31 insertions(+), 26 deletions(-) diff --git a/components/mpas-ocean/src/mode_init/Registry_parabolic_bowl.xml b/components/mpas-ocean/src/mode_init/Registry_parabolic_bowl.xml index 0541dd0fb408..37df0e3abd4f 100644 --- a/components/mpas-ocean/src/mode_init/Registry_parabolic_bowl.xml +++ b/components/mpas-ocean/src/mode_init/Registry_parabolic_bowl.xml @@ -31,7 +31,11 @@ + /> + diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F index 7b0449f06047..3ff67bd38051 100644 --- a/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F +++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F @@ -552,42 +552,43 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ call ocn_init_wet_average(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, layerThicknessEdgeAverage, & subgridUValues, subgridVValues, uVelocityAverage(iEdge), vVelocityAverage(iEdge)) + if ( .NOT. config_parabolic_bowl_subgrid_edge_bathymetry_max_pixel ) then + ! Evaluate bounds of look-up table range + !--------------------------------------------------------------- - ! Evaluate bounds of look-up table range - !--------------------------------------------------------------- + subgridSshEdgeTableRange(1,iEdge) = -maxval(subgridBathymetryValues(1:nSubgridEdge)) + config_drying_min_cell_height + eps + subgridSshEdgeTableRange(2,iEdge) = -minval(subgridBathymetryValues(1:nSubgridEdge)) -! subgridSshEdgeTableRange(1,iEdge) = -maxval(subgridBathymetryValues(1:nSubgridEdge)) + config_drying_min_cell_height + eps -! subgridSshEdgeTableRange(2,iEdge) = -minval(subgridBathymetryValues(1:nSubgridEdge)) + ! Evaluate subgrid bathymetry + !--------------------------------------------------------------- - ! Evaluate subgrid bathymetry - !--------------------------------------------------------------- -! -! subgridEdgeBathymetryMean(iEdge) = sum(subgridBathymetryValues(1:nSubgridEdge))/real(nSubgridEdge,RKIND) -! subgridEdgeBathymetryMin(iEdge) = maxval(subgridBathymetryValues(1:nSubgridEdge)) - - ! Vertical integration of wet fraction - !--------------------------------------------------------------- + subgridEdgeBathymetryMean(iEdge) = sum(subgridBathymetryValues(1:nSubgridEdge))/real(nSubgridEdge,RKIND) + subgridEdgeBathymetryMin(iEdge) = maxval(subgridBathymetryValues(1:nSubgridEdge)) -! call ocn_init_vertical_integration(iEdge,subgridSshEdgeTableRange, nSubgridEdge, subgridBathymetryValues, subgridAreas, & -! subgridWetVolumeEdgeTable, subgridWetFractionEdgeTable) + ! Vertical integration of wet fraction + !--------------------------------------------------------------- - ! Evaluate bounds of look-up table range - !--------------------------------------------------------------- + call ocn_init_vertical_integration(iEdge,subgridSshEdgeTableRange, nSubgridEdge, subgridBathymetryValues, subgridAreas, & + subgridWetVolumeEdgeTable, subgridWetFractionEdgeTable) + else + ! Evaluate bounds of look-up table range + !--------------------------------------------------------------- - subgridSshEdgeTableRange(1,iEdge) = -maxval(CellEdgeBathymetryValues(3,1:nSubgridCellEdge)) + config_drying_min_cell_height + eps - subgridSshEdgeTableRange(2,iEdge) = -minval(CellEdgeBathymetryValues(3,1:nSubgridCellEdge)) + subgridSshEdgeTableRange(1,iEdge) = -maxval(CellEdgeBathymetryValues(3,1:nSubgridCellEdge)) + config_drying_min_cell_height + eps + subgridSshEdgeTableRange(2,iEdge) = -minval(CellEdgeBathymetryValues(3,1:nSubgridCellEdge)) - ! Evaluate bounds of look-up table range - !--------------------------------------------------------------- + ! Evaluate bounds of look-up table range + !--------------------------------------------------------------- - subgridEdgeBathymetryMean(iEdge) = sum(cellEdgeBathymetryValues(3,1:nSubgridCellEdge))/real(nSubgridCellEdge,RKIND) - subgridEdgeBathymetryMin(iEdge) = maxval(cellEdgeBathymetryValues(3,1:nSubgridCellEdge)) + subgridEdgeBathymetryMean(iEdge) = sum(cellEdgeBathymetryValues(3,1:nSubgridCellEdge))/real(nSubgridCellEdge,RKIND) + subgridEdgeBathymetryMin(iEdge) = maxval(cellEdgeBathymetryValues(3,1:nSubgridCellEdge)) - ! Vertical integration of wet fraction - !--------------------------------------------------------------- - call ocn_init_vertical_integration( iEdge, subgridSshEdgeTableRange, & + ! Vertical integration of wet fraction + !--------------------------------------------------------------- + call ocn_init_vertical_integration( iEdge, subgridSshEdgeTableRange, & nSubgridCellEdge, cellEdgeBathymetryValues(3,:), dsEdge, subgridWetVolumeEdgeTable, subgridWetFractionEdgeTable ) + endif end do !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! From 3209aada84832875e39a7a9e12543b29c799bdf3 Mon Sep 17 00:00:00 2001 From: Damrongsak Wirasaet Date: Wed, 23 Aug 2023 14:40:28 -0700 Subject: [PATCH 04/16] Set the defaluft value of the edge_max_pixel registry to true --- components/mpas-ocean/src/mode_init/Registry_parabolic_bowl.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/mpas-ocean/src/mode_init/Registry_parabolic_bowl.xml b/components/mpas-ocean/src/mode_init/Registry_parabolic_bowl.xml index 37df0e3abd4f..04d03b5d55fd 100644 --- a/components/mpas-ocean/src/mode_init/Registry_parabolic_bowl.xml +++ b/components/mpas-ocean/src/mode_init/Registry_parabolic_bowl.xml @@ -32,7 +32,7 @@ description="Level of refinement to evaluate subgrid bathymety" possible_values="Any positive integer number" /> - From d6efffe516ecb534a70f6d83096ad32d8d8735b6 Mon Sep 17 00:00:00 2001 From: Damrongsak Wirasaet Date: Thu, 24 Aug 2023 13:02:29 -0700 Subject: [PATCH 05/16] Commit an option of using max subgrid bathymetry along the edge to be used for the edge look up --- .../mode_init/mpas_ocn_init_parabolic_bowl.F | 52 ++++++++++--------- 1 file changed, 27 insertions(+), 25 deletions(-) diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F index 3ff67bd38051..d1a69a95a1f9 100644 --- a/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F +++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F @@ -527,26 +527,9 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ subgridBathymetryValues, subgridAreas, subgridSshValues, subgridUValues, subgridVValues) -!!$ call ocn_init_evaluate_subgrid_data(x, y, nSubgridTriPerSlice, nSubgridEdge, rSubgridPoints, sSubgridPoints, & -!!$ subgridBathymetryValues, subgridAreas, xSubgridCell=xSubgridCell, ySubgridCell=ySubgridCell ) ; + enddo - ! DW: Use the higher values of the pair of subcells along the - ! cell edge - cellEdgeBathymetryValues(1,:) = subgridBathymetryValues(1:2*nSubgridCellEdge - 1:2) - if ( nslice > 1 ) then - cellEdgeBathymetryValues(2,:) = & - subgridBathymetryValues(nSubgridTriPerSlice+2*nsubgridCellEdge - 1:nSubgridTriPerSlice+1:-2) - else - cellEdgeBathymetryValues(2,:) = cellEdgeBathymetryValues(1,:) - endif - - do iEdgeSegment = 1, nSubgridCellEdge - cellEdgeBathymetryValues(3,iEdgeSegment) = minval(cellEdgeBathymetryValues(1:2,iEdgeSegment) ) - end do - dsEdge(:) = sqrt( (x(2) - x(1))*(x(2) - x(1)) + (y(2) - y(1))*(y(2)- y(1)) )/nSubgridCellEdge - - ! Evaluate velocity average !--------------------------------------------------------------- call ocn_init_wet_average(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, layerThicknessEdgeAverage, & @@ -572,6 +555,22 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ call ocn_init_vertical_integration(iEdge,subgridSshEdgeTableRange, nSubgridEdge, subgridBathymetryValues, subgridAreas, & subgridWetVolumeEdgeTable, subgridWetFractionEdgeTable) else + ! DW: Use the higher values of the pair of subcells along the + ! cell edge + cellEdgeBathymetryValues(1,:) = subgridBathymetryValues(1:2*nSubgridCellEdge - 1:2) + if ( nslice > 1 ) then + cellEdgeBathymetryValues(2,:) = & + subgridBathymetryValues(nSubgridTriPerSlice+2*nsubgridCellEdge - 1:nSubgridTriPerSlice+1:-2) + else + cellEdgeBathymetryValues(2,:) = cellEdgeBathymetryValues(1,:) + endif + + do iEdgeSegment = 1, nSubgridCellEdge + cellEdgeBathymetryValues(3,iEdgeSegment) = minval(cellEdgeBathymetryValues(1:2,iEdgeSegment) ) + end do + dsEdge(:) = sqrt( (x(2) - x(1))*(x(2) - x(1)) + (y(2) - y(1))*(y(2)- y(1)) )/nSubgridCellEdge + + ! Evaluate bounds of look-up table range !--------------------------------------------------------------- @@ -589,6 +588,7 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ call ocn_init_vertical_integration( iEdge, subgridSshEdgeTableRange, & nSubgridCellEdge, cellEdgeBathymetryValues(3,:), dsEdge, subgridWetVolumeEdgeTable, subgridWetFractionEdgeTable ) endif + end do !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -1178,6 +1178,7 @@ subroutine ocn_init_vertical_integration(iCV, subgridSshTableRange, nSubgridCV, end subroutine ocn_init_vertical_integration!}}} + !*********************************************************************** ! ! routine ocn_init_wet_average @@ -1190,7 +1191,7 @@ end subroutine ocn_init_vertical_integration!}}} ! !----------------------------------------------------------------------- - subroutine ocn_init_wet_average(nSubgridCV, subgridBathymetryValues,subgridSshValues, subgridAreas, subgridThicknessAverage, & + subroutine ocn_init_wet_average(nSubgridCV, subgridBathymetryValues, subgridSshValues, subgridAreas, subgridThicknessAverage, & subgridUValues, subgridVValues, subgridUAverage, subgridVAverage)!{{{ implicit none @@ -1210,10 +1211,9 @@ subroutine ocn_init_wet_average(nSubgridCV, subgridBathymetryValues,subgridSshVa integer :: lev, tri - !------------------------------------------------------------------- - ! averageDepth = sum(subgridBathymetryValues(1:nSubgridCV))/real(nSubgridCV,RKIND) - averageDepth = sum(subgridBathymetryValues(1:nSubgridCV)*subgridAreas(1:nSubgridCV))/sum(subgridAreas(1:nSubgridCV)) - ! real(nSubgridCV,RKIND) + !-------------------------------------------------------------------- + + averageDepth = sum(subgridBathymetryValues(1:nSubgridCV))/real(nSubgridCV,RKIND) computeVelAverage = .false. if (present(subgridUValues) .and. present(subgridVValues) .and. & @@ -1225,7 +1225,6 @@ subroutine ocn_init_wet_average(nSubgridCV, subgridBathymetryValues,subgridSshVa wetArea = 0.0_RKIND HU_int = 0.0_RKIND HV_int = 0.0_RKIND - do tri = 1,nSubgridCV layerThicknessValue = subgridBathymetryValues(tri) + subgridSshValues(tri) @@ -1247,10 +1246,13 @@ subroutine ocn_init_wet_average(nSubgridCV, subgridBathymetryValues,subgridSshVa endif subgridThicknessAverage = H_int/wetArea + !-------------------------------------------------------------------- + end subroutine ocn_init_wet_average!}}} - + + ! DW subroutine ocn_init_wet_average_ssh( nSubgridCV, subgridBathymetryValues, subgridSshValues, subgridAreas, sshWetAverage)!{{{ implicit none From 09741d51202bb82418b43e40576d49b3d19b85e2 Mon Sep 17 00:00:00 2001 From: Damrongsak Wirasaet Date: Tue, 21 Nov 2023 12:02:39 -0800 Subject: [PATCH 06/16] commit buttermilk bay before trying to pull new compass --- components/mpas-ocean/src/mode_init/Makefile | 3 +- .../mpas-ocean/src/mode_init/Registry.xml | 2 + .../src/mode_init/Registry_Buttermilk_Bay.xml | 72 + .../mode_init/mpas_ocn_init_buttermilk_bay.F | 1969 +++++++++++++++++ .../src/mode_init/mpas_ocn_init_mode.F | 6 +- .../mode_init/mpas_ocn_init_parabolic_bowl.F | 2 +- 6 files changed, 2051 insertions(+), 3 deletions(-) create mode 100644 components/mpas-ocean/src/mode_init/Registry_Buttermilk_Bay.xml create mode 100644 components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F diff --git a/components/mpas-ocean/src/mode_init/Makefile b/components/mpas-ocean/src/mode_init/Makefile index d1aba81e14f5..3f0e84f84a4a 100644 --- a/components/mpas-ocean/src/mode_init/Makefile +++ b/components/mpas-ocean/src/mode_init/Makefile @@ -31,7 +31,8 @@ TEST_CASES = mpas_ocn_init_baroclinic_channel.o \ mpas_ocn_init_mixed_layer_eddy.o \ mpas_ocn_init_transport_tests.o \ mpas_ocn_init_test_sht.o \ - mpas_ocn_init_parabolic_bowl.o + mpas_ocn_init_parabolic_bowl.o \ + mpas_ocn_init_buttermilk_bay.o #mpas_ocn_init_TEMPLATE.o all: init_mode diff --git a/components/mpas-ocean/src/mode_init/Registry.xml b/components/mpas-ocean/src/mode_init/Registry.xml index 3c05635a9c1b..99986d5a070a 100644 --- a/components/mpas-ocean/src/mode_init/Registry.xml +++ b/components/mpas-ocean/src/mode_init/Registry.xml @@ -20,6 +20,8 @@ #include "Registry_mixed_layer_eddy.xml" #include "Registry_test_sht.xml" #include "Registry_parabolic_bowl.xml" +#include "Registry_Buttermilk_Bay.xml" + // #include "Registry_TEMPLATE.xml" diff --git a/components/mpas-ocean/src/mode_init/Registry_Buttermilk_Bay.xml b/components/mpas-ocean/src/mode_init/Registry_Buttermilk_Bay.xml new file mode 100644 index 000000000000..5724a2e8ca13 --- /dev/null +++ b/components/mpas-ocean/src/mode_init/Registry_Buttermilk_Bay.xml @@ -0,0 +1,72 @@ + + + + + + + + + + + + + + + + + + + + + + + diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F new file mode 100644 index 000000000000..36ff26239884 --- /dev/null +++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F @@ -0,0 +1,1969 @@ +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the LICENSE file +! distributed with this code, or at http://mpas-dev.github.com/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_init_parabolic_bowl +! +!> \brief MPAS ocean initialize case -- TEMPLATE +!> \author D. Wirasaet, S. Brus +!> \date May-June 2022 +!> \details +!> This module contains the routines for initializing the +!> parabolic_bowl initial condition (Thacker's problem) +!> +!> +!....................................................................... +!----------------------------------------------------------------------- +! +!> In order to add a new initial condition, do the following: +!> 1. In src/core_ocean/mode_init, copy these to your new initial condition name: +!> cp mpas_ocn_init_TEMPLATE.F mpas_ocn_init_your_new_name.F +!> cp Registry_TEMPLATE.xml Registry_ocn_your_new_name.xml +!> +!> 2. In those two new files, replace the following text: +!> TEMPLATE, FILL_IN_AUTHOR, FILL_IN_DATE +!> TEMPLATE uses underscores (subroutine names), like your_new_name. +!> +!> 3. Add a #include line for your registry to +!> src/core_ocean/mode_init/Registry.xml +!> +!> 4. Copy and change TEMPLATE lines in src/core_ocean/mode_init/mpas_ocn_init_mode.F +!> +!> 5. Add these dependency lines by following TEMPLATE examples in: +!> in src/core_ocean/mode_init/Makefile +! +!----------------------------------------------------------------------- + +module ocn_init_Buttermilk_bay + + use mpas_kind_types + use mpas_io_units + use mpas_derived_types + use mpas_pool_routines + use mpas_constants + use mpas_dmpar + + use ocn_constants + use ocn_config + use ocn_init_vertical_grids + use ocn_init_cell_markers + use ocn_subgrid + + use mpas_constants + use mpas_io + use mpas_io_streams + use mpas_stream_manager + + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + !-------------------------------------------------------------------- + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + + public :: ocn_init_setup_Buttermilk_bay, & + ocn_init_validate_Buttermilk_bay + + !-------------------------------------------------------------------- + ! + ! Private module variables + ! + !-------------------------------------------------------------------- + +! real (kind=RKIND):: LL, CC +! real (kind=RKIND):: oneMC2, sqrtOneMC2, oneMC +! real (kind=RKIND), parameter :: eps=1.0e-10 + real(kind=RKIND), dimension(:,:), pointer :: & + subgridWetVolumeCellTable, & + subgridWetVolumeEdgeTable, & + subgridWetVolumeVertexTable, & + subgridSshCellTableRange, & + subgridSshEdgeTableRange, & + subgridSshVertexTableRange, & + subgridWetFractionCellTable, & + subgridWetFractionEdgeTable, & + subgridWetFractionVertexTable + real(kind=RKIND), dimension(:), pointer :: & + subgridEdgeBathymetryMean, & + subgridVertexBathymetryMean, & + subgridCellBathymetryMin, & + subgridEdgeBathymetryMin, & + subgridVertexBathymetryMin, & + subgridLayerThicknessDebug + integer, pointer :: nSubgridTableLevels + + ! For netcdf topobathy input variables + integer :: nLatTopo, nLonTopo + + type (field1DReal) :: topoLat, topoLon + type (field2DReal) :: topoIC + + real(kind=RKIND), parameter:: eps = 1.0e-10_RKIND ; +!*********************************************************************** + +contains + +!*********************************************************************** +! +! routine ocn_init_setup_parabolic_bowl +! +!> \brief Setup for this initial condition +!> \author D. Wirasaet and S. Brus +!> \date May-June 2022 +!> \details +!> This routine sets up the initial conditions for this case. +!> To be run in sigma vertical coordinates and single-layer +! +!----------------------------------------------------------------------- + + subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ + use mpas_vector_operations ! To calcutate edgeNormalVector + + implicit none + !-------------------------------------------------------------------- + + type (domain_type), intent(inout) :: domain + integer, intent(out) :: iErr + + type (block_type), pointer :: block_ptr + type (mpas_pool_type), pointer :: meshPool + type (mpas_pool_type), pointer :: statePool + type (mpas_pool_type), pointer :: tracersPool + type (mpas_pool_type), pointer :: verticalMeshPool + + ! local variables + integer :: iCell, iEdge, iVertex, k, idx + real (kind=RKIND) :: yMin, yMax, xMin, xMax, dcEdgeMin, dcEdgeMinGlobal + real (kind=RKIND) :: yMinGlobal, yMaxGlobal, yMidGlobal, xMinGlobal, xMaxGlobal + real (kind=RKIND) :: localVar1, localVar2 + real (kind=RKIND), dimension(:), pointer :: interfaceLocations + + ! Define dimension pointers + integer, pointer :: nCellsSolve, nEdgesSolve, nVerticesSolve, nVertLevels, nVertLevelsP1 + integer, pointer :: index_temperature, index_salinity + integer, pointer :: maxEdges + + ! Define variable pointers + logical, pointer :: on_a_sphere + integer, dimension(:), pointer :: minLevelCell, maxLevelCell + integer, dimension(:), pointer :: nEdgesOnCell + integer, dimension(:,:), pointer :: verticesOnCell, verticesOnEdge + integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex + real (kind=RKIND), dimension(:), pointer :: xCell, yCell, refBottomDepth, refZMid, & + vertCoordMovementWeights, bottomDepth, fCell, fEdge, fVertex, dcEdge + real (kind=RKIND), dimension(:,:), pointer:: zMid + + real (kind=RKIND), dimension(:), pointer:: xEdge, yEdge, xVertex, yVertex + real (kind=RKIND) :: minBottomDepth, maxBottomDepth, globalMaxBottomDepth, globalMinBottomDepth + real (kind=RKIND), dimension(:,:), pointer :: layerThickness, restingThickness + real (kind=RKIND), dimension(:,:,:), pointer :: activeTracers + + real (kind=RKIND), dimension(:), pointer :: ssh + real (kind=RKIND), dimension(:), pointer :: areaCell + real (kind=RKIND), dimension(:,:), pointer :: edgeNormalVectors + real (kind=RKIND), dimension(:,:), pointer :: normalVelocity + + + real (kind=RKIND):: HH, uu, vv + real (kind=RKIND):: RR, num, den + real (kind=RKIND):: xshift = 0.0, yshift = 0.0 + real (kind=RKIND) :: layerThicknessEdgeAverage + real (kind=RKIND), dimension(:,:), allocatable :: rSubgridPoints, sSubgridPoints + real (kind=RKIND), dimension(:), allocatable :: subgridBathymetryValues, subgridAreas + real (kind=RKIND), dimension(:), allocatable :: subgridSshValues + real (kind=RKIND), dimension(:), allocatable :: subgridUValues, subgridVValues + real (kind=RKIND), dimension(:), allocatable :: uVelocityAverage, vVelocityAverage + integer :: nSubgridCell, nSubgridEdge, nSubgridVertex + integer :: nSubgridTriPerSlice + integer :: v1, v2 + integer :: c1, c2 + real (kind=RKIND) :: x(3), y(3) + integer :: slice, nSlice + integer :: i,j + real (kind=RKIND) :: deltaZ +! DW + integer:: jj + integer:: nsubgridCellEdge, iEdgeSegment + real (kind=RKIND):: pi + real (kind=RKIND), dimension(:,:), allocatable :: cellEdgeBathymetryValues + real (kind=RKIND), dimension(:), allocatable:: dsEdge + real (kind=RKIND), dimension(:), allocatable:: xSubgridCell, ySubgridCell + + iErr = 0 + + call ocn_subgrid_init(domain,iErr) + + if(config_init_configuration .ne. trim('buttermilk_bay')) return + + ! Determine vertical grid for configuration + call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', meshPool) + call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + call mpas_pool_get_dimension(meshPool, 'nVertLevelsP1', nVertLevelsP1) + call mpas_pool_get_config(meshPool, 'on_a_sphere', on_a_sphere) + + nVertLevels = config_Buttermilk_bay_vert_levels ; + nVertLevelsP1 = nVertLevels + 1 + + allocate(interfaceLocations(nVertLevelsP1)) + call ocn_generate_vertical_grid( config_vertical_grid, interfaceLocations, ocnConfigs ) ; + !! Mental note: interfaceLocatons = (k-1)/N ; + + ! Initalize min/max values to large positive and negative values + yMin = 1.0E10_RKIND + yMax = -1.0E10_RKIND + xMin = 1.0E10_RKIND + xMax = -1.0E10_RKIND + dcEdgeMin = 1.0E10_RKIND + + ! Determine local min and max values. + block_ptr => domain % blocklist + do while(associated(block_ptr)) + call mpas_pool_get_subpool(block_ptr % structs, 'mesh', meshPool) + + call mpas_pool_get_dimension( meshPool, 'nCellsSolve', nCellsSolve ) + call mpas_pool_get_dimension( meshPool, 'nEdgesSolve', nEdgesSolve ) + + call mpas_pool_get_array(meshPool, 'xCell', xCell) + call mpas_pool_get_array(meshPool, 'yCell', yCell) + call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge) + + yMin = min( yMin, minval(yCell(1:nCellsSolve))) + yMax = max( yMax, maxval(yCell(1:nCellsSolve))) + xMin = min( xMin, minval(xCell(1:nCellsSolve))) + xMax = max( xMax, maxval(xCell(1:nCellsSolve))) + dcEdgeMin = min( dcEdgeMin, minval(dcEdge(1:nEdgesSolve))) + + block_ptr => block_ptr % next + end do + + ! Determine global min and max values. + call mpas_dmpar_min_real(domain % dminfo, yMin, yMinGlobal) + call mpas_dmpar_max_real(domain % dminfo, yMax, yMaxGlobal) + call mpas_dmpar_min_real(domain % dminfo, xMin, xMinGlobal) + call mpas_dmpar_max_real(domain % dminfo, xMax, xMaxGlobal) + call mpas_dmpar_min_real(domain % dminfo, dcEdgeMin, dcEdgeMinGlobal) + + pi = acos(-1.0_RKIND) ; + xshift = 0.0_RKIND ; + yshift = -(YMin + dcEdgeMin*sin(pi/3.0_RKIND)) ; + + !*********************************************************************** + ! + ! Topography + ! + !*********************************************************************** + + call mpas_log_write( 'Reading bathymetry from a NetCDF file') + + if (config_Buttermilk_bay_topography_source == 'latlon_file' .or. & + config_Buttermilk_bay_topography_source == 'xy_file' ) then + call mpas_log_write( 'Reading topography data from file.') + call ocn_init_setup_Buttermilk_bay_read_topo(domain, iErr) + endif + + !-------------------------------------------------------------------- + ! Use this section to set initial values + !-------------------------------------------------------------------- + + block_ptr => domain % blocklist + do while(associated(block_ptr)) + call mpas_pool_get_subpool(block_ptr % structs, 'mesh', meshPool) + call mpas_pool_get_subpool(block_ptr % structs, 'state', statePool) + call mpas_pool_get_subpool(block_ptr % structs, 'verticalMesh', verticalMeshPool) + call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) + + call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) ; + call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) ; + call mpas_pool_get_dimension(meshPool, 'nEdgesSolve', nEdgesSolve) ; + call mpas_pool_get_dimension(meshPool, 'nVerticesSolve', nVerticesSolve) ; + call mpas_pool_get_dimension(meshPool, 'maxEdges', maxEdges) + + call mpas_pool_get_dimension(tracersPool, 'index_temperature', index_temperature) + call mpas_pool_get_dimension(tracersPool, 'index_salinity', index_salinity) + + call mpas_pool_get_array(meshPool, 'xCell', xCell) + call mpas_pool_get_array(meshPool, 'yCell', yCell) + call mpas_pool_get_array(meshPool, 'refBottomDepth', refBottomDepth) + call mpas_pool_get_array(meshPool, 'vertCoordMovementWeights', vertCoordMovementWeights) + call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) + call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) + call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) + call mpas_pool_get_array(meshPool, 'areaCell', areaCell) + + call mpas_pool_get_array(meshPool, 'verticesOnCell', verticesOnCell) + call mpas_pool_get_array(meshPool, 'verticesOnEdge', verticesOnEdge) + call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) + call mpas_pool_get_array(meshPool, 'cellsOnVertex', cellsOnVertex) + + call mpas_pool_get_array(meshPool, 'fCell', fCell) + call mpas_pool_get_array(meshPool, 'fEdge', fEdge) + call mpas_pool_get_array(meshPool, 'fVertex', fVertex) + + call mpas_pool_get_array(meshPool, 'xEdge', xEdge ) + call mpas_pool_get_array(meshPool, 'yEdge', yEdge ) + call mpas_pool_get_array(meshPool, 'xVertex', xVertex ) + call mpas_pool_get_array(meshPool, 'yVertex', yVertex ) + + call mpas_pool_get_array(statePool, 'zMid', zMid, 1) ; + + call mpas_pool_get_array(statePool, 'ssh', ssh, 1) + call mpas_pool_get_array(meshPool, 'edgeNormalVectors', edgeNormalVectors ) ; + call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity ) ; + + call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracers, 1) + call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 1) + + call mpas_pool_get_array(verticalMeshPool, 'refZMid', refZMid) + call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness) + + call mpas_pool_get_array(meshPool, 'subgridWetVolumeCellTable', & + subgridWetVolumeCellTable) + call mpas_pool_get_array(meshPool, 'subgridWetVolumeEdgeTable', & + subgridWetVolumeEdgeTable) + call mpas_pool_get_array(meshPool, 'subgridWetVolumeVertexTable', & + subgridWetVolumeVertexTable) + call mpas_pool_get_array(meshPool, 'subgridWetFractionCellTable', & + subgridWetFractionCellTable) + call mpas_pool_get_array(meshPool, 'subgridWetFractionEdgeTable', & + subgridWetFractionEdgeTable) + call mpas_pool_get_array(meshPool,'subgridWetFractionVertexTable',& + subgridWetFractionVertexTable) + call mpas_pool_get_array(meshPool, 'subgridSshCellTableRange', & + subgridSshCellTableRange) + call mpas_pool_get_array(meshPool, 'subgridSshEdgeTableRange', & + subgridSshEdgeTableRange) + call mpas_pool_get_array(meshPool, 'subgridSshVertexTableRange', & + subgridSshVertexTableRange) + call mpas_pool_get_array(meshPool, 'subgridEdgeBathymetryMean', & + subgridEdgeBathymetryMean) + call mpas_pool_get_array(meshPool, 'subgridVertexBathymetryMean', & + subgridVertexBathymetryMean) + call mpas_pool_get_array(meshPool, 'subgridCellBathymetryMin', & + subgridCellBathymetryMin) + call mpas_pool_get_array(meshPool, 'subgridEdgeBathymetryMin', & + subgridEdgeBathymetryMin) + call mpas_pool_get_array(meshPool, 'subgridVertexBathymetryMin', & + subgridVertexBathymetryMin) + call mpas_pool_get_dimension(meshPool, 'nSubgridTableLevels', & + nSubgridTableLevels) + call mpas_pool_get_array(meshPool, 'subgridLayerThicknessDebug', & + subgridLayerThicknessDebug) + + + ! if config_parabolic_bowl_adjust_domain_center == .true., + ! Adjust center of the mesh so that its center is located at (0,0) + if ( config_Buttermilk_bay_adjust_domain ) then + xCell = xCell - xshift ; + yCell = yCell - yshift ; + + xEdge = xEdge - xshift ; + yEdge = yEdge - yshift ; + + xVertex = xVertex - xshift ; + yVertex = yVertex - yshift ; + end if + + ! Initlialze vector + call mpas_initialize_vectors(meshPool) ; + + minLevelCell(:) = 1 + do iCell = 1, nCellsSolve + ! Set up vertical grid + maxLevelCell(iCell) = nVertLevels ; ! sigma coordinates + end do + + + do iCell = 1, nCellsSolve + ! Set temperature + activeTracers(index_temperature, :, iCell) = 10.0_RKIND + + ! Set salinity + activeTracers(index_salinity, :, iCell) = 30.0_RKIND + + ! Set Coriolis parameters, if other than zero + fCell(iCell) = config_Buttermilk_bay_coriolis_parameter ; + end do + + do iEdge = 1, nEdgesSolve + fEdge(iEdge) = config_Buttermilk_bay_coriolis_parameter ; + end do + + do iVertex = 1, nVerticesSolve + fVertex(iVertex) = config_Buttermilk_bay_coriolis_parameter ; + end do + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Cells + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + print*, "Begin cells", nCellsSolve + + nSubgridTriPerSlice = config_Buttermilk_bay_subgrid_level**2 + allocate(rSubgridPoints(3,maxEdges*nSubgridTriPerSlice), sSubgridPoints(3,maxEdges*nSubgridTriPerSlice)) + call ocn_init_define_subgrid_points(config_Buttermilk_bay_subgrid_level, rSubgridPoints, sSubgridPoints) + allocate(subgridBathymetryValues(maxEdges*nSubgridTriPerSlice), subgridAreas(maxEdges*nSubgridTriPerSlice)) + allocate(subgridSshValues(maxEdges*nSubgridTriPerSlice)) + allocate(subgridUValues(maxEdges*nSubgridTriPerSlice)) + allocate(subgridVValues(maxEdges*nSubgridTriPerSlice)) + + + open( unit = 101, file = 'BM_bathy.dat', action = 'write' ) ; + + do iCell = 1,nCellsSolve + ! + ! Evaluate subgrid bathymetry at centers of sub-triangles for cell slices + ! (all subdivided triangles for each cell slice are gathered into + ! subgridBathymetryValues and subgridAreas) + !--------------------------------------------------------------- + + nSubgridCell = 0 ! Counter for all subgrid triangles over cell slices + do slice = 1,nEdgesOnCell(iCell) ! Loop over cell slices + + v1 = verticesOnCell(slice,iCell) + if (slice+1 <= nEdgesOnCell(iCell)) then + v2 = verticesOnCell(slice+1,iCell) + else + v2 = verticesOnCell(1,iCell) + endif + + ! Cell slice coordinates + x(1) = xCell(iCell) + y(1) = yCell(iCell) + + x(2) = xVertex(v1) + y(2) = yVertex(v1) + + x(3) = xVertex(v2) + y(3) = yVertex(v2) + + call ocn_init_evaluate_subgrid_data(x, y, nSubgridTriPerSlice, nSubgridCell, rSubgridPoints, sSubgridPoints, & + subgridBathymetryValues, subgridAreas, subgridSshValues) + + enddo + + do jj = 1, nSubgridCell + WRITE( 101, '(I6,2F16.8)' ) iCell, subgridBathymetryValues(jj) ; + end do + + ! Evaluate bounds of look-up table range + !--------------------------------------------------------------- + + subgridSshCellTableRange(1,iCell) = -maxval(subgridBathymetryValues(1:nSubgridCell)) + config_drying_min_cell_height + eps + subgridSshCellTableRange(2,iCell) = -minval(subgridBathymetryValues(1:nSubgridCell)) + + ! Evaluate subgrid bathymetry + !--------------------------------------------------------------- + ! DW + bottomDepth(iCell) = sum(subgridBathymetryValues(1:nSubgridCell)*subgridAreas(1:nSubgridCell))/sum(subgridAreas(1:nSubgridCell)) + subgridCellBathymetryMin(iCell) = maxval(subgridBathymetryValues(1:nSubgridCell)) + + ! Vertical integration of wet fraction + !--------------------------------------------------------------- + + call ocn_init_vertical_integration(iCell,subgridSshCellTableRange, nSubgridCell, subgridBathymetryValues, subgridAreas, & + subgridWetVolumeCellTable, subgridWetFractionCellTable) + + ! Evaluate wet layerThickness average + !--------------------------------------------------------------- + call ocn_init_wet_average(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, layerThickness(1,iCell)) + + +!!$ call ocn_init_wet_average_ssh(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, ssh(iCell))!{{{ + enddo + close( unit = 101 ) ; + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Edges + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + print*, "Begin edges", nEdgesSolve + + allocate(uVelocityAverage(nEdgesSolve)) + allocate(vVelocityAverage(nEdgesSolve)) + + nSubgridCellEdge=config_Buttermilk_bay_subgrid_level + allocate( cellEdgeBathymetryValues(3,nSubgridCellEdge) ) + allocate( dsEdge(nSubgridCellEdge) ) + cellEdgeBathymetryValues = -99999 ; + + allocate(xSubgridCell(maxEdges*nSubgridTriPerSlice)) + allocate(ySubgridCell(maxEdges*nSubgridTriPErSlice)) + + + do iEdge = 1,nEdgesSolve + !-------------------------------------------------------------- + ! Evaluate subgrid bathymetry at centers of sub-triangles for edge slices + ! (all subdivided triangles for each edge slice are gathered into + ! subgridBathymetryValues and subgridAreas) + !--------------------------------------------------------------- + + nSlice = 0 + do slice = 1,2 + if (cellsOnEdge(slice,iEdge) <= nCellsSolve) then + nSlice = nSlice + 1 + endif + enddo + + nSubgridEdge = 0 ! Counter for all subgrid triangles over edge slices + do slice = 1,nSlice ! Loop over edge slices + + ! Edge slice coordinates + x(1) = xVertex(verticesOnEdge(1,iEdge)) + y(1) = yVertex(verticesOnEdge(1,iEdge)) + + x(2) = xVertex(verticesOnEdge(2,iEdge)) + y(2) = yVertex(verticesOnEdge(2,iEdge)) + + x(3) = xCell(cellsOnEdge(slice,iEdge)) + y(3) = yCell(cellsOnEdge(slice,iEdge)) + + + call ocn_init_evaluate_subgrid_data(x, y, nSubgridTriPerSlice, nSubgridEdge, rSubgridPoints, sSubgridPoints, & + subgridBathymetryValues, subgridAreas, subgridSshValues, subgridUValues, subgridVValues, & + xSubgridCell, ySubgridCell ) + + enddo + + ! Evaluate velocity average + !--------------------------------------------------------------- + call ocn_init_wet_average(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, layerThicknessEdgeAverage, & + subgridUValues, subgridVValues, uVelocityAverage(iEdge), vVelocityAverage(iEdge)) + + if ( .NOT. config_Buttermilk_bay_subgrid_edge_bathymetry_max_pixel ) then + ! Evaluate bounds of look-up table range + !--------------------------------------------------------------- + + subgridSshEdgeTableRange(1,iEdge) = -maxval(subgridBathymetryValues(1:nSubgridEdge)) + config_drying_min_cell_height + eps + subgridSshEdgeTableRange(2,iEdge) = -minval(subgridBathymetryValues(1:nSubgridEdge)) + + + ! Evaluate subgrid bathymetry + !--------------------------------------------------------------- + + subgridEdgeBathymetryMean(iEdge) = sum(subgridBathymetryValues(1:nSubgridEdge))/real(nSubgridEdge,RKIND) + subgridEdgeBathymetryMin(iEdge) = maxval(subgridBathymetryValues(1:nSubgridEdge)) + + ! Vertical integration of wet fraction + !--------------------------------------------------------------- + + call ocn_init_vertical_integration(iEdge,subgridSshEdgeTableRange, nSubgridEdge, subgridBathymetryValues, subgridAreas, & + subgridWetVolumeEdgeTable, subgridWetFractionEdgeTable) + else + ! DW: Use the higher values of the pair of subcells along the + ! cell edge + cellEdgeBathymetryValues(1,:) = subgridBathymetryValues(1:2*nSubgridCellEdge - 1:2) + if ( nslice > 1 ) then + cellEdgeBathymetryValues(2,:) = & + subgridBathymetryValues(nSubgridTriPerSlice+2*nsubgridCellEdge - 1:nSubgridTriPerSlice+1:-2) + else + cellEdgeBathymetryValues(2,:) = cellEdgeBathymetryValues(1,:) + endif + + do iEdgeSegment = 1, nSubgridCellEdge + cellEdgeBathymetryValues(3,iEdgeSegment) = minval(cellEdgeBathymetryValues(1:2,iEdgeSegment) ) + end do + dsEdge(:) = sqrt( (x(2) - x(1))*(x(2) - x(1)) + (y(2) - y(1))*(y(2)- y(1)) )/nSubgridCellEdge + + + ! Evaluate bounds of look-up table range + !--------------------------------------------------------------- + + subgridSshEdgeTableRange(1,iEdge) = -maxval(CellEdgeBathymetryValues(3,1:nSubgridCellEdge)) + config_drying_min_cell_height + eps + subgridSshEdgeTableRange(2,iEdge) = -minval(CellEdgeBathymetryValues(3,1:nSubgridCellEdge)) + + ! Evaluate bounds of look-up table range + !--------------------------------------------------------------- + + subgridEdgeBathymetryMean(iEdge) = sum(cellEdgeBathymetryValues(3,1:nSubgridCellEdge))/real(nSubgridCellEdge,RKIND) + subgridEdgeBathymetryMin(iEdge) = maxval(cellEdgeBathymetryValues(3,1:nSubgridCellEdge)) + + ! Vertical integration of wet fraction + !--------------------------------------------------------------- + call ocn_init_vertical_integration( iEdge, subgridSshEdgeTableRange, & + nSubgridCellEdge, cellEdgeBathymetryValues(3,:), dsEdge, subgridWetVolumeEdgeTable, subgridWetFractionEdgeTable ) + endif + + end do + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Vertex + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + print*, "Begin vertex", nVerticesSolve + +vertex:do iVertex = 1,nVerticesSolve + + ! Evaluate subgrid bathymetry at centers of sub-triangles for vertex triangle + !--------------------------------------------------------------- + + nSlice = 0 + do slice = 1,3 + if (cellsOnVertex(slice,iVertex) <= nCellsSolve) then + nSlice = nSlice + 1 + endif + enddo + + if (nSlice < 3) then + cycle vertex + endif + + nSubgridVertex = 0 ! Counter for all subgrid triangles over edge slices + do slice = 1,nSlice + + c1 = cellsOnVertex(slice,iVertex) + if (slice < 3) then + c2 = cellsOnVertex(slice+1,iVertex) + else + c2 = cellsOnVertex(1,iVertex) + endif + + ! Vertex slice coordinates + x(1) = xCell(c1) + y(1) = yCell(c1) + + x(2) = xCell(c2) + y(2) = yCell(c2) + + x(3) = xVertex(iVertex) + y(3) = yVertex(iVertex) + + call ocn_init_evaluate_subgrid_data(x, y, nSubgridTriPerSlice, nSubgridVertex, rSubgridPoints, sSubgridPoints, & + subgridBathymetryValues, subgridAreas) + enddo + + ! Evaluate bounds of look-up table range + !--------------------------------------------------------------- + + subgridSshVertexTableRange(1,iVertex) = -maxval(subgridBathymetryValues(1:nSubgridVertex)) + config_drying_min_cell_height + eps + subgridSshVertexTableRange(2,iVertex) = -minval(subgridBathymetryValues(1:nSubgridVertex)) + + ! Evaluate subgrid bathymetry + !--------------------------------------------------------------- + subgridVertexBathymetryMean(iVertex) = sum(subgridBathymetryValues(1:nSubgridVertex)*subgridAreas(1:nSubgridVertex))/sum(subgridAreas(1:nSubgridVertex)) + subgridVertexBathymetryMin(iVertex) = maxval(subgridBathymetryValues(1:nSubgridVertex)) + + + ! Vertical integration of wet fraction + !--------------------------------------------------------------- + + call ocn_init_vertical_integration(iVertex,subgridSshVertexTableRange, nSubgridVertex, subgridBathymetryValues, subgridAreas, & + subgridWetVolumeVertexTable, subgridWetFractionVertexTable) + + enddo vertex + + ! Find max bottom depth + maxBottomDepth = maxval( bottomDepth ) ; + minBottomDepth = minval( bottomDepth ) ; + call mpas_dmpar_max_real( domain % dminfo, maxBottomDepth, globalMaxBottomDepth ) ; + call mpas_dmpar_min_real( domain % dminfo, minBottomDepth, globalMinBottomDepth ) ; + + ! Set refBottomDepth and refZMid + do k = 1, nVertLevels + refBottomDepth(k) = globalMaxBottomDepth*interfaceLocations(k+1) ; + refZMid(k) = -0.5_RKIND*( interfaceLocations(k+1) + interfaceLocations(k))*globalMaxBottomDepth ; + end do + + ! Set vertCoordMovementWeights + vertCoordMovementWeights(:) = 1.0_RKIND + + ! Set velocity + do iEdge = 1, nEdgesSolve + + if (config_use_subgrid_wetting_drying) then + do k = 1, nVertLevels + normalVelocity(k,iEdge) = uVelocityAverage(iEdge)*edgeNormalVectors(1,iEdge) & + + vVelocityAverage(iEdge)*edgeNormalVectors(2,iEdge) ; + end do + else + call ocn_init_Buttermilk_bay_velocity(xEdge(iEdge), yEdge(iEdge), uu, vv) + do k = 1, nVertLevels + normalVelocity(k,iEdge) = uu*edgeNormalVectors(1,iEdge) + vv*edgeNormalVectors(2,iEdge) ; + end do + end if + end do + + ! Set layer thickness and ssh + if (config_use_wetting_drying) then + + do iCell = 1, nCellsSolve + ! Set up vertical grid + maxLevelCell(iCell) = nVertLevels ; ! sigma coordinates + end do + + do iCell = 1, nCellsSolve + ! + ! make sure depth is thick enough via ssh = TOTAL_DEPTH - bottomDepth + ! add a thin layer of nlayer*config_drying_min_cellhight + ! + + if (config_use_subgrid_wetting_drying) then + + call ocn_subgrid_ssh_lookup(layerThickness(1,iCell),& + subgridWetVolumeCellTable(:,iCell),& + subgridSshCellTableRange(:,iCell),& + bottomDepth(iCell),& + subgridCellBathymetryMin(iCell),& + ssh(iCell)) +! call ocn_subgrid_layer_thickness_lookup(ssh(iCell), & +! subgridWetVolumeCellTable(:,iCell), & +! subgridSshCellTableRange(:,iCell),& +! bottomDepth(iCell),& +! LayerThickness(1,iCell)) + + else + call ocn_init_Buttermilk_bay_bathymetry(xCell(iCell),yCell(iCell),bottomDepth(iCell)) + call ocn_init_Buttermilk_bay_ssh(xCell(iCell),yCell(iCell),bottomDepth(iCell),ssh(iCell)) + ssh(iCell) = - bottomDepth(iCell) + & + max(ssh(iCell) + bottomDepth(iCell), & + maxLevelCell(iCell)*(config_drying_min_cell_height + eps)) + + do k = 1, maxLevelCell(iCell) + layerThickness(k,iCell) = max(config_drying_min_cell_height + eps, & + (ssh(iCell) + bottomDepth(iCell))/real(maxLevelCell(iCell),RKIND)) + + if (layerThickness(k,iCell) < config_drying_min_cell_height) then + call mpas_log_write('layerThickness($i,$i)=$r', MPAS_LOG_CRIT, & + intArgs=(/k,iCell/), & + realArgs=(/layerThickness(k,iCell)/)) + end if + end do + endif + + + do k = 1, maxLevelCell(iCell) + restingThickness(k,iCell) = bottomDepth(iCell)/maxLevelCell(iCell) + end do + end do + + end if + + block_ptr => block_ptr % next + end do + + do iCell = 1,nCellsSolve + call ocn_subgrid_layer_thickness_lookup(ssh(iCell), & + subgridWetVolumeCellTable(:,iCell), & + subgridSshCellTableRange(:,iCell),& + bottomDepth(iCell),& + subgridLayerThicknessDebug(iCell)) + enddo + ! + + + deallocate(interfaceLocations) + if (config_global_ocean_topography_source == 'latlon_file') then + call mpas_log_write( 'Cleaning up topography IC fields') + call ocn_init_Buttermilk_bay_destroy_topo_fields() + endif + !-------------------------------------------------------------------- + + print*, "****** End Butter milk bay init *****" ; + + return ; + end subroutine ocn_init_setup_Buttermilk_bay!}}} + +!*********************************************************************** +! +! routine ocn_init_define_subgrid_points +! +!> \brief Define subgrid points on reference triangle +!> \author Steven Brus +!> \date November 2022 +!> \details Gives the r and s coordinates of the subgrid triangle vertices +!> +! +!----------------------------------------------------------------------- + + subroutine ocn_init_define_subgrid_points(nSubgridLevel, rSubgridPoints, sSubgridPoints)!{{{ + + implicit none + + integer, intent(in) :: nSubgridLevel + real (kind=RKIND), dimension(:,:), intent(inout) :: rSubgridPoints, sSubgridPoints + + integer :: i,j,k + integer :: nSubgridTri + real (kind=RKIND) :: dx + real (kind=RKIND), dimension(:), allocatable :: xPoints + + !-------------------------------------------------------------------- + + ! Equi-spaced nodes on -1,+1 + allocate(xPoints(nSubgridLevel+1)) + dx = 2.0_RKIND/real(nSubgridLevel,RKIND) + xPoints(1) = -1.0_RKIND + do i = 2,nSubgridLevel+1 + xPoints(i) = xPoints(i-1)+dx + enddo + + ! * + ! |\ + ! | \ + ! | \ + ! | \ + ! s *----* + ! |\ u |\ + ! | \ | \ + ! | \ | \ + ! | l \| l \ + ! *----*----* + ! r + + ! Trianglulate tensor product of equi-spaced nodes on reference triangle + ! r - horizontal coordinate (i index), s - vertical coordinate (j index) + nSubgridTri = 0 + do j = 1,nSubgridLevel + do i = 1,nSubgridLevel+1 - j + + nSubgridTri = nSubgridTri+1 + + ! lower triangle in pair (see l in triangle above) + rSubgridPoints(1,nSubgridTri) = xPoints(i) + sSubgridPoints(1,nSubgridTri) = xPoints(j) + + rSubgridPoints(2,nSubgridTri) = xPoints(i+1) + sSubgridPoints(2,nSubgridTri) = xPoints(j) + + rSubgridPoints(3,nSubgridTri) = xPoints(i) + sSubgridPoints(3,nSubgridTri) = xPoints(j+1) + + ! upper triangle in pair (see u in triangle above, doesn't occur next to hypotenuse) + if (i < nSubgridLevel+1 - j) then + + nSubgridTri = nSubgridTri + 1 + + rSubgridPoints(1,nSubgridTri) = xPoints(i+1) + sSubgridPoints(1,nSubgridTri) = xPoints(j) + + rSubgridPoints(2,nSubgridTri) = xPoints(i+1) + sSubgridPoints(2,nSubgridTri) = xPoints(j+1) + + rSubgridPoints(3,nSubgridTri) = xPoints(i) + sSubgridPoints(3,nSubgridTri) = xPoints(j+1) + + endif + + enddo + enddo + + print*, "exit ocn_init_define_subgrid_points" + + + !-------------------------------------------------------------------- + + end subroutine ocn_init_define_subgrid_points!}}} + +!*********************************************************************** +! +! routine ocn_init_evaluate_subgrid_data +! +!> \brief Evaluate subgrid infromation +!> \author Steven Brus +!> \date November 2022 +!> \details Evaluate data on subgrid triangles for a given triangluar grid cell region +!> +! +!----------------------------------------------------------------------- + + +! subroutine ocn_init_evaluate_subgrid_data(xTri, yTri, nSubgridTri, nSubgridCV, rSubgridPoints, sSubgridPoints, & +! subgridBathymetryValues, subgridAreas, & +! subgridSshValues, subgridUValues, subgridVValues)!{{{ + subroutine ocn_init_evaluate_subgrid_data( xTri, yTri, nSubgridTri, nSubgridCV, rSubgridPoints, sSubgridPoints, & + subgridBathymetryValues, subgridAreas, & + subgridSshValues, subgridUValues, subgridVValues, & + xSubgridCell, ySubgridCell ) !{{{ + + + implicit none + real (kind=RKIND), intent(inout) :: xTri(3), yTri(3) + integer, intent(in) :: nSubgridTri + integer, intent(inout) :: nSubgridCV + real (kind=RKIND), dimension(:,:), intent(in) :: rSubgridPoints, sSubgridPoints + real (kind=RKIND), dimension(:), intent(inout) :: subgridBathymetryValues + real (kind=RKIND), dimension(:), intent(inout) :: subgridAreas + real (kind=RKIND), dimension(:), intent(inout),optional :: subgridSshValues + real (kind=RKIND), dimension(:), intent(inout),optional :: subgridUValues, subgridVValues + + real (kind=RKIND), dimension(:), intent(inout), optional :: xSubgridCell, ySubgridCell + + real (kind=RKIND) :: rCenter, sCenter + real (kind=RKIND) :: xCenter, yCenter + real (kind=RKIND) :: xSubgridPoints(3), ySubgridPoints(3) + real (kind=RKIND) :: x, y + real (kind=RKIND) :: area + + + integer :: iPt,i + + !-------------------------------------------------------------------- + + ! Coordinates of physical triangle + ! (ensure counter-clockwise numering) + call ocn_init_tri_area(xTri(:), yTri(:), area) + if (area < 0.0_RKIND) then + x = xTri(1) + y = yTri(1) + xTri(1) = xTri(2) + yTri(1) = yTri(2) + xTri(2) = x ; + yTri(2) = y ; + area = abs(area) ; + endif + + do iPt = 1,nSubgridTri + + ! Counter over all subcells within cell/edge/vertex control volume + nSubgridCV = nSubgridCV + 1 + + ! Center sub-triangle (on reference triangle) + rCenter = sum(rSubgridPoints(:,iPt))/3.0_RKIND + sCenter = sum(sSubgridPoints(:,iPt))/3.0_RKIND + + ! Transformation of sub-triangle center to physical coordinates + call ocn_init_tri_coordinate_transform(rCenter, sCenter, xTri, yTri, xCenter, yCenter) + + ! Evaluate bathymetry + call ocn_init_Buttermilk_bay_bathymetry(xCenter, yCenter, subgridBathymetryValues(nSubgridCV)) + + ! Transformation of sub-triangle vertices to physical coordinates + do i = 1,3 + call ocn_init_tri_coordinate_transform(rSubgridPoints(i,iPt), sSubgridPoints(i,iPt), xTri, yTri, xSubgridPoints(i), ySubgridPoints(i)) + enddo + + ! Calculate area of sub-triangle + call ocn_init_tri_area(xSubgridPoints(:), ySubgridPoints(:), subgridAreas(nSubgridCV)) + + ! Optionally evalulate ssh + if (present(subgridsshValues)) then + call ocn_init_Buttermilk_bay_ssh(xCenter, yCenter, subgridBathymetryValues(nSubgridCV), subgridSshValues(nSubgridCV)) + endif + + if (present(subgridUValues).and.present(subgridVValues)) then + call ocn_init_Buttermilk_bay_velocity(xCenter, yCenter, subgridUValues(nSubgridCV), subgridVValues(nSubgridCV)) + endif + + if ( present(xSubgridCell) ) then + xSubgridCell(nSubgridCV) = xCenter ; + endif + + if ( present(ySubgridCell) ) then + ySubgridCell(nSubgridCV) = yCenter ; + end if + enddo + + !-------------------------------------------------------------------- + + end subroutine ocn_init_evaluate_subgrid_data!}}} + +!*********************************************************************** +! +! routine ocn_init_parabolic_bowl_bathymetry +! +!> \brief Compute bathymetry +!> \author Steven Brus, D. Wirasaet +!> \date November 2022 +!> \details Return the value of the bathymetry at a given x,y point +!> +! +!----------------------------------------------------------------------- + + subroutine ocn_init_Buttermilk_bay_bathymetry(x, y, depth)!{{{ + + implicit none + + real (kind=RKIND), intent(in) :: x, y + real (kind=RKIND), intent(out) :: depth + + integer:: ix(2) + real (kind=RKIND) :: xc(2) + real (kind=RKIND) :: dx(2), x0(2), xN(2), bath(4), f(4), den + real (kind=RKIND) :: xv(2), yv(2) + + !-------------------------------------------------------------------- + ! + ! DEM must be on a uniform grid + x0 = (/ topoLon%array(1), topoLat%array(1) /) ; + xN = (/ topoLon%array(nLonTopo), topoLat%array(nLatTopo) /) ; + dx = (/ topoLon%array(2) - topoLon%array(1), & + topoLat%array(2) - topoLat%array(1) /) ; + + xc = (/ x, y /) ; + ix = floor( (xc - x0)/dx ) + 1 ; + + if ( ( ix(1) >= 1 .and. ix(1) < nLonTopo ) .and. & + ( ix(2) >= 1 .and. ix(2) < nLatTopo ) ) then + ! include the west and soutth bourdaries of + ! a given uniform raster DEM + bath(1) = TopoIC%array(ix(1),ix(2)) ; + bath(2) = TopoIC%array(ix(1)+1,ix(2)) ; + bath(3) = TopoIC%array(ix(1)+1,ix(2)+1) ; + bath(4) = TopoIC%array(ix(1),ix(2)) ; + + xv(1) = topoLon%array(ix(1)) ; + xv(2) = topoLon%array(ix(1)+1) ; + + yv(1) = topoLat%array(ix(2)) ; + yv(2) = topoLat%array(ix(2)+1) ; + + den = dx(1)*dx(2) ; + + f(1) = ( xc(1) - xv(2) )*( xc(2) - yv(2) )/den ; + f(2) = ( xc(1) - xv(1) )*( xc(2) - yv(2) )/(-den) ; + f(3) = ( xc(1) - xv(1) )*( xc(2) - yv(1) )/den ; + f(4) = ( xc(1) - xv(2) )*( xc(2) - yv(1) )/(-den) ; + + depth = sum( bath*f ) ; + else + ! nearest extrapolation + ix(1) = merge( 1, ix(1), ix(1) < 1 ) ; + ix(1) = merge( nLonTopo, ix(1), ix(1) >= nLonTopo ) ; + ix(2) = merge( 1, ix(2), ix(2) < 1 ) ; + ix(2) = merge( nLatTopo, ix(2), ix(2) >= nLatTopo ) ; + + depth = TopoIC%array(ix(1),ix(2)) ; + endif + + return ; + end subroutine ocn_init_Buttermilk_bay_bathymetry!}}} + +!*********************************************************************** +! +! routine ocn_init_parabolic_bowl_ssh +! +!> \brief Compute initial ssh field +!> \author Steven Brus, D. Wirasaet +!> \date November 2022 +!> \details Use exact solution to compute ssh field for initial conditions +!> +! +!----------------------------------------------------------------------- + + subroutine ocn_init_Buttermilk_bay_ssh(x, y, bottomDepth, ssh)!{{{ + + implicit none + + real (kind=RKIND), intent(in) :: x, y + real (kind=RKIND), intent(in) :: bottomDepth + real (kind=RKIND), intent(out) :: ssh + real (kind=RKIND) :: RR + + !-------------------------------------------------------------------- + +! RR = sqrt(x**2 + y**2) ; +! +! ssh = (sqrtOneMC2/oneMC) - 1.0_RKIND - ((RR**2)/(LL**2))*( (oneMC2/(oneMC**2)) - 1.0_RKIND ) ; +! ssh = config_parabolic_bowl_b0*ssh ; +! !ssh = - bottomDepth + max(ssh + bottomDepth, config_drying_min_cell_height + eps) + + ssh = 0.0_RKIND ; + ssh = merge( -bottomDepth + eps, ssh, ssh + bottomDepth < eps ) ; + + + !-------------------------------------------------------------------- + + end subroutine ocn_init_Buttermilk_bay_ssh!}}} + +!*********************************************************************** +! +! routine ocn_init_parabolic_bowl_velocity +! +!> \brief Compute initial velocity field +!> \author Steven Brus, D. Wirasaet +!> \date November 2022 +!> \details Use exact solution to comupte velocity field for initial conditions +!> +! +!----------------------------------------------------------------------- + + subroutine ocn_init_Buttermilk_bay_velocity(x, y, u, v)!{{{ + + implicit none + + real (kind=RKIND), intent(in) :: x, y + real (kind=RKIND), intent(out) :: u, v + real (kind=RKIND) :: RR, HH + + !-------------------------------------------------------------------- + + u = 0.0_RKIND ; + v = 0.0_RKIND ; + !-------------------------------------------------------------------- + + end subroutine ocn_init_Buttermilk_bay_velocity!}}} + +!*********************************************************************** +! +! routine ocn_init_tri_area +! +!> \brief Compute triangle area +!> \author Steven Brus +!> \date November 2022 +!> \details Compute area of triangle given vertex coordinates +!> +! +!----------------------------------------------------------------------- + + subroutine ocn_init_tri_area(x, y, area)!{{{ + + implicit none + + real (kind=RKIND), intent(in) :: x(3), y(3) + real (kind=RKIND), intent(out) :: area + + !-------------------------------------------------------------------- + + area = 0.5_RKIND*(x(2)*y(3) - y(2)*x(3) - x(1)*y(3) + y(1)*x(3) + x(1)*y(2) - y(1)*x(2)) + + !-------------------------------------------------------------------- + + end subroutine ocn_init_tri_area!}}} + +!*********************************************************************** +! +! routine ocn_init_tri_coordinate_transform +! +!> \brief Transform reference triangle coordinates to mesh coordinates +!> \author Steven Brus +!> \date November 2022 +!> \details Evaluate r,s reference coordinates in x,y mesh +!> +! +!----------------------------------------------------------------------- + + subroutine ocn_init_tri_coordinate_transform(r, s, xTri, yTri, x, y)!{{{ + + implicit none + + real (kind=RKIND), intent(in) :: r, s + real (kind=RKIND), intent(in) :: xTri(3), yTri(3) + real (kind=RKIND), intent(out) :: x, y + + !-------------------------------------------------------------------- + + ! Transformation of sub-triangle center to physical coordinates + x = 0.5_RKIND*(-(r+s)*xTri(1) + (1.0_RKIND+r)*xTri(2) + (1.0_RKIND+s)*xTri(3)) + y = 0.5_RKIND*(-(r+s)*yTri(1) + (1.0_RKIND+r)*yTri(2) + (1.0_RKIND+s)*yTri(3)) + + !-------------------------------------------------------------------- + + end subroutine ocn_init_tri_coordinate_transform!}}} + +!*********************************************************************** +! +! routine ocn_init_vertical_integration +! +!> \brief Compute the wet volume per unit area lookup table +!> \author Steven Brus +!> \date November 2022 +!> \details Integrate the wet fraction over discrete ssh values +!> +! +!----------------------------------------------------------------------- + + subroutine ocn_init_vertical_integration(iCV, subgridSshTableRange, nSubgridCV, subgridBathymetryValues, subgridAreas, & + subgridWetVolumeTable, subgridWetFractionTable)!{{{ + + implicit none + + integer, intent(in) :: iCV + real (kind=RKIND), dimension(:,:), intent(in) :: subgridSshTableRange + integer, intent(in) :: nSubgridCV + real (kind=RKIND), dimension(:), intent(in) :: subgridBathymetryValues, subgridAreas + real (kind=RKIND), dimension(:,:), intent(inout) :: subgridWetVolumeTable + real (kind=RKIND), dimension(:,:), intent(inout) :: subgridWetFractionTable + + real (kind=RKIND) :: deltaZ, ssh, pVal + integer :: lev, tri + + + !-------------------------------------------------------------------- + + deltaZ = (subgridSshTableRange(2,iCV)-subgridSshTableRange(1,iCV))/real(nSubgridTableLevels-1,RKIND); + + subgridWetVolumeTable(1,iCV) = config_drying_min_cell_height + eps + ssh = subgridSshTableRange(1,iCV) + deltaZ + subgridWetVolumeTable(2,iCV) = subgridWetVolumeTable(1,iCV)*sum(subgridAreas(1:nSubgridCV)) + subgridWetFractionTable(1,iCV) = 0.0_RKIND + do lev = 2,nSubgridTableLevels + do tri = 1,nSubgridCV + + pVal = 0.0_RKIND + if (subgridBathymetryValues(tri) + ssh >= 0.0_RKIND) then + pVal = 1.0_RKIND + endif + + subgridWetVolumeTable(lev,iCV) = subgridWetVolumeTable(lev,iCV) + pVal*deltaZ*subgridAreas(tri) + subgridWetFractionTable(lev,iCV) = subgridWetFractionTable(lev,iCV) + pVal*subgridAreas(tri) + + enddo + + if (lev < nSubgridTableLevels) then + subgridWetVolumeTable(lev+1,iCV) = subgridWetVolumeTable(lev,iCV) + endif + + subgridWetVolumeTable(lev,iCV) = subgridWetVolumeTable(lev,iCV)/sum(subgridAreas(1:nSubgridCV)) + subgridWetFractionTable(lev,iCV) = subgridWetFractionTable(lev,iCV)/sum(subgridAreas(1:nSubgridCV)) + + ssh = ssh + deltaZ + enddo + + !-------------------------------------------------------------------- + + end subroutine ocn_init_vertical_integration!}}} + + +!*********************************************************************** +! +! routine ocn_init_wet_average +! +!> \brief Compute thickness and velocity averages over wet area +!> \author Steven Brus +!> \date November 2022 +!> \details Compute the thickness and velocity averages over an area +!> based on the subgrid wet area +! +! DW: Shall we rename it to 'ocn_init_grid_average'? +!----------------------------------------------------------------------- + + subroutine ocn_init_wet_average(nSubgridCV, subgridBathymetryValues, subgridSshValues, subgridAreas, subgridThicknessAverage, & + subgridUValues, subgridVValues, subgridUAverage, subgridVAverage)!{{{ + + implicit none + + integer, intent(in) :: nSubgridCV + real (kind=RKIND), dimension(:), intent(in) :: subgridBathymetryValues, subgridAreas + real (kind=RKIND), dimension(:), intent(in) :: subgridSshValues + real (kind=RKIND), intent(inout) :: subgridThicknessAverage + real (kind=RKIND), dimension(:), intent(in), optional :: subgridUValues, subgridVValues + real (kind=RKIND), intent(inout), optional :: subgridUAverage, subgridVAverage + + real (kind=RKIND) :: deltaZ, ssh, pVal + real (kind=RKIND) :: averageDepth, wetArea, layerThicknessValue + real (kind=RKIND) :: H_int + real (kind=RKIND) :: HU_int, HV_int + logical :: computeVelAverage + integer :: lev, tri + + + !-------------------------------------------------------------------- + + averageDepth = sum(subgridBathymetryValues(1:nSubgridCV))/real(nSubgridCV,RKIND) + + computeVelAverage = .false. + if (present(subgridUValues) .and. present(subgridVValues) .and. & + present(subgridUAverage) .and. present(subgridVAverage)) then + computeVelAverage = .true. + endif + + H_int = 0.0_RKIND + wetArea = 0.0_RKIND + HU_int = 0.0_RKIND + HV_int = 0.0_RKIND + do tri = 1,nSubgridCV + + layerThicknessValue = subgridBathymetryValues(tri) + subgridSshValues(tri) + if (layerThicknessValue < config_drying_min_cell_height + eps) then + layerThicknessValue = config_drying_min_cell_height + eps + endif + H_int = H_int + layerThicknessValue*subgridAreas(tri) + wetArea = wetArea + subgridAreas(tri) + + if (computeVelAverage) then + HU_int = HU_int + layerThicknessValue*subgridUValues(tri)*subgridAreas(tri) + HV_int = HV_int + layerThicknessValue*subgridVValues(tri)*subgridAreas(tri) + endif + enddo + + if (computeVelAverage) then + subgridUAverage = HU_int/H_int + subgridVAverage = HV_int/H_int + endif + subgridThicknessAverage = H_int/wetArea + + + !-------------------------------------------------------------------- + + end subroutine ocn_init_wet_average!}}} + + ! + subroutine ocn_init_wet_average_ssh( nSubgridCV, subgridBathymetryValues, subgridSshValues, subgridAreas, sshWetAverage)!{{{ + implicit none + + integer, intent(in) :: nSubgridCV + real (kind=RKIND), dimension(:), intent(in) :: subgridBathymetryValues, subgridAreas + real (kind=RKIND), dimension(:), intent(in) :: subgridSshValues + real (kind=RKIND), intent(inout) :: sshWetAverage + + real (kind=RKIND) :: deltaZ, ssh, pVal + real (kind=RKIND) :: averageDepth, wetArea, layerThicknessValue + integer :: lev, tri + + + sshWetAverage = 0.0_RKIND + wetArea = 0.0_RKIND + + do tri = 1, nsubgridCV + layerThicknessValue = subgridBathymetryValues(tri) + subgridSshValues(tri) + + if (layerThicknessValue > config_drying_min_cell_height + eps ) then + sshWetAverage = sshWetAverage + subgridSshValues(tri)*subgridAreas(tri) + wetArea = wetArea + subgridAreas(tri) + endif + end do + + if ( WetArea > 0.0_RKIND ) then + sshWetAverage = sshWetAverage/WetArea ; + else + sshWetAverage = -maxval(subgridBathymetryValues(1:nsubgridCV)) & + + config_drying_min_cell_height ; + endif + + return ; + end subroutine ocn_init_wet_average_ssh + +!*********************************************************************** +! +! routine ocn_init_setup_global_ocean_read_topo +! +!> \brief Read the topography IC file +!> \author D. Wirasaet, S. Brus +!> \date Sep 2023 +!> \details +!> This routine reads the topography IC file, including latitude and longitude +!> information for topography data. +!> +!> Adapted from ocn_init_setup_global_ocean_read_topo +!----------------------------------------------------------------------- + subroutine ocn_init_setup_Buttermilk_bay_read_topo(domain, iErr)!{{{ + type (domain_type), intent(inout) :: domain + integer, intent(out) :: iErr + + type (MPAS_Stream_type) :: topographyStream + + iErr = 0 + + ! Define stream for depth levels + call MPAS_createStream(topographyStream, domain % iocontext, & + config_Buttermilk_bay_topography_file, MPAS_IO_NETCDF, & + MPAS_IO_READ, ierr=iErr) + + ! Setup topoLat, topoLon, and topoIC fields for stream to be read in + topoLat % fieldName = trim(config_Buttermilk_bay_topography_lat_varname) + topoLat % dimSizes(1) = nLatTopo + topoLat % dimNames(1) = trim(config_Buttermilk_bay_topography_nlat_dimname) + topoLat % isVarArray = .false. + topoLat % isPersistent = .true. + topoLat % isActive = .true. + topoLat % hasTimeDimension = .false. + topoLat % block => domain % blocklist + allocate(topoLat % attLists(1)) + allocate(topoLat % array(nLatTopo)) + + topoLon % fieldName = trim(config_Buttermilk_bay_topography_lon_varname) + topoLon % dimSizes(1) = nLonTopo + topoLon % dimNames(1) = trim(config_Buttermilk_bay_topography_nlon_dimname) + topoLon % isVarArray = .false. + topoLon % isPersistent = .true. + topoLon % isActive = .true. + topoLon % hasTimeDimension = .false. + topoLon % block => domain % blocklist + allocate(topoLon % attLists(1)) + allocate(topoLon % array(nLonTopo)) + + topoIC % fieldName = trim(config_Buttermilk_bay_topography_varname) + topoIC % dimSizes(1) = nLonTopo + topoIC % dimSizes(2) = nLatTopo + topoIC % dimNames(1) = trim(config_Buttermilk_bay_topography_nlon_dimname) + topoIC % dimNames(2) = trim(config_Buttermilk_bay_topography_nlat_dimname) + topoIC % isVarArray = .false. + topoIC % isPersistent = .true. + topoIC % isActive = .true. + topoIC % hasTimeDimension = .false. + topoIC % block => domain % blocklist + allocate(topoIC % attLists(1)) + allocate(topoIC % array(nLonTopo, nLatTopo)) + + ! Add topoLat, topoLon, and topoIC fields to stream + call MPAS_streamAddField(topographyStream, topoLat, iErr) + call MPAS_streamAddField(topographyStream, topoLon, iErr) + call MPAS_streamAddField(topographyStream, topoIC, iErr) + + + ! Read stream + call MPAS_readStream(topographyStream, 1, iErr) + + ! Close stream + call MPAS_closeStream(topographyStream) + + + end subroutine ocn_init_setup_Buttermilk_bay_read_topo!}}} + + +!*********************************************************************** +! +! routine ocn_init_Buttermilk_bay_destroy_topo_fields +! +!> \brief Topography field cleanup routine +!> \author D. Wirasaet and S. Brus +!> \date Sep 2023 +!> \details +!> This routine destroys the fields that were created to hold topography +!> initial condition information +!> +!> NOTE: adapteed from ocn_init_global_ocaen_destroy_topo_fileds +!----------------------------------------------------------------------- + + subroutine ocn_init_Buttermilk_bay_destroy_topo_fields()!{{{ + implicit none + + deallocate(topoIC % array) + deallocate(topoLat % array) + deallocate(topoLon % array) + end subroutine ocn_init_Buttermilk_bay_destroy_topo_fields!}}} + + +!*********************************************************************** +! +! routine ocn_init_Buttermilk_bay +! +!> \brief Validation for this initial condition +!> \author D. Wirasaet and Steven Brus +!> \date Sep 2022 +!> \details +!> This routine validates the configuration options for this case. +!> +!----------------------------------------------------------------------- + subroutine ocn_init_validate_Buttermilk_bay(configPool, packagePool, iocontext, iErr)!{{{ + implicit none + + !-------------------------------------------------------------------- + type (mpas_pool_type), intent(inout) :: configPool, packagePool + type (mpas_io_context_type), intent(inout) :: iocontext + + integer, intent(out) :: iErr + + ! character (len=StrKIND), pointer :: config_init_configuration + integer, pointer :: config_vert_levels, config_Buttermilk_bay_vert_levels + + + type (mpas_io_context_type), pointer :: iocontext_ptr + type (MPAS_IO_Handle_type) :: inputFile + character (len=StrKIND), pointer :: config_init_configuration, & + config_Buttermilk_bay_topography_source, & + config_Buttermilk_bay_topography_file, & + config_Buttermilk_bay_topography_nlat_dimname, & + config_Buttermilk_bay_topography_nlon_dimname + + iErr = 0 + + call mpas_pool_get_config(configPool, 'config_init_configuration', config_init_configuration) + + if(config_init_configuration .ne. trim('buttermilk_bay')) return + + call mpas_pool_get_config(configPool, 'config_vert_levels', config_vert_levels) + call mpas_pool_get_config(configPool, & + & 'config_Buttermilk_bay_vert_levels', config_Buttermilk_bay_vert_levels) + + if(config_vert_levels <= 0 .and. config_Buttermilk_bay_vert_levels > 0) then + config_vert_levels = config_Buttermilk_bay_vert_levels + else if (config_vert_levels <= 0) then + call mpas_log_write( 'Validation failed for Buttermilk bay.'// & + & 'Not given a usable value for vertical levels.', MPAS_LOG_CRIT) + iErr = 1 + end if + + + !---------- adapted from ocn_init_validate_global_ocean + call mpas_pool_get_config(configPool, 'config_Buttermilk_bay_topography_source', & + config_Buttermilk_bay_topography_source) + call mpas_pool_get_config(configPool, 'config_Buttermilk_bay_topography_file', & + config_Buttermilk_bay_topography_file) + call mpas_pool_get_config(configPool, 'config_Buttermilk_bay_topography_nlat_dimname', & + config_Buttermilk_bay_topography_nlat_dimname) + call mpas_pool_get_config(configPool, 'config_Buttermilk_bay_topography_nlon_dimname', & + config_Buttermilk_bay_topography_nlon_dimname) + + + if (config_Buttermilk_bay_topography_source /= 'latlon_file' .and. & + config_Buttermilk_bay_topography_source /= 'xy_file') then + call mpas_log_write( 'Unexpected value for & + & config_Buttermilk_bay_topography_source: ' & + // trim(config_Buttermilk_bay_topography_source), MPAS_LOG_CRIT) + iErr = 1 + return + end if + + ! print*, "message ", config_Buttermilk_bay_topography_file, config_Buttermilk_bay_topography_source + if (config_Buttermilk_bay_topography_file == 'none' .and. & + (config_Buttermilk_bay_topography_source == 'latlon_file' .or. & + config_Buttermilk_bay_topography_source == 'xy_file') ) then + call mpas_log_write( 'Validation failed for Buttermilk bay test case. ' & + // 'Invalid filename for config_Buttermilk_bay_topography_file ' & + // config_Buttermilk_bay_topography_file // ' ' & + // config_Buttermilk_bay_topography_source , MPAS_LOG_CRIT) + iErr = 1 + return + end if + +!! call mpas_log_write( ' in ocn_init_validate_buttermilk_bay '//config_Buttermilk_bay_topography_source, MPAS_LOG_CRIT ) + + if (config_Buttermilk_bay_topography_source == 'latlon_file' .or. & + config_Buttermilk_bay_topography_source == 'xy_file' ) then + +!! call mpas_log_write( ' before io_open', MPAS_LOG_CRIT ) + + inputFile = MPAS_io_open(config_Buttermilk_bay_topography_file, & + MPAS_IO_READ, MPAS_IO_NETCDF, iocontext_ptr, ierr=iErr) + if (iErr /= 0) then + call mpas_log_write( 'could not open file '// & + & trim(config_Buttermilk_bay_topography_file), MPAS_LOG_CRIT) + return + end if + + call mpas_log_write( 'after io_open', MPAS_LOG_CRIT ) + + + call MPAS_io_inq_dim(inputFile, config_Buttermilk_bay_topography_nlat_dimname, nLatTopo, iErr) + call MPAS_io_inq_dim(inputFile, config_Buttermilk_bay_topography_nlon_dimname, nLonTopo, iErr) + + call MPAS_io_close(inputFile, iErr) + + call mpas_log_write( ' in ocn_init_validate_buttermilk_bay & + after io_close', MPAS_LOG_CRIT ) + + if ( config_Buttermilk_bay_topography_latlon_degrees .and. & + config_Buttermilk_bay_topography_source == 'latlon_file' ) then + topoLat % array(:) = topoLat % array(:) * pii / 180.0_RKIND + topoLon % array(:) = topoLon % array(:) * pii / 180.0_RKIND + end if + + end if + !---------- + + call mpas_log_write( ' Done ocn_init_validate_buttermilk_bay ', MPAS_LOG_CRIT ) + !-------------------------------------------------------------------- + + end subroutine ocn_init_validate_ButterMilk_bay!}}} + + +!!*********************************************************************** +!! +!! routine ocn_init_validate_global_ocean +!! +!!> \brief Validation for global ocean test case +!!> \author Doug Jacobsen +!!> \date 03/04/2014 +!!> \details +!!> This routine validates the configuration options for the global ocean test case. +!! +!!----------------------------------------------------------------------- +! +! subroutine ocn_init_validate_global_ocean(configPool, packagePool, iocontext, iErr)!{{{ +! +! !-------------------------------------------------------------------- +! +! type (mpas_pool_type), intent(inout) :: configPool, packagePool +! type (mpas_io_context_type), intent(inout), target :: iocontext +! integer, intent(out) :: iErr +! +! type (mpas_io_context_type), pointer :: iocontext_ptr +! type (MPAS_IO_Handle_type) :: inputFile +! character (len=StrKIND), pointer :: config_init_configuration, & +! config_global_ocean_topography_source, & +! config_global_ocean_depth_file, & +! config_global_ocean_depth_dimname, & +! config_global_ocean_temperature_file, & +! config_global_ocean_salinity_file, & +! config_global_ocean_tracer_nlat_dimname, & +! config_global_ocean_tracer_nlon_dimname, & +! config_global_ocean_tracer_ndepth_dimname, & +! config_global_ocean_topography_file, & +! config_global_ocean_topography_nlat_dimname, & +! config_global_ocean_topography_nlon_dimname, & +! config_global_ocean_windstress_file, & +! config_global_ocean_windstress_nlat_dimname, & +! config_global_ocean_windstress_nlon_dimname, & +! config_global_ocean_land_ice_topo_file, & +! config_global_ocean_land_ice_topo_nlat_dimname, & +! config_global_ocean_land_ice_topo_nlon_dimname, & +! config_global_ocean_swData_file, & +! config_global_ocean_swData_nlon_dimname, & +! config_global_ocean_swData_nlat_dimname, & +! config_global_ocean_ecosys_file, & +! config_global_ocean_ecosys_nlat_dimname, & +! config_global_ocean_ecosys_nlon_dimname, & +! config_global_ocean_ecosys_ndepth_dimname +! +! integer, pointer :: config_vert_levels, config_global_ocean_tracer_vert_levels, & +! config_global_ocean_ecosys_vert_levels +! logical, pointer :: config_use_ecosysTracers +! logical, pointer :: landIceInitActive, config_global_ocean_depress_by_land_ice +! logical, pointer :: criticalPassagesActive, config_global_ocean_deepen_critical_passages +! +! iocontext_ptr => iocontext +! iErr = 0 +! +! call mpas_pool_get_config(configPool, 'config_init_configuration', config_init_configuration) +! +! if(config_init_configuration .ne. trim('global_ocean')) return +! +! call mpas_pool_get_config(configPool, 'config_vert_levels', config_vert_levels) +! +! call mpas_pool_get_config(configPool, 'config_global_ocean_topography_source', & +! config_global_ocean_topography_source) +! call mpas_pool_get_config(configPool, 'config_global_ocean_depth_file', & +! config_global_ocean_depth_file) +! call mpas_pool_get_config(configPool, 'config_global_ocean_depth_dimname', & +! config_global_ocean_depth_dimname) +! call mpas_pool_get_config(configPool, 'config_global_ocean_temperature_file', & +! config_global_ocean_temperature_file) +! call mpas_pool_get_config(configPool, 'config_global_ocean_salinity_file', & +! config_global_ocean_salinity_file) +! call mpas_pool_get_config(configPool, 'config_global_ocean_tracer_nlat_dimname', & +! config_global_ocean_tracer_nlat_dimname) +! call mpas_pool_get_config(configPool, 'config_global_ocean_tracer_nlon_dimname', & +! config_global_ocean_tracer_nlon_dimname) +! call mpas_pool_get_config(configPool, 'config_global_ocean_tracer_ndepth_dimname', & +! config_global_ocean_tracer_ndepth_dimname) +! call mpas_pool_get_config(configPool, 'config_global_ocean_tracer_vert_levels', & +! config_global_ocean_tracer_vert_levels) +! call mpas_pool_get_config(configPool, 'config_global_ocean_topography_file', & +! config_global_ocean_topography_file) +! call mpas_pool_get_config(configPool, 'config_global_ocean_topography_nlat_dimname', & +! config_global_ocean_topography_nlat_dimname) +! call mpas_pool_get_config(configPool, 'config_global_ocean_topography_nlon_dimname', & +! config_global_ocean_topography_nlon_dimname) +! call mpas_pool_get_config(configPool, 'config_global_ocean_depress_by_land_ice', & +! config_global_ocean_depress_by_land_ice) +! call mpas_pool_get_config(configPool, 'config_use_ecosysTracers', & +! config_use_ecosysTracers) +! call mpas_pool_get_config(configPool, 'config_global_ocean_land_ice_topo_file', & +! config_global_ocean_land_ice_topo_file) +! call mpas_pool_get_config(configPool, 'config_global_ocean_land_ice_topo_nlat_dimname', & +! config_global_ocean_land_ice_topo_nlat_dimname) +! call mpas_pool_get_config(configPool, 'config_global_ocean_land_ice_topo_nlon_dimname', & +! config_global_ocean_land_ice_topo_nlon_dimname) +! call mpas_pool_get_config(configPool, 'config_global_ocean_deepen_critical_passages', & +! config_global_ocean_deepen_critical_passages) +! call mpas_pool_get_config(configPool, 'config_global_ocean_windstress_file', & +! config_global_ocean_windstress_file) +! call mpas_pool_get_config(configPool, 'config_global_ocean_windstress_nlat_dimname', & +! config_global_ocean_windstress_nlat_dimname) +! call mpas_pool_get_config(configPool, 'config_global_ocean_windstress_nlon_dimname', & +! config_global_ocean_windstress_nlon_dimname) +! call mpas_pool_get_config(configPool, 'config_global_ocean_swData_file', & +! config_global_ocean_swData_file) +! call mpas_pool_get_config(configPool, 'config_global_ocean_swData_nlat_dimname', & +! config_global_ocean_swData_nlat_dimname) +! call mpas_pool_get_config(configPool, 'config_global_ocean_swData_nlon_dimname', & +! config_global_ocean_swData_nlon_dimname) +! call mpas_pool_get_config(configPool, 'config_global_ocean_ecosys_file', & +! config_global_ocean_ecosys_file) +! call mpas_pool_get_config(configPool, 'config_global_ocean_ecosys_nlat_dimname', & +! config_global_ocean_ecosys_nlat_dimname) +! call mpas_pool_get_config(configPool, 'config_global_ocean_ecosys_nlon_dimname', & +! config_global_ocean_ecosys_nlon_dimname) +! call mpas_pool_get_config(configPool, 'config_global_ocean_ecosys_ndepth_dimname', & +! config_global_ocean_ecosys_ndepth_dimname) +! call mpas_pool_get_config(configPool, 'config_global_ocean_ecosys_vert_levels', & +! config_global_ocean_ecosys_vert_levels) +! +! call mpas_pool_get_package(packagePool, 'landIceInitActive', landIceInitActive) +! if ( config_global_ocean_depress_by_land_ice) then +! landIceInitActive = .true. +! end if +! +! call mpas_pool_get_package(packagePool, 'criticalPassagesActive', criticalPassagesActive) +! if ( config_global_ocean_deepen_critical_passages) then +! criticalPassagesActive = .true. +! end if +! +! if (trim(config_global_ocean_depth_file) == 'none') then +! call mpas_log_write( 'Validation failed for global ocean. ' & +! // 'Invalid filename for config_global_ocean_depth_file', MPAS_LOG_CRIT) +! iErr = 1 +! return +! end if +! +! inputFile = MPAS_io_open(config_global_ocean_depth_file, MPAS_IO_READ, MPAS_IO_NETCDF, iocontext_ptr, ierr=iErr) +! if (iErr .ne. 0) then +! call mpas_log_write( 'could not open file '// trim(config_global_ocean_depth_file), MPAS_LOG_CRIT) +! return +! end if +! +! call MPAS_io_inq_dim(inputFile, config_global_ocean_depth_dimname, nDepthOutput, iErr) +! +! call MPAS_io_close(inputFile, iErr) +! +! if (trim(config_global_ocean_temperature_file) == 'none') then +! call mpas_log_write( 'Validation failed for global ocean. ' & +! // 'Invalid filename for config_global_ocean_temperature_file', MPAS_LOG_CRIT) +! iErr = 1 +! return +! end if +! +! if (trim(config_global_ocean_salinity_file) == 'none') then +! call mpas_log_write( 'Validation failed for global ocean. ' & +! // 'Invalid filename for config_global_ocean_salinity_file', MPAS_LOG_CRIT) +! iErr = 1 +! return +! end if +! +! inputFile = MPAS_io_open(config_global_ocean_temperature_file, MPAS_IO_READ, MPAS_IO_NETCDF, iocontext_ptr, ierr=iErr) +! if (iErr .ne. 0) then +! call mpas_log_write( 'could not open file '// trim(config_global_ocean_temperature_file), MPAS_LOG_CRIT) +! return +! end if +! +! call MPAS_io_inq_dim(inputFile, config_global_ocean_tracer_nlat_dimname, nLatTracer, iErr) +! call MPAS_io_inq_dim(inputFile, config_global_ocean_tracer_nlon_dimname, nLonTracer, iErr) +! call MPAS_io_inq_dim(inputFile, config_global_ocean_tracer_ndepth_dimname, nDepthTracer, iErr) +! +! call MPAS_io_close(inputFile, iErr) +! +! if (config_global_ocean_tracer_vert_levels <= 0 .and. nDepthTracer > 0) then +! config_global_ocean_tracer_vert_levels = nDepthTracer +! else if(config_global_ocean_tracer_vert_levels <= 0) then +! call mpas_log_write( 'Validation failed for global ocean. ' & +! // 'Value of config_global_ocean_tracer_vert_levels=-1 ' & +! // 'but nDepthTracer was not correctly read from input file.', MPAS_LOG_CRIT) +! iErr = 1 +! end if +! +! if (trim(config_global_ocean_windstress_file) == 'none') then +! call mpas_log_write( 'Validation failed for global ocean. ' & +! // 'Invalid filename for config_global_ocean_windstress_file', MPAS_LOG_CRIT) +! iErr = 1 +! return +! end if +! +! inputFile = MPAS_io_open(config_global_ocean_swData_file, MPAS_IO_READ, MPAS_IO_NETCDF, iocontext_ptr, ierr=iErr) +! +! call MPAS_io_inq_dim(inputFile, config_global_ocean_swData_nlat_dimname, nLatSW, iErr) +! call MPAS_io_inq_dim(inputFile, config_global_ocean_swData_nlon_dimname, nLonSW, iErr) +! +! call MPAS_io_close(inputFile, iErr) +! +! if (config_global_ocean_topography_source /= 'latlon_file' .and. & +! config_global_ocean_topography_source /= 'mpas_variable') then +! call mpas_log_write( 'Unexpected value for config_global_ocean_topography_source: ' & +! // trim(config_global_ocean_topography_source), MPAS_LOG_CRIT) +! iErr = 1 +! return +! end if +! +! if (config_global_ocean_topography_file == 'none' .and. & +! config_global_ocean_topography_source == 'latlon_file') then +! call mpas_log_write( 'Validation failed for global ocean. ' & +! // 'Invalid filename for config_global_ocean_topography_file', MPAS_LOG_CRIT) +! iErr = 1 +! return +! end if +! +! if (config_global_ocean_topography_source == 'latlon_file') then +! inputFile = MPAS_io_open(config_global_ocean_topography_file, MPAS_IO_READ, MPAS_IO_NETCDF, iocontext_ptr, ierr=iErr) +! if (iErr /= 0) then +! call mpas_log_write( 'could not open file '// trim(config_global_ocean_topography_file), MPAS_LOG_CRIT) +! return +! end if +! +! call MPAS_io_inq_dim(inputFile, config_global_ocean_topography_nlat_dimname, nLatTopo, iErr) +! call MPAS_io_inq_dim(inputFile, config_global_ocean_topography_nlon_dimname, nLonTopo, iErr) +! +! call MPAS_io_close(inputFile, iErr) +! end if +! +! inputFile = MPAS_io_open(config_global_ocean_windstress_file, MPAS_IO_READ, MPAS_IO_NETCDF, iocontext_ptr, ierr=iErr) +! if (iErr .ne. 0) then +! call mpas_log_write( 'could not open file '// trim(config_global_ocean_windstress_file), MPAS_LOG_CRIT) +! return +! end if +! +! call MPAS_io_inq_dim(inputFile, config_global_ocean_windstress_nlat_dimname, nLatWind, iErr) +! call MPAS_io_inq_dim(inputFile, config_global_ocean_windstress_nlon_dimname, nLonWind, iErr) +! +! call MPAS_io_close(inputFile, iErr) +! +! if (config_vert_levels <= 0 .and. nDepthOutput > 0) then +! config_vert_levels = nDepthOutput +! else if(config_vert_levels <= 0) then +! call mpas_log_write( 'Validation failed for global ocean. Not given a usable value for vertical levels.', MPAS_LOG_CRIT) +! iErr = 1 +! end if +! +! if ( config_use_ecosysTracers ) then +! if (trim(config_global_ocean_ecosys_file) == 'none') then +! call mpas_log_write( & +! 'Validation failed for global ocean. Invalid filename for config_global_ocean_windstress_file', MPAS_LOG_CRIT) +! iErr = 1 +! return +! end if +! +! inputFile = MPAS_io_open(config_global_ocean_ecosys_file, MPAS_IO_READ, MPAS_IO_NETCDF, iocontext_ptr, ierr=iErr) +! +! call MPAS_io_inq_dim(inputFile, config_global_ocean_ecosys_nlat_dimname, nLatEcosys, iErr) +! call MPAS_io_inq_dim(inputFile, config_global_ocean_ecosys_nlon_dimname, nLonEcosys, iErr) +! call MPAS_io_inq_dim(inputFile, config_global_ocean_ecosys_ndepth_dimname, nDepthEcosys, iErr) +! +! call MPAS_io_close(inputFile, iErr) +! +! if (config_global_ocean_ecosys_vert_levels <= 0 .and. nDepthEcosys > 0) then +! config_global_ocean_ecosys_vert_levels = nDepthEcosys +! else if(config_global_ocean_ecosys_vert_levels <= 0) then +! call mpas_log_write( 'Validation failed for global ocean. ' & +! // 'Value of config_global_ocean_ecosys_vert_levels=-1, ' & +! // 'but nDepthEcosys was not correctly read from input file.', MPAS_LOG_CRIT) +! iErr = 1 +! end if +! +! end if +! +! if (config_global_ocean_depress_by_land_ice) then +! if (config_global_ocean_land_ice_topo_file == 'none' .and. & +! config_global_ocean_topography_source == 'latlon_file') then +! call mpas_log_write( 'Validation failed for global ocean. '// & +! 'Invalid filename for config_global_ocean_land_ice_topo_file', MPAS_LOG_CRIT) +! iErr = 1 +! return +! end if +! +! if(config_global_ocean_topography_source == 'latlon_file') then +! inputFile = MPAS_io_open(config_global_ocean_land_ice_topo_file, MPAS_IO_READ, MPAS_IO_NETCDF, iocontext_ptr, ierr=iErr) +! if (iErr /= 0) then +! call mpas_log_write( 'could not open file '// trim(config_global_ocean_land_ice_topo_file), MPAS_LOG_CRIT) +! return +! end if +! +! call MPAS_io_inq_dim(inputFile, config_global_ocean_land_ice_topo_nlat_dimname, nLatLandIceThk, iErr) +! call MPAS_io_inq_dim(inputFile, config_global_ocean_land_ice_topo_nlon_dimname, nLonLandIceThk, iErr) +! +! call MPAS_io_close(inputFile, iErr) +! end if +! end if +! +! !-------------------------------------------------------------------- +! +! end subroutine ocn_init_validate_global_ocean!}}} +! + + + +!!*********************************************************************** +!! +!! routine ocn_init_setup_global_ocean_read_topo +!! +!!> \brief Read the topography IC file +!!> \author Doug Jacobsen +!!> \date 03/04/2014 +!!> \details +!!> This routine reads the topography IC file, including latitude and longitude +!!> information for topography data. +!! +!!----------------------------------------------------------------------- +! +! subroutine ocn_init_setup_global_ocean_read_topo(domain, iErr)!{{{ +! type (domain_type), intent(inout) :: domain +! integer, intent(out) :: iErr +! +! type (MPAS_Stream_type) :: topographyStream +! +! iErr = 0 +! +! ! Define stream for depth levels +! call MPAS_createStream(topographyStream, domain % iocontext, config_global_ocean_topography_file, MPAS_IO_NETCDF, & +! MPAS_IO_READ, ierr=iErr) +! +! ! Setup topoLat, topoLon, and topoIC fields for stream to be read in +! topoLat % fieldName = trim(config_global_ocean_topography_lat_varname) +! topoLat % dimSizes(1) = nLatTopo +! topoLat % dimNames(1) = trim(config_global_ocean_topography_nlat_dimname) +! topoLat % isVarArray = .false. +! topoLat % isPersistent = .true. +! topoLat % isActive = .true. +! topoLat % hasTimeDimension = .false. +! topoLat % block => domain % blocklist +! allocate(topoLat % attLists(1)) +! allocate(topoLat % array(nLatTopo)) +! +! topoLon % fieldName = trim(config_global_ocean_topography_lon_varname) +! topoLon % dimSizes(1) = nLonTopo +! topoLon % dimNames(1) = trim(config_global_ocean_topography_nlon_dimname) +! topoLon % isVarArray = .false. +! topoLon % isPersistent = .true. +! topoLon % isActive = .true. +! topoLon % hasTimeDimension = .false. +! topoLon % block => domain % blocklist +! allocate(topoLon % attLists(1)) +! allocate(topoLon % array(nLonTopo)) +! +! topoIC % fieldName = trim(config_global_ocean_topography_varname) +! topoIC % dimSizes(1) = nLonTopo +! topoIC % dimSizes(2) = nLatTopo +! topoIC % dimNames(1) = trim(config_global_ocean_topography_nlon_dimname) +! topoIC % dimNames(2) = trim(config_global_ocean_topography_nlat_dimname) +! topoIC % isVarArray = .false. +! topoIC % isPersistent = .true. +! topoIC % isActive = .true. +! topoIC % hasTimeDimension = .false. +! topoIC % block => domain % blocklist +! allocate(topoIC % attLists(1)) +! allocate(topoIC % array(nLonTopo, nLatTopo)) +! +! ! Add topoLat, topoLon, and topoIC fields to stream +! call MPAS_streamAddField(topographyStream, topoLat, iErr) +! call MPAS_streamAddField(topographyStream, topoLon, iErr) +! call MPAS_streamAddField(topographyStream, topoIC, iErr) +! +! if(config_global_ocean_topography_has_ocean_frac) then +! oceanFracIC % fieldName = trim(config_global_ocean_topography_ocean_frac_varname) +! oceanFracIC % dimSizes(1) = nLonTopo +! oceanFracIC % dimSizes(2) = nLatTopo +! oceanFracIC % dimNames(1) = trim(config_global_ocean_topography_nlon_dimname) +! oceanFracIC % dimNames(2) = trim(config_global_ocean_topography_nlat_dimname) +! oceanFracIC % isVarArray = .false. +! oceanFracIC % isPersistent = .true. +! oceanFracIC % isActive = .true. +! oceanFracIC % hasTimeDimension = .false. +! oceanFracIC % block => domain % blocklist +! allocate(oceanFracIC % attLists(1)) +! allocate(oceanFracIC % array(nLonTopo, nLatTopo)) +! +! call MPAS_streamAddField(topographyStream, oceanFracIC, iErr) +! end if +! +! ! Read stream +! call MPAS_readStream(topographyStream, 1, iErr) +! +! ! Close stream +! call MPAS_closeStream(topographyStream) +! +! if (config_global_ocean_topography_latlon_degrees) then +! topoLat % array(:) = topoLat % array(:) * pii / 180.0_RKIND +! topoLon % array(:) = topoLon % array(:) * pii / 180.0_RKIND +! end if +! +! end subroutine ocn_init_setup_global_ocean_read_topo!}}} + + + +!*********************************************************************** + +end module ocn_init_Buttermilk_bay + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! vim: foldmethod=marker diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_mode.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_mode.F index 5d93ef354763..37f1cc9c58d7 100644 --- a/components/mpas-ocean/src/mode_init/mpas_ocn_init_mode.F +++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_mode.F @@ -64,8 +64,9 @@ module ocn_init_mode use ocn_init_transport_tests use ocn_init_test_sht - ! DW + use ocn_init_parabolic_bowl + use ocn_init_Buttermilk_bay implicit none private @@ -311,6 +312,7 @@ function ocn_init_mode_run(domain) result(iErr)!{{{ call ocn_init_setup_transport_tests(domain, ierr) call ocn_init_setup_test_sht(domain, ierr) call ocn_init_setup_parabolic_bowl(domain, iErr) + call ocn_init_setup_Buttermilk_bay(domain, iErr) !call ocn_init_setup_TEMPLATE(domain, ierr) call mpas_log_write( ' Completed setup of: ' // trim(config_init_configuration)) @@ -430,6 +432,8 @@ subroutine ocn_init_mode_validate_configuration(configPool, packagePool, ioconte call ocn_init_validate_parabolic_bowl(configPool, packagePool, iocontext, iErr=err_tmp) iErr = ior(iErr, err_tmp) + call ocn_init_validate_Buttermilk_bay(configPool, packagePool, iocontext, iErr=err_tmp) + iErr = ior(iErr, err_tmp) ! call ocn_init_validate_TEMPLATE(configPool, packagePool, iocontext, iErr=err_tmp) ! iErr = ior(iErr, err_tmp) diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F index d1a69a95a1f9..c69c2419bb10 100644 --- a/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F +++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F @@ -7,7 +7,7 @@ ! !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! -! ocn_init_parabolic_bowl +! ocn_init_Buttermilk_bay ! !> \brief MPAS ocean initialize case -- TEMPLATE !> \author D. Wirasaet, S. Brus From 98d6e2bd90d5aa6f9ff9168adb44f9347810b372 Mon Sep 17 00:00:00 2001 From: "D. Wirasaet" Date: Wed, 13 Dec 2023 10:15:38 -0800 Subject: [PATCH 07/16] commit buttermilk bay in progress material --- .../mpas_ocn_time_integration_rk4.F | 5 + .../mode_init/mpas_ocn_init_buttermilk_bay.F | 202 ++++++++++++------ .../mpas-ocean/src/shared/mpas_ocn_subgrid.F | 3 +- .../src/shared/mpas_ocn_tidal_forcing.F | 5 + 4 files changed, 151 insertions(+), 64 deletions(-) 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 58755c871d12..1f0cdb75a8bd 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 @@ -692,6 +692,9 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ !$omp end parallel end if +! print*, config_use_tidal_forcing +! print*, trim(config_tidal_forcing_type) + ! direct application of tidal boundary condition if (config_use_tidal_forcing .and. trim(config_tidal_forcing_type) == 'direct') then do iCell=1, nCells @@ -708,6 +711,8 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ 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 + + print*, iCell, totalDepth, layerThicknessNew(1, iCell) ; end if end do end if diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F index 36ff26239884..b46f9106d24c 100644 --- a/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F +++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F @@ -145,6 +145,7 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ type (mpas_pool_type), pointer :: statePool type (mpas_pool_type), pointer :: tracersPool type (mpas_pool_type), pointer :: verticalMeshPool + type (mpas_pool_type), pointer :: forcingPool ! local variables integer :: iCell, iEdge, iVertex, k, idx @@ -177,8 +178,10 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ real (kind=RKIND), dimension(:), pointer :: areaCell real (kind=RKIND), dimension(:,:), pointer :: edgeNormalVectors real (kind=RKIND), dimension(:,:), pointer :: normalVelocity + ! Elevation Bcs + real (kind=RKIND), dimension(:), pointer :: tidalInputMask - + real (kind=RKIND):: HH, uu, vv real (kind=RKIND):: RR, num, den real (kind=RKIND):: xshift = 0.0, yshift = 0.0 @@ -194,16 +197,15 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ integer :: c1, c2 real (kind=RKIND) :: x(3), y(3) integer :: slice, nSlice - integer :: i,j real (kind=RKIND) :: deltaZ ! DW - integer:: jj + integer:: i, j, jj integer:: nsubgridCellEdge, iEdgeSegment real (kind=RKIND):: pi real (kind=RKIND), dimension(:,:), allocatable :: cellEdgeBathymetryValues real (kind=RKIND), dimension(:), allocatable:: dsEdge real (kind=RKIND), dimension(:), allocatable:: xSubgridCell, ySubgridCell - + real (kind=RKIND):: bathymetryMin, bathymetryMax iErr = 0 call ocn_subgrid_init(domain,iErr) @@ -258,9 +260,12 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ call mpas_dmpar_max_real(domain % dminfo, xMax, xMaxGlobal) call mpas_dmpar_min_real(domain % dminfo, dcEdgeMin, dcEdgeMinGlobal) - pi = acos(-1.0_RKIND) ; - xshift = 0.0_RKIND ; - yshift = -(YMin + dcEdgeMin*sin(pi/3.0_RKIND)) ; + pi = acos(-1.0_RKIND) + + xshift = xMin + yshift = (3.0_RKIND*yMin + dcEdgeMin*sin(pi/3.0_RKIND))/3.0_RKIND + ! print*, "xMin, yMin = ", xMin, yMin, dcEdgeMin + ! print*, "xshift, yshift = ", xshift, yshift !*********************************************************************** ! @@ -281,12 +286,15 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ !-------------------------------------------------------------------- block_ptr => domain % blocklist - do while(associated(block_ptr)) +blockptr_doloop:do while(associated(block_ptr)) call mpas_pool_get_subpool(block_ptr % structs, 'mesh', meshPool) call mpas_pool_get_subpool(block_ptr % structs, 'state', statePool) call mpas_pool_get_subpool(block_ptr % structs, 'verticalMesh', verticalMeshPool) call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) + call mpas_pool_get_subpool(block_ptr % structs, 'forcing', forcingPool) + call mpas_pool_get_array(forcingPool, 'tidalInputMask', tidalInputMask) + call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) ; call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) ; call mpas_pool_get_dimension(meshPool, 'nEdgesSolve', nEdgesSolve) ; @@ -377,6 +385,19 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ xVertex = xVertex - xshift ; yVertex = yVertex - yshift ; + + ! get min, max coordinates of model domain ! + ! after adjusting the coordinates ! + yMin = min( yMin, minval(yCell(1:nCellsSolve))) + yMax = max( yMax, maxval(yCell(1:nCellsSolve))) + xMin = min( xMin, minval(xCell(1:nCellsSolve))) + xMax = max( xMax, maxval(xCell(1:nCellsSolve))) + + ! Determine global min and max values. + call mpas_dmpar_min_real(domain % dminfo, yMin, yMinGlobal) + call mpas_dmpar_max_real(domain % dminfo, yMax, yMaxGlobal) + call mpas_dmpar_min_real(domain % dminfo, xMin, xMinGlobal) + call mpas_dmpar_max_real(domain % dminfo, xMax, xMaxGlobal) end if ! Initlialze vector @@ -423,7 +444,14 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ allocate(subgridVValues(maxEdges*nSubgridTriPerSlice)) + allocate(xSubgridCell(maxEdges*nSubgridTriPerSlice)) + allocate(ySubgridCell(maxEdges*nSubgridTriPErSlice)) + xSubgridCell = 0.0_RKIND + ySubgridCell = 0.0_RKIND + open( unit = 101, file = 'BM_bathy.dat', action = 'write' ) ; + open( unit = 102, file = 'BM_ssh.dat', action = 'write' ) ; + do iCell = 1,nCellsSolve ! @@ -452,21 +480,30 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ x(3) = xVertex(v2) y(3) = yVertex(v2) - call ocn_init_evaluate_subgrid_data(x, y, nSubgridTriPerSlice, nSubgridCell, rSubgridPoints, sSubgridPoints, & - subgridBathymetryValues, subgridAreas, subgridSshValues) - + call ocn_init_evaluate_subgrid_data( x, y, nSubgridTriPerSlice, nSubgridCell, rSubgridPoints, sSubgridPoints, & + subgridBathymetryValues, subgridAreas, subgridSshValues, & + subgridUValues, subgridVValues, xSubgridCell, ysubgridCell ) !{{{ enddo do jj = 1, nSubgridCell - WRITE( 101, '(I6,2F16.8)' ) iCell, subgridBathymetryValues(jj) ; + WRITE( 101, '(I6,3F16.8)' ) iCell, xSubgridCell(jj), ysubgridCell(jj), subgridBathymetryValues(jj) ; + WRITE( 102, '(I6,3F16.8)' ) iCell, xSubgridCell(jj), ysubgridCell(jj), subgridSshValues(jj) ; end do ! Evaluate bounds of look-up table range - !--------------------------------------------------------------- + !--------------------------------------------------------------- + bathymetryMin = maxval( subgridBathymetryValues(1:nSubgridCell) ) + bathymetryMax = minval( subgridBathymetryValues(1:nSubgridCell) ) + + if ( abs(bathymetryMin - bathymetryMax) > 200.0*eps ) then + subgridSshCellTableRange(1,iCell) = -maxval(subgridBathymetryValues(1:nSubgridCell)) + config_drying_min_cell_height + eps + subgridSshCellTableRange(2,iCell) = -minval(subgridBathymetryValues(1:nSubgridCell)) + else + ! flat bathy ! + subgridSshCellTableRange(1,iCell) = -bathymetryMin + config_drying_min_cell_height + eps + subgridSshCellTableRange(2,iCell) = -bathymetryMin + 2.0*config_drying_min_cell_height + eps + endif - subgridSshCellTableRange(1,iCell) = -maxval(subgridBathymetryValues(1:nSubgridCell)) + config_drying_min_cell_height + eps - subgridSshCellTableRange(2,iCell) = -minval(subgridBathymetryValues(1:nSubgridCell)) - ! Evaluate subgrid bathymetry !--------------------------------------------------------------- ! DW @@ -481,12 +518,12 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ ! Evaluate wet layerThickness average !--------------------------------------------------------------- - call ocn_init_wet_average(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, layerThickness(1,iCell)) - +!!$ call ocn_init_wet_average(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, layerThickness(1,iCell)) -!!$ call ocn_init_wet_average_ssh(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, ssh(iCell))!{{{ + call ocn_init_wet_average_ssh(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, ssh(iCell))!{{{ enddo close( unit = 101 ) ; + close( unit = 102 ) ; !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Edges @@ -501,9 +538,6 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ allocate( dsEdge(nSubgridCellEdge) ) cellEdgeBathymetryValues = -99999 ; - allocate(xSubgridCell(maxEdges*nSubgridTriPerSlice)) - allocate(ySubgridCell(maxEdges*nSubgridTriPErSlice)) - do iEdge = 1,nEdgesSolve !-------------------------------------------------------------- @@ -704,6 +738,7 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ maxLevelCell(iCell) = nVertLevels ; ! sigma coordinates end do + open( unit = 103, file = 'ssh_cell.dat', action = 'write') do iCell = 1, nCellsSolve ! ! make sure depth is thick enough via ssh = TOTAL_DEPTH - bottomDepth @@ -712,17 +747,17 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ if (config_use_subgrid_wetting_drying) then - call ocn_subgrid_ssh_lookup(layerThickness(1,iCell),& - subgridWetVolumeCellTable(:,iCell),& - subgridSshCellTableRange(:,iCell),& - bottomDepth(iCell),& - subgridCellBathymetryMin(iCell),& - ssh(iCell)) -! call ocn_subgrid_layer_thickness_lookup(ssh(iCell), & -! subgridWetVolumeCellTable(:,iCell), & -! subgridSshCellTableRange(:,iCell),& -! bottomDepth(iCell),& -! LayerThickness(1,iCell)) +! call ocn_subgrid_ssh_lookup(layerThickness(1,iCell),& +! subgridWetVolumeCellTable(:,iCell),& +! subgridSshCellTableRange(:,iCell),& +! bottomDepth(iCell),& +! subgridCellBathymetryMin(iCell),& +! ssh(iCell)) + call ocn_subgrid_layer_thickness_lookup(ssh(iCell), & + subgridWetVolumeCellTable(:,iCell), & + subgridSshCellTableRange(:,iCell),& + bottomDepth(iCell),& + LayerThickness(1,iCell)) else call ocn_init_Buttermilk_bay_bathymetry(xCell(iCell),yCell(iCell),bottomDepth(iCell)) @@ -741,18 +776,22 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ realArgs=(/layerThickness(k,iCell)/)) end if end do - endif + endif do k = 1, maxLevelCell(iCell) restingThickness(k,iCell) = bottomDepth(iCell)/maxLevelCell(iCell) end do + + write(103,*) iCell, ssh(iCell), LayerThickness(1,iCell) ; + end do - + close( unit = 103 ) ; + end if block_ptr => block_ptr % next - end do + end do blockptr_doloop do iCell = 1,nCellsSolve call ocn_subgrid_layer_thickness_lookup(ssh(iCell), & @@ -763,6 +802,26 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ enddo ! + open( unit = 101, file = 'bcmask.dat', action = 'write' ) ; + ! Set tidal boundary mask + do iCell = 1, nCellsSolve + tidalInputMask(iCell) = 0.0_RKIND + if ( yCell(iCell) < (yMin+(dcEdgeMin*sin(pi/3.0_RKIND)/2.0_RKIND)) & + .and. yCell(iCell) > (yMin-(dcEDgeMin*sin(pi/3.0_RKIND)/2.0_RKIND)) ) then + + if ( (xCell(iCell) - dcEdgeMin/2.0_RKIND) > 2048.0_RKIND & + .and. (xCell(iCell) + dcEdgeMin/2.0_RKIND) < 3072.0_RKIND ) then + tidalInputMask(iCell) = 1.0_RKIND + end if + ! spread it over multiple cells + ! if (yCell(iCell) > (25.0e3 - 3*dcEdgeMinGlobal)) then + ! tidalInputMask(iCell) = exp(-((yCell(iCell)-25.0e3)/dcEdgeMinGlobal)**2.0) + end if + write(101,'(I10,F16.8)') iCell, tidalInputMask(iCell) ; + end do + + close( 101 ) ; + ! deallocate(interfaceLocations) if (config_global_ocean_topography_source == 'latlon_file') then @@ -948,7 +1007,7 @@ subroutine ocn_init_evaluate_subgrid_data( xTri, yTri, nSubgridTri, nSubgridCV, ! Calculate area of sub-triangle call ocn_init_tri_area(xSubgridPoints(:), ySubgridPoints(:), subgridAreas(nSubgridCV)) - ! Optionally evalulate ssh + ! Optionally evalulate ssh ! if (present(subgridsshValues)) then call ocn_init_Buttermilk_bay_ssh(xCenter, yCenter, subgridBathymetryValues(nSubgridCV), subgridSshValues(nSubgridCV)) endif @@ -1002,6 +1061,7 @@ subroutine ocn_init_Buttermilk_bay_bathymetry(x, y, depth)!{{{ dx = (/ topoLon%array(2) - topoLon%array(1), & topoLat%array(2) - topoLat%array(1) /) ; + ! Bilienar interpolation xc = (/ x, y /) ; ix = floor( (xc - x0)/dx ) + 1 ; @@ -1029,7 +1089,7 @@ subroutine ocn_init_Buttermilk_bay_bathymetry(x, y, depth)!{{{ depth = sum( bath*f ) ; else - ! nearest extrapolation + ! nearest extrapolation ! ix(1) = merge( 1, ix(1), ix(1) < 1 ) ; ix(1) = merge( nLonTopo, ix(1), ix(1) >= nLonTopo ) ; ix(2) = merge( 1, ix(2), ix(2) < 1 ) ; @@ -1037,7 +1097,7 @@ subroutine ocn_init_Buttermilk_bay_bathymetry(x, y, depth)!{{{ depth = TopoIC%array(ix(1),ix(2)) ; endif - + return ; end subroutine ocn_init_Buttermilk_bay_bathymetry!}}} @@ -1071,9 +1131,12 @@ subroutine ocn_init_Buttermilk_bay_ssh(x, y, bottomDepth, ssh)!{{{ ! !ssh = - bottomDepth + max(ssh + bottomDepth, config_drying_min_cell_height + eps) ssh = 0.0_RKIND ; - ssh = merge( -bottomDepth + eps, ssh, ssh + bottomDepth < eps ) ; + ! ssh = merge( -bottomDepth + eps, ssh, ssh + bottomDepth < eps ) ; + ! ssh = - bottomDepth + max( ssh + bottomDepth, config_drying_min_cell_height + eps ) ; + ssh = - bottomDepth + max( ssh + bottomDepth, 0.0_RKIND ) ; + return !-------------------------------------------------------------------- end subroutine ocn_init_Buttermilk_bay_ssh!}}} @@ -1397,19 +1460,26 @@ subroutine ocn_init_setup_Buttermilk_bay_read_topo(domain, iErr)!{{{ topoIC % block => domain % blocklist allocate(topoIC % attLists(1)) allocate(topoIC % array(nLonTopo, nLatTopo)) + ! Add topoLat, topoLon, and topoIC fields to stream call MPAS_streamAddField(topographyStream, topoLat, iErr) call MPAS_streamAddField(topographyStream, topoLon, iErr) call MPAS_streamAddField(topographyStream, topoIC, iErr) - - + ! Read stream call MPAS_readStream(topographyStream, 1, iErr) + topoIC%array = -topoIC%array ; + ! Close stream call MPAS_closeStream(topographyStream) + if ( config_Buttermilk_bay_topography_latlon_degrees .and. & + config_Buttermilk_bay_topography_source == 'latlon_file' ) then + topoLat % array(:) = topoLat % array(:) * pii / 180.0_RKIND + topoLon % array(:) = topoLon % array(:) * pii / 180.0_RKIND + end if end subroutine ocn_init_setup_Buttermilk_bay_read_topo!}}} @@ -1499,6 +1569,11 @@ subroutine ocn_init_validate_Buttermilk_bay(configPool, packagePool, iocontext, config_Buttermilk_bay_topography_nlon_dimname) + call mpas_log_write( config_Buttermilk_bay_topography_source ) + call mpas_log_write( config_Buttermilk_bay_topography_file ) + call mpas_log_write( config_Buttermilk_bay_topography_nlat_dimname ) + call mpas_log_write( config_Buttermilk_bay_topography_nlon_dimname ) + if (config_Buttermilk_bay_topography_source /= 'latlon_file' .and. & config_Buttermilk_bay_topography_source /= 'xy_file') then call mpas_log_write( 'Unexpected value for & @@ -1520,42 +1595,43 @@ subroutine ocn_init_validate_Buttermilk_bay(configPool, packagePool, iocontext, return end if -!! call mpas_log_write( ' in ocn_init_validate_buttermilk_bay '//config_Buttermilk_bay_topography_source, MPAS_LOG_CRIT ) + call mpas_log_write( ' in ocn_init_validate_buttermilk_bay '//config_Buttermilk_bay_topography_source ) if (config_Buttermilk_bay_topography_source == 'latlon_file' .or. & config_Buttermilk_bay_topography_source == 'xy_file' ) then !! call mpas_log_write( ' before io_open', MPAS_LOG_CRIT ) - inputFile = MPAS_io_open(config_Buttermilk_bay_topography_file, & - MPAS_IO_READ, MPAS_IO_NETCDF, iocontext_ptr, ierr=iErr) - if (iErr /= 0) then - call mpas_log_write( 'could not open file '// & - & trim(config_Buttermilk_bay_topography_file), MPAS_LOG_CRIT) - return - end if +!! inputFile = MPAS_io_open(trim(config_Buttermilk_bay_topography_file), & +!! MPAS_IO_READ, MPAS_IO_NETCDF, iocontext_ptr, ierr=iErr) +!! if (iErr /= 0) then +!! call mpas_log_write( 'could not open file '// & +!! & trim(config_Buttermilk_bay_topography_file), MPAS_LOG_CRIT) +!! return +!! end if - call mpas_log_write( 'after io_open', MPAS_LOG_CRIT ) +!! call mpas_log_write( 'after io_open', MPAS_LOG_CRIT ) + nLatTopo = 3585 ; + nLonTopo = 3585 ; +!! call MPAS_io_inq_dim(inputFile, config_Buttermilk_bay_topography_nlat_dimname, nLatTopo, iErr) +!! call MPAS_io_inq_dim(inputFile, config_Buttermilk_bay_topography_nlon_dimname, nLonTopo, iErr) - call MPAS_io_inq_dim(inputFile, config_Buttermilk_bay_topography_nlat_dimname, nLatTopo, iErr) - call MPAS_io_inq_dim(inputFile, config_Buttermilk_bay_topography_nlon_dimname, nLonTopo, iErr) +!! call MPAS_io_close(inputFile, iErr) - call MPAS_io_close(inputFile, iErr) +!! call mpas_log_write( ' in ocn_init_validate_buttermilk_bay & +!! after io_close', MPAS_LOG_CRIT ) - call mpas_log_write( ' in ocn_init_validate_buttermilk_bay & - after io_close', MPAS_LOG_CRIT ) - - if ( config_Buttermilk_bay_topography_latlon_degrees .and. & - config_Buttermilk_bay_topography_source == 'latlon_file' ) then - topoLat % array(:) = topoLat % array(:) * pii / 180.0_RKIND - topoLon % array(:) = topoLon % array(:) * pii / 180.0_RKIND - end if +!! if ( config_Buttermilk_bay_topography_latlon_degrees .and. & +!! config_Buttermilk_bay_topography_source == 'latlon_file' ) then +!! topoLat % array(:) = topoLat % array(:) * pii / 180.0_RKIND +!! topoLon % array(:) = topoLon % array(:) * pii / 180.0_RKIND +!! end if end if !---------- - call mpas_log_write( ' Done ocn_init_validate_buttermilk_bay ', MPAS_LOG_CRIT ) + call mpas_log_write( ' Done ocn_init_validate_buttermilk_bay ' ) !-------------------------------------------------------------------- end subroutine ocn_init_validate_ButterMilk_bay!}}} diff --git a/components/mpas-ocean/src/shared/mpas_ocn_subgrid.F b/components/mpas-ocean/src/shared/mpas_ocn_subgrid.F index 6c2a5abcd345..fdda3b99e55f 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_subgrid.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_subgrid.F @@ -134,7 +134,8 @@ subroutine ocn_subgrid_layer_thickness_lookup(zeta, & !layerThick = 0.0_RKIND layerThick = config_drying_min_cell_height + eps !call mpas_log_write('thickness lookup, below table: $r $r $r',realArgs=(/layerThick,zeta,bathymetry/)) - else + else + ! do lev = 1, nSubgridTableLevels-1 zeta0 = (real(lev,RKIND)-1.0_RKIND)*deltaZ + tableMin zeta1 = zeta0 + deltaZ 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..3c7f1d046a59 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tidal_forcing.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tidal_forcing.F @@ -253,6 +253,9 @@ subroutine ocn_tidal_forcing_build_array(domain, meshPool, forcingPool, statePoo end if ! ensure tidalHeight can't push flow below anticipated minimum + ! + ! - DW: will need some mod here + ! if (config_use_wetting_drying .and. tidalInputMask(iCell) == 1.0_RKIND) then ! ensure that tidal height can't force below total minimum thickness ! condition wrong to ensure that there isn't any drying according to criteria @@ -282,6 +285,8 @@ subroutine ocn_tidal_forcing_build_array(domain, meshPool, forcingPool, statePoo end do ! do iCell = 1, nCells +! print*, "Tidal heigh = ", tidalHeight ; + end subroutine ocn_tidal_forcing_build_array !}}} From a11ebce3374eec58d474563a660934eba744d107 Mon Sep 17 00:00:00 2001 From: "D. Wirasaet" Date: Thu, 4 Jan 2024 11:29:55 -0800 Subject: [PATCH 08/16] 2nd attemp comiit --- .../mode_init/mpas_ocn_init_buttermilk_bay.F | 17 ++++++++++++----- .../src/shared/mpas_ocn_init_routines.F | 19 +++++++++++++++++++ .../mpas-ocean/src/shared/mpas_ocn_subgrid.F | 17 ++++++++++++----- 3 files changed, 43 insertions(+), 10 deletions(-) diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F index b46f9106d24c..f4c3cfeb9c43 100644 --- a/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F +++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F @@ -495,13 +495,17 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ bathymetryMin = maxval( subgridBathymetryValues(1:nSubgridCell) ) bathymetryMax = minval( subgridBathymetryValues(1:nSubgridCell) ) - if ( abs(bathymetryMin - bathymetryMax) > 200.0*eps ) then - subgridSshCellTableRange(1,iCell) = -maxval(subgridBathymetryValues(1:nSubgridCell)) + config_drying_min_cell_height + eps + if ( abs(bathymetryMin - bathymetryMax) > 100.0*eps ) then +! subgridSshCellTableRange(1,iCell) = -maxval(subgridBathymetryValues(1:nSubgridCell)) + config_drying_min_cell_height + eps +! subgridSshCellTableRange(2,iCell) = -minval(subgridBathymetryValues(1:nSubgridCell)) + subgridSshCellTableRange(1,iCell) = -maxval(subgridBathymetryValues(1:nSubgridCell)) subgridSshCellTableRange(2,iCell) = -minval(subgridBathymetryValues(1:nSubgridCell)) else ! flat bathy ! - subgridSshCellTableRange(1,iCell) = -bathymetryMin + config_drying_min_cell_height + eps - subgridSshCellTableRange(2,iCell) = -bathymetryMin + 2.0*config_drying_min_cell_height + eps +! subgridSshCellTableRange(1,iCell) = -bathymetryMin + config_drying_min_cell_height + eps +! subgridSshCellTableRange(2,iCell) = -bathymetryMin + 2.0*config_drying_min_cell_height + eps + subgridSshCellTableRange(1,iCell) = -bathymetryMin + config_drying_min_cell_height + subgridSshCellTableRange(2,iCell) = -bathymetryMin + 2.0*config_drying_min_cell_height endif ! Evaluate subgrid bathymetry @@ -519,9 +523,12 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ ! Evaluate wet layerThickness average !--------------------------------------------------------------- !!$ call ocn_init_wet_average(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, layerThickness(1,iCell)) - call ocn_init_wet_average_ssh(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, ssh(iCell))!{{{ enddo + + + + close( unit = 101 ) ; close( unit = 102 ) ; diff --git a/components/mpas-ocean/src/shared/mpas_ocn_init_routines.F b/components/mpas-ocean/src/shared/mpas_ocn_init_routines.F index a037b558cc9c..5c4e3379a313 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_init_routines.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_init_routines.F @@ -122,6 +122,8 @@ subroutine ocn_init_routines_block(block, dt, err)!{{{ type (mpas_pool_iterator_type) :: groupItr + logical, save:: first = .true. + call mpas_pool_get_dimension(block % dimensions, 'nCells', nCells) call mpas_pool_get_dimension(block % dimensions, 'nEdges', nEdges) call mpas_pool_get_dimension(block % dimensions, 'nVertices', nVertices) @@ -153,6 +155,15 @@ subroutine ocn_init_routines_block(block, dt, err)!{{{ end do end if + if ( first ) then + open( unit = 109, file = 'ssh0_before.dat' ) + do iCell = 1, nCells + write(109,'(I7,E17.9)') iCell, ssh( iCell ) ; + end do + close(109) ; + end if + + #ifdef MPAS_OPENACC !$acc enter data copyin(layerThickness, normalVelocity) !$acc update device (normalTransportVelocity, & @@ -268,7 +279,15 @@ subroutine ocn_init_routines_block(block, dt, err)!{{{ err = ior(err, err1) + if ( first ) then + open( unit = 109, file = 'ssh0_afer.dat' ) + do iCell = 1, nCells + write(109,'(I7,E17.9)') iCell, ssh( iCell ) ; + end do + close(109) ; + end if + first = .false. end subroutine ocn_init_routines_block!}}} diff --git a/components/mpas-ocean/src/shared/mpas_ocn_subgrid.F b/components/mpas-ocean/src/shared/mpas_ocn_subgrid.F index fdda3b99e55f..41a97b2cb85e 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_subgrid.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_subgrid.F @@ -271,6 +271,8 @@ subroutine ocn_subgrid_ssh_lookup(layerThick, & real (kind=RKIND) :: zeta0, zeta1 real (kind=RKIND) :: phi0, phi1 + real (kind=RKIND):: zetatmp + tableMin = subgridTableRange(1) tableMax = subgridTableRange(2) nSubgridTableLevels = config_subgrid_table_levels_init @@ -285,17 +287,18 @@ subroutine ocn_subgrid_ssh_lookup(layerThick, & ! call mpas_log_write('ssh, layerThick: $r $r',realArgs=(/zeta0,layerThick0/)) !enddo + zetatmp = zeta ; if (layerThick >= layerThickMax) then zeta = layerThick - bathymetryMean - !call mpas_log_write('ssh lookup, above table: $r $r $r',realArgs=(/layerThick,zeta,bathymetryMean/)) +!S call mpas_log_write('ssh lookup, above table: $r $r $r',realArgs=(/layerThick,zeta,bathymetryMean/)) else if (layerThick <= layerThickMin) then zeta = config_drying_min_cell_height + eps - bathymetryMin - !call mpas_log_write('ssh lookup, below table: $r $r $r',realArgs=(/layerThick,zeta,bathymetryMin/)) +!S call mpas_log_write('ssh lookup, below table: $r $r $r',realArgs=(/layerThick,zeta,bathymetryMin/)) else - do lev = 1, nSubgridTableLevels-1 + levloop: do lev = 1, nSubgridTableLevels-1 zeta0 = (real(lev,RKIND)-1.0_RKIND)*deltaZ + tableMin zeta1 = zeta0 + deltaZ @@ -310,13 +313,17 @@ subroutine ocn_subgrid_ssh_lookup(layerThick, & zeta = phi0*zeta0 + phi1*zeta1 !call mpas_log_write('ssh lookup, in table: $r $r $r $r',realArgs=(/layerThick,zeta,bathymetryMean,zeta+bathymetryMean/)) - return + !return + exit levloop end if - end do + end do levloop end if !call mpas_log_write('$r',realArgs=(/zeta/)) + call mpas_log_write('$r $r $r $r $r', realArgs=(/ & + layerThick,layerThickMin,layerThickMax,zetatmp,zeta /) ) ; + end subroutine ocn_subgrid_ssh_lookup!}}} From bd5cb4ab6cfe5c653bbc3ba314a047948490cb6e Mon Sep 17 00:00:00 2001 From: "D. Wirasaet" Date: Thu, 4 Jan 2024 11:30:23 -0800 Subject: [PATCH 09/16] 2nd attempt --- components/mpas-ocean/src/Registry.xml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index d0ad68d8da35..32bfbe064a43 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -2561,15 +2561,15 @@ description="Pre-computed wet volume per unit area for each edge at given ssh value" packages="subgridWetDryPKG" /> - - - From ce91a35551a80cbe5fa3cf8d4f33bb33b7635974 Mon Sep 17 00:00:00 2001 From: "D. Wirasaet" Date: Thu, 18 Jan 2024 11:42:27 -0800 Subject: [PATCH 10/16] Mile stone for buttermilk bay test case --- .../mpas_ocn_time_integration_rk4.F | 29 ++++- .../mode_init/mpas_ocn_init_buttermilk_bay.F | 84 +++++++++++---- .../src/shared/mpas_ocn_diagnostics.F | 75 ++++++++++--- .../src/shared/mpas_ocn_init_routines.F | 34 +++--- .../mpas-ocean/src/shared/mpas_ocn_subgrid.F | 14 +-- .../src/shared/mpas_ocn_tidal_forcing.F | 20 +++- .../mpas_ocn_vel_forcing_surface_stress.F | 4 + .../src/shared/mpas_ocn_wetting_drying.F | 101 ++++++++++++++++++ 8 files changed, 298 insertions(+), 63 deletions(-) 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 1f0cdb75a8bd..5d553e52f77d 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 @@ -44,6 +44,9 @@ module ocn_time_integration_rk4 use ocn_surface_land_ice_fluxes use ocn_transport_tests + use ocn_subgrid + + implicit none private save @@ -195,6 +198,10 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ real (kind=RKIND), dimension(:,:,:), pointer :: activeTracersCur, activeTracersNew + ! DW + real(kind=RKIND), dimension(:,:), pointer :: subgridSshCellTableRange, subgridWetVolumeCellTable + + ! Get config options call mpas_pool_get_config(domain % configs, 'config_mom_del4', config_mom_del4) call mpas_pool_get_config(domain % configs, 'config_filter_btr_mode', config_filter_btr_mode) @@ -249,6 +256,13 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) + call mpas_pool_get_array(meshPool, 'subgridWetVolumeCellTable', & + subgridWetVolumeCellTable) + + call mpas_pool_get_array(meshPool, 'subgridSshCellTableRange', & + subgridSshCellTableRange) + + ! Lower k-loop limit of 1 rather than minLevel* needed in *New = *Cur ! assignments below are needed to maintain bit-for-bit results @@ -692,9 +706,6 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ !$omp end parallel end if -! print*, config_use_tidal_forcing -! print*, trim(config_tidal_forcing_type) - ! direct application of tidal boundary condition if (config_use_tidal_forcing .and. trim(config_tidal_forcing_type) == 'direct') then do iCell=1, nCells @@ -712,7 +723,17 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ !(1.0_RKIND - tidalInputMask(iCell))*layerThicknessNew(k, iCell) ! generalized tappered assumption code end do - print*, iCell, totalDepth, layerThicknessNew(1, iCell) ; + ! DW + if ( config_use_subgrid_wetting_drying ) then + call ocn_subgrid_layer_thickness_lookup( tidalBCValue(iCell), & + subgridWetVolumeCellTable(:,iCell), & + subgridSshCellTableRange(:,iCell),& + bottomDepth(iCell), & + LayerThicknessNew(1,iCell) ) + + end if + ! + end if end do end if diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F index f4c3cfeb9c43..a33252157800 100644 --- a/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F +++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F @@ -198,6 +198,9 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ real (kind=RKIND) :: x(3), y(3) integer :: slice, nSlice real (kind=RKIND) :: deltaZ + + + ! DW integer:: i, j, jj integer:: nsubgridCellEdge, iEdgeSegment @@ -443,7 +446,6 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ allocate(subgridUValues(maxEdges*nSubgridTriPerSlice)) allocate(subgridVValues(maxEdges*nSubgridTriPerSlice)) - allocate(xSubgridCell(maxEdges*nSubgridTriPerSlice)) allocate(ySubgridCell(maxEdges*nSubgridTriPErSlice)) xSubgridCell = 0.0_RKIND @@ -504,7 +506,7 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ ! flat bathy ! ! subgridSshCellTableRange(1,iCell) = -bathymetryMin + config_drying_min_cell_height + eps ! subgridSshCellTableRange(2,iCell) = -bathymetryMin + 2.0*config_drying_min_cell_height + eps - subgridSshCellTableRange(1,iCell) = -bathymetryMin + config_drying_min_cell_height + subgridSshCellTableRange(1,iCell) = -bathymetryMin subgridSshCellTableRange(2,iCell) = -bathymetryMin + 2.0*config_drying_min_cell_height endif @@ -526,8 +528,18 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ call ocn_init_wet_average_ssh(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, ssh(iCell))!{{{ enddo - - + do iCell = 1, nCellsSolve + call ocn_subgrid_ssh_lookup( config_drying_min_cell_height, & + subgridWetVolumeCellTable(:,iCell), & + subgridSshCellTableRange(:,iCell),& + bottomDepth(iCell),& + subgridCellBathymetryMin(iCell),& + subgridSShCellTableRange(3,iCell) ) + + if ( ssh(iCell) < subgridSshCellTableRange(3,iCell) ) then + ssh(iCell) = subgridSshCellTableRange(3,iCell) ; + end if + end do close( unit = 101 ) ; close( unit = 102 ) ; @@ -589,10 +601,11 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ ! Evaluate bounds of look-up table range !--------------------------------------------------------------- - subgridSshEdgeTableRange(1,iEdge) = -maxval(subgridBathymetryValues(1:nSubgridEdge)) + config_drying_min_cell_height + eps +! subgridSshEdgeTableRange(1,iEdge) = -maxval(subgridBathymetryValues(1:nSubgridEdge)) + config_drying_min_cell_height + eps +! subgridSshEdgeTableRange(2,iEdge) = -minval(subgridBathymetryValues(1:nSubgridEdge)) + subgridSshEdgeTableRange(1,iEdge) = -maxval(subgridBathymetryValues(1:nSubgridEdge)) subgridSshEdgeTableRange(2,iEdge) = -minval(subgridBathymetryValues(1:nSubgridEdge)) - ! Evaluate subgrid bathymetry !--------------------------------------------------------------- @@ -624,7 +637,9 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ ! Evaluate bounds of look-up table range !--------------------------------------------------------------- - subgridSshEdgeTableRange(1,iEdge) = -maxval(CellEdgeBathymetryValues(3,1:nSubgridCellEdge)) + config_drying_min_cell_height + eps +! subgridSshEdgeTableRange(1,iEdge) = -maxval(CellEdgeBathymetryValues(3,1:nSubgridCellEdge)) + config_drying_min_cell_height + eps +! subgridSshEdgeTableRange(2,iEdge) = -minval(CellEdgeBathymetryValues(3,1:nSubgridCellEdge)) + subgridSshEdgeTableRange(1,iEdge) = -maxval(CellEdgeBathymetryValues(3,1:nSubgridCellEdge)) subgridSshEdgeTableRange(2,iEdge) = -minval(CellEdgeBathymetryValues(3,1:nSubgridCellEdge)) ! Evaluate bounds of look-up table range @@ -641,6 +656,19 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ end do + ! + ! find an ssh value corresponding to drying_min_cell-height + ! of each edge + do iEdge = 1,nEdgesSolve + call ocn_subgrid_ssh_lookup( config_drying_min_cell_height, & + subgridWetVolumeEdgeTable(:,iEdge),& + subgridSshEdgeTableRange(:,iEdge),& + subgridEdgeBathymetryMean(iEdge),& + subgridEdgeBathymetryMin(iEdge),& + subgridSshEdgeTableRange(3,iEdge) ) + end do + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Vertex !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -688,10 +716,12 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ ! Evaluate bounds of look-up table range !--------------------------------------------------------------- - - subgridSshVertexTableRange(1,iVertex) = -maxval(subgridBathymetryValues(1:nSubgridVertex)) + config_drying_min_cell_height + eps +! subgridSshVertexTableRange(1,iVertex) = -maxval(subgridBathymetryValues(1:nSubgridVertex)) + config_drying_min_cell_height + eps +! subgridSshVertexTableRange(2,iVertex) = -minval(subgridBathymetryValues(1:nSubgridVertex)) + subgridSshVertexTableRange(1,iVertex) = -maxval(subgridBathymetryValues(1:nSubgridVertex)) subgridSshVertexTableRange(2,iVertex) = -minval(subgridBathymetryValues(1:nSubgridVertex)) + ! Evaluate subgrid bathymetry !--------------------------------------------------------------- subgridVertexBathymetryMean(iVertex) = sum(subgridBathymetryValues(1:nSubgridVertex)*subgridAreas(1:nSubgridVertex))/sum(subgridAreas(1:nSubgridVertex)) @@ -706,6 +736,19 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ enddo vertex + ! + ! find an ssh value corresponding to drying_min_cell-height + ! of each edge + do iVertex = 1,nVerticesSolve + + call ocn_subgrid_ssh_lookup( config_drying_min_cell_height, & + subgridWetVolumeVertexTable(:,iVertex),& + subgridSshVertexTableRange(:,iVertex),& + subgridVertexBathymetryMean(iVertex),& + subgridVertexBathymetryMin(iVertex),& + subgridSshVertexTableRange(3,iVertex) ) + end do + ! Find max bottom depth maxBottomDepth = maxval( bottomDepth ) ; minBottomDepth = minval( bottomDepth ) ; @@ -761,14 +804,14 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ ! subgridCellBathymetryMin(iCell),& ! ssh(iCell)) call ocn_subgrid_layer_thickness_lookup(ssh(iCell), & - subgridWetVolumeCellTable(:,iCell), & - subgridSshCellTableRange(:,iCell),& - bottomDepth(iCell),& - LayerThickness(1,iCell)) + subgridWetVolumeCellTable(:,iCell), & + subgridSshCellTableRange(:,iCell),& + bottomDepth(iCell),& + LayerThickness(1,iCell)) else - call ocn_init_Buttermilk_bay_bathymetry(xCell(iCell),yCell(iCell),bottomDepth(iCell)) - call ocn_init_Buttermilk_bay_ssh(xCell(iCell),yCell(iCell),bottomDepth(iCell),ssh(iCell)) +! call ocn_init_Buttermilk_bay_bathymetry(xCell(iCell),yCell(iCell),bottomDepth(iCell)) +! call ocn_init_Buttermilk_bay_ssh(xCell(iCell),yCell(iCell),bottomDepth(iCell),ssh(iCell)) ssh(iCell) = - bottomDepth(iCell) + & max(ssh(iCell) + bottomDepth(iCell), & maxLevelCell(iCell)*(config_drying_min_cell_height + eps)) @@ -1265,7 +1308,8 @@ subroutine ocn_init_vertical_integration(iCV, subgridSshTableRange, nSubgridCV, deltaZ = (subgridSshTableRange(2,iCV)-subgridSshTableRange(1,iCV))/real(nSubgridTableLevels-1,RKIND); - subgridWetVolumeTable(1,iCV) = config_drying_min_cell_height + eps + ! subgridWetVolumeTable(1,iCV) = config_drying_min_cell_height + eps + subgridWetVolumeTable(1,iCV) = 0.0_RKIND ssh = subgridSshTableRange(1,iCV) + deltaZ subgridWetVolumeTable(2,iCV) = subgridWetVolumeTable(1,iCV)*sum(subgridAreas(1:nSubgridCV)) subgridWetFractionTable(1,iCV) = 0.0_RKIND @@ -1390,7 +1434,8 @@ subroutine ocn_init_wet_average_ssh( nSubgridCV, subgridBathymetryValues, subgri do tri = 1, nsubgridCV layerThicknessValue = subgridBathymetryValues(tri) + subgridSshValues(tri) - if (layerThicknessValue > config_drying_min_cell_height + eps ) then +! if (layerThicknessValue > config_drying_min_cell_height + eps ) then + if ( layerThicknessValue > 0.0_RKIND ) then sshWetAverage = sshWetAverage + subgridSshValues(tri)*subgridAreas(tri) wetArea = wetArea + subgridAreas(tri) endif @@ -1399,8 +1444,9 @@ subroutine ocn_init_wet_average_ssh( nSubgridCV, subgridBathymetryValues, subgri if ( WetArea > 0.0_RKIND ) then sshWetAverage = sshWetAverage/WetArea ; else - sshWetAverage = -maxval(subgridBathymetryValues(1:nsubgridCV)) & - + config_drying_min_cell_height ; +! sshWetAverage = -maxval(subgridBathymetryValues(1:nsubgridCV)) & +! + config_drying_min_cell_height ; + sshWetAverage = -maxval(subgridBathymetryValues(1:nsubgridCV)) ; endif return ; diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F index 099dc698dce7..14936ac1c68a 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F @@ -719,6 +719,9 @@ subroutine ocn_diagnostic_solve_layerThicknessEdge(normalVelocity, & layerThicknessCell1, & layerThicknessCell2 + real(kind=RKIND):: eps = 1.0e-10_RKIND + + ! End preamble !----------------------------------------------------------------- ! Begin code @@ -738,12 +741,15 @@ subroutine ocn_diagnostic_solve_layerThicknessEdge(normalVelocity, & do iEdge = 1, nEdgesAll cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) - sshMean = 0.5_RKIND * (ssh(cell1) + ssh(cell2)) - call ocn_subgrid_layer_thickness_lookup(sshMean, & - subgridWetVolumeEdgeTable(:,iEdge), & - subgridSshEdgeTableRange(:,iEdge), & - subgridEdgeBathymetryMean(iEdge), & - layerThickEdgeMean(1,iEdge)) + !sshMean = 0.5_RKIND * (ssh(cell1) + ssh(cell2)) + !call ocn_subgrid_layer_thickness_lookup(sshMean, & + ! subgridWetVolumeEdgeTable(:,iEdge), & + ! subgridSshEdgeTableRange(:,iEdge), & + ! subgridEdgeBathymetryMean(iEdge), & + ! layerThickEdgeMean(1,iEdge)) + layerThickEdgeMean(1,iEdge) = 0.5_RKIND * & + ( layerThickness(1,cell1) + & + layerThickness(1,cell2) ) end do #ifndef MPAS_OPENACC !$omp end do @@ -796,12 +802,33 @@ subroutine ocn_diagnostic_solve_layerThicknessEdge(normalVelocity, & !$omp parallel !$omp do schedule(runtime) private(k) #endif - do iEdge = 1, nEdgesAll - do k = 1,nVertLevels - layerThickEdgeFlux(k,iEdge) = & - layerThickEdgeMean(k,iEdge) - end do - end do + if ( config_use_subgrid_wetting_drying ) then + ! + do iEdge = 1, nEdgesAll + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + sshMean = 0.5_RKIND * (ssh(cell1) + ssh(cell2)) + call ocn_subgrid_layer_thickness_lookup( sshMean, & + subgridWetVolumeEdgeTable(:,iEdge), & + subgridSshEdgeTableRange(:,iEdge), & + subgridEdgeBathymetryMean(iEdge), & + layerThickEdgeFlux(1,iEdge) ) + + if ( layerThickEdgeFlux(1,iEdge) < eps ) then + layerThickEdgeFlux(1,iEdge) = layerThickEdgeMean(1,iEdge) ; + end if + end do + ! + else + + do iEdge = 1, nEdgesAll + do k = 1,nVertLevels + layerThickEdgeFlux(k,iEdge) = & + layerThickEdgeMean(k,iEdge) + end do + end do + + end if #ifndef MPAS_OPENACC !$omp end do !$omp end parallel @@ -858,6 +885,16 @@ subroutine ocn_diagnostic_solve_layerThicknessEdge(normalVelocity, & subgridEdgeBathymetryMean(iEdge), & layerThickEdgeFlux(1,iedge)) end if + + if ( layerThickEdgeFlux(1,iedge) < eps ) then + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + + layerThickEdgeFlux(1,iEdge) = 0.5_RKIND * & + ( layerThickness(1,cell1) + & + layerThickness(1,cell2) ) + end if + ! print*, iEdge, layerThickEdgeFlux(1,iedge), layerThickEdgeMean(1,iedge) end do #ifndef MPAS_OPENACC !$omp end do @@ -2501,6 +2538,8 @@ subroutine ocn_diagnostic_solve_z_coordinates(layerThickness, zMid, zTop, ssh)!{ real (kind=RKIND), dimension(:), intent(out) :: & ssh + logical, save:: first = .true. + !----------------------------------------------------------------- ! ! local variables @@ -2549,14 +2588,24 @@ subroutine ocn_diagnostic_solve_z_coordinates(layerThickness, zMid, zTop, ssh)!{ bottomDepth(iCell), & subgridCellBathymetryMin(iCell), & ssh(iCell)) + +! if ( first ) then +! if ( abs(ssh(iCell)) < 1.0e-10_RKIND ) then +! ssh(iCell) = 0.0_RKIND +! end if +! endif + zTop(1,iCell) = ssh(iCell) - zMid(1,iCell) = -bottomDepth(iCell) + 0.5_RKIND*layerThickness(1,iCell) + zMid(1,iCell) = -bottomDepth(iCell) + 0.5_RKIND*layerThickness(1,iCell) else ! copy zTop(1,iCell) into sea-surface height array ssh(iCell) = zTop(minLevelCell(iCell),iCell) end if end do + + first = .false. ; + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel diff --git a/components/mpas-ocean/src/shared/mpas_ocn_init_routines.F b/components/mpas-ocean/src/shared/mpas_ocn_init_routines.F index 5c4e3379a313..d177e459e883 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_init_routines.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_init_routines.F @@ -122,7 +122,7 @@ subroutine ocn_init_routines_block(block, dt, err)!{{{ type (mpas_pool_iterator_type) :: groupItr - logical, save:: first = .true. +! logical, save:: first = .true. call mpas_pool_get_dimension(block % dimensions, 'nCells', nCells) call mpas_pool_get_dimension(block % dimensions, 'nEdges', nEdges) @@ -155,13 +155,13 @@ subroutine ocn_init_routines_block(block, dt, err)!{{{ end do end if - if ( first ) then - open( unit = 109, file = 'ssh0_before.dat' ) - do iCell = 1, nCells - write(109,'(I7,E17.9)') iCell, ssh( iCell ) ; - end do - close(109) ; - end if +! if ( first ) then +! open( unit = 109, file = 'ssh0_before.dat' ) +! do iCell = 1, nCells +! write(109,'(I7,E17.9)') iCell, ssh( iCell ) ; +! end do +! close(109) ; +! end if #ifdef MPAS_OPENACC @@ -279,15 +279,15 @@ subroutine ocn_init_routines_block(block, dt, err)!{{{ err = ior(err, err1) - if ( first ) then - open( unit = 109, file = 'ssh0_afer.dat' ) - do iCell = 1, nCells - write(109,'(I7,E17.9)') iCell, ssh( iCell ) ; - end do - close(109) ; - end if - - first = .false. +! if ( first ) then +! open( unit = 109, file = 'ssh0_afer.dat' ) +! do iCell = 1, nCells +! write(109,'(I7,E17.9)') iCell, ssh( iCell ) ; +! end do +! close(109) ; +! end if +! +! first = .false. end subroutine ocn_init_routines_block!}}} diff --git a/components/mpas-ocean/src/shared/mpas_ocn_subgrid.F b/components/mpas-ocean/src/shared/mpas_ocn_subgrid.F index 41a97b2cb85e..cc231a168c4e 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_subgrid.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_subgrid.F @@ -131,8 +131,8 @@ subroutine ocn_subgrid_layer_thickness_lookup(zeta, & layerThick = zeta + bathymetry !call mpas_log_write('thickness lookup, above table: $r $r $r',realArgs=(/layerThick,zeta,bathymetry/)) else if (zeta <= tableMin) then - !layerThick = 0.0_RKIND - layerThick = config_drying_min_cell_height + eps + layerThick = 0.0_RKIND + ! layerThick = config_drying_min_cell_height + eps !call mpas_log_write('thickness lookup, below table: $r $r $r',realArgs=(/layerThick,zeta,bathymetry/)) else ! @@ -188,7 +188,7 @@ subroutine ocn_subgrid_wet_fraction_lookup(zeta, & !----------------------------------------------------------------- integer :: lev - integer :: nSubgridTableLevels + integer :: nSubgrsdTableLevels real (kind=RKIND) :: tableMin real (kind=RKIND) :: tableMax @@ -293,9 +293,9 @@ subroutine ocn_subgrid_ssh_lookup(layerThick, & !S call mpas_log_write('ssh lookup, above table: $r $r $r',realArgs=(/layerThick,zeta,bathymetryMean/)) else if (layerThick <= layerThickMin) then - zeta = config_drying_min_cell_height + eps - bathymetryMin + zeta = - bathymetryMin ! prevent_drying likely fails to ensure positive water columbn +!$ zeta = config_drying_min_cell_height + eps - bathymetryMin !S call mpas_log_write('ssh lookup, below table: $r $r $r',realArgs=(/layerThick,zeta,bathymetryMin/)) - else levloop: do lev = 1, nSubgridTableLevels-1 @@ -321,8 +321,8 @@ subroutine ocn_subgrid_ssh_lookup(layerThick, & end if !call mpas_log_write('$r',realArgs=(/zeta/)) - call mpas_log_write('$r $r $r $r $r', realArgs=(/ & - layerThick,layerThickMin,layerThickMax,zetatmp,zeta /) ) ; + !call mpas_log_write('$r $r $r $r $r', realArgs=(/ & + ! layerThick,layerThickMin,layerThickMax,zetatmp,zeta /) ) ; end subroutine ocn_subgrid_ssh_lookup!}}} 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 3c7f1d046a59..5208405340f2 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tidal_forcing.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tidal_forcing.F @@ -218,6 +218,8 @@ subroutine ocn_tidal_forcing_build_array(domain, meshPool, forcingPool, statePoo real (kind=RKIND) :: totalDepth, tidalHeight character (len=StrKIND), pointer :: xtime + + real (kind=RKIND), dimension(:), pointer:: subgridSShCekkTableRange if ( .not. tidalfluxOn ) return @@ -233,6 +235,12 @@ subroutine ocn_tidal_forcing_build_array(domain, meshPool, forcingPool, statePoo call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) + + ! DW + call mpas_pool_get_array(meshPool, 'subgridSshCellTableRange', & + subgridSshCellTableRange) + + ! tidal fields are needed only over 0 and 1 halos nCells = nCellsArray( 2 ) @@ -254,12 +262,19 @@ subroutine ocn_tidal_forcing_build_array(domain, meshPool, forcingPool, statePoo ! ensure tidalHeight can't push flow below anticipated minimum ! - ! - DW: will need some mod here + ! - DW: will need some modfication here ! if (config_use_wetting_drying .and. tidalInputMask(iCell) == 1.0_RKIND) then ! ensure that tidal height can't force below total minimum thickness ! condition wrong to ensure that there isn't any drying according to criteria - tidalHeight = max(-bottomDepth(iCell) + (float(maxLevelCell(iCell))+1.0_RKIND)*config_drying_min_cell_height, tidalHeight) + ! tidalHeight = max(-bottomDepth(iCell) + (float(maxLevelCell(iCell))+1.0_RKIND)*config_drying_min_cell_height, tidalHeight) + + ! DW + if ( config_use_subgrid_wetting_drying ) then + tidalHeight = max(tidalHeight, subgridSShCellTableRange(3,iCell) ) ; + else + tidalHeight = max(-bottomDepth(iCell) + (float(maxLevelCell(iCell))+1.0_RKIND)*config_drying_min_cell_height, tidalHeight) + end if end if ! compute total depth for relative thickness contribution @@ -285,7 +300,6 @@ subroutine ocn_tidal_forcing_build_array(domain, meshPool, forcingPool, statePoo end do ! do iCell = 1, nCells -! print*, "Tidal heigh = ", tidalHeight ; end subroutine ocn_tidal_forcing_build_array !}}} diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_forcing_surface_stress.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_forcing_surface_stress.F index 216252b32c6d..762a42669be3 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vel_forcing_surface_stress.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_forcing_surface_stress.F @@ -166,10 +166,14 @@ subroutine ocn_vel_forcing_surface_stress_tend(sfcFlxAttCoeff, & absorb = transCoeffTop - transCoeffBot remainingStress = remainingStress - absorb + ! + ! print*, iEdge, layerThickEdgeMean(k,iEdge) ; + ! tend(k,iEdge) = tend(k,iEdge) + & edgeMask(k,iEdge)*sfcStress(iEdge) * & absorb/rho_sw/layerThickEdgeMean(k,iEdge) + zTop = zBot transCoeffTop = transCoeffBot enddo 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 aeb1dcdb4b7d..b6dc8eb4c192 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F @@ -241,6 +241,9 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying integer :: iEdge, k + logical, save:: first = .true. + integer:: iCell, ii + err = 0 call mpas_pool_get_subpool(block % structs, 'tend', tendPool) @@ -264,9 +267,44 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying ! ensure cells stay wet by selectively damping cells with a damping tendency to make ! sure tendency doesn't dry cells + +! if ( first ) then +! open( unit = 100, file = 'edoncell.dat' ) ; +! do iCell = 1, nCellsAll +! write( 100, '(I10)', advance = 'no' ) iCell +! do ii = 1, nEdgesOnCell(iCell) +! iEdge = edgesOnCell(ii, iCell) ; +! +! write( 100, '(I10)', advance = 'no' ) iEdge +! end do +! write( 100, * ) '' ; +! end do +! close(100) +! +! ! +! open( unit = 100, file = 'cellonedge.dat' ) ; +! do iEdge = 1, nEdgesAll +! write( 100, '(3I10)', advance = 'no' ) iEdge, cellsOnEdge(1:2,iEdge) ; +! end do +! write(100,*) '' +! close(100) +! ! +! end if + call ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, & normalTransportVelocity, rkSubstepWeight, wettingVelocityFactor, err) +! if ( first ) then +! open( unit = 101, file = 'wdfactor.dat' ) +! +! do iEdge = 1, nEdgesAll +! write( 101, '(I10,2F16.8,2I10)' ) iEdge, wettingVelocityFactor(1,iEdge), normalVelocity(1,iEdge), minLevelEdgeTop(iEdge), maxLevelEdgeBot(iEdge) +! end do +! +! close( 101 ) ; +! end if +! first = .false. ; + ! prevent drying from happening with selective wettingVelocityFactor if (config_zero_drying_velocity) then !$omp parallel @@ -302,6 +340,7 @@ end subroutine ocn_prevent_drying_rk4 !}}} !> to prevent cells from drying. ! !----------------------------------------------------------------------- + subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, & normalVelocity, dt, wettingVelocityFactor, err)!{{{ @@ -357,17 +396,27 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness real (kind=RKIND) :: hCrit, hRampMin, hEdgeTotal character (len=100) :: log_string + +! logical, save:: first = .true. +! integer:: dryflag, ii err = 0 hCrit = config_drying_min_cell_height + +! if ( first ) then +! open( unit = 102, file = 'wddiag.dat' ) ; +! open( unit = 103, file = 'wdflaged.dat') ; +! end if + if (.not. config_zero_drying_velocity) return ! need predicted transport velocity to limit drying flux !$omp parallel !$omp do schedule(runtime) private(i, iEdge, k, divOutFlux, layerThickness) do iCell = 1, nCellsAll +! dryflag = 0 ; do k = minLevelCell(iCell), maxLevelCell(iCell) divOutFlux = 0.0_RKIND layerThickness = min(layerThicknessProvis(k, iCell), layerThicknessCur(k, iCell)) @@ -384,6 +433,7 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness end if end do layerThickness = layerThickness + dt * divOutFlux + ! if layer thickness is too small, limit divergence flux outwards with ! opposite velocity @@ -397,6 +447,17 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness end if end if end do +! dryflag = 1 ; + +! if ( first ) then +! write( 103, '(2I10)', advance = 'no') iCell, dryflag ; +! do ii = 1, nEdgesOnCell(iCell) +! iEdge = edgesOnCell(ii, iCell) +! write( 103, '(I8,F16.8)', advance = 'no' ) iEdge, wettingVelocityFactor(1, iEdge) ; +! end do +! write( 103, *) '' ; +! end if + elseif (config_zero_drying_velocity_ramp .and. & (layerThickness > & hCrit + config_drying_safety_height) .and. & @@ -414,14 +475,54 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness end if end if end do +! dryflag = -1 +! if ( first ) then +! write( 103, '(2I10)', advance = 'no') iCell, dryflag ; +! do ii = 1, nEdgesOnCell(iCell) +! iEdge = edgesOnCell(ii, iCell) +! write( 103, '(I8,F16.8)', advance = 'no' ) iEdge, wettingVelocityFactor(1, iEdge) ; +! end do +! write( 103, *) '' ; +! end if end if +! else +! +! if ( first ) then +! write( 103, '(2I10)', advance = 'no') iCell, dryflag ; +! do ii = 1, nEdgesOnCell(iCell) +! iEdge = edgesOnCell(ii, iCell) +! write( 103, '(I8,F16.8)', advance = 'no' ) iEdge, wettingVelocityFactor(1, iEdge) ; +! end do +! write( 103, *) '' ; +! end if +! +! end if + +! if ( first ) then +! write( 102, '(I10,2F16.8,I5)') iCell, layerThickness, dt*divOutFlux, dryflag +! end if end do end do !$omp end do !$omp end parallel +! if ( first ) close(102) ; +! if ( first ) close(103) ; +! +! if ( first ) then +! open( unit = 101, file = 'wdfactor1.dat' ) +! +! do iEdge = 1, nEdgesAll +! write( 101, '(F16.8)' ) wettingVelocityFactor(1,iEdge) +! end do +! +! close( 101 ) ; +! end if +! +! first = .false. ; + end subroutine ocn_wetting_drying_wettingVelocity !}}} From cea952e44f1a8de6c5f153b0243597c1f1c440df Mon Sep 17 00:00:00 2001 From: "D. Wirasaet" Date: Thu, 25 Jan 2024 10:53:51 -0800 Subject: [PATCH 11/16] check point alternative wet/dry implemenation --- components/mpas-ocean/src/shared/Makefile | 4 +- .../src/shared/mpas_ocn_wetting_drying.F | 223 +++++++++++++----- 2 files changed, 167 insertions(+), 60 deletions(-) diff --git a/components/mpas-ocean/src/shared/Makefile b/components/mpas-ocean/src/shared/Makefile index 9f867d138d53..3977c1af36f2 100644 --- a/components/mpas-ocean/src/shared/Makefile +++ b/components/mpas-ocean/src/shared/Makefile @@ -69,14 +69,14 @@ OBJS = mpas_ocn_init_routines.o \ mpas_ocn_time_average_coupled.o \ mpas_ocn_framework_forcing.o \ mpas_ocn_time_varying_forcing.o \ + mpas_ocn_subgrid.o \ mpas_ocn_wetting_drying.o \ mpas_ocn_vel_tidal_potential.o \ mpas_ocn_vel_forcing_topographic_wave_drag.o \ mpas_ocn_transport_tests.o \ mpas_ocn_vel_self_attraction_loading.o \ mpas_ocn_vertical_advection.o \ - mpas_ocn_stokes_drift.o \ - mpas_ocn_subgrid.o + mpas_ocn_stokes_drift.o all: $(OBJS) 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 b6dc8eb4c192..2f6939d2e5ee 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F @@ -34,6 +34,9 @@ module ocn_wetting_drying use ocn_gm use ocn_mesh + ! DW + use ocn_subgrid + implicit none private save @@ -291,9 +294,15 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying ! ! ! end if - call ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, & - normalTransportVelocity, rkSubstepWeight, wettingVelocityFactor, err) + if ( config_use_subgrid_wetting_drying .and. & + config_pressure_gradient_type == 'ssh_gradient' ) then + call ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, & + normalTransportVelocity, rkSubstepWeight, wettingVelocityFactor, err, block) + else + call ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, & + normalTransportVelocity, rkSubstepWeight, wettingVelocityFactor, err) + end if ! if ( first ) then ! open( unit = 101, file = 'wdfactor.dat' ) ! @@ -343,7 +352,7 @@ end subroutine ocn_prevent_drying_rk4 !}}} subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, & - normalVelocity, dt, wettingVelocityFactor, err)!{{{ + normalVelocity, dt, wettingVelocityFactor, err, block )!{{{ !----------------------------------------------------------------- ! @@ -366,6 +375,10 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness real (kind=RKIND), intent(in) :: & dt !< Input: time step + ! optional + type (block_type), optional, intent(in) :: block + + !----------------------------------------------------------------- ! ! input/output variables @@ -396,7 +409,19 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness real (kind=RKIND) :: hCrit, hRampMin, hEdgeTotal character (len=100) :: log_string + + ! + type (mpas_pool_type), pointer :: meshPool + ! array pointers for subgrid look-up tables + real (kind=RKIND), dimension(:), pointer :: bottomDepth + real(kind=RKIND), dimension(:,:), pointer :: subgridSshCellTableRange, subgridWetVolumeCellTable + real(kind=RKIND), dimension(:), pointer :: subgridCellBathymetryMin + + ! + integer:: cellDummy(2), cellCur, CellNei, sshCur, sshNei + real (kind=RKIND), allocatable:: layerThicknessCell(:), sshCell(:) + ! logical, save:: first = .true. ! integer:: dryflag, ii @@ -405,11 +430,27 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness hCrit = config_drying_min_cell_height -! if ( first ) then -! open( unit = 102, file = 'wddiag.dat' ) ; -! open( unit = 103, file = 'wdflaged.dat') ; -! end if + ! for 1-layer subgrid wetting/dring ! + if ( config_use_subgrid_wetting_drying .and. & + config_pressure_gradient_type == 'ssh_gradient' ) then + allocate( layerThicknessCell(nCellsAll) ) + allocate( sshCell(nCellsAll) ) + layerThickness = 0.0_RKIND + sshCell = 0.0_RKIND + + call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) + + call mpas_pool_get_array( meshPool, 'bottomDepth', bottomDepth) + call mpas_pool_get_array( meshPool, 'subgridWetVolumeCellTable', & + subgridWetVolumeCellTable ) + call mpas_pool_get_array( meshPool, 'subgridSshCellTableRange', & + subgridSshCellTableRange ) + call mpas_pool_get_array( meshPool, 'subgridCellBathymetryMin', & + subgridCellBathymetryMin ) + + end if + if (.not. config_zero_drying_velocity) return ! need predicted transport velocity to limit drying flux @@ -433,7 +474,21 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness end if end do layerThickness = layerThickness + dt * divOutFlux - + + ! + if ( config_use_subgrid_wetting_drying .and. & + config_pressure_gradient_type == 'ssh_gradient' ) then + layerThicknessCell(iCell) = layerThickness + + call ocn_subgrid_ssh_lookup( layerThicknessCell(iCell), & + subgridWetVolumeCellTable(:,iCell), & + subgridSshCellTableRange(:,iCell), & + bottomDepth(iCell), & + subgridCellBathymetryMin(iCell), & + sshCell(iCell) ) + + end if + ! ! if layer thickness is too small, limit divergence flux outwards with ! opposite velocity @@ -447,16 +502,6 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness end if end if end do -! dryflag = 1 ; - -! if ( first ) then -! write( 103, '(2I10)', advance = 'no') iCell, dryflag ; -! do ii = 1, nEdgesOnCell(iCell) -! iEdge = edgesOnCell(ii, iCell) -! write( 103, '(I8,F16.8)', advance = 'no' ) iEdge, wettingVelocityFactor(1, iEdge) ; -! end do -! write( 103, *) '' ; -! end if elseif (config_zero_drying_velocity_ramp .and. & (layerThickness > & @@ -475,54 +520,116 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness end if end if end do -! dryflag = -1 -! if ( first ) then -! write( 103, '(2I10)', advance = 'no') iCell, dryflag ; -! do ii = 1, nEdgesOnCell(iCell) -! iEdge = edgesOnCell(ii, iCell) -! write( 103, '(I8,F16.8)', advance = 'no' ) iEdge, wettingVelocityFactor(1, iEdge) ; -! end do -! write( 103, *) '' ; -! end if - end if -! else -! -! if ( first ) then -! write( 103, '(2I10)', advance = 'no') iCell, dryflag ; -! do ii = 1, nEdgesOnCell(iCell) -! iEdge = edgesOnCell(ii, iCell) -! write( 103, '(I8,F16.8)', advance = 'no' ) iEdge, wettingVelocityFactor(1, iEdge) ; -! end do -! write( 103, *) '' ; -! end if -! -! end if - -! if ( first ) then -! write( 102, '(I10,2F16.8,I5)') iCell, layerThickness, dt*divOutFlux, dryflag -! end if end do end do !$omp end do !$omp end parallel -! if ( first ) close(102) ; -! if ( first ) close(103) ; -! -! if ( first ) then -! open( unit = 101, file = 'wdfactor1.dat' ) -! -! do iEdge = 1, nEdgesAll -! write( 101, '(F16.8)' ) wettingVelocityFactor(1,iEdge) -! end do -! -! close( 101 ) ; -! end if -! -! first = .false. ; + ! DW implementation ! + if ( config_use_modify_wetting_drying_implementation ) then + ! + if ( config_use_subgrid_wetting_drying .and. & + config_pressure_gradient_type == 'ssh_gradient' ) then + ! + wettingVelocityFactor = 0.0_RKIND + ! + cellLoop: do iCell = 1, nCellsAll + ! if layer thickness is too small, limit divergence flux outwards with + ! opposite velocity + ifDryCell: if ( layerThicknessCell(iCell) <= & + hCrit + config_drying_safety_height ) then + + edgeLoop: do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + ! + if ( normalVelocity(1,iEdge)*edgeSignOnCell(i, iCell) < 0.0_RKIND ) then + ! + wettingVelocityFactor(1,iEdge) = 1.0_RKIND + + else if ( normalVelocity(1,iEdge)*edgeSignOnCell(i, iCell) == 0.0_RKIND ) then + ! + wettingVelocityFactor(1,iEdge) = 1.0_RKIND + + cellDummy(1:2) = cellsOnEdge(1:2,iEdge) + + cellCur = iCell + cellNei = merge( cellDummy(2), cellDummy(1), iCell == cellDummy(1) ) + if ( cellNei > nCellsAll ) cellNei = cellCur ! check with Steven on how to handle this index + + sshCur = sshCell(cellCur) ; + sshNei = sshCell(cellNei) ; + + ! if the neigbor cell is anticiapted wet ! + if ( layerThicknessCell(cellNei) > & + hCrit + config_drying_safety_height ) then + if ( sshCur < sshNei ) then + wettingVelocityFactor(1,iEdge) = 0.0_RKIND + end if + end if + ! + end if + end do edgeLoop + + end if ifDryCell + end do cellLoop + + ! Ramping + hRampMin = config_zero_drying_velocity_ramp_hmin + cellLoopRamp: do iCell = 1, nCellsAll + ! + ifRampCell: if ( config_zero_drying_velocity_ramp .and. & + (layerThicknessCell(iCell) > & + hCrit + config_drying_safety_height) .and. & + (layerThicknessCell(iCell) <= config_zero_drying_velocity_ramp_hmax) ) then + + ! + ! Following O'Dea et al. (2021), if total upwinded wct is less than + ! 2*critical thickness, apply damping at each edge + edgeLoopRamp: do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + if ( normalVelocity(1, iEdge) * edgeSignOnCell(i, iCell) <= 0.0_RKIND ) then + wettingVelocityFactor(1, iEdge) = 1.0_RKIND - & + tanh(50.0_RKIND * (layerThicknessCell(iCell) - hRampMin)/hRampMin) + end if + ! + if ( normalVelocity(1, iEdge) * edgeSignOnCell(i,iCell) == 0.0_RKIND ) then + cellDummy(1:2) = cellsOnEdge(1:2,iEdge) + + cellCur = iCell ; + cellNei = merge( cellDummy(2), cellDummy(1), iCell == cellDummy(1) ) + if ( cellNei > nCellsAll ) cellNei = cellCur ! check + ! with Steven on how to handle this + + sshCur = sshCell(cellCur) ; + sshNei = sshCell(cellNei) ; + ! if the neigbor cell is anticiapted dry ! + if ( layerThicknessCell(cellNei) <= & + hCrit + config_drying_safety_height ) then + if ( sshCur < sshNei ) then + wettingVelocityFactor(1,iEdge) = 1.0_RKIND + end if + end if + end if + ! + end do edgeLoopRamp + ! + end if ifRampCell + ! + end do cellLoopRamp + ! + end if + ! + end if + + if ( config_use_subgrid_wetting_drying .and. & + config_pressure_gradient_type == 'ssh_gradient' ) then + deallocate( layerThicknessCell ) ; + deallocate( sshCell ) ; + end if + end subroutine ocn_wetting_drying_wettingVelocity !}}} From 9c4aa8378256de81182966a3cd5bff27921a79e4 Mon Sep 17 00:00:00 2001 From: "D. Wirasaet" Date: Mon, 29 Jan 2024 11:05:33 -0800 Subject: [PATCH 12/16] add non subgrid set up to the buttermilk bay test case --- components/mpas-ocean/src/Registry.xml | 4 + .../src/mode_init/Registry_Buttermilk_Bay.xml | 5 +- .../mode_init/mpas_ocn_init_buttermilk_bay.F | 43 ++---- .../src/shared/mpas_ocn_wetting_drying.F | 140 ++++++++++-------- 4 files changed, 104 insertions(+), 88 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 32bfbe064a43..b7eabd6baf47 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -1164,6 +1164,10 @@ + diff --git a/components/mpas-ocean/src/mode_init/Registry_Buttermilk_Bay.xml b/components/mpas-ocean/src/mode_init/Registry_Buttermilk_Bay.xml index 5724a2e8ca13..2b7e46a2ce64 100644 --- a/components/mpas-ocean/src/mode_init/Registry_Buttermilk_Bay.xml +++ b/components/mpas-ocean/src/mode_init/Registry_Buttermilk_Bay.xml @@ -22,11 +22,14 @@ + 100.0*eps ) then -! subgridSshCellTableRange(1,iCell) = -maxval(subgridBathymetryValues(1:nSubgridCell)) + config_drying_min_cell_height + eps -! subgridSshCellTableRange(2,iCell) = -minval(subgridBathymetryValues(1:nSubgridCell)) subgridSshCellTableRange(1,iCell) = -maxval(subgridBathymetryValues(1:nSubgridCell)) subgridSshCellTableRange(2,iCell) = -minval(subgridBathymetryValues(1:nSubgridCell)) else ! flat bathy ! -! subgridSshCellTableRange(1,iCell) = -bathymetryMin + config_drying_min_cell_height + eps -! subgridSshCellTableRange(2,iCell) = -bathymetryMin + 2.0*config_drying_min_cell_height + eps subgridSshCellTableRange(1,iCell) = -bathymetryMin subgridSshCellTableRange(2,iCell) = -bathymetryMin + 2.0*config_drying_min_cell_height endif @@ -537,7 +531,7 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ subgridSShCellTableRange(3,iCell) ) if ( ssh(iCell) < subgridSshCellTableRange(3,iCell) ) then - ssh(iCell) = subgridSshCellTableRange(3,iCell) ; + ssh(iCell) = subgridSshCellTableRange(3,iCell) ; end if end do @@ -795,14 +789,9 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ ! add a thin layer of nlayer*config_drying_min_cellhight ! - if (config_use_subgrid_wetting_drying) then - -! call ocn_subgrid_ssh_lookup(layerThickness(1,iCell),& -! subgridWetVolumeCellTable(:,iCell),& -! subgridSshCellTableRange(:,iCell),& -! bottomDepth(iCell),& -! subgridCellBathymetryMin(iCell),& -! ssh(iCell)) +! if (config_use_subgrid_wetting_drying) then + if ( config_Buttermilk_bay_subgrid_setup ) then + ! Intical contion for a subgrid run call ocn_subgrid_layer_thickness_lookup(ssh(iCell), & subgridWetVolumeCellTable(:,iCell), & subgridSshCellTableRange(:,iCell),& @@ -810,25 +799,25 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ LayerThickness(1,iCell)) else -! call ocn_init_Buttermilk_bay_bathymetry(xCell(iCell),yCell(iCell),bottomDepth(iCell)) -! call ocn_init_Buttermilk_bay_ssh(xCell(iCell),yCell(iCell),bottomDepth(iCell),ssh(iCell)) - ssh(iCell) = - bottomDepth(iCell) + & - max(ssh(iCell) + bottomDepth(iCell), & - maxLevelCell(iCell)*(config_drying_min_cell_height + eps)) - + ! Initial condition for a standard run + ssh(iCell) = -bottomDepth(iCell) ; do k = 1, maxLevelCell(iCell) - layerThickness(k,iCell) = max(config_drying_min_cell_height + eps, & - (ssh(iCell) + bottomDepth(iCell))/real(maxLevelCell(iCell),RKIND)) - + ! + layerThickness(k,iCell) = max( config_drying_min_cell_height, & + bottomDepth(iCell)/real(maxLevelCell(iCell),RKIND) ) + + if (layerThickness(k,iCell) < config_drying_min_cell_height) then call mpas_log_write('layerThickness($i,$i)=$r', MPAS_LOG_CRIT, & intArgs=(/k,iCell/), & realArgs=(/layerThickness(k,iCell)/)) end if + + ssh(iCell) = ssh(iCell) + layerThickness(k,iCell) ; end do - endif - - + ! + endif + do k = 1, maxLevelCell(iCell) restingThickness(k,iCell) = bottomDepth(iCell)/maxLevelCell(iCell) end do 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 2f6939d2e5ee..ffb2df230912 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F @@ -294,8 +294,8 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying ! ! ! end if - if ( config_use_subgrid_wetting_drying .and. & - config_pressure_gradient_type == 'ssh_gradient' ) then + if ( config_use_modify_wetting_drying_implementation .and. & + config_pressure_gradient_type == 'ssh_gradient' ) then call ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, & normalTransportVelocity, rkSubstepWeight, wettingVelocityFactor, err, block) @@ -430,26 +430,6 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness hCrit = config_drying_min_cell_height - ! for 1-layer subgrid wetting/dring ! - if ( config_use_subgrid_wetting_drying .and. & - config_pressure_gradient_type == 'ssh_gradient' ) then - allocate( layerThicknessCell(nCellsAll) ) - allocate( sshCell(nCellsAll) ) - layerThickness = 0.0_RKIND - sshCell = 0.0_RKIND - - call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) - - call mpas_pool_get_array( meshPool, 'bottomDepth', bottomDepth) - call mpas_pool_get_array( meshPool, 'subgridWetVolumeCellTable', & - subgridWetVolumeCellTable ) - call mpas_pool_get_array( meshPool, 'subgridSshCellTableRange', & - subgridSshCellTableRange ) - call mpas_pool_get_array( meshPool, 'subgridCellBathymetryMin', & - subgridCellBathymetryMin ) - - end if - if (.not. config_zero_drying_velocity) return @@ -457,7 +437,7 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness !$omp parallel !$omp do schedule(runtime) private(i, iEdge, k, divOutFlux, layerThickness) do iCell = 1, nCellsAll -! dryflag = 0 ; + do k = minLevelCell(iCell), maxLevelCell(iCell) divOutFlux = 0.0_RKIND layerThickness = min(layerThicknessProvis(k, iCell), layerThicknessCur(k, iCell)) @@ -475,20 +455,6 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness end do layerThickness = layerThickness + dt * divOutFlux - ! - if ( config_use_subgrid_wetting_drying .and. & - config_pressure_gradient_type == 'ssh_gradient' ) then - layerThicknessCell(iCell) = layerThickness - - call ocn_subgrid_ssh_lookup( layerThicknessCell(iCell), & - subgridWetVolumeCellTable(:,iCell), & - subgridSshCellTableRange(:,iCell), & - bottomDepth(iCell), & - subgridCellBathymetryMin(iCell), & - sshCell(iCell) ) - - end if - ! ! if layer thickness is too small, limit divergence flux outwards with ! opposite velocity @@ -528,18 +494,72 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness !$omp end parallel - ! DW implementation ! - if ( config_use_modify_wetting_drying_implementation ) then + ! DW implementation + ! - For onle layer only + if ( config_use_modify_wetting_drying_implementation .and. & + config_pressure_gradient_type == 'ssh_gradient') then + ! + allocate( layerThicknessCell(nCellsAll) ) + allocate( sshCell(nCellsAll) ) + layerThickness = 0.0_RKIND + sshCell = 0.0_RKIND + + call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) + call mpas_pool_get_array( meshPool, 'bottomDepth', bottomDepth) + + ! for 1-layer subgrid wetting/dring ! + if ( config_use_subgrid_wetting_drying ) then + call mpas_pool_get_array( meshPool, & + 'subgridWetVolumeCellTable', subgridWetVolumeCellTable ) + call mpas_pool_get_array( meshPool, & + 'subgridSshCellTableRange', subgridSshCellTableRange ) + call mpas_pool_get_array( meshPool, & + 'subgridCellBathymetryMin', subgridCellBathymetryMin ) + + end if + + wettingVelocityFactor = 0.0_RKIND ! - if ( config_use_subgrid_wetting_drying .and. & - config_pressure_gradient_type == 'ssh_gradient' ) then - ! - wettingVelocityFactor = 0.0_RKIND - ! - cellLoop: do iCell = 1, nCellsAll - ! if layer thickness is too small, limit divergence flux outwards with - ! opposite velocity - ifDryCell: if ( layerThicknessCell(iCell) <= & + layerLoop: do iCell = 1, nCellsAll + do k = minLevelCell(iCell), maxLevelCell(iCell) + divOutFlux = 0.0_RKIND + layerThickness = min(layerThicknessProvis(k, iCell), layerThicknessCur(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 + layerThickness = layerThickness + dt * divOutFlux + end do + + layerThicknessCell(iCell) = layerThickness + ! + if ( config_use_subgrid_wetting_drying ) then + + call ocn_subgrid_ssh_lookup( layerThicknessCell(iCell), & + subgridWetVolumeCellTable(:,iCell), & + subgridSshCellTableRange(:,iCell), & + bottomDepth(iCell), & + subgridCellBathymetryMin(iCell), & + sshCell(iCell) ) + else + sshCell(iCell) = layerThicknessCell(iCell) - bottomDepth(iCell) + end if + ! + ! Check with Steven if updating hao is needed + end do layerLoop + + cellLoop: do iCell = 1, nCellsAll + ! if layer thickness is too small, limit divergence flux outwards with + ! opposite velocity + ifDryCell: if ( layerThicknessCell(iCell) <= & hCrit + config_drying_safety_height ) then edgeLoop: do i = 1, nEdgesOnCell(iCell) @@ -573,14 +593,14 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness end if end do edgeLoop - end if ifDryCell - end do cellLoop + end if ifDryCell + end do cellLoop - ! Ramping - hRampMin = config_zero_drying_velocity_ramp_hmin - cellLoopRamp: do iCell = 1, nCellsAll - ! - ifRampCell: if ( config_zero_drying_velocity_ramp .and. & + ! Ramping + hRampMin = config_zero_drying_velocity_ramp_hmin + cellLoopRamp: do iCell = 1, nCellsAll + ! + ifRampCell: if ( config_zero_drying_velocity_ramp .and. & (layerThicknessCell(iCell) > & hCrit + config_drying_safety_height) .and. & (layerThicknessCell(iCell) <= config_zero_drying_velocity_ramp_hmax) ) then @@ -616,15 +636,15 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness ! end do edgeLoopRamp ! - end if ifRampCell - ! - end do cellLoopRamp + end if ifRampCell + ! + end do cellLoopRamp ! - end if + ! end if ! end if - if ( config_use_subgrid_wetting_drying .and. & + if ( config_use_modify_wetting_drying_implementation .and. & config_pressure_gradient_type == 'ssh_gradient' ) then deallocate( layerThicknessCell ) ; deallocate( sshCell ) ; From ab33659de87add9af6f17d0e73d0a8ef5e571d4e Mon Sep 17 00:00:00 2001 From: "D. Wirasaet" Date: Wed, 31 Jan 2024 11:52:24 -0800 Subject: [PATCH 13/16] Parabolic bowl setup --- .../mode_init/mpas_ocn_init_buttermilk_bay.F | 2 +- .../mode_init/mpas_ocn_init_parabolic_bowl.F | 132 +++++++++++++----- 2 files changed, 96 insertions(+), 38 deletions(-) diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F index f3d60b97ad37..32eae7feb828 100644 --- a/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F +++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F @@ -495,7 +495,7 @@ subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{ bathymetryMin = maxval( subgridBathymetryValues(1:nSubgridCell) ) bathymetryMax = minval( subgridBathymetryValues(1:nSubgridCell) ) - if ( abs(bathymetryMin - bathymetryMax) > 100.0*eps ) then + if ( abs(bathymetrymin - bathymetrymax) > 100.0*eps ) then subgridSshCellTableRange(1,iCell) = -maxval(subgridBathymetryValues(1:nSubgridCell)) subgridSshCellTableRange(2,iCell) = -minval(subgridBathymetryValues(1:nSubgridCell)) else diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F index c69c2419bb10..fe60cd85290e 100644 --- a/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F +++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F @@ -189,6 +189,7 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ real (kind=RKIND), dimension(:), allocatable:: dsEdge integer:: jj real (kind=RKIND), dimension(:), allocatable:: xSubgridCell, ySubgridCell + real (kind=RKIND):: bathymetryMin, bathymetryMax @@ -449,9 +450,20 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ ! Evaluate bounds of look-up table range !--------------------------------------------------------------- + bathymetryMin = maxval( subgridBathymetryValues(1:nSubgridCell) ) + bathymetryMax = minval( subgridBathymetryValues(1:nSubgridCell) ) - subgridSshCellTableRange(1,iCell) = -maxval(subgridBathymetryValues(1:nSubgridCell)) + config_drying_min_cell_height + eps - subgridSshCellTableRange(2,iCell) = -minval(subgridBathymetryValues(1:nSubgridCell)) + if ( abs(bathymetryMin - bathymetryMax) > 100.0*eps ) then + subgridSshCellTableRange(1,iCell) = -maxval(subgridBathymetryValues(1:nSubgridCell)) + subgridSshCellTableRange(2,iCell) = -minval(subgridBathymetryValues(1:nSubgridCell)) + else + ! flat bathy ! + subgridSshCellTableRange(1,iCell) = -bathymetryMin + subgridSshCellTableRange(2,iCell) = -bathymetryMin + 2.0*config_drying_min_cell_height + endif + +! subgridSshCellTableRange(1,iCell) = -maxval(subgridBathymetryValues(1:nSubgridCell)) + config_drying_min_cell_height + eps +! subgridSshCellTableRange(2,iCell) = -minval(subgridBathymetryValues(1:nSubgridCell)) ! Evaluate subgrid bathymetry !--------------------------------------------------------------- @@ -467,15 +479,26 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ ! Evaluate wet layerThickness average !--------------------------------------------------------------- - call ocn_init_wet_average(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, layerThickness(1,iCell)) +!!$ call ocn_init_wet_average(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, layerThickness(1,iCell)) -!!$ call ocn_init_wet_average_ssh(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, ssh(iCell))!{{{ - - + call ocn_init_wet_average_ssh(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, ssh(iCell))!{{{ enddo - + if (config_use_subgrid_wetting_drying) then + do iCell = 1, nCellsSolve + call ocn_subgrid_ssh_lookup( config_drying_min_cell_height, & + subgridWetVolumeCellTable(:,iCell), & + subgridSshCellTableRange(:,iCell), & + bottomDepth(iCell), & + subgridCellBathymetryMin(iCell), & + subgridSShCellTableRange(3,iCell) ) + + if ( ssh(iCell) < subgridSshCellTableRange(3,iCell) ) then + ssh(iCell) = subgridSshCellTableRange(3,iCell) ; + end if + end do + end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Edges @@ -538,10 +561,18 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ if ( .NOT. config_parabolic_bowl_subgrid_edge_bathymetry_max_pixel ) then ! Evaluate bounds of look-up table range !--------------------------------------------------------------- + bathymetryMin = maxval( subgridBathymetryValues(1:nSubgridCell) ) + bathymetryMax = minval( subgridBathymetryValues(1:nSubgridCell) ) - subgridSshEdgeTableRange(1,iEdge) = -maxval(subgridBathymetryValues(1:nSubgridEdge)) + config_drying_min_cell_height + eps - subgridSshEdgeTableRange(2,iEdge) = -minval(subgridBathymetryValues(1:nSubgridEdge)) + if ( abs(bathymetrymin - bathymetrymax) > 100.0*eps ) then + subgridSshEdgeTableRange(1,iEdge) = -maxval(subgridBathymetryValues(1:nSubgridEdge)) + subgridSshEdgeTableRange(2,iEdge) = -minval(subgridBathymetryValues(1:nSubgridEdge)) + else + ! flat bathy ! + subgridSshEdgeTableRange(1,iCell) = -bathymetryMin + subgridSshEdgeTableRange(2,iCell) = -bathymetryMin + 2.0*config_drying_min_cell_height + endif ! Evaluate subgrid bathymetry !--------------------------------------------------------------- @@ -573,9 +604,17 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ ! Evaluate bounds of look-up table range !--------------------------------------------------------------- - - subgridSshEdgeTableRange(1,iEdge) = -maxval(CellEdgeBathymetryValues(3,1:nSubgridCellEdge)) + config_drying_min_cell_height + eps - subgridSshEdgeTableRange(2,iEdge) = -minval(CellEdgeBathymetryValues(3,1:nSubgridCellEdge)) + bathymetryMin = maxval( cellEdgeBathymetryValues(3,1:nSubgridCellEdge) ) + bathymetryMax = minval( cellEdgeBathymetryValues(3,1:nSubgridCellEdge) ) + + if ( abs(bathymetrymin - bathymetrymax) > 100.0*eps ) then + subgridSshEdgeTableRange(1,iEdge) = -maxval(CellEdgeBathymetryValues(3,1:nSubgridCellEdge)) + subgridSshEdgeTableRange(2,iEdge) = -minval(CellEdgeBathymetryValues(3,1:nSubgridCellEdge)) + else + ! flat bathy ! + subgridSshEdgeTableRange(1,iCell) = -bathymetryMin + subgridSshEdgeTableRange(2,iCell) = -bathymetryMin + 2.0*config_drying_min_cell_height + endif ! Evaluate bounds of look-up table range !--------------------------------------------------------------- @@ -586,11 +625,26 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ ! Vertical integration of wet fraction !--------------------------------------------------------------- call ocn_init_vertical_integration( iEdge, subgridSshEdgeTableRange, & - nSubgridCellEdge, cellEdgeBathymetryValues(3,:), dsEdge, subgridWetVolumeEdgeTable, subgridWetFractionEdgeTable ) + nSubgridCellEdge, cellEdgeBathymetryValues(3,:), dsEdge, & + subgridWetVolumeEdgeTable, subgridWetFractionEdgeTable ) endif end do + ! + ! find an ssh value corresponding to drying_min_cell-height + ! of each edge + if (config_use_subgrid_wetting_drying) then + do iEdge = 1,nEdgesSolve + call ocn_subgrid_ssh_lookup( config_drying_min_cell_height, & + subgridWetVolumeEdgeTable(:,iEdge), & + subgridSshEdgeTableRange(:,iEdge), & + subgridEdgeBathymetryMean(iEdge), & + subgridEdgeBathymetryMin(iEdge), & + subgridSshEdgeTableRange(3,iEdge) ) + end do + end if + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Vertex !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -690,9 +744,9 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ ! Set layer thickness and ssh if (config_use_wetting_drying) then - do iCell = 1, nCellsSolve - ! Set up vertical grid - maxLevelCell(iCell) = nVertLevels ; ! sigma coordinates + do iCell = 1, nCellsSolve + ! Set up vertical grid + maxLevelCell(iCell) = nVertLevels ; ! sigma coordinates end do do iCell = 1, nCellsSolve @@ -703,27 +757,28 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ if (config_use_subgrid_wetting_drying) then - call ocn_subgrid_ssh_lookup(layerThickness(1,iCell),& - subgridWetVolumeCellTable(:,iCell),& - subgridSshCellTableRange(:,iCell),& - bottomDepth(iCell),& - subgridCellBathymetryMin(iCell),& - ssh(iCell)) -! call ocn_subgrid_layer_thickness_lookup(ssh(iCell), & -! subgridWetVolumeCellTable(:,iCell), & -! subgridSshCellTableRange(:,iCell),& -! bottomDepth(iCell),& -! LayerThickness(1,iCell)) +! call ocn_subgrid_ssh_lookup(layerThickness(1,iCell),& +! subgridWetVolumeCellTable(:,iCell),& +! subgridSshCellTableRange(:,iCell),& +! bottomDepth(iCell),& +! subgridCellBathymetryMin(iCell),& +! ssh(iCell)) + + call ocn_subgrid_layer_thickness_lookup(ssh(iCell), & + subgridWetVolumeCellTable(:,iCell), & + subgridSshCellTableRange(:,iCell),& + bottomDepth(iCell),& + LayerThickness(1,iCell)) else - call ocn_init_parabolic_bowl_bathymetry(xCell(iCell),yCell(iCell),bottomDepth(iCell)) - call ocn_init_parabolic_bowl_ssh(xCell(iCell),yCell(iCell),bottomDepth(iCell),ssh(iCell)) + !call ocn_init_parabolic_bowl_bathymetry(xCell(iCell),yCell(iCell),bottomDepth(iCell)) + !call ocn_init_parabolic_bowl_ssh(xCell(iCell),yCell(iCell),bottomDepth(iCell),ssh(iCell)) ssh(iCell) = - bottomDepth(iCell) + & max(ssh(iCell) + bottomDepth(iCell), & - maxLevelCell(iCell)*(config_drying_min_cell_height + eps)) - + maxLevelCell(iCell)*(config_drying_min_cell_height+eps*0.0)) + do k = 1, maxLevelCell(iCell) - layerThickness(k,iCell) = max(config_drying_min_cell_height + eps, & + layerThickness(k,iCell) = max(config_drying_min_cell_height + eps*0.0, & (ssh(iCell) + bottomDepth(iCell))/real(maxLevelCell(iCell),RKIND)) if (layerThickness(k,iCell) < config_drying_min_cell_height) then @@ -745,14 +800,14 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{ block_ptr => block_ptr % next end do - do iCell = 1,nCellsSolve + do iCell = 1,nCellsSolve call ocn_subgrid_layer_thickness_lookup(ssh(iCell), & subgridWetVolumeCellTable(:,iCell), & subgridSshCellTableRange(:,iCell),& bottomDepth(iCell),& subgridLayerThicknessDebug(iCell)) - enddo - ! + enddo + ! deallocate(interfaceLocations) @@ -1147,7 +1202,9 @@ subroutine ocn_init_vertical_integration(iCV, subgridSshTableRange, nSubgridCV, deltaZ = (subgridSshTableRange(2,iCV)-subgridSshTableRange(1,iCV))/real(nSubgridTableLevels-1,RKIND); - subgridWetVolumeTable(1,iCV) = config_drying_min_cell_height + eps +! subgridWetVolumeTable(1,iCV) = config_drying_min_cell_height + eps + + subgridWetVolumeTable(1,iCV) = 0.0_RKIND ssh = subgridSshTableRange(1,iCV) + deltaZ subgridWetVolumeTable(2,iCV) = subgridWetVolumeTable(1,iCV)*sum(subgridAreas(1:nSubgridCV)) subgridWetFractionTable(1,iCV) = 0.0_RKIND @@ -1273,7 +1330,8 @@ subroutine ocn_init_wet_average_ssh( nSubgridCV, subgridBathymetryValues, subgri do tri = 1, nsubgridCV layerThicknessValue = subgridBathymetryValues(tri) + subgridSshValues(tri) - if (layerThicknessValue > config_drying_min_cell_height + eps ) then +! if (layerThicknessValue > config_drying_min_cell_height + eps ) then + if ( layerThicknessValue > config_drying_min_cell_height ) then sshWetAverage = sshWetAverage + subgridSshValues(tri)*subgridAreas(tri) wetArea = wetArea + subgridAreas(tri) endif From 00d2dda783c2e427172d96adba457e9ac51db2d1 Mon Sep 17 00:00:00 2001 From: "D. Wirasaet" Date: Sun, 11 Feb 2024 14:34:14 -0800 Subject: [PATCH 14/16] Check point: updating wetting/drying (add halo update) and parabolic bowl --- components/mpas-ocean/src/Registry.xml | 11 +- .../mpas_ocn_time_integration_rk4.F | 13 ++- .../mode_init/mpas_ocn_init_buttermilk_bay.F | 29 ++--- .../mode_init/mpas_ocn_init_parabolic_bowl.F | 41 +++++-- .../shared/mpas_ocn_diagnostics_variables.F | 11 ++ .../src/shared/mpas_ocn_wetting_drying.F | 101 ++++++++++++------ 6 files changed, 144 insertions(+), 62 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index b7eabd6baf47..06ec49b50949 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -1554,7 +1554,6 @@ - @@ -2930,7 +2929,7 @@ missing_value="FILLVAL" missing_value_mask="cellMask" packages="forwardMode;analysisMode" /> - @@ -3424,6 +3423,14 @@ description="gradient of sea surface height perturbation from self-attraction and loading" packages="tidalPotentialForcingPKG" /> + +