adiabat: add optional mode-2 pass control and skip redundant solves in RCE loop

Author: Nicholaswogan

clima/cython/AdiabatClimate.pyx

@@ -830,6 +830,15 @@ cdef class AdiabatClimate:
     def __set__(self, bint val):
       wa_pxd.adiabatclimate_prevent_overconvection_set(self._ptr, &val)
 
+  property require_mode2:
+    "bool. If true, require passing through mode 2 when mode 1 converges before optional mode 3 polishing."
+    def __get__(self):
+      cdef bint val
+      wa_pxd.adiabatclimate_require_mode2_get(self._ptr, &val)
+      return bool(val)
+    def __set__(self, bint val):
+      wa_pxd.adiabatclimate_require_mode2_set(self._ptr, &val)
+
   property rtol:
     "float. Relative tolerance of integration."
     def __get__(self):

clima/cython/AdiabatClimate_pxd.pxd

@@ -159,6 +159,8 @@ cdef extern void adiabatclimate_convective_max_boundary_shift_set(AdiabatClimate
 
 cdef extern void adiabatclimate_prevent_overconvection_get(AdiabatClimate *ptr, bint *val)
 cdef extern void adiabatclimate_prevent_overconvection_set(AdiabatClimate *ptr, bint *val)
+cdef extern void adiabatclimate_require_mode2_get(AdiabatClimate *ptr, bint *val)
+cdef extern void adiabatclimate_require_mode2_set(AdiabatClimate *ptr, bint *val)
 
 cdef extern void adiabatclimate_rtol_get(AdiabatClimate *ptr, double *val)
 cdef extern void adiabatclimate_rtol_set(AdiabatClimate *ptr, double *val)

clima/fortran/AdiabatClimate.f90

@@ -996,6 +996,24 @@ subroutine adiabatclimate_prevent_overconvection_set(ptr, val) bind(c)
   c%prevent_overconvection = val
 end subroutine
 
+subroutine adiabatclimate_require_mode2_get(ptr, val) bind(c)
+  use clima, only: AdiabatClimate
+  type(c_ptr), value, intent(in) :: ptr
+  logical(c_bool), intent(out) :: val
+  type(AdiabatClimate), pointer :: c
+  call c_f_pointer(ptr, c)
+  val = c%require_mode2
+end subroutine
+
+subroutine adiabatclimate_require_mode2_set(ptr, val) bind(c)
+  use clima, only: AdiabatClimate
+  type(c_ptr), value, intent(in) :: ptr
+  logical(c_bool), intent(in) :: val
+  type(AdiabatClimate), pointer :: c
+  call c_f_pointer(ptr, c)
+  c%require_mode2 = val
+end subroutine
+
 subroutine adiabatclimate_rtol_get(ptr, val) bind(c)
   use clima, only: AdiabatClimate
   type(c_ptr), value, intent(in) :: ptr

src/adiabat/clima_adiabat.f90

@@ -120,6 +120,9 @@ module clima_adiabat
     integer :: convective_max_boundary_shift = -1
     !> If true, shrink convective tops when a strong inversion exists just above.
     logical :: prevent_overconvection = .true.
+    !> If true, require passing through mode 2 when mode 1 converges
+    !> before optional mode 3 polishing.
+    logical :: require_mode2 = .false.
     !> Per-layer lockout counter used by prevent_overconvection polishing.
     !> A positive value temporarily prevents immediate re-demotion of the same
     !> convective top boundary to avoid ABAB oscillations.

src/adiabat/clima_adiabat_solve.f90

@@ -186,7 +186,7 @@ module function AdiabatClimate_RCE(self, P_i_surf, T_surf_guess, T_guess, convec
     logical, allocatable :: convecting_with_below_save(:,:)
     real(dp), allocatable :: T_in(:), x_init(:), dFdt(:), x_stage(:), x_sol(:), f_sol(:)
     integer :: i, j, k, mode_update
-    logical :: solver_ok, converged1
+    logical :: solver_ok, perform_solve
     logical :: mask_changed
 
     converged = .false.
@@ -226,14 +226,21 @@ module function AdiabatClimate_RCE(self, P_i_surf, T_surf_guess, T_guess, convec
       if (allocated(err)) return
     endif
 
-    converged1 = .false.
+    perform_solve = .true.
+    mode_update = 1
+    if (self%max_rc_iters_convection <= 1) then
+      mode_update = 2
+    endif
+
     do i = 1,self%max_rc_iters
       j = i
 
       if (self%verbose) then
         print"(1x,'Iteration =',i3)", i
       endif
 
+      if (perform_solve) then
+
       if (allocated(x_init)) deallocate(x_init)
       allocate(x_init(size(self%inds_Tx)))
       x_init(1) = self%T_surf
@@ -299,55 +306,74 @@ module function AdiabatClimate_RCE(self, P_i_surf, T_surf_guess, T_guess, convec
       call AdiabatClimate_objective(self, P_i_surf, x_sol, dFdt, f_sol, err)
       if (allocated(err)) return
 
+      endif
+      perform_solve = .true.
+
       ! Save the current convective zones
       convecting_with_below_save = reshape(convecting_with_below_save,shape=[self%nz,i],pad=self%convecting_with_below)
 
-      if (.not.converged1) then
-        if (i < self%max_rc_iters_convection) then
-          ! permit both convective<->radiative flips
-          mode_update = 1
-        else
-          ! permit only radiative->convective growth
-          mode_update = 2
-        endif
-      else
-        ! final polish: trim over-convective tops only
-        mode_update = 3
-      endif
-
       call AdiabatClimate_update_convecting_zones(self, P_i_surf, [self%T_surf, self%T], mode_update, err)
       if (allocated(err)) return
       mask_changed = .not. all(convecting_with_below_save(:,i) .eqv. self%convecting_with_below)
 
-      if (.not.converged1) then
+      ! Change modes and check convergence
+      select case (mode_update)
+      case (1)
         if (.not.mask_changed) then
-          converged1 = .true.
-          if (.not.self%prevent_overconvection) then
-            converged = .true.
-          else
+          if (self%require_mode2) then
+            ! We must pass through mode 2
+            mode_update = 2
+            perform_solve = .false.
+            cycle
+          endif
+          if (self%prevent_overconvection) then
+            ! We must skip to mode 3
+            mode_update = 3
+            perform_solve = .false.
             if (self%verbose) then
               print'(1x,A)','Preventing overconvection'
             endif
-            ! Transition directly into polish once stage-1/2 converges.
-            call AdiabatClimate_update_convecting_zones(self, P_i_surf, [self%T_surf, self%T], 3, err)
-            if (allocated(err)) return
-            mask_changed = .not. all(convecting_with_below_save(:,i) .eqv. self%convecting_with_below)
-            if (.not.mask_changed) converged = .true.
+            cycle
+          endif
+          ! If we get here, then we are converged
+          converged = .true.
+          exit
+        else
+          ! Mask is still changing
+          if (i >= self%max_rc_iters_convection - 1) then
+            ! We skip to mode 2
+            mode_update = 2
+            cycle
           endif
         endif
-      else
-        if (.not.mask_changed) converged = .true.
-      endif
-
-      if (converged) then
-        if (self%verbose) then
-          print'(1x,A)','CONVERGED'
+      case (2)
+        if (.not.mask_changed) then
+          if (self%prevent_overconvection) then
+            ! We move on to mode 3
+            mode_update = 3
+            perform_solve = .false.
+            if (self%verbose) then
+              print'(1x,A)','Preventing overconvection'
+            endif
+            cycle
+          endif
+          ! Otherwise we are converged
+          converged = .true.
+          exit
         endif
-        exit
-      endif
+      case (3)
+        if (.not.mask_changed) then
+          converged = .true.
+          exit
+        endif
+      end select
 
     enddo
 
+    if (converged .and. self%verbose) then
+      print'(1x,A)','CONVERGED'
+    endif
+
     ! Return all information to what is was prior to checking for the root
     call AdiabatClimate_set_convecting_zones(self, convecting_with_below_save(:,j), err)
     if (allocated(err)) return