diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 9ce37b56fdfc..ddb169e32538 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -776,6 +776,7 @@ + @@ -1537,6 +1538,9 @@ is the value of that variable from the *previous* time level! + diff --git a/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.cpp b/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.cpp index 8ecff9880441..0cb56962f6b4 100644 --- a/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.cpp +++ b/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.cpp @@ -69,7 +69,7 @@ std::vector indexToTriangleID, std::vector indexToVertexID, vertexToFCell, vertexProcIDs, triangleToFVertex, indexToEdgeID, edgeToFEdge, fVertexToTriangleID, fCellToVertex, iceMarginEdgesLIds, dirichletNodesIDs; std::vector dissipationHeatOnPrisms, velocityOnVertices, velocityOnCells, - elevationData, thicknessData, betaData, bedTopographyData, stiffnessFactorData, effecPressData, muFrictionData, temperatureDataOnPrisms, smbData, thicknessOnCells, bodyForceOnBasalCell; + elevationData, thicknessData, betaData, bedTopographyData, stiffnessFactorData, effecPressData, muFrictionData, bedRoughnessData, bulkFrictionData, basalDebrisData, temperatureDataOnPrisms, smbData, thicknessOnCells, bodyForceOnBasalCell; std::vector isVertexBoundary, isBoundaryEdge; // only needed for creating ASCII mesh @@ -240,7 +240,8 @@ void velocity_solver_init_fo(double const *levelsRatio_F) { void velocity_solver_solve_fo(double const* bedTopography_F, double const* lowerSurface_F, double const* thickness_F, double * beta_F, double const* smb_F, double const* temperature_F, double const* stiffnessFactor_F, - double const* effecPress_F, double const* muFriction_F, + double* effecPress_F, double const* muFriction_F, + double const* bedRoughnessRC_F, double const* bulkFriction_F, double const* basalDebris_F, double* const dirichletVelocityXValue, double* const dirichletVelocitYValue, double* u_normal_F, double* bodyForce_F, double* dissipation_heat_F, double* xVelocityOnCell, double* yVelocityOnCell, double const* deltat, @@ -269,7 +270,7 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower if (!isDomainEmpty) { std::vector > marineBdyExtensionMap; - importFields(marineBdyExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, stiffnessFactor_F, effecPress_F, muFriction_F, temperature_F, smb_F, minThickness); + importFields(marineBdyExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, stiffnessFactor_F, effecPress_F, muFriction_F, bedRoughnessRC_F, bulkFriction_F, basalDebris_F, temperature_F, smb_F, minThickness); std::vector regulThk(thicknessData); for (int index = 0; index < nVertices; index++) @@ -286,7 +287,7 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower Ordering, first_time_step, indexToVertexID, indexToTriangleID, minBeta, regulThk, levelsNormalizedThickness, elevationData, thicknessData, betaData, bedTopographyData, smbData, - stiffnessFactorData, effecPressData, muFrictionData, + stiffnessFactorData, effecPressData, muFrictionData, bedRoughnessData, bulkFrictionData, basalDebrisData, temperatureDataOnPrisms, bodyForceOnBasalCell, dissipationHeatOnPrisms, velocityOnVertices, albany_error, dt); *error=albany_error; @@ -299,6 +300,8 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower } exportBeta(beta_F); + exportEffectivePressure(effecPress_F); + mapVerticesToCells(velocityOnVertices, &velocityOnCells[0], 2, nLayers, Ordering); @@ -1117,6 +1120,7 @@ double signedTriangleAreaOnSphere(const double* x, const double* y, void importFields(std::vector >& marineBdyExtensionMap, double const* bedTopography_F, double const * lowerSurface_F, double const * thickness_F, double const * beta_F, double const* stiffnessFactor_F, double const* effecPress_F, double const* muFriction_F, + double const* bedRoughnessRC_F, double const* bulkFriction_F, double const* basalDebris_F, double const * temperature_F, double const * smb_F, double eps) { int vertexLayerShift = (Ordering == 0) ? 1 : nLayers + 1; @@ -1131,6 +1135,12 @@ void importFields(std::vector >& marineBdyExtensionMap, dou effecPressData.assign(nVertices, 1e10); if (muFriction_F!= 0) muFrictionData.assign(nVertices, 1e10); + if (bedRoughnessRC_F!= 0) + bedRoughnessData.assign(nVertices, 1e10); + if (bulkFriction_F != 0) + bulkFrictionData.assign(nVertices, 1e10); + if (basalDebris_F != 0) + basalDebrisData.assign(nVertices, 1e10); if(temperature_F != 0) temperatureDataOnPrisms.assign(nLayers * nTriangles, 1e10); if (smb_F != 0) @@ -1153,6 +1163,12 @@ void importFields(std::vector >& marineBdyExtensionMap, dou effecPressData[index] = effecPress_F[iCell] / unit_length; if (muFriction_F != 0) muFrictionData[index] = muFriction_F[iCell]; + if (bedRoughnessRC_F != 0) + bedRoughnessData[index] = bedRoughnessRC_F[iCell] / unit_length; + if (bulkFriction_F != 0) + bulkFrictionData[index] = bulkFriction_F[iCell]; + if (basalDebris_F != 0) + basalDebrisData[index] = basalDebris_F[iCell]; } int lElemColumnShift = (Ordering == 1) ? 1 : nTriangles; @@ -1361,6 +1377,18 @@ void exportBeta(double * beta_F) { allToAll (beta_F, sendCellsList_F, recvCellsList_F, 1); } +void exportEffectivePressure(double * effecPress_F) { + std::fill(effecPress_F, effecPress_F + nCells_F, 0.); + for (int index = 0; index < nVertices; index++) { + int fCell = vertexToFCell[index]; + effecPress_F[fCell] = effecPressData[index] * unit_length; + } +#ifdef changeTrianglesOwnership + allToAll (effecPress_F, &sendCellsListReversed, &recvCellsListReversed, 1); +#endif + allToAll (effecPress_F, sendCellsList_F, recvCellsList_F, 1); +} + void createReducedMPI(int nLocalEntities, MPI_Comm& reduced_comm_id) { int numProcs, me; if (reduced_comm_id != MPI_COMM_NULL) @@ -1630,6 +1658,7 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep double const* surfaceAirTemperature_F, double const* basalHeatFlux_F, double const* stiffnessFactor_F, double const* effecPress_F, double const* muFriction_F, + double const* bedRoughnessRC_F, double const* bulkFriction_F, double const* basalDebris_F, double const* thickness_F, double const* thicknessUncertainty_F, double const* smb_F, double const* smbUncertainty_F, double const* bmb_F, double const* bmbUncertainty_F, @@ -1715,7 +1744,8 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep std::vector > marineBdyExtensionMap; // local map to be created by importFields importFields(marineBdyExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, - stiffnessFactor_F, effecPress_F, muFriction_F, temperature_F, smb_F, minThickness); + stiffnessFactor_F, effecPress_F, muFriction_F, bedRoughnessRC_F, bulkFriction_F, basalDebris_F, + temperature_F, smb_F, minThickness); import2DFieldsObservations(marineBdyExtensionMap, thicknessUncertainty_F, @@ -1760,6 +1790,7 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep write_ascii_mesh_field(betaData, "beta"); write_ascii_mesh_field(muFrictionData, "mu"); write_ascii_mesh_field(muLogData, "mu_log"); + write_ascii_mesh_field(bedRoughnessData, "bed_roughness"); write_ascii_mesh_field(stiffnessFactorData, "stiffening_factor"); write_ascii_mesh_field(stiffnessLogData, "stiffening_factor_log"); @@ -1780,6 +1811,7 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep write_ascii_mesh_field(effecPressData, "effective_pressure"); + // These two fields are more complicated than the others so cannot use the function to write out std::cout << "Writing temperature.ascii." << std::endl; diff --git a/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.hpp b/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.hpp index fb8a907216ac..31db53c1686e 100644 --- a/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.hpp +++ b/components/mpas-albany-landice/src/mode_forward/Interface_velocity_solver.hpp @@ -90,7 +90,7 @@ void velocity_solver_solve_l1l2(double const* lowerSurface_F, void velocity_solver_solve_fo(double const* bedTopography_F, double const* lowerSurface_F, double const* thickness_F, double * beta_F, double const* smb_F, double const* temperature_F, double const* stiffnessFactor_F, - double const* effecPress_F, double const* muFriction_F, + double* effecPress_F, double const* muFriction_F, double const* bedRoughnessRC_F, double const* bulkFriction_F, double const* basalDebris_F, double* const dirichletVelocityXValue = 0, double* const dirichletVelocitYValue = 0, double* u_normal_F = 0, double* bodyForce_F = 0, double* dissipation_heat_F = 0, double* xVelocityOnCell = 0, double* yVelocityOnCell = 0, double const * deltat = 0, @@ -134,6 +134,7 @@ void write_ascii_mesh(int const* indexToCellID_F, double const* surfaceAirTemperature_F, double const* basalHeatFlux_F, double const* stiffnessFactor_F, double const* effecPress_F, double const* muFriction_F, + double const* bedRoughnessRC_F, double const* bulkFriction_F, double const* basalDebris_F, double const* thickness_F, double const* thicknessUncertainty_F, double const* smb_F, double const* smbUncertainty_F, double const* bmb_F, double const* bmbUncertainty_F, @@ -161,8 +162,11 @@ extern void velocity_solver_solve_fo__(int nLayers, int nGlobalVertices, const std::vector& bedTopographyData, const std::vector& smbData, const std::vector& stiffnessFactorData, - const std::vector& effecPressData, + std::vector& effecPressData, const std::vector& muFrictionData, + const std::vector& bedRoughnessData, + const std::vector& bulkFrictionData, + const std::vector& basalDebrisData, const std::vector& temperatureDataOnPrisms, std::vector& bodyForceOnBasalCell, std::vector& dissipationHeatOnPrisms, @@ -207,7 +211,8 @@ void createReducedMPI(int nLocalEntities, MPI_Comm& reduced_comm_id); void importFields(std::vector >& marineBdyExtensionMap, double const* bedTopography_F, double const* lowerSurface_F, double const* thickness_F, - double const* beta_F = 0, double const* stiffnessFactor_F = 0, double const* effecPress_F = 0, double const* muFriction_F = 0, double const* temperature_F = 0, double const* smb_F = 0, double eps = 0); + double const* beta_F = 0, double const* stiffnessFactor_F = 0, double const* effecPress_F = 0, double const* muFriction_F = 0, double const* bedRoughness_F = 0, + double const* bulkFriction_F = 0, double const* basalDebris_F = 0, double const* temperature_F = 0, double const* smb_F = 0, double eps = 0); void import2DFieldsObservations(std::vector >& marineBdyExtensionMap, double const * lowerSurface_F, @@ -231,6 +236,7 @@ std::vector extendMaskByOneLayer(int const* verticesMask_F); void exportDissipationHeat(double * dissipationHeat_F); void exportBodyForce(double * bodyForce_F); void exportBeta(double * beta_F); +void exportEffectivePressure(double * effecPress_F); void get_prism_velocity_on_FEdges(double* uNormal, const std::vector& velocityOnCells, diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F index 077ff4c1773f..b5eb319655cd 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F @@ -452,6 +452,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro real (kind=RKIND), dimension(:), pointer :: effectivePressure real (kind=RKIND), dimension(:), pointer :: effectivePressureLimited real (kind=RKIND), dimension(:), pointer :: muFriction + real (kind=RKIND), dimension(:), pointer :: bedRoughnessRC real (kind=RKIND), pointer :: deltat integer, dimension(:), pointer :: vertexMask, cellMask, edgeMask integer, dimension(:,:), pointer :: dirichletVelocityMask @@ -516,6 +517,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro call mpas_pool_get_array(velocityPool, 'drivingStressVert', drivingStressVert) call mpas_pool_get_array(velocityPool, 'drivingStress', drivingStress) call mpas_pool_get_array(velocityPool, 'muFriction', muFriction) + call mpas_pool_get_array(velocityPool, 'bedRoughnessRC', bedRoughnessRC) call mpas_pool_get_array(velocityPool, 'anyDynamicVertexMaskChanged', anyDynamicVertexMaskChanged) call mpas_pool_get_array(velocityPool, 'dirichletMaskChanged', dirichletMaskChanged) call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMask, timeLevel = 1) @@ -598,6 +600,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro call velocity_solver_solve_FO(bedTopography, lowerSurface, thickness, & betaSolve, sfcMassBal, temperature, stiffnessFactor, & effectivePressureLimited, muFriction, & + bedRoughnessRC, & uReconstructX, uReconstructY, & ! Dirichlet boundary values to apply where dirichletVelocityMask=1 normalVelocity, drivingStressVert, dissipationVertexField % array, uReconstructX, uReconstructY, & ! return values deltat, albanyVelocityError) ! return values @@ -787,6 +790,7 @@ subroutine li_velocity_external_write_albany_mesh(domain) real (kind=RKIND), dimension(:), pointer :: stiffnessFactor real (kind=RKIND), dimension(:), pointer :: effectivePressure real (kind=RKIND), dimension(:), pointer :: muFriction + real (kind=RKIND), dimension(:), pointer :: bedRoughnessRC integer, dimension(:,:), pointer :: dirichletVelocityMask type (mpas_pool_type), pointer :: meshPool, geometryPool, thermalPool, observationsPool, velocityPool, scratchPool, hydroPool real (kind=RKIND), pointer :: config_sea_level, config_ice_density, config_ocean_density @@ -848,6 +852,7 @@ subroutine li_velocity_external_write_albany_mesh(domain) ! Velocity variables call mpas_pool_get_array(velocityPool, 'beta', beta) call mpas_pool_get_array(velocityPool, 'muFriction', muFriction) + call mpas_pool_get_array(velocityPool, 'bedRoughnessRC', bedRoughnessRC) call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMask, timeLevel = 1) call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor) @@ -901,6 +906,7 @@ subroutine li_velocity_external_write_albany_mesh(domain) surfaceAirTemperature, basalHeatFlux, & stiffnessFactor, & effectivePressure, muFriction, & + bedRoughnessRC, & thickness, thicknessUncertainty, & sfcMassBal, sfcMassBalUncertainty, & floatingBasalMassBal, floatingBasalMassBalUncertainty, &