From 57c9dfa96c124ccfea53bc2670e6e23039f65aad Mon Sep 17 00:00:00 2001 From: Elizabeth Hunke Date: Sat, 23 Sep 2023 22:00:15 -0500 Subject: [PATCH 1/5] optionally call icepack routines; fix threading bug --- components/mpas-seaice/driver/ice_comp_mct.F | 54 +++++++++++++++---- .../src/shared/mpas_seaice_icepack.F | 12 ++--- 2 files changed, 48 insertions(+), 18 deletions(-) diff --git a/components/mpas-seaice/driver/ice_comp_mct.F b/components/mpas-seaice/driver/ice_comp_mct.F index 214fca9caeb8..7a7acaa3595d 100644 --- a/components/mpas-seaice/driver/ice_comp_mct.F +++ b/components/mpas-seaice/driver/ice_comp_mct.F @@ -48,6 +48,8 @@ module ice_comp_mct use seaice_analysis_driver use seaice_column, only: seaice_column_reinitialize_fluxes, & !colpkg seaice_column_coupling_prep + use seaice_icepack, only: seaice_icepack_reinitialize_fluxes, & + seaice_icepack_coupling_prep use seaice_constants, only: coupleAlarmID, & seaiceFreshWaterFreezingPoint, & seaiceDensityIce, & @@ -69,6 +71,7 @@ module ice_comp_mct use seaice_error, only: seaice_check_critical_error use ice_colpkg, only: colpkg_sea_freezing_temperature + use icepack_intfc, only: icepack_sea_freezing_temperature ! ! !PUBLIC MEMBER FUNCTIONS: @@ -745,8 +748,12 @@ end subroutine xml_stream_get_attributes call mpas_reset_clock_alarm(domain % clock, coupleAlarmID, ierr=ierr) ! Coupling prep - call seaice_column_coupling_prep(domain) - + call MPAS_pool_get_config(domain % configs, "config_column_physics_type", tempCharConfig) + if (trim(tempCharConfig) == "icepack") then + call seaice_icepack_coupling_prep(domain) + else if (trim(tempCharConfig) == "column_package") then + call seaice_column_coupling_prep(domain) + endif ! config_column_physics_type !----------------------------------------------------------------------- ! @@ -1053,7 +1060,9 @@ subroutine ice_run_mct( EClock, cdata_i, x2i_i, i2x_i)!{{{ logical :: streamActive, debugOn logical, pointer :: config_write_output_on_startup logical, save :: first=.true. - character (len=StrKIND), pointer :: config_restart_timestamp_name + character (len=StrKIND), pointer :: & + config_restart_timestamp_name, & + config_column_physics_type real(kind=RKIND), pointer :: & dayOfNextShortwaveCalculation ! needed for CESM like coupled simulations @@ -1069,6 +1078,7 @@ subroutine ice_run_mct( EClock, cdata_i, x2i_i, i2x_i)!{{{ if (debugOn) call mpas_log_write("=== Beginning ice_run_mct ===") call mpas_pool_get_config(domain % configs, 'config_restart_timestamp_name', config_restart_timestamp_name) + call MPAS_pool_get_config(domain % configs, "config_column_physics_type", config_column_physics_type) ! Setup log information. call shr_file_getLogUnit (shrlogunit) @@ -1076,7 +1086,11 @@ subroutine ice_run_mct( EClock, cdata_i, x2i_i, i2x_i)!{{{ call shr_file_setLogUnit (iceLogUnit) ! reinitialize fluxes - call seaice_column_reinitialize_fluxes(domain) + if (trim(config_column_physics_type) == "icepack") then + call seaice_icepack_reinitialize_fluxes(domain) + else if (trim(config_column_physics_type) == "column_package") then + call seaice_column_reinitialize_fluxes(domain) + endif ! config_column_physics_type ! Import state from coupler call ice_import_mct(x2i_i, ierr) @@ -1803,6 +1817,7 @@ subroutine ice_import_mct(x2i_i, errorCode)!{{{ config_use_column_biogeochemistry character(len=strKIND), pointer :: & + config_column_physics_type, & config_thermodynamics_type, & config_ocean_surface_type @@ -1913,6 +1928,7 @@ subroutine ice_import_mct(x2i_i, errorCode)!{{{ do while(associated(block_ptr)) configs => block_ptr % configs + call mpas_pool_get_config(configs, "config_column_physics_type", config_column_physics_type) call mpas_pool_get_config(configs, "config_thermodynamics_type", config_thermodynamics_type) call mpas_pool_get_config(configs, "config_ocean_surface_type", config_ocean_surface_type) call mpas_pool_get_config(configs, "config_use_aerosols", config_use_aerosols) @@ -1988,7 +2004,11 @@ subroutine ice_import_mct(x2i_i, errorCode)!{{{ seaSurfaceTemperature(i) = x2i_i % rAttr(index_x2i_So_t, n) seaSurfaceSalinity(i) = x2i_i % rAttr(index_x2i_So_s, n) - seaFreezingTemperature(i) = colpkg_sea_freezing_temperature(seaSurfaceSalinity(i)) + if (trim(config_column_physics_type) == "icepack") then + seaFreezingTemperature(i) = icepack_sea_freezing_temperature(seaSurfaceSalinity(i)) + else if (trim(config_column_physics_type) == "column_package") then + seaFreezingTemperature(i) = colpkg_sea_freezing_temperature(seaSurfaceSalinity(i)) + endif ! config_column_physics_type uOceanVelocity(i) = x2i_i % rAttr(index_x2i_So_u, n) vOceanVelocity(i) = x2i_i % rAttr(index_x2i_So_v, n) @@ -2012,8 +2032,7 @@ subroutine ice_import_mct(x2i_i, errorCode)!{{{ ! coupling step as a freshwater and salt flux. This step is required to balance mass ! and heat with the ocean. - call frazil_mass(freezingMeltingPotential(i), frazilMassFluxRev, seaSurfaceSalinity(i), & - config_thermodynamics_type) + call frazil_mass(freezingMeltingPotential(i), frazilMassFluxRev, seaSurfaceSalinity(i)) frazilMassAdjust(i) = frazilMassFlux-frazilMassFluxRev @@ -2754,8 +2773,7 @@ end subroutine basal_pressure!}}} ! !IROUTINE: frazil_mass ! ! !INTERFACE - subroutine frazil_mass(freezingPotential, frazilMassFlux, seaSurfaceSalinity, & - config_thermodynamics_type) + subroutine frazil_mass(freezingPotential, frazilMassFlux, seaSurfaceSalinity) ! ! !DESCRIPTION: ! Calculate frazil mass based on on the sea surface salinity, and frazil heat flux @@ -2769,14 +2787,21 @@ subroutine frazil_mass(freezingPotential, frazilMassFlux, seaSurfaceSalinity, & liquidus_temperature_mush, & enthalpy_mush +!echmod - not implemented yet +! use icepack_intfc, only: & +! icepack_liquidus_temperature_mush, & +! icepack_enthalpy_mush + ! !INPUT PARAMETERS: real (kind=RKIND), intent(in) :: freezingPotential real (kind=RKIND), intent(in) :: seaSurfaceSalinity - character(len=strKIND), intent(in) :: config_thermodynamics_type ! !OUTPUT PARAMETERS: real (kind=RKIND), intent(out) :: frazilMassFlux + character(len=strKIND), pointer :: config_thermodynamics_type + character(len=strKIND), pointer :: config_column_physics_type + !EOP !BOC !----------------------------------------------------------------------- @@ -2790,6 +2815,8 @@ subroutine frazil_mass(freezingPotential, frazilMassFlux, seaSurfaceSalinity, & qi0new, & vi0new + call MPAS_pool_get_config(domain % configs, "config_thermodynamics_type", config_thermodynamics_type) + call MPAS_pool_get_config(domain % configs, "config_column_physics_type", config_column_physics_type) if (freezingPotential > 0.0_RKIND) then @@ -2799,8 +2826,15 @@ subroutine frazil_mass(freezingPotential, frazilMassFlux, seaSurfaceSalinity, & else Si0new = seaSurfaceSalinity**2 / (4.0_RKIND*seaiceFrazilSalinityReduction) endif +!echmod - not yet implemented in icepack - add these functions to icepack_intfc with prefix icepack_ +! if (trim(config_column_physics_type) == "icepack") then +! Ti = icepack_liquidus_temperature_mush(Si0new/seaiceFrazilIcePorosity) +! qi0new = icepack_enthalpy_mush(Ti, Si0new) +! else if (trim(config_column_physics_type) == "column_package") then Ti = liquidus_temperature_mush(Si0new/seaiceFrazilIcePorosity) qi0new = enthalpy_mush(Ti, Si0new) +! endif ! config_column_physics_type + else qi0new = -seaiceDensityIce*seaiceLatentHeatMelting endif ! ktherm diff --git a/components/mpas-seaice/src/shared/mpas_seaice_icepack.F b/components/mpas-seaice/src/shared/mpas_seaice_icepack.F index 2ba7e93222cf..10e971670e18 100644 --- a/components/mpas-seaice/src/shared/mpas_seaice_icepack.F +++ b/components/mpas-seaice/src/shared/mpas_seaice_icepack.F @@ -23,9 +23,6 @@ module seaice_icepack use seaice_error -! use ice_kinds_mod, only: & !colpkg -! char_len_long - implicit none private @@ -195,8 +192,6 @@ module seaice_icepack real(kind=RKIND), dimension(:), allocatable :: & tracerArrayCell - ! warnings string kind -! integer, parameter :: strKINDWarnings = char_len_long integer, parameter :: strKINDWarnings = strKIND contains @@ -3160,8 +3155,9 @@ subroutine column_radiation(domain, clock, lInitialization) ! snow grain radius array allocate(snow_grain_radius(nSnowLayers,nCategories)) - !$omp parallel do default(shared) firstprivate(aerosolsArray,index_shortwaveAerosol) & - !$omp& private(iCategory,iAerosol,lonCellColumn) + !$omp parallel do default(shared) & + !$omp& firstprivate(aerosolsArray,index_shortwaveAerosol,snow_grain_radius) & + !$omp& private(iCategory,iAerosol,lonCellColumn,iSnowLayer) do iCell = 1, nCellsSolve ! set aerosols array @@ -5887,7 +5883,7 @@ subroutine seaice_init_icepack_constants() call icepack_configure() call seaice_icepack_write_warnings(icepack_warnings_aborted()) -!echmod - BFB but 20% slower when this is set in seaice_constants instead??? +!echmod - set this in seaice_constants instead seaicePi = pi iceAreaMinimum = seaicePuny From 34fc863a62952a97d7bf0bc51aa54cb403e69f72 Mon Sep 17 00:00:00 2001 From: Elizabeth Hunke Date: Sun, 24 Sep 2023 13:42:11 -0500 Subject: [PATCH 2/5] optionally call icepack routines; clean up constants --- .../mpas_seaice_temperatures.F | 44 +++++++++++++++++++ .../src/shared/mpas_seaice_constants.F | 15 ++----- .../src/shared/mpas_seaice_prescribed.F | 28 ++++++++++-- 3 files changed, 73 insertions(+), 14 deletions(-) diff --git a/components/mpas-seaice/src/analysis_members/mpas_seaice_temperatures.F b/components/mpas-seaice/src/analysis_members/mpas_seaice_temperatures.F index 587b78bbe6f6..9c6136e807d8 100644 --- a/components/mpas-seaice/src/analysis_members/mpas_seaice_temperatures.F +++ b/components/mpas-seaice/src/analysis_members/mpas_seaice_temperatures.F @@ -222,6 +222,10 @@ subroutine seaice_compute_temperatures(domain, instance, timeLevel, err)!{{{ colpkg_ice_temperature, & colpkg_snow_temperature + use icepack_intfc, only: & + icepack_ice_temperature, & + icepack_snow_temperature + !----------------------------------------------------------------- ! ! input variables @@ -261,6 +265,9 @@ subroutine seaice_compute_temperatures(domain, instance, timeLevel, err)!{{{ tracersPool, & temperaturesAMPool + character(len=strKIND), pointer :: & + config_column_physics_type + real(kind=RKIND), dimension(:,:,:), pointer :: & iceAreaCategory, & snowVolumeCategory, & @@ -287,6 +294,8 @@ subroutine seaice_compute_temperatures(domain, instance, timeLevel, err)!{{{ block => domain % blocklist do while (associated(block)) + call MPAS_pool_get_config(block % configs, "config_column_physics_type", config_column_physics_type) + call MPAS_pool_get_subpool(block % structs, 'temperaturesAM', temperaturesAMPool) call MPAS_pool_get_subpool(block % structs, 'tracers', tracersPool) @@ -309,6 +318,39 @@ subroutine seaice_compute_temperatures(domain, instance, timeLevel, err)!{{{ snowTemperature = 0.0_RKIND ! compute temperatures + + if (trim(config_column_physics_type) == "icepack") then + + do iCell = 1, nCellsSolve + do iCategory = 1, nCategories + + ! check if ice present + if (iceAreaCategory(1,iCategory,iCell) > 1e-11_RKIND) then + + ! ice layers + do iIceLayer = 1, nIceLayers + iceTemperature(iIceLayer, iCategory, iCell) = & + icepack_ice_temperature(iceEnthalpy(iIceLayer,iCategory,iCell), & + iceSalinity(iIceLayer,iCategory,iCell)) + enddo ! iIceLayer + + ! snow layers + if (snowVolumeCategory(1,iCategory,iCell) > 1e-11_RKIND) then + + do iSnowLayer = 1, nSnowLayers + snowTemperature(iSnowLayer, iCategory, iCell) = & + icepack_snow_temperature(snowEnthalpy(iSnowLayer,iCategory,iCell)) + enddo ! iIceLayer + + endif ! snowVolumeCategory + + endif ! iceAreaCategory + + enddo ! iCategory + enddo ! iCell + + else if (trim(config_column_physics_type) == "column_package") then + do iCell = 1, nCellsSolve do iCategory = 1, nCategories @@ -337,6 +379,8 @@ subroutine seaice_compute_temperatures(domain, instance, timeLevel, err)!{{{ enddo ! iCategory enddo ! iCell + endif ! config_column_physics_type + block => block % next enddo diff --git a/components/mpas-seaice/src/shared/mpas_seaice_constants.F b/components/mpas-seaice/src/shared/mpas_seaice_constants.F index 4060e8fe6939..9fe7e151f5c9 100644 --- a/components/mpas-seaice/src/shared/mpas_seaice_constants.F +++ b/components/mpas-seaice/src/shared/mpas_seaice_constants.F @@ -50,13 +50,7 @@ module seaice_constants seaiceLatentHeatVaporization = SHR_CONST_LATVAP ,&! latent heat, vaporization freshwater (J/kg) seaiceLatentHeatMelting = SHR_CONST_LATICE ,&! latent heat of melting of fresh ice (J/kg) seaiceReferenceSalinity = SHR_CONST_ICE_REF_SAL,&! ice reference salinity (ppt) - seaiceSnowPatchiness = 0.005_RKIND ,&! parameter for fractional snow area (m) - -#ifdef RASM_MODS - seaiceIceOceanDragCoefficient = 0.00962_RKIND ! ice-ocn drag coefficient for RASM as temporary measure -#else - seaiceIceOceanDragCoefficient = 0.00536_RKIND ! ice-ocn drag coefficient -#endif + seaiceSnowPatchiness = 0.005_RKIND ! parameter for fractional snow area (m) ! R_gC2molC = SHR_CONST_MWC ,&! molar mass of carbon @@ -92,9 +86,7 @@ module seaice_constants seaiceLatentHeatMelting = seaiceLatentHeatSublimation & ! latent heat of melting of fresh ice (J/kg) - seaiceLatentHeatVaporization, & seaiceReferenceSalinity = 4._RKIND ,&! ice reference salinity (ppt) - seaiceSnowPatchiness = 0.02_RKIND ,&! parameter for fractional snow area (m) - - seaiceIceOceanDragCoefficient = 0.00536_RKIND ! ice-ocn drag coefficient + seaiceSnowPatchiness = 0.02_RKIND ! parameter for fractional snow area (m) ! R_gC2molC = 12.0107_RKIND, & ! molar mass of carbon #endif @@ -144,7 +136,8 @@ module seaice_constants ! dynamics constants real(kind=RKIND), public :: & - seaiceIceStrengthConstantHiblerP = 2.75e4_RKIND ,&! P* constant in Hibler strength formula + seaiceIceOceanDragCoefficient = 0.00536_RKIND, &! ice-ocn drag coefficient + seaiceIceStrengthConstantHiblerP = 2.75e4_RKIND , &! P* constant in Hibler strength formula seaiceIceStrengthConstantHiblerC = 20._RKIND ! C* constant in Hibler strength formula ! minimum sea ice area diff --git a/components/mpas-seaice/src/shared/mpas_seaice_prescribed.F b/components/mpas-seaice/src/shared/mpas_seaice_prescribed.F index fe6c16dca6f3..57d1daba2eeb 100644 --- a/components/mpas-seaice/src/shared/mpas_seaice_prescribed.F +++ b/components/mpas-seaice/src/shared/mpas_seaice_prescribed.F @@ -135,10 +135,19 @@ subroutine seaice_run_prescribed_ice(domain) colpkg_enthalpy_ice, & colpkg_salinity_profile - use seaice_column, only: & +! use icepack_intfc, only: & +! icepack_enthalpy_snow, & +! colpkg_enthalpy_ice, & !echmod - in icepack_init_trcr +! colpkg_salinity_profile !echmod - in icepack_init_thermo + + use seaice_column, only: & ! colpkg seaice_column_reinitialize_fluxes, & seaice_column_aggregate + use seaice_icepack, only: & + seaice_icepack_reinitialize_fluxes, & + seaice_icepack_aggregate + type (domain_type), intent(inout) :: & domain !< Input/Output: @@ -146,6 +155,9 @@ subroutine seaice_run_prescribed_ice(domain) config_use_prescribed_ice, & config_use_prescribed_ice_forcing + character(len=strKIND), pointer :: & + config_column_physics_type + type(block_type), pointer :: & blockPtr @@ -199,6 +211,7 @@ subroutine seaice_run_prescribed_ice(domain) iIceLayer, & iSnowLayer + call MPAS_pool_get_config(domain % configs, "config_column_physics_type", config_column_physics_type) call mpas_pool_get_config(domain % configs, "config_use_prescribed_ice", config_use_prescribed_ice) if (config_use_prescribed_ice) then @@ -288,6 +301,7 @@ subroutine seaice_run_prescribed_ice(domain) temperatureGradient = seaFreezingTemperature(iCell) - surfaceTemperature(1,iCategory,iCell) +!echmod - add icepack option for ice and snow quantities ! ice quantities do iIceLayer = 1, nIceLayers @@ -338,7 +352,11 @@ subroutine seaice_run_prescribed_ice(domain) enddo ! aggregate tracers - call seaice_column_aggregate(domain) + if (trim(config_column_physics_type) == "icepack") then + call seaice_icepack_aggregate(domain) + else if (trim(config_column_physics_type) == "column_package") then + call seaice_column_aggregate(domain) + endif ! config_column_physics_type ! set non-computed fluxes, ice velocities, ice-ocn stresses to zero blockPtr => domain % blocklist @@ -363,7 +381,11 @@ subroutine seaice_run_prescribed_ice(domain) enddo ! reinitialize fluxes - call seaice_column_reinitialize_fluxes(domain) + if (trim(config_column_physics_type) == "icepack") then + call seaice_icepack_reinitialize_fluxes(domain) + else if (trim(config_column_physics_type) == "column_package") then + call seaice_column_reinitialize_fluxes(domain) + endif ! config_column_physics_type endif ! prescribed ice mode From 3e94f62222dc973678dec3bff363a679f3d29d54 Mon Sep 17 00:00:00 2001 From: Elizabeth Hunke Date: Mon, 25 Sep 2023 16:52:04 -0500 Subject: [PATCH 3/5] simplify icepack interface in MPAS-SI driver modules --- components/mpas-seaice/driver/ice_comp_mct.F | 23 +++--- .../src/shared/mpas_seaice_forcing.F | 6 +- .../src/shared/mpas_seaice_icepack.F | 77 +------------------ .../src/shared/mpas_seaice_initialize.F | 32 ++++---- .../src/shared/mpas_seaice_prescribed.F | 31 ++++++-- 5 files changed, 61 insertions(+), 108 deletions(-) diff --git a/components/mpas-seaice/driver/ice_comp_mct.F b/components/mpas-seaice/driver/ice_comp_mct.F index 7a7acaa3595d..0e0ce90a224d 100644 --- a/components/mpas-seaice/driver/ice_comp_mct.F +++ b/components/mpas-seaice/driver/ice_comp_mct.F @@ -2787,10 +2787,9 @@ subroutine frazil_mass(freezingPotential, frazilMassFlux, seaSurfaceSalinity) liquidus_temperature_mush, & enthalpy_mush -!echmod - not implemented yet -! use icepack_intfc, only: & -! icepack_liquidus_temperature_mush, & -! icepack_enthalpy_mush + use icepack_intfc, only: & + icepack_liquidus_temperature, & + icepack_enthalpy_mush ! !INPUT PARAMETERS: real (kind=RKIND), intent(in) :: freezingPotential @@ -2826,14 +2825,14 @@ subroutine frazil_mass(freezingPotential, frazilMassFlux, seaSurfaceSalinity) else Si0new = seaSurfaceSalinity**2 / (4.0_RKIND*seaiceFrazilSalinityReduction) endif -!echmod - not yet implemented in icepack - add these functions to icepack_intfc with prefix icepack_ -! if (trim(config_column_physics_type) == "icepack") then -! Ti = icepack_liquidus_temperature_mush(Si0new/seaiceFrazilIcePorosity) -! qi0new = icepack_enthalpy_mush(Ti, Si0new) -! else if (trim(config_column_physics_type) == "column_package") then - Ti = liquidus_temperature_mush(Si0new/seaiceFrazilIcePorosity) - qi0new = enthalpy_mush(Ti, Si0new) -! endif ! config_column_physics_type + + if (trim(config_column_physics_type) == "icepack") then + Ti = icepack_liquidus_temperature(Si0new/seaiceFrazilIcePorosity) + qi0new = icepack_enthalpy_mush(Ti, Si0new) + else if (trim(config_column_physics_type) == "column_package") then + Ti = liquidus_temperature_mush(Si0new/seaiceFrazilIcePorosity) + qi0new = enthalpy_mush(Ti, Si0new) + endif ! config_column_physics_type else qi0new = -seaiceDensityIce*seaiceLatentHeatMelting diff --git a/components/mpas-seaice/src/shared/mpas_seaice_forcing.F b/components/mpas-seaice/src/shared/mpas_seaice_forcing.F index 823f549f08a3..4de0435ce92a 100644 --- a/components/mpas-seaice/src/shared/mpas_seaice_forcing.F +++ b/components/mpas-seaice/src/shared/mpas_seaice_forcing.F @@ -2495,8 +2495,8 @@ end subroutine oceanic_forcing subroutine prepare_oceanic_coupling_variables_ncar(block, firstTimeStep) - use seaice_icepack, only: & - seaice_icepack_sea_freezing_temperature + use icepack_intfc, only: & + icepack_sea_freezing_temperature use seaice_column, only: & seaice_column_sea_freezing_temperature @@ -2552,7 +2552,7 @@ subroutine prepare_oceanic_coupling_variables_ncar(block, firstTimeStep) oceanMixedLayerDepth(iCell) = max(oceanMixedLayerDepth(iCell), 0.0_RKIND) ! sea freezing temperature - seaFreezingTemperature(iCell) = seaice_icepack_sea_freezing_temperature(seaSurfaceSalinity(iCell)) + seaFreezingTemperature(iCell) = icepack_sea_freezing_temperature(seaSurfaceSalinity(iCell)) enddo ! iCell else if (trim(config_column_physics_type) == "column_package") then diff --git a/components/mpas-seaice/src/shared/mpas_seaice_icepack.F b/components/mpas-seaice/src/shared/mpas_seaice_icepack.F index 10e971670e18..4e2929e9ef50 100644 --- a/components/mpas-seaice/src/shared/mpas_seaice_icepack.F +++ b/components/mpas-seaice/src/shared/mpas_seaice_icepack.F @@ -48,11 +48,7 @@ module seaice_icepack seaice_icepack_init_itd, & seaice_icepack_init_ocean_conc, & seaice_icepack_ice_strength, & - seaice_icepack_sea_freezing_temperature, & - seaice_icepack_liquidus_temperature, & - seaice_icepack_enthalpy_snow, & seaice_icepack_enthalpy_ice, & - seaice_icepack_salinity_profile, & seaice_init_icepack_constants ! tracer object @@ -5722,54 +5718,9 @@ end subroutine seaice_icepack_ice_strength !----------------------------------------------------------------------- - function seaice_icepack_sea_freezing_temperature(seaSurfaceSalinity) result(seaFreezingTemperature) - - use icepack_intfc, only: & - icepack_sea_freezing_temperature - - real(kind=RKIND), intent(in) :: seaSurfaceSalinity - real(kind=RKIND) :: seaFreezingTemperature - - seaFreezingTemperature = icepack_sea_freezing_temperature(seaSurfaceSalinity) - - end function seaice_icepack_sea_freezing_temperature - - !----------------------------------------------------------------------- - - function seaice_icepack_liquidus_temperature(salinity) result(liquidusTemperature) - - use icepack_intfc, only: & - icepack_liquidus_temperature - - real(kind=RKIND), intent(in) :: salinity - real(kind=RKIND) :: liquidusTemperature - - liquidusTemperature = icepack_liquidus_temperature(salinity) - - end function seaice_icepack_liquidus_temperature - - !----------------------------------------------------------------------- - - function seaice_icepack_enthalpy_snow(snowTemperature) result(snowEnthalpy) - - use icepack_intfc, only: & - icepack_enthalpy_snow - - real(kind=RKIND), intent(in) :: snowTemperature - real(kind=RKIND) :: snowEnthalpy - - snowEnthalpy = icepack_enthalpy_snow(snowTemperature) - - end function seaice_icepack_enthalpy_snow - - !----------------------------------------------------------------------- - function seaice_icepack_enthalpy_ice(iceTemperature, iceSalinity) result(iceEnthalpy) -!echmod - this function will be needed for the prescribed ice configuration - not tested -!echmod - call enthalpy_mush from icepack_mushy_physics via icepack_intfc (not yet available) - - use ice_mushy_physics, only: enthalpy_mush + use icepack_intfc, only: icepack_enthalpy_mush use seaice_constants, only: & seaiceMeltingTemperatureDepression, & ! melting temperature depression factor (C/ppt) @@ -5793,7 +5744,7 @@ function seaice_icepack_enthalpy_ice(iceTemperature, iceSalinity) result(iceEnth if (trim(config_thermodynamics_type) == "mushy") then - qin = enthalpy_mush(iceTemperature, iceSalinity) + qin = icepack_enthalpy_mush(iceTemperature, iceSalinity) else @@ -5805,30 +5756,6 @@ function seaice_icepack_enthalpy_ice(iceTemperature, iceSalinity) result(iceEnth end function seaice_icepack_enthalpy_ice - !----------------------------------------------------------------------- - - function seaice_icepack_salinity_profile(depth) result(iceSalinity) - -!echmod - this function will be needed for the prescribed ice configuration - not tested - - use seaice_constants, only: & - pii, & !echmod - use SHR_CONST_PI - seaiceMaximumSalinity ! ice maximum salinity (ppt) - - real(kind=RKIND), intent(in) :: & - depth ! depth - - real(kind=RKIND) :: & - iceSalinity ! initial salinity profile - - real (kind=RKIND), parameter :: & - nsal = 0.407_RKIND, & - msal = 0.573_RKIND - - iceSalinity = (seaiceMaximumSalinity/2.0_RKIND)*(1.0_RKIND-cos(pii*depth**(nsal/(msal+depth)))) - - end function seaice_icepack_salinity_profile - !----------------------------------------------------------------------- ! initialize constants !----------------------------------------------------------------------- diff --git a/components/mpas-seaice/src/shared/mpas_seaice_initialize.F b/components/mpas-seaice/src/shared/mpas_seaice_initialize.F index 8dd2e1000c92..80fac35089c4 100644 --- a/components/mpas-seaice/src/shared/mpas_seaice_initialize.F +++ b/components/mpas-seaice/src/shared/mpas_seaice_initialize.F @@ -22,6 +22,9 @@ module seaice_initialize use mpas_log, only: mpas_log_write use mpas_timer, only: mpas_timer_start, mpas_timer_stop + use icepack_intfc, only: & + icepack_warnings_aborted + implicit none private @@ -800,9 +803,10 @@ subroutine init_ice_cice_default(& use seaice_constants, only: & seaiceDegreesToRadians - use seaice_icepack, only: & - seaice_icepack_init_trcr, & - seaice_icepack_enthalpy_snow + use icepack_intfc, only: & + icepack_init_trcr, & + icepack_enthalpy_snow + use seaice_column, only: & seaice_column_init_trcr, & seaice_column_enthalpy_snow @@ -928,7 +932,7 @@ subroutine init_ice_cice_default(& snowVolumeCategory(1,iCategory,iCell) = min(iceAreaCategory(1,iCategory,iCell) * initialCategorySnowThickness, & 0.2_RKIND * iceVolumeCategory(1,iCategory,iCell)) - call seaice_icepack_init_trcr(& + call icepack_init_trcr(& airTemperature(iCell), & seaFreezingTemperature(iCell), & initialSalinityProfile(:,iCell), & @@ -972,7 +976,7 @@ subroutine init_ice_cice_default(& surfaceTemperature(1,:,iCell) = seaFreezingTemperature(iCell) if (trim(config_column_physics_type) == "icepack") then do iSnowLayer = 1, nSnowLayers - snowEnthalpy(iSnowLayer,:,iCell) = seaice_icepack_enthalpy_snow(0.0_RKIND) + snowEnthalpy(iSnowLayer,:,iCell) = icepack_enthalpy_snow(0.0_RKIND) end do else if (trim(config_column_physics_type) == "column_package") then do iSnowLayer = 1, nSnowLayers @@ -1306,9 +1310,10 @@ subroutine init_ice_ridging(& use seaice_constants, only: & seaiceDegreesToRadians - use seaice_icepack, only: & - seaice_icepack_init_trcr, & - seaice_icepack_enthalpy_snow + use icepack_intfc, only: & + icepack_init_trcr, & + icepack_enthalpy_snow + use seaice_column, only: & seaice_column_init_trcr, & seaice_column_enthalpy_snow @@ -1449,7 +1454,7 @@ subroutine init_ice_ridging(& snowVolumeCategory(1,iCategory,iCell) = min(iceAreaCategory(1,iCategory,iCell) * initialCategorySnowThickness, & 0.2_RKIND * iceVolumeCategory(1,iCategory,iCell)) - call seaice_icepack_init_trcr(& + call icepack_init_trcr(& airTemperature(iCell), & seaFreezingTemperature(iCell), & initialSalinityProfile(:,iCell), & @@ -1471,7 +1476,7 @@ subroutine init_ice_ridging(& surfaceTemperature(1,:,iCell) = seaFreezingTemperature(iCell) do iSnowLayer = 1, nSnowLayers - snowEnthalpy(iSnowLayer,:,iCell) = seaice_icepack_enthalpy_snow(0.0_RKIND) + snowEnthalpy(iSnowLayer,:,iCell) = icepack_enthalpy_snow(0.0_RKIND) end do endif @@ -2400,8 +2405,9 @@ end subroutine initial_halo_exchanges subroutine initialize_coupler_fields(domain) - use seaice_icepack, only: & - seaice_icepack_liquidus_temperature, & + use icepack_intfc, only: & + icepack_liquidus_temperature + use seaice_icepack, only: & !echmod - use icepack_intfc directly seaice_icepack_init_ocean_conc, & seaice_icepack_initial_air_drag_coefficient use seaice_column, only: & @@ -2524,7 +2530,7 @@ subroutine initialize_coupler_fields(domain) if (trim(config_column_physics_type) == "icepack") then do iCell = 1, nCells - seaFreezingTemperature(iCell) = seaice_icepack_liquidus_temperature(seaSurfaceSalinity(iCell)) + seaFreezingTemperature(iCell) = icepack_liquidus_temperature(seaSurfaceSalinity(iCell)) enddo ! iCell else if (trim(config_column_physics_type) == "column_package") then do iCell = 1, nCells diff --git a/components/mpas-seaice/src/shared/mpas_seaice_prescribed.F b/components/mpas-seaice/src/shared/mpas_seaice_prescribed.F index 57d1daba2eeb..73580b09949f 100644 --- a/components/mpas-seaice/src/shared/mpas_seaice_prescribed.F +++ b/components/mpas-seaice/src/shared/mpas_seaice_prescribed.F @@ -135,16 +135,16 @@ subroutine seaice_run_prescribed_ice(domain) colpkg_enthalpy_ice, & colpkg_salinity_profile -! use icepack_intfc, only: & -! icepack_enthalpy_snow, & -! colpkg_enthalpy_ice, & !echmod - in icepack_init_trcr -! colpkg_salinity_profile !echmod - in icepack_init_thermo - use seaice_column, only: & ! colpkg seaice_column_reinitialize_fluxes, & seaice_column_aggregate + use icepack_intfc, only: & + icepack_enthalpy_snow, & + icepack_salinity_profile + use seaice_icepack, only: & + seaice_icepack_enthalpy_ice, & seaice_icepack_reinitialize_fluxes, & seaice_icepack_aggregate @@ -302,6 +302,25 @@ subroutine seaice_run_prescribed_ice(domain) temperatureGradient = seaFreezingTemperature(iCell) - surfaceTemperature(1,iCategory,iCell) !echmod - add icepack option for ice and snow quantities + if (trim(config_column_physics_type) == "icepack") then + + ! ice quantities + do iIceLayer = 1, nIceLayers + + depth = (real(iIceLayer,kind=RKIND) - 0.5_RKIND) / real(nIceLayers,kind=RKIND) + iceTemperature = surfaceTemperature(1,iCategory,iCell) + temperatureGradient * depth + iceSalinity(iIceLayer,iCategory,iCell) = icepack_salinity_profile(depth) + iceEnthalpy(iIceLayer,iCategory,iCell) = seaice_icepack_enthalpy_ice(iceTemperature,iceSalinity(iIceLayer,iCategory,iCell)) + + enddo ! iIceLayer + + ! snow quantities + do iSnowLayer = 1, nSnowLayers + snowEnthalpy(iSnowLayer,iCategory,iCell) = icepack_enthalpy_snow(surfaceTemperature(1,iCategory,iCell)) + enddo ! iSnowLayer + + else if (trim(config_column_physics_type) == "column_package") then + ! ice quantities do iIceLayer = 1, nIceLayers @@ -319,6 +338,8 @@ subroutine seaice_run_prescribed_ice(domain) endif + endif ! config_column_physics_type + else surfaceTemperature(1,iCategory,iCell) = seaFreezingTemperature(iCell) From a3b23ec30e592d0d45f2e8cf55f5353ab35388b4 Mon Sep 17 00:00:00 2001 From: Elizabeth Hunke Date: Mon, 25 Sep 2023 17:49:51 -0500 Subject: [PATCH 4/5] call icepack interface routines directly with warnings call; clean up --- .../src/shared/mpas_seaice_icepack.F | 78 ++----------------- .../src/shared/mpas_seaice_initialize.F | 27 +++++-- .../src/shared/mpas_seaice_velocity_solver.F | 10 ++- 3 files changed, 34 insertions(+), 81 deletions(-) diff --git a/components/mpas-seaice/src/shared/mpas_seaice_icepack.F b/components/mpas-seaice/src/shared/mpas_seaice_icepack.F index 4e2929e9ef50..29a159771f03 100644 --- a/components/mpas-seaice/src/shared/mpas_seaice_icepack.F +++ b/components/mpas-seaice/src/shared/mpas_seaice_icepack.F @@ -45,11 +45,10 @@ module seaice_icepack seaice_icepack_coupling_prep, & seaice_icepack_finalize, & seaice_icepack_init_trcr, & - seaice_icepack_init_itd, & seaice_icepack_init_ocean_conc, & - seaice_icepack_ice_strength, & seaice_icepack_enthalpy_ice, & - seaice_init_icepack_constants + seaice_init_icepack_constants, & + seaice_icepack_write_warnings ! tracer object type, private :: ciceTracerObjectType @@ -1213,7 +1212,7 @@ subroutine seaice_icepack_postdynamics_time_integration(domain, clock) call mpas_timer_start("Column snow") if (config_use_column_snow_tracers) & - call icepack_snow(domain) + call column_snow(domain) call mpas_timer_stop("Column snow") !----------------------------------------------------------------- @@ -2650,7 +2649,7 @@ end subroutine column_prep_radiation !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! -! icepack_snow +! column_snow ! !> \brief Enable snow grain aging, effective snow density, wind compaction and redistribution !> @@ -2669,7 +2668,7 @@ end subroutine column_prep_radiation ! !----------------------------------------------------------------------- - subroutine icepack_snow(domain) + subroutine column_snow(domain) use icepack_intfc, only: & icepack_step_snow, & @@ -2870,7 +2869,7 @@ subroutine icepack_snow(domain) block => block % next end do - end subroutine icepack_snow + end subroutine column_snow !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! @@ -5577,29 +5576,6 @@ end subroutine seaice_icepack_init_trcr !----------------------------------------------------------------------- - subroutine seaice_icepack_init_itd(& - nCategories, & - categoryThicknessLimits) - - use icepack_intfc, only: & - icepack_init_itd - - integer, intent(in) :: & - nCategories ! number of thickness categories - - real(kind=RKIND), intent(out) :: & - categoryThicknessLimits(0:nCategories) ! category limits (m) - - call icepack_init_itd(& - nCategories, & - categoryThicknessLimits) - - call seaice_icepack_write_warnings(icepack_warnings_aborted()) - - end subroutine seaice_icepack_init_itd - - !----------------------------------------------------------------------- - subroutine seaice_icepack_init_ocean_conc(& oceanAmmoniumConc, & oceanDMSPConc, & @@ -5676,48 +5652,6 @@ end subroutine seaice_icepack_init_ocean_conc !----------------------------------------------------------------------- - subroutine seaice_icepack_ice_strength(& - nCategories, & - iceAreaCell, & - iceVolumeCell, & - openWaterArea, & - iceAreaCategory, & - iceVolumeCategory, & - icePressure) - - use icepack_intfc, only: & - icepack_ice_strength - - integer, intent(in) :: & - nCategories ! number of thickness categories - - real(kind=RKIND), intent(in) :: & - iceAreaCell, & ! concentration of ice - iceVolumeCell, & ! volume per unit area of ice (m) - openWaterArea ! concentration of open water - - real(kind=RKIND), dimension(:), intent(in) :: & - iceAreaCategory, & ! concentration of ice - iceVolumeCategory ! volume per unit area of ice (m) - - real(kind=RKIND), intent(inout) :: & - icePressure ! ice strength (N/m) - - call icepack_ice_strength(& - nCategories, & - iceAreaCell, & - iceVolumeCell, & - openWaterArea, & - iceAreaCategory, & - iceVolumeCategory, & - icePressure) - - call seaice_icepack_write_warnings(icepack_warnings_aborted()) - - end subroutine seaice_icepack_ice_strength - - !----------------------------------------------------------------------- - function seaice_icepack_enthalpy_ice(iceTemperature, iceSalinity) result(iceEnthalpy) use icepack_intfc, only: icepack_enthalpy_mush diff --git a/components/mpas-seaice/src/shared/mpas_seaice_initialize.F b/components/mpas-seaice/src/shared/mpas_seaice_initialize.F index 80fac35089c4..b05347d3269a 100644 --- a/components/mpas-seaice/src/shared/mpas_seaice_initialize.F +++ b/components/mpas-seaice/src/shared/mpas_seaice_initialize.F @@ -22,9 +22,6 @@ module seaice_initialize use mpas_log, only: mpas_log_write use mpas_timer, only: mpas_timer_start, mpas_timer_stop - use icepack_intfc, only: & - icepack_warnings_aborted - implicit none private @@ -805,7 +802,11 @@ subroutine init_ice_cice_default(& use icepack_intfc, only: & icepack_init_trcr, & - icepack_enthalpy_snow + icepack_enthalpy_snow, & + icepack_warnings_aborted + + use seaice_icepack, only: & + seaice_icepack_write_warnings use seaice_column, only: & seaice_column_init_trcr, & @@ -942,6 +943,7 @@ subroutine init_ice_cice_default(& nSnowLayers, & iceEnthalpy(:,iCategory,iCell), & snowEnthalpy(:,iCategory,iCell)) + call seaice_icepack_write_warnings(icepack_warnings_aborted()) enddo ! iCategory else if (trim(config_column_physics_type) == "column_package") then @@ -1026,8 +1028,13 @@ subroutine initial_category_areas_and_volumes(& use seaice_constants, only: & seaicePuny + use icepack_intfc, only: & + icepack_init_itd, & + icepack_warnings_aborted + use seaice_icepack, only: & - seaice_icepack_init_itd + seaice_icepack_write_warnings + use seaice_column, only: & seaice_column_init_itd @@ -1089,9 +1096,10 @@ subroutine initial_category_areas_and_volumes(& if (.not. config_use_column_package) then if (trim(config_column_physics_type) == "icepack") then - call seaice_icepack_init_itd(& + call icepack_init_itd(& nCategories, & categoryThicknessLimits) + call seaice_icepack_write_warnings(icepack_warnings_aborted()) else if (trim(config_column_physics_type) == "column_package") then call seaice_column_init_itd(& nCategories, & @@ -1312,7 +1320,11 @@ subroutine init_ice_ridging(& use icepack_intfc, only: & icepack_init_trcr, & - icepack_enthalpy_snow + icepack_enthalpy_snow, & + icepack_warnings_aborted + + use seaice_icepack, only: & + seaice_icepack_write_warnings use seaice_column, only: & seaice_column_init_trcr, & @@ -1464,6 +1476,7 @@ subroutine init_ice_ridging(& nSnowLayers, & iceEnthalpy(:,iCategory,iCell), & snowEnthalpy(:,iCategory,iCell)) + call seaice_icepack_write_warnings(icepack_warnings_aborted()) enddo ! iCategory diff --git a/components/mpas-seaice/src/shared/mpas_seaice_velocity_solver.F b/components/mpas-seaice/src/shared/mpas_seaice_velocity_solver.F index f1215793aebc..d2be3dd81109 100644 --- a/components/mpas-seaice/src/shared/mpas_seaice_velocity_solver.F +++ b/components/mpas-seaice/src/shared/mpas_seaice_velocity_solver.F @@ -1307,8 +1307,13 @@ end subroutine new_ice_velocities!}}} subroutine ice_strength(domain) + use icepack_intfc, only: & + icepack_ice_strength, & + icepack_warnings_aborted + use seaice_icepack, only: & - seaice_icepack_ice_strength + seaice_icepack_write_warnings + use seaice_column, only: & seaice_column_ice_strength @@ -1417,7 +1422,7 @@ subroutine ice_strength(domain) if (solveStress(iCell) == 1) then ! this routine doesnt reset icePressure - call seaice_icepack_ice_strength(& + call icepack_ice_strength(& nCategories, & iceAreaCell(iCell), & iceVolumeCell(iCell), & @@ -1425,6 +1430,7 @@ subroutine ice_strength(domain) iceAreaCategory(1,:,iCell), & iceVolumeCategory(1,:,iCell), & icePressure(iCell)) + call seaice_icepack_write_warnings(icepack_warnings_aborted()) endif ! solveStress From c3640952ad44aba93f49120cad710ef5953d434c Mon Sep 17 00:00:00 2001 From: Elizabeth Hunke Date: Tue, 26 Sep 2023 11:49:43 -0500 Subject: [PATCH 5/5] updating Icepack to hash 78286bf --- components/mpas-seaice/src/icepack | 2 +- components/mpas-seaice/src/shared/mpas_seaice_icepack.F | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/components/mpas-seaice/src/icepack b/components/mpas-seaice/src/icepack index cb0d9469f784..78286bfa0b2f 160000 --- a/components/mpas-seaice/src/icepack +++ b/components/mpas-seaice/src/icepack @@ -1 +1 @@ -Subproject commit cb0d9469f784f49fe5acaef8b320bbb4507a1b74 +Subproject commit 78286bfa0b2f8a9dc446f61a76e7a57499c96715 diff --git a/components/mpas-seaice/src/shared/mpas_seaice_icepack.F b/components/mpas-seaice/src/shared/mpas_seaice_icepack.F index 29a159771f03..40cfb6aa9060 100644 --- a/components/mpas-seaice/src/shared/mpas_seaice_icepack.F +++ b/components/mpas-seaice/src/shared/mpas_seaice_icepack.F @@ -218,7 +218,7 @@ subroutine seaice_init_icepack_physics_package_parameters(domain) if (config_use_column_package) then ! set non activated variable pointers to other memory - call init_column_non_activated_pointers(domain) + call init_icepack_non_activated_pointers(domain) ! initialize the column package tracer object call init_column_tracer_object(domain, ciceTracerObject) @@ -13518,7 +13518,7 @@ end function config_cice_int !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! -! init_column_non_activated_pointers +! init_icepack_non_activated_pointers ! !> \brief !> \author Adrian K. Turner, LANL @@ -13528,7 +13528,7 @@ end function config_cice_int ! !----------------------------------------------------------------------- - subroutine init_column_non_activated_pointers(domain) + subroutine init_icepack_non_activated_pointers(domain) type(domain_type) :: domain @@ -13854,7 +13854,7 @@ subroutine init_column_non_activated_pointers(domain) block => block % next end do - end subroutine init_column_non_activated_pointers + end subroutine init_icepack_non_activated_pointers !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !