From c52896d77d3307623a43503b270747f9423724b1 Mon Sep 17 00:00:00 2001
From: Paulo Kubota <paulo.kubota@inpe.br>
Date: Thu, 25 Jun 2026 22:28:42 +0000
Subject: [PATCH] check control  numerical stability

---
 namelist_monan/namelist.atmosphere            |  2 +
 namelist_monan/namelist.atmosphere.TEMPLATE   |  2 +
 src/core_atmosphere/Registry.xml              | 10 ++++
 .../dynamics/mpas_atm_time_integration.F      |  2 +-
 src/core_atmosphere/mpas_atm_core.F           | 46 ++++++++++++++++++-
 5 files changed, 60 insertions(+), 2 deletions(-)

diff --git a/namelist_monan/namelist.atmosphere b/namelist_monan/namelist.atmosphere
index 9b786b74..24c4d7d7 100644
--- a/namelist_monan/namelist.atmosphere
+++ b/namelist_monan/namelist.atmosphere
@@ -6,6 +6,8 @@
     config_split_dynamics_transport = true
     config_number_of_sub_steps = 2
     config_dynamics_split_steps = 3
+    config_coldstart_substeps = 3
+    config_coldstart_steps_to_substep = 30
     config_h_mom_eddy_visc2 = 0.0
     config_h_mom_eddy_visc4 = 0.0
     config_v_mom_eddy_visc2 = 0.0
diff --git a/namelist_monan/namelist.atmosphere.TEMPLATE b/namelist_monan/namelist.atmosphere.TEMPLATE
index 2a102cb9..c060a3e7 100755
--- a/namelist_monan/namelist.atmosphere.TEMPLATE
+++ b/namelist_monan/namelist.atmosphere.TEMPLATE
@@ -6,6 +6,8 @@
     config_split_dynamics_transport = true
     config_number_of_sub_steps = 2
     config_dynamics_split_steps = 3
+    config_coldstart_substeps = 3
+    config_coldstart_steps_to_substep = 30
     config_h_mom_eddy_visc2 = 0.0
     config_h_mom_eddy_visc4 = 0.0
     config_v_mom_eddy_visc2 = 0.0
diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml
index 1f778472..3ce4a832 100644
--- a/src/core_atmosphere/Registry.xml
+++ b/src/core_atmosphere/Registry.xml
@@ -107,6 +107,16 @@
                      description="When config_split_dynamics_transport = T, the number of RK steps per transport step"
                      possible_values="Positive integer values"/>
 
+                <nml_option name="config_coldstart_substeps" type="integer" default_value="1"
+                     units="-"
+                     description="Number of dynamics sub-steps per full timestep during cold-start stabilization window. Each full timestep dt is advanced as substeps of dt/config_coldstart_substeps. Physics timestep is unchanged. Set to 1 to disable (default). Recommended value: 3 for cold-starts from GFS analysis at 10 km resolution."
+                     possible_values="Positive integer values (1 = disabled, 3 recommended for cold-start)"/>
+
+                <nml_option name="config_coldstart_steps_to_substep" type="integer" default_value="0"
+                     units="-"
+                     description="Number of full model timesteps (counting from timestep 1) during which cold-start dynamics substepping is active. Set to 0 to disable (default). Recommended value: 30 for 10 km cold-starts from GFS (covers first 30 min at dt=60s). Must be combined with config_coldstart_substeps > 1 to have any effect."
+                     possible_values="Non-negative integer values (0 = disabled, 20-30 recommended for cold-start)"/>
+
                 <nml_option name="config_h_mom_eddy_visc2" type="real" default_value="0.0" in_defaults="false"
                      units="m^2 s^{-1}"
                      description="$\nabla^2$ eddy viscosity for horizontal diffusion of momentum"
diff --git a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F
index 42c5059a..7cb43548 100644
--- a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F
+++ b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F
@@ -1175,7 +1175,7 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group)
          !requires that the subroutine atm_advance_scalars_mono was called on the third Runge Kutta step, so that a halo
          !update for the scalars at time_levs(1) is applied. A halo update for the scalars at time_levs(2) is done above. 
          if (config_monotonic) then
-            rqvdynten(:,:) = ( scalars_2(index_qv,:,:) - scalars_1(index_qv,:,:) ) / config_dt
+            rqvdynten(:,:) = ( scalars_2(index_qv,:,:) - scalars_1(index_qv,:,:) ) / dt
          else
             rqvdynten(:,:) = 0._RKIND
          end if
diff --git a/src/core_atmosphere/mpas_atm_core.F b/src/core_atmosphere/mpas_atm_core.F
index 4093873b..df0cdc1c 100644
--- a/src/core_atmosphere/mpas_atm_core.F
+++ b/src/core_atmosphere/mpas_atm_core.F
@@ -982,6 +982,17 @@ subroutine atm_do_timestep(domain, dt, itimestep)
       real (kind=RKIND) :: xtime_s
       integer :: ierr
 
+      ! Cold-start substepping — controlled via namelist (nhyd_model record):
+      !   config_coldstart_substeps          : number of dynamics sub-steps per full step (1 = disabled)
+      !   config_coldstart_steps_to_substep  : number of full steps in which substepping is active (0 = disabled)
+      ! Physics always uses the full dt; only dynamics is sub-stepped.
+      integer, pointer :: config_coldstart_substeps
+      integer, pointer :: config_coldstart_steps_to_substep
+      integer :: step, integration_steps, itimestep_local
+      real (kind=RKIND) :: integration_dt
+      type (mpas_pool_type), pointer :: state_substep
+      logical :: substep_active
+
       clock => domain % clock
       mpas_log_info => domain % logInfo
 
@@ -1001,7 +1012,40 @@ subroutine atm_do_timestep(domain, dt, itimestep)
       endif
 #endif
 
-      call atm_timestep(domain, dt, currTime, itimestep, exchange_halo_group)
+      ! --- Cold-start substepping logic (Skamarock/JCSDA, ported for MONAN 1.4.5-rc) ---
+      ! Determine whether this full timestep falls within the substepping window
+      call mpas_pool_get_config(domain % blocklist % configs, 'config_coldstart_substeps',         config_coldstart_substeps)
+      call mpas_pool_get_config(domain % blocklist % configs, 'config_coldstart_steps_to_substep', config_coldstart_steps_to_substep)
+
+      substep_active    = .false.
+      integration_dt    = dt
+      integration_steps = 1
+      itimestep_local   = itimestep
+
+      if (itimestep .le. config_coldstart_steps_to_substep .and. config_coldstart_substeps .gt. 1) then
+         substep_active    = .true.
+         integration_dt    = dt / real(config_coldstart_substeps, kind=RKIND)
+         integration_steps = config_coldstart_substeps
+         itimestep_local   = (itimestep - 1) * config_coldstart_substeps + 1
+      end if
+
+      do step = 1, integration_steps
+
+         if (substep_active) &
+            call mpas_log_write('  [cold-start] substep $i of $i (effective dt = $r s)', &
+                                intArgs=(/step, integration_steps/), realArgs=(/integration_dt/))
+
+         call atm_timestep(domain, integration_dt, currTime, itimestep_local, exchange_halo_group)
+
+         ! Between sub-steps: promote time level 2 → time level 1 for state fields
+         if (step .lt. integration_steps) then
+            call mpas_pool_get_subpool(domain % blocklist % structs, 'state', state_substep)
+            call mpas_pool_shift_time_levels(state_substep)
+            itimestep_local = itimestep_local + 1
+         end if
+
+      end do
+      ! --- End cold-start substepping ---
 
    end subroutine atm_do_timestep