diff --git a/config_src/external/stochastic_physics/stochastic_physics.F90 b/config_src/external/stochastic_physics/stochastic_physics.F90 index 97bcff70d4..196468e317 100644 --- a/config_src/external/stochastic_physics/stochastic_physics.F90 +++ b/config_src/external/stochastic_physics/stochastic_physics.F90 @@ -58,7 +58,7 @@ end subroutine init_stochastic_physics_ocn !> Determines the stochastic physics perturbations. subroutine run_stochastic_physics_ocn(sppt_wts, skeb_wts, t_rp1, t_rp2) real, intent(inout) :: sppt_wts(:,:) !< array containing random weights for SPPT range [0,2] - real, intent(inout) :: skeb_wts(:,:) !< array containing random weights for SKEB + real, intent(inout) :: skeb_wts(:,:) !< array containing random weights for SKEB with units of a length [m] real, intent(inout) :: t_rp1(:,:) !< array containing random weights for ePBL !! perturbations (KE generation) range [0,2] real, intent(inout) :: t_rp2(:,:) !< array containing random weights for ePBL diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 156a397ff6..cbfb9f2d50 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1676,7 +1676,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & fluxes%fluxes_used = .true. if (CS%stoch_CS%do_skeb) then - call apply_skeb(CS%G,CS%GV,CS%stoch_CS,CS%u,CS%v,CS%h,CS%tv,dtdia,Time_end_thermo) + call apply_skeb(G, GV, US, CS%stoch_CS, u, v, h, tv, dtdia, Time_end_thermo) endif if (showCallTree) call callTree_waypoint("finished diabatic (step_MOM_thermo)") @@ -3665,7 +3665,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & endif ! initialize stochastic physics - call stochastics_init(CS%dt_therm, CS%G, CS%GV, CS%stoch_CS, param_file, diag, Time) + call stochastics_init(CS%dt_therm, CS%G, CS%GV, US, CS%stoch_CS, param_file, diag, Time) call callTree_leave("initialize_MOM()") call cpu_clock_end(id_clock_init) diff --git a/src/parameterizations/stochastic/MOM_stochastics.F90 b/src/parameterizations/stochastic/MOM_stochastics.F90 index ddc34fdbaa..886ae1697e 100644 --- a/src/parameterizations/stochastic/MOM_stochastics.F90 +++ b/src/parameterizations/stochastic/MOM_stochastics.F90 @@ -8,22 +8,23 @@ module MOM_stochastics ! particular version wraps all of the calls for MOM6 in the calls that had ! been used for MOM4. ! +use MOM_coms, only : Get_PElist use MOM_debugging, only : hchksum, uvchksum, qchksum use MOM_diag_mediator, only : register_diag_field, diag_ctrl, time_type, post_data use MOM_diag_mediator, only : register_static_field, enable_averages, disable_averaging -use MOM_grid, only : ocean_grid_type -use MOM_variables, only : thermo_var_ptrs use MOM_domains, only : pass_var, pass_vector, CORNER, SCALAR_PAIR -use MOM_verticalGrid, only : verticalGrid_type +use MOM_domains, only : root_PE, num_PEs use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave use MOM_file_parser, only : get_param, log_version, close_param_file, param_file_type -use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain -use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain -use MOM_domains, only : root_PE, num_PEs -use MOM_coms, only : Get_PElist +use MOM_grid, only : ocean_grid_type +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : thermo_var_ptrs +use MOM_verticalGrid, only : verticalGrid_type use MOM_EOS, only : calculate_density, EOS_domain use stochastic_physics, only : init_stochastic_physics_ocn, run_stochastic_physics_ocn +use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain +use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain #include @@ -55,12 +56,18 @@ module MOM_stochastics !! dissipation rate used to set the amplitude of SKEBS [nondim] real, allocatable :: skeb_diss(:,:,:) !< Dissipation rate used to set amplitude of SKEBS [L2 T-3 ~> m2 s-2] !! Index into this at h points. + integer :: answer_date !< The vintage of the order of arithmetic in the stochastics + !! calculations. Values below 20250701 recover the answers from + !! early in 2025, while higher values use expressions that have been + !! refactored for rotational symmetry, including with FMAs enabled. + ! stochastic patterns real, allocatable :: sppt_wts(:,:) !< Random pattern for ocean SPPT - !! tendencies with a number between 0 and 2 - real, allocatable :: skeb_wts(:,:) !< Random pattern for ocean SKEB - real, allocatable :: epbl1_wts(:,:) !< Random pattern for K.E. generation - real, allocatable :: epbl2_wts(:,:) !< Random pattern for K.E. dissipation + !! tendencies with a number between 0 and 2 [nondim] + real, allocatable :: skeb_wts(:,:) !< Random pattern of lengthscales for ocean SKEB in mks units [m] + ! Note that SKEB_wts is set via external code in mks units. + real, allocatable :: epbl1_wts(:,:) !< Random pattern for K.E. generation [nondim] + real, allocatable :: epbl2_wts(:,:) !< Random pattern for K.E. dissipation [nondim] type(time_type), pointer :: Time !< Pointer to model time (needed for sponges) type(diag_ctrl), pointer :: diag=>NULL() !< A structure that is used to regulate the @@ -77,10 +84,11 @@ module MOM_stochastics contains !! This subroutine initializes the stochastics physics control structure. -subroutine stochastics_init(dt, grid, GV, CS, param_file, diag, Time) +subroutine stochastics_init(dt, grid, GV, US, CS, param_file, diag, Time) real, intent(in) :: dt !< time step [T ~> s] type(ocean_grid_type), intent(in) :: grid !< horizontal grid information type(verticalGrid_type), intent(in) :: GV !< vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(stochastic_CS), pointer, intent(inout) :: CS !< stochastic control structure type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters type(diag_ctrl), target, intent(inout) :: diag !< structure to regulate diagnostic output @@ -94,6 +102,7 @@ subroutine stochastics_init(dt, grid, GV, CS, param_file, diag, Time) integer :: pe_zero ! root pe integer :: nxT, nxB ! number of x-points including halo integer :: nyT, nyB ! number of y-points including halo + integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. integer :: i, j, k ! loop indices real :: tmp(grid%isdB:grid%iedB,grid%jsdB:grid%jedB) ! Used to construct tapers integer :: taper_width ! Width (in cells) of the taper that brings the stochastic velocity @@ -153,6 +162,14 @@ subroutine stochastics_init(dt, grid, GV, CS, param_file, diag, Time) "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, "DEFAULT_ANSWER_DATE", default_answer_date, & + "This sets the default value for the various _ANSWER_DATE parameters.", & + default=99991231, do_not_log=.true.) + call get_param(param_file, mdl, "STOCHASTICS_ANSWER_DATE", CS%answer_date, & + "The vintage of the order of arithmetic in the stochastics calculations. "//& + "Values below 20250701 recover the answers from early in 2025, while higher "//& + "values use expressions that have been refactored for rotational symmetry.", & + default=20250101) !### Change to: default=default_answer_date) if (CS%do_sppt .OR. CS%pert_epbl .OR. CS%do_skeb) then num_procs = num_PEs() @@ -163,7 +180,7 @@ subroutine stochastics_init(dt, grid, GV, CS, param_file, diag, Time) nyT = grid%jed - grid%jsd + 1 nxB = grid%iedB - grid%isdB + 1 nyB = grid%jedB - grid%jsdB + 1 - call init_stochastic_physics_ocn(dt, grid%geoLonT, grid%geoLatT, nxT, nyT, GV%ke, & + call init_stochastic_physics_ocn(dt*US%T_to_s, grid%geoLonT, grid%geoLatT, nxT, nyT, GV%ke, & grid%geoLonBu, grid%geoLatBu, nxB, nyB, & CS%pert_epbl, CS%do_sppt, CS%do_skeb, pe_zero, mom_comm, iret) if (iret/=0) then @@ -183,23 +200,23 @@ subroutine stochastics_init(dt, grid, GV, CS, param_file, diag, Time) CS%id_sppt_wts = register_diag_field('ocean_model', 'sppt_pattern', CS%diag%axesT1, Time, & 'random pattern for sppt', 'None') CS%id_skeb_wts = register_diag_field('ocean_model', 'skeb_pattern', CS%diag%axesB1, Time, & - 'random pattern for skeb', 'None') + 'random pattern for skeb', 'm', conversion=1.0) ! SKEB_wts is set in external code in mks units of [m] CS%id_epbl1_wts = register_diag_field('ocean_model', 'epbl1_wts', CS%diag%axesT1, Time, & 'random pattern for KE generation', 'None') CS%id_epbl2_wts = register_diag_field('ocean_model', 'epbl2_wts', CS%diag%axesT1, Time, & 'random pattern for KE dissipation', 'None') CS%id_skebu = register_diag_field('ocean_model', 'skebu', CS%diag%axesCuL, Time, & - 'zonal current perts', 'None') + 'zonal current perts', 'm s-1', conversion=US%L_T_to_m_s) CS%id_skebv = register_diag_field('ocean_model', 'skebv', CS%diag%axesCvL, Time, & - 'zonal current perts', 'None') + 'zonal current perts', 'm s-1', conversion=US%L_T_to_m_s) CS%id_diss = register_diag_field('ocean_model', 'skeb_amp', CS%diag%axesTL, Time, & - 'SKEB amplitude', 'm s-1') + 'SKEB amplitude', 'm s-1', conversion=US%L_T_to_m_s) CS%id_psi = register_diag_field('ocean_model', 'psi', CS%diag%axesBL, Time, & - 'stream function', 'None') + 'stream function', 'm2 s-1', conversion=US%L_T_to_m_s*US%L_to_m) CS%id_skeb_taperu = register_static_field('ocean_model', 'skeb_taper_u', CS%diag%axesCu1, & - 'SKEB taper u', 'None', interp_method='none') + 'SKEB taper u', 'None', conversion=1.0, interp_method='none') CS%id_skeb_taperv = register_static_field('ocean_model', 'skeb_taper_v', CS%diag%axesCv1, & - 'SKEB taper v', 'None', interp_method='none') + 'SKEB taper v', 'None', conversion=1.0, interp_method='none') ! Initialize the "taper" fields. These fields multiply the components of the stochastic ! velocity increment in such a way as to smoothly taper them to zero at land boundaries. @@ -259,37 +276,45 @@ subroutine update_stochastics(CS) call callTree_enter("update_stochastics(), MOM_stochastics.F90") ! update stochastic physics patterns before running next time-step - call run_stochastic_physics_ocn(CS%sppt_wts,CS%skeb_wts,CS%epbl1_wts,CS%epbl2_wts) + call run_stochastic_physics_ocn(CS%sppt_wts, CS%skeb_wts, CS%epbl1_wts, CS%epbl2_wts) call callTree_leave("update_stochastics(), MOM_stochastics.F90") end subroutine update_stochastics -subroutine apply_skeb(grid,GV,CS,uc,vc,thickness,tv,dt,Time_end) +subroutine apply_skeb(grid, GV, US, CS, uc, vc, thickness, tv, dt, Time_end) type(ocean_grid_type), intent(in) :: grid !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(stochastic_CS), intent(inout) :: CS !< stochastic control structure - real, dimension(SZIB_(grid),SZJ_(grid),SZK_(GV)), intent(inout) :: uc !< zonal velocity [L T-1 ~> m s-1] real, dimension(SZI_(grid),SZJB_(grid),SZK_(GV)), intent(inout) :: vc !< meridional velocity [L T-1 ~> m s-1] real, dimension(SZI_(grid),SZJ_(grid),SZK_(GV)), intent(in) :: thickness !< thickness [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< points to thermodynamic fields real, intent(in) :: dt !< time increment [T ~> s] type(time_type), intent(in) :: Time_end !< Time at the end of the interval -! locals - - real, dimension(SZIB_(grid),SZJB_(grid),SZK_(GV)) :: psi !< Streamfunction for stochastic velocity increments - !! [L2 T-1 ~> m2 s-1] - real, dimension(SZIB_(grid),SZJ_(grid) ,SZK_(GV)) :: ustar !< Stochastic u velocity increment [L T-1 ~> m s-1] - real, dimension(SZI_(grid) ,SZJB_(grid),SZK_(GV)) :: vstar !< Stochastic v velocity increment [L T-1 ~> m s-1] - real, dimension(SZI_(grid),SZJ_(grid)) :: diss_tmp !< Temporary array used in smoothing skeb_diss - !! [L2 T-3 ~> m2 s-2] - real, dimension(3,3) :: local_weights !< 3x3 stencil weights used in smoothing skeb_diss - !! [L2 ~> m2] - - real :: shr,ten,tot,kh - integer :: i,j,k,iter + + ! local variables + real, dimension(SZIB_(grid),SZJB_(grid),SZK_(GV)) :: psi !< Streamfunction for stochastic velocity increments + !! [L2 T-1 ~> m2 s-1] + real, dimension(SZIB_(grid),SZJ_(grid) ,SZK_(GV)) :: ustar !< Stochastic u velocity increment [L T-1 ~> m s-1] + real, dimension(SZI_(grid) ,SZJB_(grid),SZK_(GV)) :: vstar !< Stochastic v velocity increment [L T-1 ~> m s-1] + real, dimension(SZI_(grid),SZJ_(grid)) :: diss_tmp !< Temporary array used in smoothing skeb_diss + !! [L2 T-3 ~> m2 s-2] + real, dimension(SZI_(grid),SZJ_(grid)) :: area_wt !< Masked cell areas used in spatial filter [L2 ~> m2] + real, dimension(3,3) :: local_weights !< 3x3 stencil weights used in smoothing skeb_diss + !! [L2 ~> m2] + + real :: shr ! Horizonal shear [T-1 ~> s-1] + real :: ten ! Horizonal tension of the flow [T-1 ~> s-1] + real :: tot ! The magnitude of the combined shear and tension [T-1 ~> s-1] + real :: kh ! A smooothing factor [nondim] + real :: sum_wtd_skeb_diss ! The rotationally symmetric sum of the surrounding values of skeb times + ! the area weights used to filter skeb_diss [L4 T-3 ~> m4 s-3] + real :: sum_area_wts ! A rotationally symmetric sum of the surrounding area weights + ! that are used to filter skeb_diss [L2 ~> m2] + integer :: i, j, k, iter integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state call callTree_enter("apply_skeb(), MOM_stochastics.F90") @@ -302,53 +327,90 @@ subroutine apply_skeb(grid,GV,CS,uc,vc,thickness,tv,dt,Time_end) enddo ; enddo enddo - !kh needs to be scaled + ! kh needs to be scaled + kh = 1.0 !(120*111)**2 + if (CS%answer_date < 20250701) then + do k=1,GV%ke + do j=grid%jsc,grid%jec ; do i=grid%isc,grid%iec + ! Shear in [T-1 ~> s-1] + shr = (vc(i,J,k)-vc(i-1,J,k)) * grid%mask2dCv(i,J)*grid%mask2dCv(i-1,J)*grid%IdxCv(i,J) + & + (uc(I,j,k)-uc(I,j-1,k)) * grid%mask2dCu(I,j)*grid%mask2dCu(I,j-1)*grid%IdyCu(I,j) + ! Tension in [T-1 ~> s-1] + ten = (vc(i,J,k)-vc(i-1,J,k)) * grid%mask2dCv(i,J)*grid%mask2dCv(i-1,J)*grid%IdyCv(i,J) + & + (uc(I,j,k)-uc(I,j-1,k)) * grid%mask2dCu(I,j)*grid%mask2dCu(I,j-1)*grid%IdxCu(I,j) + + tot = sqrt( shr**2 + ten**2 ) * grid%mask2dT(i,j) + CS%skeb_diss(i,j,k) = tot**3 * kh * grid%areaT(i,j) !!**2 + enddo ; enddo + enddo + else ! This version has parentheses to preserve rotational symmetry when FMAs are enabled. + do k=1,GV%ke + do j=grid%jsc,grid%jec ; do i=grid%isc,grid%iec + ! Shear in [T-1 ~> s-1] + shr = ((vc(i,J,k)-vc(i-1,J,k)) * grid%mask2dCv(i,J)*grid%mask2dCv(i-1,J)*grid%IdxCv(i,J)) + & + ((uc(I,j,k)-uc(I,j-1,k)) * grid%mask2dCu(I,j)*grid%mask2dCu(I,j-1)*grid%IdyCu(I,j)) + ! Tension in [T-1 ~> s-1] + ten = ((vc(i,J,k)-vc(i-1,J,k)) * grid%mask2dCv(i,J)*grid%mask2dCv(i-1,J)*grid%IdyCv(i,J)) + & + ((uc(I,j,k)-uc(I,j-1,k)) * grid%mask2dCu(I,j)*grid%mask2dCu(I,j-1)*grid%IdxCu(I,j)) + + tot = sqrt( shr**2 + ten**2 ) * grid%mask2dT(i,j) + CS%skeb_diss(i,j,k) = tot**3 * kh * grid%areaT(i,j) !!**2 + enddo ; enddo + enddo + endif + endif ! Sets CS%skeb_diss in [L2 T-3 ~> m2 s-3] without GM or FrictWork - kh=1!(120*111)**2 - do k=1,GV%ke - do j=grid%jsc,grid%jec ; do i=grid%isc,grid%iec - ! Shear - shr = (vc(i,J,k)-vc(i-1,J,k))*grid%mask2dCv(i,J)*grid%mask2dCv(i-1,J)*grid%IdxCv(i,J)+& - (uc(I,j,k)-uc(I,j-1,k))*grid%mask2dCu(I,j)*grid%mask2dCu(I,j-1)*grid%IdyCu(I,j) - ! Tension - ten = (vc(i,J,k)-vc(i-1,J,k))*grid%mask2dCv(i,J)*grid%mask2dCv(i-1,J)*grid%IdyCv(i,J)+& - (uc(I,j,k)-uc(I,j-1,k))*grid%mask2dCu(I,j)*grid%mask2dCu(I,j-1)*grid%IdxCu(I,j) - - tot = sqrt( shr**2 + ten**2 ) * grid%mask2dT(i,j) - CS%skeb_diss(i,j,k) = tot**3 * kh * grid%areaT(i,j)!!**2 - enddo ; enddo - enddo - endif ! Sets CS%skeb_diss without GM or FrictWork + if (CS%skeb_npass >= 1) then + do j=grid%jsc-2,grid%jec+2 ; do i=grid%isc-2,grid%iec+2 + area_wt(i,j) = grid%mask2dT(i,j)*grid%areaT(i,j) + enddo ; enddo + endif ! smooth dissipation skeb_npass times do iter=1,CS%skeb_npass if (mod(iter,2) == 1) call pass_var(CS%skeb_diss, grid%domain) do k=1,GV%ke + if (CS%answer_date < 20250701) then + ! Do the filter with expressions that do not preserve rotational symmetry. + do j=grid%jsc-1,grid%jec+1 ; do i=grid%isc-1,grid%iec+1 + local_weights(:,:) = area_wt(i-1:i+1,j-1:j+1) + diss_tmp(i,j) = sum(local_weights(:,:)*CS%skeb_diss(i-1:i+1,j-1:j+1,k)) / & + (sum(local_weights) + 1.e-16*US%m_to_L**2) + enddo ; enddo + else + ! This spatial filter preserves rotational symmeetry (including with FMAs), but is + ! mathematically equivalent to the older sum-based form above + do j=grid%jsc-1,grid%jec+1 ; do i=grid%isc-1,grid%iec+1 + sum_area_wts = area_wt(i,j) + & + (((area_wt(i-1,j) + area_wt(i+1,j)) + (area_wt(i,j-1) + area_wt(i,j+1))) + & + ((area_wt(i-1,j-1) + area_wt(i+1,j+1)) + (area_wt(i-1,j+1) + area_wt(i+1,j-1)))) + sum_wtd_skeb_diss = CS%skeb_diss(i,j,k) * area_wt(i+1,j) + & + ((( (CS%skeb_diss(i-1,j,k) * area_wt(i-1,j)) + (CS%skeb_diss(i+1,j,k) * area_wt(i+1,j)) ) + & + ( (CS%skeb_diss(i,j-1,k) * area_wt(i,j-1)) + (CS%skeb_diss(i,j+1,k) * area_wt(i,j+1)) )) + & + (( (CS%skeb_diss(i-1,j-1,k) * area_wt(i-1,j-1)) + (CS%skeb_diss(i-1,j-1,k) * area_wt(i+1,j+1)) ) + & + ( (CS%skeb_diss(i-1,j+1,k) * area_wt(i-1,j+1)) + (CS%skeb_diss(i+1,j-1,k) * area_wt(i+1,j-1)) ))) + diss_tmp(i,j) = sum_wtd_skeb_diss / (sum_area_wts + 1.e-16*US%m_to_L**2) + enddo ; enddo + endif do j=grid%jsc-1,grid%jec+1 ; do i=grid%isc-1,grid%iec+1 - ! This does not preserve rotational symmetry - local_weights = grid%mask2dT(i-1:i+1,j-1:j+1)*grid%areaT(i-1:i+1,j-1:j+1) - diss_tmp(i,j) = sum(local_weights*CS%skeb_diss(i-1:i+1,j-1:j+1,k)) / & - (sum(local_weights) + 1.E-16) - enddo ; enddo - do j=grid%jsc-1,grid%jec+1 ; do i=grid%isc-1,grid%iec+1 - if (grid%mask2dT(i,j)==0.) cycle - CS%skeb_diss(i,j,k) = diss_tmp(i,j) + CS%skeb_diss(i,j,k) = grid%mask2dT(i,j) * diss_tmp(i,j) enddo ; enddo enddo enddo call pass_var(CS%skeb_diss, grid%domain) - ! call hchksum(CS%skeb_diss, "SKEB DISS", grid%HI, haloshift=2) - ! call qchksum(CS%skeb_wts, "SKEB WTS", grid%HI, haloshift=1) + ! call hchksum(CS%skeb_diss, "SKEB DISS", grid%HI, haloshift=2, unscale=US%L_T_to_m_s**2*US%s_to_T) + ! call qchksum(CS%skeb_wts, "SKEB WTS", grid%HI, haloshift=1) ! SKEB_wts comes in from external code in mks units. do k=1,GV%ke do J=grid%jscB-1,grid%jecB ; do I=grid%iscB-1,grid%iecB + ! psi has units of [L2 T-1 ~> m2 s-1] because skeb_wts is in mks units of [m]. psi(I,J,k) = sqrt(0.25 * dt * max((CS%skeb_diss(i ,j ,k) + CS%skeb_diss(i+1,j+1,k)) + & (CS%skeb_diss(i ,j+1,k) + CS%skeb_diss(i+1,j ,k)), 0.) ) & - * CS%skeb_wts(I,J) + * US%m_to_L*CS%skeb_wts(I,J) enddo ; enddo enddo - !call qchksum(psi,"SKEB PSI", grid%HI, haloshift=1) + !call qchksum(psi,"SKEB PSI", grid%HI, haloshift=1, unscale=US%L_T_to_m_s*US%L_to_m) !call pass_var(psi, grid%domain, position=CORNER) do k=1,GV%ke do j=grid%jsc,grid%jec ; do I=grid%iscB,grid%iecB @@ -361,7 +423,7 @@ subroutine apply_skeb(grid,GV,CS,uc,vc,thickness,tv,dt,Time_end) enddo ; enddo enddo - !call uvchksum("SKEB increment [uv]", ustar, vstar, grid%HI) + !call uvchksum("SKEB increment [uv]", ustar, vstar, grid%HI, unscale=US%L_T_to_m_s) call enable_averages(dt, Time_end, CS%diag) if (CS%id_diss > 0) then @@ -393,7 +455,7 @@ end subroutine apply_skeb !! input fields have valid values in the first two halo points upon entry. subroutine smooth_x9_uv(G, field_u, field_v, zero_land) type(ocean_grid_type), intent(in) :: G !< Ocean grid - real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: field_u !< u-point field to be smoothed[arbitrary] + real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: field_u !< u-point field to be smoothed [arbitrary] real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: field_v !< v-point field to be smoothed [arbitrary] logical, optional, intent(in) :: zero_land !< If present and false, return the average !! of the surrounding ocean points when