NCAR to main (2026-04-22)

config_src/drivers/nuopc_cap/mom_cap.F90

@@ -8,6 +8,7 @@ module MOM_cap_mod
 
 use MOM_domains,              only: get_domain_extent
 use MOM_io,                   only: stdout, io_infra_end
+use MOM_io,                   only: insert_ensemble_appendix
 use mpp_domains_mod,          only: mpp_get_compute_domains
 use mpp_domains_mod,          only: mpp_get_ntile_count, mpp_get_pelist, mpp_get_global_domain
 use mpp_domains_mod,          only: mpp_get_domain_npes
@@ -28,6 +29,7 @@ module MOM_cap_mod
 use MOM_ocean_model_nuopc,    only: ocean_model_init_sfc, ocean_model_flux_init
 use MOM_ocean_model_nuopc,    only: ocean_model_init, update_ocean_model, ocean_model_end
 use MOM_ocean_model_nuopc,    only: get_ocean_grid, get_eps_omesh, query_ocean_state
+use MOM_ocean_model_nuopc,    only: stoch_restart_needed
 use MOM_cap_time,             only: AlarmInit
 use MOM_cap_methods,          only: mom_import, mom_export, mom_set_geomtype, mod2med_areacor
 use MOM_cap_methods,          only: med2mod_areacor, state_diagnose
@@ -1775,7 +1777,7 @@ subroutine ModelAdvance(gcomp, rc)
   character(240)                         :: msgString
   character(ESMF_MAXSTR)                 :: casename
   integer                                :: iostat
-  integer                                :: writeunit
+  integer                                :: rpointer_unit
   type(ESMF_VM)                          :: vm
   integer                                :: n, i
   character(240)                         :: import_timestr, export_timestr
