diff --git a/components/mpas-albany-landice/bld/build-namelist b/components/mpas-albany-landice/bld/build-namelist index 619f567f7d60..068363fb98c6 100755 --- a/components/mpas-albany-landice/bld/build-namelist +++ b/components/mpas-albany-landice/bld/build-namelist @@ -456,6 +456,8 @@ if ($MALI_DYNAMIC eq 'TRUE') { } else { add_default($nl, 'config_tracer_advection', 'val'=>"none"); } +add_default($nl, 'config_horiz_tracer_adv_order'); +add_default($nl, 'config_advection_coef_3rd_order'); add_default($nl, 'config_restore_thickness_after_advection'); add_default($nl, 'config_zero_sfcMassBalApplied_over_bare_land'); @@ -509,6 +511,7 @@ add_default($nl, 'config_damagecalvingParameter'); add_default($nl, 'config_ismip6_retreat_k'); add_default($nl, 'config_calving_error_threshold'); add_default($nl, 'config_distribute_unablatedVolumeDynCell'); +add_default($nl, 'config_update_velocity_before_calving'); ################################## # Namelist group: thermal_solver # @@ -574,6 +577,8 @@ add_default($nl, 'config_dynamic_thickness'); add_default($nl, 'config_dt'); add_default($nl, 'config_time_integration'); +add_default($nl, 'config_rk_order'); +add_default($nl, 'config_rk3_stages'); add_default($nl, 'config_adaptive_timestep'); add_default($nl, 'config_min_adaptive_timestep'); add_default($nl, 'config_max_adaptive_timestep'); @@ -635,6 +640,7 @@ add_default($nl, 'config_print_calving_info'); add_default($nl, 'config_print_thermal_info'); add_default($nl, 'config_always_compute_fem_grid'); add_default($nl, 'config_print_velocity_cleanup_details'); +add_default($nl, 'config_check_tracer_monotonicity'); #################################### # Namelist group: subglacial_hydro # diff --git a/components/mpas-albany-landice/bld/build-namelist-section b/components/mpas-albany-landice/bld/build-namelist-section index 533d13a9b1f5..8f4bc139c98a 100644 --- a/components/mpas-albany-landice/bld/build-namelist-section +++ b/components/mpas-albany-landice/bld/build-namelist-section @@ -22,6 +22,8 @@ add_default($nl, 'config_effective_pressure_max'); add_default($nl, 'config_thickness_advection'); add_default($nl, 'config_tracer_advection'); +add_default($nl, 'config_horiz_tracer_adv_order'); +add_default($nl, 'config_advection_coef_3rd_order'); add_default($nl, 'config_restore_thickness_after_advection'); add_default($nl, 'config_zero_sfcMassBalApplied_over_bare_land'); @@ -71,6 +73,7 @@ add_default($nl, 'config_damagecalvingParameter'); add_default($nl, 'config_ismip6_retreat_k'); add_default($nl, 'config_calving_error_threshold'); add_default($nl, 'config_distribute_unablatedVolumeDynCell'); +add_default($nl, 'config_update_velocity_before_calving'); ################################## # Namelist group: thermal_solver # @@ -136,6 +139,8 @@ add_default($nl, 'config_dynamic_thickness'); add_default($nl, 'config_dt'); add_default($nl, 'config_time_integration'); +add_default($nl, 'config_rk_order'); +add_default($nl, 'config_rk3_stages'); add_default($nl, 'config_adaptive_timestep'); add_default($nl, 'config_min_adaptive_timestep'); add_default($nl, 'config_max_adaptive_timestep'); @@ -197,6 +202,7 @@ add_default($nl, 'config_print_calving_info'); add_default($nl, 'config_print_thermal_info'); add_default($nl, 'config_always_compute_fem_grid'); add_default($nl, 'config_print_velocity_cleanup_details'); +add_default($nl, 'config_check_tracer_monotonicity'); #################################### # Namelist group: subglacial_hydro # diff --git a/components/mpas-albany-landice/bld/namelist_files/namelist_defaults_mali.xml b/components/mpas-albany-landice/bld/namelist_files/namelist_defaults_mali.xml index 0db22a1d0d35..5cebbf9cd136 100644 --- a/components/mpas-albany-landice/bld/namelist_files/namelist_defaults_mali.xml +++ b/components/mpas-albany-landice/bld/namelist_files/namelist_defaults_mali.xml @@ -19,6 +19,8 @@ 'fo' 'none' +3 +0.25 .false. .true. @@ -62,6 +64,7 @@ -170.0 0.1 .false. +.false. 'none' @@ -115,6 +118,8 @@ '0000-00-01_00:00:00' 'forward_euler' +2 +3 .false. 3600.0 3.15e9 @@ -129,8 +134,8 @@ .false. 'rpointer.glc' -'0000-01-01_00:00:00' -'0000-01-01_00:00:00' +'0001-01-01_00:00:00' +'0001-01-01_00:00:00' 'none' 'noleap' 'gregorian' @@ -159,6 +164,7 @@ .false. .true. .false. +.false. .false. diff --git a/components/mpas-albany-landice/bld/namelist_files/namelist_definition_mali.xml b/components/mpas-albany-landice/bld/namelist_files/namelist_definition_mali.xml index 32c4085b868b..c04d52840a6c 100644 --- a/components/mpas-albany-landice/bld/namelist_files/namelist_definition_mali.xml +++ b/components/mpas-albany-landice/bld/namelist_files/namelist_definition_mali.xml @@ -139,7 +139,7 @@ Default: Defined in namelist_defaults.xml category="advection" group="advection"> Selection of the method for advecting thickness ('fo' = first-order upwinding). -Valid values: 'fo', 'none' +Valid values: 'fo', 'fct', 'none' Default: Defined in namelist_defaults.xml @@ -147,7 +147,23 @@ Default: Defined in namelist_defaults.xml category="advection" group="advection"> Selection of the method for advecting tracers. -Valid values: 'fo', 'none' +Valid values: 'fo', 'fct', 'none' +Default: Defined in namelist_defaults.xml + + + +Order of polynomial used for tracer reconstruction at cell edges + +Valid values: 2, 3 and 4 +Default: Defined in namelist_defaults.xml + + + +Reconstruction of 3rd-order reconstruction to blend with 4th-order reconstuction. Equivalent to beta in Skamarock and Gassmann (2011) eq. 7. 0 is fully 4th order, 1 is fully 3rd order. + +Valid values: any real between 0 and 1 Default: Defined in namelist_defaults.xml @@ -469,6 +485,14 @@ Valid values: .true. or .false. Default: Defined in namelist_defaults.xml + +If true, add an additional velocity solve between advection and calving. If false, use velocity field from beginning of time step to calculate calving rate. For certain calving laws, like damage threshold calving, it is not necessary to update the velocity before calving, while for von Mises stress and eigencalving, it is more accurate to have an updated velocity state before solving for calvingThickness. + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + @@ -835,9 +859,25 @@ Default: Defined in namelist_defaults.xml -Time integration method (currently, only forward Euler is supported). +Time integration method. -Valid values: 'forward_euler' +Valid values: 'forward_euler' or 'runge_kutta' +Default: Defined in namelist_defaults.xml + + + +Order of Runge-Kutta time integration to use. A value of 1 would be equivalent to forward euler, but will cause an error to avoid unnecessary redundancy. Values of 2 and 3 indicate strong-stability preserving RK2 and RK3. There is currently no support for classical RK2 or RK4 methods. + +Valid values: 2 or 3 +Default: Defined in namelist_defaults.xml + + + +Number of stages for strong stability preserving RK3 time integration. If set to 3, this involves 3 velocity solves and a maximum CFL fraction of 1. If set to 4, this involves 4 velocity solves, but the maximum CFL fraction is 2. + +Valid values: 3 or 4 Default: Defined in namelist_defaults.xml @@ -1133,6 +1173,14 @@ Valid values: .true. or .false. Default: Defined in namelist_defaults.xml + +Check tracer monotonicity at the end of the monotonic advection routine and write warnings to log file if not monotonic. + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 14b8527e6225..279d7fb3c0a0 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -42,7 +42,7 @@ - + possible_values="'fo', 'fct', 'none'" + /> + + + @@ -477,8 +489,16 @@ possible_values="Any time interval of the format 'YYYY-MM-DD_HH:MM:SS', but limited by CFL condition." /> + + - - @@ -638,6 +658,10 @@ description="After velocity is calculated there are a few checks for appropriate values in certain geometric configurations. Setting this option to .true. will cause detailed information about those adjustments to be printed." possible_values=".true. or .false." /> + @@ -741,6 +765,7 @@ + @@ -850,6 +875,7 @@ + @@ -1207,9 +1233,12 @@ is the value of that variable from the *previous* time level! - + @@ -1270,6 +1299,9 @@ is the value of that variable from the *previous* time level! + @@ -1413,6 +1445,9 @@ is the value of that variable from the *previous* time level! + diff --git a/components/mpas-albany-landice/src/landice.cmake b/components/mpas-albany-landice/src/landice.cmake index 3b550fd34067..bed29521388d 100644 --- a/components/mpas-albany-landice/src/landice.cmake +++ b/components/mpas-albany-landice/src/landice.cmake @@ -43,6 +43,8 @@ list(APPEND RAW_SOURCES core_landice/shared/mpas_li_constants.F core_landice/shared/mpas_li_mask.F core_landice/shared/mpas_li_setup.F + core_landice/shared/mpas_li_mesh.F + core_landice/shared/mpas_li_config.F ) # analysis members @@ -57,9 +59,11 @@ list(APPEND RAW_SOURCES core_landice/mode_forward/mpas_li_core.F core_landice/mode_forward/mpas_li_core_interface.F core_landice/mode_forward/mpas_li_time_integration.F - core_landice/mode_forward/mpas_li_time_integration_fe.F + core_landice/mode_forward/mpas_li_time_integration_fe_rk.F core_landice/mode_forward/mpas_li_diagnostic_vars.F core_landice/mode_forward/mpas_li_advection.F + core_landice/mode_forward/mpas_li_advection_fct_shared.F + core_landice/mode_forward/mpas_li_advection_fct.F core_landice/mode_forward/mpas_li_calving.F core_landice/mode_forward/mpas_li_statistics.F core_landice/mode_forward/mpas_li_velocity.F diff --git a/components/mpas-albany-landice/src/mode_forward/Makefile b/components/mpas-albany-landice/src/mode_forward/Makefile index 71cae3b4d1d1..380be8cb0202 100644 --- a/components/mpas-albany-landice/src/mode_forward/Makefile +++ b/components/mpas-albany-landice/src/mode_forward/Makefile @@ -4,9 +4,11 @@ OBJS = mpas_li_core.o \ mpas_li_core_interface.o \ mpas_li_time_integration.o \ - mpas_li_time_integration_fe.o \ + mpas_li_time_integration_fe_rk.o \ mpas_li_diagnostic_vars.o \ mpas_li_advection.o \ + mpas_li_advection_fct_shared.o \ + mpas_li_advection_fct.o \ mpas_li_calving.o \ mpas_li_statistics.o \ mpas_li_velocity.o \ @@ -33,9 +35,9 @@ mpas_li_core.o: mpas_li_time_integration.o \ mpas_li_statistics.o \ mpas_li_calving.o -mpas_li_time_integration.o: mpas_li_time_integration_fe.o +mpas_li_time_integration.o: mpas_li_time_integration_fe_rk.o -mpas_li_time_integration_fe.o: mpas_li_advection.o \ +mpas_li_time_integration_fe_rk.o: mpas_li_advection.o \ mpas_li_calving.o \ mpas_li_thermal.o \ mpas_li_iceshelf_melt.o \ @@ -45,8 +47,12 @@ mpas_li_time_integration_fe.o: mpas_li_advection.o \ mpas_li_bedtopo.o mpas_li_advection.o: mpas_li_thermal.o \ + mpas_li_advection_fct.o \ + mpas_li_advection_fct_shared.o \ mpas_li_diagnostic_vars.o +mpas_li_advection_fct.o: mpas_li_advection_fct_shared.o + mpas_li_calving.o: mpas_li_thermal.o \ mpas_li_advection.o diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F index 4b84ab7206fd..56b018cdd399 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F @@ -11,8 +11,8 @@ ! li_advection ! !> \brief MPAS land ice advection solver -!> \author William Lipscomb -!> \date December 2015 +!> \author William Lipscomb, Trevor Hillebrand +!> \date December 2015, updated Sept 2023 to include FCT !> \details !> This module contains the routines for advecting ice thickness !> and tracers for land ice. @@ -27,6 +27,8 @@ module li_advection use li_setup use li_mask use li_constants + use li_config + use li_mesh implicit none private @@ -44,7 +46,8 @@ module li_advection !-------------------------------------------------------------------- public :: li_advection_thickness_tracers, & li_layer_normal_velocity, & - li_update_geometry + li_update_geometry, & + li_grounded_to_floating !-------------------------------------------------------------------- ! @@ -60,8 +63,8 @@ module li_advection ! subroutine li_advection_thickness_tracers ! !> \brief Advection for ice thickness and tracers -!> \author William Lipscomb -!> \date November 2015 +!> \author William Lipscomb, Trevor Hillebrand +!> \date November 2015, updated Sept 2023 to include FCT !> \details !> This routine (1) computes new values of ice thickness and tracers under !> horizontal advection, (2) applies surface and basal mass balance @@ -72,6 +75,7 @@ module li_advection !----------------------------------------------------------------------- subroutine li_advection_thickness_tracers(& + domain, & dt, & meshPool, & velocityPool, & @@ -82,7 +86,8 @@ subroutine li_advection_thickness_tracers(& advectTracersIn) use li_setup, only: li_calculate_layerThickness - + use li_tracer_advection_fct, only: li_tracer_advection_fct_tend + use li_tracer_advection_fct_shared !----------------------------------------------------------------- ! ! input variables @@ -104,6 +109,7 @@ subroutine li_advection_thickness_tracers(& ! input/output variables ! !----------------------------------------------------------------- + type (domain_type), intent(inout) :: domain !< Input/Output: domain object type (mpas_pool_type), intent(inout) :: & velocityPool !< Input/output: velocity information @@ -162,7 +168,9 @@ subroutine li_advection_thickness_tracers(& layerNormalVelocity, & ! normal component of velocity on cell edges at layer midpoints layerThickness, & ! thickness of each layer layerThicknessOld, & ! old layer thickness - layerThicknessEdge ! layer thickness on upstream edge of cell + layerThicknessEdge, & ! layer thickness on upstream edge of cell + normalVelocity, & ! horizontal velocity on interfaces + layerThicknessEdgeFlux ! higher order thickness flux on layers and edges real (kind=RKIND), dimension(:,:,:), pointer :: & advectedTracers, & ! values of advected tracers @@ -187,6 +195,9 @@ subroutine li_advection_thickness_tracers(& cellMaskTemporaryField, & ! scratch field containing old values of cellMask thermalCellMaskField + ! Allocatable arrays need for flux-corrected transport advection + real (kind=RKIND), dimension(:,:,:), allocatable :: tend + integer, dimension(:), pointer :: & cellMask, & ! integer bitmask for cells edgeMask, & ! integer bitmask for edges @@ -208,15 +219,22 @@ subroutine li_advection_thickness_tracers(& config_ice_density, & ! rhoi config_ocean_density, & ! rhoo config_sea_level, & ! sea level relative to z = 0 - config_thermal_thickness + config_thermal_thickness, & + config_advection_coef_3rd_order + + integer, pointer :: config_num_halos logical :: advectTracers ! if true, then advect tracers as well as thickness integer :: iCell1, iCell2, theGroundedCell + integer, parameter :: horizAdvOrder = 4 + integer, parameter :: vertAdvOrder = 1 real(kind=RKIND) :: GLfluxSign, thicknessFluxEdge integer :: err_tmp + integer :: nTracers + integer, pointer :: maxEdges2 !WHL - debug integer, dimension(:), pointer :: indexToCellID, indexToEdgeID @@ -277,6 +295,8 @@ subroutine li_advection_thickness_tracers(& ! get arrays from the velocity pool call mpas_pool_get_array(velocityPool, 'layerNormalVelocity', layerNormalVelocity) + call mpas_pool_get_array(velocityPool, 'layerThicknessEdgeFlux', layerThicknessEdgeFlux) + call mpas_pool_get_array(velocityPool, 'normalVelocity', normalVelocity) call mpas_pool_get_array(velocityPool, 'fluxAcrossGroundingLine', fluxAcrossGroundingLine) call mpas_pool_get_array(velocityPool, 'fluxAcrossGroundingLineOnCells', fluxAcrossGroundingLineOnCells) @@ -295,11 +315,12 @@ subroutine li_advection_thickness_tracers(& call mpas_pool_get_config(liConfigs, 'config_ice_density', config_ice_density) call mpas_pool_get_config(liConfigs, 'config_ocean_density', config_ocean_density) call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level) + call mpas_pool_get_config(liConfigs, 'config_num_halos', config_num_halos) - !WHL - debug call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges) call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + call mpas_pool_get_dimension(meshPool, 'maxEdges2', maxEdges2) call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell) call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) call mpas_pool_get_array(meshPool, 'indexToCellID', indexToCellID) @@ -339,6 +360,9 @@ subroutine li_advection_thickness_tracers(& call mpas_allocate_scratch_field(thermalCellMaskField, .true.) thermalCellMask => thermalCellMaskField % array + ! define max number of advection cells as max number of edges*2 + nAdvCellsMax = maxEdges2 + ! given the old thickness, compute the thickness in each layer call li_calculate_layerThickness(meshPool, thickness, layerThickness) @@ -348,14 +372,14 @@ subroutine li_advection_thickness_tracers(& ! save old copycellMask for determining cells changing from grounded to floating and vice versa cellMaskTemporaryField % array(:) = cellMask(:) + layerThicknessEdgeFlux(:,:) = 0.0_RKIND !----------------------------------------------------------------- ! Horizontal transport of thickness and tracers !----------------------------------------------------------------- - select case (trim(config_thickness_advection)) - - case ('fo') + if ( (trim(config_thickness_advection) .eq. 'fo') .or. & + (trim(config_thickness_advection) .eq. 'fct') ) then if (config_print_thickness_advection_info) then call mpas_log_write('Using first-order upwind for thickness advection') @@ -367,7 +391,7 @@ subroutine li_advection_thickness_tracers(& endif - if (advectTracers) then + if ( (advectTracers) .or. (trim(config_thickness_advection) .eq. 'fct') ) then ! Copy the old tracer values into the advectedTracers array. ! This requires setting a tracer pointer to either temperature or enthalpy, @@ -381,7 +405,8 @@ subroutine li_advection_thickness_tracers(& thermalPool, & advectedTracers, & surfaceTracers, & - basalTracers) + basalTracers, & + nTracers) else ! no tracer advection; pass a tracer array full of zeroes @@ -389,7 +414,13 @@ subroutine li_advection_thickness_tracers(& endif - ! Transport thickness and tracers using a first-order upwind scheme + if ((trim(config_thickness_advection) .eq. 'fct') .or. & + (trim(config_tracer_advection) .eq. 'fct' ) ) then + allocate(tend(nTracers,nVertLevels,nCells+1)) + tend(:,:,:) = 0.0_RKIND + endif + + ! Transport thickness and tracers ! Note: For the enthalpy scheme, temperature and waterFrac are the primary prognostic ! variables to be updated, but enthalpy is the advected tracer (for reasons of ! energy conservation). @@ -402,8 +433,10 @@ subroutine li_advection_thickness_tracers(& intArgs=(/iCell, indexToCellID(iCell)/), realArgs=(/thickness(iCell), advectedTracers(1,1,iCell)/)) do iEdgeOnCell = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(iEdgeOnCell,iCell) - call mpas_log_write('iEdge (local=$i, global=$i), layerNormalVelocity=$r:', & - intArgs=(/iEdge, indexToEdgeID(iEdge)/), realArgs=(/layerNormalVelocity(1,iEdge)/)) + call mpas_log_write('iEdge (local=$i, global=$i), & + layerNormalVelocity=$r:', & + intArgs=(/iEdge, indexToEdgeID(iEdge)/), & + realArgs=(/layerNormalVelocity(1,iEdge)/)) enddo endif enddo @@ -414,18 +447,80 @@ subroutine li_advection_thickness_tracers(& advectedTracersOld(:,:,:) = advectedTracers(:,:,:) ! compute new values of layer thickness and tracers - call advect_thickness_tracers_upwind(& - dt, & - meshPool, & - layerNormalVelocity, & - layerThicknessEdge, & - layerThicknessOld, & - advectedTracersOld, & - cellMask, & - edgeMask, & - layerThickness, & - advectedTracers, & - err) + if ((trim(config_thickness_advection) .eq. 'fo') .and. & + ( (trim(config_tracer_advection) .eq. 'fo') .or. & + (trim(config_tracer_advection) .eq. 'none') ) ) then + call advect_thickness_tracers_upwind(& + dt, & + meshPool, & + layerNormalVelocity, & + layerThicknessEdge, & + layerThicknessOld, & + advectedTracersOld, & + cellMask, & + edgeMask, & + layerThickness, & + advectedTracers, & + err) + elseif ((trim(config_thickness_advection) .eq. 'fo') .and. & + trim(config_tracer_advection) .eq. 'fct') then + call advect_thickness_tracers_upwind(& + dt, & + meshPool, & + layerNormalVelocity, & + layerThicknessEdge, & + layerThicknessOld, & + advectedTracersOld, & + cellMask, & + edgeMask, & + layerThickness, & + advectedTracers, & + err) + ! Reset tracers after fo advection for fct advection + advectedTracers(:,:,:) = advectedTracersOld(:,:,:) + + ! Pass FO upwind normalThicknessFlux (layerThicknessEdge * layerNormalVelocity) + ! to fct tracer routine + call li_tracer_advection_fct_tend(& + tend, advectedTracers, layerThicknessOld, & + layerThicknessEdge * layerNormalVelocity, 0 * normalVelocity, dt, & + nTracers, layerThicknessEdgeFlux, computeBudgets=.false.)!{{{ + elseif (trim(config_thickness_advection) .eq. 'fct') then + ! Call fct routine for thickness first, and use layerThicknessEdgeFlux + ! returned by that call as normalThicknessFlux for call to tracer fct + call li_tracer_advection_fct_tend(& + tend(nTracers:,:,:), advectedTracers(nTracers:,:,:), layerThicknessOld * 0.0_RKIND + 1.0_RKIND, & + layerNormalVelocity, 0.0_RKIND * normalVelocity, dt, & + 1, layerThicknessEdgeFlux, computeBudgets=.false.) + ! layerThickness is last tracer. However, for some reason + ! this: layerThickness(:,:) = advectedTracers(nTracers,:,:) does not conserve mass! + ! This does conserve mass: + layerThickness(:,:) = layerThickness(:,:) + tend(nTracers,:,:) * dt + + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'layerThicknessEdgeFlux') + call mpas_timer_stop("halo updates") + + if (trim(config_tracer_advection) .eq. 'fct') then + ! Call fct for tracers, using layerThicknessEdgeFlux + ! from fct thickness advection as normalThicknessFlux + call li_tracer_advection_fct_tend(& + tend(1:nTracers-1,:,:), advectedTracers(1:nTracers-1,:,:), layerThicknessOld, & + layerThicknessEdgeFlux, 0.0_RKIND * normalVelocity, dt, & + nTracers-1, computeBudgets=.false.) + elseif (trim(config_tracer_advection) .eq. 'none') then + ! do nothing + else + err_tmp = 1 + call mpas_log_write(trim(config_tracer_advection) // & + ' tracer advection is not currently supported with fct thickness advection.', MPAS_LOG_ERR) + endif + else + err_tmp = 1 + call mpas_log_write("config_thickness_advection = " // trim(config_thickness_advection) // & + ", config_tracer_advection = " // trim(config_tracer_advection) // & + " is not a supported combination.", MPAS_LOG_ERR) + endif if (config_print_thickness_advection_info) then do iCell = 1, nCells @@ -535,10 +630,6 @@ subroutine li_advection_thickness_tracers(& call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) err = ior(err, err_tmp) - ! Calculate volume converted from grounded to floating - ! This needs to be determined after SMB/BMB are applied because those can change floating/grounded state - call grounded_to_floating(cellMaskTemporaryField % array, cellMask, thickness, groundedToFloatingThickness, nCells) - ! Calculate flux across grounding line ! Do this after new thickness & mask have been calculated, including SMB/BMB fluxAcrossGroundingLine(:) = 0.0_RKIND @@ -554,10 +645,16 @@ subroutine li_advection_thickness_tracers(& GLfluxSign = -1.0_RKIND theGroundedCell = iCell2 endif - do k = 1, nVertLevels - thicknessFluxEdge = layerNormalVelocity(k,iEdge) * dvEdge(iEdge) * layerThicknessEdge(k,iEdge) - fluxAcrossGroundingLine(iEdge) = fluxAcrossGroundingLine(iEdge) + GLfluxSign * thicknessFluxEdge / dvEdge(iEdge) - enddo + if (trim(config_thickness_advection) == 'fct') then + do k = 1, nVertLevels + fluxAcrossGroundingLine(iEdge) = fluxAcrossGroundingLine(iEdge) + GLfluxSign * layerThicknessEdgeFlux(k,iEdge) + enddo + else + do k = 1, nVertLevels + layerThicknessEdgeFlux(k,iEdge) = layerNormalVelocity(k,iEdge) * layerThicknessEdge(k,iEdge) + fluxAcrossGroundingLine(iEdge) = fluxAcrossGroundingLine(iEdge) + GLfluxSign * layerThicknessEdgeFlux(k,iEdge) + enddo + endif ! assign to grounded cell in fluxAcrossGroundingLineOnCells if (thickness(theGroundedCell) <= 0.0_RKIND) then ! This should never be the case, but checking to avoid possible divide by zero @@ -601,23 +698,34 @@ subroutine li_advection_thickness_tracers(& geometryPool, & thermalPool, & advectedTracers) - endif - case ('none') + elseif (trim(config_thickness_advection) .eq. 'none') then ! Do nothing - case default + else call mpas_log_write(trim(config_thickness_advection) // ' is not a valid option for thickness/tracer advection.', & MPAS_LOG_ERR) err_tmp = 1 - end select + end if ! config_thickness_advection err = ior(err,err_tmp) + ! Deallocate arrays for fct + if ( (trim(config_thickness_advection) .eq. 'fct') .or. & + (trim(config_tracer_advection) .eq. 'fct') ) then + deallocate( nAdvCellsForEdge, & + advCellsForEdge, & + advCoefs, & + advCoefs3rd, & + advMaskHighOrder, & + advMask2ndOrder, & + tend) + endif + ! clean up call mpas_deallocate_scratch_field(advectedTracersField, .true.) call mpas_deallocate_scratch_field(advectedTracersOldField, .true.) @@ -966,8 +1074,8 @@ end subroutine apply_mass_balance ! subroutine tracer_setup ! !> \brief Assemble a 3D array for tracer advection -!> \author William Lipscomb -!> \date November 2015 +!> \author William Lipscomb, Trevor Hillebrand +!> \date November 2015, updated Sept 2023 for FCT !> \details !> This routine assembles a 3D array for tracer advection. !> Each tracer in the array is transported conservatively, along @@ -981,10 +1089,12 @@ subroutine tracer_setup(& thermalPool, & advectedTracers, & surfaceTracers, & - basalTracers) + basalTracers, & + nTracers) use li_thermal, only: li_temperature_to_enthalpy_kelvin - + use li_tracer_advection_fct_shared + use li_tracer_advection_fct !----------------------------------------------------------------- ! ! input variables @@ -1015,6 +1125,7 @@ subroutine tracer_setup(& basalTracers ! tracer values for new ice at lower surface ! dimension 1 = maxTracers, 2 = nCells + integer, intent(out) :: nTracers !----------------------------------------------------------------- ! ! local variables @@ -1022,7 +1133,8 @@ subroutine tracer_setup(& !----------------------------------------------------------------- character (len=StrKIND), pointer :: & - config_thermal_solver + config_thermal_solver, & + config_thickness_advection logical, pointer :: & config_calculate_damage @@ -1036,7 +1148,8 @@ subroutine tracer_setup(& real (kind=RKIND), dimension(:,:), pointer :: & temperature, & ! interior ice temperature waterFrac, & ! interior water fraction - enthalpy ! interior ice enthalpy + enthalpy, & ! interior ice enthalpy + layerThickness real (kind=RKIND), dimension(:), pointer :: & surfaceAirTemperature, & ! surface air temperature @@ -1049,12 +1162,17 @@ subroutine tracer_setup(& thickness ! ice thickness real (kind=RKIND), dimension(:), pointer :: & - damage ! damage + damage, & ! damage + passiveTracer2d real (kind=RKIND), pointer :: & config_ice_density ! ice density - integer :: iCell, iTracer, k + integer :: iCell, iTracer, k, err, err1, err2 + + err = 0 + err1 = 0 + err2 = 0 ! get dimensions call mpas_pool_get_dimension(meshPool, 'nCells', nCells) @@ -1065,7 +1183,9 @@ subroutine tracer_setup(& ! get arrays from geometry pool call mpas_pool_get_array(geometryPool, 'thickness', thickness) + call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) call mpas_pool_get_array(geometryPool, 'damage', damage) + call mpas_pool_get_array(geometryPool, 'passiveTracer2d', passiveTracer2d) ! get arrays from thermal pool call mpas_pool_get_array(thermalPool, 'temperature', temperature) @@ -1078,6 +1198,7 @@ subroutine tracer_setup(& call mpas_pool_get_config(liConfigs, 'config_thermal_solver', config_thermal_solver) call mpas_pool_get_config(liConfigs, 'config_calculate_damage', config_calculate_damage) call mpas_pool_get_config(liConfigs, 'config_ice_density', config_ice_density) + call mpas_pool_get_config(liConfigs, 'config_thickness_advection', config_thickness_advection) ! initialize advectedTracers(:,:,:) = 0.0_RKIND @@ -1155,8 +1276,41 @@ subroutine tracer_setup(& endif + ! always advect passiveTracer2d + iTracer = iTracer + 1 + ! copy 2d passiveTracer2d to all vertical levels in the tracer array + do k = 1, nVertLevels + advectedTracers(iTracer,k,:) = passiveTracer2d(:) + enddo + surfaceTracers(iTracer,:) = 0.0_RKIND + basalTracers(iTracer,:) = 0.0_RKIND + + if (trim(config_thickness_advection) == 'fct') then + ! increment the tracer index + iTracer = iTracer + 1 + + ! copy 2d thickness to all vertical levels in the tracer array + advectedTracers(iTracer,:,:) = layerThickness(:,:) + surfaceTracers(iTracer,:) = 0.0_RKIND + basalTracers(iTracer,:) = 0.0_RKIND + + endif + + nTracers = iTracer ! Note: Other tracers (e.g., ice age) can be added here as needed. ! May need to increase maxTracers in the Registry. + if ( (trim(config_tracer_advection) == 'fct') .or. & + (trim(config_thickness_advection) == 'fct') ) then + call li_tracer_advection_fct_shared_init(geometryPool, err1) + call li_tracer_advection_fct_init(err2) + + if (err1 /= 0 .or. err2 /= 0) then + err = 1 + call mpas_log_write( & + 'Error encountered during fct tracer advection init', & + MPAS_LOG_ERR, masterOnly=.true.) + endif + endif end subroutine tracer_setup @@ -1228,7 +1382,8 @@ subroutine tracer_finish(& thickness ! ice thickness real (kind=RKIND), dimension(:), pointer :: & - damage ! damage + damage, & ! damage + passiveTracer2d real (kind=RKIND), dimension(:,:), pointer :: & temperature, & ! interior ice temperature @@ -1249,6 +1404,7 @@ subroutine tracer_finish(& ! get arrays from geometry pool call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'damage', damage) + call mpas_pool_get_array(geometryPool, 'passiveTracer2d', passiveTracer2d) ! get arrays from thermal pool call mpas_pool_get_array(thermalPool, 'temperature', temperature) @@ -1299,6 +1455,12 @@ subroutine tracer_finish(& enddo endif + ! Always advect passiveTracer2d + iTracer = iTracer + 1 + do iCell = 1, nCellsSolve + passiveTracer2d(iCell) = sum(advectedTracers(iTracer,:,iCell) * layerThicknessFractions) ! thickness-weighted average + enddo + ! Note: Other tracers (e.g., ice age) can be added here as needed. ! May need to increase maxTracers in the Registry. @@ -1942,7 +2104,7 @@ subroutine vertical_remap(thickness, cellMask, meshPool, layerThickness, tracers end subroutine vertical_remap - subroutine grounded_to_floating(cellMaskOrig, cellMaskNew, thicknessNew, groundedToFloatingThickness, nCells) + subroutine li_grounded_to_floating(cellMaskOrig, cellMaskNew, thicknessNew, groundedToFloatingThickness, nCells) !----------------------------------------------------------------- ! input variables !----------------------------------------------------------------- @@ -1986,7 +2148,7 @@ subroutine grounded_to_floating(cellMaskOrig, cellMaskNew, thicknessNew, grounde endif enddo - end subroutine grounded_to_floating + end subroutine li_grounded_to_floating !*********************************************************************** diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_advection_fct.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_advection_fct.F new file mode 100644 index 000000000000..52c15eef720d --- /dev/null +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_advection_fct.F @@ -0,0 +1,606 @@ +! and the University Corporation for Atmospheric Research (UCAR). +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! +! 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 +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! li_tracer_advection_fct +! +!> \brief MPAS monotonic tracer advection with FCT +!> \author Mark Petersen, David Lee, Doug Jacobsen, Phil Jones, Trevor Hillebrand +!> \date October 2017, updated May 2019, ported to MALI Sept 2023 +!> \details +!> This module contains routines for monotonic advection of tracers +!> using a Flux Corrected Transport (FCT) algorithm +! +!----------------------------------------------------------------------- + +module li_tracer_advection_fct + + ! module includes +#ifdef _ADV_TIMERS + use mpas_timer +#endif + use mpas_kind_types + + use li_config + use li_mesh + use li_tracer_advection_fct_shared + + implicit none + private + save + + ! module private variables + real (kind=RKIND) :: & + coef3rdOrder !< high-order horizontal coefficient + logical :: & + monotonicityCheck !< flag to check monotonicity + + ! public method interfaces + public :: li_tracer_advection_fct_tend, & + li_tracer_advection_fct_init + +!********************************************************************** + + contains + +!********************************************************************** +! +! routine li_tracer_advection_fct_tend +! +!> \brief MPAS monotonic tracer horizontal advection tendency with FCT +!> \author Mark Petersen, David Lee, Doug Jacobsen, Phil Jones, Trevor Hillebrand +!> \date October 2017, updated May 2019, ported to MALI Sept 2023 +!> \details +!> This routine computes the monotonic tracer horizontal advection +!> tendency using a flux-corrected transport (FCT) algorithm. +! +!----------------------------------------------------------------------- + + subroutine li_tracer_advection_fct_tend( & + tend, tracers, layerThickness, & + normalThicknessFlux, w, dt, & + nTracers, & + layerThicknessEdgeFlux, & + computeBudgets)!{{{ + use li_mesh + !----------------------------------------------------------------- + ! Input/Output parameters + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:,:), intent(inout) :: & + tend !< [inout] Tracer tendency to which advection added + real (kind=RKIND), dimension(:,:), intent(inout), optional :: & + layerThicknessEdgeFlux !< [inout] used to compute higher order normalThicknessFlux + !----------------------------------------------------------------- + ! Input parameters + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:,:), intent(inout) :: & + tracers !< [in] Current tracer values + + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThickness, &!< [in] Thickness + normalThicknessFlux, &!< [in] Thichness weighted velocitiy + w !< [in] Vertical velocity + + real (kind=RKIND), intent(in) :: & + dt !< [in] Timestep + + integer, intent(in) :: nTracers + + logical, intent(in) :: & + computeBudgets !< [in] Flag to compute active tracer budgets + !----------------------------------------------------------------- + ! Local variables + !----------------------------------------------------------------- + + integer :: & + i, iCell, iEdge, &! horz indices + cell1, cell2, &! neighbor cell indices + k,kmin, kmax, &! vert index variants + kmin1, kmax1, &! vert index variants + iTracer ! tracer index + + real (kind=RKIND) :: & + signedFactor, &! temp factor including flux sign + tracerMinNew, &! updated tracer minimum + tracerMaxNew, &! updated tracer maximum + tracerUpwindNew, &! tracer updated with upwind flx + scaleFactor, &! factor for normalizing fluxes + flux, &! flux temporary + tracerWeight, &! tracer weighting temporary + invAreaCell1, &! inverse cell area + coef1, coef3 ! temporary coefficients + + real (kind=RKIND), dimension(:), allocatable :: & + wgtTmp, &! vertical temporaries for + flxTmp, sgnTmp ! high-order flux computation + + real (kind=RKIND), dimension(:,:), allocatable :: & + tracerCur, &! reordered current tracer + tracerMax, &! max tracer in neighbors for limiting + tracerMin, &! min tracer in neighbors for limiting + hNewInv, &! inverse of new layer thickness + hProv, &! provisional layer thickness + hProvInv, &! inverse of provisional layer thickness + flxIn, &! flux coming into each cell + flxOut, &! flux going out of each cell + workTend, &! temp for holding some tendency values + lowOrderFlx, &! low order flux for FCT + highOrderFlx ! high order flux for FCT + + real (kind=RKIND), parameter :: & + eps = 1.e-10_RKIND ! small number to avoid numerical difficulties + + ! end of preamble + !---------------- + ! begin code + +#ifdef _ADV_TIMERS + call mpas_timer_start('allocates') +#endif + + ! allocate temporary arrays + allocate(wgtTmp (nVertLevels), & + flxTmp (nVertLevels), & + sgnTmp (nVertLevels), & + tracerCur (nVertLevels ,nCells+1), & + tracerMin (nVertLevels ,nCells), & + tracerMax (nVertLevels ,nCells), & + hNewInv (nVertLevels ,nCells), & + hProv (nVertLevels ,nCells), & + hProvInv (nVertLevels ,nCells), & + flxIn (nVertLevels ,nCells+1), & + flxOut (nVertLevels ,nCells+1), & + workTend (nVertLevels ,nCells+1), & + lowOrderFlx (nVertLevels+1,max(nCells,nEdges)+1), & + highOrderFlx(nVertLevels+1,max(nCells,nEdges)+1)) + + ! Initialize variables so you don't get floating point exceptions + wgtTmp(:) = 0.0_RKIND + flxTmp(:) = 0.0_RKIND + sgnTmp(:) = 0.0_RKIND + tracerCur(:,:) = 0.0_RKIND + tracerMin(:,:) = 0.0_RKIND + tracerMax(:,:) = 0.0_RKIND + hNewInv(:,:) = 0.0_RKIND + hProv(:,:) = 0.0_RKIND + hProvInv(:,:) = 0.0_RKIND + flxIn(:,:) = 0.0_RKIND + flxOut(:,:) = 0.0_RKIND + workTend(:,:) = 0.0_RKIND + lowOrderFlx(:,:) = 0.0_RKIND + highOrderFlx(:,:) = 0.0_RKIND + +#ifdef _ADV_TIMERS + call mpas_timer_stop('allocates') + call mpas_timer_start('prov thickness') +#endif + + ! Compute some provisional layer thicknesses + ! Note: This assumes we are in the first part of the horizontal/ + ! vertical operator splitting, which is true because currently + ! we dont flip order and horizontal is always first. + ! See notes in commit 2cd4a89d. + + do iCell = 1, nCells + invAreaCell1 = dt/areaCell(iCell) + kmin = 1 !minLevelCell(iCell) + kmax = nVertLevels !maxLevelCell(iCell) + do k = kmin, kmax + hProv(k, iCell) = layerThickness(k, iCell) + end do + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i,iCell) + signedFactor = invAreaCell1*dvEdge(iEdge)* & + edgeSignOnCell(i,iCell) + ! Provisional layer thickness is after horizontal + ! thickness flux only + do k = kmin, kmax + hProv(k,iCell) = hProv(k,iCell) & + + signedFactor*normalThicknessFlux(k,iEdge) + end do + end do + ! New layer thickness is after horizontal and vertical + ! thickness flux + do k = kmin, kmax + if (hProv(k, iCell) > 0.0_RKIND) then + hProvInv(k,iCell) = 1.0_RKIND/ hProv(k,iCell) + hNewInv (k,iCell) = 1.0_RKIND/(hProv(k,iCell) - & + dt*w(k,iCell) + dt*w(k+1, iCell)) + else + hProvInv(k, iCell) = 0.0_RKIND + hNewInv(k, iCell) = 0.0_RKIND + endif + end do + end do + +#ifdef _ADV_TIMERS + call mpas_timer_stop('prov thickness') +#endif + + ! Loop over tracers. One tracer is advected at a time. + do iTracer = 1, nTracers + +#ifdef _ADV_TIMERS + call mpas_timer_start('cell init') +#endif + + ! Extract current tracer and change index order to improve locality + do iCell = 1, nCells+1 + do k=1, nVertLevels + tracerCur(k,iCell) = tracers(iTracer,k,iCell) + end do ! k loop + end do ! iCell loop + + ! Compute the high and low order horizontal fluxes. +#ifdef _ADV_TIMERS + call mpas_timer_stop('cell init') + call mpas_timer_start('tracer bounds') +#endif + + ! set nCells to first halo level + ! nCells = nCellsHalo(1) + + ! Determine bounds on tracer (tracerMin and tracerMax) from + ! surrounding cells for later limiting. + + do iCell = 1, nCellsHalo(1) + kmin = 1 !minLevelCell(iCell) + kmax = nVertLevels !maxLevelCell(iCell) + do k=kmin,kmax + tracerMin(k,iCell) = tracerCur(k,iCell) + tracerMax(k,iCell) = tracerCur(k,iCell) + end do + ! TODO: Determine how/if this translates to MALI + do i = 1, nEdgesOnCell(iCell) + cell2 = cellsOnCell(i,iCell) + !kmin1 = max(kmin, minLevelCell(cell2)) + !kmax1 = min(kmax, maxLevelCell(cell2)) + do k=kmin, kmax !kmin1,kmax1 + tracerMax(k,iCell) = max(tracerMax(k,iCell), & + tracerCur(k,cell2)) + tracerMin(k,iCell) = min(tracerMin(k,iCell), & + tracerCur(k,cell2)) + end do ! k loop + end do ! i loop over nEdgesOnCell + end do ! cell loop + +#ifdef _ADV_TIMERS + call mpas_timer_stop('tracer bounds') + call mpas_timer_start('horiz flux') +#endif + ! Need all the edges around the 1 halo cells and owned cells + ! nEdges = nEdgesHalo(2) + + ! Compute the high order horizontal flux + do iEdge = 1, nEdgesHalo(2) + cell1 = cellsOnEdge(1, iEdge) + cell2 = cellsOnEdge(2, iEdge) + + ! compute some common intermediate factors + do k = 1, nVertLevels + wgtTmp(k) = normalThicknessFlux(k,iEdge)* & + advMaskHighOrder(k,iEdge) + sgnTmp(k) = sign(1.0_RKIND, & + normalThicknessFlux(k,iEdge)) + flxTmp(k) = 0.0_RKIND + end do + + ! Compute 3rd or 4th fluxes where requested. + do i = 1, nAdvCellsForEdge(iEdge) + iCell = advCellsForEdge(i,iEdge) + coef1 = advCoefs (i,iEdge) + coef3 = advCoefs3rd (i,iEdge)*coef3rdOrder + do k = 1, nVertLevels !minLevelCell(iCell), maxLevelCell(iCell) + flxTmp(k) = flxTmp(k) + tracerCur(k,iCell)* & + wgtTmp(k)*(coef1 + coef3*sgnTmp(k)) + end do ! k loop + end do ! i loop over nAdvCellsForEdge + + do k=1,nVertLevels + highOrderFlx(k,iEdge) = flxTmp(k) + end do + + ! Compute 2nd order fluxes where needed. + ! Also compute low order upwind horizontal flux (monotonic) + ! Remove low order flux from the high order flux + ! Store left over high order flux in highOrderFlx array + do k = 1, nVertLevels !minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + tracerWeight = advMask2ndOrder(k,iEdge) & + * (1.0_RKIND - advMaskHighOrder(k, iEdge)) & + * (dvEdge(iEdge) * 0.5_RKIND) & + * normalThicknessFlux(k, iEdge) + + lowOrderFlx(k,iEdge) = dvEdge(iEdge) * & + (max(0.0_RKIND,normalThicknessFlux(k,iEdge))*tracerCur(k,cell1) & + + min(0.0_RKIND,normalThicknessFlux(k,iEdge))*tracerCur(k,cell2)) + + highOrderFlx(k,iEdge) = highOrderFlx(k,iedge) & + + tracerWeight*(tracerCur(k,cell1) & + + tracerCur(k,cell2)) + + ! only remove low order flux where high order flux is valid to + ! avoid introducing nonzero values to highOrderFlx where it is invalid + highOrderFlx(k,iEdge) = highOrderFlx(k,iEdge) & + - lowOrderFlx(k,iEdge) * & + (advMaskHighOrder(k,iEdge) + advMask2ndOrder(k,iEdge) & + * (1.0_RKIND - advMaskHighOrder(k, iEdge)) ) + end do ! k loop + end do ! iEdge loop + +#ifdef _ADV_TIMERS + call mpas_timer_stop('horiz flux') + call mpas_timer_start('scale factor build') +#endif + + ! Initialize flux arrays for all cells + do iCell = 1, nCells+1 + do k=1, nVertLevels + workTend(k, iCell) = 0.0_RKIND + flxIn (k, iCell) = 0.0_RKIND + flxOut (k, iCell) = 0.0_RKIND + end do ! k loop + end do ! iCell loop + + ! Need one halo of cells around owned cells + ! nCells = nCellsHalo(1) + + do iCell = 1, nCellsHalo(1) + invAreaCell1 = 1.0_RKIND / areaCell(iCell) + + ! Finish computing the low order horizontal fluxes + ! Upwind fluxes are accumulated in workTend + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + signedFactor = edgeSignOnCell(i, iCell) * invAreaCell1 + + do k = 1, nVertLevels !minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + + ! Here workTend is the advection tendency due to the + ! upwind (low order) fluxes. + workTend(k,iCell) = workTend(k,iCell) & + + signedFactor*lowOrderFlx(k,iEdge) + + ! Accumulate remaining high order fluxes + flxOut(k,iCell) = flxOut(k,iCell) + min(0.0_RKIND, & + signedFactor*highOrderFlx(k,iEdge)) + flxIn (k,iCell) = flxIn (k,iCell) + max(0.0_RKIND, & + signedFactor*highOrderFlx(k,iEdge)) + + end do + end do + + ! Build the factors for the FCT + ! Computed using the bounds that were computed previously, + ! and the bounds on the newly updated value + ! Factors are placed in the flxIn and flxOut arrays + do k = 1, nVertLevels !minLevelCell(iCell), maxLevelCell(iCell) + ! Here workTend is the upwind tendency + tracerUpwindNew = (tracerCur(k,iCell)*layerThickness(k,iCell) & + + dt*workTend(k,iCell))*hProvInv(k,iCell) + tracerMinNew = tracerUpwindNew & + + dt*flxOut(k,iCell)*hProvInv(k,iCell) + tracerMaxNew = tracerUpwindNew & + + dt*flxIn (k,iCell)*hProvInv(k,iCell) + + scaleFactor = (tracerMax(k,iCell) - tracerUpwindNew)/ & + (tracerMaxNew - tracerUpwindNew + eps) + flxIn (k,iCell) = min(1.0_RKIND, & + max(0.0_RKIND, scaleFactor)) + scaleFactor = (tracerUpwindNew - tracerMin(k,iCell))/ & + (tracerUpwindNew - tracerMinNew + eps) + flxOut(k,iCell) = min(1.0_RKIND, & + max(0.0_RKIND, scaleFactor)) + end do ! k loop + end do ! iCell loop + +#ifdef _ADV_TIMERS + call mpas_timer_stop('scale factor build') + call mpas_timer_start('rescale horiz fluxes') +#endif + ! Need all of the edges around owned cells + ! nEdges = nEdgesHalo(1) + ! rescale the high order horizontal fluxes + do iEdge = 1, nEdgesHalo(1) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + do k = 1, nVertLevels !minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + highOrderFlx(k,iEdge) = max(0.0_RKIND,highOrderFlx(k,iEdge))* & + min(flxOut(k,cell1), flxIn (k,cell2)) & + + min(0.0_RKIND,highOrderFlx(k,iEdge))* & + min(flxIn (k,cell1), flxOut(k,cell2)) + end do ! k loop + end do ! iEdge loop + +#ifdef _ADV_TIMERS + call mpas_timer_stop('rescale horiz fluxes') + call mpas_timer_start('flux accumulate') +#endif + + ! Accumulate the scaled high order vertical tendencies + ! and the upwind tendencies + do iCell = 1, nCellsSolve + invAreaCell1 = 1.0_RKIND / areaCell(iCell) + + ! Accumulate the scaled high order horizontal tendencies + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + signedFactor = invAreaCell1*edgeSignOnCell(i,iCell) + do k = 1, nVertLevels !minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + ! workTend on RHS is upwind tendency + ! workTend on LHS is total horiz advect tendency + workTend(k,iCell) = workTend(k,iCell) & + + signedFactor*highOrderFlx(k,iEdge) + end do + end do + + do k = 1, nVertLevels !minLevelCell(iCell), maxLevelCell(iCell) + ! workTend on RHS is total horiz advection tendency + ! tracerCur on LHS is provisional tracer after + ! horizontal fluxes only. + tracerCur(k,iCell) = (tracerCur(k,iCell)* & + layerThickness(k,iCell) & + + dt*workTend(k,iCell)) & + * hProvInv(k,iCell) + tend(iTracer,k,iCell) = tend(iTracer,k,iCell) & + + workTend(k,iCell) + end do + + end do ! iCell loop + +#ifdef _ADV_TIMERS + call mpas_timer_stop('flux accumulate') + call mpas_timer_start('advect diags horiz') +#endif + ! Compute budget and monotonicity diagnostics if needed + + ! Use layerThicknessEdgeFlux from the call to fct for thickness + ! advection as the higher-order normalThicknessFlux for fct tracer advection. + if (present(layerThicknessEdgeFlux)) then + do iEdge = 1,nEdges + do k = 1, nVertLevels + ! Save u*h*T flux on edge for analysis. This variable will be + ! divided by h at the end of the time step. + ! average normal velocities from layer interfaces to layer midpoints + if (dvEdge(iEdge) > 0.0_RKIND) then + layerThicknessEdgeFlux(k,iEdge) = & + (lowOrderFlx(k,iEdge) + highOrderFlx(k,iEdge))/dvEdge(iEdge) + else + layerThicknessEdgeFlux(k,iEdge) = 0.0_RKIND + endif + enddo + enddo + endif + +! ! TODO: Determine if it is necessary to define activeTracerHorizontalAdvectionTendency +! ! for MALI +! !do iCell = 1, nCellsSolve +! !do k = 1, nVertLevels !minLevelCell(iCell), maxLevelCell(iCell) +! ! activeTracerHorizontalAdvectionTendency(iTracer,k,iCell) = & +! ! workTend(k,iCell) +! !end do +! !end do ! iCell loop + + if (monotonicityCheck) then + ! Check tracer values against local min,max to detect + ! non-monotone values and write warning if found + + do iCell = 1, nCellsSolve + do k = 1, nVertLevels !minLevelCell(iCell), maxLevelCell(iCell) + if(tracerCur(k,iCell) < tracerMin(k, iCell)-eps) then + call mpas_log_write( & + 'Horizontal minimum out of bounds on tracer: $i $r $r ', & + MPAS_LOG_WARN, intArgs=(/iTracer/), & + realArgs=(/ tracerMin(k, iCell), tracerCur(k,iCell) /) ) + end if + + if(tracerCur(k,iCell) > tracerMax(k,iCell)+eps) then + call mpas_log_write( & + 'Horizontal maximum out of bounds on tracer: $i $r $r ', & + MPAS_LOG_WARN, intArgs=(/iTracer/), & + realArgs=(/ tracerMax(k, iCell), tracerCur(k,iCell) /) ) + end if + end do + end do + end if ! monotonicity check +#ifdef _ADV_TIMERS + call mpas_timer_stop('advect diags horiz') +#endif + +!TODO: implement vertical advection? + ! Update tracer array + tracers(iTracer,:,:) = tracerCur(:,:) + end do ! iTracer loop + +#ifdef _ADV_TIMERS + call mpas_timer_start('deallocates') +#endif + + deallocate(wgtTmp, & + flxTmp, & + sgnTmp, & + tracerCur, & + tracerMin, & + tracerMax, & + hNewInv, & + hProv, & + hProvInv, & + flxIn, & + flxOut, & + workTend, & + lowOrderFlx, & + highOrderFlx) + +#ifdef _ADV_TIMERS + call mpas_timer_stop('deallocates') +#endif + + end subroutine li_tracer_advection_fct_tend!}}} + +!********************************************************************** +! +! routine li_tracer_advection_fct_init +! +!> \brief MPAS initialize monotonic tracer advection tendency with FCT +!> \author Mark Petersen, David Lee, Doug Jacobsen, Phil Jones +!> \date October 2017, updated May 2019 +!> \details +!> This routine initializes monotonic tracer advection quantities for +!> the flux-corrected transport (FCT) algorithm. +! +!----------------------------------------------------------------------- + + subroutine li_tracer_advection_fct_init(err)!{{{ + + !*** output parameters + + integer, intent(out) :: & + err !< [out] error flag + + ! end of preamble + !---------------- + ! begin code + + err = 0 ! initialize error code to success + + ! Check that the halo is wide enough for FCT + if (config_num_halos < 3) then + call mpas_log_write( & + 'Monotonic advection cannot be used with less than 3 halos.', & + MPAS_LOG_CRIT) + err = -1 + end if + + ! Set blending coefficient if 3rd order horizontal advection chosen + select case (config_horiz_tracer_adv_order) + case (2) + coef3rdOrder = 0.0_RKIND + case (3) + coef3rdOrder = config_advection_coef_3rd_order + case (4) + coef3rdOrder = 0.0_RKIND + case default + coef3rdOrder = 0.0_RKIND + call mpas_log_write( & + 'Invalid value for horizontal advection order, defaulting to 2',& + MPAS_LOG_WARN) + end select ! horizontal advection order + + ! Set flag for checking monotonicity + monotonicityCheck = config_check_tracer_monotonicity + + end subroutine li_tracer_advection_fct_init!}}} + +!*********************************************************************** + +end module li_tracer_advection_fct + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_advection_fct_shared.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_advection_fct_shared.F new file mode 100644 index 000000000000..31c6c10472ce --- /dev/null +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_advection_fct_shared.F @@ -0,0 +1,781 @@ +! 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 +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! li_tracer_advection_fct_shared +! +!> \brief MPAS advection common coefficients and variables +!> \author Doug Jacobsen, Trevor Hillebrand +!> \date 03/09/12, ported to MALI Sept 2023 +!> \details +!> This module contains the routines and arrays for some common +!> tracer advection quantities. +! +!----------------------------------------------------------------------- + +module li_tracer_advection_fct_shared + + use mpas_kind_types + use mpas_derived_types + use mpas_hash + use mpas_sort + use mpas_geometry_utils + use mpas_log + + use li_config + use li_mesh + use li_mask + + implicit none + private + save + + !-------------------------------------------------------------------- + ! Public variables (should only be used by other advection modules) + !-------------------------------------------------------------------- + + integer, public :: & + nAdvCellsMax ! largest number of advection cells for any edge + + integer, public, dimension(:), allocatable :: & + nAdvCellsForEdge ! number of cells contrib to advection at edge + + integer, public, dimension(:,:), allocatable :: & + advCellsForEdge ! index of cells contributing to advection at edge + + real (kind=RKIND), public, dimension(:,:), allocatable :: & + advCoefs, &! common advection coefficients + advCoefs3rd, &! common advection coeffs for high order + advMaskHighOrder, & ! mask where high order advection terms are valid + advMask2ndOrder ! mask where 2nd order advection terms are valid + + !-------------------------------------------------------------------- + ! Public member functions + !-------------------------------------------------------------------- + + public :: li_tracer_advection_fct_shared_init + +!*********************************************************************** + + contains + +!*********************************************************************** +! +! routine li_tracer_advection_fct_shared_init +! +!> \brief MPAS tracer advection coefficients +!> \author Doug Jacobsen, Bill Skamarock, +!> adapted for MALI by Trevor Hillebrand and Matt Hoffman +!> \date 03/09/12, adapted for MALI 09/2023 +!> \details +!> This routine precomputes advection coefficients for horizontal +!> advection of tracers. +! +!----------------------------------------------------------------------- + + subroutine li_tracer_advection_fct_shared_init(geometryPool, err)!{{{ + use li_config + use li_mesh + !----------------------------------------------------------------- + ! Output variables + !----------------------------------------------------------------- + + integer, intent(out) :: err !< [out] Error flag + + !----------------------------------------------------------------- + ! Input variables + !----------------------------------------------------------------- + + type (mpas_pool_type), intent(in) :: geometryPool + + !----------------------------------------------------------------- + ! Local variables + !----------------------------------------------------------------- + + integer :: & + n, i, k, &! loop indices for neighbor, vertical loops + iEdge, iCell, &! loop indices for edge, cell loops + cell1, cell2, & ! neighbor cell indices across edge + iNeighbor, jCell + + integer, dimension(:), allocatable :: & + cellIndx ! cell indices gathered from neighbors + + integer, dimension(:), pointer :: cellMask + + integer, dimension(:,:), allocatable :: & + cellIndxSorted ! nbr cell indices sorted + + real (kind=RKIND), dimension(:,:,:), allocatable :: & + derivTwo !< 2nd derivative values for polynomial fit to tracers + + type (hashtable) :: cell_hash + + ! End preamble + !------------- + ! Begin code + + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + + err = 0 ! initialize error code + + ! define max number of advection cells as max number of edges*2 + nAdvCellsMax = maxEdges2 + + ! Allocate common variables + allocate( nAdvCellsForEdge ( nEdges), & + advCellsForEdge (nAdvCellsMax,nEdges), & + advCoefs (nAdvCellsMax,nEdges), & + advCoefs3rd (nAdvCellsMax,nEdges), & + advMaskHighOrder(nVertLevels, nEdges), & + advMask2ndOrder(nVertLevels, nEdges), & + derivTwo (nAdvCellsMax,2,nEdges)) + + ! Initialize all these allocatable arrays + nAdvCellsForEdge(:) = 0 + advCellsForEdge(:,:) = 0 + advCoefs(:,:) = 0.0_RKIND + advCoefs3rd(:,:) = 0.0_RKIND + advMaskHighOrder(:,:) = 0.0_RKIND + advMask2ndOrder(:,:) = 0.0_RKIND + derivTwo(:,:,:) = 0.0_RKIND + + ! Compute derivTwo array + call computeDerivTwo(derivTwo, err) + if (err /= 0) then + call mpas_log_write( & + 'Error computing derivTwo in ocn advect shared init ', & + MPAS_LOG_CRIT) + deallocate (derivTwo) + return + endif + + allocate(cellIndx ( maxEdges2 + 2), & + cellIndxSorted(2, maxEdges2 + 2)) + + cellIndx(:) = 0 + cellIndxSorted(:,:) = 0 + + ! calculate boundaryCell. A boundary cell is one that + ! has at least one empty or non-dynamic neighbor + boundaryCell(:) = 0 + do iCell = 1, nCells + do iNeighbor = 1, nEdgesOnCell(iCell) + jCell = cellsOnCell(iNeighbor, iCell) + if (.not. li_mask_is_dynamic_ice(cellMask(jCell))) then + boundaryCell(iCell) = 1 + exit ! no need to loop over more neigbors + endif + enddo + enddo + + call mpas_log_write('Marked $i boundary cells out of $i cells belonging to this processor', & + intArgs=(/sum(boundaryCell(1:nCellsSolve)), nCellsSolve/)) + + ! Compute masks and coefficients + do iEdge = 1, nEdges + nAdvCellsForEdge(iEdge) = 0 + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + + ! at boundaries, must stay at low order + if (boundaryCell(cell1) == 1 .or. & + boundaryCell(cell2) == 1) then + advMaskHighOrder(:,iEdge) = 0.0_RKIND + else + advMaskHighOrder(:,iEdge) = 1.0_RKIND + end if + + ! second order mask must have ice on both sides + if ( li_mask_is_dynamic_ice(cellMask(cell1)) .and. li_mask_is_dynamic_ice(cellMask(cell2)) ) then + advMask2ndOrder(:,iEdge) = 1.0_RKIND + else + advMask2ndOrder(:,iEdge) = 0.0_RKIND + end if + + ! do only if this edge flux is needed to update owned cells + if (cell1 <= nCells .and. cell2 <= nCells) then + ! Insert cellsOnEdge to list of advection cells + call mpas_hash_init(cell_hash) + call mpas_hash_insert(cell_hash, cell1) + call mpas_hash_insert(cell_hash, cell2) + cellIndx(1) = cell1 + cellIndx(2) = cell2 + cellIndxSorted(1,1) = indexToCellID(cell1) + cellIndxSorted(2,1) = cell1 + cellIndxSorted(1,2) = indexToCellID(cell2) + cellIndxSorted(2,2) = cell2 + n = 2 + + ! Build unique list of cells used for advection on edge + ! by expanding to the extended neighbor cells + do i = 1, nEdgesOnCell(cell1) + if (.not. mpas_hash_search(cell_hash, & + cellsOnCell(i,cell1))) then + n = n + 1 + cellIndx(n) = cellsOnCell(i, cell1) + cellIndxSorted(1,n) = indexToCellID(cellsOnCell(i,cell1)) + cellIndxSorted(2,n) = cellsOnCell(i,cell1) + call mpas_hash_insert(cell_hash,cellsOnCell(i,cell1)) + end if + end do + + do i = 1, nEdgesOnCell(cell2) + if(.not. mpas_hash_search(cell_hash, & + cellsOnCell(i, cell2))) then + n = n + 1 + cellIndx(n) = cellsOnCell(i, cell2) + cellIndxSorted(1,n) = indexToCellID(cellsOnCell(i,cell2)) + cellIndxSorted(2,n) = cellsOnCell(i,cell2) + call mpas_hash_insert(cell_hash, cellsOnCell(i,cell2)) + end if + end do + + call mpas_hash_destroy(cell_hash) + + ! sort the cell indices by cellID + call mpas_quicksort(n, cellIndxSorted) + + ! store local cell indices for high-order calculations + nAdvCellsForEdge(iEdge) = n + do iCell = 1, nAdvCellsForEdge(iEdge) + advCellsForEdge(iCell,iEdge) = cellIndxSorted(2,iCell) + end do + + ! equation 7 in Skamarock, W. C., & Gassmann, A. (2011): + ! F(u,psi)_{i+1/2} = u_{i+1/2} * + ! [ 1/2 (psi_{i+1} + psi_i) term 1 + ! - 1/12(dx^2psi_{i+1} + dx^2psi_i) term 2 + ! + sign(u) beta/12 (dx^2psi_{i+1} - dx^2psi_i)] term 3 + ! (note minus sign) + ! + ! advCoefs accounts for terms 1 and 2 in SG11 equation 7. + ! Term 1 is the 2nd-order flux-function term. advCoefs + ! accounts for this with the "+ 0.5" lines below. In the + ! advection routines that use these coefficients, the + ! 2nd-order flux loop is then skipped. Term 2 is the + ! 4th-order flux-function term. advCoefs accounts for + ! term 3, the beta term. beta > 0 corresponds to the + ! third-order flux function. The - sign in the derivTwo + ! accumulation is for the i+1 part of term 3, while + ! the + sign is for the i part. + + do i=1,nAdvCellsMax + advCoefs (i,iEdge) = 0.0_RKIND + advCoefs3rd(i,iEdge) = 0.0_RKIND + end do + + ! pull together third and fourth order contributions to the + ! flux first from cell1 + i = mpas_binary_search(cellIndxSorted, 2, 1, & + nAdvCellsForEdge(iEdge), indexToCellID(cell1)) + if (i <= nAdvCellsForEdge(iEdge)) then + advCoefs (i,iEdge) = advCoefs (i,iEdge) & + + derivTwo(1,1,iEdge) + advCoefs3rd(i,iEdge) = advCoefs3rd(i,iEdge) & + + derivTwo(1,1,iEdge) + end if + + do iCell = 1, nEdgesOnCell(cell1) + i = mpas_binary_search(cellIndxSorted, 2, 1, & + nAdvCellsForEdge(iEdge), & + indexToCellID(cellsOnCell(iCell,cell1))) + if (i <= nAdvCellsForEdge(iEdge)) then + advCoefs (i,iEdge) = advCoefs (i,iEdge) & + + derivTwo(iCell+1,1,iEdge) + advCoefs3rd(i,iEdge) = advCoefs3rd(i,iEdge) & + + derivTwo(iCell+1,1,iEdge) + end if + end do + + ! pull together third and fourth order contributions to the + ! flux now from cell2 + i = mpas_binary_search(cellIndxSorted, 2, 1, & + nAdvCellsForEdge(iEdge), indexToCellID(cell2)) + if (i <= nAdvCellsForEdge(iEdge)) then + advCoefs (i,iEdge) = advCoefs (i,iEdge) & + + derivTwo(1,2,iEdge) + advCoefs3rd(i,iEdge) = advCoefs3rd(i,iEdge) & + - derivTwo(1,2,iEdge) + end if + + do iCell = 1, nEdgesOnCell(cell2) + i = mpas_binary_search(cellIndxSorted, 2, 1, & + nAdvCellsForEdge(iEdge), & + indexToCellID(cellsOnCell(iCell,cell2))) + if (i <= nAdvCellsForEdge(iEdge)) then + advCoefs (i,iEdge) = advCoefs (i,iEdge) & + + derivTwo(iCell+1,2,iEdge) + advCoefs3rd(i,iEdge) = advCoefs3rd(i,iEdge) & + - derivTwo(iCell+1,2,iEdge) + end if + end do + + do iCell = 1,nAdvCellsForEdge(iEdge) + advCoefs (iCell,iEdge) = - (dcEdge(iEdge)**2)* & + advCoefs(iCell,iEdge) / 12. + advCoefs3rd(iCell,iEdge) = - (dcEdge(iEdge)**2)* & + advCoefs3rd(iCell,iEdge) / 12. + end do + + ! 2nd order centered contribution + ! place this in the main flux weights + i = mpas_binary_search(cellIndxSorted, 2, 1, & + nAdvCellsForEdge(iEdge), & + indexToCellID(cell1)) + if (i <= nAdvCellsForEdge(iEdge)) then + advCoefs(i,iEdge) = advCoefs(i, iEdge) + 0.5 + end if + + i = mpas_binary_search(cellIndxSorted, 2, 1, & + nAdvCellsForEdge(iEdge), & + indexToCellID(cell2)) + if (i <= nAdvCellsForEdge(iEdge)) then + advCoefs(i,iEdge) = advCoefs(i, iEdge) + 0.5 + end if + + ! multiply by edge length - thus the flux is just dt*ru + ! times the results of the vector-vector multiply + do iCell=1,nAdvCellsForEdge(iEdge) + advCoefs (iCell,iEdge) = dvEdge(iEdge)* & + advCoefs (iCell,iEdge) + advCoefs3rd(iCell,iEdge) = dvEdge(iEdge)* & + advCoefs3rd(iCell,iEdge) + end do + end if ! only do for edges of owned-cells + end do ! end loop over edges + + deallocate(cellIndx, & + cellIndxSorted, & + derivTwo) + + ! If 2nd order advection, disable high-order terms by + ! setting mask to zero. + if (config_horiz_tracer_adv_order == 2) then + advMaskHighOrder(:,:) = 0.0_RKIND + endif + + ! Copy module variables to device + !$acc enter data copyin(nAdvCellsForEdge, advCellsForEdge, & + !$acc advCoefs, advCoefs3rd, advMaskHighOrder) + + !-------------------------------------------------------------------- + + end subroutine li_tracer_advection_fct_shared_init !}}} + +!*********************************************************************** +! +! routine computeDerivTwo +! +!> \brief MPAS deriv two computation +!> \author Doug Jacobsen, Bill Skamarock +!> \date 03/09/12 +!> \details +!> This routine precomputes the second derivative values for tracer +!> advection. It computes cell coefficients for the polynomial fit +!> as described in: +!> Skamarock, W. C., & Gassmann, A. (2011). +!> Conservative Transport Schemes for Spherical Geodesic Meshs: +!> High-Order Flux Operators for ODE-Based Time Integration. +!> Monthly Weather Review, 139(9), 2962-2975. +!> doi:10.1175/MWR-D-10-05056.1 +!> This is performed during model initialization. +! +!----------------------------------------------------------------------- + + subroutine computeDerivTwo(derivTwo, err)!{{{ + + !----------------------------------------------------------------- + ! Output variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:,:), intent(out) :: & + derivTwo !< [out] 2nd deriv values of polynomial for tracer fit + + integer, intent(out) :: err !< [out] returned error flag + + !----------------------------------------------------------------- + ! Local variables + !----------------------------------------------------------------- + + logical :: & + addCell, &! flag for adding cell add to list + doCell ! flag for whether to process cell + + integer, parameter :: & + polynomial_order = 2 ! option for polynomial order, forced to 2 + + integer :: & + i, j, n, &! loop counters + ip1, ip2, &! temps for i+1, i+2 + iCell, iEdge, &! loop indices for cell, edge loops + ma, na, cell_add, mw + + real (kind=RKIND) :: & + xec, yec, zec, &! arc bisection coords + thetae_tmp, &! angle + xv1, xv2, yv1, yv2, zv1, zv2, &! vertex cart coords + pii, &! pi math constant + length_scale, &! length scale + cos2t, costsint, sin2t ! trig function temps + + real (kind=RKIND), dimension(:,:), allocatable :: & + thetae ! sphere angle between vectors + + real (kind=RKIND), dimension(:), allocatable :: & + theta_abs, &! sphere angle + angle_2d ! 2d angle + + integer, dimension(25) :: cell_list + + real (kind=RKIND), dimension(25) :: & + xc, yc, zc, &! cell center coordinates + xp, yp, &! + thetav, thetat, dl_sphere + + real (kind=RKIND) :: & + amatrix(25,25), bmatrix(25,25), wmatrix(25,25) + + logical, parameter :: onSphere=.false. + + ! End preamble + !------------- + ! Begin code + + ! set return code and check proper poly order + err = 0 + + if (polynomial_order > 2) then + call mpas_log_write( & + 'Polynomial for second derivitave can only be 2', & + MPAS_LOG_ERR) + err = 1 + return + end if + + ! allocate local arrays + allocate(thetae(2, nEdges)) + allocate(theta_abs(nCells)) + allocate(angle_2d(maxEdges)) + + ! Initialize derivTwo and pi + + pii = 2.*asin(1.0) + derivTwo(:,:,:) = 0.0_RKIND + + do iCell = 1, nCells + + cell_list(1) = iCell + do i=2, nEdgesOnCell(iCell)+1 + cell_list(i) = cellsOnCell(i-1,iCell) + end do + n = nEdgesOnCell(iCell) + 1 + + !if ( polynomial_order > 2 ) then + ! do i=2, nEdgesOnCell(iCell) + 1 + ! do j=1, nEdgesOnCell( cell_list(i) ) + ! cell_add = cellsOnCell(j,cell_list(i)) + ! addCell = .true. + ! do k=1,n + ! if ( cell_add == cell_list(k) ) addCell = .false. + ! end do + ! if (addCell) then + ! n = n+1 + ! cell_list(n) = cell_add + ! end if + ! end do + ! end do + !end if + + ! check to see if we are reaching outside the halo + + doCell = .true. + do i=1,n + if (cell_list(i) > nCells) doCell = .false. + end do + + if ( .not. doCell ) cycle + + ! compute poynomial fit for this cell if all needed + ! neighbors exist + if ( onSphere ) then + + do i=1,n + j = cell_list(i) + xc(i) = xCell(j) / sphereRadius + yc(i) = yCell(j) / sphereRadius + zc(i) = zCell(j) / sphereRadius + end do + + if ( zc(1) == 1.0_RKIND) then + theta_abs(iCell) = pii/2.0_RKIND + else + theta_abs(iCell) = pii/2.0_RKIND & + - mpas_sphere_angle( xc(1), yc(1), zc(1), & + xc(2), yc(2), zc(2), & + 0.0_RKIND, 0.0_RKIND, 1.0_RKIND) + end if + + ! angles from cell center to neighbor centers (thetav) + + do i=1,n-1 + + ip2 = i+2 + if (ip2 > n) ip2 = 2 + + thetav(i) = mpas_sphere_angle(xc( 1), yc( 1), zc( 1),& + xc(i+1), yc(i+1), zc(i+1),& + xc(ip2), yc(ip2), zc(ip2)) + + dl_sphere(i) = sphereRadius* & + mpas_arc_length(xc( 1),yc( 1),zc( 1), & + xc(i+1),yc(i+1),zc(i+1)) + end do + + length_scale = 1.0_RKIND + do i=1,n-1 + dl_sphere(i) = dl_sphere(i)/length_scale + end do + + ! thetat(1) = 0. ! this defines the x direction, cell center 1 -> + ! this defines the x direction, longitude line + thetat(1) = theta_abs(iCell) + do i=2,n-1 + thetat(i) = thetat(i-1) + thetav(i-1) + end do + + do i=1,n-1 + xp(i) = cos(thetat(i)) * dl_sphere(i) + yp(i) = sin(thetat(i)) * dl_sphere(i) + end do + + else ! On an x-y plane + + do i=1,n-1 + + angle_2d(i) = angleEdge(edgesOnCell(i,iCell)) + iEdge = edgesOnCell(i,iCell) + if ( iCell .ne. cellsOnEdge(1,iEdge)) & + angle_2d(i) = angle_2d(i) - pii + + xp(i) = dcEdge(edgesOnCell(i,iCell)) * cos(angle_2d(i)) + yp(i) = dcEdge(edgesOnCell(i,iCell)) * sin(angle_2d(i)) + + end do + + end if + + ma = n-1 + mw = nEdgesOnCell(iCell) + + bmatrix = 0. + amatrix = 0. + wmatrix = 0. + + if (polynomial_order == 2) then + na = 6 + ma = ma+1 + + amatrix(1,1) = 1. + wmatrix(1,1) = 1. + do i=2,ma + amatrix(i,1) = 1. + amatrix(i,2) = xp(i-1) + amatrix(i,3) = yp(i-1) + amatrix(i,4) = xp(i-1)**2 + amatrix(i,5) = xp(i-1) * yp(i-1) + amatrix(i,6) = yp(i-1)**2 + + wmatrix(i,i) = 1. + end do + + else if (polynomial_order == 3) then + na = 10 + ma = ma+1 + + amatrix(1,1) = 1. + wmatrix(1,1) = 1. + do i=2,ma + amatrix(i,1) = 1. + amatrix(i,2) = xp(i-1) + amatrix(i,3) = yp(i-1) + + amatrix(i,4) = xp(i-1)**2 + amatrix(i,5) = xp(i-1) * yp(i-1) + amatrix(i,6) = yp(i-1)**2 + + amatrix(i,7) = xp(i-1)**3 + amatrix(i,8) = yp(i-1) * (xp(i-1)**2) + amatrix(i,9) = xp(i-1) * (yp(i-1)**2) + amatrix(i,10) = yp(i-1)**3 + + wmatrix(i,i) = 1. + + end do + + else + na = 15 + ma = ma+1 + + amatrix(1,1) = 1. + wmatrix(1,1) = 1. + do i=2,ma + amatrix(i,1) = 1. + amatrix(i,2) = xp(i-1) + amatrix(i,3) = yp(i-1) + + amatrix(i,4) = xp(i-1)**2 + amatrix(i,5) = xp(i-1) * yp(i-1) + amatrix(i,6) = yp(i-1)**2 + + amatrix(i,7) = xp(i-1)**3 + amatrix(i,8) = yp(i-1) * (xp(i-1)**2) + amatrix(i,9) = xp(i-1) * (yp(i-1)**2) + amatrix(i,10) = yp(i-1)**3 + + amatrix(i,11) = xp(i-1)**4 + amatrix(i,12) = yp(i-1) * (xp(i-1)**3) + amatrix(i,13) = (xp(i-1)**2)*(yp(i-1)**2) + amatrix(i,14) = xp(i-1) * (yp(i-1)**3) + amatrix(i,15) = yp(i-1)**4 + + wmatrix(i,i) = 1. + + end do + + do i=1,mw + wmatrix(i,i) = 1. + end do + + end if + + call mpas_poly_fit_2( amatrix, bmatrix, wmatrix, ma, na, 25 ) + + do i=1, nEdgesOnCell(iCell) + ip1 = i+1 + if (ip1 > n-1) ip1 = 1 + + iEdge = edgesOnCell(i,iCell) + + if ( onSphere ) then + xv1 = xVertex(verticesOnEdge(1,iedge)) / sphereRadius + yv1 = yVertex(verticesOnEdge(1,iedge)) / sphereRadius + zv1 = zVertex(verticesOnEdge(1,iedge)) / sphereRadius + xv2 = xVertex(verticesOnEdge(2,iedge)) / sphereRadius + yv2 = yVertex(verticesOnEdge(2,iedge)) / sphereRadius + zv2 = zVertex(verticesOnEdge(2,iedge)) / sphereRadius + else + xv1 = xVertex(verticesOnEdge(1,iedge)) + yv1 = yVertex(verticesOnEdge(1,iedge)) + zv1 = zVertex(verticesOnEdge(1,iedge)) + xv2 = xVertex(verticesOnEdge(2,iedge)) + yv2 = yVertex(verticesOnEdge(2,iedge)) + zv2 = zVertex(verticesOnEdge(2,iedge)) + end if + + if ( onSphere ) then + call mpas_arc_bisect( xv1, yv1, zv1, & + xv2, yv2, zv2, & + xec, yec, zec ) + + thetae_tmp = mpas_sphere_angle( & + xc( 1), yc( 1), zc( 1), & + xc(i+1), yc(i+1), zc(i+1), & + xec, yec, zec ) + thetae_tmp = thetae_tmp + thetat(i) + if (iCell == cellsOnEdge(1,iEdge)) then + thetae(1, edgesOnCell(i,iCell)) = thetae_tmp + else + thetae(2, edgesOnCell(i,iCell)) = thetae_tmp + end if + end if + + end do + + ! fill second derivative stencil for rk advection + + do i=1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i,iCell) + + if ( onSphere ) then + if (iCell == cellsOnEdge(1,iEdge)) then + + cos2t = cos(thetae(1, edgesOnCell(i,iCell))) + sin2t = sin(thetae(1, edgesOnCell(i,iCell))) + costsint = cos2t*sin2t + cos2t = cos2t**2 + sin2t = sin2t**2 + + do j=1,n + derivTwo(j,1,iEdge) = 2.*cos2t *bmatrix(4,j) & + + 2.*costsint*bmatrix(5,j) & + + 2.*sin2t *bmatrix(6,j) + end do + + else + + cos2t = cos(thetae(2, edgesOnCell(i,iCell))) + sin2t = sin(thetae(2, edgesOnCell(i,iCell))) + costsint = cos2t*sin2t + cos2t = cos2t**2 + sin2t = sin2t**2 + + do j=1,n + derivTwo(j,2,iEdge) = 2.*cos2t *bmatrix(4,j) & + + 2.*costsint*bmatrix(5,j) & + + 2.*sin2t *bmatrix(6,j) + end do + end if + + else + + cos2t = cos(angle_2d(i)) + sin2t = sin(angle_2d(i)) + costsint = cos2t*sin2t + cos2t = cos2t**2 + sin2t = sin2t**2 + +! do j=1,n +! +! derivTwo(j,1,iEdge) = 2.*xe(iEdge)*xe(iEdge)*bmatrix(4,j) & +! + 2.*xe(iEdge)*ye(iEdge)*bmatrix(5,j) & +! + 2.*ye(iEdge)*ye(iEdge)*bmatrix(6,j) +! end do + + if (iCell == cellsOnEdge(1,iEdge)) then + do j=1,n + derivTwo(j,1,iEdge) = 2.*cos2t *bmatrix(4,j) & + + 2.*costsint*bmatrix(5,j) & + + 2.*sin2t *bmatrix(6,j) + end do + else + do j=1,n + derivTwo(j,2,iEdge) = 2.*cos2t *bmatrix(4,j) & + + 2.*costsint*bmatrix(5,j) & + + 2.*sin2t *bmatrix(6,j) + end do + end if + + end if + end do + + end do ! end of loop over cells + + deallocate(thetae) + deallocate(theta_abs) + deallocate(angle_2d) + + !-------------------------------------------------------------------- + + end subroutine computeDerivTwo !}}} + +!*********************************************************************** + +end module li_tracer_advection_fct_shared + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F index 042d1b2b5207..304ed0b0349c 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F @@ -74,7 +74,7 @@ module li_calving !> (3) Calve ice based on an ice thickness threshold !----------------------------------------------------------------------- - subroutine li_calve_ice(domain, err) + subroutine li_calve_ice(domain, err, solveVeloAfterCalving) use li_advection @@ -92,7 +92,7 @@ subroutine li_calve_ice(domain, err) ! output variables !----------------------------------------------------------------- integer, intent(out) :: err !< Output: error flag - + logical, intent(out), optional :: solveVeloAfterCalving !----------------------------------------------------------------- ! local variables !----------------------------------------------------------------- @@ -111,13 +111,14 @@ subroutine li_calve_ice(domain, err) logical, pointer :: config_apply_calving_mask real(kind=RKIND), pointer :: config_calving_timescale - integer, pointer :: nCells + integer, pointer :: nCells, nCellsSolve integer, pointer :: config_number_of_blocks real (kind=RKIND), pointer :: deltat !< time step (s) integer, dimension(:), pointer :: & - indexToCellID ! list of global cell IDs + indexToCellID, & ! list of global cell IDs + cellMask real (kind=RKIND) :: & calvingFraction ! fraction of ice that calves in each column; depends on calving_timescale @@ -136,14 +137,22 @@ subroutine li_calve_ice(domain, err) real (kind=RKIND), dimension(:), pointer :: calvingVelocity + real (kind=RKIND), dimension(:), pointer :: areaCell + type (field1dReal), pointer :: originalThicknessField real (kind=RKIND), dimension(:), pointer :: originalThickness + type (field1dInteger), pointer :: cellMaskTemporaryField + integer :: iCell + real (kind=RKIND), parameter :: smallNumber = 1.0e-6_RKIND + integer :: err_tmp + real (kind=RKIND), dimension(1) :: calvingSumLocal, calvingSumGlobal + err = 0 err_tmp = 0 @@ -155,6 +164,15 @@ subroutine li_calve_ice(domain, err) call mpas_pool_get_config(liConfigs, 'config_data_calving', config_data_calving) call mpas_pool_get_config(liConfigs, 'config_number_of_blocks', config_number_of_blocks) + if (present(solveVeloAfterCalving)) then + call mpas_pool_get_subpool(domain % blocklist % structs, 'geometry', geometryPool) + call mpas_pool_get_field(geometryPool, 'cellMaskTemporary', cellMaskTemporaryField) + call mpas_allocate_scratch_field(cellMaskTemporaryField, .true.) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + ! Store cellMask prior to calving + cellMaskTemporaryField % array(:) = cellMask(:) + endif + ! Zero calvingThickness here instead of or in addition to in individual subroutines. ! This is necessary when using damage threshold calving with other calving ! routines. Some individual routines still set calvingThickness to zero, but @@ -206,7 +224,7 @@ subroutine li_calve_ice(domain, err) ! However, the eigencalving method requires multiple applications of the calvingThickness ! to the thickness. So the simplest method to apply data calving is to store the old ! thickness and then set it back when we are done. - if (config_data_calving) then + if ( config_data_calving .or. present(solveVeloAfterCalving) ) then call mpas_pool_get_subpool(domain % blocklist % structs, 'scratch', scratchPool) call mpas_pool_get_field(scratchPool, 'workCell2', originalThicknessField) call mpas_allocate_scratch_field(originalThicknessField, single_block_in = .false.) @@ -337,8 +355,22 @@ subroutine li_calve_ice(domain, err) block => block % next end do + if (present(solveVeloAfterCalving)) then + call mpas_pool_get_array(meshPool, 'areaCell', areaCell) + call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) + call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) + calvingSumLocal(1) = sum(calvingThickness(1:nCellsSolve) * areaCell(1:nCellsSolve) * & + real(li_mask_is_dynamic_ice_int(cellMaskTemporaryField % array(1:nCellsSolve)), RKIND)) + call mpas_dmpar_sum_real_array(domain % dminfo, 1, calvingSumLocal(1), calvingSumGlobal(1)) + if (calvingSumGlobal(1) > smallNumber) then + solveVeloAfterCalving = .true. + else + solveVeloAfterCalving = .false. + endif + call mpas_deallocate_scratch_field(cellMaskTemporaryField, .true.) + endif - if (config_data_calving) then + if ( config_data_calving .or. present(solveVeloAfterCalving) ) then call mpas_deallocate_scratch_field(originalThicknessField, single_block_in=.false.) endif @@ -2322,6 +2354,7 @@ subroutine li_apply_front_ablation_velocity(meshPool, geometryPool, velocityPool integer, dimension(:,:), pointer :: cellsOnEdge integer, dimension(:), pointer :: cellMask, edgeMask integer, dimension(:), pointer :: frontAblationMask + integer, dimension(:), pointer :: openOceanMask real (kind=RKIND), dimension(:), pointer :: calvingThickness, faceMeltingThickness real (kind=RKIND), pointer :: config_sea_level real (kind=RKIND), pointer :: config_calving_error_threshold @@ -2390,6 +2423,7 @@ subroutine li_apply_front_ablation_velocity(meshPool, geometryPool, velocityPool call mpas_pool_get_array(geometryPool, 'ablatedVolumeNonDynCell', ablatedVolumeNonDynCell) call mpas_pool_get_array(geometryPool, 'ablatedVolumeDynCell', ablatedVolumeDynCell) call mpas_pool_get_array(geometryPool, 'frontAblationMask', frontAblationMask) + call mpas_pool_get_array(geometryPool, 'openOceanMask', openOceanMask) call mpas_pool_get_array(geometryPool, 'requiredAblationVolumeNonDynEdge', requiredAblationVolumeNonDynEdge) call mpas_pool_get_array(geometryPool, 'requiredAblationVolumeDynEdge', requiredAblationVolumeDynEdge) call mpas_pool_get_array(geometryPool, 'requiredAblationVolumeNonDynCell', requiredAblationVolumeNonDynCell) @@ -2460,15 +2494,28 @@ subroutine li_apply_front_ablation_velocity(meshPool, geometryPool, velocityPool endif if ( applyToFloating ) then + call find_open_ocean(domain, openOceanMask, err_tmp) + err = ior(err, err_tmp) do iCell = 1,nCells ! If applyToFloating: marine marginal cells are dynamic floating margin cells if ( li_mask_is_floating_ice(cellMask(iCell)) .and. & li_mask_is_dynamic_margin(cellMask(iCell)) ) then - frontAblationMask(iCell) = 1 + do iEdge = 1, nEdgesOnCell(iCell) + iNeighbor = cellsOnCell(iEdge, iCell) + ! only mark this cell if it neighbors the open ocean (don't apply calving to holes in ice shelf) + if (openOceanMask(iNeighbor) == 1) then + frontAblationMask(iCell) = 1 + exit ! no need to keep searching for neighbors + endif + enddo endif enddo endif + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'frontAblationMask') + call mpas_timer_stop("halo updates") + ! Now extend the mask forward to include non-dynamic floating neighbors (appropriate for any of the mask types) ! Use value 2 to indicate the non-dynamic cells do iCell = 1,nCells @@ -3011,6 +3058,93 @@ function scale_edge_length(angleEdgeHere, u, v) end function scale_edge_length +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! routine find_open_ocean +! +!> \brief Generate a mask of open ocean (i.e. ignore 'holes' in ice shelf) +!----------------------------------------------------------------------- + subroutine find_open_ocean(domain, openOceanMask, err) + + !----------------------------------------------------------------- + ! input variables + !----------------------------------------------------------------- + + !----------------------------------------------------------------- + ! input/output variables + !----------------------------------------------------------------- + type (domain_type), intent(inout) :: & + domain !< Input/Output: domain object + + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + integer, dimension(:), intent(out) :: openOceanMask + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + type (block_type), pointer :: block + type (mpas_pool_type), pointer :: geometryPool, meshPool, scratchPool + integer, dimension(:,:), pointer :: cellsOnCell ! list of cells that neighbor each cell + integer, dimension(:), pointer :: nEdgesOnCell + type (field1dInteger), pointer :: seedMaskField, growMaskField + integer, dimension(:), pointer :: seedMask, growMask, cellMask + real (kind=RKIND), dimension(:), pointer :: bedTopography, thickness + real (kind=RKIND), pointer :: config_sea_level + integer, pointer :: nCells + integer :: iCell, iNeighbor, jCell + + err = 0 + + block => domain % blocklist + + call mpas_pool_get_subpool(domain % blocklist % structs, 'geometry', geometryPool) + call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', meshPool) + call mpas_pool_get_subpool(domain % blocklist % structs, 'scratch', scratchPool) + call mpas_pool_get_dimension(meshPool, 'nCells', nCells) + call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) + call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell) + call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) + call mpas_pool_get_array(geometryPool, 'thickness', thickness) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + + call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level) + call mpas_pool_get_field(scratchPool, 'seedMask', seedMaskField) + call mpas_allocate_scratch_field(seedMaskField, single_block_in = .true.) + seedMask => seedMaskField % array + call mpas_pool_get_field(scratchPool, 'growMask', growMaskField) + call mpas_allocate_scratch_field(growMaskField, single_block_in = .true.) + growMask => growMaskField % array + + seedMask(:) = 0 + growMask(:) = 0 + do iCell = 1, nCells + if ((bedTopography(iCell) < config_sea_level) .and. ( .not. li_mask_is_dynamic_ice(cellMask(iCell)) ) ) then + growMask(iCell) = 1 ! includes floating non-dynamic cells + do iNeighbor = 1, nEdgesOnCell(iCell) + jCell = cellsOnCell(iNeighbor, iCell) + if (jCell == nCells + 1) then + seedMask(iCell) = 1 ! this is the seed locations - open ocean along domain boundary + exit + endif + enddo + endif + enddo + + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'seedMask') + call mpas_timer_stop("halo updates") + + call li_flood_fill(seedMask, growMask, domain) + openOceanMask = seedMask + + call mpas_deallocate_scratch_field(seedMaskField, single_block_in=.true.) + call mpas_deallocate_scratch_field(growMaskField, single_block_in=.true.) + end subroutine find_open_ocean + + !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! ! routine damagecalving diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F index a4fb31c371a6..10ee1e94b305 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F @@ -15,6 +15,8 @@ module li_core use li_velocity use li_setup use li_mask + use li_mesh + use li_config implicit none private @@ -122,6 +124,7 @@ function li_core_init(domain, startTimeStamp) result(err) err_tmp = 0 globalErr = 0 + call li_config_init(domain % configs) call mpas_pool_get_config(domain % configs, 'config_do_restart', config_do_restart) ! @@ -207,6 +210,11 @@ function li_core_init(domain, startTimeStamp) result(err) call mpas_timer_stop('reset_io_alarms') call mpas_log_write("Finished processing recurring input streams at initial time.") + ! === + ! Initialize land ice mesh public paramters + call li_meshCreate(domain) + ! === + ! === ! Initialize some time stuff on each block ! === @@ -255,6 +263,11 @@ function li_core_init(domain, startTimeStamp) result(err) block => block % next end do + ! === + ! This is where li_updateMeshFIelds would be called if we wanted to be consistent with mpas-ocean. + ! We don't currently think this needs to be called, so leaving it out for now. + ! === + ! === ! === Initialize modules === ! === @@ -325,8 +338,6 @@ function li_core_init(domain, startTimeStamp) result(err) ! === call li_velocity_external_write_albany_mesh(domain) - ! check for errors and exit - call mpas_dmpar_max_int(domain % dminfo, err, globalErr) ! Find out if any blocks got an error if (globalErr > 0) then call mpas_log_write("An error has occurred in li_core_init. Aborting...", MPAS_LOG_CRIT) @@ -664,6 +675,9 @@ function li_core_finalize(domain) result(err) err = ior(err, err_tmp) + call li_meshDestroy(err_tmp) + err = ior(err, err_tmp) + ! === error check and exit call mpas_dmpar_max_int(domain % dminfo, err, globalErr) ! Find out if any blocks got an error if (globalErr > 0) then diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration.F index c1fe24dfc1d6..d8bd4ae94a67 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration.F @@ -26,7 +26,7 @@ module li_time_integration use mpas_timekeeping use mpas_log - use li_time_integration_fe + use li_time_integration_fe_rk use li_setup use li_constants @@ -166,10 +166,9 @@ subroutine li_timestep(domain, err) !call mpas_log_write('Using ' // trim(config_time_integration) // ' time integration.') select case (config_time_integration) case ('forward_euler') - call li_time_integrator_forwardeuler(domain, err_tmp) - case ('rk4') - call mpas_log_write(trim(config_time_integration) // ' is not currently supported.', MPAS_LOG_ERR) - err_tmp = 1 + call li_time_integrator_forwardeuler_rungekutta(domain, err_tmp) + case ('runge_kutta') + call li_time_integrator_forwardeuler_rungekutta(domain, err_tmp) case default call mpas_log_write(trim(config_time_integration) // ' is not a valid land ice time integration option.', MPAS_LOG_ERR) err_tmp = 1 @@ -177,51 +176,6 @@ subroutine li_timestep(domain, err) err = ior(err,err_tmp) - ! === - ! === Adaptive timestep: update clock information - ! === - ! Set time step in clock object since the time step could have changed - ! Need to get value out of a block, but all blocks should have the same value, so just use first block - if (config_adaptive_timestep) then - block => domain % blocklist - call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) - call mpas_pool_get_array(meshPool, 'deltat', deltat) - ! convert dt in seconds to timeInterval type - call mpas_set_timeInterval(timeStepInterval, dt=deltat, ierr=err_tmp) - err = ior(err,err_tmp) - ! update the clock with the timeInterval - call mpas_set_clock_timestep(domain % clock, timeStepInterval, err_tmp) - err = ior(err,err_tmp) - endif - - ! === - ! === Update clock information - ! === - ! Advance clock - needed to wait until after time step is completed in case the dt has changed! - call mpas_advance_clock(domain % clock) - currTime = mpas_get_clock_time(domain % clock, MPAS_NOW, err_tmp) - call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp, ierr=err_tmp) - err = ior(err, err_tmp) - call mpas_log_write(' Completed timestep. New time is: ' // trim(timeStamp), flushNow=.true.) - - - block => domain % blocklist - do while (associated(block)) - ! Assign the time stamp for this time step - call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) - call mpas_pool_get_array(meshPool, 'xtime', xtime) - xtime = timeStamp - - ! compute time since start of simulation, in days - call mpas_pool_get_array(meshPool, 'simulationStartTime', simulationStartTime) - call mpas_pool_get_array(meshPool, 'daysSinceStart',daysSinceStart) - call mpas_set_time(simulationStartTime_timeType, dateTimeString=simulationStartTime) - call mpas_get_timeInterval(currTime - simulationStartTime_timeType, dt=daysSinceStart) - daysSinceStart = daysSinceStart / seconds_per_day - - block => block % next - end do - ! === error check if (err > 0) then call mpas_log_write("An error has occurred in li_timestep.", MPAS_LOG_ERR) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F similarity index 61% rename from components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F rename to components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F index 2bbc58d3f6f8..6ca13f8deb0d 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F @@ -9,17 +9,17 @@ !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! -! li_time_integration_fe +! li_time_integration_fe_rk ! -!> \brief MPAS land ice Forward Euler time integration scheme -!> \author Matt Hoffman -!> \date 17 April 2011 +!> \brief MPAS land ice Forward Euler and Runge-Kutta time integration schemes +!> \author Matt Hoffman, Trevor Hillebrand +!> \date 17 April 2011, Runge-Kutta added Sept 2023 !> \details -!> This module contains the Forward Euler time integration scheme +!> This module contains the Forward Euler and Runge-Kutta time integration schemes ! !----------------------------------------------------------------------- -module li_time_integration_fe +module li_time_integration_fe_rk use mpas_derived_types use mpas_pool_routines @@ -30,11 +30,13 @@ module li_time_integration_fe use li_advection use li_calving, only: li_calve_ice, li_restore_calving_front, li_calculate_damage, li_finalize_damage_after_advection - use li_thermal, only: li_thermal_solver + use li_thermal, only: li_thermal_solver, li_enthalpy_to_temperature_kelvin use li_iceshelf_melt use li_diagnostic_vars use li_setup use li_constants + use li_mesh + use li_mask implicit none private @@ -51,7 +53,7 @@ module li_time_integration_fe ! !-------------------------------------------------------------------- - public :: li_time_integrator_forwardeuler + public :: li_time_integrator_forwardeuler_rungekutta !-------------------------------------------------------------------- ! @@ -67,20 +69,22 @@ module li_time_integration_fe !*********************************************************************** ! -! routine li_time_integrator_forwardeuler +! routine li_time_integrator_forwardeuler_rungekutta ! -!> \brief Forward Euler time integration scheme -!> \author Matthew Hoffman -!> \date 10 January 2012 +!> \brief Forward Euler and Runge-Kutta time integration schemes +!> \author Matthew Hoffman, Trevor Hillebrand +!> \date 10 January 2012, updated Sept 2023 for Runge-Kutta !> \details -!> This routine performs Forward Euler time integration. +!> This routine performs Forward Euler and Runge-Kutta time integration. ! !----------------------------------------------------------------------- - subroutine li_time_integrator_forwardeuler(domain, err) + subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) use li_subglacial_hydro use li_velocity use li_bedtopo + use li_mask + use li_advection, only: li_grounded_to_floating !----------------------------------------------------------------- ! input variables @@ -92,28 +96,134 @@ subroutine li_time_integrator_forwardeuler(domain, err) type (domain_type), intent(inout) :: & domain !< Input/Output: domain object - !----------------------------------------------------------------- ! output variables !----------------------------------------------------------------- integer, intent(out) :: err !< Output: error flag - !----------------------------------------------------------------- ! local variables !----------------------------------------------------------------- type (block_type), pointer :: block integer :: err_tmp - - logical, pointer :: config_restore_calving_front + type (mpas_pool_type), pointer :: geometryPool, thermalPool, meshPool, velocityPool + + logical, pointer :: config_restore_calving_front, & + config_restore_calving_front_prevent_retreat logical, pointer :: config_calculate_damage logical, pointer :: config_finalize_damage_after_advection + logical, pointer :: config_update_velocity_before_calving + character (len=StrKIND), pointer :: config_thickness_advection + character (len=StrKIND), pointer :: config_thermal_solver + character (len=StrKIND), pointer :: config_time_integration + integer, pointer :: config_rk_order, config_rk3_stages + integer :: rkStage, iCell, iTracer, k + real (kind=RKIND), dimension(:,:), pointer :: layerThickness, & + temperature, & + waterFrac + real (kind=RKIND), dimension(:), pointer :: thickness, damage, passiveTracer2d + real (kind=RKIND), dimension(:), pointer :: sfcMassBalApplied, & + basalMassBalApplied, & + groundedSfcMassBalApplied, & + groundedBasalMassBalApplied, & + floatingBasalMassBalApplied, & + fluxAcrossGroundingLine, & + fluxAcrossGroundingLineOnCells, & + groundedToFloatingThickness + real (kind=RKIND), dimension(:), allocatable :: sfcMassBalAccum, & + basalMassBalAccum, & + groundedSfcMassBalAccum, & + groundedBasalMassBalAccum, & + floatingBasalMassBalAccum, & + fluxAcrossGroundingLineAccum, & + fluxAcrossGroundingLineOnCellsAccum + real (kind=RKIND), dimension(:), pointer :: layerCenterSigma, layerThicknessFractions + integer, dimension(:), pointer :: cellMaskPrev ! cell mask before advection + real (kind=RKIND), dimension(:,:), allocatable :: layerThicknessPrev, & + layerThicknessTmp, & + temperaturePrev, & + waterFracPrev, & + passiveTracer3d, & + passiveTracer3dPrev, & + damage3d, & + damage3dPrev + integer, pointer :: nVertLevels + integer, pointer :: nCells, nEdges + integer :: iCell1, iCell2, iEdge, theGroundedCell + integer, dimension(:), pointer :: edgeMask, cellMask + real (kind=RKIND), pointer :: deltat, config_ice_density + real (kind=RKIND) :: deltatFull + real (kind=RKIND), dimension(4) :: rkSubstepWeights + real (kind=RKIND), dimension(4) :: rkSSPweights + real (kind=RKIND), dimension(4) :: rkTendWeights ! Weights used for calculating budget terms + integer :: nRKstages + logical :: solveVeloAfterCalving err = 0 err_tmp = 0 + solveVeloAfterCalving = .false. call mpas_pool_get_config(liConfigs, 'config_restore_calving_front', config_restore_calving_front) + call mpas_pool_get_config(liConfigs, 'config_restore_calving_front_prevent_retreat', config_restore_calving_front_prevent_retreat) call mpas_pool_get_config(liConfigs, 'config_calculate_damage',config_calculate_damage) - call mpas_pool_get_config(liConfigs, 'config_finalize_damage_after_advection',config_finalize_damage_after_advection) + call mpas_pool_get_config(liConfigs, 'config_finalize_damage_after_advection', config_finalize_damage_after_advection) + call mpas_pool_get_config(liConfigs, 'config_thickness_advection', config_thickness_advection) + call mpas_pool_get_config(liConfigs, 'config_thermal_solver', config_thermal_solver) + call mpas_pool_get_config(liConfigs, 'config_rk_order', config_rk_order) + call mpas_pool_get_config(liConfigs, 'config_rk3_stages', config_rk3_stages) + call mpas_pool_get_config(liConfigs, 'config_time_integration', config_time_integration) + call mpas_pool_get_config(liConfigs, 'config_ice_density', config_ice_density) + call mpas_pool_get_config(liConfigs, 'config_update_velocity_before_calving', config_update_velocity_before_calving) + + call mpas_pool_get_subpool(domain % blocklist % structs, 'geometry', geometryPool) + call mpas_pool_get_subpool(domain % blocklist % structs, 'thermal', thermalPool) + call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', meshPool) + call mpas_pool_get_subpool(domain % blocklist % structs, 'velocity', velocityPool) + + call mpas_pool_get_array(meshPool, 'layerThicknessFractions', layerThicknessFractions) + call mpas_pool_get_array(meshPool, 'deltat', deltat) + call mpas_pool_get_array(meshPool, 'layerCenterSigma', layerCenterSigma) + call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + call mpas_pool_get_dimension(meshPool, 'nCells', nCells) + call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges) + + call mpas_pool_get_array(velocityPool, 'fluxAcrossGroundingLineOnCells', fluxAcrossGroundingLineOnCells) + call mpas_pool_get_array(geometryPool, 'groundedToFloatingThickness', groundedToFloatingThickness) + + allocate(temperaturePrev(nVertLevels, nCells+1)) + allocate(waterFracPrev(nVertLevels, nCells+1)) + allocate(layerThicknessPrev(nVertLevels, nCells+1)) + allocate(layerThicknessTmp(nVertLevels, nCells+1)) + allocate(damage3dPrev(nVertLevels, nCells+1)) + allocate(damage3d(nVertLevels, nCells+1)) + allocate(passiveTracer3dPrev(nVertLevels, nCells+1)) + allocate(passiveTracer3d(nVertLevels, nCells+1)) + allocate(cellMaskPrev(nCells+1)) + + allocate(sfcMassBalAccum(nCells+1)) + allocate(basalMassBalAccum(nCells+1)) + allocate(groundedSfcMassBalAccum(nCells+1)) + allocate(groundedBasalMassBalAccum(nCells+1)) + allocate(floatingBasalMassBalAccum(nCells+1)) + allocate(fluxAcrossGroundingLineAccum(nEdges+1)) + allocate(fluxAcrossGroundingLineOnCellsAccum(nCells+1)) + + temperaturePrev(:,:) = 0.0_RKIND + waterFracPrev(:,:) = 0.0_RKIND + layerThicknessPrev(:,:) = 0.0_RKIND + layerThicknessTmp(:,:) = 0.0_RKIND + damage3dPrev(:,:) = 0.0_RKIND + damage3d(:,:) = 0.0_RKIND + passiveTracer3dPrev(:,:) = 0.0_RKIND + cellMaskPrev(:) = 0 + passiveTracer3d(:,:) = 0.0_RKIND + + sfcMassBalAccum(:) = 0.0_RKIND + basalMassBalAccum(:) = 0.0_RKIND + groundedSfcMassBalAccum(:) = 0.0_RKIND + groundedBasalMassBalAccum(:) = 0.0_RKIND + floatingBasalMassBalAccum(:) = 0.0_RKIND + fluxAcrossGroundingLineAccum(:) = 0.0_RKIND + fluxAcrossGroundingLineOnCellsAccum(:) = 0.0_RKIND ! === Prepare for advection (including CFL checks) =========== ! This has to come first currently, because it sets the time step! @@ -122,6 +232,12 @@ subroutine li_time_integrator_forwardeuler(domain, err) err = ior(err, err_tmp) call mpas_timer_stop("advection prep") +! === Advance the clock before all other physics happen =========== + call mpas_timer_start("advancing clock") + call advance_clock(domain, err_tmp) + err = ior(err, err_tmp) + call mpas_timer_stop("advancing clock") + !TODO: Determine whether grounded melting should in fact be called first ! === Face melting for grounded ice =========== call mpas_timer_start("face melting for grounded ice") @@ -129,6 +245,8 @@ subroutine li_time_integrator_forwardeuler(domain, err) err = ior(err, err_tmp) call mpas_timer_stop("face melting for grounded ice") +! *** TODO: Should basal melt rate calculation and column physics go inside RK loop? *** + ! === Basal melting for floating ice =========== call mpas_timer_start("basal melting for floating ice") call li_basal_melt_floating_ice(domain, err_tmp) @@ -141,27 +259,263 @@ subroutine li_time_integrator_forwardeuler(domain, err) err = ior(err, err_tmp) call mpas_timer_stop("vertical therm") -! === calculate damage =========== - if (config_calculate_damage) then - call mpas_timer_start("damage") - call li_calculate_damage(domain, err_tmp) - err = ior(err, err_tmp) - call mpas_timer_stop("damage") - endif +! *** TODO: Should basal melt rate calculation, column physics, and hydrology go inside RK loop? *** + call mpas_pool_get_array(geometryPool, 'thickness', thickness) + call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) + call mpas_pool_get_array(geometryPool, 'passiveTracer2d', passiveTracer2d) + call mpas_pool_get_array(geometryPool, 'damage', damage) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + + call mpas_pool_get_array(thermalPool, 'temperature', temperature) + call mpas_pool_get_array(thermalPool, 'waterFrac', waterFrac) + ! Save relevant fields before RK loop, to be used in update at the end + temperaturePrev(:,:) = temperature(:,:) + waterFracPrev(:,:) = waterFrac(:,:) + layerThicknessPrev(:,:) = layerThickness(:,:) + cellMaskPrev(:) = cellMask(:) + do k = 1, nVertLevels + damage3dPrev(k,:) = damage(:) + passiveTracer3dPrev(k,:) = passiveTracer2d(:) + enddo + deltatFull = deltat ! Save deltat in order to reset it at end of RK loop + + ! Set RK weights based on desired time integration method. Note + ! that rkSubstepWeights are used to update at each sub-step, and + ! are thus offset from the typical writing of the coefficients + ! by one index. e.g., the coefficients for RK4 are usually written + ! (0, 1/2, 1/2, 1), while we use (1/2, 1/2, 1). The last entry of 1.0 + ! is simply for ease of implementation. + rkSSPweights(:) = -9999.0_RKIND ! Initialized to this value to make it obvious if + ! a weight is used that should not be. Appropriate + ! weights are updated for each case below + rkTendWeights(:) = -9999.0_RKIND + rkSubstepWeights(:) = -9999.0_RKIND + if (trim(config_time_integration) == 'forward_euler') then + nRKstages = 1 + rkSSPweights(1) = 1.0_RKIND + rkTendWeights(1) = 1.0_RKIND + rkSubstepWeights(1) = 1.0_RKIND + elseif ( (trim(config_time_integration) == 'runge_kutta') .and. & + (config_rk_order == 2) ) then + ! use Strong Stability Preserving RK2. Could also + ! add config option to use standard endpoint or midpoint methods, but + ! these are algorithmically more complex and less suitable for our domains + nRKstages = 2 + rkSubstepWeights(1:2) = 1.0_RKIND + + rkSSPweights(1) = 1.0_RKIND + rkSSPweights(2) = 0.5_RKIND + + rkTendWeights(1) = 0.5_RKIND + rkTendWeights(2) = 0.5_RKIND + elseif ( (trim(config_time_integration) == 'runge_kutta') .and. & + (config_rk_order == 3) ) then + if (config_rk3_stages == 3) then + ! use three-stage Strong Stability Preserving RK3 + nRKstages = 3 + rkSubstepWeights(1:3) = 1.0_RKIND + + rkSSPweights(1) = 1.0_RKIND + rkSSPweights(2) = 3.0_RKIND / 4.0_RKIND + rkSSPweights(3) = 1.0_RKIND / 3.0_RKIND + + rkTendWeights(1) = 1.0_RKIND / 6.0_RKIND + rkTendWeights(2) = 1.0_RKIND / 6.0_RKIND + rkTendWeights(3) = 2.0_RKIND / 3.0_RKIND + elseif (config_rk3_stages == 4) then + ! use four-stage Strong Stability Preserving RK3 + nRKstages = 4 + rkSubstepWeights(:) = 0.5_RKIND + + rkSSPweights(1) = 1.0_RKIND + rkSSPweights(2) = 0.0_RKIND + rkSSPweights(3) = 2.0_RKIND / 3.0_RKIND + rkSSPweights(4) = 0.0_RKIND + + rkTendWeights(1) = 1.0_RKIND / 6.0_RKIND + rkTendWeights(2) = 1.0_RKIND / 6.0_RKIND + rkTendWeights(3) = 1.0_RKIND / 6.0_RKIND + rkTendWeights(4) = 1.0_RKIND / 2.0_RKIND + else + err = 1 + call mpas_log_write('config_rk3_stages must 3 or 4', & + messageType=MPAS_LOG_ERR) + return + endif + else + err = 1 + call mpas_log_write('config_time_integration = ' // trim(config_time_integration) & + // ' is not supported with config_rk_order = $i', & + intArgs=(/config_rk_order/), messageType=MPAS_LOG_ERR) + return + endif +! *** Start RK loop *** + do rkStage = 1, nRKstages + call mpas_log_write('beginning rk stage $i of $i', & + intArgs=(/rkStage, nRKstages/)) + deltat = deltatFull * rkSubstepWeights(rkStage) + + ! === calculate damage =========== + if (config_calculate_damage) then + call mpas_timer_start("damage") + call li_calculate_damage(domain, err_tmp) + err = ior(err, err_tmp) + call mpas_timer_stop("damage") + endif + + ! === Compute new state for prognostic variables === + call mpas_timer_start("advect thickness and tracers") + call advection_solver(domain, err_tmp) + err = ior(err, err_tmp) + call mpas_timer_stop("advect thickness and tracers") + + ! If using SSP RK, then update thickness and tracers incrementally. + ! For first RK stage, thickness and tracer updates above are sufficient. + ! Likewise, for the 4-stage SSP RK3 the last stage is just a forward euler update. + if ( (rkStage > 1) .and. (rkStage < 4) ) then + call mpas_pool_get_array(geometryPool, 'thickness', thickness) + call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) + call mpas_pool_get_array(geometryPool, 'passiveTracer2d', passiveTracer2d) + call mpas_pool_get_array(geometryPool, 'damage', damage) + + call mpas_pool_get_array(thermalPool, 'temperature', temperature) + call mpas_pool_get_array(thermalPool, 'waterFrac', waterFrac) + + layerThicknessTmp(:,:) = layerThickness(:,:) + layerThickness(:,:) = rkSSPweights(rkStage) * layerThicknessPrev(:,:) + & + (1.0_RKIND - rkSSPweights(rkStage)) * layerThickness(:,:) + thickness = sum(layerThickness, 1) + ! Calculate masks after updating thickness + call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + err = ior(err, err_tmp) + + if (trim(config_thermal_solver) .ne. 'none') then + do iCell = 1, nCells + do k = 1, nVertLevels + if (layerThickness(k,iCell) > 0.0_RKIND) then + temperature(k,iCell) = ( rkSSPweights(rkStage) * temperaturePrev(k,iCell) * layerThicknessPrev(k,iCell) + & + (1.0_RKIND - rkSSPweights(rkStage)) * temperature(k,iCell) * & + layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) + waterFrac(k,iCell) = ( rkSSPweights(rkStage) * waterFracPrev(k,iCell) * layerThicknessPrev(k,iCell) + & + (1.0_RKIND - rkSSPweights(rkStage)) * waterFrac(k,iCell) * & + layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) + endif + enddo + enddo + endif + + if (config_calculate_damage) then + do iCell = 1, nCells + do k = 1, nVertLevels + if (layerThickness(k,iCell) > 0.0_RKIND) then + damage3d(k,iCell) = ( rkSSPweights(rkStage) * damage3dPrev(k,iCell) * layerThicknessPrev(k,iCell) + & + (1.0_RKIND - rkSSPweights(rkStage)) * damage(iCell) * & + layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) + else + damage3d(k,iCell) = 0.0_RKIND + endif + enddo + damage(iCell) = sum(damage3d(:, iCell) * layerThicknessFractions) + enddo + endif + + do iCell = 1, nCells + do k = 1, nVertLevels + if (layerThickness(k,iCell) > 0.0_RKIND) then + passiveTracer3d(k,iCell) = ( rkSSPweights(rkStage) * passiveTracer3dPrev(k,iCell) * layerThicknessPrev(k,iCell) + & + (1.0_RKIND - rkSSPweights(rkStage)) * passiveTracer2d(iCell) * & + layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) + else + passiveTracer3d(k,iCell) = 0.0_RKIND + endif + enddo + passiveTracer2d(iCell) = sum(passiveTracer3d(:,iCell) * layerThicknessFractions) + enddo + + endif + + ! === Ensure damage is within bounds before velocity solve === + if ( config_finalize_damage_after_advection ) then + call mpas_timer_start("finalize damage") + call li_finalize_damage_after_advection(domain, err_tmp) + err = ior(err, err_tmp) + call mpas_timer_stop("finalize damage") + endif + + call mpas_pool_get_array(geometryPool, 'sfcMassBalApplied', sfcMassBalApplied) + call mpas_pool_get_array(geometryPool, 'groundedSfcMassBalApplied', groundedSfcMassBalApplied) + call mpas_pool_get_array(geometryPool, 'basalMassBalApplied', basalMassBalApplied) + call mpas_pool_get_array(geometryPool, 'floatingBasalMassBalApplied', floatingBasalMassBalApplied) + call mpas_pool_get_array(geometryPool, 'groundedBasalMassBalApplied', groundedBasalMassBalApplied) + call mpas_pooL_get_array(velocityPool, 'fluxAcrossGroundingLine', fluxAcrossGroundingLine) + call mpas_pooL_get_array(velocityPool, 'fluxAcrossGroundingLineOnCells', fluxAcrossGroundingLineOnCells) + call mpas_pool_get_array(geometryPool, 'thickness', thickness) + + ! update budgets + sfcMassBalAccum = sfcMassBalAccum + rkTendWeights(rkStage) * sfcMassBalApplied + groundedSfcMassBalAccum = groundedSfcMassBalAccum + rkTendWeights(rkStage) * groundedSfcMassBalApplied + basalMassBalAccum = basalMassBalAccum + rkTendWeights(rkStage) * basalMassBalApplied + groundedBasalMassBalAccum = groundedBasalMassBalAccum + rkTendWeights(rkStage) * groundedBasalMassBalApplied + floatingBasalMassBalAccum = floatingBasalMassBalAccum + rkTendWeights(rkStage) * floatingBasalMassBalApplied + fluxAcrossGroundingLineAccum = fluxAcrossGroundingLineAccum + rkTendWeights(rkStage) * fluxAcrossGroundingLine + fluxAcrossGroundingLineOnCellsAccum = fluxAcrossGroundingLineOnCellsAccum + rkTendWeights(rkStage) * fluxAcrossGroundingLineOnCells + + ! Halo updates + call mpas_timer_start("halo updates") + + call mpas_dmpar_field_halo_exch(domain, 'thickness') + call mpas_dmpar_field_halo_exch(domain, 'temperature') + call mpas_dmpar_field_halo_exch(domain, 'waterFrac') + call mpas_dmpar_field_halo_exch(domain, 'damage') + call mpas_dmpar_field_halo_exch(domain, 'passiveTracer2d') + + call mpas_timer_stop("halo updates") + + ! Update velocity for each RK step + ! === Solve Velocity ===================== + if ( config_update_velocity_before_calving .or. ( (.not. config_update_velocity_before_calving) & + .and. (rkStage < nRKstages) ) ) then -! === Compute new state for prognostic variables ================================== - call mpas_timer_start("advect thickness and tracers") - call advection_solver(domain, err_tmp) + if (config_restore_calving_front) then + ! restore the calving front to its initial position before velocity solve. + call li_restore_calving_front(domain, err_tmp) + err = ior(err, err_tmp) + endif + + call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) + err = ior(err, err_tmp) + endif +! *** end RK loop *** + enddo + +! Finalize budget updates + ! Update masks after RK integration + call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) err = ior(err, err_tmp) - call mpas_timer_stop("advect thickness and tracers") - -! === finalize damage after advection =========== - if (config_finalize_damage_after_advection) then - call mpas_timer_start("finalize damage") - call li_finalize_damage_after_advection(domain, err_tmp) - err = ior(err, err_tmp) - call mpas_timer_stop("finalize damage") - endif + + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) + + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'cellMask') + call mpas_dmpar_field_halo_exch(domain, 'edgeMask') + call mpas_dmpar_field_halo_exch(domain, 'vertexMask') + call mpas_timer_stop("halo updates") + + ! Calculate volume converted from grounded to floating + ! This needs to be determined after SMB/BMB are applied because those can change floating/grounded state + call li_grounded_to_floating(cellMaskPrev, cellMask, thickness, groundedToFloatingThickness, nCells) + + sfcMassBalApplied(:) = sfcMassBalAccum(:) + groundedSfcMassBalApplied(:) = groundedSfcMassBalAccum(:) + basalMassBalApplied(:) = basalMassBalAccum(:) + groundedBasalMassBalApplied(:) = groundedBasalMassBalAccum(:) + floatingBasalMassBalApplied(:) = floatingBasalMassBalAccum(:) + fluxAcrossGroundingLine(:) = fluxAcrossGroundingLineAccum(:) + fluxAcrossGroundingLineOnCells(:) = fluxAcrossGroundingLineOnCellsAccum(:) + +! Reset time step to full length after RK loop + deltat = deltatFull ! === Update subglacial hydrology =========== ! It's not clear where the best place to call this should be. @@ -177,15 +531,19 @@ subroutine li_time_integrator_forwardeuler(domain, err) call mpas_timer_start("calve_ice") ! ice calving - call li_calve_ice(domain, err_tmp) - err = ior(err, err_tmp) + if ( config_update_velocity_before_calving ) then + call li_calve_ice(domain, err_tmp, solveVeloAfterCalving) + else + call li_calve_ice(domain, err_tmp) + solveVeloAfterCalving = .true. + endif + err = ior(err, err_tmp) if (config_restore_calving_front) then - ! restore the calving front to its initial position; calving options are ignored + ! restore the calving front to its initial position before velocity solve. call li_restore_calving_front(domain, err_tmp) err = ior(err, err_tmp) endif - call mpas_timer_stop("calve_ice") call mpas_timer_start("halo updates") @@ -201,10 +559,11 @@ subroutine li_time_integrator_forwardeuler(domain, err) call li_bedtopo_solve(domain, err=err_tmp) err = ior(err, err_tmp) -! === Solve Velocity ===================== - ! During time-stepping, we always solveVelo - call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) - err = ior(err, err_tmp) +! === Solve velocity for final state ===================== + if (solveVeloAfterCalving) then + call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) + err = ior(err, err_tmp) + endif ! === Calculate diagnostic variables for new state ===================== @@ -216,11 +575,26 @@ subroutine li_time_integrator_forwardeuler(domain, err) ! === error check if (err == 1) then - call mpas_log_write("An error has occurred in li_time_integrator_forwardeuler.", MPAS_LOG_ERR) + call mpas_log_write("An error has occurred in li_time_integrator_forwardeuler_rungekutta.", MPAS_LOG_ERR) endif + deallocate(temperaturePrev) + deallocate(waterFracPrev) + deallocate(layerThicknessPrev) + deallocate(layerThicknessTmp) + deallocate(damage3dPrev) + deallocate(damage3d) + deallocate(passiveTracer3dPrev) + deallocate(passiveTracer3d) + deallocate(cellMaskPrev) + deallocate(sfcMassBalAccum) + deallocate(basalMassBalAccum) + deallocate(groundedSfcMassBalAccum) + deallocate(groundedBasalMassBalAccum) + deallocate(floatingBasalMassBalAccum) + deallocate(fluxAcrossGroundingLineAccum) !-------------------------------------------------------------------- - end subroutine li_time_integrator_forwardeuler + end subroutine li_time_integrator_forwardeuler_rungekutta @@ -280,7 +654,7 @@ subroutine prepare_advection(domain, err) real (kind=RKIND), dimension(:,:), pointer :: normalVelocity real (kind=RKIND), dimension(:,:), pointer :: layerNormalVelocity real (kind=RKIND), pointer :: calvingCFLdt, faceMeltingCFLdt - integer, pointer :: processLimitingTimestep + integer, pointer :: processLimitingTimestep, config_rk_order, config_rk3_stages integer, dimension(:), pointer :: edgeMask logical, pointer :: config_print_thickness_advection_info @@ -288,7 +662,8 @@ subroutine prepare_advection(domain, err) logical, pointer :: config_adaptive_timestep_include_DCFL character (len=StrKIND), pointer :: & - config_thickness_advection ! method for advecting thickness and tracers + config_thickness_advection, & ! method for advecting thickness and tracers + config_time_integration integer :: & allowableAdvecDtProcNumberHere, & @@ -340,6 +715,9 @@ subroutine prepare_advection(domain, err) call mpas_pool_get_config(liConfigs, 'config_adaptive_timestep', config_adaptive_timestep) call mpas_pool_get_config(liConfigs, 'config_adaptive_timestep_include_DCFL', config_adaptive_timestep_include_DCFL) call mpas_pool_get_config(liConfigs, 'config_thickness_advection', config_thickness_advection) + call mpas_pool_get_config(liConfigs, 'config_time_integration', config_time_integration) + call mpas_pool_get_config(liConfigs, 'config_rk_order', config_rk_order) + call mpas_pool_get_config(liConfigs, 'config_rk3_stages', config_rk3_stages) if (trim(config_thickness_advection) == 'none') then if (config_adaptive_timestep) then @@ -385,7 +763,6 @@ subroutine prepare_advection(domain, err) err = ior(err, err_tmp) allowableAdvecDtOnProc = min(allowableAdvecDtOnProc, allowableAdvecDt) - ! Calculate diffusive CFL timestep, if needed ! This used to be only calculated if (config_adaptive_timestep_include_DCFL) but for simplicity, ! now it is always calculated. That allows assessment of the DCFL even when it is not being obeyed @@ -410,6 +787,12 @@ subroutine prepare_advection(domain, err) block => block % next end do + ! If using 4-stage SSPRK3, CFL number of 2 is theoretically allowed + if ( (trim(config_time_integration) == 'runge_kutta') .and. & + (config_rk_order == 3) .and. (config_rk3_stages == 4) ) then + allowableAdvecDtOnProc = allowableAdvecDtOnProc * 2.0_RKIND + allowableDiffDtOnProc = allowableDiffDtOnProc * 2.0_RKIND + endif ! Local advective CFL info call mpas_set_timeInterval(allowableAdvecDtOnProcInterval, dt=allowableAdvecDtOnProc, ierr=err_tmp) @@ -598,7 +981,6 @@ subroutine advection_solver(domain, err) ! input/output variables !----------------------------------------------------------------- type (domain_type), intent(inout) :: domain !< Input/Output: domain object - !----------------------------------------------------------------- ! output variables !----------------------------------------------------------------- @@ -620,6 +1002,7 @@ subroutine advection_solver(domain, err) real (kind=RKIND), dimension(:), pointer :: thicknessOld real (kind=RKIND), dimension(:), pointer :: thicknessNew real (kind=RKIND), dimension(:), pointer :: thickness + real (kind=RKIND), dimension(:,:), pointer :: layerThickness real (kind=RKIND), dimension(:,:), pointer :: temperature real (kind=RKIND), dimension(:,:), pointer :: waterFrac @@ -661,7 +1044,8 @@ subroutine advection_solver(domain, err) ! === ! === Calculate layerThicknessEdge, which is needed for advection ! === - if (trim(config_thickness_advection) == 'fo' .or. trim(config_tracer_advection) == 'fo') then + if (trim(config_thickness_advection) == 'fo' .or. trim(config_tracer_advection) == 'fo' .or. & + trim(config_thickness_advection) == 'fct' .or. trim(config_tracer_advection) == 'fct') then block => domain % blocklist do while (associated(block)) @@ -697,7 +1081,9 @@ subroutine advection_solver(domain, err) call mpas_pool_get_subpool(block % structs, 'scratch', scratchPool) call mpas_pool_get_array(meshPool, 'deltat', deltat) - if (trim(config_thickness_advection) == 'fo' .and. trim(config_tracer_advection) == 'fo') then + if ( (trim(config_thickness_advection) == 'fo' .and. trim(config_tracer_advection) == 'fo') .or. & + (trim(config_thickness_advection) == 'fct') .or. & + (trim(config_thickness_advection) == 'fo' .and. trim(config_tracer_advection) == 'fct') ) then ! Note: This subroutine requires that thickness and tracers are correct in halos @@ -706,6 +1092,7 @@ subroutine advection_solver(domain, err) endif call li_advection_thickness_tracers(& + domain, & deltat, & meshPool, & velocityPool, & @@ -724,6 +1111,7 @@ subroutine advection_solver(domain, err) endif call li_advection_thickness_tracers(& + domain, & deltat, & meshPool, & velocityPool, & @@ -752,7 +1140,6 @@ subroutine advection_solver(domain, err) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_array(geometryPool, 'thickness', thickness, timeLevel=1) - allocate( masktmp(nCells + 1) ) masktmp = 0 @@ -789,6 +1176,7 @@ subroutine advection_solver(domain, err) call mpas_dmpar_field_halo_exch(domain, 'waterFrac') call mpas_dmpar_field_halo_exch(domain, 'enthalpy') call mpas_dmpar_field_halo_exch(domain, 'damage') + call mpas_dmpar_field_halo_exch(domain, 'passiveTracer2d') call mpas_timer_stop("halo updates") @@ -803,7 +1191,10 @@ subroutine advection_solver(domain, err) if ( config_restore_thickness_after_advection ) then call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'thicknessOld', thicknessOld) + call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) thickness(:) = thicknessOld(:) + + call li_calculate_layerThickness(meshPool, thickness, layerThickness) endif call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) @@ -1080,6 +1471,115 @@ subroutine calculate_layerThicknessEdge(meshPool, geometryPool, velocityPool, er !-------------------------------------------------------------------- end subroutine calculate_layerThicknessEdge +!*********************************************************************** +! +! routine advance_clock +! +!> \brief Update the clock +!> \author Holly Han +!> \date July 2023 +!> \details +!> This routine advances the clock and is called right after the +!> subroutine 'prepare_advection' such that the clock is updated +!> earlier in the timestep before other physics are calculated +!----------------------------------------------------------------------- + + subroutine advance_clock(domain, err) + + use mpas_timekeeping + + !----------------------------------------------------------------- + ! input variables + !----------------------------------------------------------------- + + !----------------------------------------------------------------- + ! input/output variables + !----------------------------------------------------------------- + + type (domain_type), intent(inout) :: domain !< Input/Output: domain object + + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + type (block_type), pointer :: block + + type (mpas_pool_type), pointer :: meshPool + type (MPAS_Time_Type) :: currTime !< current time as time type + type (MPAS_TimeInterval_type) :: timeStepInterval !< the current time step as an interval + type (MPAS_Time_type) :: simulationStartTime_timeType + + logical, pointer :: config_adaptive_timestep + + character (len=StrKIND), pointer :: xtime + character (len=StrKIND), pointer :: simulationStartTime + character(len=StrKIND) :: timeStamp !< current time as a string + + real (kind=RKIND), pointer :: daysSinceStart + real (kind=RKIND), pointer :: deltat ! variable in blocks + + integer :: err_tmp + + err = 0 + + call mpas_pool_get_config(liConfigs, 'config_adaptive_timestep', config_adaptive_timestep) + + ! === + ! === Adaptive timestep: update clock information + ! === + ! Set time step in clock object since the time step could have changed + ! Need to get value out of a block, but all blocks should have the same value, so just use first block + if (config_adaptive_timestep) then + block => domain % blocklist + call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) + call mpas_pool_get_array(meshPool, 'deltat', deltat) + ! convert dt in seconds to timeInterval type + call mpas_set_timeInterval(timeStepInterval, dt=deltat, ierr=err_tmp) + err = ior(err,err_tmp) + ! update the clock with the timeInterval + call mpas_set_clock_timestep(domain % clock, timeStepInterval, err_tmp) + err = ior(err,err_tmp) + endif + + ! === + ! === Update clock information + ! === + ! Advance clock - needed to wait until after time step is completed in case the dt has changed! + call mpas_advance_clock(domain % clock) + currTime = mpas_get_clock_time(domain % clock, MPAS_NOW, err_tmp) + call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp, ierr=err_tmp) + err = ior(err, err_tmp) + call mpas_log_write(' Completed timestep. New time is: ' // trim(timeStamp), flushNow=.true.) + + block => domain % blocklist + do while (associated(block)) + ! Assign the time stamp for this time step + call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) + call mpas_pool_get_array(meshPool, 'xtime', xtime) + xtime = timeStamp + + ! compute time since start of simulation, in days + call mpas_pool_get_array(meshPool, 'simulationStartTime', simulationStartTime) + call mpas_pool_get_array(meshPool, 'daysSinceStart',daysSinceStart) + call mpas_set_time(simulationStartTime_timeType, dateTimeString=simulationStartTime) + call mpas_get_timeInterval(currTime - simulationStartTime_timeType, dt=daysSinceStart) + daysSinceStart = daysSinceStart / seconds_per_day + + block => block % next + end do + + ! === error check + if (err > 0) then + call mpas_log_write("An error has occurred in advance_clock.", MPAS_LOG_ERR) + endif + + !-------------------------------------------------------------------- + end subroutine advance_clock + -end module li_time_integration_fe +end module li_time_integration_fe_rk 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 4cbe4da25bc8..077ff4c1773f 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 @@ -616,8 +616,12 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro ! Now interpolate from vertices to cell centers - call li_interpolate_vertex_to_cell_2d(meshPool, dissipationVertexField % array, heatDissipation) - heatDissipation = heatDissipation / (config_ice_density * cp_ice) + if (albanyVelocityError == 0) then + ! the dissipationVertexField can have garbage if the solver didn't converge, + ! so keep previous timestep's field if nonconvergence + call li_interpolate_vertex_to_cell_2d(meshPool, dissipationVertexField % array, heatDissipation) + heatDissipation = heatDissipation / (config_ice_density * cp_ice) + endif call mpas_deallocate_scratch_field(dissipationVertexField, .true.) call li_interpolate_vertex_to_cell_1d(meshPool, drivingStressVert, drivingStress) diff --git a/components/mpas-albany-landice/src/shared/Makefile b/components/mpas-albany-landice/src/shared/Makefile index b9e1f6820813..e0f1828298d9 100644 --- a/components/mpas-albany-landice/src/shared/Makefile +++ b/components/mpas-albany-landice/src/shared/Makefile @@ -3,6 +3,8 @@ OBJS = mpas_li_constants.o \ mpas_li_mask.o \ + mpas_li_mesh.o \ + mpas_li_config.o \ mpas_li_setup.o all: $(OBJS) @@ -13,7 +15,9 @@ mpas_li_setup.o: mpas_li_mask.o: mpas_li_setup.o +mpas_li_mesh.o: +mpas_li_config.o: clean: $(RM) *.o *.mod *.f90 diff --git a/components/mpas-albany-landice/src/shared/mpas_li_config.F b/components/mpas-albany-landice/src/shared/mpas_li_config.F new file mode 100644 index 000000000000..771c89bc6afc --- /dev/null +++ b/components/mpas-albany-landice/src/shared/mpas_li_config.F @@ -0,0 +1,55 @@ +! 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.io/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! li_config +! +!> \brief MALI specific config +!> \details +!> This module contains config specific to MALI. +! +!----------------------------------------------------------------------- + +module li_config + + use mpas_derived_types + use mpas_pool_routines + use mpas_kind_types + + implicit none + public + save + +#include "../inc/config_declare.inc" + +!*********************************************************************** + +contains + +!*********************************************************************** +! +! routine li_config_init +! +!> \brief Initializes the MALI config +!> \details +!> This routine sets up config for use in MALI. +! +!----------------------------------------------------------------------- + subroutine li_config_init(configPool)!{{{ + type (mpas_pool_type), pointer :: configPool + +#include "../inc/config_get.inc" + + end subroutine li_config_init!}}} + +!*********************************************************************** + +end module li_config + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! vim: foldmethod=marker diff --git a/components/mpas-albany-landice/src/shared/mpas_li_mesh.F b/components/mpas-albany-landice/src/shared/mpas_li_mesh.F new file mode 100644 index 000000000000..0c7c78533de0 --- /dev/null +++ b/components/mpas-albany-landice/src/shared/mpas_li_mesh.F @@ -0,0 +1,589 @@ +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! \file mpas_li_mesh.F +! +! 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.io/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! li_mesh +! +!> \brief MALI mesh structure with GPU support +!> \author Rob Aulwes and Phil Jones, modified for MALI +!> by Trevor Hillebrand and Matthew Hoffman +!> \date 14 Jan 2020, added to MALI 2023 +!> \details +!> This module creates and maintains a primary land ice mesh structure +!> and ensures all mesh variables are copied to an accelerator device +!> if needed. Currently it consists of pointers to the existing MPAS mesh pool +!> variables, but is intended to eventually replace the mesh pool later. +! +!------------------------------------------------------------------------------- + +module li_mesh + + ! module dependencies + use mpas_dmpar + use mpas_derived_types + use mpas_pool_routines + use mpas_constants + use mpas_log + + implicit none + private + + !---------------------------------------------------------------------------- + ! + ! Public parameters + ! + !---------------------------------------------------------------------------- + !{{{ + + logical, public :: & + onSphere ! this mesh is on the sphere + + real (kind=RKIND), public :: & + sphereRadius ! radius of sphere for spherical meshes + + integer, public :: &! mesh, array sizes + nCells, &! total number of local (owned+halo) cells in primary + nEdges, &! total number of local edge midpoints + nVertices, &! total number of local cells in dual (cell vertices) + nCellsSolve, &! number of cells owned by the local domain + nEdgesSolve, &! number of edges owned by the local domain + nVerticesSolve, &! number of vertices owned by the local domain + maxEdges, &! largest number of edges any polygon has + maxEdges2, &! 2x the largest number of edges any polygon has + vertexDegree, &! number of cells or edges touching each vertex + nVertLevels ! number of vertical levels + + integer, public, dimension(:), allocatable :: & + nCellsHalo, &! number of owned+halo(n) cells in local domain + nEdgesHalo, &! number of owned+halo(n) edges in local domain + nVerticesHalo ! number of owned+halo(n) vertices in local domain + + integer, public, dimension(:), pointer :: & + nEdgesOnEdge, &! number of edges connected to each edge point + nEdgesOnCell, & ! number of edges associated with each cell center + indexToCellID, &! global ID of each local cell + indexToEdgeID, &! global ID of each local edge + indexToVertexID ! global ID of each local vertex + + integer, public, dimension(:, :), pointer :: & + edgesOnEdge, &! index of edges connected to each edge + cellsOnEdge, &! index of cells connected to each edge + verticesOnEdge, &! index of vertices connected to each edge + cellsOnCell, &! index of cells connected to each cell + edgesOnCell, &! index of edges connected to each cell + verticesOnCell, &! index of vertices connected to each cell + cellsOnVertex, &! index of cells connected to each vertex + edgesOnVertex, & ! index of edges connected to each vertex + edgeSignOnCell, &! sign of edge contributions to a cell + edgeSignOnVertex ! sign of edge contributions to a vertex + + real(kind=RKIND), public, dimension(:), pointer :: & + latCell, &! latitude of cell centers + lonCell, &! longitude of cell centers + xCell, &! Cartesian x coord of cell center + yCell, &! Cartesian y coord of cell center + zCell, &! Cartesian z coord of cell center + latEdge, &! latitude of edge + lonEdge, &! longitude of edge + xEdge, &! Cartesian x coord of edge + yEdge, &! Cartesian y coord of edge + zEdge, &! Cartesian z coord of edge + latVertex, &! latitude of vertex + lonVertex, &! longitude of vertex + xVertex, &! Cartesian coord of vertex + yVertex, &! Cartesian y coord of vertex + zVertex, &! Cartesian z coord of vertex + dcEdge, &! length of edge = dist between cells across edge + dvEdge, &! length of edge = dist between vertices along edge + areaCell, &! area of each cell + areaTriangle, &! area of each cell on dual grid + meshDensity, &! density of mesh + angleEdge ! angle the edge normal makes with local east + + real(kind=RKIND), public, dimension(:), allocatable :: & + invAreaCell ! inverse of area of each cell + + ! Multiplicative masks and vectors for various conditions + integer, public, dimension(:), allocatable :: & + boundaryCell ! mask for boundary cells at each level + + real(kind=RKIND), public, dimension(:, :), pointer :: & + weightsOnEdge, &! weights on each edge + kiteAreasOnVertex, &! real (vertexDegree nVertices) + edgeNormalVectors, &! normal unit vector at edge + localVerticalUnitVectors ! local unit vector iin vertical + + real(kind=RKIND), public, dimension(:, :, :), pointer :: & + cellTangentPlane, &! two vectors defining tangent plane at cell center + coeffs_reconstruct ! coeffs for reconstructing vectors at cell centers + + !}}} + + !---------------------------------------------------------------------------- + ! + ! Public member functions + ! + !---------------------------------------------------------------------------- + !{{{ + + public :: & + li_meshCreate, & + li_meshDestroy + !}}} + +!*********************************************************************** + +contains + +!*********************************************************************** +! +! li_meshCreate +! +!> \brief Creates the ocean mesh data structure on both host and device +!> \author Rob Aulwes and Phil Jones +!> \date 14 Jan 2020 +!> \details +!> This module creates and maintains public ocean mesh data +!> and ensures all mesh variables are copied to an accelerator device +!> if needed. +! +!----------------------------------------------------------------------- + + subroutine li_meshCreate(domain) !{{{ + + ! Input arguments + + type(domain_type) :: & + domain !< [in] MPAS type to describe domain + + ! Local variables + + integer :: & + blockCount ! counter for number of blocks + + type(block_type), pointer :: & + block ! variables in current subblock + + type(mpas_pool_type), pointer :: & + meshPool ! mesh variables in MPAS pool structure + + real (kind=RKIND) :: & + maxDensityLocal, maxDensityGlobal ! temps for mesh density + + ! scalar pointers for retrieval, but convert to actual scalars in struct + logical, pointer :: & + on_a_sphere + + real (kind=RKIND), pointer :: & + sphere_radius + + integer, pointer :: &! mesh dimensions + nCellsTmp, &! + nEdgesTmp, &! + nVerticesTmp, &! + maxEdgesTmp, &! + maxEdges2Tmp, &! + vertexDegreeTmp, &! + nVertLevelsTmp ! + !nVertLevelsP1Tmp ! + + ! temporary pointers for converting index arrays + integer, dimension(:), pointer :: & + nCellsArrayTmp, & + nEdgesArrayTmp, & + nVerticesArrayTmp + + ! temporary pointers for converting masks + integer i, k, n ! loop indices + integer, dimension(:, :), pointer :: & + edgeSignOnCellTmp, & + edgeSignOnVertexTmp + + !*** + !*** end of preamble, begin code + !*** + + ! We only support one block so test for condition here + blockCount = 0 + block => domain%blocklist + do while (associated(block)) + blockCount = blockCount + 1 + if (blockCount > 1) then + call mpas_log_write( & + 'li_meshCreate: more than one block no longer supported', & + MPAS_LOG_CRIT) + endif + block => block%next + end do + + ! Reset to original block + block => domain%blocklist + + ! retrieve the mpas mesh pool + call mpas_pool_get_subpool(block%structs, 'mesh', meshPool) + + !----------------------------------------------------------------- + ! first set pointers/values for all mesh variables + ! many variables already initialized based on read of mesh file + !----------------------------------------------------------------- + + ! set all mesh properties + call mpas_pool_get_config(meshPool, 'on_a_sphere', & + on_a_sphere) + call mpas_pool_get_config(meshPool, 'sphere_radius', & + sphere_radius) + + ! set all mesh dimensions + call mpas_pool_get_dimension(meshPool, 'nCells', & + nCellsTmp) + call mpas_pool_get_dimension(meshPool, 'nEdges', & + nEdgesTmp) + call mpas_pool_get_dimension(meshPool, 'nVertices', & + nVerticesTmp) + call mpas_pool_get_dimension(meshPool, 'maxEdges', & + maxEdgesTmp) + call mpas_pool_get_dimension(meshPool, 'maxEdges2', & + maxEdges2Tmp) + call mpas_pool_get_dimension(meshPool, 'vertexDegree', & + vertexDegreeTmp) + call mpas_pool_get_dimension(meshPool, 'nVertLevels', & + nVertLevelsTmp) + call mpas_pool_get_dimension(meshPool, 'nCellsArray', & + nCellsArrayTmp) + call mpas_pool_get_dimension(meshPool, 'nEdgesArray', & + nEdgesArrayTmp) + call mpas_pool_get_dimension(meshPool, 'nVerticesArray', & + nVerticesArrayTmp) + + ! translate scalar pointers to scalars in new mesh structure + onSphere = on_a_sphere + sphereRadius = sphere_radius + maxEdges = maxEdgesTmp + maxEdges2 = maxEdges2Tmp + vertexDegree = vertexDegreeTmp + nVertLevels = nVertLevelsTmp + + ! convert previous index limits into new halo definitions + nCells = nCellsTmp + nEdges = nEdgesTmp + nVertices = nVerticesTmp + + n = size(nCellsArrayTmp) + allocate (nCellsHalo(n - 1)) + nCellsSolve = nCellsArrayTmp(1) + do i = 2, n + nCellsHalo(i - 1) = nCellsArrayTmp(i) + end do + + n = size(nEdgesArrayTmp) + allocate (nEdgesHalo(n - 1)) + nEdgesSolve = nEdgesArrayTmp(1) + do i = 2, n + nEdgesHalo(i - 1) = nEdgesArrayTmp(i) + end do + + n = size(nVerticesArrayTmp) + allocate (nVerticesHalo(n - 1)) + nVerticesSolve = nVerticesArrayTmp(1) + do i = 2, n + nVerticesHalo(i - 1) = nVerticesArrayTmp(i) + end do + + ! set pointers for a lot of connectivity info + call mpas_pool_get_array(meshPool, 'nEdgesOnEdge', & + nEdgesOnEdge) + call mpas_pool_get_array(meshPool, 'nEdgesOnCell', & + nEdgesOnCell) + call mpas_pool_get_array(meshPool, 'indexToCellID', & + indexToCellID) + call mpas_pool_get_array(meshPool, 'indexToEdgeID', & + indexToEdgeID) + call mpas_pool_get_array(meshPool, 'indexToVertexID', & + indexToVertexID) + call mpas_pool_get_array(meshPool, 'edgesOnEdge', & + edgesOnEdge) + call mpas_pool_get_array(meshPool, 'cellsOnEdge', & + cellsOnEdge) + call mpas_pool_get_array(meshPool, 'verticesOnEdge', & + verticesOnEdge) + call mpas_pool_get_array(meshPool, 'cellsOnCell', & + cellsOnCell) + call mpas_pool_get_array(meshPool, 'edgesOnCell', & + edgesOnCell) + call mpas_pool_get_array(meshPool, 'verticesOnCell', & + verticesOnCell) + call mpas_pool_get_array(meshPool, 'cellsOnVertex', & + cellsOnVertex) + call mpas_pool_get_array(meshPool, 'edgesOnVertex', & + edgesOnVertex) + + ! now set a number of physics and numerical properties of mesh + call mpas_pool_get_array(meshPool, 'latCell', & + latCell) + call mpas_pool_get_array(meshPool, 'lonCell', & + lonCell) + call mpas_pool_get_array(meshPool, 'xCell', & + xCell) + call mpas_pool_get_array(meshPool, 'yCell', & + yCell) + call mpas_pool_get_array(meshPool, 'zCell', & + zCell) + call mpas_pool_get_array(meshPool, 'latEdge', & + latEdge) + call mpas_pool_get_array(meshPool, 'lonEdge', & + lonEdge) + call mpas_pool_get_array(meshPool, 'xEdge', & + xEdge) + call mpas_pool_get_array(meshPool, 'yEdge', & + yEdge) + call mpas_pool_get_array(meshPool, 'zEdge', & + zEdge) + call mpas_pool_get_array(meshPool, 'latVertex', & + latVertex) + call mpas_pool_get_array(meshPool, 'lonVertex', & + lonVertex) + call mpas_pool_get_array(meshPool, 'xVertex', & + xVertex) + call mpas_pool_get_array(meshPool, 'yVertex', & + yVertex) + call mpas_pool_get_array(meshPool, 'zVertex', & + zVertex) + call mpas_pool_get_array(meshPool, 'dcEdge', & + dcEdge) + call mpas_pool_get_array(meshPool, 'dvEdge', & + dvEdge) + call mpas_pool_get_array(meshPool, 'areaCell', & + areaCell) + call mpas_pool_get_array(meshPool, 'areaTriangle', & + areaTriangle) + call mpas_pool_get_array(meshPool, 'weightsOnEdge', & + weightsOnEdge) + call mpas_pool_get_array(meshPool, 'meshDensity', & + meshDensity) + call mpas_pool_get_array(meshPool, 'angleEdge', & + angleEdge) + call mpas_pool_get_array(meshPool, 'weightsOnEdge', & + weightsOnEdge) + call mpas_pool_get_array(meshPool, 'kiteAreasOnVertex', & + kiteAreasOnVertex) + call mpas_pool_get_array(meshPool, 'edgeNormalVectors', & + edgeNormalVectors) + call mpas_pool_get_array(meshPool, 'localVerticalUnitVectors', & + localVerticalUnitVectors) + call mpas_pool_get_array(meshPool, 'cellTangentPlane', & + cellTangentPlane) + call mpas_pool_get_array(meshPool, 'coeffs_reconstruct', & + coeffs_reconstruct) + + ! For masks, we wish to convert to real multiplicative masks + ! so retrieve integer version pointers and allocate real masks + ! Once these are converted in Registry, we can eliminate this. + call mpas_pool_get_array(meshPool, 'edgeSignOnCell', & + edgeSignOnCell) + call mpas_pool_get_array(meshPool, 'edgeSignOnVertex', & + edgeSignOnVertex) + + allocate ( & + boundaryCell(nCells+1), & + invAreaCell(nCells+1)) + + !----------------------------------------------------------------- + ! Now that all pointers are set and mesh variables allocated + ! we initialize other mesh quantities + !----------------------------------------------------------------- + + ! Start by initializing vertical mesh, min/max cells and + ! sign/index fields + call meshSignIndexFields() + + areaCell(nCells+1) = -1.0e34_RKIND + + ! Compute the inverse of areaCell + do n = 1, nCells + invAreaCell(n) = 1.0_RKIND / areaCell(n) + end do + invAreaCell(nCells+1) = 0.0_RKIND + +!------------------------------------------------------------------------------- + + end subroutine li_meshCreate !}}} + +!******************************************************************************* +! +! li_meshDestroy +! +!> \brief Destroy mesh structure and removes from device +!> \author Rob Aulwes and Phil Jones +!> \date 14 Jan 2020 +!> \details +!> This module removes the mesh variables from the device and invalidates +!> all pointers in the mesh structure. +! +!------------------------------------------------------------------------------- + + subroutine li_meshDestroy(err) !{{{ + + ! Input variables + + ! Since the ocnMesh is currently a public module variable, no inputs + ! here, but eventually may want to treat ocnMesh as a specific + ! instantiation instead and pass via args everywhere. If so, need an + ! input mesh here + + ! Output variables + + integer, intent(out) :: & + err ! returned error flag + + ! Local variables + + !*** + !*** end of preamble, begin code + !*** + + err = 0 + + ! Reset all scalars to zero + onSphere = .false. + sphereRadius = 0.0_RKIND + nCells = 0 + nEdges = 0 + nVertices = 0 + nCellsSolve = 0 + nEdgesSolve = 0 + nVerticesSolve = 0 + maxEdges = 0 + maxEdges2 = 0 + vertexDegree = 0 + nVertLevels = 0 + + ! Now nullify all pointers to invalidate fields + ! If this becomes the only mesh structure and mesh pool is eliminated, + ! then we will want to deallocate here instead of nullify. + + nullify (nEdgesOnEdge, & + nEdgesOnCell, & + indexToCellID, & + indexToEdgeID, & + indexToVertexID, & + edgesOnEdge, & + cellsOnEdge, & + verticesOnEdge, & + cellsOnCell, & + edgesOnCell, & + verticesOnCell, & + cellsOnVertex, & + edgesOnVertex, & + latCell, & + lonCell, & + xCell, & + yCell, & + zCell, & + latEdge, & + lonEdge, & + xEdge, & + yEdge, & + zEdge, & + latVertex, & + lonVertex, & + xVertex, & + yVertex, & + zVertex, & + dcEdge, & + dvEdge, & + areaCell, & + areaTriangle, & + meshDensity, & + angleEdge, & + weightsOnEdge, & + kiteAreasOnVertex, & + edgeNormalVectors, & + localVerticalUnitVectors, & + cellTangentPlane, & + edgeSignOnCell, & + edgeSignOnVertex, & + coeffs_reconstruct) + + deallocate (nEdgesHalo, & + nCellsHalo, & + nVerticesHalo, & + boundaryCell, & + invAreaCell) + +!------------------------------------------------------------------------------- + + end subroutine li_meshDestroy !}}} + +!*********************************************************************** +! +! routine li_meshSignIndexFields +! +!> \brief set up sign and index fields +!> \author Doug Jacobsen, Mark Petersen, Todd Ringler +!> \date September 2011 +!> \details +!> This routine initializes edgeSignOnCell, edgeSignOnVertex, and +!> kiteIndexOnCell. +! +!----------------------------------------------------------------------- + + subroutine meshSignIndexFields()!{{{ + + !----------------------------------------------------------------- + ! Local variables + !----------------------------------------------------------------- + integer :: iCell, iEdge, iVertex, i, j + + ! End preamble + !------------- + ! Begin code + + ! Initialize to zero + edgeSignOnCell = 0.0_RKIND + edgeSignOnVertex = 0.0_RKIND + + do iCell = 1, nCells + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + iVertex = verticesOnCell(i, iCell) + + ! Vector points from cell 1 to cell 2 + if (iCell == cellsOnEdge(1, iEdge)) then + edgeSignOnCell(i, iCell) = -1.0_RKIND + else + edgeSignOnCell(i, iCell) = 1.0_RKIND + end if + + end do + end do + + do iVertex = 1, nVertices + do i = 1, vertexDegree + iEdge = edgesOnVertex(i, iVertex) + + ! Vector points from vertex 1 to vertex 2 + if (iVertex == verticesOnEdge(1,iEdge)) then + edgeSignOnVertex(i,iVertex) = -1.0_RKIND + else + edgeSignOnVertex(i,iVertex) = 1.0_RKIND + end if + end do + end do + + !-------------------------------------------------------------------- + + end subroutine meshSignIndexFields !}}} +!*********************************************************************** + +end module li_mesh + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! vim: foldmethod=marker diff --git a/components/mpas-albany-landice/src/shared/mpas_li_setup.F b/components/mpas-albany-landice/src/shared/mpas_li_setup.F index dd6b2f2431c5..e231b4733bbd 100644 --- a/components/mpas-albany-landice/src/shared/mpas_li_setup.F +++ b/components/mpas-albany-landice/src/shared/mpas_li_setup.F @@ -121,7 +121,7 @@ subroutine li_setup_config_options( domain, err ) call mpas_pool_get_config(liConfigs, 'config_tracer_advection', config_tracer_advection) if (( (trim(config_thermal_solver) == 'temperature') .or. & (trim(config_thermal_solver) == 'enthalpy') ) .and. & - (trim(config_tracer_advection) /= 'fo') )then + ( (trim(config_tracer_advection) /= 'fo') .and. (trim(config_tracer_advection) /= 'fct') ) )then call mpas_log_write("Setting config_tracer_advection='fo' because a thermal solver has been selected. " // & "(config_tracer_advection was set to: " // trim(config_tracer_advection) // ")", MPAS_LOG_WARN) config_tracer_advection = 'fo'