diff --git a/.testing/tools/parse_perf.py b/.testing/tools/parse_perf.py index b86b1cc106..7cbffd995d 100755 --- a/.testing/tools/parse_perf.py +++ b/.testing/tools/parse_perf.py @@ -58,9 +58,16 @@ def parse_perf_report(perf_data_path): # get per-symbol count else: - tokens = line.split() - symbol = tokens[2] - period = int(tokens[3]) + try: + tokens = line.split() + symbol = tokens[2] + period = int(tokens[3]) + except ValueError: + print("parse_perf.py: Error extracting symbol count", + file=sys.stderr) + print("line:", repr(line), file=sys.stderr) + print("tokens:", tokens, file=sys.stderr) + raise profile[event_name]['symbol'][symbol] = period diff --git a/ac/makedep b/ac/makedep index 4c9cc9229b..e0d350857e 100755 --- a/ac/makedep +++ b/ac/makedep @@ -34,6 +34,8 @@ re_procedure = re.compile( # Preprocessor expression tokenization +# NOTE: Labels and attributes could be assigned here, but for now we just use +# the token string as the label. cpp_scanner = re.Scanner([ (r'defined', lambda scanner, token: token), (r'[_A-Za-z][_0-9a-zA-Z]*', lambda scanner, token: token), @@ -56,13 +58,15 @@ cpp_scanner = re.Scanner([ (r'&', lambda scanner, token: token), (r'\|\|', lambda scanner, token: token), (r'\|', lambda scanner, token: token), - (r'^\#if', None), + (r'^ *\# *if', None), (r'\s+', None), ]) cpp_operate = { + '(': lambda x: x, '!': lambda x: not x, + 'defined': lambda x, y: x in y, '*': lambda x, y: x * y, '/': lambda x, y: x // y, '+': lambda x, y: x + y, @@ -85,6 +89,7 @@ cpp_operate = { cpp_op_rank = { '(': 13, '!': 12, + 'defined': 12, '*': 11, '/': 11, '+': 10, @@ -527,7 +532,7 @@ def cpp_expr_eval(expr, macros=None): if macros is None: macros = {} - results, remainder = cpp_scanner.scan(expr) + results, remainder = cpp_scanner.scan(expr.strip()) # Abort if any characters are not tokenized if remainder: @@ -545,72 +550,59 @@ def cpp_expr_eval(expr, macros=None): tokens = iter(results) for tok in tokens: - # Evaluate "defined()" statements - if tok == 'defined': - tok = next(tokens) - - parens = tok == '(' - if parens: - tok = next(tokens) + if tok in cpp_op_rank.keys(): + while cpp_op_rank[tok] <= cpp_op_rank[prior_op]: - # NOTE: Any key in `macros` is considered to be set, even if the - # value is None. - value = tok in macros + # Unary operators are "look ahead" so we always skip them. + # (However, `op` below could be a unary operator.) + if tok in ('!', 'defined', '('): + break - # Negation - while prior_op == '!': + second = stack.pop() op = stack.pop() - assert op == '!' - value = cpp_operate[op](value) - prior_op = stack[-1] if stack else None - - stack.append(value) - if parens: - tok = next(tokens) - assert tok == ')' + if op == '(': + value = second - elif tok.isdigit(): - value = int(tok) - stack.append(value) + elif op == '!': + if isinstance(second, str): + if second.isidentifier(): + second = macros.get(second, '0') + if second.isdigit(): + second = int(second) - elif tok.isidentifier(): - # "Identifiers that are not macros, which are all considered to be - # the number zero." (CPP manual, 4.2.2) - value = macros.get(tok, '0') - if value.isdigit(): - value = int(value) - stack.append(value) + value = cpp_operate[op](second) - elif tok in cpp_op_rank.keys(): - while cpp_op_rank[tok] <= cpp_op_rank[prior_op]: + elif op == 'defined': + value = cpp_operate[op](second, macros) - # Skip unary prefix operators (only '!' at the moment) - if tok == '!': - break + else: + first = stack.pop() - second = stack.pop() - op = stack.pop() - first = stack.pop() + if isinstance(first, str): + if first.isidentifier(): + first = macros.get(first, '0') + if first.isdigit(): + first = int(first) - value = cpp_operate[op](first, second) - prior_op = stack[-1] if stack else None + if isinstance(second, str): + if second.isidentifier(): + second = macros.get(second, '0') + if second.isdigit(): + second = int(second) - if prior_op == '(': - prior_op = None - if tok == ')': - stack.pop() + value = cpp_operate[op](first, second) + prior_op = stack[-1] if stack else None stack.append(value) - if tok == ')': - prior_op = stack[-2] if stack and len(stack) > 1 else None - else: + # The ) "operator" has already been applied, so it can be dropped. + if tok != ')': stack.append(tok) prior_op = tok - if prior_op in ('(',): - prior_op = None + elif tok.isdigit() or tok.isidentifier(): + stack.append(tok) else: print("Unsupported token:", tok) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 543d77a0f3..a083402fde 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -40,7 +40,7 @@ module MOM_ALE use MOM_remapping, only : remapping_core_h, remapping_core_w use MOM_remapping, only : remappingSchemesDoc, remappingDefaultScheme use MOM_remapping, only : interpolate_column, reintegrate_column -use MOM_remapping, only : remapping_CS, dzFromH1H2 +use MOM_remapping, only : remapping_CS, dzFromH1H2, remapping_set_param use MOM_string_functions, only : uppercase, extractWord, extract_integer use MOM_tracer_registry, only : tracer_registry_type, tracer_type, MOM_tracer_chkinv use MOM_unit_scaling, only : unit_scale_type @@ -147,6 +147,7 @@ module MOM_ALE public pre_ALE_adjustments public ALE_remap_init_conds public ALE_register_diags +public ALE_set_extrap_boundaries ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with @@ -176,6 +177,7 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS) logical :: force_bounds_in_subcell logical :: local_logical logical :: remap_boundary_extrap + logical :: init_boundary_extrap type(hybgen_regrid_CS), pointer :: hybgen_regridCS => NULL() ! Control structure for hybgen regridding ! for sharing parameters. @@ -225,6 +227,10 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS) call get_param(param_file, mdl, "REMAP_BOUNDARY_EXTRAP", remap_boundary_extrap, & "If true, values at the interfaces of boundary cells are "//& "extrapolated instead of piecewise constant", default=.false.) + call get_param(param_file, mdl, "INIT_BOUNDARY_EXTRAP", init_boundary_extrap, & + "If true, values at the interfaces of boundary cells are "//& + "extrapolated instead of piecewise constant during initialization."//& + "Defaults to REMAP_BOUNDARY_EXTRAP.", default=remap_boundary_extrap) 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) @@ -237,13 +243,13 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS) if (.not.GV%Boussinesq) CS%answer_date = max(CS%answer_date, 20230701) call initialize_remapping( CS%remapCS, string, & - boundary_extrapolation=remap_boundary_extrap, & + boundary_extrapolation=init_boundary_extrap, & check_reconstruction=check_reconstruction, & check_remapping=check_remapping, & force_bounds_in_subcell=force_bounds_in_subcell, & answer_date=CS%answer_date) call initialize_remapping( CS%vel_remapCS, vel_string, & - boundary_extrapolation=remap_boundary_extrap, & + boundary_extrapolation=init_boundary_extrap, & check_reconstruction=check_reconstruction, & check_remapping=check_remapping, & force_bounds_in_subcell=force_bounds_in_subcell, & @@ -308,6 +314,18 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS) if (CS%show_call_tree) call callTree_leave("ALE_init()") end subroutine ALE_init +!> Sets the boundary extrapolation set for the remapping type. +subroutine ALE_set_extrap_boundaries( param_file, CS) + type(param_file_type), intent(in) :: param_file !< Parameter file + type(ALE_CS), pointer :: CS !< Module control structure + + logical :: remap_boundary_extrap + call get_param(param_file, "MOM_ALE", "REMAP_BOUNDARY_EXTRAP", remap_boundary_extrap, & + "If true, values at the interfaces of boundary cells are "//& + "extrapolated instead of piecewise constant", default=.false.) + call remapping_set_param(CS%remapCS, boundary_extrapolation=remap_boundary_extrap) +end subroutine ALE_set_extrap_boundaries + !> Initialize diagnostics for the ALE module. subroutine ALE_register_diags(Time, G, GV, US, diag, CS) type(time_type),target, intent(in) :: Time !< Time structure diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 3521131431..0b602be944 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -56,6 +56,7 @@ module MOM use MOM_ALE, only : ALE_remap_tracers, ALE_remap_velocities use MOM_ALE, only : ALE_remap_set_h_vel, ALE_remap_set_h_vel_via_dz use MOM_ALE, only : ALE_update_regrid_weights, pre_ALE_diagnostics, ALE_register_diags +use MOM_ALE, only : ALE_set_extrap_boundaries use MOM_ALE_sponge, only : rotate_ALE_sponge, update_ALE_sponge_field use MOM_barotropic, only : Barotropic_CS use MOM_boundary_update, only : call_OBC_register, OBC_register_end, update_OBC_CS @@ -3120,8 +3121,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & endif endif endif - if ( CS%use_ALE_algorithm ) call ALE_updateVerticalGridType( CS%ALE_CSp, GV ) - + if ( CS%use_ALE_algorithm ) then + call ALE_set_extrap_boundaries (param_file, CS%ALE_CSp) + call callTree_waypoint("returned from ALE_init() (initialize_MOM)") + call ALE_updateVerticalGridType( CS%ALE_CSp, GV ) + endif ! The basic state variables have now been fully initialized, so update their halos and ! calculate any derived thermodynmics quantities. diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90 index 056b171ba8..00a289ab9a 100644 --- a/src/core/MOM_CoriolisAdv.F90 +++ b/src/core/MOM_CoriolisAdv.F90 @@ -902,20 +902,20 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS, pbv, Wav if (associated(AD%rv_x_u)) then do J=Jsq,Jeq ; do i=is,ie AD%rv_x_u(i,J,k) = -G%IdyCv(i,J) * C1_12 * & - ((q2(I,J) + q2(I-1,J) + q2(I-1,J-1)) * uh(I-1,j,k) + & - (q2(I-1,J) + q2(I,J) + q2(I,J-1)) * uh(I,j,k) + & - (q2(I-1,J) + q2(I,J+1) + q2(I,J)) * uh(I,j+1,k) + & - (q2(I,J) + q2(I-1,J+1) + q2(I-1,J)) * uh(I-1,j+1,k)) + (((((q2(I,J) + q2(I-1,J-1)) + q2(I-1,J)) * uh(I-1,j,k)) + & + (((q2(I-1,J) + q2(I,J+1)) + q2(I,J)) * uh(I,j+1,k))) + & + ((((q2(I-1,J) + q2(I,J-1)) + q2(I,J)) * uh(I,j,k))+ & + (((q2(I,J) + q2(I-1,J+1)) + q2(I-1,J)) * uh(I-1,j+1,k)))) enddo ; enddo endif if (associated(AD%rv_x_v)) then do j=js,je ; do I=Isq,Ieq AD%rv_x_v(I,j,k) = G%IdxCu(I,j) * C1_12 * & - ((q2(I+1,J) + q2(I,J) + q2(I,J-1)) * vh(i+1,J,k) + & - (q2(I-1,J) + q2(I,J) + q2(I,J-1)) * vh(i,J,k) + & - (q2(I-1,J-1) + q2(I,J) + q2(I,J-1)) * vh(i,J-1,k) + & - (q2(I+1,J-1) + q2(I,J) + q2(I,J-1)) * vh(i+1,J-1,k)) + (((((q2(I+1,J) + q2(I,J-1)) + q2(I,J)) * vh(i+1,J,k)) + & + (((q2(I-1,J-1) + q2(I,J)) + q2(I,J-1)) * vh(i,J-1,k))) + & + ((((q2(I-1,J) + q2(I,J-1)) + q2(I,J)) * vh(i,J,k)) + & + (((q2(I+1,J-1) + q2(I,J)) + q2(I,J-1)) * vh(i+1,J-1,k)))) enddo ; enddo endif endif diff --git a/src/core/MOM_density_integrals.F90 b/src/core/MOM_density_integrals.F90 index 9fed528e71..b73524e362 100644 --- a/src/core/MOM_density_integrals.F90 +++ b/src/core/MOM_density_integrals.F90 @@ -141,14 +141,21 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m] ! Local variables - real :: T5(5), S5(5) ! Temperatures and salinities at five quadrature points [C ~> degC] and [S ~> ppt] - real :: p5(5) ! Pressures at five quadrature points [R L2 T-2 ~> Pa] - real :: r5(5) ! Densities at five quadrature points [R ~> kg m-3] + real :: T5((5*HI%iscB+1):(5*(HI%iecB+2))) ! Temperatures along a line of subgrid locations [C ~> degC] + real :: S5((5*HI%iscB+1):(5*(HI%iecB+2))) ! Salinities along a line of subgrid locations [S ~> ppt] + real :: p5((5*HI%iscB+1):(5*(HI%iecB+2))) ! Pressures along a line of subgrid locations [R L2 T-2 ~> Pa] + real :: r5((5*HI%iscB+1):(5*(HI%iecB+2))) ! Densities anomalies along a line of subgrid locations [R ~> kg m-3] + real :: T15((15*HI%iscB+1):(15*(HI%iecB+1))) ! Temperatures at an array of subgrid locations [C ~> degC] + real :: S15((15*HI%iscB+1):(15*(HI%iecB+1))) ! Salinities at an array of subgrid locations [S ~> ppt] + real :: p15((15*HI%iscB+1):(15*(HI%iecB+1))) ! Pressures at an array of subgrid locations [R L2 T-2 ~> Pa] + real :: r15((15*HI%iscB+1):(15*(HI%iecB+1))) ! Densities at an array of subgrid locations [R ~> kg m-3] real :: rho_anom ! The depth averaged density anomaly [R ~> kg m-3] real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim] real :: GxRho ! The product of the gravitational acceleration and reference density [R L2 Z-1 T-2 ~> Pa m-1] real :: I_Rho ! The inverse of the Boussinesq density [R-1 ~> m3 kg-1] real :: dz ! The layer thickness [Z ~> m] + real :: dz_x(5,HI%iscB:HI%iecB) ! Layer thicknesses along an x-line of subgrid locations [Z ~> m] + real :: dz_y(5,HI%isc:HI%iec) ! Layer thicknesses along a y-line of subgrid locations [Z ~> m] real :: z0pres ! The height at which the pressure is zero [Z ~> m] real :: hWght ! A pressure-thickness below topography [Z ~> m] real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [Z ~> m] @@ -162,7 +169,10 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & logical :: do_massWeight ! Indicates whether to do mass weighting. logical :: use_rho_ref ! Pass rho_ref to the equation of state for more accurate calculation ! of density anomalies. - integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, i, j, m, n + integer, dimension(2) :: EOSdom_h5 ! The 5-point h-point i-computational domain for the equation of state + integer, dimension(2) :: EOSdom_q15 ! The 3x5-point q-point i-computational domain for the equation of state + integer, dimension(2) :: EOSdom_h15 ! The 3x5-point h-point i-computational domain for the equation of state + integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, i, j, m, n, pos ! These array bounds work for the indexing convention of the input arrays, but ! on the computational domain defined for the output arrays. @@ -188,123 +198,169 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & "dz_neglect must be present if useMassWghtInterp is present and true.") endif ; endif - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - dz = z_t(i,j) - z_b(i,j) - do n=1,5 - T5(n) = T(i,j) ; S5(n) = S(i,j) - p5(n) = -GxRho*((z_t(i,j) - z0pres) - 0.25*real(n-1)*dz) + ! Set the loop ranges for equation of state calculations at various points. + EOSdom_h5(1) = 1 ; EOSdom_h5(2) = 5*(Ieq-Isq+2) + EOSdom_q15(1) = 1 ; EOSdom_q15(2) = 15*(Ieq-Isq+1) + EOSdom_h15(1) = 1 ; EOSdom_h15(2) = 15*(HI%iec-HI%isc+1) + + do j=Jsq,Jeq+1 + do i=Isq,Ieq+1 + dz = z_t(i,j) - z_b(i,j) + do n=1,5 + T5(i*5+n) = T(i,j) ; S5(i*5+n) = S(i,j) + p5(i*5+n) = -GxRho*((z_t(i,j) - z0pres) - 0.25*real(n-1)*dz) + enddo enddo + if (use_rho_ref) then - call calculate_density(T5, S5, p5, r5, EOS, rho_ref=rho_ref) - ! Use Boole's rule to estimate the pressure anomaly change. - rho_anom = C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) + call calculate_density(T5, S5, p5, r5, EOS, EOSdom_h5, rho_ref=rho_ref) else - call calculate_density(T5, S5, p5, r5, EOS) - ! Use Boole's rule to estimate the pressure anomaly change. - rho_anom = C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) - rho_ref + call calculate_density(T5, S5, p5, r5, EOS, EOSdom_h5) endif - dpa(i,j) = G_e*dz*rho_anom - ! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of - ! the pressure anomaly. - if (present(intz_dpa)) intz_dpa(i,j) = 0.5*G_e*dz**2 * & - (rho_anom - C1_90*(16.0*(r5(4)-r5(2)) + 7.0*(r5(5)-r5(1))) ) - enddo ; enddo + do i=Isq,Ieq+1 + ! Use Boole's rule to estimate the pressure anomaly change. + rho_anom = C1_90*(7.0*(r5(i*5+1)+r5(i*5+5)) + 32.0*(r5(i*5+2)+r5(i*5+4)) + 12.0*r5(i*5+3)) + if (.not.use_rho_ref) rho_anom = rho_anom - rho_ref + dz = z_t(i,j) - z_b(i,j) + dpa(i,j) = G_e*dz*rho_anom + ! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of + ! the pressure anomaly. + if (present(intz_dpa)) intz_dpa(i,j) = 0.5*G_e*dz**2 * & + (rho_anom - C1_90*(16.0*(r5(i*5+4)-r5(i*5+2)) + 7.0*(r5(i*5+5)-r5(i*5+1))) ) + enddo + enddo - if (present(intx_dpa)) then ; do j=js,je ; do I=Isq,Ieq - ! hWght is the distance measure by which the cell is violation of - ! hydrostatic consistency. For large hWght we bias the interpolation of - ! T & S along the top and bottom integrals, akin to thickness weighting. - hWght = 0.0 - if (do_massWeight) & - hWght = max(0., -bathyT(i,j)-z_t(i+1,j), -bathyT(i+1,j)-z_t(i,j)) - if (hWght > 0.) then - hL = (z_t(i,j) - z_b(i,j)) + dz_neglect - hR = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect - hWght = hWght * ( (hL-hR)/(hL+hR) )**2 - iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR ) - hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom - hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom + if (present(intx_dpa)) then ; do j=js,je + do I=Isq,Ieq + ! hWght is the distance measure by which the cell is violation of + ! hydrostatic consistency. For large hWght we bias the interpolation of + ! T & S along the top and bottom integrals, akin to thickness weighting. + hWght = 0.0 + if (do_massWeight) & + hWght = max(0., -bathyT(i,j)-z_t(i+1,j), -bathyT(i+1,j)-z_t(i,j)) + if (hWght > 0.) then + hL = (z_t(i,j) - z_b(i,j)) + dz_neglect + hR = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect + hWght = hWght * ( (hL-hR)/(hL+hR) )**2 + iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR ) + hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom + hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom + else + hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0 + endif + + do m=2,4 + ! T, S, and z are interpolated in the horizontal. The z interpolation + ! is linear, but for T and S it may be thickness weighted. + wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L + wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR + dz_x(m,i) = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i+1,j) - z_b(i+1,j)) + pos = i*15+(m-2)*5 + T15(pos+1) = wtT_L*T(i,j) + wtT_R*T(i+1,j) + S15(pos+1) = wtT_L*S(i,j) + wtT_R*S(i+1,j) + p15(pos+1) = -GxRho*((wt_L*z_t(i,j) + wt_R*z_t(i+1,j)) - z0pres) + do n=2,5 + T15(pos+n) = T15(pos+1) ; S15(pos+n) = S15(pos+1) + p15(pos+n) = p15(pos+n-1) + GxRho*0.25*dz_x(m,i) + enddo + enddo + enddo + + if (use_rho_ref) then + call calculate_density(T15, S15, p15, r15, EOS, EOSdom_q15, rho_ref=rho_ref) else - hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0 + call calculate_density(T15, S15, p15, r15, EOS, EOSdom_q15) endif - intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j) - do m=2,4 - ! T, S, and z are interpolated in the horizontal. The z interpolation - ! is linear, but for T and S it may be thickness weighted. - wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L - wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR - dz = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i+1,j) - z_b(i+1,j)) - T5(1) = wtT_L*T(i,j) + wtT_R*T(i+1,j) - S5(1) = wtT_L*S(i,j) + wtT_R*S(i+1,j) - p5(1) = -GxRho*((wt_L*z_t(i,j) + wt_R*z_t(i+1,j)) - z0pres) - do n=2,5 - T5(n) = T5(1) ; S5(n) = S5(1) ; p5(n) = p5(n-1) + GxRho*0.25*dz - enddo + do I=Isq,Ieq + intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j) + ! Use Boole's rule to estimate the pressure anomaly change. if (use_rho_ref) then - call calculate_density(T5, S5, p5, r5, EOS, rho_ref=rho_ref) - ! Use Boole's rule to estimate the pressure anomaly change. - intz(m) = G_e*dz*( C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3))) + do m=2,4 + pos = i*15+(m-2)*5 + intz(m) = G_e*dz_x(m,i)*( C1_90*( 7.0*(r15(pos+1)+r15(pos+5)) + & + 32.0*(r15(pos+2)+r15(pos+4)) + & + 12.0*r15(pos+3))) + enddo else - call calculate_density(T5, S5, p5, r5, EOS) - ! Use Boole's rule to estimate the pressure anomaly change. - intz(m) = G_e*dz*( C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) - rho_ref ) + do m=2,4 + pos = i*15+(m-2)*5 + intz(m) = G_e*dz_x(m,i)*( C1_90*( 7.0*(r15(pos+1)+r15(pos+5)) + & + 32.0*(r15(pos+2)+r15(pos+4)) + & + 12.0*r15(pos+3)) - rho_ref ) + enddo endif + ! Use Boole's rule to integrate the bottom pressure anomaly values in x. + intx_dpa(i,j) = C1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + & + 12.0*intz(3)) + enddo + enddo ; endif + if (present(inty_dpa)) then ; do J=Jsq,Jeq + do i=is,ie + ! hWght is the distance measure by which the cell is violation of + ! hydrostatic consistency. For large hWght we bias the interpolation of + ! T & S along the top and bottom integrals, akin to thickness weighting. + hWght = 0.0 + if (do_massWeight) & + hWght = max(0., -bathyT(i,j)-z_t(i,j+1), -bathyT(i,j+1)-z_t(i,j)) + if (hWght > 0.) then + hL = (z_t(i,j) - z_b(i,j)) + dz_neglect + hR = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect + hWght = hWght * ( (hL-hR)/(hL+hR) )**2 + iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR ) + hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom + hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom + else + hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0 + endif + + do m=2,4 + ! T, S, and z are interpolated in the horizontal. The z interpolation + ! is linear, but for T and S it may be thickness weighted. + wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L + wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR + dz_y(m,i) = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i,j+1) - z_b(i,j+1)) + pos = i*15+(m-2)*5 + T15(pos+1) = wtT_L*T(i,j) + wtT_R*T(i,j+1) + S15(pos+1) = wtT_L*S(i,j) + wtT_R*S(i,j+1) + p15(pos+1) = -GxRho*((wt_L*z_t(i,j) + wt_R*z_t(i,j+1)) - z0pres) + do n=2,5 + T15(pos+n) = T15(pos+1) ; S15(pos+n) = S15(pos+1) + p15(pos+n) = p15(pos+n-1) + GxRho*0.25*dz_y(m,i) + enddo + enddo enddo - ! Use Boole's rule to integrate the bottom pressure anomaly values in x. - intx_dpa(i,j) = C1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + & - 12.0*intz(3)) - enddo ; enddo ; endif - if (present(inty_dpa)) then ; do J=Jsq,Jeq ; do i=is,ie - ! hWght is the distance measure by which the cell is violation of - ! hydrostatic consistency. For large hWght we bias the interpolation of - ! T & S along the top and bottom integrals, akin to thickness weighting. - hWght = 0.0 - if (do_massWeight) & - hWght = max(0., -bathyT(i,j)-z_t(i,j+1), -bathyT(i,j+1)-z_t(i,j)) - if (hWght > 0.) then - hL = (z_t(i,j) - z_b(i,j)) + dz_neglect - hR = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect - hWght = hWght * ( (hL-hR)/(hL+hR) )**2 - iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR ) - hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom - hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom + if (use_rho_ref) then + call calculate_density(T15(15*HI%isc+1:), S15(15*HI%isc+1:), p15(15*HI%isc+1:), & + r15(15*HI%isc+1:), EOS, EOSdom_h15, rho_ref=rho_ref) else - hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0 + call calculate_density(T15(15*HI%isc+1:), S15(15*HI%isc+1:), p15(15*HI%isc+1:), & + r15(15*HI%isc+1:), EOS, EOSdom_h15) endif - intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1) - do m=2,4 - ! T, S, and z are interpolated in the horizontal. The z interpolation - ! is linear, but for T and S it may be thickness weighted. - wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L - wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR - dz = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i,j+1) - z_b(i,j+1)) - T5(1) = wtT_L*T(i,j) + wtT_R*T(i,j+1) - S5(1) = wtT_L*S(i,j) + wtT_R*S(i,j+1) - p5(1) = -GxRho*((wt_L*z_t(i,j) + wt_R*z_t(i,j+1)) - z0pres) - do n=2,5 - T5(n) = T5(1) ; S5(n) = S5(1) - p5(n) = p5(n-1) + GxRho*0.25*dz + do i=is,ie + intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1) + ! Use Boole's rule to estimate the pressure anomaly change. + do m=2,4 + pos = i*15+(m-2)*5 + if (use_rho_ref) then + intz(m) = G_e*dz_y(m,i)*( C1_90*(7.0*(r15(pos+1)+r15(pos+5)) + & + 32.0*(r15(pos+2)+r15(pos+4)) + & + 12.0*r15(pos+3))) + else + intz(m) = G_e*dz_y(m,i)*( C1_90*(7.0*(r15(pos+1)+r15(pos+5)) + & + 32.0*(r15(pos+2)+r15(pos+4)) + & + 12.0*r15(pos+3)) - rho_ref ) + endif enddo - if (use_rho_ref) then - call calculate_density(T5, S5, p5, r5, EOS, rho_ref=rho_ref) - ! Use Boole's rule to estimate the pressure anomaly change. - intz(m) = G_e*dz*( C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3))) - else - call calculate_density(T5, S5, p5, r5, EOS) - ! Use Boole's rule to estimate the pressure anomaly change. - intz(m) = G_e*dz*( C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) - rho_ref ) - endif - + ! Use Boole's rule to integrate the values. + inty_dpa(i,j) = C1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + & + 12.0*intz(3)) enddo - ! Use Boole's rule to integrate the values. - inty_dpa(i,j) = C1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + & - 12.0*intz(3)) - enddo ; enddo ; endif + enddo ; endif end subroutine int_density_dz_generic_pcm @@ -414,10 +470,9 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & logical :: use_rho_ref ! Pass rho_ref to the equation of state for more accurate calculation ! of density anomalies. logical :: use_varT, use_varS, use_covarTS ! Logicals for SGS variances fields - integer, dimension(2) :: EOSdom_q5 ! The 5-point q-point i-computational domain for the equation of state + integer, dimension(2) :: EOSdom_h5 ! The 5-point h-point i-computational domain for the equation of state integer, dimension(2) :: EOSdom_q15 ! The 3x5-point q-point i-computational domain for the equation of state integer, dimension(2) :: EOSdom_h15 ! The 3x5-point h-point i-computational domain for the equation of state - integer :: Isq, Ieq, Jsq, Jeq, i, j, m, n, pos Isq = HI%IscB ; Ieq = HI%IecB ; Jsq = HI%JscB ; Jeq = HI%JecB @@ -456,8 +511,8 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & enddo ! Set the loop ranges for equation of state calculations at various points. - EOSdom_q5(1) = 1 ; EOSdom_q5(2) = (ieq-isq+2)*5 - EOSdom_q15(1) = 1 ; EOSdom_q15(2) = 15*(ieq-isq+1) + EOSdom_h5(1) = 1 ; EOSdom_h5(2) = 5*(Ieq-Isq+2) + EOSdom_q15(1) = 1 ; EOSdom_q15(2) = 15*(Ieq-Isq+1) EOSdom_h15(1) = 1 ; EOSdom_h15(2) = 15*(HI%iec-HI%isc+1) ! 1. Compute vertical integrals @@ -475,12 +530,12 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & if (use_varS) S25(i*5+1:i*5+5) = tv%varS(i,j,k) enddo if (use_Stanley_eos) then - call calculate_density(T5, S5, p5, T25, TS5, S25, r5, EOS, EOSdom_q5, rho_ref=rho_ref) + call calculate_density(T5, S5, p5, T25, TS5, S25, r5, EOS, EOSdom_h5, rho_ref=rho_ref) else if (use_rho_ref) then - call calculate_density(T5, S5, p5, r5, EOS, EOSdom_q5, rho_ref=rho_ref) + call calculate_density(T5, S5, p5, r5, EOS, EOSdom_h5, rho_ref=rho_ref) else - call calculate_density(T5, S5, p5, r5, EOS, EOSdom_q5) + call calculate_density(T5, S5, p5, r5, EOS, EOSdom_h5) u5(:) = r5(:) - rho_ref endif endif @@ -491,8 +546,8 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & rho_anom = C1_90*(7.0*(r5(i*5+1)+r5(i*5+5)) + 32.0*(r5(i*5+2)+r5(i*5+4)) + 12.0*r5(i*5+3)) dpa(i,j) = G_e*dz(i)*rho_anom if (present(intz_dpa)) then - ! Use a Boole's-rule-like fifth-order accurate estimate of - ! the double integral of the pressure anomaly. + ! Use a Boole's-rule-like fifth-order accurate estimate of + ! the double integral of the pressure anomaly. intz_dpa(i,j) = 0.5*G_e*dz(i)**2 * & (rho_anom - C1_90*(16.0*(r5(i*5+4)-r5(i*5+2)) + 7.0*(r5(i*5+5)-r5(i*5+1))) ) endif @@ -504,8 +559,8 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & - rho_ref dpa(i,j) = G_e*dz(i)*rho_anom if (present(intz_dpa)) then - ! Use a Boole's-rule-like fifth-order accurate estimate of - ! the double integral of the pressure anomaly. + ! Use a Boole's-rule-like fifth-order accurate estimate of + ! the double integral of the pressure anomaly. intz_dpa(i,j) = 0.5*G_e*dz(i)**2 * & (rho_anom - C1_90*(16.0*(u5(i*5+4)-u5(i*5+2)) + 7.0*(u5(i*5+5)-u5(i*5+1))) ) endif @@ -774,13 +829,26 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & ! a parabolic interpolation is used to compute intermediate values. ! Local variables - real :: T5(5) ! Temperatures along a line of subgrid locations [C ~> degC] - real :: S5(5) ! Salinities along a line of subgrid locations [S ~> ppt] - real :: T25(5) ! SGS temperature variance along a line of subgrid locations [C2 ~> degC2] - real :: TS5(5) ! SGS temperature-salinity covariance along a line of subgrid locations [C S ~> degC ppt] - real :: S25(5) ! SGS salinity variance along a line of subgrid locations [S2 ~> ppt2] - real :: p5(5) ! Pressures at five quadrature points [R L2 T-2 ~> Pa] - real :: r5(5) ! Density anomalies from rho_ref at quadrature points [R ~> kg m-3] + real :: T5((5*HI%iscB+1):(5*(HI%iecB+2))) ! Temperatures along a line of subgrid locations [C ~> degC] + real :: S5((5*HI%iscB+1):(5*(HI%iecB+2))) ! Salinities along a line of subgrid locations [S ~> ppt] + real :: T25((5*HI%iscB+1):(5*(HI%iecB+2))) ! SGS temperature variance along a line of subgrid + ! locations [C2 ~> degC2] + real :: TS5((5*HI%iscB+1):(5*(HI%iecB+2))) ! SGS temp-salt covariance along a line of subgrid + ! locations [C S ~> degC ppt] + real :: S25((5*HI%iscB+1):(5*(HI%iecB+2))) ! SGS salinity variance along a line of subgrid locations [S2 ~> ppt2] + real :: p5((5*HI%iscB+1):(5*(HI%iecB+2))) ! Pressures along a line of subgrid locations [R L2 T-2 ~> Pa] + real :: r5((5*HI%iscB+1):(5*(HI%iecB+2))) ! Densities anomalies along a line of subgrid + ! locations [R ~> kg m-3] + real :: T15((15*HI%iscB+1):(15*(HI%iecB+1))) ! Temperatures at an array of subgrid locations [C ~> degC] + real :: S15((15*HI%iscB+1):(15*(HI%iecB+1))) ! Salinities at an array of subgrid locations [S ~> ppt] + real :: T215((15*HI%iscB+1):(15*(HI%iecB+1))) ! SGS temperature variance along a line of subgrid + ! locations [C2 ~> degC2] + real :: TS15((15*HI%iscB+1):(15*(HI%iecB+1))) ! SGS temp-salt covariance along a line of subgrid + ! locations [C S ~> degC ppt] + real :: S215((15*HI%iscB+1):(15*(HI%iecB+1))) ! SGS salinity variance along a line of subgrid + ! locations [S2 ~> ppt2] + real :: p15((15*HI%iscB+1):(15*(HI%iecB+1))) ! Pressures at an array of subgrid locations [R L2 T-2 ~> Pa] + real :: r15((15*HI%iscB+1):(15*(HI%iecB+1))) ! Densities at an array of subgrid locations [R ~> kg m-3] real :: wt_t(5), wt_b(5) ! Top and bottom weights [nondim] real :: rho_anom ! The integrated density anomaly [R ~> kg m-3] real :: w_left, w_right ! Left and right weights [nondim] @@ -790,6 +858,8 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & real :: GxRho ! The gravitational acceleration times density [R L2 Z-1 T-2 ~> kg m-2 s-2] real :: I_Rho ! The inverse of the Boussinesq density [R-1 ~> m3 kg-1] real :: dz ! Layer thicknesses at tracer points [Z ~> m] + real :: dz_x(5,HI%iscB:HI%iecB) ! Layer thicknesses along an x-line of subgrid locations [Z ~> m] + real :: dz_y(5,HI%isc:HI%iec) ! Layer thicknesses along a y-line of subgrid locations [Z ~> m] real :: massWeightToggle ! A non-dimensional toggle factor (0 or 1) [nondim] real :: Ttl, Tbl, Tml, Ttr, Tbr, Tmr ! Temperatures at the velocity cell corners [C ~> degC] real :: Stl, Sbl, Sml, Str, Sbr, Smr ! Salinities at the velocity cell corners [S ~> ppt] @@ -801,9 +871,12 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & real :: hWght ! A topographically limited thickness weight [Z ~> m] real :: hL, hR ! Thicknesses to the left and right [Z ~> m] real :: iDenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2] - integer :: Isq, Ieq, Jsq, Jeq, i, j, m, n + integer, dimension(2) :: EOSdom_h5 ! The 5-point h-point i-computational domain for the equation of state + integer, dimension(2) :: EOSdom_q15 ! The 3x5-point q-point i-computational domain for the equation of state + integer, dimension(2) :: EOSdom_h15 ! The 3x5-point h-point i-computational domain for the equation of state + integer :: Isq, Ieq, Jsq, Jeq, i, j, m, n, pos logical :: use_PPM ! If false, assume zero curvature in reconstruction, i.e. PLM - logical :: use_varT, use_varS, use_covarTS + logical :: use_varT, use_varS, use_covarTS ! Logicals for SGS variances fields Isq = HI%IscB ; Ieq = HI%IecB ; Jsq = HI%JscB ; Jeq = HI%JecB @@ -824,226 +897,277 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & use_covarTS = .false. use_varS = .false. if (use_stanley_eos) then - use_varT = associated(tv%varT) - use_covarTS = associated(tv%covarTS) - use_varS = associated(tv%varS) + use_varT = associated(tv%varT) + use_covarTS = associated(tv%covarTS) + use_varS = associated(tv%varS) endif T25(:) = 0. TS5(:) = 0. S25(:) = 0. + T215(:) = 0. + TS15(:) = 0. + S215(:) = 0. do n = 1, 5 wt_t(n) = 0.25 * real(5-n) wt_b(n) = 1.0 - wt_t(n) enddo + ! Set the loop ranges for equation of state calculations at various points. + EOSdom_h5(1) = 1 ; EOSdom_h5(2) = 5*(Ieq-Isq+2) + EOSdom_q15(1) = 1 ; EOSdom_q15(2) = 15*(Ieq-Isq+1) + EOSdom_h15(1) = 1 ; EOSdom_h15(2) = 15*(HI%iec-HI%isc+1) + ! 1. Compute vertical integrals - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - if (use_PPM) then - ! Curvature coefficient of the parabolas - s6 = 3.0 * ( 2.0*tv%S(i,j,k) - ( S_t(i,j,k) + S_b(i,j,k) ) ) - t6 = 3.0 * ( 2.0*tv%T(i,j,k) - ( T_t(i,j,k) + T_b(i,j,k) ) ) - endif - dz = e(i,j,K) - e(i,j,K+1) - do n=1,5 - p5(n) = -GxRho*((e(i,j,K) - z0pres) - 0.25*real(n-1)*dz) - ! Salinity and temperature points are reconstructed with PPM - S5(n) = wt_t(n) * S_t(i,j,k) + wt_b(n) * ( S_b(i,j,k) + s6 * wt_t(n) ) - T5(n) = wt_t(n) * T_t(i,j,k) + wt_b(n) * ( T_b(i,j,k) + t6 * wt_t(n) ) + do j=Jsq,Jeq+1 + do i=Isq,Ieq+1 + if (use_PPM) then + ! Curvature coefficient of the parabolas + s6 = 3.0 * ( 2.0*tv%S(i,j,k) - ( S_t(i,j,k) + S_b(i,j,k) ) ) + t6 = 3.0 * ( 2.0*tv%T(i,j,k) - ( T_t(i,j,k) + T_b(i,j,k) ) ) + endif + dz = e(i,j,K) - e(i,j,K+1) + do n=1,5 + p5(I*5+n) = -GxRho*((e(i,j,K) - z0pres) - 0.25*real(n-1)*dz) + ! Salinity and temperature points are reconstructed with PPM + S5(I*5+n) = wt_t(n) * S_t(i,j,k) + wt_b(n) * ( S_b(i,j,k) + s6 * wt_t(n) ) + T5(I*5+n) = wt_t(n) * T_t(i,j,k) + wt_b(n) * ( T_b(i,j,k) + t6 * wt_t(n) ) + enddo + if (use_stanley_eos) then + if (use_varT) T25(I*5+1:I*5+5) = tv%varT(i,j,k) + if (use_covarTS) TS5(I*5+1:I*5+5) = tv%covarTS(i,j,k) + if (use_varS) S25(I*5+1:I*5+5) = tv%varS(i,j,k) + endif enddo + if (use_stanley_eos) then - if (use_varT) T25(:) = tv%varT(i,j,k) - if (use_covarTS) TS5(:) = tv%covarTS(i,j,k) - if (use_varS) S25(:) = tv%varS(i,j,k) - call calculate_density(T5, S5, p5, T25, TS5, S25, r5, EOS, rho_ref=rho_ref) + call calculate_density(T5, S5, p5, T25, TS5, S25, r5, EOS, EOSdom_h5, rho_ref=rho_ref) else - call calculate_density(T5, S5, p5, r5, EOS, rho_ref=rho_ref) + call calculate_density(T5, S5, p5, r5, EOS, EOSdom_h5, rho_ref=rho_ref) endif - ! Use Boole's rule to estimate the pressure anomaly change. - rho_anom = C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) - dpa(i,j) = G_e*dz*rho_anom - if (present(intz_dpa)) then - ! Use a Boole's-rule-like fifth-order accurate estimate of - ! the double integral of the pressure anomaly. - intz_dpa(i,j) = 0.5*G_e*dz**2 * & - (rho_anom - C1_90*(16.0*(r5(4)-r5(2)) + 7.0*(r5(5)-r5(1))) ) - endif - enddo ; enddo ! end loops on j and i + do i=Isq,Ieq+1 + dz = e(i,j,K) - e(i,j,K+1) + ! Use Boole's rule to estimate the pressure anomaly change. + rho_anom = C1_90*(7.0*(r5(i*5+1)+r5(i*5+5)) + 32.0*(r5(i*5+2)+r5(i*5+4)) + 12.0*r5(i*5+3)) + dpa(i,j) = G_e*dz*rho_anom + if (present(intz_dpa)) then + ! Use a Boole's-rule-like fifth-order accurate estimate of + ! the double integral of the pressure anomaly. + intz_dpa(i,j) = 0.5*G_e*dz**2 * & + (rho_anom - C1_90*(16.0*(r5(i*5+4)-r5(i*5+2)) + 7.0*(r5(i*5+5)-r5(i*5+1))) ) + endif + enddo ! end loop on i + enddo ! end loop on j ! 2. Compute horizontal integrals in the x direction - if (present(intx_dpa)) then ; do j=HI%jsc,HI%jec ; do I=Isq,Ieq - ! Corner values of T and S - ! hWght is the distance measure by which the cell is violation of - ! hydrostatic consistency. For large hWght we bias the interpolation - ! of T,S along the top and bottom integrals, almost like thickness - ! weighting. - ! Note: To work in terrain following coordinates we could offset - ! this distance by the layer thickness to replicate other models. - hWght = massWeightToggle * & - max(0., -bathyT(i,j)-e(i+1,j,K), -bathyT(i+1,j)-e(i,j,K)) - if (hWght > 0.) then - hL = (e(i,j,K) - e(i,j,K+1)) + dz_subroundoff - hR = (e(i+1,j,K) - e(i+1,j,K+1)) + dz_subroundoff - hWght = hWght * ( (hL-hR)/(hL+hR) )**2 - iDenom = 1./( hWght*(hR + hL) + hL*hR ) - Ttl = ( (hWght*hR)*T_t(i+1,j,k) + (hWght*hL + hR*hL)*T_t(i,j,k) ) * iDenom - Tbl = ( (hWght*hR)*T_b(i+1,j,k) + (hWght*hL + hR*hL)*T_b(i,j,k) ) * iDenom - Tml = ( (hWght*hR)*tv%T(i+1,j,k)+ (hWght*hL + hR*hL)*tv%T(i,j,k) ) * iDenom - Ttr = ( (hWght*hL)*T_t(i,j,k) + (hWght*hR + hR*hL)*T_t(i+1,j,k) ) * iDenom - Tbr = ( (hWght*hL)*T_b(i,j,k) + (hWght*hR + hR*hL)*T_b(i+1,j,k) ) * iDenom - Tmr = ( (hWght*hL)*tv%T(i,j,k) + (hWght*hR + hR*hL)*tv%T(i+1,j,k) ) * iDenom - Stl = ( (hWght*hR)*S_t(i+1,j,k) + (hWght*hL + hR*hL)*S_t(i,j,k) ) * iDenom - Sbl = ( (hWght*hR)*S_b(i+1,j,k) + (hWght*hL + hR*hL)*S_b(i,j,k) ) * iDenom - Sml = ( (hWght*hR)*tv%S(i+1,j,k) + (hWght*hL + hR*hL)*tv%S(i,j,k) ) * iDenom - Str = ( (hWght*hL)*S_t(i,j,k) + (hWght*hR + hR*hL)*S_t(i+1,j,k) ) * iDenom - Sbr = ( (hWght*hL)*S_b(i,j,k) + (hWght*hR + hR*hL)*S_b(i+1,j,k) ) * iDenom - Smr = ( (hWght*hL)*tv%S(i,j,k) + (hWght*hR + hR*hL)*tv%S(i+1,j,k) ) * iDenom - else - Ttl = T_t(i,j,k); Tbl = T_b(i,j,k); Ttr = T_t(i+1,j,k); Tbr = T_b(i+1,j,k) - Tml = tv%T(i,j,k); Tmr = tv%T(i+1,j,k) - Stl = S_t(i,j,k); Sbl = S_b(i,j,k); Str = S_t(i+1,j,k); Sbr = S_b(i+1,j,k) - Sml = tv%S(i,j,k); Smr = tv%S(i+1,j,k) - endif - - do m=2,4 - w_left = wt_t(m) ; w_right = wt_b(m) - - ! Salinity and temperature points are linearly interpolated in - ! the horizontal. The subscript (1) refers to the top value in - ! the vertical profile while subscript (5) refers to the bottom - ! value in the vertical profile. - T_top = w_left*Ttl + w_right*Ttr - T_mn = w_left*Tml + w_right*Tmr - T_bot = w_left*Tbl + w_right*Tbr - - S_top = w_left*Stl + w_right*Str - S_mn = w_left*Sml + w_right*Smr - S_bot = w_left*Sbl + w_right*Sbr - - ! Pressure - dz = w_left*(e(i,j,K) - e(i,j,K+1)) + w_right*(e(i+1,j,K) - e(i+1,j,K+1)) - p5(1) = -GxRho*((w_left*e(i,j,K) + w_right*e(i+1,j,K)) - z0pres) - do n=2,5 - p5(n) = p5(n-1) + GxRho*0.25*dz - enddo - - ! Parabolic reconstructions in the vertical for T and S - if (use_PPM) then - ! Coefficients of the parabolas - s6 = 3.0 * ( 2.0*S_mn - ( S_top + S_bot ) ) - t6 = 3.0 * ( 2.0*T_mn - ( T_top + T_bot ) ) - endif - do n=1,5 - S5(n) = wt_t(n) * S_top + wt_b(n) * ( S_bot + s6 * wt_t(n) ) - T5(n) = wt_t(n) * T_top + wt_b(n) * ( T_bot + t6 * wt_t(n) ) - enddo - if (use_stanley_eos) then - if (use_varT) T25(:) = w_left*tv%varT(i,j,k) + w_right*tv%varT(i+1,j,k) - if (use_covarTS) TS5(:) = w_left*tv%covarTS(i,j,k) + w_right*tv%covarTS(i+1,j,k) - if (use_varS) S25(:) = w_left*tv%varS(i,j,k) + w_right*tv%varS(i+1,j,k) - call calculate_density(T5, S5, p5, T25, TS5, S25, r5, EOS, rho_ref=rho_ref) + if (present(intx_dpa)) then ; do j=HI%jsc,HI%jec + do I=Isq,Ieq + ! Corner values of T and S + ! hWght is the distance measure by which the cell is violation of + ! hydrostatic consistency. For large hWght we bias the interpolation + ! of T,S along the top and bottom integrals, almost like thickness + ! weighting. + ! Note: To work in terrain following coordinates we could offset + ! this distance by the layer thickness to replicate other models. + hWght = massWeightToggle * & + max(0., -bathyT(i,j)-e(i+1,j,K), -bathyT(i+1,j)-e(i,j,K)) + if (hWght > 0.) then + hL = (e(i,j,K) - e(i,j,K+1)) + dz_subroundoff + hR = (e(i+1,j,K) - e(i+1,j,K+1)) + dz_subroundoff + hWght = hWght * ( (hL-hR)/(hL+hR) )**2 + iDenom = 1./( hWght*(hR + hL) + hL*hR ) + Ttl = ( (hWght*hR)*T_t(i+1,j,k) + (hWght*hL + hR*hL)*T_t(i,j,k) ) * iDenom + Tbl = ( (hWght*hR)*T_b(i+1,j,k) + (hWght*hL + hR*hL)*T_b(i,j,k) ) * iDenom + Tml = ( (hWght*hR)*tv%T(i+1,j,k)+ (hWght*hL + hR*hL)*tv%T(i,j,k) ) * iDenom + Ttr = ( (hWght*hL)*T_t(i,j,k) + (hWght*hR + hR*hL)*T_t(i+1,j,k) ) * iDenom + Tbr = ( (hWght*hL)*T_b(i,j,k) + (hWght*hR + hR*hL)*T_b(i+1,j,k) ) * iDenom + Tmr = ( (hWght*hL)*tv%T(i,j,k) + (hWght*hR + hR*hL)*tv%T(i+1,j,k) ) * iDenom + Stl = ( (hWght*hR)*S_t(i+1,j,k) + (hWght*hL + hR*hL)*S_t(i,j,k) ) * iDenom + Sbl = ( (hWght*hR)*S_b(i+1,j,k) + (hWght*hL + hR*hL)*S_b(i,j,k) ) * iDenom + Sml = ( (hWght*hR)*tv%S(i+1,j,k) + (hWght*hL + hR*hL)*tv%S(i,j,k) ) * iDenom + Str = ( (hWght*hL)*S_t(i,j,k) + (hWght*hR + hR*hL)*S_t(i+1,j,k) ) * iDenom + Sbr = ( (hWght*hL)*S_b(i,j,k) + (hWght*hR + hR*hL)*S_b(i+1,j,k) ) * iDenom + Smr = ( (hWght*hL)*tv%S(i,j,k) + (hWght*hR + hR*hL)*tv%S(i+1,j,k) ) * iDenom else - call calculate_density(T5, S5, p5, r5, EOS, rho_ref=rho_ref) + Ttl = T_t(i,j,k); Tbl = T_b(i,j,k); Ttr = T_t(i+1,j,k); Tbr = T_b(i+1,j,k) + Tml = tv%T(i,j,k); Tmr = tv%T(i+1,j,k) + Stl = S_t(i,j,k); Sbl = S_b(i,j,k); Str = S_t(i+1,j,k); Sbr = S_b(i+1,j,k) + Sml = tv%S(i,j,k); Smr = tv%S(i+1,j,k) endif - ! Use Boole's rule to estimate the pressure anomaly change. - intz(m) = G_e*dz*( C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) ) - enddo ! m - intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j) + do m=2,4 + w_left = wt_t(m) ; w_right = wt_b(m) - ! Use Boole's rule to integrate the bottom pressure anomaly values in x. - intx_dpa(I,j) = C1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + 12.0*intz(3)) + ! Salinity and temperature points are linearly interpolated in + ! the horizontal. The subscript (1) refers to the top value in + ! the vertical profile while subscript (5) refers to the bottom + ! value in the vertical profile. + T_top = w_left*Ttl + w_right*Ttr + T_mn = w_left*Tml + w_right*Tmr + T_bot = w_left*Tbl + w_right*Tbr - enddo ; enddo ; endif + S_top = w_left*Stl + w_right*Str + S_mn = w_left*Sml + w_right*Smr + S_bot = w_left*Sbl + w_right*Sbr - ! 3. Compute horizontal integrals in the y direction - if (present(inty_dpa)) then ; do J=Jsq,Jeq ; do i=HI%isc,HI%iec - ! Corner values of T and S - ! hWght is the distance measure by which the cell is violation of - ! hydrostatic consistency. For large hWght we bias the interpolation - ! of T,S along the top and bottom integrals, almost like thickness - ! weighting. - ! Note: To work in terrain following coordinates we could offset - ! this distance by the layer thickness to replicate other models. - hWght = massWeightToggle * & - max(0., -bathyT(i,j)-e(i,j+1,K), -bathyT(i,j+1)-e(i,j,K)) - if (hWght > 0.) then - hL = (e(i,j,K) - e(i,j,K+1)) + dz_subroundoff - hR = (e(i,j+1,K) - e(i,j+1,K+1)) + dz_subroundoff - hWght = hWght * ( (hL-hR)/(hL+hR) )**2 - iDenom = 1./( hWght*(hR + hL) + hL*hR ) - Ttl = ( (hWght*hR)*T_t(i,j+1,k) + (hWght*hL + hR*hL)*T_t(i,j,k) ) * iDenom - Tbl = ( (hWght*hR)*T_b(i,j+1,k) + (hWght*hL + hR*hL)*T_b(i,j,k) ) * iDenom - Tml = ( (hWght*hR)*tv%T(i,j+1,k)+ (hWght*hL + hR*hL)*tv%T(i,j,k) ) * iDenom - Ttr = ( (hWght*hL)*T_t(i,j,k) + (hWght*hR + hR*hL)*T_t(i,j+1,k) ) * iDenom - Tbr = ( (hWght*hL)*T_b(i,j,k) + (hWght*hR + hR*hL)*T_b(i,j+1,k) ) * iDenom - Tmr = ( (hWght*hL)*tv%T(i,j,k) + (hWght*hR + hR*hL)*tv%T(i,j+1,k) ) * iDenom - Stl = ( (hWght*hR)*S_t(i,j+1,k) + (hWght*hL + hR*hL)*S_t(i,j,k) ) * iDenom - Sbl = ( (hWght*hR)*S_b(i,j+1,k) + (hWght*hL + hR*hL)*S_b(i,j,k) ) * iDenom - Sml = ( (hWght*hR)*tv%S(i,j+1,k)+ (hWght*hL + hR*hL)*tv%S(i,j,k) ) * iDenom - Str = ( (hWght*hL)*S_t(i,j,k) + (hWght*hR + hR*hL)*S_t(i,j+1,k) ) * iDenom - Sbr = ( (hWght*hL)*S_b(i,j,k) + (hWght*hR + hR*hL)*S_b(i,j+1,k) ) * iDenom - Smr = ( (hWght*hL)*tv%S(i,j,k) + (hWght*hR + hR*hL)*tv%S(i,j+1,k) ) * iDenom + ! Pressure + dz_x(m,i) = w_left*(e(i,j,K) - e(i,j,K+1)) + w_right*(e(i+1,j,K) - e(i+1,j,K+1)) + + pos = i*15+(m-2)*5 + p15(pos+1) = -GxRho*((w_left*e(i,j,K) + w_right*e(i+1,j,K)) - z0pres) + do n=2,5 + p15(pos+n) = p15(pos+n-1) + GxRho*0.25*dz_x(m,i) + enddo + + ! Parabolic reconstructions in the vertical for T and S + if (use_PPM) then + ! Coefficients of the parabolas + s6 = 3.0 * ( 2.0*S_mn - ( S_top + S_bot ) ) + t6 = 3.0 * ( 2.0*T_mn - ( T_top + T_bot ) ) + endif + do n=1,5 + S15(pos+n) = wt_t(n) * S_top + wt_b(n) * ( S_bot + s6 * wt_t(n) ) + T15(pos+n) = wt_t(n) * T_top + wt_b(n) * ( T_bot + t6 * wt_t(n) ) + enddo + if (use_stanley_eos) then + if (use_varT) T215(pos+1:pos+5) = w_left*tv%varT(i,j,k) + w_right*tv%varT(i+1,j,k) + if (use_covarTS) TS15(pos+1:pos+5) = w_left*tv%covarTS(i,j,k) + w_right*tv%covarTS(i+1,j,k) + if (use_varS) S215(pos+1:pos+5) = w_left*tv%varS(i,j,k) + w_right*tv%varS(i+1,j,k) + endif + if (use_stanley_eos) then + call calculate_density(T5, S5, p5, T25, TS5, S25, r5, EOS, rho_ref=rho_ref) + else + call calculate_density(T5, S5, p5, r5, EOS, rho_ref=rho_ref) + endif + enddo + enddo + + if (use_stanley_eos) then + call calculate_density(T15, S15, p15, T215, TS15, S215, r15, EOS, EOSdom_q15, rho_ref=rho_ref) else - Ttl = T_t(i,j,k); Tbl = T_b(i,j,k); Ttr = T_t(i,j+1,k); Tbr = T_b(i,j+1,k) - Tml = tv%T(i,j,k); Tmr = tv%T(i,j+1,k) - Stl = S_t(i,j,k); Sbl = S_b(i,j,k); Str = S_t(i,j+1,k); Sbr = S_b(i,j+1,k) - Sml = tv%S(i,j,k); Smr = tv%S(i,j+1,k) + call calculate_density(T15, S15, p15, r15, EOS, EOSdom_q15, rho_ref=rho_ref) endif - do m=2,4 - w_left = wt_t(m) ; w_right = wt_b(m) - - ! Salinity and temperature points are linearly interpolated in - ! the horizontal. The subscript (1) refers to the top value in - ! the vertical profile while subscript (5) refers to the bottom - ! value in the vertical profile. - T_top = w_left*Ttl + w_right*Ttr - T_mn = w_left*Tml + w_right*Tmr - T_bot = w_left*Tbl + w_right*Tbr - - S_top = w_left*Stl + w_right*Str - S_mn = w_left*Sml + w_right*Smr - S_bot = w_left*Sbl + w_right*Sbr - - ! Pressure - dz = w_left*(e(i,j,K) - e(i,j,K+1)) + w_right*(e(i,j+1,K) - e(i,j+1,K+1)) - p5(1) = -GxRho*((w_left*e(i,j,K) + w_right*e(i,j+1,K)) - z0pres) - do n=2,5 - p5(n) = p5(n-1) + GxRho*0.25*dz - enddo + do I=Isq,Ieq + do m=2,4 + pos = i*15+(m-2)*5 + ! Use Boole's rule to estimate the pressure anomaly change. + intz(m) = G_e*dz_x(m,i)*( C1_90*( 7.0*(r15(pos+1)+r15(pos+5)) + & + 32.0*(r15(pos+2)+r15(pos+4)) + & + 12.0*r15(pos+3)) ) + enddo ! m + intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j) - ! Parabolic reconstructions in the vertical for T and S - if (use_PPM) then - ! Coefficients of the parabolas - s6 = 3.0 * ( 2.0*S_mn - ( S_top + S_bot ) ) - t6 = 3.0 * ( 2.0*T_mn - ( T_top + T_bot ) ) - endif - do n=1,5 - S5(n) = wt_t(n) * S_top + wt_b(n) * ( S_bot + s6 * wt_t(n) ) - T5(n) = wt_t(n) * T_top + wt_b(n) * ( T_bot + t6 * wt_t(n) ) - enddo + ! Use Boole's rule to integrate the bottom pressure anomaly values in x. + intx_dpa(I,j) = C1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + 12.0*intz(3)) - if (use_stanley_eos) then - if (use_varT) T25(:) = w_left*tv%varT(i,j,k) + w_right*tv%varT(i,j+1,k) - if (use_covarTS) TS5(:) = w_left*tv%covarTS(i,j,k) + w_right*tv%covarTS(i,j+1,k) - if (use_varS) S25(:) = w_left*tv%varS(i,j,k) + w_right*tv%varS(i,j+1,k) - call calculate_density(T5, S5, p5, T25, TS5, S25, r5, EOS, rho_ref=rho_ref) + enddo + enddo ; endif + + ! 3. Compute horizontal integrals in the y direction + if (present(inty_dpa)) then ; do J=Jsq,Jeq + do i=HI%isc,HI%iec + ! Corner values of T and S + ! hWght is the distance measure by which the cell is violation of + ! hydrostatic consistency. For large hWght we bias the interpolation + ! of T,S along the top and bottom integrals, almost like thickness + ! weighting. + ! Note: To work in terrain following coordinates we could offset + ! this distance by the layer thickness to replicate other models. + hWght = massWeightToggle * & + max(0., -bathyT(i,j)-e(i,j+1,K), -bathyT(i,j+1)-e(i,j,K)) + if (hWght > 0.) then + hL = (e(i,j,K) - e(i,j,K+1)) + dz_subroundoff + hR = (e(i,j+1,K) - e(i,j+1,K+1)) + dz_subroundoff + hWght = hWght * ( (hL-hR)/(hL+hR) )**2 + iDenom = 1./( hWght*(hR + hL) + hL*hR ) + Ttl = ( (hWght*hR)*T_t(i,j+1,k) + (hWght*hL + hR*hL)*T_t(i,j,k) ) * iDenom + Tbl = ( (hWght*hR)*T_b(i,j+1,k) + (hWght*hL + hR*hL)*T_b(i,j,k) ) * iDenom + Tml = ( (hWght*hR)*tv%T(i,j+1,k)+ (hWght*hL + hR*hL)*tv%T(i,j,k) ) * iDenom + Ttr = ( (hWght*hL)*T_t(i,j,k) + (hWght*hR + hR*hL)*T_t(i,j+1,k) ) * iDenom + Tbr = ( (hWght*hL)*T_b(i,j,k) + (hWght*hR + hR*hL)*T_b(i,j+1,k) ) * iDenom + Tmr = ( (hWght*hL)*tv%T(i,j,k) + (hWght*hR + hR*hL)*tv%T(i,j+1,k) ) * iDenom + Stl = ( (hWght*hR)*S_t(i,j+1,k) + (hWght*hL + hR*hL)*S_t(i,j,k) ) * iDenom + Sbl = ( (hWght*hR)*S_b(i,j+1,k) + (hWght*hL + hR*hL)*S_b(i,j,k) ) * iDenom + Sml = ( (hWght*hR)*tv%S(i,j+1,k)+ (hWght*hL + hR*hL)*tv%S(i,j,k) ) * iDenom + Str = ( (hWght*hL)*S_t(i,j,k) + (hWght*hR + hR*hL)*S_t(i,j+1,k) ) * iDenom + Sbr = ( (hWght*hL)*S_b(i,j,k) + (hWght*hR + hR*hL)*S_b(i,j+1,k) ) * iDenom + Smr = ( (hWght*hL)*tv%S(i,j,k) + (hWght*hR + hR*hL)*tv%S(i,j+1,k) ) * iDenom else - call calculate_density(T5, S5, p5, r5, EOS, rho_ref=rho_ref) + Ttl = T_t(i,j,k); Tbl = T_b(i,j,k); Ttr = T_t(i,j+1,k); Tbr = T_b(i,j+1,k) + Tml = tv%T(i,j,k); Tmr = tv%T(i,j+1,k) + Stl = S_t(i,j,k); Sbl = S_b(i,j,k); Str = S_t(i,j+1,k); Sbr = S_b(i,j+1,k) + Sml = tv%S(i,j,k); Smr = tv%S(i,j+1,k) endif - ! Use Boole's rule to estimate the pressure anomaly change. - intz(m) = G_e*dz*( C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) ) - enddo ! m - intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1) + do m=2,4 + w_left = wt_t(m) ; w_right = wt_b(m) + + ! Salinity and temperature points are linearly interpolated in + ! the horizontal. The subscript (1) refers to the top value in + ! the vertical profile while subscript (5) refers to the bottom + ! value in the vertical profile. + T_top = w_left*Ttl + w_right*Ttr + T_mn = w_left*Tml + w_right*Tmr + T_bot = w_left*Tbl + w_right*Tbr + + S_top = w_left*Stl + w_right*Str + S_mn = w_left*Sml + w_right*Smr + S_bot = w_left*Sbl + w_right*Sbr + + ! Pressure + dz_y(m,i) = w_left*(e(i,j,K) - e(i,j,K+1)) + w_right*(e(i,j+1,K) - e(i,j+1,K+1)) + + pos = i*15+(m-2)*5 + p15(pos+1) = -GxRho*((w_left*e(i,j,K) + w_right*e(i,j+1,K)) - z0pres) + do n=2,5 + p15(pos+n) = p15(pos+n-1) + GxRho*0.25*dz_y(m,i) + enddo + + ! Parabolic reconstructions in the vertical for T and S + if (use_PPM) then + ! Coefficients of the parabolas + s6 = 3.0 * ( 2.0*S_mn - ( S_top + S_bot ) ) + t6 = 3.0 * ( 2.0*T_mn - ( T_top + T_bot ) ) + endif + do n=1,5 + S15(pos+n) = wt_t(n) * S_top + wt_b(n) * ( S_bot + s6 * wt_t(n) ) + T15(pos+n) = wt_t(n) * T_top + wt_b(n) * ( T_bot + t6 * wt_t(n) ) + enddo + + if (use_stanley_eos) then + if (use_varT) T215(pos+1:pos+5) = w_left*tv%varT(i,j,k) + w_right*tv%varT(i,j+1,k) + if (use_covarTS) TS15(pos+1:pos+5) = w_left*tv%covarTS(i,j,k) + w_right*tv%covarTS(i,j+1,k) + if (use_varS) S215(pos+1:pos+5) = w_left*tv%varS(i,j,k) + w_right*tv%varS(i,j+1,k) + endif + enddo + enddo + + if (use_stanley_eos) then + call calculate_density(T15(15*HI%isc+1:), S15(15*HI%isc+1:), p15(15*HI%isc+1:), & + T215(15*HI%isc+1:), TS15(15*HI%isc+1:), S215(15*HI%isc+1:), & + r15(15*HI%isc+1:), EOS, EOSdom_h15, rho_ref=rho_ref) + else + call calculate_density(T15(15*HI%isc+1:), S15(15*HI%isc+1:), p15(15*HI%isc+1:), & + r15(15*HI%isc+1:), EOS, EOSdom_h15, rho_ref=rho_ref) + endif - ! Use Boole's rule to integrate the bottom pressure anomaly values in y. - inty_dpa(i,J) = C1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + 12.0*intz(3)) + do i=HI%isc,HI%iec + do m=2,4 + ! Use Boole's rule to estimate the pressure anomaly change. + pos = i*15+(m-2)*5 + intz(m) = G_e*dz_y(m,i)*( C1_90*( 7.0*(r15(pos+1)+r15(pos+5)) + & + 32.0*(r15(pos+2)+r15(pos+4)) + & + 12.0*r15(pos+3)) ) + enddo ! m + intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1) - enddo ; enddo ; endif + ! Use Boole's rule to integrate the bottom pressure anomaly values in y. + inty_dpa(i,J) = C1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + 12.0*intz(3)) + enddo + enddo ; endif end subroutine int_density_dz_generic_ppm @@ -1161,12 +1285,19 @@ subroutine int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, d ! series for log(1-eps/1+eps) that assumes that |eps| < 0.34. ! Local variables - real :: T5(5) ! Temperatures at five quadrature points [C ~> degC] - real :: S5(5) ! Salinities at five quadrature points [S ~> ppt] - real :: p5(5) ! Pressures at five quadrature points [R L2 T-2 ~> Pa] - real :: a5(5) ! Specific volumes at five quadrature points [R-1 ~> m3 kg-1] + real :: T5((5*HI%isd+1):(5*(HI%ied+2))) ! Temperatures along a line of subgrid locations [C ~> degC] + real :: S5((5*HI%isd+1):(5*(HI%ied+2))) ! Salinities along a line of subgrid locations [S ~> ppt] + real :: p5((5*HI%isd+1):(5*(HI%ied+2))) ! Pressures along a line of subgrid locations [R L2 T-2 ~> Pa] + real :: a5((5*HI%isd+1):(5*(HI%ied+2))) ! Specific volumes anomalies along a line of subgrid + ! locations [R-1 ~> m3 kg-3] + real :: T15((15*HI%isd+1):(15*(HI%ied+1))) ! Temperatures at an array of subgrid locations [C ~> degC] + real :: S15((15*HI%isd+1):(15*(HI%ied+1))) ! Salinities at an array of subgrid locations [S ~> ppt] + real :: p15((15*HI%isd+1):(15*(HI%ied+1))) ! Pressures at an array of subgrid locations [R L2 T-2 ~> Pa] + real :: a15((15*HI%isd+1):(15*(HI%ied+1))) ! Specific volumes at an array of subgrid locations [R ~> kg m-3] real :: alpha_anom ! The depth averaged specific density anomaly [R-1 ~> m3 kg-1] real :: dp ! The pressure change through a layer [R L2 T-2 ~> Pa] + real :: dp_x(5,SZIB_(HI)) ! The pressure change through a layer along an x-line of subgrid locations [Z ~> m] + real :: dp_y(5,SZI_(HI)) ! The pressure change through a layer along a y-line of subgrid locations [Z ~> m] real :: hWght ! A pressure-thickness below topography [R L2 T-2 ~> Pa] real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa] real :: iDenom ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2] @@ -1178,7 +1309,10 @@ subroutine int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, d ! 5 sub-column locations [L2 T-2 ~> m2 s-2] logical :: do_massWeight ! Indicates whether to do mass weighting. real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim] - integer :: Isq, Ieq, Jsq, Jeq, ish, ieh, jsh, jeh, i, j, m, n, halo + integer, dimension(2) :: EOSdom_h5 ! The 5-point h-point i-computational domain for the equation of state + integer, dimension(2) :: EOSdom_q15 ! The 3x5-point q-point i-computational domain for the equation of state + integer, dimension(2) :: EOSdom_h15 ! The 3x5-point h-point i-computational domain for the equation of state + integer :: Isq, Ieq, Jsq, Jeq, ish, ieh, jsh, jeh, i, j, m, n, pos, halo Isq = HI%IscB ; Ieq = HI%IecB ; Jsq = HI%JscB ; Jeq = HI%JecB halo = 0 ; if (present(halo_size)) halo = MAX(halo_size,0) @@ -1195,110 +1329,146 @@ subroutine int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, d "dP_neglect must be present if useMassWghtInterp is present and true.") endif ; endif - do j=jsh,jeh ; do i=ish,ieh - dp = p_b(i,j) - p_t(i,j) - do n=1,5 - T5(n) = T(i,j) ; S5(n) = S(i,j) - p5(n) = p_b(i,j) - 0.25*real(n-1)*dp + ! Set the loop ranges for equation of state calculations at various points. + EOSdom_h5(1) = 1 ; EOSdom_h5(2) = 5*(ieh-ish+1) + EOSdom_q15(1) = 1 ; EOSdom_q15(2) = 15*(Ieq-Isq+1) + EOSdom_h15(1) = 1 ; EOSdom_h15(2) = 15*(HI%iec-HI%isc+1) + + do j=jsh,jeh + do i=ish,ieh + dp = p_b(i,j) - p_t(i,j) + pos = 5*i + do n=1,5 + T5(pos+n) = T(i,j) ; S5(pos+n) = S(i,j) + p5(pos+n) = p_b(i,j) - 0.25*real(n-1)*dp + enddo enddo - call calculate_spec_vol(T5, S5, p5, a5, EOS, spv_ref=alpha_ref) + call calculate_spec_vol(T5(5*ish+1:), S5(5*ish+1:), p5(5*ish+1:), a5(5*ish+1:), EOS, & + EOSdom_h5, spv_ref=alpha_ref) - ! Use Boole's rule to estimate the interface height anomaly change. - alpha_anom = C1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + 12.0*a5(3)) - dza(i,j) = dp*alpha_anom - ! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of - ! the interface height anomaly. - if (present(intp_dza)) intp_dza(i,j) = 0.5*dp**2 * & - (alpha_anom - C1_90*(16.0*(a5(4)-a5(2)) + 7.0*(a5(5)-a5(1))) ) - enddo ; enddo + do i=ish,ieh + dp = p_b(i,j) - p_t(i,j) + ! Use Boole's rule to estimate the interface height anomaly change. + pos = 5*i + alpha_anom = C1_90*(7.0*(a5(pos+1)+a5(pos+5)) + 32.0*(a5(pos+2)+a5(pos+4)) + 12.0*a5(pos+3)) + dza(i,j) = dp*alpha_anom + ! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of + ! the interface height anomaly. + if (present(intp_dza)) intp_dza(i,j) = 0.5*dp**2 * & + (alpha_anom - C1_90*(16.0*(a5(pos+4)-a5(pos+2)) + 7.0*(a5(pos+5)-a5(pos+1))) ) + enddo + enddo - if (present(intx_dza)) then ; do j=HI%jsc,HI%jec ; do I=Isq,Ieq - ! hWght is the distance measure by which the cell is violation of - ! hydrostatic consistency. For large hWght we bias the interpolation of - ! T & S along the top and bottom integrals, akin to thickness weighting. - hWght = 0.0 - if (do_massWeight) & - hWght = max(0., bathyP(i,j)-p_t(i+1,j), bathyP(i+1,j)-p_t(i,j)) - if (hWght > 0.) then - hL = (p_b(i,j) - p_t(i,j)) + dP_neglect - hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect - hWght = hWght * ( (hL-hR)/(hL+hR) )**2 - iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR ) - hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom - hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom - else - hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0 - endif + if (present(intx_dza)) then ; do j=HI%jsc,HI%jec + do I=Isq,Ieq + ! hWght is the distance measure by which the cell is violation of + ! hydrostatic consistency. For large hWght we bias the interpolation of + ! T & S along the top and bottom integrals, akin to thickness weighting. + hWght = 0.0 + if (do_massWeight) & + hWght = max(0., bathyP(i,j)-p_t(i+1,j), bathyP(i+1,j)-p_t(i,j)) + if (hWght > 0.) then + hL = (p_b(i,j) - p_t(i,j)) + dP_neglect + hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect + hWght = hWght * ( (hL-hR)/(hL+hR) )**2 + iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR ) + hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom + hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom + else + hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0 + endif - intp(1) = dza(i,j) ; intp(5) = dza(i+1,j) - do m=2,4 - wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L - wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR + do m=2,4 + wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L + wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR + pos = i*15+(m-2)*5 - ! T, S, and p are interpolated in the horizontal. The p interpolation - ! is linear, but for T and S it may be thickness weighted. - p5(1) = wt_L*p_b(i,j) + wt_R*p_b(i+1,j) - dp = wt_L*(p_b(i,j) - p_t(i,j)) + wt_R*(p_b(i+1,j) - p_t(i+1,j)) - T5(1) = wtT_L*T(i,j) + wtT_R*T(i+1,j) - S5(1) = wtT_L*S(i,j) + wtT_R*S(i+1,j) + ! T, S, and p are interpolated in the horizontal. The p interpolation + ! is linear, but for T and S it may be thickness weighted. + p15(pos+1) = wt_L*p_b(i,j) + wt_R*p_b(i+1,j) + dp_x(m,I) = wt_L*(p_b(i,j) - p_t(i,j)) + wt_R*(p_b(i+1,j) - p_t(i+1,j)) + T15(pos+1) = wtT_L*T(i,j) + wtT_R*T(i+1,j) + S15(pos+1) = wtT_L*S(i,j) + wtT_R*S(i+1,j) - do n=2,5 - T5(n) = T5(1) ; S5(n) = S5(1) ; p5(n) = p5(n-1) - 0.25*dp + do n=2,5 + T15(pos+n) = T15(pos+1) ; S15(pos+n) = S15(pos+1) + p15(pos+n) = p15(pos+n-1) - 0.25*dp_x(m,I) + enddo enddo - call calculate_spec_vol(T5, S5, p5, a5, EOS, spv_ref=alpha_ref) + enddo + + call calculate_spec_vol(T15(15*Isq+1:), S15(15*Isq+1:), p15(15*Isq+1:), & + a15(15*Isq+1:), EOS, EOSdom_q15, spv_ref=alpha_ref) - ! Use Boole's rule to estimate the interface height anomaly change. - intp(m) = dp*( C1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + & - 12.0*a5(3))) + do I=Isq,Ieq + intp(1) = dza(i,j) ; intp(5) = dza(i+1,j) + ! Use Boole's rule to estimate the interface height anomaly change. + do m=2,4 + pos = i*15+(m-2)*5 + intp(m) = dp_x(m,I)*( C1_90*(7.0*(a15(pos+1)+a15(pos+5)) + 32.0*(a15(pos+2)+a15(pos+4)) + & + 12.0*a15(pos+3))) + enddo + ! Use Boole's rule to integrate the interface height anomaly values in x. + intx_dza(i,j) = C1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + & + 12.0*intp(3)) enddo - ! Use Boole's rule to integrate the interface height anomaly values in x. - intx_dza(i,j) = C1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + & - 12.0*intp(3)) - enddo ; enddo ; endif + enddo ; endif - if (present(inty_dza)) then ; do J=Jsq,Jeq ; do i=HI%isc,HI%iec - ! hWght is the distance measure by which the cell is violation of - ! hydrostatic consistency. For large hWght we bias the interpolation of - ! T & S along the top and bottom integrals, akin to thickness weighting. - hWght = 0.0 - if (do_massWeight) & - hWght = max(0., bathyP(i,j)-p_t(i,j+1), bathyP(i,j+1)-p_t(i,j)) - if (hWght > 0.) then - hL = (p_b(i,j) - p_t(i,j)) + dP_neglect - hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect - hWght = hWght * ( (hL-hR)/(hL+hR) )**2 - iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR ) - hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom - hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom - else - hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0 - endif + if (present(inty_dza)) then ; do J=Jsq,Jeq + do i=HI%isc,HI%iec + ! hWght is the distance measure by which the cell is violation of + ! hydrostatic consistency. For large hWght we bias the interpolation of + ! T & S along the top and bottom integrals, akin to thickness weighting. + hWght = 0.0 + if (do_massWeight) & + hWght = max(0., bathyP(i,j)-p_t(i,j+1), bathyP(i,j+1)-p_t(i,j)) + if (hWght > 0.) then + hL = (p_b(i,j) - p_t(i,j)) + dP_neglect + hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect + hWght = hWght * ( (hL-hR)/(hL+hR) )**2 + iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR ) + hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom + hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom + else + hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0 + endif + + do m=2,4 + wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L + wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR + pos = i*15+(m-2)*5 - intp(1) = dza(i,j) ; intp(5) = dza(i,j+1) - do m=2,4 - wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L - wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR - - ! T, S, and p are interpolated in the horizontal. The p interpolation - ! is linear, but for T and S it may be thickness weighted. - p5(1) = wt_L*p_b(i,j) + wt_R*p_b(i,j+1) - dp = wt_L*(p_b(i,j) - p_t(i,j)) + wt_R*(p_b(i,j+1) - p_t(i,j+1)) - T5(1) = wtT_L*T(i,j) + wtT_R*T(i,j+1) - S5(1) = wtT_L*S(i,j) + wtT_R*S(i,j+1) - do n=2,5 - T5(n) = T5(1) ; S5(n) = S5(1) ; p5(n) = p5(n-1) - 0.25*dp + ! T, S, and p are interpolated in the horizontal. The p interpolation + ! is linear, but for T and S it may be thickness weighted. + p15(pos+1) = wt_L*p_b(i,j) + wt_R*p_b(i,j+1) + dp_y(m,i) = wt_L*(p_b(i,j) - p_t(i,j)) + wt_R*(p_b(i,j+1) - p_t(i,j+1)) + T15(pos+1) = wtT_L*T(i,j) + wtT_R*T(i,j+1) + S15(pos+1) = wtT_L*S(i,j) + wtT_R*S(i,j+1) + do n=2,5 + T15(pos+n) = T15(pos+1) ; S15(pos+n) = S15(pos+1) + p15(pos+n) = p15(pos+n-1) - 0.25*dp_y(m,i) + enddo enddo - call calculate_spec_vol(T5, S5, p5, a5, EOS, spv_ref=alpha_ref) + enddo + + call calculate_spec_vol(T15(15*HI%isc+1:), S15(15*HI%isc+1:), p15(15*HI%isc+1:), & + a15(15*HI%isc+1:), EOS, EOSdom_h15, spv_ref=alpha_ref) + + do i=HI%isc,HI%iec - ! Use Boole's rule to estimate the interface height anomaly change. - intp(m) = dp*( C1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + & - 12.0*a5(3))) + intp(1) = dza(i,j) ; intp(5) = dza(i,j+1) + ! Use Boole's rule to estimate the interface height anomaly change. + do m=2,4 + pos = i*15+(m-2)*5 + intp(m) = dp_y(m,i)*( C1_90*(7.0*(a15(pos+1)+a15(pos+5)) + 32.0*(a15(pos+2)+a15(pos+4)) + & + 12.0*a15(pos+3))) + enddo + ! Use Boole's rule to integrate the interface height anomaly values in y. + inty_dza(i,j) = C1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + & + 12.0*intp(3)) enddo - ! Use Boole's rule to integrate the interface height anomaly values in y. - inty_dza(i,j) = C1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + & - 12.0*intp(3)) - enddo ; enddo ; endif + enddo ; endif end subroutine int_spec_vol_dp_generic_pcm @@ -1358,14 +1528,15 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, ! Boole's rule to do the horizontal integrals, and from a truncation in the ! series for log(1-eps/1+eps) that assumes that |eps| < 0.34. - real :: T5(5) ! Temperatures at five quadrature points [C ~> degC] - real :: S5(5) ! Salinities at five quadrature points [S ~> ppt] - real :: p5(5) ! Pressures at five quadrature points [R L2 T-2 ~> Pa] - real :: a5(5) ! Specific volumes at five quadrature points [R-1 ~> m3 kg-1] - real :: T15(15) ! Temperatures at fifteen interior quadrature points [C ~> degC] - real :: S15(15) ! Salinities at fifteen interior quadrature points [S ~> ppt] - real :: p15(15) ! Pressures at fifteen quadrature points [R L2 T-2 ~> Pa] - real :: a15(15) ! Specific volumes at fifteen quadrature points [R-1 ~> m3 kg-1] + real :: T5((5*HI%iscB+1):(5*(HI%iecB+2))) ! Temperatures along a line of subgrid locations [C ~> degC] + real :: S5((5*HI%iscB+1):(5*(HI%iecB+2))) ! Salinities along a line of subgrid locations [S ~> ppt] + real :: p5((5*HI%iscB+1):(5*(HI%iecB+2))) ! Pressures along a line of subgrid locations [R L2 T-2 ~> Pa] + real :: a5((5*HI%iscB+1):(5*(HI%iecB+2))) ! Specific volumes anomalies along a line of subgrid + ! locations [R-1 ~> m3 kg-3] + real :: T15((15*HI%iscB+1):(15*(HI%iecB+1))) ! Temperatures at an array of subgrid locations [C ~> degC] + real :: S15((15*HI%iscB+1):(15*(HI%iecB+1))) ! Salinities at an array of subgrid locations [S ~> ppt] + real :: p15((15*HI%iscB+1):(15*(HI%iecB+1))) ! Pressures at an array of subgrid locations [R L2 T-2 ~> Pa] + real :: a15((15*HI%iscB+1):(15*(HI%iecB+1))) ! Specific volumes at an array of subgrid locations [R ~> kg m-3] real :: wt_t(5), wt_b(5) ! Weights of top and bottom values at quadrature points [nondim] real :: T_top, T_bot ! Horizontally interpolated temperature at the cell top and bottom [C ~> degC] real :: S_top, S_bot ! Horizontally interpolated salinity at the cell top and bottom [S ~> ppt] @@ -1373,7 +1544,7 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, real :: alpha_anom ! The depth averaged specific density anomaly [R-1 ~> m3 kg-1] real :: dp ! The pressure change through a layer [R L2 T-2 ~> Pa] - real :: dp_90(2:4) ! The pressure change through a layer divided by 90 [R L2 T-2 ~> Pa] + real :: dp_90(2:4,SZIB_(HI)) ! The pressure change through a layer divided by 90 [R L2 T-2 ~> Pa] real :: hWght ! A pressure-thickness below topography [R L2 T-2 ~> Pa] real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa] real :: iDenom ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2] @@ -1385,6 +1556,9 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, ! 5 sub-column locations [L2 T-2 ~> m2 s-2] real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim] logical :: do_massWeight ! Indicates whether to do mass weighting. + integer, dimension(2) :: EOSdom_h5 ! The 5-point h-point i-computational domain for the equation of state + integer, dimension(2) :: EOSdom_q15 ! The 3x5-point q-point i-computational domain for the equation of state + integer, dimension(2) :: EOSdom_h15 ! The 3x5-point h-point i-computational domain for the equation of state integer :: Isq, Ieq, Jsq, Jeq, i, j, m, n, pos Isq = HI%IscB ; Ieq = HI%IecB ; Jsq = HI%JscB ; Jeq = HI%JecB @@ -1397,140 +1571,157 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, wt_b(n) = 1.0 - wt_t(n) enddo + ! Set the loop ranges for equation of state calculations at various points. + EOSdom_h5(1) = 1 ; EOSdom_h5(2) = 5*(Ieq-Isq+2) + EOSdom_q15(1) = 1 ; EOSdom_q15(2) = 15*(Ieq-Isq+1) + EOSdom_h15(1) = 1 ; EOSdom_h15(2) = 15*(HI%iec-HI%isc+1) + ! 1. Compute vertical integrals - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - dp = p_b(i,j) - p_t(i,j) - do n=1,5 ! T, S and p are linearly interpolated in the vertical. - p5(n) = wt_t(n) * p_t(i,j) + wt_b(n) * p_b(i,j) - S5(n) = wt_t(n) * S_t(i,j) + wt_b(n) * S_b(i,j) - T5(n) = wt_t(n) * T_t(i,j) + wt_b(n) * T_b(i,j) + do j=Jsq,Jeq+1 + do i=Isq,Ieq+1 + do n=1,5 ! T, S and p are linearly interpolated in the vertical. + p5(i*5+n) = wt_t(n) * p_t(i,j) + wt_b(n) * p_b(i,j) + S5(i*5+n) = wt_t(n) * S_t(i,j) + wt_b(n) * S_b(i,j) + T5(i*5+n) = wt_t(n) * T_t(i,j) + wt_b(n) * T_b(i,j) + enddo enddo - call calculate_spec_vol(T5, S5, p5, a5, EOS, spv_ref=alpha_ref) - - ! Use Boole's rule to estimate the interface height anomaly change. - alpha_anom = C1_90*((7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4))) + 12.0*a5(3)) - dza(i,j) = dp*alpha_anom - ! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of - ! the interface height anomaly. - if (present(intp_dza)) intp_dza(i,j) = 0.5*dp**2 * & - (alpha_anom - C1_90*(16.0*(a5(4)-a5(2)) + 7.0*(a5(5)-a5(1))) ) - enddo ; enddo + call calculate_spec_vol(T5, S5, p5, a5, EOS, EOSdom_h5, spv_ref=alpha_ref) + do i=Isq,Ieq+1 + ! Use Boole's rule to estimate the interface height anomaly change. + dp = p_b(i,j) - p_t(i,j) + alpha_anom = C1_90*((7.0*(a5(i*5+1)+a5(i*5+5)) + 32.0*(a5(i*5+2)+a5(i*5+4))) + 12.0*a5(i*5+3)) + dza(i,j) = dp*alpha_anom + ! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of + ! the interface height anomaly. + if (present(intp_dza)) intp_dza(i,j) = 0.5*dp**2 * & + (alpha_anom - C1_90*(16.0*(a5(i*5+4)-a5(i*5+2)) + 7.0*(a5(i*5+5)-a5(i*5+1))) ) + enddo + enddo ! 2. Compute horizontal integrals in the x direction - if (present(intx_dza)) then ; do j=HI%jsc,HI%jec ; do I=Isq,Ieq - ! hWght is the distance measure by which the cell is violation of - ! hydrostatic consistency. For large hWght we bias the interpolation - ! of T,S along the top and bottom integrals, almost like thickness - ! weighting. Note: To work in terrain following coordinates we could - ! offset this distance by the layer thickness to replicate other models. - hWght = 0.0 - if (do_massWeight) & - hWght = max(0., bathyP(i,j)-p_t(i+1,j), bathyP(i+1,j)-p_t(i,j)) - if (hWght > 0.) then - hL = (p_b(i,j) - p_t(i,j)) + dP_neglect - hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect - hWght = hWght * ( (hL-hR)/(hL+hR) )**2 - iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR ) - hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom - hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom - else - hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0 - endif + if (present(intx_dza)) then ; do j=HI%jsc,HI%jec + do I=Isq,Ieq + ! hWght is the distance measure by which the cell is violation of + ! hydrostatic consistency. For large hWght we bias the interpolation + ! of T,S along the top and bottom integrals, almost like thickness + ! weighting. Note: To work in terrain following coordinates we could + ! offset this distance by the layer thickness to replicate other models. + hWght = 0.0 + if (do_massWeight) & + hWght = max(0., bathyP(i,j)-p_t(i+1,j), bathyP(i+1,j)-p_t(i,j)) + if (hWght > 0.) then + hL = (p_b(i,j) - p_t(i,j)) + dP_neglect + hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect + hWght = hWght * ( (hL-hR)/(hL+hR) )**2 + iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR ) + hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom + hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom + else + hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0 + endif - do m=2,4 - wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L - wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR - - ! T, S, and p are interpolated in the horizontal. The p interpolation - ! is linear, but for T and S it may be thickness weighted. - P_top = wt_L*p_t(i,j) + wt_R*p_t(i+1,j) - P_bot = wt_L*p_b(i,j) + wt_R*p_b(i+1,j) - T_top = wtT_L*T_t(i,j) + wtT_R*T_t(i+1,j) - T_bot = wtT_L*T_b(i,j) + wtT_R*T_b(i+1,j) - S_top = wtT_L*S_t(i,j) + wtT_R*S_t(i+1,j) - S_bot = wtT_L*S_b(i,j) + wtT_R*S_b(i+1,j) - dp_90(m) = C1_90*(P_bot - P_top) - - ! Salinity, temperature and pressure with linear interpolation in the vertical. - pos = (m-2)*5 - do n=1,5 - p15(pos+n) = wt_t(n) * P_top + wt_b(n) * P_bot - S15(pos+n) = wt_t(n) * S_top + wt_b(n) * S_bot - T15(pos+n) = wt_t(n) * T_top + wt_b(n) * T_bot + do m=2,4 + wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L + wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR + + ! T, S, and p are interpolated in the horizontal. The p interpolation + ! is linear, but for T and S it may be thickness weighted. + P_top = wt_L*p_t(i,j) + wt_R*p_t(i+1,j) + P_bot = wt_L*p_b(i,j) + wt_R*p_b(i+1,j) + T_top = wtT_L*T_t(i,j) + wtT_R*T_t(i+1,j) + T_bot = wtT_L*T_b(i,j) + wtT_R*T_b(i+1,j) + S_top = wtT_L*S_t(i,j) + wtT_R*S_t(i+1,j) + S_bot = wtT_L*S_b(i,j) + wtT_R*S_b(i+1,j) + dp_90(m,I) = C1_90*(P_bot - P_top) + + ! Salinity, temperature and pressure with linear interpolation in the vertical. + pos = i*15+(m-2)*5 + do n=1,5 + p15(pos+n) = wt_t(n) * P_top + wt_b(n) * P_bot + S15(pos+n) = wt_t(n) * S_top + wt_b(n) * S_bot + T15(pos+n) = wt_t(n) * T_top + wt_b(n) * T_bot + enddo enddo enddo - call calculate_spec_vol(T15, S15, p15, a15, EOS, spv_ref=alpha_ref) + call calculate_spec_vol(T15, S15, p15, a15, EOS, EOSdom_q15, spv_ref=alpha_ref) - intp(1) = dza(i,j) ; intp(5) = dza(i+1,j) - do m=2,4 - ! Use Boole's rule to estimate the interface height anomaly change. - ! The integrals at the ends of the segment are already known. - pos = (m-2)*5 - intp(m) = dp_90(m)*((7.0*(a15(pos+1)+a15(pos+5)) + & - 32.0*(a15(pos+2)+a15(pos+4))) + 12.0*a15(pos+3)) + do I=Isq,Ieq + intp(1) = dza(i,j) ; intp(5) = dza(i+1,j) + do m=2,4 + ! Use Boole's rule to estimate the interface height anomaly change. + ! The integrals at the ends of the segment are already known. + pos = I*15+(m-2)*5 + intp(m) = dp_90(m,I)*((7.0*(a15(pos+1)+a15(pos+5)) + & + 32.0*(a15(pos+2)+a15(pos+4))) + 12.0*a15(pos+3)) + enddo + ! Use Boole's rule to integrate the interface height anomaly values in x. + intx_dza(I,j) = C1_90*((7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4))) + & + 12.0*intp(3)) enddo - ! Use Boole's rule to integrate the interface height anomaly values in x. - intx_dza(I,j) = C1_90*((7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4))) + & - 12.0*intp(3)) - enddo ; enddo ; endif + enddo ; endif ! 3. Compute horizontal integrals in the y direction - if (present(inty_dza)) then ; do J=Jsq,Jeq ; do i=HI%isc,HI%iec - ! hWght is the distance measure by which the cell is violation of - ! hydrostatic consistency. For large hWght we bias the interpolation - ! of T,S along the top and bottom integrals, like thickness weighting. - hWght = 0.0 - if (do_massWeight) & - hWght = max(0., bathyP(i,j)-p_t(i,j+1), bathyP(i,j+1)-p_t(i,j)) - if (hWght > 0.) then - hL = (p_b(i,j) - p_t(i,j)) + dP_neglect - hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect - hWght = hWght * ( (hL-hR)/(hL+hR) )**2 - iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR ) - hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom - hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom - else - hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0 - endif + if (present(inty_dza)) then ; do J=Jsq,Jeq + do i=HI%isc,HI%iec + ! hWght is the distance measure by which the cell is violation of + ! hydrostatic consistency. For large hWght we bias the interpolation + ! of T,S along the top and bottom integrals, like thickness weighting. + hWght = 0.0 + if (do_massWeight) & + hWght = max(0., bathyP(i,j)-p_t(i,j+1), bathyP(i,j+1)-p_t(i,j)) + if (hWght > 0.) then + hL = (p_b(i,j) - p_t(i,j)) + dP_neglect + hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect + hWght = hWght * ( (hL-hR)/(hL+hR) )**2 + iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR ) + hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom + hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom + else + hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0 + endif - do m=2,4 - wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L - wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR - - ! T, S, and p are interpolated in the horizontal. The p interpolation - ! is linear, but for T and S it may be thickness weighted. - P_top = wt_L*p_t(i,j) + wt_R*p_t(i,j+1) - P_bot = wt_L*p_b(i,j) + wt_R*p_b(i,j+1) - T_top = wtT_L*T_t(i,j) + wtT_R*T_t(i,j+1) - T_bot = wtT_L*T_b(i,j) + wtT_R*T_b(i,j+1) - S_top = wtT_L*S_t(i,j) + wtT_R*S_t(i,j+1) - S_bot = wtT_L*S_b(i,j) + wtT_R*S_b(i,j+1) - dp_90(m) = C1_90*(P_bot - P_top) - - ! Salinity, temperature and pressure with linear interpolation in the vertical. - pos = (m-2)*5 - do n=1,5 - p15(pos+n) = wt_t(n) * P_top + wt_b(n) * P_bot - S15(pos+n) = wt_t(n) * S_top + wt_b(n) * S_bot - T15(pos+n) = wt_t(n) * T_top + wt_b(n) * T_bot + do m=2,4 + wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L + wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR + + ! T, S, and p are interpolated in the horizontal. The p interpolation + ! is linear, but for T and S it may be thickness weighted. + P_top = wt_L*p_t(i,j) + wt_R*p_t(i,j+1) + P_bot = wt_L*p_b(i,j) + wt_R*p_b(i,j+1) + T_top = wtT_L*T_t(i,j) + wtT_R*T_t(i,j+1) + T_bot = wtT_L*T_b(i,j) + wtT_R*T_b(i,j+1) + S_top = wtT_L*S_t(i,j) + wtT_R*S_t(i,j+1) + S_bot = wtT_L*S_b(i,j) + wtT_R*S_b(i,j+1) + dp_90(m,i) = C1_90*(P_bot - P_top) + + ! Salinity, temperature and pressure with linear interpolation in the vertical. + pos = i*15+(m-2)*5 + do n=1,5 + p15(pos+n) = wt_t(n) * P_top + wt_b(n) * P_bot + S15(pos+n) = wt_t(n) * S_top + wt_b(n) * S_bot + T15(pos+n) = wt_t(n) * T_top + wt_b(n) * T_bot + enddo enddo enddo - call calculate_spec_vol(T15, S15, p15, a15, EOS, spv_ref=alpha_ref) + call calculate_spec_vol(T15(15*HI%isc+1:), S15(15*HI%isc+1:), p15(15*HI%isc+1:), & + a15(15*HI%isc+1:), EOS, EOSdom_h15, spv_ref=alpha_ref) - intp(1) = dza(i,j) ; intp(5) = dza(i,j+1) - do m=2,4 - ! Use Boole's rule to estimate the interface height anomaly change. - ! The integrals at the ends of the segment are already known. - pos = (m-2)*5 - intp(m) = dp_90(m) * ((7.0*(a15(pos+1)+a15(pos+5)) + & - 32.0*(a15(pos+2)+a15(pos+4))) + 12.0*a15(pos+3)) + do i=HI%isc,HI%iec + intp(1) = dza(i,j) ; intp(5) = dza(i,j+1) + do m=2,4 + ! Use Boole's rule to estimate the interface height anomaly change. + ! The integrals at the ends of the segment are already known. + pos = i*15+(m-2)*5 + intp(m) = dp_90(m,i) * ((7.0*(a15(pos+1)+a15(pos+5)) + & + 32.0*(a15(pos+2)+a15(pos+4))) + 12.0*a15(pos+3)) + enddo + ! Use Boole's rule to integrate the interface height anomaly values in x. + inty_dza(i,J) = C1_90*((7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4))) + & + 12.0*intp(3)) enddo - ! Use Boole's rule to integrate the interface height anomaly values in x. - inty_dza(i,J) = C1_90*((7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4))) + & - 12.0*intp(3)) - enddo ; enddo ; endif + enddo ; endif end subroutine int_spec_vol_dp_generic_plm diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index aeb25bc351..167b2e81f3 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -999,7 +999,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie KE_term(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) & - * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) + * ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1))) enddo ; enddo enddo call post_data(CS%id_dKEdt, KE_term, CS%diag) @@ -1018,7 +1018,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie KE_term(i,j,k) = 0.5 * G%IareaT(i,j) & - * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) + * ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1))) enddo ; enddo enddo if (CS%id_PE_to_KE > 0) call post_data(CS%id_PE_to_KE, KE_term, CS%diag) @@ -1037,7 +1037,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie KE_term(i,j,k) = 0.5 * G%IareaT(i,j) & - * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) + * ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1))) enddo ; enddo enddo call post_data(CS%id_KE_BT, KE_term, CS%diag) @@ -1056,13 +1056,13 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS enddo ; enddo do j=js,je ; do i=is,ie KE_h(i,j) = -KE(i,j,k) * G%IareaT(i,j) & - * (uh(I,j,k) - uh(I-1,j,k) + vh(i,J,k) - vh(i,J-1,k)) + * ((uh(I,j,k) - uh(I-1,j,k)) + (vh(i,J,k) - vh(i,J-1,k))) enddo ; enddo if (.not.G%symmetric) & call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie KE_term(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) & - * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) + * ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1))) enddo ; enddo enddo call post_data(CS%id_KE_Coradv, KE_term, CS%diag) @@ -1085,13 +1085,13 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS enddo ; enddo do j=js,je ; do i=is,ie KE_h(i,j) = -KE(i,j,k) * G%IareaT(i,j) & - * (uh(I,j,k) - uh(I-1,j,k) + vh(i,J,k) - vh(i,J-1,k)) + * ((uh(I,j,k) - uh(I-1,j,k)) + (vh(i,J,k) - vh(i,J-1,k))) enddo ; enddo if (.not.G%symmetric) & call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie KE_term(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) & - * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) + * ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1))) enddo ; enddo enddo call post_data(CS%id_KE_adv, KE_term, CS%diag) @@ -1110,7 +1110,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie KE_term(i,j,k) = 0.5 * G%IareaT(i,j) & - * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) + * ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1))) enddo ; enddo enddo call post_data(CS%id_KE_visc, KE_term, CS%diag) @@ -1167,7 +1167,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie KE_term(i,j,k) = 0.5 * G%IareaT(i,j) & - * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) + * ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1))) enddo ; enddo enddo call post_data(CS%id_KE_horvisc, KE_term, CS%diag) @@ -1189,7 +1189,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie KE_term(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) & - * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) + * ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1))) enddo ; enddo enddo call post_data(CS%id_KE_dia, KE_term, CS%diag) diff --git a/src/parameterizations/lateral/MOM_Zanna_Bolton.F90 b/src/parameterizations/lateral/MOM_Zanna_Bolton.F90 index b49d123377..f472118e7d 100644 --- a/src/parameterizations/lateral/MOM_Zanna_Bolton.F90 +++ b/src/parameterizations/lateral/MOM_Zanna_Bolton.F90 @@ -1083,7 +1083,7 @@ subroutine compute_energy_source(u, v, h, fx, fy, G, GV, CS) call do_group_pass(pass_KE_uv, G%domain) do j=js,je ; do i=is,ie KE_term(i,j,k) = 0.5 * G%IareaT(i,j) & - * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) + * ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1))) enddo ; enddo enddo @@ -1096,4 +1096,4 @@ subroutine compute_energy_source(u, v, h, fx, fy, G, GV, CS) end subroutine compute_energy_source -end module MOM_Zanna_Bolton \ No newline at end of file +end module MOM_Zanna_Bolton diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index c510acdf0d..178e6f76e2 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -82,7 +82,14 @@ module MOM_thickness_diffuse !! used for MEKE [H ~> m or kg m-2]. When the total depth is less !! than this, the diffusivity is scaled away. logical :: GM_src_alt !< If true, use the GM energy conversion form S^2*N^2*kappa rather - !! than the streamfunction for the GM source term. + !! than the streamfunction for the GM source term for MEKE. + integer :: MEKE_src_answer_date !< The vintage of the expressions in the GM energy conversion + !! calculation when MEKE_GM_SRC_ALT is true. Values below 20240601 + !! recover the answers from the original implementation, while higher + !! values use expressions that satisfy rotational symmetry. + logical :: MEKE_src_slope_bug !< If true, use a bug that limits the positive values, but not the + !! negative values, of the slopes used when MEKE_GM_SRC_ALT is true. + !! When this is true, it breaks rotational symmetry. logical :: use_GM_work_bug !< If true, use the incorrect sign for the !! top-level work tendency on the top layer. real :: Stanley_det_coeff !< The coefficient correlating SGS temperature variance with the mean @@ -635,12 +642,12 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! interface of a layer that is within a layer [nondim]. 0 m] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: & - Slope_y_PE, & ! 3D array of neutral slopes at v-points, set equal to Slope (below) [nondim] + Slope_y_PE, & ! 3D array of neutral slopes at v-points, set equal to Slope (below) [Z L-1 ~> nondim] hN2_y_PE ! Harmonic mean of thicknesses around the interfaces times the buoyancy frequency ! at v-points with unit conversion factors [H L2 Z-2 T-2 ~> m s-2 or kg m-2 s-2], ! used for calculating the potential energy release real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: & - Slope_x_PE, & ! 3D array of neutral slopes at u-points, set equal to Slope (below) [nondim] + Slope_x_PE, & ! 3D array of neutral slopes at u-points, set equal to Slope (below) [Z L-1 ~> nondim] hN2_x_PE ! Harmonic mean of thicknesses around the interfaces times the buoyancy frequency ! at u-points with unit conversion factors [H L2 Z-2 T-2 ~> m s-2 or kg m-2 s-2], ! used for calculating the potential energy release @@ -998,7 +1005,13 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV int_slope_u(I,j,K) * ((e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j)) slope2_Ratio_u(I,K) = (1.0 - int_slope_u(I,j,K)) * slope2_Ratio_u(I,K) - Slope_x_PE(I,j,k) = MIN(Slope,CS%slope_max) + if (CS%MEKE_src_slope_bug) then + Slope_x_PE(I,j,k) = MIN(Slope, CS%slope_max) + else + Slope_x_PE(I,j,k) = Slope + if (Slope > CS%slope_max) Slope_x_PE(I,j,k) = CS%slope_max + if (Slope < -CS%slope_max) Slope_x_PE(I,j,k) = -CS%slope_max + endif if (CS%id_slope_x > 0) CS%diagSlopeX(I,j,k) = Slope ! Estimate the streamfunction at each interface [H L2 T-1 ~> m3 s-1 or kg s-1]. @@ -1313,7 +1326,13 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV int_slope_v(i,J,K) * ((e(i,j+1,K)-e(i,j,K)) * G%IdyCv(i,J)) slope2_Ratio_v(i,K) = (1.0 - int_slope_v(i,J,K)) * slope2_Ratio_v(i,K) - Slope_y_PE(i,J,k) = MIN(Slope,CS%slope_max) + if (CS%MEKE_src_slope_bug) then + Slope_y_PE(i,J,k) = MIN(Slope, CS%slope_max) + else + Slope_y_PE(i,J,k) = Slope + if (Slope > CS%slope_max) Slope_y_PE(i,J,k) = CS%slope_max + if (Slope < -CS%slope_max) Slope_y_PE(i,J,k) = -CS%slope_max + endif if (CS%id_slope_y > 0) CS%diagSlopeY(I,j,k) = Slope Sfn_unlim_v(i,K) = -((KH_v(i,J,K)*G%dx_Cv(i,J))*Slope) @@ -1550,14 +1569,33 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV enddo ; enddo ; endif if (find_work .and. CS%GM_src_alt) then ; if (allocated(MEKE%GM_src)) then - do j=js,je ; do i=is,ie ; do k=nz,1,-1 - PE_release_h = -0.25 * (GV%H_to_RZ*US%L_to_Z**2) * & + if (CS%MEKE_src_answer_date >= 20240601) then + do j=js,je ; do i=is,ie ; do k=nz,1,-1 + PE_release_h = -0.25 * GV%H_to_RZ * & + ( (KH_u(I,j,k)*(Slope_x_PE(I,j,k)**2) * hN2_x_PE(I,j,k) + & + Kh_u(I-1,j,k)*(Slope_x_PE(I-1,j,k)**2) * hN2_x_PE(I-1,j,k)) + & + (Kh_v(i,J,k)*(Slope_y_PE(i,J,k)**2) * hN2_y_PE(i,J,k) + & + Kh_v(i,J-1,k)*(Slope_y_PE(i,J-1,k)**2) * hN2_y_PE(i,J-1,k)) ) + MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + PE_release_h + enddo ; enddo ; enddo + else + do j=js,je ; do i=is,ie ; do k=nz,1,-1 + PE_release_h = -0.25 * GV%H_to_RZ * & (KH_u(I,j,k)*(Slope_x_PE(I,j,k)**2) * hN2_x_PE(I,j,k) + & Kh_u(I-1,j,k)*(Slope_x_PE(I-1,j,k)**2) * hN2_x_PE(I-1,j,k) + & Kh_v(i,J,k)*(Slope_y_PE(i,J,k)**2) * hN2_y_PE(i,J,k) + & Kh_v(i,J-1,k)*(Slope_y_PE(i,J-1,k)**2) * hN2_y_PE(i,J-1,k)) - MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + PE_release_h - enddo ; enddo ; enddo + MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + PE_release_h + enddo ; enddo ; enddo + endif + if (CS%debug) then + call hchksum(MEKE%GM_src, 'MEKE%GM_src', G%HI, scale=US%RZ3_T3_to_W_m2*US%L_to_Z**2) + call uvchksum("KH_[uv]", Kh_u, Kh_v, G%HI, scale=US%L_to_m**2*US%s_to_T, & + scalar_pair=.true.) + call uvchksum("Slope_[xy]_PE", Slope_x_PE, Slope_y_PE, G%HI, scale=US%Z_to_L) + call uvchksum("hN2_[xy]_PE", hN2_x_PE, hN2_y_PE, G%HI, scale=GV%H_to_mks*US%L_to_Z**2*US%s_to_T**2, & + scalar_pair=.true.) + endif endif ; endif if (CS%id_slope_x > 0) call post_data(CS%id_slope_x, CS%diagSlopeX, CS%diag) @@ -2224,9 +2262,25 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) "If true, write out verbose debugging data.", & default=.false., debuggingParam=.true.) + 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, "MEKE_GM_SRC_ALT", CS%GM_src_alt, & "If true, use the GM energy conversion form S^2*N^2*kappa rather "//& "than the streamfunction for the GM source term.", default=.false.) + call get_param(param_file, mdl, "MEKE_GM_SRC_ANSWER_DATE", CS%MEKE_src_answer_date, & + "The vintage of the expressions in the GM energy conversion calculation when "//& + "MEKE_GM_SRC_ALT is true. Values below 20240601 recover the answers from the "//& + "original implementation, while higher values use expressions that satisfy "//& + "rotational symmetry.", & + default=20240101, do_not_log=.not.CS%GM_src_alt) ! ### Change default to default_answer_date. + call get_param(param_file, mdl, "MEKE_GM_SRC_ALT_SLOPE_BUG", CS%MEKE_src_slope_bug, & + "If true, use a bug that limits the positive values, but not the negative values, "//& + "of the slopes used when MEKE_GM_SRC_ALT is true. When this is true, it breaks "//& + "all of the symmetry rules that MOM6 is supposed to obey.", & + default=.true., do_not_log=.not.CS%GM_src_alt) ! ### Change default to False. + call get_param(param_file, mdl, "MEKE_GEOMETRIC", CS%MEKE_GEOMETRIC, & "If true, uses the GM coefficient formulation from the GEOMETRIC "//& "framework (Marshall et al., 2012).", default=.false.) @@ -2238,9 +2292,6 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) "The nondimensional coefficient governing the efficiency of the GEOMETRIC "//& "thickness diffusion.", units="nondim", default=0.05) - 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) call get_param(param_file, mdl, "MEKE_GEOMETRIC_ANSWER_DATE", CS%MEKE_GEOM_answer_date, & "The vintage of the expressions in the MEKE_GEOMETRIC calculation. "//& "Values below 20190101 recover the answers from the original implementation, "//& diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index 8e95edd563..f480c655d7 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -119,6 +119,10 @@ module MOM_CVMix_KPP !! enhancement with KPP [Z ~> m] logical :: STOKES_MIXING !< Flag if model is mixing down Stokes gradient !! This is relevant for which current to use in RiB + integer :: answer_date !< The vintage of the order of arithmetic in the CVMix KPP + !! calculations. Values below 20240501 recover the answers + !! from early in 2024, while higher values use expressions + !! that have been refactored for rotational symmetry. !> CVMix parameters type(CVMix_kpp_params_type), pointer :: KPP_params => NULL() @@ -199,6 +203,7 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) character(len=20) :: langmuir_mixing_opt = 'NONE' !< Langmuir mixing option to be passed to CVMix, e.g., LWF16 character(len=20) :: langmuir_entrainment_opt = 'NONE' !< Langmuir entrainment option to be !! passed to CVMix, e.g., LWF16 + integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. logical :: CS_IS_ONE=.false. !< Logical for setting Cs based on Non-local logical :: lnoDGat1=.false. !< True => G'(1) = 0 (shape function) !! False => compute G'(1) as in LMD94 @@ -217,6 +222,10 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) if (.not. KPP_init) return allocate(CS) + call get_param(paramFile, 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 openParameterBlock(paramFile,'KPP') call get_param(paramFile, mdl, 'PASSIVE', CS%passiveMode, & 'If True, puts KPP into a passive-diagnostic mode.', & @@ -471,6 +480,12 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) units="m", default=1.0, scale=US%m_to_Z) endif + call get_param(paramFile, mdl, "ANSWER_DATE", CS%answer_date, & + "The vintage of the order of arithmetic in the CVMix KPP calculations. Values "//& + "below 20240501 recover the answers from early in 2024, while higher values "//& + "use expressions that have been refactored for rotational symmetry.", & + default=20240101) !### Change to: default=default_answer_date) + call closeParameterBlock(paramFile) call get_param(paramFile, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) @@ -1390,11 +1405,17 @@ subroutine KPP_smooth_BLD(CS, G, GV, US, dz) wn = 0.125 * G%mask2dT(i,j+1) wc = 1.0 - (ww+we+wn+ws) - CS%OBLdepth(i,j) = wc * OBLdepth_prev(i,j) & - + ww * OBLdepth_prev(i-1,j) & - + we * OBLdepth_prev(i+1,j) & - + ws * OBLdepth_prev(i,j-1) & - + wn * OBLdepth_prev(i,j+1) + if (CS%answer_date < 20240501) then + CS%OBLdepth(i,j) = wc * OBLdepth_prev(i,j) & + + ww * OBLdepth_prev(i-1,j) & + + we * OBLdepth_prev(i+1,j) & + + ws * OBLdepth_prev(i,j-1) & + + wn * OBLdepth_prev(i,j+1) + else + CS%OBLdepth(i,j) = wc * OBLdepth_prev(i,j) & + + ((ww * OBLdepth_prev(i-1,j) + we * OBLdepth_prev(i+1,j)) & + + (ws * OBLdepth_prev(i,j-1) + wn * OBLdepth_prev(i,j+1))) + endif ! Apply OBLdepth smoothing at a cell only if the OBLdepth gets deeper via smoothing. if (CS%deepen_only) CS%OBLdepth(i,j) = max(CS%OBLdepth(i,j), OBLdepth_prev(i,j)) diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 index ee04f3f195..4c0baf3d26 100644 --- a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 +++ b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 @@ -59,6 +59,8 @@ module MOM_bkgnd_mixing real :: N0_2Omega !< ratio of the typical Buoyancy frequency to !! twice the Earth's rotation period, used with the !! Henyey scaling from the mixing [nondim] + real :: Henyey_max_lat !< A latitude poleward of which the Henyey profile + !! is returned to the minimum diffusivity [degN] real :: prandtl_bkgnd !< Turbulent Prandtl number used to convert !! vertical background diffusivity into viscosity [nondim] real :: Kd_tanh_lat_scale !< A nondimensional scaling for the range of @@ -282,6 +284,10 @@ subroutine bkgnd_mixing_init(Time, G, GV, US, param_file, diag, CS, physical_OBL call get_param(param_file, mdl, "OMEGA", CS%omega, & "The rotation rate of the earth.", & units="s-1", default=7.2921e-5, scale=US%T_to_s) + call get_param(param_file, mdl, "HENYEY_MAX_LAT", CS%Henyey_max_lat, & + "A latitude poleward of which the Henyey profile "//& + "is returned to the minimum diffusivity", & + units="degN", default=95.0) endif call get_param(param_file, mdl, "KD_TANH_LAT_FN", CS%Kd_tanh_lat_fn, & @@ -447,6 +453,7 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd_lay, Kd_int, Kv_bkgnd, j, G, I_x30 = 2.0 / invcosh(CS%N0_2Omega*2.0) ! This is evaluated at 30 deg. do i=is,ie abs_sinlat = abs(sin(G%geoLatT(i,j)*deg_to_rad)) + if (abs(G%geoLatT(i,j))>CS%Henyey_max_lat) abs_sinlat = min_sinlat Kd_sfc(i) = max(CS%Kd_min, CS%Kd * & ((abs_sinlat * invcosh(CS%N0_2Omega / max(min_sinlat, abs_sinlat))) * I_x30) ) enddo diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index e601fdf2f7..c103cbc71b 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -59,6 +59,7 @@ module MOM_set_visc real :: drag_bg_vel !< An assumed unresolved background velocity for !! calculating the bottom drag [L T-1 ~> m s-1]. !! Runtime parameter `DRAG_BG_VEL`. + !! Should not be used if BBL_USE_TIDAL_BG is True. real :: BBL_thick_min !< The minimum bottom boundary layer thickness [Z ~> m]. !! This might be Kv / (cdrag * drag_bg_vel) to give !! Kv as the minimum near-bottom viscosity. @@ -229,11 +230,8 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) real :: BBL_thick_max ! A huge upper bound on the boundary layer thickness [Z ~> m]. real :: kv_bbl ! The bottom boundary layer viscosity [H Z T-1 ~> m2 s-1 or Pa s] real :: C2f ! C2f = 2*f at velocity points [T-1 ~> s-1]. - - real :: U_bg_sq ! The square of an assumed background - ! velocity, for calculating the mean - ! magnitude near the bottom for use in the - ! quadratic bottom drag [L2 T-2 ~> m2 s-2]. + real :: u2_bg(SZIB_(G)) ! The square of an assumed background velocity, for calculating the mean + ! magnitude near the bottom for use in the quadratic bottom drag [L2 T-2 ~> m2 s-2]. real :: hwtot ! Sum of the thicknesses used to calculate ! the near-bottom velocity magnitude [H ~> m or kg m-2]. real :: I_hwtot ! The Adcroft reciprocal of hwtot [H-1 ~> m-1 or m2 kg-1]. @@ -338,7 +336,6 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) use_BBL_EOS = associated(tv%eqn_of_state) .and. CS%BBL_use_EOS OBC => CS%OBC - U_bg_sq = CS%drag_bg_vel * CS%drag_bg_vel cdrag_sqrt = sqrt(CS%cdrag) cdrag_sqrt_H = cdrag_sqrt * US%L_to_m * GV%m_to_H cdrag_sqrt_H_RL = cdrag_sqrt * US%L_to_Z * GV%RZ_to_H @@ -422,7 +419,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) !$OMP parallel do default(private) shared(u,v,h,dz,tv,visc,G,GV,US,CS,Rml,nz,nkmb,nkml,K2, & !$OMP Isq,Ieq,Jsq,Jeq,h_neglect,dz_neglect,Rho0x400_G, & - !$OMP U_bg_sq,cdrag_sqrt,cdrag_sqrt_H,cdrag_sqrt_H_RL, & + !$OMP cdrag_sqrt,cdrag_sqrt_H,cdrag_sqrt_H_RL, & !$OMP cdrag_L_to_H,cdrag_RL_to_H,use_BBL_EOS,BBL_thick_max, & !$OMP OBC,D_u,D_v,mask_u,mask_v,pbv) do j=Jsq,Jeq ; do m=1,2 @@ -591,6 +588,19 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0 dztot_vel = 0.0 ; dzwtot = 0.0 Thtot = 0.0 ; Shtot = 0.0 ; SpV_htot = 0.0 + + ! Set the "back ground" friction velocity scale to either the tidal amplitude or place-holder constant + if (CS%BBL_use_tidal_bg) then + if (m==1) then + u2_bg(I) = 0.5*( G%mask2dT(i,j)*(CS%tideamp(i,j)*CS%tideamp(i,j))+ & + G%mask2dT(i+1,j)*(CS%tideamp(i+1,j)*CS%tideamp(i+1,j)) ) + else + u2_bg(i) = 0.5*( G%mask2dT(i,j)*(CS%tideamp(i,j)*CS%tideamp(i,j))+ & + G%mask2dT(i,j+1)*(CS%tideamp(i,j+1)*CS%tideamp(i,j+1)) ) + endif + else + u2_bg(i) = CS%drag_bg_vel * CS%drag_bg_vel + endif do k=nz,1,-1 if (htot_vel>=CS%Hbbl) exit ! terminate the k loop @@ -606,18 +616,10 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) if ((.not.CS%linear_drag) .and. (hweight >= 0.0)) then ; if (m==1) then v_at_u = set_v_at_u(v, h, G, GV, i, j, k, mask_v, OBC) - if (CS%BBL_use_tidal_bg) then - U_bg_sq = 0.5*( G%mask2dT(i,j)*(CS%tideamp(i,j)*CS%tideamp(i,j))+ & - G%mask2dT(i+1,j)*(CS%tideamp(i+1,j)*CS%tideamp(i+1,j)) ) - endif - hutot = hutot + hweight * sqrt(u(I,j,k)*u(I,j,k) + v_at_u*v_at_u + U_bg_sq) + hutot = hutot + hweight * sqrt(u(I,j,k)*u(I,j,k) + v_at_u*v_at_u + u2_bg(I)) else u_at_v = set_u_at_v(u, h, G, GV, i, j, k, mask_u, OBC) - if (CS%BBL_use_tidal_bg) then - U_bg_sq = 0.5*( G%mask2dT(i,j)*(CS%tideamp(i,j)*CS%tideamp(i,j))+ & - G%mask2dT(i,j+1)*(CS%tideamp(i,j+1)*CS%tideamp(i,j+1)) ) - endif - hutot = hutot + hweight * sqrt(v(i,J,k)*v(i,J,k) + u_at_v*u_at_v + U_bg_sq) + hutot = hutot + hweight * sqrt(v(i,J,k)*v(i,J,k) + u_at_v*u_at_v + u2_bg(i)) endif ; endif if (use_BBL_EOS .and. (hweight >= 0.0)) then @@ -798,6 +800,19 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) if (m==1) then ; C2f = G%CoriolisBu(I,J-1) + G%CoriolisBu(I,J) else ; C2f = G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J) ; endif + ! Set the "back ground" friction velocity scale to either the tidal amplitude or place-holder constant + if (CS%BBL_use_tidal_bg) then + if (m==1) then + u2_bg(I) = 0.5*( G%mask2dT(i,j)*(CS%tideamp(i,j)*CS%tideamp(i,j))+ & + G%mask2dT(i+1,j)*(CS%tideamp(i+1,j)*CS%tideamp(i+1,j)) ) + else + u2_bg(i) = 0.5*( G%mask2dT(i,j)*(CS%tideamp(i,j)*CS%tideamp(i,j))+ & + G%mask2dT(i,j+1)*(CS%tideamp(i,j+1)*CS%tideamp(i,j+1)) ) + endif + else + u2_bg(i) = CS%drag_bg_vel * CS%drag_bg_vel + endif + ! The thickness of a rotation limited BBL ignoring stratification is ! h_f ~ Cn u* / f (limit of KW99 eq. 2.20 for N->0). ! The buoyancy limit of BBL thickness (h_N) is already in the variable htot from above. @@ -809,7 +824,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) ! xp = 1/2 + sqrt( 1/4 + (2 f h_N/u*)^2 ) ! To avoid dividing by zero if u*=0 then ! xp u* = 1/2 u* + sqrt( 1/4 u*^2 + (2 f h_N)^2 ) - if (CS%cdrag * U_bg_sq <= 0.0) then + if (CS%cdrag * u2_bg(i) <= 0.0) then ! This avoids NaNs and overflows, and could be used in all cases, ! but is not bitwise identical to the current code. ustH = ustar(i) ; root = sqrt(0.25*ustH**2 + (htot*C2f)**2) @@ -957,12 +972,12 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) if (m==1) then if (Rayleigh > 0.0) then v_at_u = set_v_at_u(v, h, G, GV, i, j, k, mask_v, OBC) - visc%Ray_u(I,j,k) = Rayleigh * sqrt(u(I,j,k)*u(I,j,k) + v_at_u*v_at_u + U_bg_sq) + visc%Ray_u(I,j,k) = Rayleigh * sqrt(u(I,j,k)*u(I,j,k) + v_at_u*v_at_u + u2_bg(I)) else ; visc%Ray_u(I,j,k) = 0.0 ; endif else if (Rayleigh > 0.0) then u_at_v = set_u_at_v(u, h, G, GV, i, j, k, mask_u, OBC) - visc%Ray_v(i,J,k) = Rayleigh * sqrt(v(i,J,k)*v(i,J,k) + u_at_v*u_at_v + U_bg_sq) + visc%Ray_v(i,J,k) = Rayleigh * sqrt(v(i,J,k)*v(i,J,k) + u_at_v*u_at_v + u2_bg(i)) else ; visc%Ray_v(i,J,k) = 0.0 ; endif endif @@ -1992,9 +2007,9 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS) real :: frac_used ! The fraction of the present layer that contributes to Dh and Ddz [nondim] real :: Dh ! The increment in layer thickness from the present layer [H ~> m or kg m-2]. real :: Ddz ! The increment in height change from the present layer [Z ~> m]. - real :: U_bg_sq ! The square of an assumed background velocity, for - ! calculating the mean magnitude near the top for use in - ! the quadratic surface drag [L2 T-2 ~> m2 s-2]. + real :: u2_bg(SZIB_(G)) ! The square of an assumed background velocity, for + ! calculating the mean magnitude near the top for use in + ! the quadratic surface drag [L2 T-2 ~> m2 s-2]. real :: h_tiny ! A very small thickness [H ~> m or kg m-2]. Layers that are less than ! h_tiny can not be the deepest in the viscous mixed layer. real :: absf ! The absolute value of f averaged to velocity points [T-1 ~> s-1]. @@ -2025,7 +2040,6 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS) associated(forces%frac_shelf_v)) ) return Rho0x400_G = 400.0*(GV%H_to_RZ / (US%L_to_Z**2 * GV%g_Earth)) - U_bg_sq = CS%drag_bg_vel * CS%drag_bg_vel cdrag_sqrt = sqrt(CS%cdrag) cdrag_sqrt_H = cdrag_sqrt * US%L_to_m * GV%m_to_H cdrag_sqrt_H_RL = cdrag_sqrt * US%L_to_Z * GV%RZ_to_H @@ -2099,7 +2113,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS) !$OMP parallel do default(private) shared(u,v,h,dz,tv,forces,visc,dt,G,GV,US,CS,use_EOS,dt_Rho0, & !$OMP nonBous_ML,h_neglect,dz_neglect,h_tiny,g_H_Rho0, & - !$OMP js,je,OBC,Isq,Ieq,nz,nkml,U_star_2d,U_bg_sq,mask_v, & + !$OMP js,je,OBC,Isq,Ieq,nz,nkml,U_star_2d,mask_v, & !$OMP cdrag_sqrt,cdrag_sqrt_H,cdrag_sqrt_H_RL,Rho0x400_G) do j=js,je ! u-point loop if (CS%dynamic_viscous_ML) then @@ -2251,7 +2265,14 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS) if (.not.CS%linear_drag) then v_at_u = set_v_at_u(v, h, G, GV, i, j, k, mask_v, OBC) - hutot = hutot + hweight * sqrt(u(I,j,k)**2 + v_at_u**2 + U_bg_sq) + ! Set the "back ground" friction velocity scale to either the tidal amplitude or place-holder constant + if (CS%BBL_use_tidal_bg) then + u2_bg(I) = 0.5*( G%mask2dT(i,j)*(CS%tideamp(i,j)*CS%tideamp(i,j))+ & + G%mask2dT(i+1,j)*(CS%tideamp(i+1,j)*CS%tideamp(i+1,j)) ) + else + u2_bg(I) = CS%drag_bg_vel * CS%drag_bg_vel + endif + hutot = hutot + hweight * sqrt(u(I,j,k)**2 + v_at_u**2 + u2_bg(I)) endif if (use_EOS) then Thtot(I) = Thtot(I) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k)) @@ -2369,7 +2390,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS) !$OMP parallel do default(private) shared(u,v,h,dz,tv,forces,visc,dt,G,GV,US,CS,use_EOS,dt_Rho0, & !$OMP nonBous_ML,h_neglect,dz_neglect,h_tiny,g_H_Rho0, & - !$OMP is,ie,OBC,Jsq,Jeq,nz,nkml,U_bg_sq,U_star_2d,mask_u, & + !$OMP is,ie,OBC,Jsq,Jeq,nz,nkml,U_star_2d,mask_u, & !$OMP cdrag_sqrt,cdrag_sqrt_H,cdrag_sqrt_H_RL,Rho0x400_G) do J=Jsq,Jeq ! v-point loop if (CS%dynamic_viscous_ML) then @@ -2523,7 +2544,14 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS) if (.not.CS%linear_drag) then u_at_v = set_u_at_v(u, h, G, GV, i, J, k, mask_u, OBC) - hutot = hutot + hweight * sqrt(v(i,J,k)**2 + u_at_v**2 + U_bg_sq) + ! Set the "back ground" friction velocity scale to either the tidal amplitude or place-holder constant + if (CS%BBL_use_tidal_bg) then + u2_bg(i) = 0.5*( G%mask2dT(i,j)*(CS%tideamp(i,j)*CS%tideamp(i,j))+ & + G%mask2dT(i,j+1)*(CS%tideamp(i,j+1)*CS%tideamp(i,j+1)) ) + else + u2_bg(i) = CS%drag_bg_vel * CS%drag_bg_vel + endif + hutot = hutot + hweight * sqrt(v(i,J,k)**2 + u_at_v**2 + u2_bg(i)) endif if (use_EOS) then Thtot(i) = Thtot(i) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k)) @@ -2967,6 +2995,11 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS call get_param(param_file, mdl, "TIDEAMP_VARNAME", tideamp_var, & "The name of the tidal amplitude variable in the input file.", & default="tideamp") + ! This value is here only to detect whether it is inadvertently used. CS%drag_bg_vel should + ! not be used if CS%BBL_use_tidal_bg is True. For this reason, we do not apply dimensions, + ! nor dimensional testing in this mode. If we ever detect a dimensional sensitivity to + ! this parameter, in this mode, then it means it is being used inappropriately. + CS%drag_bg_vel = 1.e30 else call get_param(param_file, mdl, "DRAG_BG_VEL", CS%drag_bg_vel, & "DRAG_BG_VEL is either the assumed bottom velocity (with "//& diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 824ae22016..c26ee4ac75 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -234,11 +234,15 @@ subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OB real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: omega_tau2w_u !< angle between mtm flux and wind at u-pts [rad] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: omega_tau2w_v !< angle between mtm flux and wind at v-pts [rad] - real :: pi, Cemp_CG, tmp, cos_tmp, sin_tmp, omega_tmp !< constants and dummy variables - real :: du, dv, depth, sigma, Wind_x, Wind_y !< intermediate variables - real :: taux, tauy, tauxDG, tauyDG, tauxDGup, tauyDGup, ustar2, tauh !< intermediate variables - real :: tauNLup, tauNLdn, tauNL_CG, tauNL_DG, tauNL_X, tauNL_Y, tau_MAG !< intermediate variables - real :: omega_w2s, omega_tau2s, omega_s2x, omega_tau2x, omega_tau2w, omega_s2w !< intermediate angles + real :: pi, Cemp_CG, tmp, cos_tmp, sin_tmp !< constants and dummy variables [nondim] + real :: omega_tmp !< A dummy angle [radians] + real :: du, dv !< Velocity increments [L T-1 ~> m s-1] + real :: depth !< Cumulative layer thicknesses [H ~> m or kg m=2] + real :: sigma !< Fractional depth in the mixed layer [nondim] + real :: Wind_x, Wind_y !< intermediate wind stress componenents [L2 T-2 ~> m2 s-2] + real :: taux, tauy, tauxDG, tauyDG, tauxDGup, tauyDGup, ustar2, tauh !< intermediate variables [L2 T-2 ~> m2 s-2] + real :: tauNLup, tauNLdn, tauNL_CG, tauNL_DG, tauNL_X, tauNL_Y, tau_MAG !< intermediate variables [L2 T-2 ~> m2 s-2] + real :: omega_w2s, omega_tau2s, omega_s2x, omega_tau2x, omega_tau2w, omega_s2w !< intermediate angles [radians] integer :: kblmin, kbld, kp1, k, nz !< vertical indices integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq ! horizontal indices @@ -321,6 +325,7 @@ subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OB enddo if (CS%debug) then + !### These checksum calls are missing necessary dimensional scaling factors. call uvchksum("surface tau[xy]_[uv] ", taux_u, tauy_v, G%HI, haloshift=1, scalar_pair=.true.) call uvchksum("ustar2", ustar2_u, ustar2_v, G%HI, haloshift=0, scalar_pair=.true.) call uvchksum(" hbl", hbl_u , hbl_v , G%HI, haloshift=0, scalar_pair=.true.) @@ -427,6 +432,7 @@ subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OB kbld = min( (kbl_u(I,j)) , (nz-2) ) if ( tau_u(I,j,kbld+2) > tau_u(I,j,kbld+1) ) kbld = kbld + 1 + !### This expression is dimensionally inconsistent. tauh = tau_u(I,j,kbld+1) + GV%H_subroundoff ! surface boundary conditions depth = 0. @@ -437,6 +443,7 @@ subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OB ! linear stress mag tau_MAG = (ustar2_u(I,j) * (1.-sigma) ) + (tauh * sigma ) + !### The following expressions are dimensionally inconsistent. cos_tmp = tauxDG_u(I,j,k+1) / (tau_u(I,j,k+1) + GV%H_subroundoff) sin_tmp = tauyDG_u(I,j,k+1) / (tau_u(I,j,k+1) + GV%H_subroundoff) @@ -457,6 +464,7 @@ subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OB tauNLdn = tauNL_X ! nonlocal increment and update to uold + !### The following expression is dimensionally inconsistent and missing parentheses. du = (tauNLup - tauNLdn) * (dt/CS%h_u(I,j,k) + GV%H_subroundoff) ui(I,j,k) = uold(I,j,k) + du uold(I,j,k) = du @@ -496,6 +504,7 @@ subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OB ! linear stress tau_MAG = (ustar2_v(i,J) * (1.-sigma)) + (tauh * sigma) + !### The following expressions are dimensionally inconsistent. cos_tmp = tauxDG_v(i,J,k+1) / (tau_v(i,J,k+1) + GV%H_subroundoff) sin_tmp = tauyDG_v(i,J,k+1) / (tau_v(i,J,k+1) + GV%H_subroundoff) @@ -514,6 +523,8 @@ subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OB tauNL_X = (tauNL_DG * cos_tmp - tauNL_CG * sin_tmp) tauNL_Y = (tauNL_DG * sin_tmp + tauNL_CG * cos_tmp) tauNLdn = tauNL_Y + !### The following expression is dimensionally inconsistent, [L T-1] vs. [L2 H-1 T-1] on the right, + ! and it is inconsistent with the counterpart expression for du. dv = (tauNLup - tauNLdn) * (dt/(CS%h_v(i,J,k)) ) vi(i,J,k) = vold(i,J,k) + dv vold(i,J,k) = dv @@ -2634,7 +2645,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & # include "version_variable.h" character(len=40) :: mdl = "MOM_vert_friction" ! This module's name. character(len=40) :: thickness_units - real :: Kv_mks ! KVML in MKS + real :: Kv_mks ! KVML in MKS [m2 s-1] if (associated(CS)) then call MOM_error(WARNING, "vertvisc_init called with an associated "// & diff --git a/src/user/Neverworld_initialization.F90 b/src/user/Neverworld_initialization.F90 index 05de663d46..6885b6881a 100644 --- a/src/user/Neverworld_initialization.F90 +++ b/src/user/Neverworld_initialization.F90 @@ -34,12 +34,12 @@ module Neverworld_initialization subroutine Neverworld_initialize_topography(D, G, param_file, max_depth) type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type real, dimension(G%isd:G%ied,G%jsd:G%jed), & - intent(out) :: D !< Ocean bottom depth in the units of depth_max + intent(out) :: D !< Ocean bottom depth in the units of depth_max [A] type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum ocean depth in arbitrary units + real, intent(in) :: max_depth !< Maximum ocean depth in arbitrary units [A] ! Local variables - real :: PI ! 3.1415926... calculated as 4*atan(1) + real :: PI ! 3.1415926... calculated as 4*atan(1) [nondim] real :: x, y ! Lateral positions normalized by the domain size [nondim] ! This include declares and sets the variable "version". # include "version_variable.h" @@ -84,9 +84,9 @@ end subroutine Neverworld_initialize_topography !> Returns the value of a cosine-bell function evaluated at x/L real function cosbell(x, L) - real , intent(in) :: x !< non-dimensional position [nondim] - real , intent(in) :: L !< non-dimensional width [nondim] - real :: PI !< 3.1415926... calculated as 4*atan(1) + real , intent(in) :: x !< Position in arbitrary units [A] + real , intent(in) :: L !< Width in arbitrary units [A] + real :: PI !< 3.1415926... calculated as 4*atan(1) [nondim] PI = 4.0*atan(1.0) cosbell = 0.5 * (1 + cos(PI*MIN(ABS(x/L),1.0))) @@ -95,9 +95,9 @@ end function cosbell !> Returns the value of a sin-spike function evaluated at x/L real function spike(x, L) - real , intent(in) :: x !< non-dimensional position [nondim] - real , intent(in) :: L !< non-dimensional width [nondim] - real :: PI !< 3.1415926... calculated as 4*atan(1) + real , intent(in) :: x !< Position in arbitrary units [A] + real , intent(in) :: L !< Width in arbitrary units [A] + real :: PI !< 3.1415926... calculated as 4*atan(1) [nondim] PI = 4.0*atan(1.0) spike = (1 - sin(PI*MIN(ABS(x/L),0.5))) @@ -108,9 +108,9 @@ end function spike !! If clip is present the top of the cone is cut off at "clip", which !! effectively defaults to 1. real function cone(x, x0, L, clip) - real, intent(in) :: x !< non-dimensional coordinate [nondim] - real, intent(in) :: x0 !< position of peak [nondim] - real, intent(in) :: L !< half-width of base of cone [nondim] + real, intent(in) :: x !< Coordinate in arbitrary units [A] + real, intent(in) :: x0 !< position of peak in arbitrary units [A] + real, intent(in) :: L !< half-width of base of cone in arbitrary units [A] real, optional, intent(in) :: clip !< clipping height of cone [nondim] cone = max( 0., 1. - abs(x - x0) / L ) @@ -119,10 +119,10 @@ end function cone !> Returns an s-curve s(x) s.t. s(x0)<=0, s(x0+L)>=1 and cubic in between. real function scurve(x, x0, L) - real, intent(in) :: x !< non-dimensional coordinate [nondim] - real, intent(in) :: x0 !< position of peak [nondim] - real, intent(in) :: L !< half-width of base of cone [nondim] - real :: s + real, intent(in) :: x !< Coordinate in arbitrary units [A] + real, intent(in) :: x0 !< position of peak in arbitrary units [A] + real, intent(in) :: L !< half-width of base of cone in arbitrary units [A] + real :: s ! A rescaled position [nondim] s = max( 0., min( 1.,( x - x0 ) / L ) ) scurve = ( 3. - 2.*s ) * ( s * s ) @@ -132,14 +132,14 @@ end function scurve !> Returns a "coastal" profile. real function cstprof(x, x0, L, lf, bf, sf, sh) - real, intent(in) :: x !< non-dimensional coordinate [nondim] - real, intent(in) :: x0 !< position of peak [nondim] - real, intent(in) :: L !< width of profile [nondim] + real, intent(in) :: x !< Coordinate in arbitrary units [A] + real, intent(in) :: x0 !< position of peak in arbitrary units [A] + real, intent(in) :: L !< width of profile in arbitrary units [A] real, intent(in) :: lf !< fraction of width that is "land" [nondim] real, intent(in) :: bf !< fraction of width that is "beach" [nondim] real, intent(in) :: sf !< fraction of width that is "continental slope" [nondim] real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim] - real :: s + real :: s ! A rescaled position [nondim] s = max( 0., min( 1.,( x - x0 ) / L ) ) cstprof = sh * scurve(s-lf,0.,bf) + (1.-sh) * scurve(s - (1.-sf),0.,sf) @@ -147,12 +147,12 @@ end function cstprof !> Distance between points x,y and a line segment (x0,y0) and (x0,y1). real function dist_line_fixed_x(x, y, x0, y0, y1) - real, intent(in) :: x !< non-dimensional x-coordinate [nondim] - real, intent(in) :: y !< non-dimensional y-coordinate [nondim] - real, intent(in) :: x0 !< x-position of line segment [nondim] - real, intent(in) :: y0 !< y-position of line segment end[nondim] - real, intent(in) :: y1 !< y-position of line segment end[nondim] - real :: dx, yr, dy + real, intent(in) :: x !< X-coordinate in arbitrary units [A] + real, intent(in) :: y !< Y-coordinate in arbitrary units [A] + real, intent(in) :: x0 !< x-position of line segment in arbitrary units [A] + real, intent(in) :: y0 !< y-position of line segment end in arbitrary units [A] + real, intent(in) :: y1 !< y-position of line segment end in arbitrary units [A] + real :: dx, yr, dy ! Relative positions in arbitrary units [A] dx = x - x0 yr = min( max(y0,y1), max( min(y0,y1), y ) ) ! bound y by y0,y1 @@ -162,11 +162,11 @@ end function dist_line_fixed_x !> Distance between points x,y and a line segment (x0,y0) and (x1,y0). real function dist_line_fixed_y(x, y, x0, x1, y0) - real, intent(in) :: x !< non-dimensional x-coordinate [nondim] - real, intent(in) :: y !< non-dimensional y-coordinate [nondim] - real, intent(in) :: x0 !< x-position of line segment end[nondim] - real, intent(in) :: x1 !< x-position of line segment end[nondim] - real, intent(in) :: y0 !< y-position of line segment [nondim] + real, intent(in) :: x !< X-coordinate in arbitrary units [A] + real, intent(in) :: y !< Y-coordinate in arbitrary units [A] + real, intent(in) :: x0 !< x-position of line segment end in arbitrary units [A] + real, intent(in) :: x1 !< x-position of line segment end in arbitrary units [A] + real, intent(in) :: y0 !< y-position of line segment in arbitrary units [A] dist_line_fixed_y = dist_line_fixed_x(y, x, y0, x0, x1) end function dist_line_fixed_y @@ -180,7 +180,7 @@ real function NS_coast(lon, lat, lon0, lat0, lat1, dlon, sh) real, intent(in) :: lat1 !< Latitude of coast end [degrees_N] real, intent(in) :: dlon !< "Radius" of coast profile [degrees] real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim] - real :: r + real :: r ! A relative position [nondim] r = dist_line_fixed_x( lon, lat, lon0, lat0, lat1 ) NS_coast = cstprof(r, 0., dlon, 0.125, 0.125, 0.5, sh) @@ -195,7 +195,7 @@ real function EW_coast(lon, lat, lon0, lon1, lat0, dlat, sh) real, intent(in) :: lat0 !< Latitude of coast [degrees_N] real, intent(in) :: dlat !< "Radius" of coast profile [degrees] real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim] - real :: r + real :: r ! A relative position [nondim] r = dist_line_fixed_y( lon, lat, lon0, lon1, lat0 ) EW_coast = cstprof(r, 0., dlat, 0.125, 0.125, 0.5, sh) @@ -210,7 +210,7 @@ real function NS_ridge(lon, lat, lon0, lat0, lat1, dlon, rh) real, intent(in) :: lat1 !< Latitude of ridge end [degrees_N] real, intent(in) :: dlon !< "Radius" of ridge profile [degrees] real, intent(in) :: rh !< depth of ridge as fraction of full depth [nondim] - real :: r + real :: r ! A distance from a point [degrees] r = dist_line_fixed_x( lon, lat, lon0, lat0, lat1 ) NS_ridge = 1. - rh * cone(r, 0., dlon) @@ -226,12 +226,13 @@ real function circ_ridge(lon, lat, lon0, lat0, ring_radius, ring_thickness, ridg real, intent(in) :: ring_radius !< Radius of ring [degrees] real, intent(in) :: ring_thickness !< Radial thickness of ring [degrees] real, intent(in) :: ridge_height !< Ridge height as fraction of full depth [nondim] - real :: r + real :: r ! A relative position [degrees] + real :: frac_ht ! The fractional height of the topography [nondim] r = sqrt( (lon - lon0)**2 + (lat - lat0)**2 ) ! Pseudo-distance from a point r = abs( r - ring_radius) ! Pseudo-distance from a circle - r = cone(r, 0., ring_thickness, ridge_height) ! 0 .. frac_ridge_height - circ_ridge = 1. - r ! Fractional depths (1-frac_ridge_height) .. 1 + frac_ht = cone(r, 0., ring_thickness, ridge_height) ! 0 .. frac_ridge_height + circ_ridge = 1. - frac_ht ! Fractional depths (1-frac_ridge_height) .. 1 end function circ_ridge !> This subroutine initializes layer thicknesses for the Neverworld test case, diff --git a/src/user/basin_builder.F90 b/src/user/basin_builder.F90 index 42083b2672..705925a97d 100644 --- a/src/user/basin_builder.F90 +++ b/src/user/basin_builder.F90 @@ -27,13 +27,13 @@ module basin_builder subroutine basin_builder_topography(D, G, param_file, max_depth) type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type real, dimension(G%isd:G%ied,G%jsd:G%jed), & - intent(out) :: D !< Ocean bottom depth in the units of depth_max + intent(out) :: D !< Ocean bottom depth in the units of depth_max [A] type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum ocean depth in arbitrary units + real, intent(in) :: max_depth !< Maximum ocean depth in arbitrary units [A] ! Local variables character(len=17) :: pname1, pname2 ! For construction of parameter names character(len=20) :: funcs ! Basin build function - real, dimension(20) :: pars ! Parameters for each function + real, dimension(20) :: pars ! Parameters for each function [various] real :: lon ! Longitude [degrees_E] real :: lat ! Latitude [degrees_N] integer :: i, j, n, n_funcs @@ -161,9 +161,9 @@ end subroutine basin_builder_topography !! If clip is present the top of the cone is cut off at "clip", which !! effectively defaults to 1. real function cone(x, x0, L, clip) - real, intent(in) :: x !< non-dimensional coordinate [nondim] - real, intent(in) :: x0 !< position of peak [nondim] - real, intent(in) :: L !< half-width of base of cone [nondim] + real, intent(in) :: x !< Coordinate in arbitrary units [A] + real, intent(in) :: x0 !< position of peak in arbitrary units [A] + real, intent(in) :: L !< half-width of base of cone in arbitrary units [A] real, optional, intent(in) :: clip !< clipping height of cone [nondim] cone = max( 0., 1. - abs(x - x0) / L ) @@ -172,10 +172,10 @@ end function cone !> Returns an s-curve s(x) s.t. s(x0)<=0, s(x0+L)>=1 and cubic in between. real function scurve(x, x0, L) - real, intent(in) :: x !< non-dimensional coordinate [nondim] - real, intent(in) :: x0 !< position of peak [nondim] - real, intent(in) :: L !< half-width of base of cone [nondim] - real :: s + real, intent(in) :: x !< Coordinate in arbitrary units [A] + real, intent(in) :: x0 !< position of peak in arbitrary units [A] + real, intent(in) :: L !< half-width of base of cone in arbitrary units [A] + real :: s ! A rescaled position [nondim] s = max( 0., min( 1.,( x - x0 ) / L ) ) scurve = ( 3. - 2.*s ) * ( s * s ) @@ -183,14 +183,14 @@ end function scurve !> Returns a "coastal" profile. real function cstprof(x, x0, L, lf, bf, sf, sh) - real, intent(in) :: x !< non-dimensional coordinate [nondim] - real, intent(in) :: x0 !< position of peak [nondim] - real, intent(in) :: L !< width of profile [nondim] + real, intent(in) :: x !< Coordinate in arbitrary units [A] + real, intent(in) :: x0 !< position of peak in arbitrary units [A] + real, intent(in) :: L !< width of profile in arbitrary units [A] real, intent(in) :: lf !< fraction of width that is "land" [nondim] real, intent(in) :: bf !< fraction of width that is "beach" [nondim] real, intent(in) :: sf !< fraction of width that is "continental slope" [nondim] real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim] - real :: s + real :: s ! A rescaled position [nondim] s = max( 0., min( 1.,( x - x0 ) / L ) ) cstprof = sh * scurve(s-lf,0.,bf) + (1.-sh) * scurve(s - (1.-sf),0.,sf) @@ -198,12 +198,12 @@ end function cstprof !> Distance between points x,y and a line segment (x0,y0) and (x0,y1). real function dist_line_fixed_x(x, y, x0, y0, y1) - real, intent(in) :: x !< non-dimensional x-coordinate [nondim] - real, intent(in) :: y !< non-dimensional y-coordinate [nondim] - real, intent(in) :: x0 !< x-position of line segment [nondim] - real, intent(in) :: y0 !< y-position of line segment end[nondim] - real, intent(in) :: y1 !< y-position of line segment end[nondim] - real :: dx, yr, dy + real, intent(in) :: x !< X-coordinate in arbitrary units [A] + real, intent(in) :: y !< Y-coordinate in arbitrary units [A] + real, intent(in) :: x0 !< x-position of line segment in arbitrary units [A] + real, intent(in) :: y0 !< y-position of line segment end in arbitrary units [A] + real, intent(in) :: y1 !< y-position of line segment end in arbitrary units [A] + real :: dx, yr, dy ! Relative positions in arbitrary units [A] dx = x - x0 yr = min( max(y0,y1), max( min(y0,y1), y ) ) ! bound y by y0,y1 @@ -213,11 +213,11 @@ end function dist_line_fixed_x !> Distance between points x,y and a line segment (x0,y0) and (x1,y0). real function dist_line_fixed_y(x, y, x0, x1, y0) - real, intent(in) :: x !< non-dimensional x-coordinate [nondim] - real, intent(in) :: y !< non-dimensional y-coordinate [nondim] - real, intent(in) :: x0 !< x-position of line segment end[nondim] - real, intent(in) :: x1 !< x-position of line segment end[nondim] - real, intent(in) :: y0 !< y-position of line segment [nondim] + real, intent(in) :: x !< X-coordinate in arbitrary units [A] + real, intent(in) :: y !< Y-coordinate in arbitrary units [A] + real, intent(in) :: x0 !< x-position of line segment end in arbitrary units [A] + real, intent(in) :: x1 !< x-position of line segment end in arbitrary units [A] + real, intent(in) :: y0 !< y-position of line segment in arbitrary units [A] dist_line_fixed_y = dist_line_fixed_x(y, x, y0, x0, x1) end function dist_line_fixed_y @@ -230,10 +230,11 @@ real function angled_coast(lon, lat, lon_eq, lat_mer, dr, sh) real, intent(in) :: lat_mer !< Latitude intersection with Prime Meridian [degrees_N] real, intent(in) :: dr !< "Radius" of coast profile [degrees] real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim] - real :: r + real :: r ! A relative position [degrees] + real :: I_dr ! The inverse of a distance [degrees-1] - r = 1/sqrt( lat_mer*lat_mer + lon_eq*lon_eq ) - r = r * ( lat_mer*lon + lon_eq*lat - lon_eq*lat_mer) + I_dr = 1/sqrt( lat_mer*lat_mer + lon_eq*lon_eq ) + r = I_dr * ( lat_mer*lon + lon_eq*lat - lon_eq*lat_mer) angled_coast = cstprof(r, 0., dr, 0.125, 0.125, 0.5, sh) end function angled_coast @@ -246,7 +247,7 @@ real function NS_coast(lon, lat, lonC, lat0, lat1, dlon, sh) real, intent(in) :: lat1 !< Latitude of coast end [degrees_N] real, intent(in) :: dlon !< "Radius" of coast profile [degrees] real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim] - real :: r + real :: r ! A relative position [degrees] r = dist_line_fixed_x( lon, lat, lonC, lat0, lat1 ) NS_coast = cstprof(r, 0., dlon, 0.125, 0.125, 0.5, sh) @@ -261,7 +262,7 @@ real function EW_coast(lon, lat, latC, lon0, lon1, dlat, sh) real, intent(in) :: lon1 !< Longitude of coast end [degrees_E] real, intent(in) :: dlat !< "Radius" of coast profile [degrees] real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim] - real :: r + real :: r ! A relative position [degrees] r = dist_line_fixed_y( lon, lat, lon0, lon1, latC ) EW_coast = cstprof(r, 0., dlat, 0.125, 0.125, 0.5, sh) @@ -276,7 +277,7 @@ real function NS_conic_ridge(lon, lat, lonC, lat0, lat1, dlon, rh) real, intent(in) :: lat1 !< Latitude of ridge end [degrees_N] real, intent(in) :: dlon !< "Radius" of ridge profile [degrees] real, intent(in) :: rh !< depth of ridge as fraction of full depth [nondim] - real :: r + real :: r ! A relative position [degrees] r = dist_line_fixed_x( lon, lat, lonC, lat0, lat1 ) NS_conic_ridge = 1. - rh * cone(r, 0., dlon) @@ -291,7 +292,7 @@ real function NS_scurve_ridge(lon, lat, lonC, lat0, lat1, dlon, rh) real, intent(in) :: lat1 !< Latitude of ridge end [degrees_N] real, intent(in) :: dlon !< "Radius" of ridge profile [degrees] real, intent(in) :: rh !< depth of ridge as fraction of full depth [nondim] - real :: r + real :: r ! A relative position [degrees] r = dist_line_fixed_x( lon, lat, lonC, lat0, lat1 ) NS_scurve_ridge = 1. - rh * (1. - scurve(r, 0., dlon) ) @@ -306,12 +307,13 @@ real function circ_conic_ridge(lon, lat, lon0, lat0, ring_radius, ring_thickness real, intent(in) :: ring_radius !< Radius of ring [degrees] real, intent(in) :: ring_thickness !< Radial thickness of ring [degrees] real, intent(in) :: ridge_height !< Ridge height as fraction of full depth [nondim] - real :: r + real :: r ! A relative position [degrees] + real :: frac_ht ! The fractional height of the topography [nondim] r = sqrt( (lon - lon0)**2 + (lat - lat0)**2 ) ! Pseudo-distance from a point r = abs( r - ring_radius) ! Pseudo-distance from a circle - r = cone(r, 0., ring_thickness, ridge_height) ! 0 .. frac_ridge_height - circ_conic_ridge = 1. - r ! nondim depths (1-frac_ridge_height) .. 1 + frac_ht = cone(r, 0., ring_thickness, ridge_height) ! 0 .. frac_ridge_height + circ_conic_ridge = 1. - frac_ht ! nondim depths (1-frac_ridge_height) .. 1 end function circ_conic_ridge !> A circular ridge with cutoff scurve profile @@ -323,13 +325,15 @@ real function circ_scurve_ridge(lon, lat, lon0, lat0, ring_radius, ring_thicknes real, intent(in) :: ring_radius !< Radius of ring [degrees] real, intent(in) :: ring_thickness !< Radial thickness of ring [degrees] real, intent(in) :: ridge_height !< Ridge height as fraction of full depth [nondim] - real :: r + real :: r ! A relative position [degrees] + real :: s ! A function of the normalized position [nondim] + real :: frac_ht ! The fractional height of the topography [nondim] r = sqrt( (lon - lon0)**2 + (lat - lat0)**2 ) ! Pseudo-distance from a point r = abs( r - ring_radius) ! Pseudo-distance from a circle - r = 1. - scurve(r, 0., ring_thickness) ! 0 .. 1 - r = r * ridge_height ! 0 .. frac_ridge_height - circ_scurve_ridge = 1. - r ! nondim depths (1-frac_ridge_height) .. 1 + s = 1. - scurve(r, 0., ring_thickness) ! 0 .. 1 + frac_ht = s * ridge_height ! 0 .. frac_ridge_height + circ_scurve_ridge = 1. - frac_ht ! nondim depths (1-frac_ridge_height) .. 1 end function circ_scurve_ridge end module basin_builder diff --git a/src/user/tidal_bay_initialization.F90 b/src/user/tidal_bay_initialization.F90 index 37a908d3a8..c2ca2565ce 100644 --- a/src/user/tidal_bay_initialization.F90 +++ b/src/user/tidal_bay_initialization.F90 @@ -110,7 +110,7 @@ subroutine tidal_bay_set_OBC_data(OBC, CS, G, GV, US, h, Time) enddo endif enddo ; enddo - total_area = reproducing_sum(my_area) + total_area = US%m_to_Z*US%m_to_L * reproducing_sum(my_area) my_flux = - CS%tide_flow * SIN(2.0*PI*time_sec / CS%tide_period) do n = 1, OBC%number_of_segments @@ -118,7 +118,7 @@ subroutine tidal_bay_set_OBC_data(OBC, CS, G, GV, US, h, Time) if (.not. segment%on_pe) cycle - segment%normal_vel_bt(:,:) = my_flux / (US%m_to_Z*US%m_to_L*total_area) + segment%normal_vel_bt(:,:) = my_flux / total_area segment%SSH(:,:) = cff_eta enddo ! end segment loop