Introduce a logical to bypass the flux conservation error in the sea ice thermodynamics update

sorc/rtofs_hycom.fd/src_2.2.99DHMTi-dist2B_relo_cice_v4.0e/source/ice_init.F90

@@ -113,7 +113,8 @@ subroutine input_data
       use ice_transport_driver, only: advection
       use ice_age, only: tr_iage, restart_age
       use ice_meltpond, only: tr_pond, restart_pond
-      use ice_therm_vertical, only: calc_Tsfc, heat_capacity
+      use ice_therm_vertical, only: calc_Tsfc, heat_capacity, &
+                                    loose_ferrmax, high_ferrmax
       use ice_restoring
 !
 ! !INPUT/OUTPUT PARAMETERS:
@@ -168,7 +169,8 @@ subroutine input_data
 !pgp        ocn_data_dir,   oceanmixed_file, restore_sst,   trestore,       &
         ocn_data_dir,   oceanmixed_file, restore_sst,   atm_netrad,     &
         trestore,       &
-        restore_ice,    insert_ssmi,     insert_sih
+        restore_ice,    insert_ssmi,     insert_sih,                    &
+        loose_ferrmax,  high_ferrmax
 !ajw
 
       namelist /tracer_nml/   &
@@ -274,6 +276,9 @@ subroutine input_data
       insert_ssmi = .false.       ! if true, update concentration using SSMI
       insert_sih  = .false.       ! if true, update thickness from a file
 
+      loose_ferrmax = .false.     ! if true, use a large tolerance for "Flux conservation error"
+      high_ferrmax  = 1.0e+3_dbl_kind ! a large value of heat flux
+
       latpnt(1) =  90._dbl_kind   ! latitude of diagnostic point 1 (deg)
       lonpnt(1) =   0._dbl_kind   ! longitude of point 1 (deg)
       latpnt(2) = -65._dbl_kind   ! latitude of diagnostic point 2 (deg)
@@ -491,6 +496,8 @@ subroutine input_data
       call broadcast_array (lonpnt(1:2),        master_task)
       call broadcast_scalar(runid,              master_task)
       call broadcast_scalar(runtype,            master_task)
+      call broadcast_scalar(loose_ferrmax,      master_task)
+      call broadcast_scalar(high_ferrmax,       master_task)
       if (dbug) & ! else only master_task writes to file
       call broadcast_scalar(nu_diag,            master_task)
       ! tracers
@@ -615,6 +622,21 @@ subroutine input_data
          write(nu_diag,1010) ' insert_ssmi               = ', insert_ssmi
 !ajw
          write(nu_diag,1010) ' insert_sih                = ', insert_sih
+         write(nu_diag,1010) ' loose_ferrmax             = ', loose_ferrmax
+         if (loose_ferrmax) then
+           write(nu_diag,1000) ' high_ferrmax              = ', high_ferrmax
+           print*, ' '
+           print*, ' ******** WARNING      *********** '
+           print*, ' '
+           print*, ' This setting is meant to be used in exceptional situations such as:'
+           print*, ' - Production run (strict delivery time).'
+           print*, ' - When the ice thickness is < 5 cm and no snow and melting rate is high.'
+           print*, ' - The amount of shortwave radiation > 100 W/m^2.'
+           print*, ' - See this issue for more details: https://github.com/NOAA-EMC/RTOFS_GLO/issues/140' 
+           print*, ' '
+           print*, ' ******** WARNING ENDS *********** '
+           print*, ' '
+         end if
          if (trim(atm_data_type) /= 'default') then
             write(nu_diag,*)    ' atm_data_dir              = ', &
                                   trim(atm_data_dir)

sorc/rtofs_hycom.fd/src_2.2.99DHMTi-dist2B_relo_cice_v4.0e/source/ice_therm_vertical.F90

@@ -70,6 +70,11 @@ module ice_therm_vertical
       real (kind=dbl_kind), parameter, private :: &
          ferrmax = 1.0e-3_dbl_kind    ! max allowed energy flux error (W m-2)
                                       ! recommend ferrmax < 0.01 W m-2
+      logical (kind=log_kind), public :: &
+         loose_ferrmax = .false.     ! toggle to bypass strict flux crash
+
+      real (kind=dbl_kind), public :: &
+         high_ferrmax = 1.0e+3_dbl_kind    ! Exceptional (for production) energy flux error (W m-2)
 
       character (char_len) :: stoplabel
 
@@ -1409,6 +1414,8 @@ subroutine temperature_changes (nx_block, ny_block, &
       logical (kind=log_kind) :: &
          all_converged  ! = true when all cells have converged
 
+      real (kind=dbl_kind) :: active_ferrmax ! which value of ferrmax is used?
+
       !-----------------------------------------------------------------
       ! Initialize
       !-----------------------------------------------------------------
@@ -1989,6 +1996,13 @@ subroutine temperature_changes (nx_block, ny_block, &
       !-----------------------------------------------------------------
       ! Check for convergence failures.
       !-----------------------------------------------------------------
+            ! Determine active tolerance on ferr
+            if (loose_ferrmax) then
+              active_ferrmax = high_ferrmax
+            else
+              active_ferrmax = ferrmax
+            end if
+
             if (.not.converged(ij)) then
                write(nu_diag,*) 'Thermo iteration does not converge,', &
                                 'istep1, my_task, i, j:', &
@@ -2012,10 +2026,17 @@ subroutine temperature_changes (nx_block, ny_block, &
                write(nu_diag,*) (Tsn(ij,k),k=1,nslyr)
                write(nu_diag,*) 'Final ice temperatures:'
                write(nu_diag,*) (Tin(ij,k),k=1,nilyr)
-               l_stop = .true.
-               istop = i
-               jstop = j
-               return
+
+               ! If we are within the loose tolerance, simply ignore the convergence failure
+               if ((ferr <= 0.9_dbl_kind * active_ferrmax) .and. (loose_ferrmax)) then
+                 write(nu_diag,*) 'WARNING: Ignored thermo convergence failure (within loose_ferrmax) at i,j:', i, j
+               else
+                 l_stop = .true.
+                 istop = i
+                 jstop = j
+                 return
+               endif
+
             endif
          enddo                  ! ij
       endif                     ! all_converged
@@ -4431,6 +4452,8 @@ subroutine conservation_check_vthermo(nx_block, ny_block, &
          einp        , & ! energy input during timestep (J m-2)
          ferr            ! energy conservation error (W m-2)
 
+      real (kind=dbl_kind) :: active_ferrmax ! which value of ferrmax is used?
+
       !----------------------------------------------------------------
       ! If energy is not conserved, print diagnostics and exit.
       !----------------------------------------------------------------
@@ -4444,7 +4467,15 @@ subroutine conservation_check_vthermo(nx_block, ny_block, &
          einp = (fsurfn(i,j) - flatn(i,j) + fswint(i,j) - fhocnn(i,j) &
                - fsnow(i,j)*Lfresh) * dt
          ferr = abs(efinal(ij)-einit(ij)-einp) / dt
-         if (ferr > ferrmax) then
+
+         ! Determine active tolerance on ferr
+         if (loose_ferrmax) then
+           active_ferrmax = high_ferrmax
+         else
+           active_ferrmax = ferrmax
+         end if
+
+         if (ferr > active_ferrmax) then
             l_stop = .true.
             istop = i
             jstop = j