@@ -117,6 +117,8 @@ subroutine li_SGH_init(domain, err)
117117 logical , pointer :: config_do_restart
118118 real (kind= RKIND), pointer :: config_sea_level
119119 integer , pointer :: config_num_halos
120+ integer , pointer :: nCells
121+ integer :: iCell
120122 integer :: err_tmp
121123
122124
@@ -202,25 +204,28 @@ subroutine li_SGH_init(domain, err)
202204 call mpas_pool_get_subpool(block % structs, ' hydro' , hydroPool)
203205 call mpas_pool_get_subpool(block % structs, ' geometry' , geometryPool)
204206
207+ call mpas_pool_get_dimension(hydroPool, ' nCells' , nCells)
205208 call mpas_pool_get_array(hydroPool, ' waterPressure' , waterPressure)
206209 call mpas_pool_get_array(hydroPool, ' hydropotential' , hydropotential)
207210 call mpas_pool_get_array(hydroPool, ' iceThicknessHydro' , iceThicknessHydro)
208211 call mpas_pool_get_array(geometryPool, ' cellMask' , cellMask)
209212 call mpas_pool_get_array(geometryPool, ' bedTopography' , bedTopography)
210213
211214 waterPressure = max(0.0_RKIND, waterPressure)
212- where (li_mask_is_grounded_ice(cellMask))
213- waterPressure = min(waterPressure, rhoi * gravity * iceThicknessHydro)
214- end where
215-
216- ! set pressure and hydropotential correctly on ice-free land and in ocean
217- where ((.not. (li_mask_is_grounded_ice(cellMask))) .and. (bedTopography > config_sea_level))
218- waterPressure = 0.0_RKIND
219- hydropotential = rho_water * gravity * bedTopography
220- elsewhere ((.not. (li_mask_is_grounded_ice(cellMask))) .and. (bedTopography <= config_sea_level))
221- waterPressure = gravity * rhoo * (config_sea_level - bedTopography)
222- hydropotential = 0.0_RKIND
223- end where
215+ do iCell = 1, nCells
216+ if (li_mask_is_grounded_ice(cellMask(iCell))) then
217+ waterPressure(iCell) = min(waterPressure(iCell), rhoi * gravity * iceThicknessHydro(iCell))
218+ else
219+ ! set pressure and hydropotential correctly on ice-free land and in ocean
220+ if (bedTopography(iCell) > config_sea_level) then
221+ waterPressure(iCell) = 0.0_RKIND
222+ hydropotential(iCell) = rho_water * gravity * bedTopography(iCell)
223+ else
224+ waterPressure(iCell) = gravity * rhoo * (config_sea_level - bedTopography(iCell))
225+ hydropotential(iCell) = 0.0_RKIND
226+ endif
227+ endif
228+ enddo
224229
225230 ! Initialize diagnostic pressure variables
226231 call calc_pressure_diag_vars(block, err_tmp)
@@ -1762,6 +1767,8 @@ subroutine calc_pressure_diag_vars(block, err)
17621767 real (kind= RKIND), dimension (:), pointer :: iceThicknessHydro
17631768 integer , dimension (:), pointer :: cellMask
17641769 real (kind= RKIND), pointer :: config_sea_level
1770+ integer , pointer :: nCells
1771+ integer :: iCell
17651772
17661773 err = 0
17671774
@@ -1773,6 +1780,7 @@ subroutine calc_pressure_diag_vars(block, err)
17731780 call mpas_pool_get_config(liConfigs, ' config_sea_level' , config_sea_level)
17741781 call mpas_pool_get_config(liConfigs, ' config_ocean_density' , rhoo)
17751782
1783+ call mpas_pool_get_dimension(hydroPool, ' nCells' , nCells)
17761784 call mpas_pool_get_array(hydroPool, ' effectivePressureSGH' , effectivePressureSGH)
17771785 call mpas_pool_get_array(hydroPool, ' waterPressure' , waterPressure)
17781786 call mpas_pool_get_array(geometryPool, ' bedTopography' , bedTopography)
@@ -1789,11 +1797,16 @@ subroutine calc_pressure_diag_vars(block, err)
17891797 end where
17901798
17911799 hydropotentialBase = rho_water * gravity * bedTopography + waterPressure
1792- where ((.not. li_mask_is_grounded_ice(cellMask)) .and. (bedTopography <= config_sea_level))
1793- hydropotentialBase = 0.0_RKIND !zero hydropotential in ocean
1794- elsewhere ((.not. li_mask_is_grounded_ice(cellMask)) .and. (bedTopography > config_sea_level))
1795- hydropotentialBase = rho_water * gravity * bedTopography ! for completeness, but won' t matter with one-side slope calculations on terrestrial boundaries
1796- end where
1800+ do iCell = 1 , nCells
1801+ if (.not. li_mask_is_grounded_ice(cellMask(iCell))) then
1802+ if (bedTopography(iCell) <= config_sea_level) then
1803+ hydropotentialBase(iCell) = 0.0_RKIND !zero hydropotential in ocean
1804+ else
1805+ hydropotentialBase(iCell) = rho_water * gravity * bedTopography(iCell)
1806+ ! for completeness, but won' t matter with one-side slope calculations on terrestrial boundaries
1807+ endif
1808+ endif
1809+ enddo
17971810
17981811 ! hydropotential with water thickness
17991812 hydropotential = hydropotentialBase + rho_water * gravity * waterThickness
0 commit comments