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.
MALI-Dev · Sep 1, 2023 · 52d9a0b · 52d9a0b
1 parent 8f218c3
commit 52d9a0b
Show file tree

Hide file tree

Showing 3 changed files with 40 additions and 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
@@ -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
@@ -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
@@ -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