From 47267e88d5a6a24a682025a93f478ad2382eeb56 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 23 Apr 2025 17:24:42 -0400 Subject: [PATCH] +Refactor stochastics for dimensional consistency Refactored MOM_stochastics for dimensional and rotational consistency. This involved correcting the dimensions of a hard-coded area, adding missing conversion factors to the register_diag_field calls for 6 stochastics diagnostics, and adding descriptions of the units in comments describing 7 variables in MOM_stochastics. This commit also adds the new runtime parameter STOCHASTICS_ANSWER_DATE that can be set to a value of 20250701 or higher to use rotationally symmetric expressions. The stochastic physics package itself is external to MOM6 and works in unscaled mks units, so fields passed to it (such as the timestep) need to be unscaled back to mks units. As a part of this change, new unit_scale_type arguments were added to stochastics_init and apply_skeb. By default the answers should be bitwise identical in cases without dimensional consistency testing, but the diagnostic conversion factors should correct the problems with dimensional rescaling. --- .../stochastic_physics/stochastic_physics.F90 | 2 +- src/core/MOM.F90 | 4 +- .../stochastic/MOM_stochastics.F90 | 194 ++++++++++++------ 3 files changed, 131 insertions(+), 69 deletions(-) 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