@@ -1994,25 +1996,27 @@ subroutine ModelAdvance(gcomp, rc)
           rpointer_filename = trim(rpointer_filename//timestamp)
         endif
 
-        write(restartname,'(A,".mom6.r",A)') &
-             trim(casename), timestamp
+        write(restartname,'(A,".mom6.r",A)') trim(casename), timestamp
+        write(stoch_restartname,'(A,".mom6.r_stoch",A,".nc")')  trim(casename), timestamp
+
+        call insert_ensemble_appendix(stoch_restartname, ".mom6")
+
         call ESMF_LogWrite("MOM_cap: Writing restart :  "//trim(restartname), ESMF_LOGMSG_INFO)
         ! write restart file(s)
-        call ocean_model_restart(ocean_state, restartname=restartname, num_rest_files=num_rest_files)
+        call ocean_model_restart(ocean_state, restartname=restartname, &
+                stoch_restartname=stoch_restartname, num_rest_files=num_rest_files)
         if (localPet == 0) then
            ! Write name of restart file in the rpointer file - this is currently hard-coded for the ocean
-          open(newunit=writeunit, file=rpointer_filename, form='formatted', status='unknown', iostat=iostat)
+          open(newunit=rpointer_unit, file=rpointer_filename, form='formatted', status='unknown', iostat=iostat)
           if (iostat /= 0) then
             call ESMF_LogSetError(ESMF_RC_FILE_OPEN, &
                  msg=subname//' ERROR opening '//rpointer_filename, line=__LINE__, file=u_FILE_u, rcToReturn=rc)
             return
           endif
-          if (len_trim(inst_suffix) == 0) then
-            write(writeunit,'(a)') trim(restartname)//'.nc'
-          else
-            write(writeunit,'(a)') trim(restartname)//'.'//trim(inst_suffix)//'.nc'
-          endif
 
+          call insert_ensemble_appendix(restartname, ".mom6")
+
+          write(rpointer_unit,'(a)') trim(restartname)//'.nc'
           if (num_rest_files > 1) then
             ! append i.th restart file name to rpointer
             do i=1, num_rest_files-1
@@ -2021,10 +2025,15 @@ subroutine ModelAdvance(gcomp, rc)
               else
                 write(suffix,'("_",I2)') i
               endif
-              write(writeunit,'(a)') trim(restartname) // trim(suffix) // '.nc'
+              write(rpointer_unit,'(a)') trim(restartname) // trim(suffix) // '.nc'
             enddo
           endif
-          close(writeunit)
+
+          if (stoch_restart_needed(ocean_state)) then
+            write(rpointer_unit,'(a)') trim(stoch_restartname)
+          endif
+
+          close(rpointer_unit)
         endif
       else  ! not cesm_coupled
         write(restartname,'(i4.4,2(i2.2),A,3(i2.2),A)') year, month, day,".", hour, minute, seconds, &
@@ -2091,6 +2100,7 @@ subroutine ModelSetRunClock(gcomp, rc)
   type(ESMF_Clock)         :: mclock, dclock
   type(ESMF_Time)          :: mcurrtime, dcurrtime
   type(ESMF_Time)          :: mstoptime, dstoptime
+  type(ESMF_Time)          :: mstoptime_prev ! model stop time before it is updated by this routine
   type(ESMF_TimeInterval)  :: mtimestep, dtimestep
   character(len=128)       :: mtimestring, dtimestring
   character(len=256)       :: cvalue
@@ -2118,7 +2128,8 @@ subroutine ModelSetRunClock(gcomp, rc)
                      stopTime=dstoptime, rc=rc)
   if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
-  call ESMF_ClockGet(mclock, currTime=mcurrtime, timeStep=mtimestep, rc=rc)
+  call ESMF_ClockGet(mclock, currTime=mcurrtime, timeStep=mtimestep, &
+                     stopTime=mstoptime_prev, rc=rc)
   if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
   !--------------------------------
@@ -2242,12 +2253,21 @@ subroutine ModelSetRunClock(gcomp, rc)
     endif
 
     ! create a 1-shot alarm at the driver stop time
-    stop_alarm = ESMF_AlarmCreate(mclock, ringtime=dstopTime, name = "stop_alarm", rc=rc)
-    call ESMF_LogWrite(subname//" Create Stop alarm", ESMF_LOGMSG_INFO)
-    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    if (cesm_coupled) then
+      stop_alarm = ESMF_AlarmCreate(mclock, ringtime=mstoptime_prev, name = "stop_alarm", rc=rc)
+      call ESMF_LogWrite(subname//" Create Stop alarm", ESMF_LOGMSG_INFO)
+      if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
-    call ESMF_TimeGet(dstoptime, timestring=timestr, rc=rc)
-    call ESMF_LogWrite("Stop Alarm will ring at : "//trim(timestr), ESMF_LOGMSG_INFO)
+      call ESMF_TimeGet(mstoptime_prev, timestring=timestr, rc=rc)
+      call ESMF_LogWrite("Stop Alarm will ring at : "//trim(timestr), ESMF_LOGMSG_INFO)
+    else
+      stop_alarm = ESMF_AlarmCreate(mclock, ringtime=dstopTime, name = "stop_alarm", rc=rc)
+      call ESMF_LogWrite(subname//" Create Stop alarm", ESMF_LOGMSG_INFO)
+      if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
+      call ESMF_TimeGet(dstoptime, timestring=timestr, rc=rc)
+      call ESMF_LogWrite("Stop Alarm will ring at : "//trim(timestr), ESMF_LOGMSG_INFO)
+    endif
 
     call outputlog_init(gcomp, mclock, rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return

config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90

@@ -64,6 +64,7 @@ module MOM_ocean_model_nuopc
 use MOM_surface_forcing_nuopc, only : ice_ocean_boundary_type, surface_forcing_CS
 use MOM_surface_forcing_nuopc, only : forcing_save_restart
 use get_stochy_pattern_mod,  only : write_stoch_restart_ocn
+use stochy_data_mod,         only : stoch_restfile
 use iso_fortran_env,           only : int64
 
 #include <MOM_memory.h>
@@ -83,6 +84,7 @@ module MOM_ocean_model_nuopc
 public ocean_public_type_chksum
 public get_ocean_grid, query_ocean_state
 public get_eps_omesh
+public stoch_restart_needed
 
 !> This type is used for communication with other components via the FMS coupler.
 !! The element names and types can be changed only with great deliberation, hence
@@ -180,11 +182,12 @@ module MOM_ocean_model_nuopc
                               !! steps can span multiple coupled time steps.
   logical :: diabatic_first   !< If true, apply diabatic and thermodynamic
                               !! processes before time stepping the dynamics.
-  logical :: do_sppt         !< If true, stochastically perturb the diabatic and
-                             !! write restarts
-  logical :: pert_epbl       !< If true, then randomly perturb the KE dissipation and
-                             !! genration termsand write restarts
-
+  logical :: do_sppt          !< If true, stochastically perturb the diabatic
+                              !! tendencies and write restarts
+  logical :: pert_epbl        !< If true, then randomly perturb the KE dissipation and
+                              !! generation terms and write restarts
+  logical :: do_skeb          !< If true, stochastically perturb the ocean lateral
+                              !! velocity and write restarts
   real :: eps_omesh           !< Max allowable difference between ESMF mesh and MOM6
                               !! domain coordinates
 
@@ -268,6 +271,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
   integer :: secs, days
   type(param_file_type) :: param_file !< A structure to parse for run-time parameters
   logical :: use_temperature
+  integer :: i, k
 
   call callTree_enter("ocean_model_init(), ocean_model_MOM.F90")
   if (associated(OS)) then
@@ -283,6 +287,19 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
   call time_interp_external_init
 
   OS%Time = Time_in
+  if(present(input_restart_file)) then
+      k = index(input_restart_file, ' ')
+      if (k==0) k = len_trim(input_restart_file)
+      i = index(input_restart_file, '.r.')
+      if (i>0) then
+         stoch_restfile = input_restart_file(1:i)//'r_stoch'//input_restart_file(i+2:k)
+
+         if (is_root_pe()) then
+           write(stdout,*) 'input_restart_file =', input_restart_file
+           write(stdout,*) 'stoch_restfile =', stoch_restfile
+         endif
+      endif
+  endif
   call initialize_MOM(OS%Time, Time_init, param_file, OS%dirs, OS%MOM_CSp, &
                       Time_in, offline_tracer_mode=OS%offline_tracer_mode, &
                       input_restart_file=input_restart_file, &
@@ -389,7 +406,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
   ! vertical integrals, since the related 3-d sums are not negligible in cost.
   call allocate_surface_state(OS%sfc_state, OS%grid, use_temperature, &
                               do_integrals=.true., gas_fields_ocn=gas_fields_ocn, &
-                              use_meltpot=use_melt_pot, use_marbl_tracers=OS%use_MARBL)
+                              use_meltpot=use_melt_pot, use_MARBL_tracers=OS%use_MARBL)
 
   call surface_forcing_init(Time_in, OS%grid, OS%US, param_file, OS%diag, &
                             OS%forcing_CSp, OS%restore_salinity, OS%restore_temp, OS%use_waves)
@@ -446,6 +463,10 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
                  "production and dissipation terms.  Amplitude and correlations are "//&
                  "controlled by the nam_stoch namelist in the UFS model only.", &
                  default=.false.)
+  call get_param(param_file, mdl, "DO_SKEB", OS%do_skeb, &
+                 "If true, then stochastically perturb the currents "//&
+                 "using the stochastic kinetic energy backscatter scheme.",&
+                 default=.false.)
 
   call close_param_file(param_file)
   call diag_mediator_close_registration(OS%diag)
@@ -545,9 +566,6 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
 
   if (OS%fluxes%fluxes_used) then
 
-    ! enable_averages() is necessary to post forcing fields to diagnostics
-    call enable_averages(dt_coupling, OS%Time + Ocean_coupling_time_step, OS%diag)
-
     if (do_thermo) &
       call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%fluxes, index_bnds, OS%Time, dt_coupling, &
                                OS%grid, OS%US, OS%forcing_CSp, OS%sfc_state, &
@@ -772,8 +790,8 @@ subroutine ocean_model_restart(OS, timestamp, restartname, stoch_restartname, nu
     endif
   endif
   if (present(stoch_restartname)) then
-    if (OS%do_sppt .OR. OS%pert_epbl) then
-      call write_stoch_restart_ocn('RESTART/'//trim(stoch_restartname))
+    if (stoch_restart_needed(OS)) then
+      call write_stoch_restart_ocn(trim(stoch_restartname))
     endif
   endif
 
@@ -1130,7 +1148,6 @@ end subroutine Ocean_stock_pe
 
 !> Write out checksums for fields from the ocean surface state
 subroutine ocean_public_type_chksum(id, timestep, ocn)
-
   character(len=*),        intent(in) :: id  !< An identifying string for this call
   integer,                 intent(in) :: timestep !< The number of elapsed timesteps
   type(ocean_public_type), intent(in) :: ocn !< A structure containing various publicly
@@ -1156,17 +1173,23 @@ end subroutine ocean_public_type_chksum
 
 subroutine get_ocean_grid(OS, Gridp)
   ! Obtain the ocean grid.
-  type(ocean_state_type) :: OS
-  type(ocean_grid_type) , pointer :: Gridp
+  type(ocean_state_type), intent(in) :: OS
+  type(ocean_grid_type) , pointer, intent(out) :: Gridp
 
   Gridp => OS%grid
   return
 end subroutine get_ocean_grid
 
 !> Returns eps_omesh read from param file
 real function get_eps_omesh(OS)
-  type(ocean_state_type) :: OS
+  type(ocean_state_type), intent(in) :: OS
   get_eps_omesh = OS%eps_omesh; return
 end function
 
+!> Returns true if a stochastic restart file is needed
+logical function stoch_restart_needed(OS)
+  type(ocean_state_type), intent(in) :: OS
+  stoch_restart_needed = OS%do_sppt .OR. OS%pert_epbl .OR. OS%do_skeb
+end function stoch_restart_needed
+
 end module MOM_ocean_model_nuopc

config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90

@@ -83,7 +83,7 @@ module MOM_surface_forcing_nuopc
                                 !! pressure limited by max_p_surf instead of the
                                 !! full atmospheric pressure.  The default is true.
   logical :: use_CFC            !< enables the MOM_CFC_cap tracer package.
-  logical :: use_marbl_tracers  !< enables the MARBL tracer package.
+  logical :: use_MARBL_tracers  !< enables the MARBL tracer package.
   logical :: enthalpy_cpl       !< Controls if enthalpy terms are provided by the coupler or computed
                                 !! internally.
   real :: gust_const            !< constant unresolved background gustiness for ustar [R L Z T-2 ~> Pa]
@@ -328,7 +328,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
   if (fluxes%dt_buoy_accum < 0) then
     call allocate_forcing_type(G, fluxes, water=.true., heat=.true., ustar=.true., &
                                press=.true., fix_accum_bug=.not.CS%ustar_gustless_bug, &
-                               cfc=CS%use_CFC, marbl=CS%use_marbl_tracers, hevap=CS%enthalpy_cpl, &
+                               cfc=CS%use_CFC, marbl=CS%use_MARBL_tracers, hevap=CS%enthalpy_cpl, &
                                tau_mag=.true., ice_ncat=IOB%ice_ncat)
     call safe_alloc_ptr(fluxes%omega_w2x,isd,ied,jsd,jed)
     call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed)
@@ -617,7 +617,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
 
   ! Copy MARBL-specific IOB fields into fluxes; also set some MARBL-specific forcings to other values
   ! (constants, values from netCDF, etc)
-  if (CS%use_marbl_tracers) &
+  if (CS%use_MARBL_tracers) &
     call convert_driver_fields_to_forcings(IOB%atm_fine_dust_flux, IOB%atm_coarse_dust_flux, &
                                            IOB%seaice_dust_flux, IOB%atm_bc_flux, IOB%seaice_bc_flux, &
                                            IOB%nhx_dep, IOB%noy_dep, IOB%atm_co2_prog, IOB%atm_co2_diag, &
@@ -1274,7 +1274,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
   call get_param(param_file, mdl, "USE_CFC_CAP", CS%use_CFC, &
                  default=.false., do_not_log=.true.)
 
-  call get_param(param_file, mdl, "USE_MARBL_TRACERS", CS%use_marbl_tracers, &
+  call get_param(param_file, mdl, "USE_MARBL_TRACERS", CS%use_MARBL_tracers, &
                  default=.false., do_not_log=.true.)
 
   call get_param(param_file, mdl, "ENTHALPY_FROM_COUPLER", CS%enthalpy_cpl, &
@@ -1470,7 +1470,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
 
   call register_forcing_type_diags(Time, diag, US, CS%use_temperature, CS%handles, &
                                    use_berg_fluxes=iceberg_flux_diags, use_waves=use_waves, &
-                                   use_cfcs=CS%use_CFC, use_glc_runoff=glc_runoff_diags)
+                                   use_cfcs=CS%use_CFC, use_MARBL_tracers=CS%use_MARBL_tracers, &
+                                    use_glc_runoff=glc_runoff_diags)
 
   call get_param(param_file, mdl, "ALLOW_FLUX_ADJUSTMENTS", CS%allow_flux_adjustments, &
                  "If true, allows flux adjustments to specified via the "//&
@@ -1485,7 +1486,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
   endif
 
   ! Set up MARBL forcing control structure
-  call MARBL_forcing_init(G, US, param_file, diag, Time, CS%inputdir, CS%use_marbl_tracers, &
+  call MARBL_forcing_init(G, US, param_file, diag, Time, CS%inputdir, CS%use_MARBL_tracers, &
       CS%marbl_forcing_CSp)
 
   if (present(restore_salt)) then ; if (restore_salt) then

config_src/drivers/solo_driver/MOM_surface_forcing.F90

@@ -118,7 +118,7 @@ module MOM_surface_forcing
                              !! rotationally invariant and more likely to be the same between compilers.
   logical :: ustar_gustless_bug   !< If true, include a bug in the time-averaging of the
                                   !! gustless wind friction velocity.
-  logical :: use_marbl_tracers              !< If true, allocate memory for forcing needed by MARBL
+  logical :: use_MARBL_tracers              !< If true, allocate memory for forcing needed by MARBL
   ! if WIND_CONFIG=='scurves' then use the following to define a piecewise scurve profile
   real :: scurves_ydata(20) = 90. !< Latitudes of scurve nodes [degreesN]
   real :: scurves_taux(20) = 0.   !< Zonal wind stress values at scurve nodes [R L Z T-2 ~> Pa]
@@ -289,7 +289,7 @@ subroutine set_forcing(sfc_state, forces, fluxes, day_start, day_interval, G, US
     ! Allocate memory for the mechanical and thermodynamic forcing fields.
     call allocate_mech_forcing(G, forces, stress=.true., ustar=.not.CS%nonBous, press=.true., tau_mag=CS%nonBous)
 
-    call allocate_forcing_type(G, fluxes, ustar=.not.CS%nonBous, marbl=CS%use_marbl_tracers, tau_mag=CS%nonBous, &
+    call allocate_forcing_type(G, fluxes, ustar=.not.CS%nonBous, marbl=CS%use_MARBL_tracers, tau_mag=CS%nonBous, &
                                fix_accum_bug=.not.CS%ustar_gustless_bug)
     if (trim(CS%buoy_config) /= "NONE") then
       if ( CS%use_temperature ) then
@@ -386,7 +386,7 @@ subroutine set_forcing(sfc_state, forces, fluxes, day_start, day_interval, G, US
     endif
   endif
 
-  if (CS%use_marbl_tracers) then
+  if (CS%use_MARBL_tracers) then
     call MARBL_forcing_from_data_override(fluxes, day_center, G, US, CS)
   endif
 
@@ -2150,7 +2150,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, tracer_flow_C
     call read_netCDF_data(filename, 'gustiness', CS%gust, G%Domain, &
                           rescale=US%Pa_to_RLZ_T2*US%L_to_Z) ! units in file should be [Pa]
   endif
-  call get_param(param_file, mdl, "USE_MARBL_TRACERS", CS%use_marbl_tracers, &
+  call get_param(param_file, mdl, "USE_MARBL_TRACERS", CS%use_MARBL_tracers, &
                   default=.false., do_not_log=.true.)
 
 !  All parameter settings are now known.
@@ -2183,7 +2183,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, tracer_flow_C
   endif
 
   ! Set up MARBL forcing control structure
-  call MARBL_forcing_init(G, US, param_file, diag, Time, CS%inputdir, CS%use_marbl_tracers, &
+  call MARBL_forcing_init(G, US, param_file, diag, Time, CS%inputdir, CS%use_MARBL_tracers, &
       CS%marbl_forcing_CSp)
 
   call register_forcing_type_diags(Time, diag, US, CS%use_temperature, CS%handles)