From fa1a7624b4350b7fd12c5e6af10ad6a8efcfcb09 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 6 Nov 2018 18:56:27 -0500 Subject: [PATCH 1/6] +Added run-time parameters for benchmark test case Added two new runtime parameters that are used in setting up the benchmark test case. By default all answers are bitwise identical, but the MOM_parameter_doc.short files change for the benchmark test case. --- src/user/benchmark_initialization.F90 | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/user/benchmark_initialization.F90 b/src/user/benchmark_initialization.F90 index b3b9662164..d50fe3fa05 100644 --- a/src/user/benchmark_initialization.F90 +++ b/src/user/benchmark_initialization.F90 @@ -38,8 +38,8 @@ subroutine benchmark_initialize_topography(D, G, param_file, max_depth) real :: D0 ! A constant to make the maximum ! ! basin depth MAXIMUM_DEPTH. ! real :: x, y -! This include declares and sets the variable "version". -#include "version_variable.h" + ! This include declares and sets the variable "version". +# include "version_variable.h" character(len=40) :: mdl = "benchmark_initialize_topography" ! This subroutine's name. integer :: i, j, is, ie, js, je, isd, ied, jsd, jed is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -107,12 +107,21 @@ subroutine benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state ! interface temperature for a given z and its derivative. real :: pi, z logical :: just_read + ! This include declares and sets the variable "version". +# include "version_variable.h" character(len=40) :: mdl = "benchmark_initialize_thickness" ! This subroutine's name. integer :: i, j, k, k1, is, ie, js, je, nz, itt is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + if (.not.just_read) call log_version(param_file, mdl, version, "") + call get_param(param_file, mdl, "BENCHMARK_ML_DEPTH_IC", ML_depth, & + "Initial mixed layer depth in the benchmark test case.", & + units='m', default=50.0, scale=US%m_to_Z, do_not_log=just_read) + call get_param(param_file, mdl, "BENCHMARK_THERMOCLINE_SCALE", thermocline_scale, & + "Initial thermocline depth scale in the benchmark test case.", & + default=500.0, units="m", scale=US%m_to_Z, do_not_log=just_read) if (just_read) return ! This subroutine has no run-time parameters. @@ -120,8 +129,6 @@ subroutine benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state k1 = GV%nk_rho_varies + 1 - ML_depth = 50.0 * US%m_to_Z - thermocline_scale = 500.0 * US%m_to_Z a_exp = 0.9 ! This block calculates T0(k) for the purpose of diagnosing where the From 6ccb6b6f7a5a94836ba2d3bbad9384c3ced3dcd5 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 6 Nov 2018 18:58:32 -0500 Subject: [PATCH 2/6] +Added a run-time parameter for circle_obcs Added a new runtime parameter, DISK_IC_AMPLITUDE, that is used in setting up the circle_obcs test case. By default all answers are bitwise identical, but the MOM_parameter_doc.short files change for the benchmark test case. --- src/user/circle_obcs_initialization.F90 | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/user/circle_obcs_initialization.F90 b/src/user/circle_obcs_initialization.F90 index f72a6e1830..7ba02a7acc 100644 --- a/src/user/circle_obcs_initialization.F90 +++ b/src/user/circle_obcs_initialization.F90 @@ -37,10 +37,11 @@ subroutine circle_obcs_initialize_thickness(h, G, GV, param_file, just_read_para ! negative because it is positive upward. real :: eta1D(SZK_(GV)+1)! Interface height relative to the sea surface ! positive upward, in in depth units (Z). + real :: IC_amp ! The amplitude of the initial height displacement, in H. real :: diskrad, rad, xCenter, xRadius, lonC, latC, xOffset logical :: just_read -! This include declares and sets the variable "version". -#include "version_variable.h" + ! This include declares and sets the variable "version". +# include "version_variable.h" character(len=40) :: mdl = "circle_obcs_initialization" ! This module's name. integer :: i, j, k, is, ie, js, je, nz @@ -61,6 +62,10 @@ subroutine circle_obcs_initialize_thickness(h, G, GV, param_file, just_read_para "The x-offset of the initially elevated disk in the \n"//& "circle_obcs test case.", units=G%x_axis_units, & default = 0.0, do_not_log=just_read) + call get_param(param_file, mdl, "DISK_IC_AMPLITUDE", IC_amp, & + "Initial amplitude of interface height displacements \n"//& + "in the circle_obcs test case.", & + units='m', default=5.0, scale=GV%m_to_H, do_not_log=just_read) if (just_read) return ! All run-time parameters have been read, so return. @@ -93,12 +98,11 @@ subroutine circle_obcs_initialize_thickness(h, G, GV, param_file, just_read_para rad = rad*(2.*asin(1.)) ! Map 0-1 to 0-pi if (nz==1) then ! The model is barotropic - h(i,j,k) = h(i,j,k) + GV%m_to_H * 1.0*0.5*(1.+cos(rad)) ! cosine bell + h(i,j,k) = h(i,j,k) + IC_amp * 0.5*(1.+cos(rad)) ! cosine bell else ! The model is baroclinic do k = 1, nz - h(i,j,k) = h(i,j,k) - GV%m_to_H * 0.5*(1.+cos(rad)) & ! cosine bell - * 5.0 * real( 2*k-nz ) + h(i,j,k) = h(i,j,k) - 0.5*(1.+cos(rad)) * IC_amp * real( 2*k-nz ) enddo endif enddo ; enddo From 0874ac4fbf447b5e2ba9ef134c145f97b58056be Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 6 Nov 2018 18:59:05 -0500 Subject: [PATCH 3/6] +Added the INTERFACE_IC_QUANTA runtime parameter Added a new parameter, INTERFACE_IC_QUANTA, that is used in the seamount and dumbbell test case to ensure that the initial interface heights can be easily represented in output files and are not sensitive to order of arithmetic changes in the setup. By default all answers are bitwise identical, but the MOM_parameter_doc.short files change for the seamount and dumbbell test cases. --- src/user/dumbbell_initialization.F90 | 20 ++++++++++++++------ src/user/seamount_initialization.F90 | 21 ++++++++++++++------- 2 files changed, 28 insertions(+), 13 deletions(-) diff --git a/src/user/dumbbell_initialization.F90 b/src/user/dumbbell_initialization.F90 index ce8fcdef39..5a5a845449 100644 --- a/src/user/dumbbell_initialization.F90 +++ b/src/user/dumbbell_initialization.F90 @@ -6,7 +6,7 @@ module dumbbell_initialization use MOM_domains, only : sum_across_PEs use MOM_dyn_horgrid, only : dyn_horgrid_type use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, is_root_pe -use MOM_file_parser, only : get_param, param_file_type +use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_get_input, only : directories use MOM_grid, only : ocean_grid_type use MOM_sponge, only : set_up_sponge_field, initialize_sponge, sponge_CS @@ -86,12 +86,14 @@ subroutine dumbbell_initialize_thickness ( h, G, GV, US, param_file, just_read_p ! negative because it is positive upward. real :: eta1D(SZK_(G)+1)! Interface height relative to the sea surface ! positive upward, in depth units (Z). - integer :: i, j, k, is, ie, js, je, nz - real :: x - real :: delta_h - real :: min_thickness, S_surf, S_range, S_ref, S_light, S_dense + real :: min_thickness ! The minimum layer thicknesses, in Z. + real :: S_surf, S_range, S_ref, S_light, S_dense ! Various salinities, in ppt. + real :: eta_IC_quanta ! The granularity of quantization of intial interface heights, in Z-1. + ! This include declares and sets the variable "version". +# include "version_variable.h" character(len=20) :: verticalCoordinate logical :: just_read ! If true, just read parameters but set nothing. + integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -100,6 +102,7 @@ subroutine dumbbell_initialize_thickness ( h, G, GV, US, param_file, just_read_p if (.not.just_read) & call MOM_mesg("MOM_initialization.F90, initialize_thickness_uniform: setting thickness") + if (.not.just_read) call log_version(param_file, mdl, version, "") call get_param(param_file, mdl,"MIN_THICKNESS", min_thickness, & 'Minimum thickness for layer',& units='m', default=1.0e-3, do_not_log=just_read, scale=US%m_to_Z) @@ -125,6 +128,10 @@ subroutine dumbbell_initialize_thickness ( h, G, GV, US, param_file, just_read_p call get_param(param_file, mdl, "S_REF", S_ref, default=35.0, do_not_log=.true.) call get_param(param_file, mdl, "TS_RANGE_S_LIGHT", S_light, default = S_Ref, do_not_log=.true.) call get_param(param_file, mdl, "TS_RANGE_S_DENSE", S_dense, default = S_Ref, do_not_log=.true.) + call get_param(param_file, mdl, "INTERFACE_IC_QUANTA", eta_IC_quanta, & + "The granularity of initial interface height values \n"//& + "per meter, to avoid sensivity to order-of-arithmetic changes.", & + default=2048.0, units="m-1", scale=US%Z_to_m, do_not_log=just_read) if (just_read) return ! All run-time parameters have been read, so return. do K=1,nz+1 @@ -137,7 +144,8 @@ subroutine dumbbell_initialize_thickness ( h, G, GV, US, param_file, just_read_p e0(K) = - G%max_depth * ( ( S_light - S_surf ) + ( S_dense - S_light ) * & ( (real(K)-1.5) / real(nz-1) ) ) / S_range ! Force round numbers ... the above expression has irrational factors ... - e0(K) = nint(2048.*US%Z_to_m*e0(K)) / (2048.*US%Z_to_m) + if (eta_IC_quanta > 0.0) & + e0(K) = nint(eta_IC_quanta*e0(K)) / eta_IC_quanta e0(K) = min(real(1-K)*GV%Angstrom_Z, e0(K)) ! Bound by surface e0(K) = max(-G%max_depth, e0(K)) ! Bound by bottom enddo diff --git a/src/user/seamount_initialization.F90 b/src/user/seamount_initialization.F90 index 42ea66d0ad..3a5d7d186f 100644 --- a/src/user/seamount_initialization.F90 +++ b/src/user/seamount_initialization.F90 @@ -85,11 +85,13 @@ subroutine seamount_initialize_thickness ( h, G, GV, US, param_file, just_read_p logical, optional, intent(in) :: just_read_params !< If present and true, this call will !! only read parameters without changing h. - real :: e0(SZK_(G)+1) ! The resting interface heights, in m, usually ! - ! negative because it is positive upward. ! - real :: eta1D(SZK_(G)+1)! Interface height relative to the sea surface ! - ! positive upward, in m. ! - real :: min_thickness, S_surf, S_range, S_ref, S_light, S_dense + real :: e0(SZK_(G)+1) ! The resting interface heights, in depth units (Z), usually + ! negative because it is positive upward. + real :: eta1D(SZK_(G)+1)! Interface height relative to the sea surface + ! positive upward, in depth units (Z). + real :: min_thickness ! The minimum layer thicknesses, in Z. + real :: S_surf, S_range, S_ref, S_light, S_dense ! Various salinities, in ppt. + real :: eta_IC_quanta ! The granularity of quantization of intial interface heights, in Z-1. character(len=20) :: verticalCoordinate logical :: just_read ! If true, just read parameters but set nothing. integer :: i, j, k, is, ie, js, je, nz @@ -102,7 +104,7 @@ subroutine seamount_initialize_thickness ( h, G, GV, US, param_file, just_read_p call MOM_mesg("MOM_initialization.F90, initialize_thickness_uniform: setting thickness") call get_param(param_file, mdl,"MIN_THICKNESS",min_thickness, & - 'Minimum thickness for layer',& + 'Minimum thickness for layer', & units='m', default=1.0e-3, do_not_log=just_read, scale=US%m_to_Z) call get_param(param_file, mdl,"REGRIDDING_COORDINATE_MODE",verticalCoordinate, & default=DEFAULT_COORDINATE_MODE, do_not_log=just_read) @@ -126,6 +128,10 @@ subroutine seamount_initialize_thickness ( h, G, GV, US, param_file, just_read_p call get_param(param_file, mdl, "S_REF", S_ref, default=35.0, do_not_log=.true.) call get_param(param_file, mdl, "TS_RANGE_S_LIGHT", S_light, default = S_Ref, do_not_log=.true.) call get_param(param_file, mdl, "TS_RANGE_S_DENSE", S_dense, default = S_Ref, do_not_log=.true.) + call get_param(param_file, mdl, "INTERFACE_IC_QUANTA", eta_IC_quanta, & + "The granularity of initial interface height values \n"//& + "per meter, to avoid sensivity to order-of-arithmetic changes.", & + default=2048.0, units="m-1", scale=US%Z_to_m, do_not_log=just_read) if (just_read) return ! All run-time parameters have been read, so return. do K=1,nz+1 @@ -138,7 +144,8 @@ subroutine seamount_initialize_thickness ( h, G, GV, US, param_file, just_read_p e0(K) = - G%max_depth * ( ( S_light - S_surf ) + ( S_dense - S_light ) * & ( (real(K)-1.5) / real(nz-1) ) ) / S_range ! Force round numbers ... the above expression has irrational factors ... - e0(K) = nint(2048.*US%Z_to_m*e0(K))/(2048.*US%Z_to_m) + if (eta_IC_quanta > 0.0) & + e0(K) = nint(eta_IC_quanta*e0(K)) / eta_IC_quanta e0(K) = min(real(1-K)*GV%Angstrom_Z, e0(K)) ! Bound by surface e0(K) = max(-G%max_depth, e0(K)) ! Bound by bottom enddo From 9259ba1ea8b11fd3a48f04aeedb645a2b047e3c9 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 6 Nov 2018 18:59:26 -0500 Subject: [PATCH 4/6] +Added run-time parameters for sloshing test case Added two new runtime parameters, SLOSHING_IC_AMPLITUDE and SLOSHING_IC_BUG, that are used in setting up the sloshing test case. By default all answers are bitwise identical, but the MOM_parameter_doc.short files change for the sloshing test case. --- src/user/sloshing_initialization.F90 | 28 +++++++++++++++++++++------- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/src/user/sloshing_initialization.F90 b/src/user/sloshing_initialization.F90 index bda2a8d895..5c4c5ec5b6 100644 --- a/src/user/sloshing_initialization.F90 +++ b/src/user/sloshing_initialization.F90 @@ -6,7 +6,7 @@ module sloshing_initialization use MOM_domains, only : sum_across_PEs use MOM_dyn_horgrid, only : dyn_horgrid_type use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, is_root_pe -use MOM_file_parser, only : get_param, param_file_type +use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_get_input, only : directories use MOM_grid, only : ocean_grid_type use MOM_sponge, only : set_up_sponge_field, initialize_sponge, sponge_CS @@ -25,8 +25,6 @@ module sloshing_initialization public sloshing_initialize_thickness public sloshing_initialize_temperature_salinity -character(len=40) :: mdl = "sloshing_initialization" !< This module's name. - contains !> Initialization of topography. @@ -73,14 +71,27 @@ subroutine sloshing_initialize_thickness ( h, G, GV, US, param_file, just_read_p real :: weight_z ! A (misused?) depth-space weighting, in inconsistent units. real :: x1, y1, x2, y2 ! Dimensonless parameters. real :: x, t ! Dimensionless depth coordinates? + logical :: use_IC_bug ! If true, set the initial conditions retaining an old bug. logical :: just_read ! If true, just read parameters but set nothing. + ! This include declares and sets the variable "version". +# include "version_variable.h" + character(len=40) :: mdl = "sloshing_initialization" !< This module's name. integer :: i, j, k, is, ie, js, je, nx, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke just_read = .false. ; if (present(just_read_params)) just_read = just_read_params - if (just_read) return ! This subroutine has no run-time parameters. + if (.not.just_read) call log_version(param_file, mdl, version, "") + call get_param(param_file, mdl, "SLOSHING_IC_AMPLITUDE", a0, & + "Initial amplitude of sloshing internal interface height \n"//& + "displacements it the sloshing test case.", & + units='m', default=75.0, scale=US%m_to_Z, do_not_log=just_read) + call get_param(param_file, mdl, "SLOSHING_IC_BUG", use_IC_bug, & + "If true, use code with a bug to set the sloshing initial conditions.", & + default=.true., do_not_log=just_read) + + if (just_read) return ! All run-time parameters have been read, so return. ! Define thicknesses do j=G%jsc,G%jec ; do i=G%isc,G%iec @@ -116,14 +127,17 @@ subroutine sloshing_initialize_thickness ( h, G, GV, US, param_file, just_read_p enddo ! 2. Define displacement - a0 = 75.0 * US%m_to_Z ! 75m Displacement amplitude in depth units. + ! a0 is set via get_param; by default a0 is a 75m Displacement amplitude in depth units. do k = 1,nz+1 weight_z = - 4.0 * ( z_unif(k) + 0.5 )**2 + 1.0 x = G%geoLonT(i,j) / G%len_lon - !### Perhaps the '+ weight_z' here should be '* weight_z' - RWH - displ(k) = a0 * cos(acos(-1.0)*x) + weight_z * US%m_to_Z + if (use_IC_bug) then + displ(k) = a0 * cos(acos(-1.0)*x) + weight_z * US%m_to_Z + else + displ(k) = a0 * cos(acos(-1.0)*x) * weight_z + endif if ( k == 1 ) then displ(k) = 0.0 From 0ace7767073a8991b442c99462472972bd462a15 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 6 Nov 2018 18:59:44 -0500 Subject: [PATCH 5/6] Added a comment in USER_initialize_thickness Added a comment in USER_initialize_thickness and removed a constant that is being multiplied by 0. All answers are bitwise identical. --- src/user/user_initialization.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/user/user_initialization.F90 b/src/user/user_initialization.F90 index 174cd0ac8f..0a4ce7ccaa 100644 --- a/src/user/user_initialization.F90 +++ b/src/user/user_initialization.F90 @@ -92,7 +92,7 @@ subroutine USER_initialize_thickness(h, G, GV, param_file, just_read_params) if (just_read) return ! All run-time parameters have been read, so return. - h(:,:,1) = 0.0 * GV%m_to_H + h(:,:,1) = 0.0 ! h should be set in units of H. if (first_call) call write_user_log(param_file) From dfb9d9fc3c715ad88d61eab276c3f35beaff17a3 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 8 Nov 2018 18:44:44 -0500 Subject: [PATCH 6/6] Corrected an openMP directive Corrected an openMP directives that omitted a recently added unit_scale_type variable. All answers are bitwise identical, and the code once again compiles with openMP enabled. --- src/diagnostics/MOM_wave_speed.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index cd71800290..3a136bcd9b 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -135,7 +135,7 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, & rescale = 1024.0**4 ; I_rescale = 1.0/rescale min_h_frac = tol1 / real(nz) -!$OMP parallel do default(none) shared(is,ie,js,je,nz,h,G,GV,min_h_frac,use_EOS,T,S,tv,& +!$OMP parallel do default(none) shared(is,ie,js,je,nz,h,G,GV,US,min_h_frac,use_EOS,T,S,tv,& !$OMP calc_modal_structure,l_use_ebt_mode,modal_structure, & !$OMP l_mono_N2_column_fraction,l_mono_N2_depth,CS, & !$OMP Z_to_Pa,cg1,g_Rho0,rescale,I_rescale,L2_to_Z2) &