Skip to content

Commit 9cb9304

Browse files
jiandewangpjpegion
andauthored
EMC stochastic candidate 20211028 (mom-ocean#1538)
The stochastic physics feature has been added in MOM6. The following are from Phil Pegion: The ocean stochastic physics has been re-coded such that there is a wrapper in config_src/external/OCEAN_stochastic_phyiscs that contains the calls to the external stochastic_physics repository. This has been added to support non-UFS applications of MOM6 where the stochastic_physics repository is not part of the build. The init and run procedures are called from src/core/MOM.F90. I have also created a new control structure stochastic_CS, which contains the logical variables, and random patterns which are then passed into src/parameterizations/vertical/MOM_diabadic_driver.F90 and src/parameterizations/vertical/MOM_energetic-PBL.F90. The writing of the ocean stochastic restarts sit in config_src/nuopc_cap/mom_cap.F90 Co-authored-by: pjpegion <Philip.Pegion@noaa.gov>
1 parent 8d80d64 commit 9cb9304

7 files changed

Lines changed: 350 additions & 23 deletions

File tree

config_src/drivers/nuopc_cap/mom_cap.F90

Lines changed: 8 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -97,6 +97,7 @@ module MOM_cap_mod
9797
use NUOPC_Model, only: model_label_SetRunClock => label_SetRunClock
9898
use NUOPC_Model, only: model_label_Finalize => label_Finalize
9999
use NUOPC_Model, only: SetVM
100+
100101
!$use omp_lib , only : omp_set_num_threads
101102

102103
implicit none; private
@@ -1524,7 +1525,7 @@ subroutine ModelAdvance(gcomp, rc)
15241525
integer :: nc
15251526
type(ESMF_Time) :: MyTime
15261527
integer :: seconds, day, year, month, hour, minute
1527-
character(ESMF_MAXSTR) :: restartname, cvalue
1528+
character(ESMF_MAXSTR) :: restartname, cvalue, stoch_restartname
15281529
character(240) :: msgString
15291530
character(ESMF_MAXSTR) :: casename
15301531
integer :: iostat
@@ -1738,14 +1739,19 @@ subroutine ModelAdvance(gcomp, rc)
17381739
! write the final restart without a timestamp
17391740
if (ESMF_AlarmIsRinging(stop_alarm, rc=rc)) then
17401741
write(restartname,'(A)')"MOM.res"
1742+
write(stoch_restartname,'(A)')"ocn_stoch.res.nc"
17411743
else
17421744
write(restartname,'(A,I4.4,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2)') &
17431745
"MOM.res.", year, month, day, hour, minute, seconds
1746+
write(stoch_restartname,'(A,I4.4,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2,A)') &
1747+
"ocn_stoch.res.", year, month, day, hour, minute, seconds,".nc"
17441748
endif
17451749
call ESMF_LogWrite("MOM_cap: Writing restart : "//trim(restartname), ESMF_LOGMSG_INFO)
17461750

17471751
! write restart file(s)
1748-
call ocean_model_restart(ocean_state, restartname=restartname)
1752+
call ocean_model_restart(ocean_state, restartname=restartname, &
1753+
stoch_restartname=stoch_restartname)
1754+
17491755
endif
17501756

17511757
if (is_root_pe()) then

config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90

Lines changed: 25 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -62,6 +62,7 @@ module MOM_ocean_model_nuopc
6262
use MOM_surface_forcing_nuopc, only : convert_IOB_to_forces, ice_ocn_bnd_type_chksum
6363
use MOM_surface_forcing_nuopc, only : ice_ocean_boundary_type, surface_forcing_CS
6464
use MOM_surface_forcing_nuopc, only : forcing_save_restart
65+
use get_stochy_pattern_mod, only : write_stoch_restart_ocn
6566
use iso_fortran_env, only : int64
6667

6768
#include <MOM_memory.h>
@@ -176,6 +177,10 @@ module MOM_ocean_model_nuopc
176177
!! steps can span multiple coupled time steps.
177178
logical :: diabatic_first !< If true, apply diabatic and thermodynamic
178179
!! processes before time stepping the dynamics.
180+
logical :: do_sppt !< If true, stochastically perturb the diabatic and
181+
!! write restarts
182+
logical :: pert_epbl !< If true, then randomly perturb the KE dissipation and
183+
!! genration termsand write restarts
179184

180185
real :: eps_omesh !< Max allowable difference between ESMF mesh and MOM6
181186
!! domain coordinates
@@ -425,6 +430,17 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
425430
endif
426431

427432
call extract_surface_state(OS%MOM_CSp, OS%sfc_state)
433+
! get number of processors and PE list for stocasthci physics initialization
434+
call get_param(param_file, mdl, "DO_SPPT", OS%do_sppt, &
435+
"If true, then stochastically perturb the thermodynamic "//&
436+
"tendencies of T,S, and h. Amplitude and correlations are "//&
437+
"controlled by the nam_stoch namelist in the UFS model only.", &
438+
default=.false.)
439+
call get_param(param_file, mdl, "PERT_EPBL", OS%pert_epbl, &
440+
"If true, then stochastically perturb the kinetic energy "//&
441+
"production and dissipation terms. Amplitude and correlations are "//&
442+
"controlled by the nam_stoch namelist in the UFS model only.", &
443+
default=.false.)
428444

429445
call close_param_file(param_file)
430446
call diag_mediator_close_registration(OS%diag)
@@ -686,14 +702,17 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
686702
end subroutine update_ocean_model
687703

688704
!> This subroutine writes out the ocean model restart file.
689-
subroutine ocean_model_restart(OS, timestamp, restartname, num_rest_files)
705+
subroutine ocean_model_restart(OS, timestamp, restartname, stoch_restartname, num_rest_files)
690706
type(ocean_state_type), pointer :: OS !< A pointer to the structure containing the
691707
!! internal ocean state being saved to a restart file
692708
character(len=*), optional, intent(in) :: timestamp !< An optional timestamp string that should be
693709
!! prepended to the file name. (Currently this is unused.)
694710
character(len=*), optional, intent(in) :: restartname !< Name of restart file to use
695711
!! This option distinguishes the cesm interface from the
696712
!! non-cesm interface
713+
character(len=*), optional, intent(in) :: stoch_restartname !< Name of restart file to use
714+
!! This option distinguishes the cesm interface from the
715+
!! non-cesm interface
697716
integer, optional, intent(out) :: num_rest_files !< number of restart files written
698717

699718
if (.not.MOM_state_is_synchronized(OS%MOM_CSp)) &
@@ -733,6 +752,11 @@ subroutine ocean_model_restart(OS, timestamp, restartname, num_rest_files)
733752
endif
734753
endif
735754
endif
755+
if (present(stoch_restartname)) then
756+
if (OS%do_sppt .OR. OS%pert_epbl) then
757+
call write_stoch_restart_ocn('RESTART/'//trim(stoch_restartname))
758+
endif
759+
endif
736760

737761
end subroutine ocean_model_restart
738762
! </SUBROUTINE> NAME="ocean_model_restart"
Lines changed: 68 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,68 @@
1+
! The are stubs for ocean stochastic physics
2+
! the fully functional code is available at
3+
! http://github.com/noaa-psd/stochastic_physics
4+
module stochastic_physics
5+
6+
implicit none
7+
8+
private
9+
10+
public :: init_stochastic_physics_ocn
11+
public :: run_stochastic_physics_ocn
12+
13+
contains
14+
15+
!!!!!!!!!!!!!!!!!!!!
16+
subroutine init_stochastic_physics_ocn(delt,geoLonT,geoLatT,nx,ny,nz,pert_epbl_in,do_sppt_in, &
17+
mpiroot, mpicomm, iret)
18+
implicit none
19+
real,intent(in) :: delt !< timestep in seconds between calls to run_stochastic_physics_ocn
20+
integer,intent(in) :: nx !< number of gridpoints in the x-direction of the compute grid
21+
integer,intent(in) :: ny !< number of gridpoints in the y-direction of the compute grid
22+
integer,intent(in) :: nz !< number of gridpoints in the z-direction of the compute grid
23+
real,intent(in) :: geoLonT(nx,ny) !< Longitude in degrees
24+
real,intent(in) :: geoLatT(nx,ny) !< Latitude in degrees
25+
logical,intent(in) :: pert_epbl_in !< logical flag, if true generate random pattern for ePBL perturbations
26+
logical,intent(in) :: do_sppt_in !< logical flag, if true generate random pattern for SPPT perturbations
27+
integer,intent(in) :: mpiroot !< root processor
28+
integer,intent(in) :: mpicomm !< mpi communicator
29+
integer, intent(out) :: iret !< return code
30+
31+
iret=0
32+
if (pert_epbl_in .EQV. .true. ) then
33+
print*,'pert_epbl needs to be false if using the stub'
34+
iret=-1
35+
endif
36+
if (do_sppt_in.EQV. .true. ) then
37+
print*,'do_sppt needs to be false if using the stub'
38+
iret=-1
39+
endif
40+
return
41+
end subroutine init_stochastic_physics_ocn
42+
43+
subroutine run_stochastic_physics_ocn(sppt_wts,t_rp1,t_rp2)
44+
implicit none
45+
real, intent(inout) :: sppt_wts(:,:) !< array containing random weights for SPPT range [0,2]
46+
real, intent(inout) :: t_rp1(:,:) !< array containing random weights for ePBL
47+
!! perturbations (KE generation) range [0,2]
48+
real, intent(inout) :: t_rp2(:,:) !< array containing random weights for ePBL
49+
!! perturbations (KE dissipation) range [0,2]
50+
return
51+
end subroutine run_stochastic_physics_ocn
52+
53+
end module stochastic_physics
54+
55+
module get_stochy_pattern_mod
56+
57+
private
58+
59+
public :: write_stoch_restart_ocn
60+
61+
contains
62+
subroutine write_stoch_restart_ocn(sfile)
63+
64+
character(len=*) :: sfile !< name of restart file
65+
return
66+
end subroutine write_stoch_restart_ocn
67+
68+
end module get_stochy_pattern_mod

src/core/MOM.F90

Lines changed: 9 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -59,6 +59,7 @@ module MOM
5959
use MOM_coord_initialization, only : MOM_initialize_coord
6060
use MOM_diabatic_driver, only : diabatic, diabatic_driver_init, diabatic_CS, extract_diabatic_member
6161
use MOM_diabatic_driver, only : adiabatic, adiabatic_driver_init, diabatic_driver_end
62+
use MOM_stochastics, only : stochastics_init, update_stochastics, stochastic_CS
6263
use MOM_diagnostics, only : calculate_diagnostic_fields, MOM_diagnostics_init
6364
use MOM_diagnostics, only : register_transport_diags, post_transport_diagnostics
6465
use MOM_diagnostics, only : register_surface_diags, write_static_fields
@@ -396,6 +397,7 @@ module MOM
396397
type(ODA_CS), pointer :: odaCS => NULL() !< a pointer to the control structure for handling
397398
!! ensemble model state vectors and data assimilation
398399
!! increments and priors
400+
type(stochastic_CS), pointer :: stoch_CS => NULL() !< a pointer to the stochastics control structure
399401
end type MOM_control_struct
400402

401403
public initialize_MOM, finish_MOM_initialization, MOM_end
@@ -640,6 +642,8 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
640642
call disable_averaging(CS%diag)
641643
endif
642644
endif
645+
! advance the random pattern if stochastic physics is active
646+
if (CS%stoch_CS%do_sppt .OR. CS%stoch_CS%pert_epbl) call update_stochastics(CS%stoch_CS)
643647

644648
if (do_dyn) then
645649
if (G%nonblocking_updates) &
@@ -790,6 +794,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
790794
enddo ; enddo
791795
endif
792796

797+
793798
call step_MOM_dynamics(forces, CS%p_surf_begin, CS%p_surf_end, dt, &
794799
dt_therm_here, bbl_time_int, CS, &
795800
Time_local, Waves=Waves)
@@ -1342,7 +1347,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
13421347
call cpu_clock_begin(id_clock_diabatic)
13431348

13441349
call diabatic(u, v, h, tv, CS%Hml, fluxes, CS%visc, CS%ADp, CS%CDp, dtdia, &
1345-
Time_end_thermo, G, GV, US, CS%diabatic_CSp, OBC=CS%OBC, Waves=Waves)
1350+
Time_end_thermo, G, GV, US, CS%diabatic_CSp, CS%stoch_CS,OBC=CS%OBC, Waves=Waves)
13461351
fluxes%fluxes_used = .true.
13471352

13481353
if (showCallTree) call callTree_waypoint("finished diabatic (step_MOM_thermo)")
@@ -2834,6 +2839,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
28342839
call init_oda(Time, G, GV, CS%diag, CS%odaCS)
28352840
endif
28362841

2842+
! initialize stochastic physics
2843+
call stochastics_init(CS%dt_therm, CS%G, CS%GV, CS%stoch_CS, param_file, diag, Time)
2844+
28372845
!### This could perhaps go here instead of in finish_MOM_initialization?
28382846
! call fix_restart_scaling(GV)
28392847
! call fix_restart_unit_scaling(US)
Lines changed: 144 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,144 @@
1+
!> Top-level module for the MOM6 ocean model in coupled mode.
2+
module MOM_stochastics
3+
4+
! This file is part of MOM6. See LICENSE.md for the license.
5+
6+
! This is the top level module for the MOM6 ocean model. It contains routines
7+
! for initialization, update, and writing restart of stochastic physics. This
8+
! particular version wraps all of the calls for MOM6 in the calls that had
9+
! been used for MOM4.
10+
!
11+
use MOM_diag_mediator, only : register_diag_field, diag_ctrl, time_type
12+
use MOM_grid, only : ocean_grid_type
13+
use MOM_verticalGrid, only : verticalGrid_type
14+
use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe
15+
use MOM_error_handler, only : callTree_enter, callTree_leave
16+
use MOM_file_parser, only : get_param, log_version, close_param_file, param_file_type
17+
use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain
18+
use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain
19+
use MOM_domains, only : root_PE,num_PEs
20+
use MOM_coms, only : Get_PElist
21+
use stochastic_physics, only : init_stochastic_physics_ocn, run_stochastic_physics_ocn
22+
23+
#include <MOM_memory.h>
24+
25+
implicit none ; private
26+
27+
public stochastics_init, update_stochastics
28+
29+
!> This control structure holds parameters for the MOM_stochastics module
30+
type, public:: stochastic_CS
31+
logical :: do_sppt !< If true, stochastically perturb the diabatic
32+
logical :: pert_epbl !< If true, then randomly perturb the KE dissipation and genration terms
33+
integer :: id_sppt_wts = -1 !< Diagnostic id for SPPT
34+
integer :: id_epbl1_wts=-1 !< Diagnostic id for epbl generation perturbation
35+
integer :: id_epbl2_wts=-1 !< Diagnostic id for epbl dissipation perturbation
36+
! stochastic patterns
37+
real, allocatable :: sppt_wts(:,:) !< Random pattern for ocean SPPT
38+
!! tendencies with a number between 0 and 2
39+
real, allocatable :: epbl1_wts(:,:) !< Random pattern for K.E. generation
40+
real, allocatable :: epbl2_wts(:,:) !< Random pattern for K.E. dissipation
41+
type(diag_ctrl), pointer :: diag !< structure used to regulate timing of diagnostic output
42+
type(time_type), pointer :: Time !< Pointer to model time (needed for sponges)
43+
end type stochastic_CS
44+
45+
contains
46+
47+
!! This subroutine initializes the stochastics physics control structure.
48+
subroutine stochastics_init(dt, grid, GV, CS, param_file, diag, Time)
49+
real, intent(in) :: dt !< time step [T ~> s]
50+
type(ocean_grid_type), intent(in) :: grid !< horizontal grid information
51+
type(verticalGrid_type), intent(in) :: GV !< vertical grid structure
52+
type(stochastic_CS), pointer, intent(inout):: CS !< stochastic control structure
53+
type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
54+
type(diag_ctrl), target, intent(inout) :: diag !< structure to regulate diagnostic output
55+
type(time_type), target :: Time !< model time
56+
! Local variables
57+
integer,allocatable :: pelist(:) ! list of pes for this instance of the ocean
58+
integer :: mom_comm ! list of pes for this instance of the ocean
59+
integer :: num_procs ! number of processors to pass to stochastic physics
60+
integer :: iret ! return code from stochastic physics
61+
integer :: me ! my pe
62+
integer :: pe_zero ! root pe
63+
integer :: nx ! number of x-points including halo
64+
integer :: ny ! number of x-points including halo
65+
66+
! This include declares and sets the variable "version".
67+
#include "version_variable.h"
68+
character(len=40) :: mdl = "ocean_stochastics_init" ! This module's name.
69+
70+
call callTree_enter("ocean_model_stochastic_init(), MOM_stochastics.F90")
71+
if (associated(CS)) then
72+
call MOM_error(WARNING, "MOM_stochastics_init called with an "// &
73+
"associated control structure.")
74+
return
75+
else ; allocate(CS) ; endif
76+
77+
CS%diag => diag
78+
CS%Time => Time
79+
80+
! Read all relevant parameters and write them to the model log.
81+
call log_version(param_file, mdl, version, "")
82+
83+
! get number of processors and PE list for stocasthci physics initialization
84+
call get_param(param_file, mdl, "DO_SPPT", CS%do_sppt, &
85+
"If true, then stochastically perturb the thermodynamic "//&
86+
"tendemcies of T,S, amd h. Amplitude and correlations are "//&
87+
"controlled by the nam_stoch namelist in the UFS model only.", &
88+
default=.false.)
89+
call get_param(param_file, mdl, "PERT_EPBL", CS%pert_epbl, &
90+
"If true, then stochastically perturb the kinetic energy "//&
91+
"production and dissipation terms. Amplitude and correlations are "//&
92+
"controlled by the nam_stoch namelist in the UFS model only.", &
93+
default=.false.)
94+
if (CS%do_sppt .OR. CS%pert_epbl) then
95+
num_procs=num_PEs()
96+
allocate(pelist(num_procs))
97+
call Get_PElist(pelist,commID = mom_comm)
98+
pe_zero=root_PE()
99+
nx = grid%ied - grid%isd + 1
100+
ny = grid%jed - grid%jsd + 1
101+
call init_stochastic_physics_ocn(dt,grid%geoLonT,grid%geoLatT,nx,ny,GV%ke, &
102+
CS%pert_epbl,CS%do_sppt,pe_zero,mom_comm,iret)
103+
if (iret/=0) then
104+
call MOM_error(FATAL, "call to init_stochastic_physics_ocn failed")
105+
return
106+
endif
107+
108+
if (CS%do_sppt) allocate(CS%sppt_wts(grid%isd:grid%ied,grid%jsd:grid%jed))
109+
if (CS%pert_epbl) then
110+
allocate(CS%epbl1_wts(grid%isd:grid%ied,grid%jsd:grid%jed))
111+
allocate(CS%epbl2_wts(grid%isd:grid%ied,grid%jsd:grid%jed))
112+
endif
113+
endif
114+
CS%id_sppt_wts = register_diag_field('ocean_model', 'sppt_pattern', CS%diag%axesT1, Time, &
115+
'random pattern for sppt', 'None')
116+
CS%id_epbl1_wts = register_diag_field('ocean_model', 'epbl1_wts', CS%diag%axesT1, Time, &
117+
'random pattern for KE generation', 'None')
118+
CS%id_epbl2_wts = register_diag_field('ocean_model', 'epbl2_wts', CS%diag%axesT1, Time, &
119+
'random pattern for KE dissipation', 'None')
120+
121+
if (is_root_pe()) &
122+
write(*,'(/12x,a/)') '=== COMPLETED MOM STOCHASTIC INITIALIZATION ====='
123+
124+
call callTree_leave("ocean_model_init(")
125+
return
126+
end subroutine stochastics_init
127+
128+
!> update_ocean_model uses the forcing in Ice_ocean_boundary to advance the
129+
!! ocean model's state from the input value of Ocean_state (which must be for
130+
!! time time_start_update) for a time interval of Ocean_coupling_time_step,
131+
!! returning the publicly visible ocean surface properties in Ocean_sfc and
132+
!! storing the new ocean properties in Ocean_state.
133+
subroutine update_stochastics(CS)
134+
type(stochastic_CS), intent(inout) :: CS !< diabatic control structure
135+
call callTree_enter("update_stochastics(), MOM_stochastics.F90")
136+
137+
! update stochastic physics patterns before running next time-step
138+
call run_stochastic_physics_ocn(CS%sppt_wts,CS%epbl1_wts,CS%epbl2_wts)
139+
140+
return
141+
end subroutine update_stochastics
142+
143+
end module MOM_stochastics
144+

0 commit comments

Comments
 (0)