From 52d9a0be50af6b6b2e2dc83c192a0e4ec6ac7a3f Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Fri, 1 Sep 2023 08:12:09 -0700 Subject: [PATCH] Enable fct thickness advection without tracer advection Enable fct thickness advection without tracer advection and throw error if using unsupported combination of thickness and tracer advection. Also initialize some previously uninitialized allocatable arrays. --- .../src/mode_forward/mpas_li_advection.F | 39 ++++++++++++------- .../mpas_li_advection_fct_shared.F | 17 +++++++- .../mpas_li_time_integration_fe.F | 2 +- 3 files changed, 40 insertions(+), 18 deletions(-) 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 f2cd48a70bd7..df8f1b95ba8a 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 @@ -386,7 +386,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, @@ -484,8 +484,7 @@ subroutine li_advection_thickness_tracers(& tend, advectedTracers, layerThicknessOld, & layerThicknessEdge * layerNormalVelocity, 0 * normalVelocity, dt, & nTracers, activeTracerHorizontalAdvectionEdgeFlux, computeBudgets=.false.)!{{{ - elseif ((trim(config_thickness_advection) .eq. 'fct') .and. & - trim(config_tracer_advection) .eq. 'fct') then + elseif (trim(config_thickness_advection) .eq. 'fct') then ! Call fct routine for thickness first, and use activeTracerHorizontalAdvectionEdgeFlux ! returned by that call as normalThicknessFlux for call to tracer fct call li_tracer_advection_fct_tend(& @@ -497,17 +496,25 @@ subroutine li_advection_thickness_tracers(& ! This does conserve mass: layerThickness(:,:) = layerThickness(:,:) + tend(nTracers,:,:) * dt - ! Call fct for tracers, using activeTracerHorizontalAdvectionEdgeFlux - ! from fct thickness advection as normalThicknessFlux - call li_tracer_advection_fct_tend(& - tend(1:nTracers-1,:,:), advectedTracers(1:nTracers-1,:,:), layerThicknessOld, & - activeTracerHorizontalAdvectionEdgeFlux(nTracers,:,:), 0.0_RKIND * normalVelocity, dt, & - nTracers-1, computeBudgets=.false.) - elseif ((trim(config_thickness_advection) .eq. 'fct') .and. & - trim(config_tracer_advection) .eq. 'fo') then + if (trim(config_tracer_advection) .eq. 'fct') then + ! Call fct for tracers, using activeTracerHorizontalAdvectionEdgeFlux + ! from fct thickness advection as normalThicknessFlux + call li_tracer_advection_fct_tend(& + tend(1:nTracers-1,:,:), advectedTracers(1:nTracers-1,:,:), layerThicknessOld, & + activeTracerHorizontalAdvectionEdgeFlux(nTracers,:,:), 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_tracer_advection = fo is not currently supposed & - with config_thickness_advection = fct", MPAS_LOG_ERR) + 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 @@ -701,7 +708,8 @@ subroutine li_advection_thickness_tracers(& err = ior(err,err_tmp) ! Deallocate arrays for fct - if (trim(config_tracer_advection) .eq. 'fct') then + if ( (trim(config_thickness_advection) .eq. 'fct') .or. & + (trim(config_tracer_advection) .eq. 'fct') ) then deallocate( nAdvCellsForEdge, & advCellsForEdge, & advCoefs, & @@ -1285,7 +1293,8 @@ subroutine tracer_setup(& 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') then + 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) 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 index ff34119dc8fb..ce6748e9ea37 100644 --- 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 @@ -132,9 +132,19 @@ subroutine li_tracer_advection_fct_shared_init(geometryPool, err)!{{{ advCellsForEdge (nAdvCellsMax,nEdges), & advCoefs (nAdvCellsMax,nEdges), & advCoefs3rd (nAdvCellsMax,nEdges), & - advMaskHighOrder(nVertLevels ,nEdges), & - advMask2ndOrder(nVertLevels ,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 @@ -148,6 +158,9 @@ subroutine li_tracer_advection_fct_shared_init(geometryPool, err)!{{{ 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 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.F index ccce755c4501..803fbaa79d5f 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.F @@ -699,7 +699,7 @@ subroutine advection_solver(domain, err) call mpas_pool_get_array(meshPool, 'deltat', deltat) if ( (trim(config_thickness_advection) == 'fo' .and. trim(config_tracer_advection) == 'fo') .or. & - (trim(config_thickness_advection) == 'fct' .and. trim(config_tracer_advection) == 'fct') .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