additions for stochy restarts

config_src/nuopc_driver/mom_cap.F90

@@ -91,7 +91,9 @@ module MOM_cap_mod
 use NUOPC_Model, only: model_label_SetRunClock    => label_SetRunClock
 use NUOPC_Model, only: model_label_Finalize       => label_Finalize
 use NUOPC_Model, only: SetVM
+#ifdef UFS
 use get_stochy_pattern_mod, only: write_stoch_restart_ocn
+#endif
 
 implicit none; private
 
@@ -1587,14 +1589,17 @@ subroutine ModelAdvance(gcomp, rc)
         call ocean_model_restart(ocean_state, restartname=restartname)
 
         ! write stochastic physics restart file if active
+#ifdef UFS
         if (ESMF_AlarmIsRinging(stop_alarm, rc=rc)) then
-           write(restartname,'(A)')"ocn_stoch.res.nc")
+           write(restartname,'(A)')"ocn_stoch.res.nc"
         else
            write(restartname,'(A,I4.4,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2,A)') &
-                "oc_stoch.res.", year, month, day, hour, minute, seconds,".nc"
+                "ocn_stoch.res.", year, month, day, hour, minute, seconds,".nc"
         endif
         call ESMF_LogWrite("MOM_cap: Writing restart :  "//trim(restartname), ESMF_LOGMSG_INFO)
-        call write_stoch_restart_ocn('RESTART/'//trim(timestamp)//'.ocn_stoch.res.nc')
+        if (is_root_pe()) print*,'calling write_stoch_restart_ocn ',trim(restartname)
+        call write_stoch_restart_ocn('RESTART/'//trim(restartname))
+#endif
      endif
 
      if (is_root_pe()) then

config_src/solo_driver/MOM_driver.F90

@@ -74,6 +74,9 @@ program MOM_main
 
   use MOM_wave_interface, only: wave_parameters_CS, MOM_wave_interface_init
   use MOM_wave_interface, only: MOM_wave_interface_init_lite, Update_Surface_Waves
+#ifdef UFS
+  use get_stochy_pattern_mod, only: write_stoch_restart_ocn
+#endif
 
   implicit none
 

src/core/MOM.F90

@@ -142,7 +142,9 @@ module MOM
 use MOM_offline_main,          only : offline_fw_fluxes_into_ocean, offline_fw_fluxes_out_ocean
 use MOM_offline_main,          only : offline_advection_layer, offline_transport_end
 use MOM_ALE,                   only : ale_offline_tracer_final, ALE_main_offline
+#ifdef UFS
 use stochastic_physics,        only : init_stochastic_physics_ocn
+#endif
 
 implicit none ; private
 
@@ -229,8 +231,6 @@ module MOM
   logical :: offline_tracer_mode = .false.
                     !< If true, step_offline() is called instead of step_MOM().
                     !! This is intended for running MOM6 in offline tracer mode
-  logical :: do_stochy = .false.
-                    !< If true, call stochastic physics pattern generator
 
   type(time_type), pointer :: Time   !< pointer to the ocean clock
   real    :: dt                      !< (baroclinic) dynamics time step [T ~> s]
@@ -2361,19 +2361,19 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   call copy_dyngrid_to_MOM_grid(dG_in, G_in, US)
   call destroy_dyn_horgrid(dG_in)
 
+  do_epbl=.false.
+  do_sppt=.false.
+#ifdef UFS
   num_procs=num_PEs()
   allocate(pelist(num_procs))
   call Get_PElist(pelist,commID = mom_comm)
   me=PE_here()
   master=root_PE()
 
   !print*,'callling init_stochastic_physics_ocn',maxval(G%geoLatT)
-  do_epbl=.false.
-  do_sppt=.false.
+  if (master) print*,'about to call init_stochastic_physics'
   call init_stochastic_physics_ocn(CS%dt_therm,G%geoLonT,G%geoLatT,G%ied-G%isd+1,G%jed-G%jsd+1,nz,do_epbl,do_sppt,master,mom_comm,iret)
-  if (do_sppt .eq. .true.) CS%do_stochy=.true.
-  if (do_epbl .eq. .true.) CS%do_stochy=.true.
-  !print*,'back from init_stochastic_physics_ocn',CS%do_stochy
+#endif
 
   ! Set a few remaining fields that are specific to the ocean grid type.
   call set_first_direction(G, first_direction)

src/parameterizations/vertical/MOM_diabatic_driver.F90

@@ -69,7 +69,9 @@ module MOM_diabatic_driver
 use MOM_verticalGrid,        only : verticalGrid_type, get_thickness_units
 use MOM_wave_speed,          only : wave_speeds
 use MOM_wave_interface,      only : wave_parameters_CS
+#ifdef UFS
 use stochastic_physics,      only : run_stochastic_physics_ocn
+#endif
 
 
 implicit none ; private
@@ -292,23 +294,28 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
 
   real, dimension(SZI_(G),SZJ_(G))            :: sppt_wts
   real, dimension(SZI_(G),SZJ_(G),2)          :: t_rp
+#ifdef UFS
   real, dimension(SZI_(G),SZJ_(G),SZK_(G))   :: h_in
   real, dimension(SZI_(G),SZJ_(G),SZK_(G))   :: t_in        !< thickness [H ~> m or kg m-2]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G))   :: s_in        !< thickness [H ~> m or kg m-2]
   real :: t_tend,s_tend,h_tend                  ! holder for tendencey needed for SPPT
   real :: t_pert,s_pert,h_pert                  ! holder for tendencey needed for SPPT
+#endif
 
   if (G%ke == 1) return
 
+#ifdef UFS
    ! save  copy of the date for SPPT
   if (CS%do_sppt) then
      h_in=h
      t_in=tv%T
      s_in=tv%S
    endif
+   print*,'calling run_stochastic_physics'
    call run_stochastic_physics_ocn(t_rp,sppt_wts)
    !print*,'in diabatic',CS%do_sppt,size(t_in,1),size(t_in,2),size(t_in,3),size(sppt_wts,1),size(sppt_wts,2)
-   !print*,'in diabatic',CS%do_sppt,minval(sppt_wts),maxval(sppt_wts)
+   print*,'in diabatic',CS%do_sppt,minval(sppt_wts),maxval(sppt_wts)
+   print*,'in diabatic',CS%do_sppt,minval(t_rp),maxval(t_rp)
   if (CS%id_t_rp1 > 0) then
      call post_data(CS%id_t_rp1, t_rp(:,:,1), CS%diag)
   endif
@@ -318,6 +325,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
   if (CS%id_sppt_wts > 0) then
      call post_data(CS%id_sppt_wts, sppt_wts, CS%diag)
   endif
+#endif
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
 
@@ -465,6 +473,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
 
   if (CS%debugConservation) call MOM_state_stats('leaving diabatic', u, v, h, tv%T, tv%S, G, GV, US)
 
+#ifdef UFS
   if (CS%do_sppt) then
     do k=1,nz
       do j=js,je
@@ -488,6 +497,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
       enddo
     enddo
   endif
+#endif
 
 end subroutine diabatic
 
@@ -891,7 +901,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, t_rp, visc, ADp, CDp, d
     endif
 
     call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
-    call energetic_PBL(h, u_h, v_h, tv, fluxes, t_rp, dt, Kd_ePBL, G, GV, US, &
+    call energetic_PBL(h, u_h, v_h, tv, fluxes, t_rp, CS%do_epbl, dt, Kd_ePBL, G, GV, US, &
                        CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves)
 
     if (associated(Hml)) then
@@ -1623,7 +1633,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, t_rp, visc, ADp, CDp, dt, Time
     endif
 
     call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
-    call energetic_PBL(h, u_h, v_h, tv, fluxes, t_rp, dt, Kd_ePBL, G, GV, US, &
+    call energetic_PBL(h, u_h, v_h, tv, fluxes, t_rp, CS%do_epbl, dt, Kd_ePBL, G, GV, US, &
                        CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves)
 
     if (associated(Hml)) then

src/parameterizations/vertical/MOM_energetic_PBL.F90

@@ -246,7 +246,7 @@ module MOM_energetic_PBL
 !!  mixed layer model.  It assumes that heating, cooling and freshwater fluxes
 !!  have already been applied.  All calculations are done implicitly, and there
 !!  is no stability limit on the time step.
-subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, t_rp, dt, Kd_int, G, GV, US, CS, &
+subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, t_rp, stoch_epbl, dt, Kd_int, G, GV, US, CS, &
                          dSV_dT, dSV_dS, TKE_forced, buoy_flux, dt_diag, last_call, &
                          dT_expected, dS_expected, Waves )
   type(ocean_grid_type),   intent(inout) :: G      !< The ocean's grid structure.
@@ -285,6 +285,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, t_rp, dt, Kd_int, G, GV,
                                                    !! call to mixedlayer_init.
   real, dimension(SZI_(G),SZJ_(G),2), &
                            intent(in)    :: t_rp !< random pattern to perturb wind
+  logical,                 intent(in)    :: stoch_epbl !< flag to pertrub production and dissipation of TKE
   real, dimension(SZI_(G),SZJ_(G)), &
                            intent(in)    :: buoy_flux !< The surface buoyancy flux [Z2 T-3 ~> m2 s-3].
   real,          optional, intent(in)    :: dt_diag   !< The diagnostic time step, which may be less
@@ -438,8 +439,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, t_rp, dt, Kd_int, G, GV,
       do K=1,nz+1 ; Kd(K) = 0.0 ; enddo
 
       ! Make local copies of surface forcing and process them.
-      !print*,'PJP EPBL',minval(t_rp),maxval(t_rp)
-      u_star = fluxes%ustar(i,j)!*t_rp(i,j)
+      u_star = fluxes%ustar(i,j)
       u_star_Mean = fluxes%ustar_gustless(i,j)
       B_flux = buoy_flux(i,j)
       if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then
@@ -463,7 +463,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, t_rp, dt, Kd_int, G, GV,
 
       call ePBL_column(h, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, TKE_forcing, B_flux, absf, &
                        u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, &
-                       US, CS, eCD, t_rp(i,j,1),t_rp(i,j,2), dt_diag=dt_diag, Waves=Waves, G=G, i=i, j=j)
+                       US, CS, eCD, t_rp(i,j,1),t_rp(i,j,2), stoch_epbl, dt_diag=dt_diag, Waves=Waves, G=G, i=i, j=j)
 
       ! applly stochastic perturbation to TKE generation
 
@@ -548,7 +548,7 @@ end subroutine energetic_PBL
 !!  mixed layer model for a single column of water.
 subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, absf, &
                        u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, US, CS, eCD, &
-                       t_rp1,t_rp2, dt_diag, Waves, G, i, j)
+                       t_rp1,t_rp2, stoch_epbl, dt_diag, Waves, G, i, j)
   type(verticalGrid_type), intent(in)    :: GV     !< The ocean's vertical grid structure.
   type(unit_scale_type),   intent(in)    :: US     !< A dimensional unit scaling type
   real, dimension(SZK_(GV)), intent(in)  :: h      !< Layer thicknesses [H ~> m or kg m-2].
@@ -587,7 +587,8 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs
                                                    !! call to mixedlayer_init.
   type(ePBL_column_diags), intent(inout) :: eCD    !< A container for passing around diagnostics.
   real,                    intent(in)    :: t_rp1  !< random value to perturb TKE production
-  real,                    intent(in)    :: t_rp2  !< random value to perturb TKE production
+  real,                    intent(in)    :: t_rp2  !< random value to perturb TKE dissipation
+  logical,                 intent(in)    :: stoch_epbl !< flag to pertrub production and dissipation of TKE
   real,          optional, intent(in)    :: dt_diag   !< The diagnostic time step, which may be less
                                                    !! than dt if there are two calls to mixedlayer [T ~> s].
   type(wave_parameters_CS), &
@@ -888,8 +889,8 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs
       else
         mech_TKE = MSTAR_total * (dt*GV%Rho0* u_star**3)
       endif
-      ! stochastically pertrub mech_TKE
-         mech_TKE=mech_TKE*t_rp1
+      ! stochastically pertrub mech_TKE in the UFS
+      if (stoch_epbl) mech_TKE=mech_TKE*t_rp1
 
       if (CS%TKE_diagnostics) then
         eCD%dTKE_conv = 0.0 ; eCD%dTKE_mixing = 0.0
@@ -973,7 +974,7 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs
         if (CS%TKE_diagnostics) &
           !eCD%dTKE_mech_decay = eCD%dTKE_mech_decay + (exp_kh-1.0) * mech_TKE * I_dtdiag
           eCD%dTKE_mech_decay = exp_kh
-        mech_TKE = mech_TKE * (1+(exp_kh-1) * t_rp2)
+        if (stoch_epbl) mech_TKE = mech_TKE * (1+(exp_kh-1) * t_rp2)
 
         !   Accumulate any convectively released potential energy to contribute
         ! to wstar and to drive penetrating convection.