diff --git a/config_src/solo_driver/MOM_surface_forcing.F90 b/config_src/solo_driver/MOM_surface_forcing.F90 index 56d7d5a846..a113d18871 100644 --- a/config_src/solo_driver/MOM_surface_forcing.F90 +++ b/config_src/solo_driver/MOM_surface_forcing.F90 @@ -381,7 +381,6 @@ subroutine wind_forcing_const(sfc_state, forces, tau_x0, tau_y0, day, G, US, CS) Pa_conversion = US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z !set steady surface wind stresses, in units of Pa. - !### mag_tau = US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z * sqrt( tau_x0**2 + tau_y0**2) mag_tau = Pa_conversion * sqrt( tau_x0**2 + tau_y0**2) do j=js,je ; do I=is-1,Ieq @@ -1674,11 +1673,12 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, tracer_flow_C call get_param(param_file, mdl, "LATENT_HEAT_VAPORIZATION", CS%latent_heat_vapor, & "The latent heat of fusion.", units="J/kg", default=hlv) if (CS%restorebuoy) then + ! These three variables use non-standard time units, but are rescaled as they are read. call get_param(param_file, mdl, "FLUXCONST", CS%Flux_const, & "The constant that relates the restoring surface fluxes "//& "to the relative surface anomalies (akin to a piston "//& "velocity). Note the non-MKS units.", & - units="m day-1", scale=US%m_to_Z*US%T_to_s, & + units="m day-1", scale=US%m_to_Z*US%T_to_s/86400.0, & fail_if_missing=.true., unscaled=flux_const_default) if (CS%use_temperature) then @@ -1686,23 +1686,16 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, tracer_flow_C "The constant that relates the restoring surface temperature "//& "flux to the relative surface anomaly (akin to a piston "//& "velocity). Note the non-MKS units.", & - units="m day-1", scale=1.0, & ! scale=US%m_to_Z*US%T_to_s, + units="m day-1", scale=1.0/86400.0, & ! scale=US%m_to_Z*US%T_to_s, default=flux_const_default) call get_param(param_file, mdl, "FLUXCONST_S", CS%Flux_const_S, & "The constant that relates the restoring surface salinity "//& "flux to the relative surface anomaly (akin to a piston "//& "velocity). Note the non-MKS units.", & - units="m day-1", scale=US%m_to_Z*US%T_to_s, & + units="m day-1", scale=US%m_to_Z*US%T_to_s/86400.0, & default=flux_const_default) endif - !### Convert flux constants from m day-1 to m s-1. Folding these into the scaling - ! factors above could change a division into a multiply by a reciprocal, which could - ! change answers at the level of roundoff. - CS%Flux_const = CS%Flux_const / 86400.0 - CS%Flux_const_T = CS%Flux_const_T / 86400.0 - CS%Flux_const_S = CS%Flux_const_S / 86400.0 - if (trim(CS%buoy_config) == "linear") then call get_param(param_file, mdl, "SST_NORTH", CS%T_north, & "With buoy_config linear, the sea surface temperature "//& diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 423cc65687..e23e740c9c 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -574,7 +574,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m if (coordinateMode(coord_mode) == REGRIDDING_ADAPTIVE) then call get_param(param_file, mdl, "ADAPT_TIME_RATIO", adaptTimeRatio, & - "Ratio of ALE timestep to grid timescale.", units="s", default=1e-1) !### Should the units be "nondim"? + "Ratio of ALE timestep to grid timescale.", units="nondim", default=1e-1) call get_param(param_file, mdl, "ADAPT_ZOOM_DEPTH", adaptZoom, & "Depth of near-surface zooming region.", units="m", default=200.0, scale=GV%m_to_H) call get_param(param_file, mdl, "ADAPT_ZOOM_COEFF", adaptZoomCoeff, & diff --git a/src/ALE/coord_slight.F90 b/src/ALE/coord_slight.F90 index 2e41d36473..92de6e1ec3 100644 --- a/src/ALE/coord_slight.F90 +++ b/src/ALE/coord_slight.F90 @@ -687,7 +687,7 @@ subroutine rho_interfaces_col(rho_col, h_col, z_col, rho_tgt, nz, z_col_new, & if (k_layer > 0) then ! The new location is inside of layer k_layer. ! Note that this is coded assuming that this layer is stably stratified. if (.not.(ppoly_i_E(k1,2) > ppoly_i_E(k1,1))) call MOM_error(FATAL, & - "build_grid_SLight: Erroneously searching for an interface in an unstratified layer.") !### COMMENT OUT LATER? + "build_grid_SLight: Erroneously searching for an interface in an unstratified layer.") ! Use the false position method to find the location (degree <= 1) or the first guess. zf = (rt - ppoly_i_E(k1,1)) / (ppoly_i_E(k1,2) - ppoly_i_E(k1,1)) @@ -698,7 +698,7 @@ subroutine rho_interfaces_col(rho_col, h_col, z_col, rho_tgt, nz, z_col_new, & ! Bracket the root. zf1 = 0.0 ; rfn1 = a(1) zf2 = 1.0 ; rfn2 = a(1) + (a(2) + (a(3) + (a(4) + a(5)))) - if (rfn1 * rfn2 > 0.0) call MOM_error(FATAL, "build_grid_SLight: Bad bracketing.") !### COMMENT OUT LATER? + if (rfn1 * rfn2 > 0.0) call MOM_error(FATAL, "build_grid_SLight: Bad bracketing.") do itt=1,max_itt rfn = a(1) + zf*(a(2) + zf*(a(3) + zf*(a(4) + zf*a(5)))) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 690e5250db..18e2c2e5b8 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -194,8 +194,8 @@ module MOM !! multiple coupling timesteps. real :: t_dyn_rel_diag !< The time of the diagnostics relative to diabatic processes and remapping !! [T ~> s]. t_dyn_rel_diag is always positive, since the diagnostics must lag. - integer :: ndyn_per_adv = 0 !< Number of calls to dynamics since the last call to advection. - !### Must be saved if thermo spans coupling? + logical :: preadv_h_stored = .false. !< If true, the thicknesses from before the advective cycle + !! have been stored for use in diagnostics. type(diag_ctrl) :: diag !< structure to regulate diagnostic output timing type(vertvisc_type) :: visc !< structure containing vertical viscosities, @@ -225,6 +225,7 @@ module MOM logical :: debug !< If true, write verbose checksums for debugging purposes. integer :: ntrunc !< number u,v truncations since last call to write_energy + integer :: cont_stencil !< The stencil for thickness from the continuity solver. ! These elements are used to control the dynamics updates. logical :: do_dynamics !< If false, does not call step_MOM_dyn_*. This is an !! undocumented run-time flag that is fragile. @@ -292,6 +293,9 @@ module MOM real :: bad_val_sst_min !< Minimum SST before triggering bad value message [degC] real :: bad_val_sss_max !< Maximum SSS before triggering bad value message [ppt] real :: bad_val_col_thick !< Minimum column thickness before triggering bad value message [m] + logical :: answers_2018 !< If true, use expressions for the surface properties that recover + !! the answers from the end of 2018. Otherwise, use more appropriate + !! expressions that differ at roundoff for non-Boussinsq cases. type(MOM_diag_IDs) :: IDs !< Handles used for diagnostics. type(transport_diag_IDs) :: transport_IDs !< Handles used for transport diagnostics. @@ -330,7 +334,8 @@ module MOM !< Pointer to the MOM along-isopycnal tracer diffusion control structure type(tracer_flow_control_CS), pointer :: tracer_flow_CSp => NULL() !< Pointer to the control structure that orchestrates the calling of tracer packages - !### update_OBC_CS might not be needed outside of initialization? + ! Although update_OBC_CS is not used directly outside of initialization, other modules + ! set pointers to this type, so it should be kept for the duration of the run. type(update_OBC_CS), pointer :: update_OBC_CSp => NULL() !< Pointer to the control structure for updating open boundary condition properties type(ocean_OBC_type), pointer :: OBC => NULL() @@ -660,12 +665,12 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_int_in, CS, & endif ! end of block "(CS%diabatic_first .and. (CS%t_dyn_rel_adv==0.0))" if (do_dyn) then - ! Store pre-dynamics grids for proper diagnostic remapping for transports - ! or advective tendencies. If there are more dynamics steps per advective - ! steps (i.e DT_THERM /= DT), this needs to be stored at the first call. - if (CS%ndyn_per_adv == 0 .and. CS%t_dyn_rel_adv == 0.) then + ! Store pre-dynamics thicknesses for proper diagnostic remapping for transports or + ! advective tendencies. If there are more than one dynamics steps per advective + ! step (i.e DT_THERM > DT), this needs to be stored at the first dynamics call. + if (.not.CS%preadv_h_stored .and. (CS%t_dyn_rel_adv == 0.)) then call diag_copy_diag_to_storage(CS%diag_pre_dyn, h, CS%diag) - CS%ndyn_per_adv = CS%ndyn_per_adv + 1 + CS%preadv_h_stored = .true. endif ! The pre-dynamics velocities might be stored for debugging truncations. @@ -719,7 +724,6 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_int_in, CS, & if (do_advection) then ! Do advective transport and lateral tracer mixing. call step_MOM_tracer_dyn(CS, G, GV, US, h, Time_local) - CS%ndyn_per_adv = 0 if (CS%diabatic_first .and. abs(CS%t_dyn_rel_thermo) > 1e-6*dt) call MOM_error(FATAL, & "step_MOM: Mismatch between the dynamics and diabatic times "//& "with DIABATIC_FIRST.") @@ -927,7 +931,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & call thickness_diffuse(h, CS%uhtr, CS%vhtr, CS%tv, dt_thermo, G, GV, US, & CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp) call cpu_clock_end(id_clock_thick_diff) - call pass_var(h, G%Domain, clock=id_clock_pass) !###, halo=max(2,cont_stensil)) + call pass_var(h, G%Domain, clock=id_clock_pass, halo=max(2,CS%cont_stencil)) call disable_averaging(CS%diag) if (showCallTree) call callTree_waypoint("finished thickness_diffuse_first (step_MOM)") @@ -1002,7 +1006,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & if (CS%debug) call hchksum(h,"Post-thickness_diffuse h", G%HI, haloshift=1, scale=GV%H_to_m) call cpu_clock_end(id_clock_thick_diff) - call pass_var(h, G%Domain, clock=id_clock_pass) !###, halo=max(2,cont_stensil)) + call pass_var(h, G%Domain, clock=id_clock_pass, halo=max(2,CS%cont_stencil)) if (showCallTree) call callTree_waypoint("finished thickness_diffuse (step_MOM)") endif @@ -1017,7 +1021,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & call mixedlayer_restrat(h, CS%uhtr, CS%vhtr, CS%tv, forces, dt, CS%visc%MLD, & CS%VarMix, G, GV, US, CS%mixedlayer_restrat_CSp) call cpu_clock_end(id_clock_ml_restrat) - call pass_var(h, G%Domain, clock=id_clock_pass) !###, halo=max(2,cont_stensil)) + call pass_var(h, G%Domain, clock=id_clock_pass, halo=max(2,CS%cont_stencil)) if (CS%debug) then call hchksum(h,"Post-mixedlayer_restrat h", G%HI, haloshift=1, scale=GV%H_to_m) call uvchksum("Post-mixedlayer_restrat [uv]htr", & @@ -1121,6 +1125,8 @@ subroutine step_MOM_tracer_dyn(CS, G, GV, US, h, Time_local) call do_group_pass(pass_T_S, G%Domain, clock=id_clock_pass) endif + CS%preadv_h_stored = .false. + end subroutine step_MOM_tracer_dyn !> MOM_step_thermo orchestrates the thermodynamic time stepping and vertical @@ -1562,6 +1568,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ! with accumulated heat deficit returned to surface ocean. logical :: bound_salinity ! If true, salt is added to keep salinity above ! a minimum value, and the deficit is reported. + logical :: default_2018_answers ! The default setting for the various 2018_ANSWERS flags. logical :: use_conT_absS ! If true, the prognostics T & S are conservative temperature ! and absolute salinity. Care should be taken to convert them ! to potential temperature and practical salinity before @@ -1875,6 +1882,13 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & "triggered, if CHECK_BAD_SURFACE_VALS is true.", units="m", & default=0.0) endif + call get_param(param_file, "MOM", "DEFAULT_2018_ANSWERS", default_2018_answers, & + "This sets the default value for the various _2018_ANSWERS parameters.", & + default=.true.) + call get_param(param_file, "MOM", "SURFACE_2018_ANSWERS", CS%answers_2018, & + "If true, use expressions for the surface properties that recover the answers "//& + "from the end of 2018. Otherwise, use more appropriate expressions that differ "//& + "at roundoff for non-Boussinsq cases.", default=default_2018_answers) call get_param(param_file, "MOM", "SAVE_INITIAL_CONDS", save_IC, & "If true, write the initial conditions to a file given "//& @@ -2314,7 +2328,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & CS%dt, CS%ADp, CS%CDp, MOM_internal_state, CS%VarMix, CS%MEKE, & CS%thickness_diffuse_CSp, & CS%OBC, CS%update_OBC_CSp, CS%ALE_CSp, CS%set_visc_CSp, & - CS%visc, dirs, CS%ntrunc, calc_dtbt=calc_dtbt) + CS%visc, dirs, CS%ntrunc, calc_dtbt=calc_dtbt, cont_stencil=CS%cont_stencil) if (CS%dtbt_reset_period > 0.0) then CS%dtbt_reset_interval = real_to_time(CS%dtbt_reset_period) ! Set dtbt_reset_time to be the next even multiple of dtbt_reset_interval. @@ -2332,13 +2346,13 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & param_file, diag, CS%dyn_unsplit_RK2_CSp, restart_CSp, & CS%ADp, CS%CDp, MOM_internal_state, CS%MEKE, CS%OBC, & CS%update_OBC_CSp, CS%ALE_CSp, CS%set_visc_CSp, CS%visc, dirs, & - CS%ntrunc) + CS%ntrunc, cont_stencil=CS%cont_stencil) else call initialize_dyn_unsplit(CS%u, CS%v, CS%h, Time, G, GV, US, & param_file, diag, CS%dyn_unsplit_CSp, restart_CSp, & CS%ADp, CS%CDp, MOM_internal_state, CS%MEKE, CS%OBC, & CS%update_OBC_CSp, CS%ALE_CSp, CS%set_visc_CSp, CS%visc, dirs, & - CS%ntrunc) + CS%ntrunc, cont_stencil=CS%cont_stencil) endif call callTree_waypoint("dynamics initialized (initialize_MOM)") @@ -2711,21 +2725,26 @@ subroutine extract_surface_state(CS, sfc_state) !! structure shared with the calling routine !! data in this structure is intent out. - ! local + ! Local variables real :: hu, hv ! Thicknesses interpolated to velocity points [H ~> m or kg m-2] - type(ocean_grid_type), pointer :: G => NULL() !< pointer to a structure containing - !! metrics and related information + type(ocean_grid_type), pointer :: G => NULL() !< pointer to a structure containing + !! metrics and related information type(verticalGrid_type), pointer :: GV => NULL() !< structure containing vertical grid info type(unit_scale_type), pointer :: US => NULL() !< structure containing various unit conversion factors real, dimension(:,:,:), pointer :: & h => NULL() !< h : layer thickness [H ~> m or kg m-2] - real :: depth(SZI_(CS%G)) !< Distance from the surface in depth units [Z ~> m] + real :: depth(SZI_(CS%G)) !< Distance from the surface in depth units [Z ~> m] or [H ~> m or kg m-2] real :: depth_ml !< Depth over which to average to determine mixed - !! layer properties [Z ~> m] - real :: dh !< Thickness of a layer within the mixed layer [Z ~> m] + !! layer properties [Z ~> m] or [H ~> m or kg m-2] + real :: dh !< Thickness of a layer within the mixed layer [Z ~> m] or [H ~> m or kg m-2] real :: mass !< Mass per unit area of a layer [kg m-2] real :: bathy_m !< The depth of bathymetry [m] (not Z), used for error checking. real :: T_freeze !< freezing temperature [degC] + real :: I_depth !< The inverse of depth [Z-1 ~> m-1] or [H-1 ~> m-1 or m2 kg-1] + real :: missing_depth !< The portion of depth_ml that can not be found in a column [H ~> m or kg m-2] + real :: H_rescale !< A conversion factor from thickness units to the units used in the + !! calculation of properties of the uppermost ocean [nondim] or [Z H-1 ~> 1 or m3 kg-1] + ! After the ANSWERS_2018 flag has been obsoleted, H_rescale will be 1. real :: delT(SZI_(CS%G)) !< T-T_freeze [degC] logical :: use_temperature !< If true, temp and saln used as state variables. integer :: i, j, k, is, ie, js, je, nz, numberOfErrors, ig, jg @@ -2777,9 +2796,9 @@ subroutine extract_surface_state(CS, sfc_state) enddo ; enddo else ! (CS%Hmix >= 0.0) - !### This calculation should work in thickness (H) units instead of Z, but that - !### would change answers at roundoff in non-Boussinesq cases. + H_rescale = 1.0 ; if (CS%answers_2018) H_rescale = GV%H_to_Z depth_ml = CS%Hmix + if (.not.CS%answers_2018) depth_ml = CS%Hmix*GV%Z_to_H ! Determine the mean tracer properties of the uppermost depth_ml fluid. !$OMP parallel do default(shared) private(depth,dh) do j=js,je @@ -2793,8 +2812,8 @@ subroutine extract_surface_state(CS, sfc_state) enddo do k=1,nz ; do i=is,ie - if (depth(i) + h(i,j,k)*GV%H_to_Z < depth_ml) then - dh = h(i,j,k)*GV%H_to_Z + if (depth(i) + h(i,j,k)*H_rescale < depth_ml) then + dh = h(i,j,k)*H_rescale elseif (depth(i) < depth_ml) then dh = depth_ml - depth(i) else @@ -2810,16 +2829,36 @@ subroutine extract_surface_state(CS, sfc_state) enddo ; enddo ! Calculate the average properties of the mixed layer depth. do i=is,ie - if (depth(i) < GV%H_subroundoff*GV%H_to_Z) & - depth(i) = GV%H_subroundoff*GV%H_to_Z - if (use_temperature) then - sfc_state%SST(i,j) = sfc_state%SST(i,j) / depth(i) - sfc_state%SSS(i,j) = sfc_state%SSS(i,j) / depth(i) + if (CS%answers_2018) then + if (depth(i) < GV%H_subroundoff*H_rescale) & + depth(i) = GV%H_subroundoff*H_rescale + if (use_temperature) then + sfc_state%SST(i,j) = sfc_state%SST(i,j) / depth(i) + sfc_state%SSS(i,j) = sfc_state%SSS(i,j) / depth(i) + else + sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) / depth(i) + endif else - sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) / depth(i) + if (depth(i) < GV%H_subroundoff*H_rescale) then + I_depth = 1.0 / (GV%H_subroundoff*H_rescale) + missing_depth = GV%H_subroundoff*H_rescale - depth(i) + if (use_temperature) then + sfc_state%SST(i,j) = (sfc_state%SST(i,j) + missing_depth*CS%tv%T(i,j,1)) * I_depth + sfc_state%SSS(i,j) = (sfc_state%SSS(i,j) + missing_depth*CS%tv%S(i,j,1)) * I_depth + else + sfc_state%sfc_density(i,j) = (sfc_state%sfc_density(i,j) + & + missing_depth*US%R_to_kg_m3*GV%Rlay(1)) * I_depth + endif + else + I_depth = 1.0 / depth(i) + if (use_temperature) then + sfc_state%SST(i,j) = sfc_state%SST(i,j) * I_depth + sfc_state%SSS(i,j) = sfc_state%SSS(i,j) * I_depth + else + sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) * I_depth + endif + endif endif - !### Verify that this is no longer needed. - ! sfc_state%Hml(i,j) = US%Z_to_m * depth(i) enddo enddo ! end of j loop @@ -2828,9 +2867,8 @@ subroutine extract_surface_state(CS, sfc_state) ! required by the speed diagnostic on the non-symmetric grid. ! This assumes that u and v halos have already been updated. if (CS%Hmix_UV>0.) then - !### This calculation should work in thickness (H) units instead of Z, but that - !### would change answers at roundoff in non-Boussinesq cases. depth_ml = CS%Hmix_UV + if (.not.CS%answers_2018) depth_ml = CS%Hmix_UV*GV%Z_to_H !$OMP parallel do default(shared) private(depth,dh,hv) do J=js-1,ie do i=is,ie @@ -2838,7 +2876,7 @@ subroutine extract_surface_state(CS, sfc_state) sfc_state%v(i,J) = 0.0 enddo do k=1,nz ; do i=is,ie - hv = 0.5 * (h(i,j,k) + h(i,j+1,k)) * GV%H_to_Z + hv = 0.5 * (h(i,j,k) + h(i,j+1,k)) * H_rescale if (depth(i) + hv < depth_ml) then dh = hv elseif (depth(i) < depth_ml) then @@ -2851,9 +2889,7 @@ subroutine extract_surface_state(CS, sfc_state) enddo ; enddo ! Calculate the average properties of the mixed layer depth. do i=is,ie - if (depth(i) < GV%H_subroundoff*GV%H_to_Z) & - depth(i) = GV%H_subroundoff*GV%H_to_Z - sfc_state%v(i,J) = sfc_state%v(i,J) / depth(i) + sfc_state%v(i,J) = sfc_state%v(i,J) / max(depth(i), GV%H_subroundoff*H_rescale) enddo enddo ! end of j loop @@ -2864,7 +2900,7 @@ subroutine extract_surface_state(CS, sfc_state) sfc_state%u(I,j) = 0.0 enddo do k=1,nz ; do I=is-1,ie - hu = 0.5 * (h(i,j,k) + h(i+1,j,k)) * GV%H_to_Z + hu = 0.5 * (h(i,j,k) + h(i+1,j,k)) * H_rescale if (depth(i) + hu < depth_ml) then dh = hu elseif (depth(I) < depth_ml) then @@ -2877,9 +2913,7 @@ subroutine extract_surface_state(CS, sfc_state) enddo ; enddo ! Calculate the average properties of the mixed layer depth. do I=is-1,ie - if (depth(I) < GV%H_subroundoff*GV%H_to_Z) & - depth(I) = GV%H_subroundoff*GV%H_to_Z - sfc_state%u(I,j) = sfc_state%u(I,j) / depth(I) + sfc_state%u(I,j) = sfc_state%u(I,j) / max(depth(I), GV%H_subroundoff*H_rescale) enddo enddo ! end of j loop else ! Hmix_UV<=0. @@ -2933,7 +2967,7 @@ subroutine extract_surface_state(CS, sfc_state) !$OMP parallel do default(shared) do j=js,je ; do i=is,ie ! Convert from gSalt to kgSalt - sfc_state%salt_deficit(i,j) = 1000.0 * US%R_to_kg_m3*US%Z_to_m*CS%tv%salt_deficit(i,j) + sfc_state%salt_deficit(i,j) = 0.001 * US%R_to_kg_m3*US%Z_to_m*CS%tv%salt_deficit(i,j) enddo ; enddo endif if (allocated(sfc_state%TempxPmE) .and. associated(CS%tv%TempxPmE)) then diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 14fc918b60..0ccf4d8f3b 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -203,8 +203,16 @@ module MOM_barotropic !! (false) is to use a predictor continuity step to !! find the pressure field, and then do a corrector !! continuity step using a weighted average of the - !! old and new velocities, with weights of (1-BEBT) - !! and BEBT. + !! old and new velocities, with weights of (1-BEBT) and BEBT. + logical :: nonlin_stress !< If true, use the full depth of the ocean at the start of the + !! barotropic step when calculating the surface stress contribution to + !! the barotropic acclerations. Otherwise use the depth based on bathyT. + real :: BT_Coriolis_scale !< A factor by which the barotropic Coriolis acceleration anomaly + !! terms are scaled. + logical :: answers_2018 !< If true, use expressions for the barotropic solver that recover + !! the answers from the end of 2018. Otherwise, use more efficient + !! or general expressions. + logical :: dynamic_psurf !< If true, add a dynamic pressure due to a viscous !! ice shelf, for instance. real :: Dmin_dyn_psurf !< The minimum depth to use in limiting the size @@ -327,6 +335,7 @@ module MOM_barotropic real :: uh_WW !< The zonal transport when ubt=ubt_WW [H L2 T-1 ~> m3 s-1 or kg s-1]. real :: uh_EE !< The zonal transport when ubt=ubt_EE [H L2 T-1 ~> m3 s-1 or kg s-1]. end type local_BT_cont_u_type + !> A desciption of the functional dependence of transport at a v-point type, private :: local_BT_cont_v_type real :: FA_v_NN !< The effective open face area for meridional barotropic transport @@ -476,7 +485,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! relative to eta_PF, with SAL effects included [H ~> m or kg m-2]. ! These are always allocated with symmetric memory and wide halos. - real :: q(SZIBW_(CS),SZJBW_(CS)) ! A pseudo potential vorticity [T-1 Z-1 ~> s-1 m-1]. + real :: q(SZIBW_(CS),SZJBW_(CS)) ! A pseudo potential vorticity [T-1 Z-1 ~> s-1 m-1] + ! or [T-1 H-1 ~> s-1 m-1 or m2 s-1 kg-1] real, dimension(SZIBW_(CS),SZJW_(CS)) :: & ubt, & ! The zonal barotropic velocity [L T-1 ~> m s-1]. bt_rem_u, & ! The fraction of the barotropic zonal velocity that remains @@ -507,7 +517,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, Rayleigh_u, & ! A Rayleigh drag timescale operating at u-points [T-1 ~> s-1]. PFu_bt_sum, & ! The summed zonal barotropic pressure gradient force [L T-2 ~> m s-2]. Coru_bt_sum, & ! The summed zonal barotropic Coriolis acceleration [L T-2 ~> m s-2]. - DCor_u, & ! A simply averaged depth at u points [Z ~> m]. + DCor_u, & ! An averaged depth or total thickness at u points [Z ~> m] or [H ~> m or kg m-2]. Datu ! Basin depth at u-velocity grid points times the y-grid ! spacing [H L ~> m2 or kg m-1]. real, dimension(SZIW_(CS),SZJBW_(CS)) :: & @@ -538,7 +548,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! [L T-2 ~> m s-2]. Corv_bt_sum, & ! The summed meridional barotropic Coriolis acceleration, ! [L T-2 ~> m s-2]. - DCor_v, & ! A simply averaged depth at v points [Z ~> m]. + DCor_v, & ! An averaged depth or total thickness at v points [Z ~> m] or [H ~> m or kg m-2]. Datv ! Basin depth at v-velocity grid points times the x-grid ! spacing [H L ~> m2 or kg m-1]. real, target, dimension(SZIW_(CS),SZJW_(CS)) :: & @@ -606,6 +616,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, time_step_end, & ! The end time of a barotropic step. time_end_in ! The end time for diagnostics when this routine started. real :: time_int_in ! The diagnostics' time interval when this routine started. + real :: Htot_avg ! The average total thickness of the tracer columns adjacent to a + ! velocity point [H ~> m or kg m-2] logical :: do_hifreq_output ! If true, output occurs every barotropic step. logical :: use_BT_cont, do_ave, find_etaav, find_PF, find_Cor logical :: ice_is_rigid, nonblock_setup, interp_eta_PF @@ -624,6 +636,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, real :: Htot ! The total thickness [H ~> m or kg m-2]. real :: eta_cor_max ! The maximum fluid that can be added as a correction to eta [H ~> m or kg m-2]. real :: accel_underflow ! An acceleration that is so small it should be zeroed out [L T-2 ~> m s-2]. + real :: h_neglect ! A thickness that is so small it is usually lost + ! in roundoff and can be neglected [H ~> m or kg m-2]. real, allocatable, dimension(:) :: wt_vel, wt_eta, wt_accel, wt_trans, wt_accel2 real :: sum_wt_vel, sum_wt_eta, sum_wt_accel, sum_wt_trans @@ -651,6 +665,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB MS%isdw = CS%isdw ; MS%iedw = CS%iedw ; MS%jsdw = CS%jsdw ; MS%jedw = CS%jedw + h_neglect = GV%H_subroundoff Idt = 1.0 / dt accel_underflow = CS%vel_underflow * Idt @@ -815,23 +830,43 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, enddo ; enddo else q(:,:) = 0.0 ; DCor_u(:,:) = 0.0 ; DCor_v(:,:) = 0.0 - ! This option has not yet been written properly. - ! ### bathyT here should be replaced with bathyT+eta(Bous) or eta(non-Bous). - !$OMP parallel do default(shared) - do j=js,je ; do I=is-1,ie - DCor_u(I,j) = 0.5 * (G%bathyT(i+1,j) + G%bathyT(i,j)) - enddo ; enddo - !$OMP parallel do default(shared) - do J=js-1,je ; do i=is,ie - DCor_v(i,J) = 0.5 * (G%bathyT(i,j+1) + G%bathyT(i,j)) - enddo ; enddo - !$OMP parallel do default(shared) - do J=js-1,je ; do I=is-1,ie - q(I,J) = 0.25 * G%CoriolisBu(I,J) * & - ((G%areaT(i,j) + G%areaT(i+1,j+1)) + (G%areaT(i+1,j) + G%areaT(i,j+1))) / & - ((G%areaT(i,j) * G%bathyT(i,j) + G%areaT(i+1,j+1) * G%bathyT(i+1,j+1)) + & - (G%areaT(i+1,j) * G%bathyT(i+1,j) + G%areaT(i,j+1) * G%bathyT(i,j+1)) ) - enddo ; enddo + if (GV%Boussinesq) then + !$OMP parallel do default(shared) + do j=js,je ; do I=is-1,ie + DCor_u(I,j) = 0.5 * (max(GV%Z_to_H*G%bathyT(i+1,j) + eta_in(i+1,j), 0.0) + & + max(GV%Z_to_H*G%bathyT(i,j) + eta_in(i,j), 0.0) ) + enddo ; enddo + !$OMP parallel do default(shared) + do J=js-1,je ; do i=is,ie + DCor_v(i,J) = 0.5 * (max(GV%Z_to_H*G%bathyT(i,j+1) + eta_in(i+1,j), 0.0) + & + max(GV%Z_to_H*G%bathyT(i,j) + eta_in(i,j), 0.0) ) + enddo ; enddo + !$OMP parallel do default(shared) + do J=js-1,je ; do I=is-1,ie + q(I,J) = 0.25 * (CS%BT_Coriolis_scale * G%CoriolisBu(I,J)) * & + ((G%areaT(i,j) + G%areaT(i+1,j+1)) + (G%areaT(i+1,j) + G%areaT(i,j+1))) / & + (max((G%areaT(i,j) * max(GV%Z_to_H*G%bathyT(i,j) + eta_in(i,j), 0.0) + & + G%areaT(i+1,j+1) * max(GV%Z_to_H*G%bathyT(i+1,j+1) + eta_in(i+1,j+1), 0.0)) + & + (G%areaT(i+1,j) * max(GV%Z_to_H*G%bathyT(i+1,j) + eta_in(i+1,j), 0.0) + & + G%areaT(i,j+1) * max(GV%Z_to_H*G%bathyT(i,j+1) + eta_in(i,j+1), 0.0)), h_neglect) ) + enddo ; enddo + else + !$OMP parallel do default(shared) + do j=js,je ; do I=is-1,ie + DCor_u(I,j) = 0.5 * (eta_in(i+1,j) + eta_in(i,j)) + enddo ; enddo + !$OMP parallel do default(shared) + do J=js-1,je ; do i=is,ie + DCor_v(i,J) = 0.5 * (eta_in(i,j+1) + eta_in(i,j)) + enddo ; enddo + !$OMP parallel do default(shared) + do J=js-1,je ; do I=is-1,ie + q(I,J) = 0.25 * (CS%BT_Coriolis_scale * G%CoriolisBu(I,J)) * & + ((G%areaT(i,j) + G%areaT(i+1,j+1)) + (G%areaT(i+1,j) + G%areaT(i,j+1))) / & + (max((G%areaT(i,j) * eta_in(i,j) + G%areaT(i+1,j+1) * eta_in(i+1,j+1)) + & + (G%areaT(i+1,j) * eta_in(i+1,j) + G%areaT(i,j+1) * eta_in(i,j+1)), h_neglect) ) + enddo ; enddo + endif ! With very wide halos, q and D need to be calculated on the available data ! domain and then updated onto the full computational domain. @@ -976,50 +1011,6 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, Datu, Datv, BTCL_u, BTCL_v) endif -! Here the vertical average accelerations due to the Coriolis, advective, -! pressure gradient and horizontal viscous terms in the layer momentum -! equations are calculated. These will be used to determine the difference -! between the accelerations due to the average of the layer equations and the -! barotropic calculation. - - !$OMP parallel do default(shared) - do j=js,je ; do I=is-1,ie - ! ### IDatu here should be replaced with 1/D+eta(Bous) or 1/eta(non-Bous). - ! ### although with BT_cont_types IDatu should be replaced by - ! ### CS%dy_Cu(I,j) / (d(uhbt)/du) (with appropriate bounds). - BT_force_u(I,j) = forces%taux(I,j) * mass_accel_to_Z * CS%IDatu(I,j)*visc_rem_u(I,j,1) - enddo ; enddo - !$OMP parallel do default(shared) - do J=js-1,je ; do i=is,ie - ! ### IDatv here should be replaced with 1/D+eta(Bous) or 1/eta(non-Bous). - ! ### although with BT_cont_types IDatv should be replaced by - ! ### CS%dx_Cv(I,j) / (d(vhbt)/dv) (with appropriate bounds). - BT_force_v(i,J) = forces%tauy(i,J) * mass_accel_to_Z * CS%IDatv(i,J)*visc_rem_v(i,J,1) - enddo ; enddo - if (present(taux_bot) .and. present(tauy_bot)) then - if (associated(taux_bot) .and. associated(tauy_bot)) then - !$OMP parallel do default(shared) - do j=js,je ; do I=is-1,ie - BT_force_u(I,j) = BT_force_u(I,j) - taux_bot(I,j) * mass_to_Z * CS%IDatu(I,j) - enddo ; enddo - !$OMP parallel do default(shared) - do J=js-1,je ; do i=is,ie - BT_force_v(i,J) = BT_force_v(i,J) - tauy_bot(i,J) * mass_to_Z * CS%IDatv(i,J) - enddo ; enddo - endif - endif - - ! bc_accel_u & bc_accel_v are only available on the potentially - ! non-symmetric computational domain. - !$OMP parallel do default(shared) - do j=js,je ; do k=1,nz ; do I=Isq,Ieq - BT_force_u(I,j) = BT_force_u(I,j) + wt_u(I,j,k) * bc_accel_u(I,j,k) - enddo ; enddo ; enddo - !$OMP parallel do default(shared) - do J=Jsq,Jeq ; do k=1,nz ; do i=is,ie - BT_force_v(i,J) = BT_force_v(i,J) + wt_v(i,J,k) * bc_accel_v(i,J,k) - enddo ; enddo ; enddo - ! Determine the difference between the sum of the layer fluxes and the ! barotropic fluxes found from the same input velocities. if (add_uh0) then @@ -1128,6 +1119,82 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ubt_first(:,:) = ubt(:,:) ; vbt_first(:,:) = vbt(:,:) endif +! Here the vertical average accelerations due to the Coriolis, advective, +! pressure gradient and horizontal viscous terms in the layer momentum +! equations are calculated. These will be used to determine the difference +! between the accelerations due to the average of the layer equations and the +! barotropic calculation. + + !$OMP parallel do default(shared) + do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) > 0.0) then + if (CS%nonlin_stress) then + if (GV%Boussinesq) then + Htot_avg = 0.5*(max(CS%bathyT(i,j)*GV%Z_to_H + eta(i,j), 0.0) + & + max(CS%bathyT(i+1,j)*GV%Z_to_H + eta(i+1,j), 0.0)) + else + Htot_avg = 0.5*(eta(i,j) + eta(i+1,j)) + endif + if (Htot_avg*CS%dy_Cu(I,j) <= 0.0) then + CS%IDatu(I,j) = 0.0 + elseif (use_BT_cont) then ! Reconsider the max and whether there should be some scaling. + CS%IDatu(I,j) = CS%dy_Cu(I,j) / (max(find_duhbt_du(ubt(I,j), BTCL_u(I,j), US), & + CS%dy_Cu(I,j)*Htot_avg) ) + else + CS%IDatu(I,j) = 1.0 / Htot_avg + endif + endif + + BT_force_u(I,j) = forces%taux(I,j) * mass_accel_to_Z * CS%IDatu(I,j)*visc_rem_u(I,j,1) + else + BT_force_u(I,j) = 0.0 + endif ; enddo ; enddo + !$OMP parallel do default(shared) + do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.0) then + if (CS%nonlin_stress) then + if (GV%Boussinesq) then + Htot_avg = 0.5*(max(CS%bathyT(i,j)*GV%Z_to_H + eta(i,j), 0.0) + & + max(CS%bathyT(i,j+1)*GV%Z_to_H + eta(i,j+1), 0.0)) + else + Htot_avg = 0.5*(eta(i,j) + eta(i,j+1)) + endif + if (Htot_avg*CS%dx_Cv(i,J) <= 0.0) then + CS%IDatv(i,J) = 0.0 + elseif (use_BT_cont) then ! Reconsider the max and whether there should be some scaling. + CS%IDatv(i,J) = CS%dx_Cv(i,J) / (max(find_dvhbt_dv(vbt(i,J), BTCL_v(i,J), US), & + CS%dx_Cv(i,J)*Htot_avg) ) + else + CS%IDatv(i,J) = 1.0 / Htot_avg + endif + endif + + BT_force_v(i,J) = forces%tauy(i,J) * mass_accel_to_Z * CS%IDatv(i,J)*visc_rem_v(i,J,1) + else + BT_force_v(i,J) = 0.0 + endif ; enddo ; enddo + if (present(taux_bot) .and. present(tauy_bot)) then + if (associated(taux_bot) .and. associated(tauy_bot)) then + !$OMP parallel do default(shared) + do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) > 0.0) then + BT_force_u(I,j) = BT_force_u(I,j) - taux_bot(I,j) * mass_to_Z * CS%IDatu(I,j) + endif ; enddo ; enddo + !$OMP parallel do default(shared) + do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.0) then + BT_force_v(i,J) = BT_force_v(i,J) - tauy_bot(i,J) * mass_to_Z * CS%IDatv(i,J) + endif ; enddo ; enddo + endif + endif + + ! bc_accel_u & bc_accel_v are only available on the potentially + ! non-symmetric computational domain. + !$OMP parallel do default(shared) + do j=js,je ; do k=1,nz ; do I=Isq,Ieq + BT_force_u(I,j) = BT_force_u(I,j) + wt_u(I,j,k) * bc_accel_u(I,j,k) + enddo ; enddo ; enddo + !$OMP parallel do default(shared) + do J=Jsq,Jeq ; do k=1,nz ; do i=is,ie + BT_force_v(i,J) = BT_force_v(i,J) + wt_v(i,J,k) * bc_accel_v(i,J,k) + enddo ; enddo ; enddo + if (CS%gradual_BT_ICs) then !$OMP parallel do default(shared) do j=js,je ; do I=is-1,ie @@ -1403,7 +1470,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, (gtot_N(i,j) * (Datv(i,J)*G%IdyCv(i,J)) + & gtot_S(i,j) * (Datv(i,J-1)*G%IdyCv(i,J-1)))) + & ((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + & - (G%CoriolisBu(I-1,J)**2 + G%CoriolisBu(I,J-1)**2))) + (G%CoriolisBu(I-1,J)**2 + G%CoriolisBu(I,J-1)**2)) * CS%BT_Coriolis_scale**2 ) H_eff_dx2 = max(H_min_dyn * ((G%IdxT(i,j))**2 + (G%IdyT(i,j))**2), & G%IareaT(i,j) * & ((Datu(I,j)*G%IdxCu(I,j) + Datu(I-1,j)*G%IdxCu(I-1,j)) + & @@ -1441,8 +1508,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, if (id_clock_calc_pre > 0) call cpu_clock_end(id_clock_calc_pre) if (id_clock_pass_pre > 0) call cpu_clock_begin(id_clock_pass_pre) - if (.not.use_BT_cont) & !### IS THIS OK HERE? - call complete_group_pass(CS%pass_Dat_uv, CS%BT_Domain) + if (.not.use_BT_cont) call complete_group_pass(CS%pass_Dat_uv, CS%BT_Domain) call complete_group_pass(CS%pass_force_hbt0_Cor_ref, CS%BT_Domain) call complete_group_pass(CS%pass_eta_bt_rem, CS%BT_Domain) @@ -1505,6 +1571,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, allocate(wt_accel2(nstep+nfilter+1)) do n=1,nstep+nfilter ! Modify this to use a different filter... + + ! This is a filter that ramps down linearly over a time dt_filt. if ( (n==nstep) .or. (dt_filt - abs(n-nstep)*dtbt >= 0.0)) then wt_vel(n) = 1.0 ; wt_eta(n) = 1.0 elseif (dtbt + dt_filt - abs(n-nstep)*dtbt > 0.0) then @@ -1512,8 +1580,9 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, else wt_vel(n) = 0.0 ; wt_eta(n) = 0.0 endif -!### if (n < nstep-nfilter) then ; wt_vel(n) = 0.0 ; else ; wt_vel(n) = 1.0 ; endif -!### if (n < nstep-nfilter) then ; wt_eta(n) = 0.0 ; else ; wt_eta(n) = 1.0 ; endif + ! This is a simple stepfunction filter. + ! if (n < nstep-nfilter) then ; wt_vel(n) = 0.0 ; else ; wt_vel(n) = 1.0 ; endif + ! wt_eta(n) = wt_vel(n) ! The rest should not be changed. sum_wt_vel = sum_wt_vel + wt_vel(n) ; sum_wt_eta = sum_wt_eta + wt_eta(n) @@ -1529,13 +1598,17 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, I_sum_wt_eta = 1.0 / sum_wt_eta ; I_sum_wt_trans = 1.0 / sum_wt_trans do n=1,nstep+nfilter wt_vel(n) = wt_vel(n) * I_sum_wt_vel - wt_accel2(n) = wt_accel(n) + if (CS%answers_2018) then + wt_accel2(n) = wt_accel(n) + ! wt_trans(n) = wt_trans(n) * I_sum_wt_trans + else + wt_accel2(n) = wt_accel(n) * I_sum_wt_accel + wt_trans(n) = wt_trans(n) * I_sum_wt_trans + endif wt_accel(n) = wt_accel(n) * I_sum_wt_accel wt_eta(n) = wt_eta(n) * I_sum_wt_eta -! wt_trans(n) = wt_trans(n) * I_sum_wt_trans enddo - sum_wt_vel = 0.0 ; sum_wt_eta = 0.0 ; sum_wt_accel = 0.0 ; sum_wt_trans = 0.0 ! The following loop contains all of the time steps. @@ -2003,7 +2076,6 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, eta(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * CS%IareaT(i,j)) * & ((uhbt(I-1,j) - uhbt(I,j)) + (vhbt(i,J-1) - vhbt(i,J))) eta_wtd(i,j) = eta_wtd(i,j) + eta(i,j) * wt_eta(n) - ! Should there be a concern if eta drops below 0 or G%bathyT? enddo ; enddo if (do_hifreq_output) then @@ -2024,6 +2096,18 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, call hchksum(eta, trim(mesg)//" eta", CS%debug_BT_HI, haloshift=iev-ie, scale=GV%H_to_m) endif + if (GV%Boussinesq) then + do j=js,je ; do i=is,ie + if (eta(i,j) < -GV%Z_to_H*G%bathyT(i,j)) & + call MOM_error(WARNING, "btstep: eta has dropped below bathyT.") + enddo ; enddo + else + do j=js,je ; do i=is,ie + if (eta(i,j) < 0.0) & + call MOM_error(WARNING, "btstep: negative eta in a non-Boussinesq barotropic solver.") + enddo ; enddo + endif + enddo ! end of do n=1,ntimestep if (id_clock_calc > 0) call cpu_clock_end(id_clock_calc) if (id_clock_calc_post > 0) call cpu_clock_begin(id_clock_calc_post) @@ -2031,8 +2115,12 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! Reset the time information in the diag type. if (do_hifreq_output) call enable_averaging(time_int_in, time_end_in, CS%diag) - I_sum_wt_vel = 1.0 / sum_wt_vel ; I_sum_wt_eta = 1.0 / sum_wt_eta - I_sum_wt_accel = 1.0 / sum_wt_accel ; I_sum_wt_trans = 1.0 / sum_wt_trans + if (CS%answers_2018) then + I_sum_wt_vel = 1.0 / sum_wt_vel ; I_sum_wt_eta = 1.0 / sum_wt_eta + I_sum_wt_accel = 1.0 / sum_wt_accel ; I_sum_wt_trans = 1.0 / sum_wt_trans + else + I_sum_wt_vel = 1.0 ; I_sum_wt_eta = 1.0 ; I_sum_wt_accel = 1.0 ; I_sum_wt_trans = 1.0 + endif if (find_etaav) then ; do j=js,je ; do i=is,ie etaav(i,j) = eta_sum(i,j) * I_sum_wt_accel @@ -2089,21 +2177,30 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, if (id_clock_pass_post > 0) call cpu_clock_end(id_clock_pass_post) if (id_clock_calc_post > 0) call cpu_clock_begin(id_clock_calc_post) - do j=js,je ; do I=is-1,ie - CS%ubtav(I,j) = ubt_sum(I,j) * I_sum_wt_trans - uhbtav(I,j) = uhbt_sum(I,j) * I_sum_wt_trans - ! The following line would do approximately nothing, as I_sum_wt_accel ~= 1. - !### u_accel_bt(I,j) = u_accel_bt(I,j) * I_sum_wt_accel - ubt_wtd(I,j) = ubt_wtd(I,j) * I_sum_wt_vel - enddo ; enddo + if (CS%answers_2018) then + do j=js,je ; do I=is-1,ie + CS%ubtav(I,j) = ubt_sum(I,j) * I_sum_wt_trans + uhbtav(I,j) = uhbt_sum(I,j) * I_sum_wt_trans + ubt_wtd(I,j) = ubt_wtd(I,j) * I_sum_wt_vel + enddo ; enddo + + do J=js-1,je ; do i=is,ie + CS%vbtav(i,J) = vbt_sum(i,J) * I_sum_wt_trans + vhbtav(i,J) = vhbt_sum(i,J) * I_sum_wt_trans + vbt_wtd(i,J) = vbt_wtd(i,J) * I_sum_wt_vel + enddo ; enddo + else + do j=js,je ; do I=is-1,ie + CS%ubtav(I,j) = ubt_sum(I,j) + uhbtav(I,j) = uhbt_sum(I,j) + enddo ; enddo + + do J=js-1,je ; do i=is,ie + CS%vbtav(i,J) = vbt_sum(i,J) + vhbtav(i,J) = vhbt_sum(i,J) + enddo ; enddo + endif - do J=js-1,je ; do i=is,ie - CS%vbtav(i,J) = vbt_sum(i,J) * I_sum_wt_trans - vhbtav(i,J) = vhbt_sum(i,J) * I_sum_wt_trans - ! The following line would do approximately nothing, as I_sum_wt_accel ~= 1. - !### v_accel_bt(i,J) = v_accel_bt(i,J) * I_sum_wt_accel - vbt_wtd(i,J) = vbt_wtd(i,J) * I_sum_wt_vel - enddo ; enddo if (id_clock_calc_post > 0) call cpu_clock_end(id_clock_calc_post) if (id_clock_pass_post > 0) call cpu_clock_begin(id_clock_pass_post) @@ -2359,7 +2456,7 @@ subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add) ((gtot_E(i,j)*Datu(I,j)*G%IdxCu(I,j) + gtot_W(i,j)*Datu(I-1,j)*G%IdxCu(I-1,j)) + & (gtot_N(i,j)*Datv(i,J)*G%IdyCv(i,J) + gtot_S(i,j)*Datv(i,J-1)*G%IdyCv(i,J-1))) + & ((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + & - (G%CoriolisBu(I-1,J)**2 + G%CoriolisBu(I,J-1)**2))) + (G%CoriolisBu(I-1,J)**2 + G%CoriolisBu(I,J-1)**2)) * CS%BT_Coriolis_scale**2 ) if (Idt_max2 * min_max_dt2 > 1.0) min_max_dt2 = 1.0 / Idt_max2 enddo ; enddo dtbt_max = sqrt(min_max_dt2 / dgeo_de) @@ -3050,6 +3147,31 @@ function find_uhbt(u, BTC, US) result(uhbt) end function find_uhbt +!> The function find_duhbt_du determines the marginal zonal face area for a given velocity. +function find_duhbt_du(u, BTC, US) result(duhbt_du) + real, intent(in) :: u !< The local zonal velocity [L T-1 ~> m s-1] + type(local_BT_cont_u_type), intent(in) :: BTC !< A structure containing various fields that + !! allow the barotropic transports to be calculated consistently + !! with the layers' continuity equations. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + + real :: duhbt_du !< The zonal barotropic face area [L H ~> m2] + + if (u == 0.0) then + duhbt_du = 0.5*(BTC%FA_u_E0 + BTC%FA_u_W0) ! Note the potential discontinuity here. + elseif (u < BTC%uBT_EE) then + duhbt_du = BTC%FA_u_EE + elseif (u < 0.0) then + duhbt_du = (BTC%FA_u_E0 + 3.0*BTC%uh_crvE * u**2) + elseif (u <= BTC%uBT_WW) then + duhbt_du = (BTC%FA_u_W0 + 3.0*BTC%uh_crvW * u**2) + else ! (u > BTC%uBT_WW) + duhbt_du = BTC%FA_u_WW + endif + +end function find_duhbt_du + + !> This function inverts the transport function to determine the barotopic !! velocity that is consistent with a given transport. function uhbt_to_ubt(uhbt, BTC, US, guess) result(ubt) @@ -3167,6 +3289,29 @@ function find_vhbt(v, BTC, US) result(vhbt) end function find_vhbt +!> The function find_vhbt determines the meridional transport for a given velocity. +function find_dvhbt_dv(v, BTC, US) result(dvhbt_dv) + real, intent(in) :: v !< The local meridional velocity [L T-1 ~> m s-1] + type(local_BT_cont_v_type), intent(in) :: BTC !< A structure containing various fields that + !! allow the barotropic transports to be calculated consistently + !! with the layers' continuity equations. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real :: dvhbt_dv !< The meridional barotropic face area [L H ~> m2] + + if (v == 0.0) then + dvhbt_dv = 0.5*(BTC%FA_v_N0 + BTC%FA_v_S0) ! Note the potential discontinuity here. + elseif (v < BTC%vBT_NN) then + dvhbt_dv = BTC%FA_v_NN + elseif (v < 0.0) then + dvhbt_dv = BTC%FA_v_N0 + 3.0*BTC%vh_crvN * v**2 + elseif (v <= BTC%vBT_SS) then + dvhbt_dv = BTC%FA_v_S0 + 3.0*BTC%vh_crvS * v**2 + else ! (v > BTC%vBT_SS) + dvhbt_dv = BTC%FA_v_SS + endif + +end function find_dvhbt_dv + !> This function inverts the transport function to determine the barotopic !! velocity that is consistent with a given transport. function vhbt_to_vbt(vhbt, BTC, US, guess) result(vbt) @@ -3752,10 +3897,13 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, ! a restart file to the internal representation in this run. real :: uH_rescale ! A rescaling factor for thickness transports from the representation in ! a restart file to the internal representation in this run. + real :: mean_SL ! The mean sea level that is used along with the bathymetry to estimate the + ! geometry when LINEARIZED_BT_CORIOLIS is true or BT_NONLIN_STRESS is false [Z ~> m]. real, allocatable, dimension(:,:) :: lin_drag_h type(memory_size_type) :: MS type(group_pass_type) :: pass_static_data, pass_q_D_Cor type(group_pass_type) :: pass_bt_hbt_btav, pass_a_polarity + logical :: default_2018_answers ! The default setting for the various 2018_ANSWERS flags. logical :: apply_bt_drag, use_BT_cont_type character(len=48) :: thickness_units, flux_units character*(40) :: hvel_str @@ -3862,6 +4010,10 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, "to do a corrector continuity step using a weighted "//& "average of the old and new velocities, with weights "//& "of (1-BEBT) and BEBT.", default=.false.) + call get_param(param_file, mdl, "BT_NONLIN_STRESS", CS%nonlin_stress, & + "If true, use the full depth of the ocean at the start of the barotropic "//& + "step when calculating the surface stress contribution to the barotropic "//& + "acclerations. Otherwise use the depth based on bathyT.", default=.false.) call get_param(param_file, mdl, "DYNAMIC_SURFACE_PRESSURE", CS%dynamic_psurf, & "If true, add a dynamic pressure due to a viscous ice "//& @@ -3882,6 +4034,16 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, "if DYNAMIC_SURFACE_PRESSURE is true. Stable values "//& "are < ~1.0.", units="nondim", default=0.9) endif + call get_param(param_file, mdl, "BT_CORIOLIS_SCALE", CS%BT_Coriolis_scale, & + "A factor by which the barotropic Coriolis anomaly terms are scaled.", & + units="nondim", default=1.0) + call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, & + "This sets the default value for the various _2018_ANSWERS parameters.", & + default=.true.) + call get_param(param_file, mdl, "BAROTROPIC_2018_ANSWERS", CS%answers_2018, & + "If true, use expressions for the barotropic solver that recover the answers "//& + "from the end of 2018. Otherwise, use more efficient or general expressions.", & + default=default_2018_answers) call get_param(param_file, mdl, "TIDES", CS%tides, & "If true, apply tidal momentum forcing.", default=.false.) @@ -3990,7 +4152,9 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, "quite large if this is true.", default=CS%debug, & debuggingParam=.true.) - CS%linearized_BT_PV = .true. + call get_param(param_file, mdl, "LINEARIZED_BT_CORIOLIS", CS%linearized_BT_PV, & + "If true use the bottom depth instead of the total water column thickness "//& + "in the barotropic Coriolis term calculations.", default=.true.) call get_param(param_file, mdl, "BEBT", CS%bebt, & "BEBT determines whether the barotropic time stepping "//& "uses the forward-backward time-stepping scheme or a "//& @@ -4106,18 +4270,22 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, ALLOC_(CS%D_u_Cor(CS%isdw-1:CS%iedw,CS%jsdw:CS%jedw)) ALLOC_(CS%D_v_Cor(CS%isdw:CS%iedw,CS%jsdw-1:CS%jedw)) CS%q_D(:,:) = 0.0 ; CS%D_u_Cor(:,:) = 0.0 ; CS%D_v_Cor(:,:) = 0.0 + + Mean_SL = G%Z_ref do j=js,je ; do I=is-1,ie - CS%D_u_Cor(I,j) = 0.5 * (G%bathyT(i+1,j) + G%bathyT(i,j)) + CS%D_u_Cor(I,j) = 0.5 * (max(Mean_SL+G%bathyT(i+1,j),0.0) + max(Mean_SL+G%bathyT(i,j),0.0)) enddo ; enddo do J=js-1,je ; do i=is,ie - CS%D_v_Cor(i,J) = 0.5 * (G%bathyT(i,j+1) + G%bathyT(i,j)) + CS%D_v_Cor(i,J) = 0.5 * (max(Mean_SL+G%bathyT(i,j+1),0.0) + max(Mean_SL+G%bathyT(i,j),0.0)) enddo ; enddo do J=js-1,je ; do I=is-1,ie if (G%mask2dT(i,j)+G%mask2dT(i,j+1)+G%mask2dT(i+1,j)+G%mask2dT(i+1,j+1)>0.) then - CS%q_D(I,J) = 0.25 * G%CoriolisBu(I,J) * & + CS%q_D(I,J) = 0.25 * (CS%BT_Coriolis_scale * G%CoriolisBu(I,J)) * & ((G%areaT(i,j) + G%areaT(i+1,j+1)) + (G%areaT(i+1,j) + G%areaT(i,j+1))) / & - ((G%areaT(i,j) * G%bathyT(i,j) + G%areaT(i+1,j+1) * G%bathyT(i+1,j+1)) + & - (G%areaT(i+1,j) * G%bathyT(i+1,j) + G%areaT(i,j+1) * G%bathyT(i,j+1)) ) + (max(((G%areaT(i,j) * max(Mean_SL+G%bathyT(i,j),0.0) + & + G%areaT(i+1,j+1) * max(Mean_SL+G%bathyT(i+1,j+1),0.0)) + & + (G%areaT(i+1,j) * max(Mean_SL+G%bathyT(i+1,j),0.0) + & + G%areaT(i,j+1) * max(Mean_SL+G%bathyT(i,j+1),0.0))), GV%H_to_Z*GV%H_subroundoff) ) else ! All four h points are masked out so q_D(I,J) will is meaningless CS%q_D(I,J) = 0. endif @@ -4329,27 +4497,30 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, ! Calculate other constants which are used for btstep. - do j=js,je ; do I=is-1,ie - if (G%mask2dCu(I,j)>0.) then - CS%IDatu(I,j) = G%mask2dCu(I,j) * 2.0 / (G%bathyT(i+1,j) + G%bathyT(i,j)) - else ! Both neighboring H points are masked out so IDatu(I,j) is meaningless - CS%IDatu(I,j) = 0. - endif - enddo ; enddo - do J=js-1,je ; do i=is,ie - if (G%mask2dCv(i,J)>0.) then - CS%IDatv(i,J) = G%mask2dCv(i,J) * 2.0 / (G%bathyT(i,j+1) + G%bathyT(i,j)) - else ! Both neighboring H points are masked out so IDatv(I,j) is meaningless - CS%IDatv(i,J) = 0. - endif - enddo ; enddo + if (.not.CS%nonlin_stress) then + Mean_SL = G%Z_ref + do j=js,je ; do I=is-1,ie + if (G%mask2dCu(I,j)>0.) then + CS%IDatu(I,j) = G%mask2dCu(I,j) * 2.0 / ((G%bathyT(i+1,j) + G%bathyT(i,j)) + 2.0*Mean_SL) + else ! Both neighboring H points are masked out so IDatu(I,j) is meaningless + CS%IDatu(I,j) = 0. + endif + enddo ; enddo + do J=js-1,je ; do i=is,ie + if (G%mask2dCv(i,J)>0.) then + CS%IDatv(i,J) = G%mask2dCv(i,J) * 2.0 / ((G%bathyT(i,j+1) + G%bathyT(i,j)) + 2.0*Mean_SL) + else ! Both neighboring H points are masked out so IDatv(I,j) is meaningless + CS%IDatv(i,J) = 0. + endif + enddo ; enddo + endif call find_face_areas(Datu, Datv, G, GV, US, CS, MS, halo=1) - if (CS%bound_BT_corr) then - ! ### Consider replacing maxvel with G%dxT(i,j) * (CS%maxCFL_BT_cont*Idt) - ! ### and G%dyT(i,j) * (CS%maxCFL_BT_cont*Idt) + if ((CS%bound_BT_corr) .and. .not.(use_BT_Cont_type .and. CS%BT_cont_bounds)) then + ! This is not used in most test cases. Were it ever to become more widely used, consider + ! replacing maxvel with min(G%dxT(i,j),G%dyT(i,j)) * (CS%maxCFL_BT_cont*Idt) . do j=js,je ; do i=is,ie - CS%eta_cor_bound(i,j) = GV%m_to_H * G%IareaT(i,j) * 0.1 * CS%maxvel * & + CS%eta_cor_bound(i,j) = G%IareaT(i,j) * 0.1 * CS%maxvel * & ((Datu(I-1,j) + Datu(I,j)) + (Datv(i,J) + Datv(i,J-1))) enddo ; enddo endif diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 8c016b11b0..005f73af11 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -232,10 +232,9 @@ module MOM_dynamics_split_RK2 contains !> RK2 splitting for time stepping MOM adiabatic dynamics -subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & - Time_local, dt, forces, p_surf_begin, p_surf_end, & - uh, vh, uhtr, vhtr, eta_av, & - G, GV, US, CS, calc_dtbt, VarMix, MEKE, thickness_diffuse_CSp, Waves) +subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_surf_begin, p_surf_end, & + uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, calc_dtbt, VarMix, & + MEKE, thickness_diffuse_CSp, Waves) type(ocean_grid_type), intent(inout) :: G !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -272,8 +271,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & logical, intent(in) :: calc_dtbt !< if true, recalculate barotropic time step type(VarMix_CS), pointer :: VarMix !< specify the spatially varying viscosities type(MEKE_type), pointer :: MEKE !< related to mesoscale eddy kinetic energy param - type(thickness_diffuse_CS), pointer :: thickness_diffuse_CSp!< Pointer to a structure containing - !! interface height diffusivities + type(thickness_diffuse_CS), pointer :: thickness_diffuse_CSp !< Pointer to a structure containing + !! interface height diffusivities type(wave_parameters_CS), optional, pointer :: Waves !< A pointer to a structure containing !! fields related to the surface wave conditions @@ -954,7 +953,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param diag, CS, restart_CS, dt, Accel_diag, Cont_diag, MIS, & VarMix, MEKE, thickness_diffuse_CSp, & OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, & - visc, dirs, ntrunc, calc_dtbt) + visc, dirs, ntrunc, calc_dtbt, cont_stencil) type(ocean_grid_type), intent(inout) :: G !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -995,6 +994,8 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param !! the number of times the velocity is !! truncated (this should be 0). logical, intent(out) :: calc_dtbt !< If true, recalculate the barotropic time step + integer, optional, intent(out) :: cont_stencil !< The stencil for thickness + !! from the continuity solver. ! local variables real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_tmp @@ -1104,6 +1105,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param grain=CLOCK_ROUTINE) call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp) + if (present(cont_stencil)) cont_stencil = continuity_stencil(CS%continuity_CSp) call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv_CSp) if (use_tides) call tidal_forcing_init(Time, G, param_file, CS%tides_CSp) call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, & diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index ed7c440010..4030d0f2da 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -77,7 +77,7 @@ module MOM_dynamics_unsplit use MOM_ALE, only : ALE_CS use MOM_barotropic, only : barotropic_CS use MOM_boundary_update, only : update_OBC_data, update_OBC_CS -use MOM_continuity, only : continuity, continuity_init, continuity_CS +use MOM_continuity, only : continuity, continuity_init, continuity_CS, continuity_stencil use MOM_CoriolisAdv, only : CorAdCalc, CoriolisAdv_init, CoriolisAdv_CS use MOM_debugging, only : check_redundant use MOM_grid, only : ocean_grid_type @@ -561,7 +561,7 @@ end subroutine register_restarts_dyn_unsplit subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS, & restart_CS, Accel_diag, Cont_diag, MIS, MEKE, & OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, & - visc, dirs, ntrunc) + visc, dirs, ntrunc, cont_stencil) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -608,6 +608,8 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS integer, target, intent(inout) :: ntrunc !< A target for the variable that !! records the number of times the velocity !! is truncated (this should be 0). + integer, optional, intent(out) :: cont_stencil !< The stencil for thickness + !! from the continuity solver. ! This subroutine initializes all of the variables that are used by this ! dynamic core, including diagnostics and the cpu clocks. @@ -651,6 +653,7 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS Accel_diag%CAu => CS%CAu ; Accel_diag%CAv => CS%CAv call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp) + if (present(cont_stencil)) cont_stencil = continuity_stencil(CS%continuity_CSp) call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv_CSp) if (use_tides) call tidal_forcing_init(Time, G, param_file, CS%tides_CSp) call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, & diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index 98de5b931c..7700507301 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -76,7 +76,7 @@ module MOM_dynamics_unsplit_RK2 use MOM_ALE, only : ALE_CS use MOM_boundary_update, only : update_OBC_data, update_OBC_CS use MOM_barotropic, only : barotropic_CS -use MOM_continuity, only : continuity, continuity_init, continuity_CS +use MOM_continuity, only : continuity, continuity_init, continuity_CS, continuity_stencil use MOM_CoriolisAdv, only : CorAdCalc, CoriolisAdv_init, CoriolisAdv_CS use MOM_debugging, only : check_redundant use MOM_grid, only : ocean_grid_type @@ -506,7 +506,7 @@ end subroutine register_restarts_dyn_unsplit_RK2 subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, US, param_file, diag, CS, & restart_CS, Accel_diag, Cont_diag, MIS, MEKE, & OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, & - visc, dirs, ntrunc) + visc, dirs, ntrunc, cont_stencil) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -551,6 +551,8 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, US, param_file, diag integer, target, intent(inout) :: ntrunc !< A target for the variable !! that records the number of times the !! velocity is truncated (this should be 0). + integer, optional, intent(out) :: cont_stencil !< The stencil for + !! thickness from the continuity solver. ! This subroutine initializes all of the variables that are used by this ! dynamic core, including diagnostics and the cpu clocks. @@ -610,6 +612,7 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, US, param_file, diag Accel_diag%CAu => CS%CAu ; Accel_diag%CAv => CS%CAv call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp) + if (present(cont_stencil)) cont_stencil = continuity_stencil(CS%continuity_CSp) call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv_CSp) if (use_tides) call tidal_forcing_init(Time, G, param_file, CS%tides_CSp) call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, & diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90 index 1a2d03bd44..10832ffe75 100644 --- a/src/core/MOM_grid.F90 +++ b/src/core/MOM_grid.F90 @@ -138,7 +138,8 @@ module MOM_grid y_axis_units !< The units that are used in labeling the y coordinate axes. real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: & - bathyT !< Ocean bottom depth at tracer points, in depth units [Z ~> m]. + bathyT !< Ocean bottom depth at tracer points, in depth units [Z ~> m]. + real :: Z_ref !< A reference value for all geometric height fields, such as bathyT [Z ~> m]. logical :: bathymetry_at_vel !< If true, there are separate values for the !! basin depths at velocity points. Otherwise the effects of @@ -194,14 +195,16 @@ subroutine MOM_grid_init(G, param_file, US, HI, global_indexing, bathymetry_at_v !! velocity points. Otherwise the effects of topography !! are entirely determined from thickness points. -! This include declares and sets the variable "version". -#include "version_variable.h" + ! Local variables + real :: mean_SeaLev_scale integer :: isd, ied, jsd, jed, nk integer :: IsdB, IedB, JsdB, JedB integer :: ied_max, jed_max integer :: niblock, njblock, nihalo, njhalo, nblocks, n, i, j logical :: local_indexing ! If false use global index values instead of having ! the data domain on each processor start at 1. + ! This include declares and sets the variable "version". +# include "version_variable.h" integer, allocatable, dimension(:) :: ibegin, iend, jbegin, jend character(len=40) :: mod_nm = "MOM_grid" ! This module's name. @@ -218,9 +221,13 @@ subroutine MOM_grid_init(G, param_file, US, HI, global_indexing, bathymetry_at_v call get_param(param_file, mod_nm, "NJBLOCK", njblock, "The number of blocks "// & "in the y-direction on each processor (for openmp).", default=1, & layoutParam=.true.) - if (present(US)) then ; if (associated(US)) G%US => US ; endif + mean_SeaLev_scale = 1.0 ; if (associated(G%US)) mean_SeaLev_scale = G%US%m_to_Z + call get_param(param_file, mod_nm, "REFERENCE_HEIGHT", G%Z_ref, & + "A reference value for geometric height fields, such as bathyT.", & + units="m", default=0.0, scale=mean_SeaLev_scale) + if (present(HI)) then G%HI = HI diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 0af2b1759b..a6cd8c048a 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -100,7 +100,7 @@ end subroutine myStats !! valid data (good=1). If no information is available, !! Then use a previous guess (prev). Optionally (smooth) !! blend the filled points to achieve a more desirable result. -subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, keep_bug, debug) +subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, debug, answers_2018) use MOM_coms, only : sum_across_PEs type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. @@ -119,9 +119,10 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, integer, optional, intent(in) :: num_pass !< The maximum number of iterations real, optional, intent(in) :: relc !< A relaxation coefficient for Laplacian (ND) real, optional, intent(in) :: crit !< A minimal value for deltas between iterations. - logical, optional, intent(in) :: keep_bug !< Use an algorithm with a bug that dates - !! to the "sienna" code release. logical, optional, intent(in) :: debug !< If true, write verbose debugging messages. + logical, optional, intent(in) :: answers_2018 !< If true, use expressions that give the same + !! answers as the code did in late 2018. Otherwise + !! add parentheses for rotational symmetry. real, dimension(SZI_(G),SZJ_(G)) :: b,r real, dimension(SZI_(G),SZJ_(G)) :: fill_pts, good_, good_new @@ -138,7 +139,7 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, integer :: npass integer :: is, ie, js, je real :: relax_coeff, acrit, ares - logical :: debug_it + logical :: debug_it, ans_2018 debug_it=.false. if (PRESENT(debug)) debug_it=debug @@ -154,12 +155,11 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, acrit = crit_default if (PRESENT(crit)) acrit = crit - siena_bug=.false. - if (PRESENT(keep_bug)) siena_bug = keep_bug - do_smooth=.false. if (PRESENT(smooth)) do_smooth=smooth + ans_2018 = .true. ; if (PRESENT(answers_2018)) ans_2018 = answers_2018 + fill_pts(:,:) = fill(:,:) nfill = sum(fill(is:ie,js:je)) @@ -189,11 +189,17 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, if (gn == 1.0) north = aout(i,j+1)*gn if (gs == 1.0) south = aout(i,j-1)*gs - ngood = ge+gw+gn+gs + if (ans_2018) then + ngood = ge+gw+gn+gs + else + ngood = (ge+gw) + (gn+gs) + endif if (ngood > 0.) then - b(i,j)=(east+west+north+south)/ngood - !### Replace this with - ! b(i,j) = ((east+west) + (north+south))/ngood + if (ans_2018) then + b(i,j)=(east+west+north+south)/ngood + else + b(i,j) = ((east+west) + (north+south))/ngood + endif fill_pts(i,j) = 0.0 good_new(i,j) = 1.0 endif @@ -230,13 +236,15 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, if (fill(i,j) == 1) then east = max(good(i+1,j),fill(i+1,j)) ; west = max(good(i-1,j),fill(i-1,j)) north = max(good(i,j+1),fill(i,j+1)) ; south = max(good(i,j-1),fill(i,j-1)) - r(i,j) = relax_coeff*(south*aout(i,j-1)+north*aout(i,j+1) + & - west*aout(i-1,j)+east*aout(i+1,j) - & - (south+north+west+east)*aout(i,j)) - !### Appropriate parentheses should be added here, but they will change answers. - ! r(i,j) = relax_coeff*( ((south*aout(i,j-1) + north*aout(i,j+1)) + & - ! (west*aout(i-1,j)+east*aout(i+1,j))) - & - ! ((south+north)+(west+east))*aout(i,j) ) + if (ans_2018) then + r(i,j) = relax_coeff*(south*aout(i,j-1)+north*aout(i,j+1) + & + west*aout(i-1,j)+east*aout(i+1,j) - & + (south+north+west+east)*aout(i,j)) + else + r(i,j) = relax_coeff*( ((south*aout(i,j-1) + north*aout(i,j+1)) + & + (west*aout(i-1,j)+east*aout(i+1,j))) - & + ((south+north)+(west+east))*aout(i,j) ) + endif else r(i,j) = 0. endif @@ -264,7 +272,7 @@ end subroutine fill_miss_2d !> Extrapolate and interpolate from a file record subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, recnum, G, tr_z, & mask_z, z_in, z_edges_in, missing_value, reentrant_x, & - tripolar_n, homogenize, m_to_Z) + tripolar_n, homogenize, m_to_Z, answers_2018) character(len=*), intent(in) :: filename !< Path to file containing tracer to be !! interpolated. @@ -285,6 +293,9 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, !! to produce perfectly "flat" initial conditions real, optional, intent(in) :: m_to_Z !< A conversion factor from meters to the units !! of depth. If missing, G%bathyT must be in m. + logical, optional, intent(in) :: answers_2018 !< If true, use expressions that give the same + !! answers as the code did in late 2018. Otherwise + !! add parentheses for rotational symmetry. ! Local variables real, dimension(:,:), allocatable :: tr_in, tr_inp ! A 2-d array for holding input data on @@ -568,7 +579,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, good2(:,:) = good(:,:) fill2(:,:) = fill(:,:) - call fill_miss_2d(tr_outf, good2, fill2, tr_prev, G, smooth=.true.) + call fill_miss_2d(tr_outf, good2, fill2, tr_prev, G, smooth=.true., answers_2018=answers_2018) call myStats(tr_outf, missing_value, is, ie, js, je, k, 'field from fill_miss_2d()') tr_z(:,:,k) = tr_outf(:,:) * G%mask2dT(:,:) @@ -587,7 +598,7 @@ end subroutine horiz_interp_and_extrap_tracer_record !> Extrapolate and interpolate using a FMS time interpolation handle subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, tr_z, mask_z, & z_in, z_edges_in, missing_value, reentrant_x, & - tripolar_n, homogenize, spongeOngrid, m_to_Z) + tripolar_n, homogenize, spongeOngrid, m_to_Z, answers_2018) integer, intent(in) :: fms_id !< A unique id used by the FMS time interpolator type(time_type), intent(in) :: Time !< A FMS time type @@ -607,6 +618,9 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t logical, optional, intent(in) :: spongeOngrid !< If present and true, the sponge data are on the model grid real, optional, intent(in) :: m_to_Z !< A conversion factor from meters to the units !! of depth. If missing, G%bathyT must be in m. + logical, optional, intent(in) :: answers_2018 !< If true, use expressions that give the same + !! answers as the code did in late 2018. Otherwise + !! add parentheses for rotational symmetry. ! Local variables real, dimension(:,:), allocatable :: tr_in,tr_inp !< A 2-d array for holding input data on @@ -841,7 +855,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t good2(:,:) = good(:,:) fill2(:,:) = fill(:,:) - call fill_miss_2d(tr_outf, good2, fill2, tr_prev, G, smooth=.true.) + call fill_miss_2d(tr_outf, good2, fill2, tr_prev, G, smooth=.true., answers_2018=answers_2018) ! if (debug) then ! call hchksum(tr_outf, 'field from fill_miss_2d ', G%HI) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index d951be33c0..ed6aa5a44d 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -2034,7 +2034,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param type(remapping_CS) :: remapCS ! Remapping parameters and work arrays logical :: homogenize, useALEremapping, remap_full_column, remap_general, remap_old_alg - logical :: answers_2018, default_2018_answers + logical :: answers_2018, default_2018_answers, hor_regrid_answers_2018 logical :: use_ice_shelf character(len=10) :: remappingScheme real :: tempAvg, saltAvg @@ -2112,15 +2112,19 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param "If false, uses the preferred remapping algorithm for initialization. "//& "If true, use an older, less robust algorithm for remapping.", & default=.true., do_not_log=just_read) - if (useALEremapping) then - call get_param(PF, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, & + call get_param(PF, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, & "This sets the default value for the various _2018_ANSWERS parameters.", & default=.true.) + if (useALEremapping) then call get_param(PF, mdl, "REMAPPING_2018_ANSWERS", answers_2018, & "If true, use the order of arithmetic and expressions that recover the "//& "answers from the end of 2018. Otherwise, use updated and more robust "//& "forms of the same expressions.", default=default_2018_answers) endif + call get_param(PF, mdl, "HOR_REGRID_2018_ANSWERS", hor_regrid_answers_2018, & + "If true, use the order of arithmetic for horizonal regridding that recovers "//& + "the answers from the end of 2018. Otherwise, use rotationally symmetric "//& + "forms of the same expressions.", default=default_2018_answers) call get_param(PF, mdl, "ICE_SHELF", use_ice_shelf, default=.false.) if (use_ice_shelf) then call get_param(PF, mdl, "ICE_THICKNESS_FILE", ice_shelf_file, & @@ -2149,8 +2153,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param return ! All run-time parameters have been read, so return. endif - !### Change this to GV%Angstrom_Z - eps_z = 1.0e-10*US%m_to_Z + eps_z = GV%Angstrom_Z eps_rho = 1.0e-10*US%kg_m3_to_R ! Read input grid coordinates for temperature and salinity field @@ -2170,11 +2173,11 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param call horiz_interp_and_extrap_tracer(tfilename, potemp_var, 1.0, 1, & G, temp_z, mask_z, z_in, z_edges_in, missing_value_temp, reentrant_x, & - tripolar_n, homogenize, m_to_Z=US%m_to_Z) + tripolar_n, homogenize, m_to_Z=US%m_to_Z, answers_2018=hor_regrid_answers_2018) call horiz_interp_and_extrap_tracer(sfilename, salin_var, 1.0, 1, & G, salt_z, mask_z, z_in, z_edges_in, missing_value_salt, reentrant_x, & - tripolar_n, homogenize, m_to_Z=US%m_to_Z) + tripolar_n, homogenize, m_to_Z=US%m_to_Z, answers_2018=hor_regrid_answers_2018) kd = size(z_in,1) diff --git a/src/initialization/MOM_tracer_initialization_from_Z.F90 b/src/initialization/MOM_tracer_initialization_from_Z.F90 index bbe61892b2..5d585466c8 100644 --- a/src/initialization/MOM_tracer_initialization_from_Z.F90 +++ b/src/initialization/MOM_tracer_initialization_from_Z.F90 @@ -90,7 +90,7 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ real :: missing_value integer :: nPoints integer :: id_clock_routine, id_clock_ALE - logical :: answers_2018, default_2018_answers + logical :: answers_2018, default_2018_answers, hor_regrid_answers_2018 logical :: reentrant_x, tripolar_n id_clock_routine = cpu_clock_id('(Initialize tracer from Z)', grain=CLOCK_ROUTINE) @@ -112,15 +112,19 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ call get_param(PF, mdl, "Z_INIT_REMAPPING_SCHEME", remapScheme, & "The remapping scheme to use if using Z_INIT_ALE_REMAPPING is True.", & default="PLM") - if (useALE) then - call get_param(PF, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, & + call get_param(PF, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, & "This sets the default value for the various _2018_ANSWERS parameters.", & default=.true.) + if (useALE) then call get_param(PF, mdl, "REMAPPING_2018_ANSWERS", answers_2018, & "If true, use the order of arithmetic and expressions that recover the "//& "answers from the end of 2018. Otherwise, use updated and more robust "//& "forms of the same expressions.", default=default_2018_answers) endif + call get_param(PF, mdl, "HOR_REGRID_2018_ANSWERS", hor_regrid_answers_2018, & + "If true, use the order of arithmetic for horizonal regridding that recovers "//& + "the answers from the end of 2018. Otherwise, use rotationally symmetric "//& + "forms of the same expressions.", default=default_2018_answers) ! These are model grid properties, but being applied to the data grid for now. ! need to revisit this (mjh) @@ -137,7 +141,7 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ call horiz_interp_and_extrap_tracer(src_file, src_var_nam, convert, recnum, & G, tr_z, mask_z, z_in, z_edges_in, missing_value, reentrant_x, tripolar_n, & - homog, m_to_Z=US%m_to_Z) + homog, m_to_Z=US%m_to_Z, answers_2018=hor_regrid_answers_2018) kd = size(z_edges_in,1)-1 call pass_var(tr_z,G%Domain) diff --git a/src/initialization/midas_vertmap.F90 b/src/initialization/midas_vertmap.F90 index f33d476cf0..e6c586ed23 100644 --- a/src/initialization/midas_vertmap.F90 +++ b/src/initialization/midas_vertmap.F90 @@ -20,13 +20,10 @@ module MIDAS_vertmap !> Fill grid edges interface fill_boundaries - module procedure fill_boundaries_real - module procedure fill_boundaries_int + module procedure fill_boundaries_real + module procedure fill_boundaries_int end interface -! real, parameter :: epsln=1.e-10 !< A hard-wired constant! - !! \todo Get rid of this constant - contains #ifdef PY_SOLO @@ -239,7 +236,7 @@ function tracer_z_init(tr_in, z_edges, e, nkml, nkbl, land_fill, wet, nlay, nlev if (kz /= k_bot_prev) then ! Calculate the intra-cell profile. if ((kz < nlevs_data(i,j)) .and. (kz > 1)) then - sl_tr = find_limited_slope(tr_1d, z_edges, kz) + sl_tr = find_limited_slope(tr_1d, z_edges, kz) endif endif if (kz > nlevs_data(i,j)) kz = nlevs_data(i,j) @@ -264,7 +261,7 @@ function tracer_z_init(tr_in, z_edges, e, nkml, nkbl, land_fill, wet, nlay, nlev ! Calculate the intra-cell profile. sl_tr = 0.0 ! ; cur_tr = 0.0 if ((kz < nlevs_data(i,j)) .and. (kz > 1)) then - sl_tr = find_limited_slope(tr_1d, z_edges, kz) + sl_tr = find_limited_slope(tr_1d, z_edges, kz) endif ! This is the piecewise linear form. tr(i,j,k) = tr(i,j,k) + wt(kz) * & @@ -378,14 +375,22 @@ subroutine determine_temperature(temp, salt, R, p_ref, niter, land_fill, h, k_st real(kind=8), dimension(size(temp,1),size(temp,3)) :: drho_dT, drho_dS real(kind=8), dimension(size(temp,1)) :: press integer :: nx, ny, nz, nt, i, j, k, n, itt - real :: dT_dS + real :: dT_dS_gauge ! The relative penalizing of temperature to salinity changes when + ! minimizing property changes while correcting density [degC ppt-1]. + real :: I_denom ! The inverse of the magnitude squared of the density gradient in + ! T-S space streched with dT_dS_gauge [m6 kg-2 ppt-1] logical :: adjust_salt, old_fit real, parameter :: S_min = 0.5, S_max=65.0 - real, parameter :: tol=1.e-4, max_t_adj=1.0, max_s_adj = 0.5 + real, parameter :: tol_T=1.e-4, tol_S=1.e-4, tol_rho=1.e-4 + real, parameter :: max_t_adj=1.0, max_s_adj = 0.5 old_fit = .true. ! reproduces siena behavior - ! will switch to the newer method which simultaneously adjusts - ! temp and salt based on the ratio of the thermal and haline coefficients. + + ! ### The whole determine_temperature subroutine needs to be reexamined, both the algorithms + ! and the extensive use of hard-coded dimensional parameters. + + ! We will switch to the newer method which simultaneously adjusts + ! temp and salt based on the ratio of the thermal and haline coefficients, once it is tested. nx=size(temp,1) ; ny=size(temp,2) ; nz=size(temp,3) @@ -411,23 +416,22 @@ subroutine determine_temperature(temp, salt, R, p_ref, niter, land_fill, h, k_st do k=k_start,nz ; do i=1,nx ! if (abs(rho(i,k)-R(k))>tol .and. hin(i,k)>epsln .and. abs(T(i,k)-land_fill) < epsln) then - if (abs(rho(i,k)-R(k))>tol) then + if (abs(rho(i,k)-R(k))>tol_rho) then if (old_fit) then - dT(i,k) = max(min((R(k)-rho(i,k)) / drho_dT(i,k), max_t_adj), -max_t_adj) - T(i,k) = max(min(T(i,k)+dT(i,k), T_max), T_min) + dT(i,k) = max(min((R(k)-rho(i,k)) / drho_dT(i,k), max_t_adj), -max_t_adj) + T(i,k) = max(min(T(i,k)+dT(i,k), T_max), T_min) else - dT_dS = 10.0 - min(-drho_dT(i,k)/drho_dS(i,k),10.) - !### RWH: Based on the dimensions alone, the expression above should be: - ! dT_dS = 10.0 - min(-drho_dS(i,k)/drho_dT(i,k),10.) - dS(i,k) = (R(k)-rho(i,k)) / (drho_dS(i,k) - drho_dT(i,k)*dT_dS ) - dT(i,k) = -dT_dS*dS(i,k) - ! dT(i,k) = max(min(dT(i,k), max_t_adj), -max_t_adj) - T(i,k) = max(min(T(i,k)+dT(i,k), T_max), T_min) - S(i,k) = max(min(S(i,k)+dS(i,k), S_max), S_min) + dT_dS_gauge = 10.0 ! 10 degC is weighted equivalently to 1 ppt. + I_denom = 1.0 / (drho_dS(i,k)**2 + dT_dS_gauge**2*drho_dT(i,k)**2) + dS(i,k) = (R(k)-rho(i,k)) * drho_dS(i,k) * I_denom + dT(i,k) = (R(k)-rho(i,k)) * dT_dS_gauge**2*drho_dT(i,k) * I_denom + + T(i,k) = max(min(T(i,k)+dT(i,k), T_max), T_min) + S(i,k) = max(min(S(i,k)+dS(i,k), S_max), S_min) endif endif enddo ; enddo - if (maxval(abs(dT)) < tol) then + if (maxval(abs(dT)) < tol_T) then adjust_salt = .false. exit iter_loop endif @@ -445,12 +449,12 @@ subroutine determine_temperature(temp, salt, R, p_ref, niter, land_fill, h, k_st #endif do k=k_start,nz ; do i=1,nx ! if (abs(rho(i,k)-R(k))>tol .and. hin(i,k)>epsln .and. abs(T(i,k)-land_fill) < epsln ) then - if (abs(rho(i,k)-R(k)) > tol) then + if (abs(rho(i,k)-R(k)) > tol_rho) then dS(i,k) = max(min((R(k)-rho(i,k)) / drho_dS(i,k), max_s_adj), -max_s_adj) S(i,k) = max(min(S(i,k)+dS(i,k), S_max), S_min) endif enddo ; enddo - if (maxval(abs(dS)) < tol) exit + if (maxval(abs(dS)) < tol_S) exit enddo ; endif temp(:,j,:)=T(:,:) @@ -480,13 +484,12 @@ subroutine find_overlap(e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z real :: Ih, e_c, tot_wt, I_totwt integer :: k - wt(:)=0.0 ; z1(:)=0.0 ; z2(:)=0.0 - k_top = k_start ; k_bot = k_start ; wt(1) = 1.0 ; z1(1) = -0.5 ; z2(1) = 0.5 + wt(:) = 0.0 ; z1(:) = 0.0 ; z2(:) = 0.0 ; k_bot = k_max + wt(1) = 1.0 ; z1(1) = -0.5 ; z2(1) = 0.5 do k=k_start,k_max ; if (e(K+1) < Z_top) exit ; enddo k_top = k - - if (k>k_max) return + if (k_top > k_max) return ! Determine the fractional weights of each layer. ! Note that by convention, e and Z_int decrease with increasing k. @@ -498,7 +501,6 @@ subroutine find_overlap(e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z z2(k) = (e_c - Z_bot) * Ih else wt(k) = MIN(e(K),Z_top) - e(K+1) ; tot_wt = wt(k) ! These are always > 0. - ! Ih = 0.0 ; if (e(K) /= e(K+1)) Ih = 1.0 / (e(K)-e(K+1)) if (e(K) /= e(K+1)) then z1(k) = (0.5*(e(K)+e(K+1)) - MIN(e(K), Z_top)) / (e(K)-e(K+1)) else ; z1(k) = -0.5 ; endif @@ -524,36 +526,27 @@ subroutine find_overlap(e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z end subroutine find_overlap -!> This subroutine determines a limited slope for val to be advected with +!> This function determines a limited slope for val to be advected with !! a piecewise limited scheme. function find_limited_slope(val, e, k) result(slope) - real, dimension(:), intent(in) :: val !< An column the values that are being interpolated. + real, dimension(:), intent(in) :: val !< A column of values that are being interpolated, in arbitrary units [A]. real, dimension(:), intent(in) :: e !< A column's interface heights [Z ~> m] or other units. integer, intent(in) :: k !< The layer whose slope is being determined. - real :: slope !< The normalized slope in the intracell distribution of val. + real :: slope !< The normalized slope in the intracell distribution of val [A Z-1 ~> A m-1] or other units. ! Local variables - real :: amn, cmn - real :: d1, d2 + real :: d1, d2 ! Thicknesses in the units of e [Z ~> m]. - if ((val(k)-val(k-1)) * (val(k)-val(k+1)) >= 0.0) then - slope = 0.0 ! ; curvature = 0.0 + d1 = 0.5*(e(K-1)-e(K+1)) ; d2 = 0.5*(e(K)-e(K+2)) + if (((val(k)-val(k-1)) * (val(k)-val(k+1)) >= 0.0) .or. (d1*d2 <= 0.0)) then + slope = 0.0 else - d1 = 0.5*(e(K-1)-e(K+1)) ; d2 = 0.5*(e(K)-e(K+2)) - if (d1*d2 > 0.0) then - slope = ((d1**2)*(val(k+1) - val(k)) + (d2**2)*(val(k) - val(k-1))) * & - (e(K) - e(K+1)) / (d1*d2*(d1+d2)) - ! slope = 0.5*(val(k+1) - val(k-1)) - ! This is S.J. Lin's form of the PLM limiter. - amn = min(abs(slope), 2.0*(max(val(k-1), val(k), val(k+1)) - val(k))) - cmn = 2.0*(val(k) - min(val(k-1), val(k), val(k+1))) - slope = sign(1.0, slope) * min(amn, cmn) - - ! min(abs(slope), 2.0*(max(val(k-1),val(k),val(k+1)) - val(k)), & - ! 2.0*(val(k) - min(val(k-1),val(k),val(k+1)))) - ! curvature = 0.0 - else - slope = 0.0 ! ; curvature = 0.0 - endif + slope = ((d1**2)*(val(k+1) - val(k)) + (d2**2)*(val(k) - val(k-1))) * & + (e(K) - e(K+1)) / (d1*d2*(d1+d2)) + ! slope = 0.5*(val(k+1) - val(k-1)) + ! This is S.J. Lin's form of the PLM limiter. + slope = sign(1.0, slope) * min(abs(slope), & + 2.0*(max(val(k-1), val(k), val(k+1)) - val(k)), & + 2.0*(val(k) - min(val(k-1), val(k), val(k+1)))) endif end function find_limited_slope diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 63811e14d7..0298bac5ab 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -67,6 +67,8 @@ module MOM_hor_visc !! viscosity is modified to include a term that !! scales quadratically with the velocity shears. logical :: use_Kh_bg_2d !< Read 2d background viscosity from a file. + logical :: Kh_bg_2d_bug !< If true, retain an answer-changing horizontal indexing bug + !! in setting the corner-point viscosities when USE_KH_BG_2D=True. real :: Kh_bg_min !< The minimum value allowed for Laplacian horizontal !! viscosity [L2 T-1 ~> m2 s-1]. The default is 0.0. logical :: use_land_mask !< Use the land mask for the computation of thicknesses @@ -1651,6 +1653,11 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) "If true, read a file containing 2-d background harmonic "//& "viscosities. The final viscosity is the maximum of the other "//& "terms and this background value.", default=.false.) + if (CS%use_Kh_bg_2d) then + call get_param(param_file, mdl, "KH_BG_2D_BUG", CS%Kh_bg_2d_bug, & + "If true, retain an answer-changing horizontal indexing bug in setting "//& + "the corner-point viscosities when USE_KH_BG_2D=True.", default=.true.) + endif call get_param(param_file, mdl, "USE_GME", CS%use_GME, & "If true, use the GM+E backscatter scheme in association \n"//& @@ -1872,8 +1879,14 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) CS%Kh_bg_xy(I,J) = MAX(Kh, Kh_vel_scale * sqrt(grid_sp_q2)) ! Use the larger of the above and values read from a file - !### This expression uses inconsistent staggering - if (CS%use_Kh_bg_2d) CS%Kh_bg_xy(I,J) = MAX(CS%Kh_bg_2d(i,j), CS%Kh_bg_xy(I,J)) + if (CS%use_Kh_bg_2d) then ; if (CS%Kh_bg_2d_bug) then + ! This option is unambiguously wrong, and should be obsoleted as soon as possible. + CS%Kh_bg_xy(I,J) = MAX(CS%Kh_bg_2d(i,j), CS%Kh_bg_xy(I,J)) + else + CS%Kh_bg_xy(I,J) = MAX(CS%Kh_bg_xy(I,J), & + 0.25*((CS%Kh_bg_2d(i,j) + CS%Kh_bg_2d(i+1,j+1)) + & + (CS%Kh_bg_2d(i+1,j) + CS%Kh_bg_2d(i,j+1))) ) + endif ; endif ! Use the larger of the above and a function of sin(latitude) if (Kh_sin_lat>0.) then diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 9c236fd937..450d6fea9b 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -69,6 +69,9 @@ module MOM_thickness_diffuse !! the GEOMETRIC thickness difussion [nondim] real :: MEKE_GEOMETRIC_epsilon !< Minimum Eady growth rate for the GEOMETRIC thickness !! diffusivity [T-1 ~> s-1]. + logical :: MEKE_GEOM_answers_2018 !< If true, use expressions in the MEKE_GEOMETRIC calculation + !! that recover the answers from the original implementation. + !! Otherwise, use expressions that satisfy rotational symmetry. logical :: Use_KH_in_MEKE !< If true, uses the thickness diffusivity calculated here to diffuse MEKE. 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. @@ -151,7 +154,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp real :: KH_u_lay(SZI_(G), SZJ_(G)) ! layer ave thickness diffusivities [L2 T-1 ~> m2 s-1] real :: KH_v_lay(SZI_(G), SZJ_(G)) ! layer ave thickness diffusivities [L2 T-1 ~> m2 s-1] - if (.not. associated(CS)) call MOM_error(FATAL, "MOM_thickness_diffuse:"// & + if (.not. associated(CS)) call MOM_error(FATAL, "MOM_thickness_diffuse: "//& "Module must be initialized before it is used.") if ((.not.CS%thickness_diffuse) .or. & @@ -366,13 +369,25 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp if (associated(MEKE)) then ; if (associated(MEKE%Kh)) then if (CS%MEKE_GEOMETRIC) then -!$OMP do - do j=js,je ; do I=is,ie - !### This will not give bitwise rotational symmetry. Add parentheses. - MEKE%Kh(i,j) = CS%MEKE_GEOMETRIC_alpha * MEKE%MEKE(i,j) / & - (0.25*(VarMix%SN_u(I,j)+VarMix%SN_u(I-1,j)+VarMix%SN_v(i,J)+VarMix%SN_v(i,J-1)) + & - CS%MEKE_GEOMETRIC_epsilon) - enddo ; enddo + if (CS%MEKE_GEOM_answers_2018) then + !$OMP do + do j=js,je ; do I=is,ie + ! This does not give bitwise rotational symmetry. + MEKE%Kh(i,j) = CS%MEKE_GEOMETRIC_alpha * MEKE%MEKE(i,j) / & + (0.25*(VarMix%SN_u(I,j)+VarMix%SN_u(I-1,j) + & + VarMix%SN_v(i,J)+VarMix%SN_v(i,J-1)) + & + CS%MEKE_GEOMETRIC_epsilon) + enddo ; enddo + else + !$OMP do + do j=js,je ; do I=is,ie + ! With the additional parentheses this gives bitwise rotational symmetry. + MEKE%Kh(i,j) = CS%MEKE_GEOMETRIC_alpha * MEKE%MEKE(i,j) / & + (0.25*((VarMix%SN_u(I,j)+VarMix%SN_u(I-1,j)) + & + (VarMix%SN_v(i,J)+VarMix%SN_v(i,J-1))) + & + CS%MEKE_GEOMETRIC_epsilon) + enddo ; enddo + endif endif endif ; endif @@ -1768,7 +1783,10 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) #include "version_variable.h" character(len=40) :: mdl = "MOM_thickness_diffuse" ! This module's name. real :: omega ! The Earth's rotation rate [T-1 ~> s-1] - real :: strat_floor + real :: strat_floor ! A floor for Brunt-Vasaila frequency in the Ferrari et al. 2010, + ! streamfunction formulation, expressed as a fraction of planetary + ! rotation [nondim]. + logical :: default_2018_answers ! The default setting for the various 2018_ANSWERS flags. if (associated(CS)) then call MOM_error(WARNING, & @@ -1852,32 +1870,38 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) default=.false., debuggingParam=.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 \n"//& + "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_GEOMETRIC", CS%MEKE_GEOMETRIC, & - "If true, uses the GM coefficient formulation \n"//& - "from the GEOMETRIC framework (Marshall et al., 2012).", default=.false.) + "If true, uses the GM coefficient formulation from the GEOMETRIC "//& + "framework (Marshall et al., 2012).", default=.false.) if (CS%MEKE_GEOMETRIC) then - call get_param(param_file, mdl, "MEKE_GEOMETRIC_EPSILON", CS%MEKE_GEOMETRIC_epsilon, & - "Minimum Eady growth rate used in the calculation of \n"//& - "GEOMETRIC thickness diffusivity.", units="s-1", default=1.0e-7, scale=US%T_to_s) - + "Minimum Eady growth rate used in the calculation of GEOMETRIC "//& + "thickness diffusivity.", units="s-1", default=1.0e-7, scale=US%T_to_s) call get_param(param_file, mdl, "MEKE_GEOMETRIC_ALPHA", CS%MEKE_GEOMETRIC_alpha, & - "The nondimensional coefficient governing the efficiency of the GEOMETRIC \n"//& + "The nondimensional coefficient governing the efficiency of the GEOMETRIC "//& "thickness diffusion.", units="nondim", default=0.05) + + call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, & + "This sets the default value for the various _2018_ANSWERS parameters.", & + default=.true.) + call get_param(param_file, mdl, "MEKE_GEOMETRIC_2018_ANSWERS", CS%MEKE_GEOM_answers_2018, & + "If true, use expressions in the MEKE_GEOMETRIC calculation that recover the "//& + "answers from the original implementation. Otherwise, use expressions that "//& + "satisfy rotational symmetry.", default=default_2018_answers) endif call get_param(param_file, mdl, "USE_KH_IN_MEKE", CS%Use_KH_in_MEKE, & - "If true, uses the thickness diffusivity calculated here to diffuse \n"//& - "MEKE.", default=.false.) + "If true, uses the thickness diffusivity calculated here to diffuse MEKE.", & + default=.false.) call get_param(param_file, mdl, "USE_GME", CS%use_GME_thickness_diffuse, & - "If true, use the GM+E backscatter scheme in association \n"//& + "If true, use the GM+E backscatter scheme in association "//& "with the Gent and McWilliams parameterization.", default=.false.) call get_param(param_file, mdl, "USE_GM_WORK_BUG", CS%use_GM_work_bug, & - "If true, compute the top-layer work tendency on the u-grid " // & + "If true, compute the top-layer work tendency on the u-grid "//& "with the incorrect sign, for legacy reproducibility.", & default=.true.) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index b19f81f8d1..ccd85280f5 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -130,6 +130,9 @@ module MOM_ALE_sponge logical :: remap_answers_2018 !< If true, use the order of arithmetic and expressions that !! recover the answers for remapping from the end of 2018. !! Otherwise, use more robust forms of the same expressions. + logical :: hor_regrid_answers_2018 !< If true, use the order of arithmetic for horizonal regridding + !! that recovers the answers from the end of 2018. Otherwise, use + !! rotationally symmetric forms of the same expressions. logical :: time_varying_sponges !< True if using newer sponge code logical :: spongeDataOngrid !< True if the sponge data are on the model horizontal grid @@ -204,6 +207,10 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ "If true, use the order of arithmetic and expressions that recover the "//& "answers from the end of 2018. Otherwise, use updated and more robust "//& "forms of the same expressions.", default=default_2018_answers) + call get_param(param_file, mdl, "HOR_REGRID_2018_ANSWERS", CS%hor_regrid_answers_2018, & + "If true, use the order of arithmetic for horizonal regridding that recovers "//& + "the answers from the end of 2018. Otherwise, use rotationally symmetric "//& + "forms of the same expressions.", default=default_2018_answers) CS%time_varying_sponges = .false. CS%nz = G%ke @@ -742,15 +749,17 @@ subroutine set_up_ALE_sponge_vel_field_varying(filename_u, fieldname_u, filename ! I am hard-wiring this call to assume that the input grid is zonally re-entrant ! In the future, this should be generalized using an interface to return the ! modulo attribute of the zonal axis (mjh). - call horiz_interp_and_extrap_tracer(CS%Ref_val_u%id,Time, 1.0,G,u_val,mask_u,z_in,z_edges_in,& - missing_value,.true.,.false.,.false., m_to_Z=US%m_to_Z) + call horiz_interp_and_extrap_tracer(CS%Ref_val_u%id, Time, 1.0, G, u_val, mask_u, z_in, z_edges_in, & + missing_value, .true., .false., .false., m_to_Z=US%m_to_Z, & + answers_2018=CS%hor_regrid_answers_2018) !!! TODO: add a velocity interface! (mjh) ! Interpolate external file data to the model grid ! I am hard-wiring this call to assume that the input grid is zonally re-entrant ! In the future, this should be generalized using an interface to return the ! modulo attribute of the zonal axis (mjh). - call horiz_interp_and_extrap_tracer(CS%Ref_val_v%id,Time, 1.0,G,v_val,mask_v,z_in,z_edges_in, & - missing_value,.true.,.false.,.false., m_to_Z=US%m_to_Z) + call horiz_interp_and_extrap_tracer(CS%Ref_val_v%id, Time, 1.0, G, v_val, mask_v, z_in, z_edges_in, & + missing_value, .true., .false., .false., m_to_Z=US%m_to_Z, & + answers_2018=CS%hor_regrid_answers_2018) ! stores the reference profile allocate(CS%Ref_val_u%p(fld_sz(3),CS%num_col_u)) CS%Ref_val_u%p(:,:) = 0.0 @@ -824,8 +833,10 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) allocate(mask_z(G%isd:G%ied,G%jsd:G%jed,1:nz_data)) sp_val(:,:,:)=0.0 mask_z(:,:,:)=0.0 - call horiz_interp_and_extrap_tracer(CS%Ref_val(m)%id,Time, 1.0,G,sp_val,mask_z,z_in,z_edges_in, & - missing_value,.true., .false.,.false.,spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z) + call horiz_interp_and_extrap_tracer(CS%Ref_val(m)%id, Time, 1.0, G, sp_val, mask_z, z_in, & + z_edges_in, missing_value, .true., .false., .false., & + spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z, & + answers_2018=CS%hor_regrid_answers_2018) allocate( hsrc(nz_data) ) allocate( tmpT1d(nz_data) ) do c=1,CS%num_col @@ -906,8 +917,9 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) allocate(sp_val(G%isdB:G%iedB,G%jsd:G%jed,1:nz_data)) allocate(mask_z(G%isdB:G%iedB,G%jsd:G%jed,1:nz_data)) ! Interpolate from the external horizontal grid and in time - call horiz_interp_and_extrap_tracer(CS%Ref_val_u%id,Time, 1.0,G,sp_val,mask_z,z_in,z_edges_in, & - missing_value, .true., .false., .false., m_to_Z=US%m_to_Z) + call horiz_interp_and_extrap_tracer(CS%Ref_val_u%id, Time, 1.0, G, sp_val, mask_z, z_in, & + z_edges_in, missing_value, .true., .false., .false., & + m_to_Z=US%m_to_Z, answers_2018=CS%hor_regrid_answers_2018) ! call pass_var(sp_val,G%Domain) ! call pass_var(mask_z,G%Domain) do c=1,CS%num_col @@ -923,8 +935,9 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) allocate(sp_val(G%isd:G%ied,G%jsdB:G%jedB,1:nz_data)) allocate(mask_z(G%isd:G%ied,G%jsdB:G%jedB,1:nz_data)) ! Interpolate from the external horizontal grid and in time - call horiz_interp_and_extrap_tracer(CS%Ref_val_v%id,Time, 1.0,G,sp_val,mask_z,z_in,z_edges_in, & - missing_value, .true., .false., .false., m_to_Z=US%m_to_Z) + call horiz_interp_and_extrap_tracer(CS%Ref_val_v%id, Time, 1.0, G, sp_val, mask_z, z_in, & + z_edges_in, missing_value, .true., .false., .false., & + m_to_Z=US%m_to_Z, answers_2018=CS%hor_regrid_answers_2018) ! call pass_var(sp_val,G%Domain) ! call pass_var(mask_z,G%Domain) diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 index 0cbe700518..57199f38d0 100644 --- a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 +++ b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 @@ -464,8 +464,7 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd_lay, Kv, j, G, GV, US, CS) do i=is,ie bckgrnd_vdc_psis = CS%bckgrnd_vdc_psim * exp(-(0.4*(G%geoLatT(i,j)+28.9))**2) bckgrnd_vdc_psin = CS%bckgrnd_vdc_psim * exp(-(0.4*(G%geoLatT(i,j)-28.9))**2) - !### Add parentheses. - CS%Kd_bkgnd(i,j,:) = CS%bckgrnd_vdc_eq + bckgrnd_vdc_psin + bckgrnd_vdc_psis + CS%Kd_bkgnd(i,j,:) = (CS%bckgrnd_vdc_eq + bckgrnd_vdc_psin) + bckgrnd_vdc_psis if (G%geoLatT(i,j) < -10.0) then CS%Kd_bkgnd(i,j,:) = CS%Kd_bkgnd(i,j,:) + CS%bckgrnd_vdc1 diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index 2625867849..c910433172 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -111,6 +111,8 @@ module MOM_bulk_mixed_layer !! using SST for temperature of liq_runoff logical :: use_calving_heat_content !< Use SST for temperature of froz_runoff logical :: salt_reject_below_ML !< It true, add salt below mixed layer (layer mode only) + logical :: convect_mom_bug !< If true, use code with a bug that causes a loss of momentum + !! conservation during mixedlayer convection. type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to regulate the !! timing of diagnostic output. @@ -1123,6 +1125,9 @@ subroutine mixedlayer_convection(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, & T_precip * netMassIn(i) * GV%H_to_RZ * fluxes%C_p * Idt if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + & T_precip * netMassIn(i) * GV%H_to_RZ + else ! This is a massless column, but zero out the summed variables anyway for safety. + htot(i) = 0.0 ; Ttot(i) = 0.0 ; Stot(i) = 0.0 ; R0_tot(i) = 0.0 ; Rcv_tot = 0.0 + uhtot(i) = 0.0 ; vhtot(i) = 0.0 ; Conv_En(i) = 0.0 ; dKE_FC(i) = 0.0 endif ; enddo ! Now do netMassOut case in this block. @@ -1288,9 +1293,11 @@ subroutine mixedlayer_convection(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, & htot(i) = htot(i) + h_ent h(i,k) = h(i,k) - h_ent d_eb(i,k) = d_eb(i,k) - h_ent - uhtot(i) = u(i,k)*h_ent ; vhtot(i) = v(i,k)*h_ent - !### I think that the line above should instead be: - ! uhtot(i) = uhtot(i) + h_ent*u(i,k) ; vhtot(i) = vhtot(i) + h_ent*v(i,k) + if (CS%convect_mom_bug) then + uhtot(i) = u(i,k)*h_ent ; vhtot(i) = v(i,k)*h_ent + else + uhtot(i) = uhtot(i) + h_ent*u(i,k) ; vhtot(i) = vhtot(i) + h_ent*v(i,k) + endif endif @@ -3568,6 +3575,9 @@ subroutine bulkmixedlayer_init(Time, G, GV, US, param_file, diag, CS) "If true, use the fluxes%calving_Hflx field to set the "//& "heat carried by runoff, instead of using SST*CP*froz_runoff.", & default=.false.) + call get_param(param_file, mdl, "BULKML_CONV_MOMENTUM_BUG", CS%convect_mom_bug, & + "If true, use code with a bug that causes a loss of momentum conservation "//& + "during mixedlayer convection.", default=.true.) call get_param(param_file, mdl, "ALLOW_CLOCKS_IN_OMP_LOOPS", & CS%allow_clocks_in_omp_loops, & diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index f8c20682ee..962dcb455e 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -709,9 +709,6 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs real :: dMKE_src_dK ! The partial derivative of MKE_src with Kddt_h(K) [R Z3 T-2 H-1 ~> J m-3 or J kg-1]. real :: Kd_guess0 ! A first guess of the diapycnal diffusivity [Z2 T-1 ~> m2 s-1]. real :: PE_chg_g0 ! The potential energy change when Kd is Kd_guess0 [R Z3 T-2 ~> J m-2] - !### The following might be unused. - real :: dPEa_dKd_g0 ! The derivative of the change in the potential energy of the column above an interface - ! with the diffusivity when the Kd is Kd_guess0 [R Z T-1 ~> J s m-4] real :: Kddt_h_g0 ! The first guess diapycnal diffusivity times a timestep divided ! by the average thicknesses around a layer [H ~> m or kg m-2]. real :: PE_chg_max ! The maximum PE change for very large values of Kddt_h(K) [R Z3 T-2 ~> J m-2]. @@ -1135,16 +1132,14 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), & pres_Z(K), dT_to_dColHt(k), dS_to_dColHt(k), & dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & - PE_chg=PE_chg_g0, dPEc_dKd=dPEa_dKd_g0, dPE_max=PE_chg_max, & - dPEc_dKd_0=dPEc_dKd_Kd0 ) + PE_chg=PE_chg_g0, dPE_max=PE_chg_max, dPEc_dKd_0=dPEc_dKd_Kd0 ) else call find_PE_chg(0.0, Kddt_h_g0, hp_a, h(k), & Th_a(k-1), Sh_a(k-1), Th_b(k), Sh_b(k), & dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE(k), dS_to_dPE(k), & pres_Z(K), dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & dT_to_dColHt(k), dS_to_dColHt(k), & - PE_chg=PE_chg_g0, dPEc_dKd=dPEa_dKd_g0, dPE_max=PE_chg_max, & - dPEc_dKd_0=dPEc_dKd_Kd0 ) + PE_chg=PE_chg_g0, dPE_max=PE_chg_max, dPEc_dKd_0=dPEc_dKd_Kd0 ) endif MKE_src = dMKE_max*(1.0 - exp(-Kddt_h_g0 * MKE2_Hharm)) @@ -1821,7 +1816,7 @@ subroutine find_mstar(CS, US, Buoyancy_Flux, UStar, UStar_Mean,& MStar = MStar * MStar_Conv_Red if (present(Langmuir_Number)) then - !### In this call, ustar was previously ustar_mean. Is this change deliberate? + !### In this call, ustar was previously ustar_mean. Is this change deliberate, Brandon? -RWH call mstar_Langmuir(CS, US, Abs_Coriolis, Buoyancy_Flux, UStar, BLD, Langmuir_Number, MStar, & MStar_LT, Convect_Langmuir_Number) endif diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90 index 9349cf06d7..b6eba22e14 100644 --- a/src/parameterizations/vertical/MOM_kappa_shear.F90 +++ b/src/parameterizations/vertical/MOM_kappa_shear.F90 @@ -78,6 +78,15 @@ module MOM_kappa_shear ! I can think of no good reason why this should be false. - RWH real :: vel_underflow !< Velocity components smaller than vel_underflow !! are set to 0 [L T-1 ~> m s-1]. + real :: kappa_src_max_chg !< The maximum permitted increase in the kappa source within an + !! iteration relative to the local source [nondim]. This must be + !! greater than 1. The lower limit for the permitted fractional + !! decrease is (1 - 0.5/kappa_src_max_chg). These limits could + !! perhaps be made dynamic with an improved iterative solver. + logical :: all_layer_TKE_bug !< If true, report back the latest estimate of TKE instead of the + !! time average TKE when there is mass in all layers. Otherwise always + !! report the time-averaged TKE, as is currently done when there + !! are some massless layers. ! logical :: layer_stagger = .false. ! If true, do the calculations centered at ! layers, rather than the interfaces. logical :: debug = .false. !< If true, write verbose debugging messages. @@ -301,8 +310,11 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & if (nz == nzc) then do K=1,nz+1 kappa_2d(i,K) = kappa_avg(K) - !### Should this be tke_avg? - tke_2d(i,K) = tke(K) + if (CS%all_layer_TKE_bug) then + tke_2d(i,K) = tke(K) + else + tke_2d(i,K) = tke_avg(K) + endif enddo else do K=1,nz+1 @@ -599,8 +611,11 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ if (nz == nzc) then do K=1,nz+1 kappa_2d(I,K,J2) = kappa_avg(K) - !### Should this be tke_avg? - tke_2d(I,K) = tke(K) + if (CS%all_layer_TKE_bug) then + tke_2d(i,K) = tke(K) + else + tke_2d(i,K) = tke_avg(K) + endif enddo else do K=1,nz+1 @@ -752,7 +767,7 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & real :: g_R0 ! g_R0 is a rescaled version of g/Rho [Z R-1 T-2 ~> m4 kg-1 s-2]. real :: Norm ! A factor that normalizes two weights to 1 [Z-2 ~> m-2]. real :: tol_dksrc ! Tolerance for the change in the kappa source within an iteration - ! relative to the local source [nondim]. + ! relative to the local source [nondim]. This must be greater than 1. real :: tol2 ! The tolerance for the change in the kappa source within an iteration ! relative to the average local source over previous iterations [nondim]. real :: tol_dksrc_low ! The tolerance for the fractional decrease in ksrc @@ -801,9 +816,15 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & gR0 = GV%z_to_H*GV%H_to_Pa g_R0 = (US%L_to_Z**2 * GV%g_Earth) / (GV%Rho0) k0dt = dt*CS%kappa_0 - !### These 3 tolerances are hard-coded and fixed for now. Perhaps these could be made dynamic later? - ! tol_dksrc = 0.5*tol_ksrc_chg ; tol_dksrc_low = 1.0 - 1.0/tol_ksrc_chg ? - tol_dksrc = 10.0 ; tol_dksrc_low = 0.95 ; tol2 = 2.0*CS%kappa_tol_err + + tol_dksrc = CS%kappa_src_max_chg + if (tol_dksrc == 10.0) then + ! This is equivalent to the expression below, but avoids changes at roundoff for the default value. + tol_dksrc_low = 0.95 + else + tol_dksrc_low = (tol_dksrc - 0.5)/tol_dksrc + endif + tol2 = 2.0*CS%kappa_tol_err dt_refinements = 5 ! Selected so that 1/2^dt_refinements < 1-tol_dksrc_low use_temperature = .false. ; if (associated(tv%T)) use_temperature = .true. @@ -2062,6 +2083,12 @@ function kappa_shear_init(Time, G, GV, US, param_file, diag, CS) "components are set to 0. A reasonable value might be "//& "1e-30 m/s, which is less than an Angstrom divided by "//& "the age of the universe.", units="m s-1", default=0.0, scale=US%m_s_to_L_T) + call get_param(param_file, mdl, "KAPPA_SHEAR_MAX_KAP_SRC_CHG", CS%kappa_src_max_chg, & + "The maximum permitted increase in the kappa source within an iteration relative "//& + "to the local source; this must be greater than 1. The lower limit for the "//& + "permitted fractional decrease is (1 - 0.5/kappa_src_max_chg). These limits "//& + "could perhaps be made dynamic with an improved iterative solver.", & + default=10.0, units="nondim") call get_param(param_file, mdl, "DEBUG_KAPPA_SHEAR", CS%debug, & "If true, write debugging data for the kappa-shear code. \n"//& @@ -2072,6 +2099,11 @@ function kappa_shear_init(Time, G, GV, US, param_file, diag, CS) "If true. use an older, dimensionally inconsistent estimate of the "//& "derivative of diffusivity with energy in the Newton's method iteration. "//& "The bug causes undercorrections when dz > 1m.", default=.true.) + call get_param(param_file, mdl, "KAPPA_SHEAR_ALL_LAYER_TKE_BUG", CS%all_layer_TKE_bug, & + "If true, report back the latest estimate of TKE instead of the time average "//& + "TKE when there is mass in all layers. Otherwise always report the time "//& + "averaged TKE, as is currently done when there are some massless layers.", & + default=.true.) ! id_clock_KQ = cpu_clock_id('Ocean KS kappa_shear', grain=CLOCK_ROUTINE) ! id_clock_avg = cpu_clock_id('Ocean KS avg', grain=CLOCK_ROUTINE) ! id_clock_project = cpu_clock_id('Ocean KS project', grain=CLOCK_ROUTINE) diff --git a/src/tracer/MOM_tracer_Z_init.F90 b/src/tracer/MOM_tracer_Z_init.F90 index 02275d7ad9..b7e00b4eba 100644 --- a/src/tracer/MOM_tracer_Z_init.F90 +++ b/src/tracer/MOM_tracer_Z_init.F90 @@ -3,7 +3,6 @@ module MOM_tracer_Z_init ! This file is part of MOM6. See LICENSE.md for the license. -!use MOM_diag_to_Z, only : find_overlap, find_limited_slope use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe ! use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type @@ -73,7 +72,7 @@ function tracer_Z_init(tr, h, filename, tr_name, G, US, missing_val, land_val) logical :: has_edges, use_missing, zero_surface character(len=80) :: loc_msg - integer :: k_top, k_bot, k_bot_prev + integer :: k_top, k_bot, k_bot_prev, k_start integer :: i, j, k, kz, is, ie, js, je, nz, nz_in is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -140,7 +139,7 @@ function tracer_Z_init(tr, h, filename, tr_name, G, US, missing_val, land_val) e(nz+1) = -G%bathyT(i,j) do k=nz,1,-1 ; e(K) = e(K+1) + dilate * h(i,j,k) ; enddo - ! Create a single-column copy of tr_in. ### CHANGE THIS LATER? + ! Create a single-column copy of tr_in. Efficiency is not an issue here. do k=1,nz_in ; tr_1d(k) = tr_in(i,j,k) ; enddo k_bot = 1 ; k_bot_prev = -1 do k=1,nz @@ -149,18 +148,18 @@ function tracer_Z_init(tr, h, filename, tr_name, G, US, missing_val, land_val) elseif (e(K) < z_edges(nz_in+1)) then tr(i,j,k) = tr_1d(nz_in) else + k_start = k_bot ! The starting point for this search call find_overlap(z_edges, e(K), e(K+1), nz_in, & - k_bot, k_top, k_bot, wt, z1, z2) + k_start, k_top, k_bot, wt, z1, z2) kz = k_top if (kz /= k_bot_prev) then ! Calculate the intra-cell profile. sl_tr = 0.0 ! ; cur_tr = 0.0 - if ((kz < nz_in) .and. (kz > 1)) call & - find_limited_slope(tr_1d, z_edges, sl_tr, kz) + if ((kz < nz_in) .and. (kz > 1)) & + sl_tr = find_limited_slope(tr_1d, z_edges, kz) endif ! This is the piecewise linear form. - tr(i,j,k) = wt(kz) * & - (tr_1d(kz) + 0.5*sl_tr*(z2(kz) + z1(kz))) + tr(i,j,k) = wt(kz) * (tr_1d(kz) + 0.5*sl_tr*(z2(kz) + z1(kz))) ! For the piecewise parabolic form add the following... ! + C1_3*cur_tr*(z2(kz)**2 + z2(kz)*z1(kz) + z1(kz)**2)) do kz=k_top+1,k_bot-1 @@ -170,8 +169,8 @@ function tracer_Z_init(tr, h, filename, tr_name, G, US, missing_val, land_val) kz = k_bot ! Calculate the intra-cell profile. sl_tr = 0.0 ! ; cur_tr = 0.0 - if ((kz < nz_in) .and. (kz > 1)) call & - find_limited_slope(tr_1d, z_edges, sl_tr, kz) + if ((kz < nz_in) .and. (kz > 1)) & + sl_tr = find_limited_slope(tr_1d, z_edges, kz) ! This is the piecewise linear form. tr(i,j,k) = tr(i,j,k) + wt(kz) * & (tr_1d(kz) + 0.5*sl_tr*(z2(kz) + z1(kz))) @@ -215,7 +214,7 @@ function tracer_Z_init(tr, h, filename, tr_name, G, US, missing_val, land_val) e(nz+1) = -G%bathyT(i,j) do k=nz,1,-1 ; e(K) = e(K+1) + dilate * h(i,j,k) ; enddo - ! Create a single-column copy of tr_in. ### CHANGE THIS LATER? + ! Create a single-column copy of tr_in. Efficiency is not an issue here. do k=1,nz_in ; tr_1d(k) = tr_in(i,j,k) ; enddo k_bot = 1 do k=1,nz @@ -224,8 +223,9 @@ function tracer_Z_init(tr, h, filename, tr_name, G, US, missing_val, land_val) elseif (z_edges(nz_in) > e(K)) then tr(i,j,k) = tr_1d(nz_in) else + k_start = k_bot ! The starting point for this search call find_overlap(z_edges, e(K), e(K+1), nz_in-1, & - k_bot, k_top, k_bot, wt, z1, z2) + k_start, k_top, k_bot, wt, z1, z2) kz = k_top if (k_top < nz_in) then @@ -410,20 +410,16 @@ end subroutine read_Z_edges !! with the depth range between Z_top and Z_bot, and the fractional weights !! of each layer. It also calculates the normalized relative depths of the range !! of each layer that overlaps that depth range. - -! ### TODO: Merge with midas_vertmap.F90:find_overlap() subroutine find_overlap(e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z2) - real, dimension(:), intent(in) :: e !< Column interface heights, in arbitrary units. - real, intent(in) :: Z_top !< Top of range being mapped to, in the units of e. - real, intent(in) :: Z_bot !< Bottom of range being mapped to, in the units of e. - integer, intent(in) :: k_max !< Number of valid layers. - integer, intent(in) :: k_start !< Layer at which to start searching. - integer, intent(inout) :: k_top !< Indices of top layers that overlap with the depth - !! range. - integer, intent(inout) :: k_bot !< Indices of bottom layers that overlap with the - !! depth range. - real, dimension(:), intent(out) :: wt !< Relative weights of each layer from k_top to k_bot. - real, dimension(:), intent(out) :: z1 !< Depth of the top limits of the part of + real, dimension(:), intent(in) :: e !< Column interface heights, [Z ~> m] or other units. + real, intent(in) :: Z_top !< Top of range being mapped to, in the units of e [Z ~> m]. + real, intent(in) :: Z_bot !< Bottom of range being mapped to, in the units of e [Z ~> m]. + integer, intent(in) :: k_max !< Number of valid layers. + integer, intent(in) :: k_start !< Layer at which to start searching. + integer, intent(out) :: k_top !< Indices of top layers that overlap with the depth range. + integer, intent(out) :: k_bot !< Indices of bottom layers that overlap with the depth range. + real, dimension(:), intent(out) :: wt !< Relative weights of each layer from k_top to k_bot [nondim]. + real, dimension(:), intent(out) :: z1 !< Depth of the top limits of the part of !! a layer that contributes to a depth level, relative to the cell center and normalized !! by the cell thickness [nondim]. Note that -1/2 <= z1 < z2 <= 1/2. real, dimension(:), intent(out) :: z2 !< Depths of the bottom limit of the part of @@ -433,17 +429,19 @@ subroutine find_overlap(e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z real :: Ih, e_c, tot_wt, I_totwt integer :: k - do k=k_start,k_max ; if (e(K+1)k_max) return + if (k_top > k_max) return ! Determine the fractional weights of each layer. ! Note that by convention, e and Z_int decrease with increasing k. - if (e(K+1)<=Z_bot) then + if (e(K+1) <= Z_bot) then wt(k) = 1.0 ; k_bot = k Ih = 0.0 ; if (e(K) /= e(K+1)) Ih = 1.0 / (e(K)-e(K+1)) e_c = 0.5*(e(K)+e(K+1)) - z1(k) = (e_c - MIN(e(K),Z_top)) * Ih + z1(k) = (e_c - MIN(e(K), Z_top)) * Ih z2(k) = (e_c - Z_bot) * Ih else wt(k) = MIN(e(K),Z_top) - e(K+1) ; tot_wt = wt(k) ! These are always > 0. @@ -453,7 +451,7 @@ subroutine find_overlap(e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z z2(k) = 0.5 k_bot = k_max do k=k_top+1,k_max - if (e(K+1)<=Z_bot) then + if (e(K+1) <= Z_bot) then k_bot = k wt(k) = e(K) - Z_bot ; z1(k) = -0.5 if (e(K) /= e(K+1)) then @@ -466,38 +464,39 @@ subroutine find_overlap(e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z if (k>=k_bot) exit enddo - I_totwt = 1.0 / tot_wt + I_totwt = 0.0 ; if (tot_wt > 0.0) I_totwt = 1.0 / tot_wt do k=k_top,k_bot ; wt(k) = I_totwt*wt(k) ; enddo endif end subroutine find_overlap -!> This subroutine determines a limited slope for val to be advected with +!> This function determines a limited slope for val to be advected with !! a piecewise limited scheme. -! ### TODO: Merge with midas_vertmap.F90:find_limited_slope() -subroutine find_limited_slope(val, e, slope, k) - real, dimension(:), intent(in) :: val !< A column of values that are being interpolated. - real, dimension(:), intent(in) :: e !< Column interface heights in arbitrary units - real, intent(out) :: slope !< Normalized slope in the intracell distribution of val. - integer, intent(in) :: k !< Layer whose slope is being determined. +function find_limited_slope(val, e, k) result(slope) + real, dimension(:), intent(in) :: val !< A column of values that are being interpolated, in arbitrary units [A]. + real, dimension(:), intent(in) :: e !< A column's interface heights [Z ~> m] or other units. + integer, intent(in) :: k !< The layer whose slope is being determined. + real :: slope !< The normalized slope in the intracell distribution of + !! val [A Z-1 ~> A m-1] or other units. ! Local variables - real :: d1, d2 ! Thicknesses in the units of e. + real :: d1, d2 ! Thicknesses in the units of e [Z ~> m]. d1 = 0.5*(e(K-1)-e(K+1)) ; d2 = 0.5*(e(K)-e(K+2)) if (((val(k)-val(k-1)) * (val(k)-val(k+1)) >= 0.0) .or. (d1*d2 <= 0.0)) then - slope = 0.0 ! ; curvature = 0.0 + slope = 0.0 else - slope = (d1**2*(val(k+1) - val(k)) + d2**2*(val(k) - val(k-1))) * & + ! This line has an extra set of parentheses on the second line, so it gives slightly + ! different answers than the version of find_limited_slope in midas_vertmap.F90. + slope = ((d1**2)*(val(k+1) - val(k)) + (d2**2)*(val(k) - val(k-1))) * & ((e(K) - e(K+1)) / (d1*d2*(d1+d2))) ! slope = 0.5*(val(k+1) - val(k-1)) ! This is S.J. Lin's form of the PLM limiter. - slope = sign(1.0,slope) * min(abs(slope), & - 2.0*(max(val(k-1),val(k),val(k+1)) - val(k)), & - 2.0*(val(k) - min(val(k-1),val(k),val(k+1)))) - ! curvature = 0.0 + slope = sign(1.0, slope) * min(abs(slope), & + 2.0*(max(val(k-1), val(k), val(k+1)) - val(k)), & + 2.0*(val(k) - min(val(k-1), val(k), val(k+1)))) endif -end subroutine find_limited_slope +end function find_limited_slope end module MOM_tracer_Z_init diff --git a/src/user/Idealized_Hurricane.F90 b/src/user/Idealized_Hurricane.F90 index b4cbb32401..ff2a533d99 100644 --- a/src/user/Idealized_Hurricane.F90 +++ b/src/user/Idealized_Hurricane.F90 @@ -45,28 +45,30 @@ module Idealized_hurricane type, public :: idealized_hurricane_CS ; private ! Parameters used to compute Holland radial wind profile - real :: rho_a !< Mean air density [kg m-3] - real :: pressure_ambient !< Pressure at surface of ambient air [Pa] - real :: pressure_central !< Pressure at surface at hurricane center [Pa] - real :: rad_max_wind !< Radius of maximum winds [m] - real :: max_windspeed !< Maximum wind speeds [m s-1] - real :: hurr_translation_spd !< Hurricane translation speed [m s-1] - real :: hurr_translation_dir !< Hurricane translation speed [m s-1] - real :: gustiness !< Gustiness (optional, used in u*) [R L Z T-1 ~> Pa] + real :: rho_a !< Mean air density [R ~> kg m-3] + real :: pressure_ambient !< Pressure at surface of ambient air [R L2 T-2 ~> Pa] + real :: pressure_central !< Pressure at surface at hurricane center [R L2 T-2 ~> Pa] + real :: rad_max_wind !< Radius of maximum winds [L ~> m] + real :: max_windspeed !< Maximum wind speeds [L T-1 ~> m s-1] + real :: hurr_translation_spd !< Hurricane translation speed [L T-1 ~> m s-1] + real :: hurr_translation_dir !< Hurricane translation direction [radians] + real :: gustiness !< Gustiness (optional, used in u*) [R L Z T-2 ~> Pa] real :: Rho0 !< A reference ocean density [R ~> kg m-3] real :: Hurr_cen_Y0 !< The initial y position of the hurricane !! This experiment is conducted in a Cartesian - !! grid and this is assumed to be in meters [m] + !! grid and this is assumed to be in meters [L ~> m] real :: Hurr_cen_X0 !< The initial x position of the hurricane !! This experiment is conducted in a Cartesian - !! grid and this is assumed to be in meters [m] - real :: Holland_A !< Parameter 'A' from the Holland formula - real :: Holland_B !< Parameter 'B' from the Holland formula + !! grid and this is assumed to be in meters [L ~> m] + real :: Holland_A !< Parameter 'A' from the Holland formula [nondim] + real :: Holland_B !< Parameter 'B' from the Holland formula [nondim] real :: Holland_AxBxDP !< 'A' x 'B' x (Pressure Ambient-Pressure central) - !! for the Holland prorfile calculation + !! for the Holland prorfile calculation [R L2 T-2 ~> Pa] logical :: relative_tau !< A logical to take difference between wind - !! and surface currents to compute the stress - + !! and surface currents to compute the stress + logical :: answers_2018 !< If true, use expressions driving the idealized hurricane test + !! case that recover the answers from the end of 2018. Otherwise use + !! expressions that are rescalable and respect rotational symmetry. ! Parameters used if in SCM (single column model) mode logical :: SCM_mode !< If true this being used in Single Column Model mode @@ -74,7 +76,7 @@ module Idealized_hurricane !! provide identical wind to reproduce a previous !! experiment, where that wind formula contained !! an error) - real :: DY_from_center !< (Fixed) distance in y from storm center path [m] + real :: dy_from_center !< (Fixed) distance in y from storm center path [L ~> m] ! Par real :: PI !< Mathematical constant @@ -97,10 +99,13 @@ subroutine idealized_hurricane_wind_init(Time, G, US, param_file, CS) type(param_file_type), intent(in) :: param_file !< Input parameter structure type(idealized_hurricane_CS), pointer :: CS !< Parameter container for this module - real :: DP, C + ! Local variables + real :: dP ! The pressure difference across the hurricane [R L2 T-2 ~> Pa] + real :: C + logical :: default_2018_answers ! The default setting for the various 2018_ANSWERS flags. -! This include declares and sets the variable "version". -#include "version_variable.h" + ! This include declares and sets the variable "version". +# include "version_variable.h" if (associated(CS)) then call MOM_error(FATAL, "idealized_hurricane_wind_init called "// & @@ -118,37 +123,34 @@ subroutine idealized_hurricane_wind_init(Time, G, US, param_file, CS) ! Parameters for computing a wind profile call get_param(param_file, mdl, "IDL_HURR_RHO_AIR", CS%rho_a, & - "Air density used to compute the idealized hurricane "//& - "wind profile.", units='kg/m3', default=1.2) - call get_param(param_file, mdl, "IDL_HURR_AMBIENT_PRESSURE", & - CS%pressure_ambient, "Ambient pressure used in the "//& - "idealized hurricane wind profile.", units='Pa', & - default=101200.) - call get_param(param_file, mdl, "IDL_HURR_CENTRAL_PRESSURE", & - CS%pressure_central, "Central pressure used in the "//& - "idealized hurricane wind profile.", units='Pa', & - default=96800.) + "Air density used to compute the idealized hurricane wind profile.", & + units='kg/m3', default=1.2, scale=US%kg_m3_to_R) + call get_param(param_file, mdl, "IDL_HURR_AMBIENT_PRESSURE", CS%pressure_ambient, & + "Ambient pressure used in the idealized hurricane wind profile.", & + units='Pa', default=101200., scale=US%m_s_to_L_T**2*US%kg_m3_to_R) + call get_param(param_file, mdl, "IDL_HURR_CENTRAL_PRESSURE", CS%pressure_central, & + "Central pressure used in the idealized hurricane wind profile.", & + units='Pa', default=96800., scale=US%m_s_to_L_T**2*US%kg_m3_to_R) call get_param(param_file, mdl, "IDL_HURR_RAD_MAX_WIND", & CS%rad_max_wind, "Radius of maximum winds used in the "//& "idealized hurricane wind profile.", units='m', & - default=50.e3) + default=50.e3, scale=US%m_to_L) call get_param(param_file, mdl, "IDL_HURR_MAX_WIND", CS%max_windspeed, & "Maximum wind speed used in the idealized hurricane"// & - "wind profile.", units='m/s', default=65.) + "wind profile.", units='m/s', default=65., scale=US%m_s_to_L_T) call get_param(param_file, mdl, "IDL_HURR_TRAN_SPEED", CS%hurr_translation_spd, & "Translation speed of hurricane used in the idealized "//& - "hurricane wind profile.", units='m/s', default=5.0) + "hurricane wind profile.", units='m/s', default=5.0, scale=US%m_s_to_L_T) call get_param(param_file, mdl, "IDL_HURR_TRAN_DIR", CS%hurr_translation_dir, & "Translation direction (towards) of hurricane used in the "//& "idealized hurricane wind profile.", units='degrees', & - default=180.0) - CS%hurr_translation_dir = CS%hurr_translation_dir * CS%Deg2Rad + default=180.0, scale=CS%Deg2Rad) call get_param(param_file, mdl, "IDL_HURR_X0", CS%Hurr_cen_X0, & "Idealized Hurricane initial X position", & - units='m', default=0.) + units='m', default=0., scale=US%m_to_L) call get_param(param_file, mdl, "IDL_HURR_Y0", CS%Hurr_cen_Y0, & "Idealized Hurricane initial Y position", & - units='m', default=0.) + units='m', default=0., scale=US%m_to_L) call get_param(param_file, mdl, "IDL_HURR_TAU_CURR_REL", CS%relative_tau, & "Current relative stress switch "//& "used in the idealized hurricane wind profile.", & @@ -163,9 +165,16 @@ subroutine idealized_hurricane_wind_init(Time, G, US, param_file, CS) "Single Column mode switch "//& "used in the SCM idealized hurricane wind profile.", & units='', default=.false.) - call get_param(param_file, mdl, "IDL_HURR_SCM_LOCY", CS%DY_from_center, & + call get_param(param_file, mdl, "IDL_HURR_SCM_LOCY", CS%dy_from_center, & "Y distance of station used in the SCM idealized hurricane "//& - "wind profile.", units='m', default=50.e3) + "wind profile.", units='m', default=50.e3, scale=US%m_to_L) + call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, & + "This sets the default value for the various _2018_ANSWERS parameters.", & + default=.true.) + call get_param(param_file, mdl, "IDL_HURR_2018_ANSWERS", CS%answers_2018, & + "If true, use expressions driving the idealized hurricane test case that recover "//& + "the answers from the end of 2018. Otherwise use expressions that are rescalable "//& + "and respect rotational symmetry.", default=default_2018_answers) ! The following parameters are model run-time parameters which are used ! and logged elsewhere and so should not be logged here. The default @@ -182,13 +191,17 @@ subroutine idealized_hurricane_wind_init(Time, G, US, param_file, CS) if (CS%BR_BENCH) then - CS%rho_a = 1.2 + CS%rho_a = 1.2*US%kg_m3_to_R + endif + dP = CS%pressure_ambient - CS%pressure_central + if (CS%answers_2018) then + C = CS%max_windspeed / sqrt( US%R_to_kg_m3 * dP ) + CS%Holland_B = C**2 * US%R_to_kg_m3*CS%rho_a * exp(1.0) + else + CS%Holland_B = CS%max_windspeed**2 * CS%rho_a * exp(1.0) / dP endif - DP = CS%pressure_ambient - CS%pressure_central - C = CS%max_windspeed / sqrt( DP ) - CS%Holland_B = C**2 * CS%rho_a * exp(1.0) - CS%Holland_A = (CS%rad_max_wind)**CS%Holland_B - CS%Holland_AxBxDP = CS%Holland_A*CS%Holland_B*DP + CS%Holland_A = (US%L_to_m*CS%rad_max_wind)**CS%Holland_B + CS%Holland_AxBxDP = CS%Holland_A*CS%Holland_B*dP end subroutine idealized_hurricane_wind_init @@ -205,17 +218,16 @@ subroutine idealized_hurricane_wind_forcing(state, forces, day, G, US, CS) integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB - real :: TX,TY !< wind stress - real :: Uocn, Vocn !< Surface ocean velocity components - real :: LAT, LON !< Grid location - real :: YY, XX !< storm relative position - real :: XC, YC !< Storm center location - real :: f !< Coriolis - real :: fbench !< The benchmark 'f' value + real :: TX, TY !< wind stress components [R L Z T-2 ~> Pa] + real :: Uocn, Vocn !< Surface ocean velocity components [L T-1 ~> m s-1] + real :: YY, XX !< storm relative position [L ~> m] + real :: XC, YC !< Storm center location [L ~> m] + real :: f_local !< Local Coriolis parameter [T-1 ~> s-1] + real :: fbench !< The benchmark 'f' value [T-1 ~> s-1] real :: fbench_fac !< A factor that is set to 0 to use the - !! benchmark 'f' value + !! benchmark 'f' value [nondim] real :: rel_tau_fac !< A factor that is set to 0 to disable - !! current relative stress calculation + !! current relative stress calculation [nondim] ! Bounds for loops and memory allocation is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -233,61 +245,67 @@ subroutine idealized_hurricane_wind_forcing(state, forces, day, G, US, CS) endif !> Compute storm center location - XC = CS%Hurr_cen_X0 + (time_type_to_real(day)*CS%hurr_translation_spd*& + XC = CS%Hurr_cen_X0 + (time_type_to_real(day)*US%s_to_T * CS%hurr_translation_spd * & cos(CS%hurr_translation_dir)) - YC = CS%Hurr_cen_Y0 + (time_type_to_real(day)*CS%hurr_translation_spd*& + YC = CS%Hurr_cen_Y0 + (time_type_to_real(day)*US%s_to_T * CS%hurr_translation_spd * & sin(CS%hurr_translation_dir)) if (CS%BR_Bench) then - ! f reset to value used in generated wind for benchmark test - fbench = 5.5659e-05 - fbench_fac = 0.0 + ! f reset to value used in generated wind for benchmark test + fbench = 5.5659e-05 * US%T_to_s + fbench_fac = 0.0 else - fbench = 0.0 - fbench_fac = 1.0 + fbench = 0.0 + fbench_fac = 1.0 endif !> Computes taux do j=js,je do I=is-1,Ieq - Uocn = state%u(I,j)*REL_TAU_FAC - Vocn = 0.25*(state%v(i,J)+state%v(i+1,J-1)& - +state%v(i+1,J)+state%v(i,J-1))*REL_TAU_FAC - f = abs(0.5*US%s_to_T*(G%CoriolisBu(I,J)+G%CoriolisBu(I,J-1)))*fbench_fac + fbench + Uocn = US%m_s_to_L_T * state%u(I,j)*REL_TAU_FAC + if (CS%answers_2018) then + Vocn = US%m_s_to_L_T * 0.25*(state%v(i,J)+state%v(i+1,J-1)& + +state%v(i+1,J)+state%v(i,J-1))*REL_TAU_FAC + else + Vocn = US%m_s_to_L_T * 0.25*((state%v(i,J)+state%v(i+1,J-1)) +& + (state%v(i+1,J)+state%v(i,J-1))) * REL_TAU_FAC + endif + f_local = abs(0.5*(G%CoriolisBu(I,J)+G%CoriolisBu(I,J-1)))*fbench_fac + fbench ! Calculate position as a function of time. if (CS%SCM_mode) then YY = YC + CS%dy_from_center XX = XC else - LAT = G%geoLatCu(I,j)*1000. ! Convert Lat from km to m. - LON = G%geoLonCu(I,j)*1000. ! Convert Lon from km to m. - YY = LAT - YC - XX = LON - XC + YY = G%geoLatCu(I,j)*1000.*US%m_to_L - YC + XX = G%geoLonCu(I,j)*1000.*US%m_to_L - XC endif - call idealized_hurricane_wind_profile(CS,f,YY,XX,Uocn,Vocn,TX,TY) - forces%taux(I,j) = G%mask2dCu(I,j) * US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z * TX + call idealized_hurricane_wind_profile(CS, US, f_local, YY, XX, Uocn, Vocn, TX, TY) + forces%taux(I,j) = G%mask2dCu(I,j) * TX enddo enddo !> Computes tauy do J=js-1,Jeq do i=is,ie - Uocn = 0.25*(state%u(I,j)+state%u(I-1,j+1)& - +state%u(I-1,j)+state%u(I,j+1))*REL_TAU_FAC - Vocn = state%v(i,J)*REL_TAU_FAC - f = abs(0.5*US%s_to_T*(G%CoriolisBu(I-1,J)+G%CoriolisBu(I,J)))*fbench_fac + fbench + if (CS%answers_2018) then + Uocn = US%m_s_to_L_T * 0.25*(state%u(I,j)+state%u(I-1,j+1) + & + state%u(I-1,j)+state%u(I,j+1))*REL_TAU_FAC + else + Uocn = US%m_s_to_L_T * 0.25*((state%u(I,j)+state%u(I-1,j+1)) + & + (state%u(I-1,j)+state%u(I,j+1))) * REL_TAU_FAC + endif + Vocn = US%m_s_to_L_T * state%v(i,J)*REL_TAU_FAC + f_local = abs(0.5*(G%CoriolisBu(I-1,J)+G%CoriolisBu(I,J)))*fbench_fac + fbench ! Calculate position as a function of time. if (CS%SCM_mode) then YY = YC + CS%dy_from_center XX = XC else - LAT = G%geoLatCv(i,J)*1000. ! Convert Lat from km to m. - LON = G%geoLonCv(i,J)*1000. ! Convert Lon from km to m. - YY = LAT - YC - XX = LON - XC + YY = G%geoLatCv(i,J)*1000.*US%m_to_L - YC + XX = G%geoLonCv(i,J)*1000.*US%m_to_L - XC endif - call idealized_hurricane_wind_profile(CS, f, YY, XX, Uocn, Vocn, TX, TY) - forces%tauy(i,J) = G%mask2dCv(i,J) * US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z * TY + call idealized_hurricane_wind_profile(CS, US, f_local, YY, XX, Uocn, Vocn, TX, TY) + forces%tauy(i,J) = G%mask2dCv(i,J) * TY enddo enddo @@ -305,34 +323,34 @@ subroutine idealized_hurricane_wind_forcing(state, forces, day, G, US, CS) end subroutine idealized_hurricane_wind_forcing !> Calculate the wind speed at a location as a function of time. -subroutine idealized_hurricane_wind_profile(CS, absf, YY, XX, UOCN, VOCN, Tx, Ty) +subroutine idealized_hurricane_wind_profile(CS, US, absf, YY, XX, UOCN, VOCN, Tx, Ty) ! Author: Brandon Reichl ! Date: Nov-20-2014 ! Aug-14-2018 Generalized for non-SCM configuration ! Input parameters - type(idealized_hurricane_CS), & - pointer :: CS !< Container for SCM parameters - real, intent(in) :: absf ! s-1] + real, intent(in) :: YY !< Location in m relative to center y [L ~> m] + real, intent(in) :: XX !< Location in m relative to center x [L ~> m] + real, intent(in) :: UOCN !< X surface current [L T-1 ~> m s-1] + real, intent(in) :: VOCN !< Y surface current [L T-1 ~> m s-1] + real, intent(out) :: Tx !< X stress [R L Z T-2 ~> Pa] + real, intent(out) :: Ty !< Y stress [R L Z T-2 ~> Pa] ! Local variables ! Wind profile terms - real :: U10 - real :: radius - real :: radius10 - real :: radius_km + real :: U10 ! The 10 m wind speed [L T-1 ~> m s-1] + real :: radius ! The distance from the hurricane center [L ~> m] + real :: radius10 ! 10 times the distance from the hurricane center [L ~> m] + real :: radius_km ! The distance from the hurricane center, perhaps in km [L ~> m] or [1000 L ~> km] real :: radiusB - real :: fcor - real :: du10 - real :: du - real :: dv + real :: tmp ! A temporary variable [R L T-1 ~> kg m-2 s-1] + real :: du10 ! The magnitude of the difference between the 10 m wind and the ocean flow [L T-1 ~> m s-1] + real :: du ! The difference between the zonal 10 m wind and the zonal ocean flow [L T-1 ~> m s-1] + real :: dv ! The difference between the meridional 10 m wind and the zonal ocean flow [L T-1 ~> m s-1] real :: CD !Wind angle variables @@ -342,12 +360,12 @@ subroutine idealized_hurricane_wind_profile(CS, absf, YY, XX, UOCN, VOCN, Tx, Ty real :: A1 real :: P1 real :: Adir - real :: V_TS - real :: U_TS + real :: V_TS ! Meridional hurricane translation speed [L T-1 ~> m s-1] + real :: U_TS ! Zonal hurricane translation speed [L T-1 ~> m s-1] ! Implementing Holland (1980) parameteric wind profile - Radius = SQRT(XX**2 + YY**2) + radius = SQRT(XX**2 + YY**2) !/ BGR ! rkm - r converted to km for Holland prof. @@ -361,72 +379,91 @@ subroutine idealized_hurricane_wind_profile(CS, absf, YY, XX, UOCN, VOCN, Tx, Ty ! if not comparing to benchmark, then use correct Holland prof. radius_km = radius endif - radiusB = (radius)**CS%Holland_B + radiusB = (US%L_to_m*radius)**CS%Holland_B !/ ! Calculate U10 in the interior (inside of 10x radius of maximum wind), ! while adjusting U10 to 0 outside of 12x radius of maximum wind. - if ( (radius/CS%rad_max_wind .gt. 0.001) .and. & - (radius/CS%rad_max_wind .lt. 10.) ) then - U10 = sqrt(CS%Holland_AxBxDP*exp(-CS%Holland_A/radiusB)/(CS%rho_A*radiusB)& - +0.25*(radius_km*absf)**2) - 0.5*radius_km*absf - elseif ( (radius/CS%rad_max_wind .gt. 10.) .and. & - (radius/CS%rad_max_wind .lt. 15.) ) then - - radius10 = CS%rad_max_wind*10. + if (CS%answers_2018) then + if ( (radius > 0.001*CS%rad_max_wind) .and. (radius < 10.*CS%rad_max_wind) ) then + U10 = sqrt(CS%Holland_AxBxDP*exp(-CS%Holland_A/radiusB) / (CS%rho_a*radiusB) + & + 0.25*(radius_km*absf)**2) - 0.5*radius_km*absf + elseif ( (radius > 10.*CS%rad_max_wind) .and. (radius < 15.*CS%rad_max_wind) ) then + radius10 = CS%rad_max_wind*10. + if (CS%BR_Bench) then + radius_km = radius10/1000. + else + radius_km = radius10 + endif + radiusB = (US%L_to_m*radius10)**CS%Holland_B - if (CS%BR_Bench) then - radius_km = radius10/1000. + U10 = (sqrt(CS%Holland_AxBxDp*exp(-CS%Holland_A/radiusB) / (CS%rho_a*radiusB) + & + 0.25*(radius_km*absf)**2) - 0.5*radius_km*absf) & + * (15. - radius/CS%rad_max_wind)/5. else - radius_km = radius10 + U10 = 0. + endif + else ! This is mathematically equivalent to that is above but more accurate. + if ( (radius > 0.001*CS%rad_max_wind) .and. (radius < 10.*CS%rad_max_wind) ) then + tmp = ( 0.5*radius_km*absf) * (CS%rho_a*radiusB) + U10 = (CS%Holland_AxBxDP * exp(-CS%Holland_A/radiusB)) / & + ( tmp + sqrt(CS%Holland_AxBxDP*exp(-CS%Holland_A/radiusB) * (CS%rho_a*radiusB) + tmp**2) ) + elseif ( (radius > 10.*CS%rad_max_wind) .and. (radius < 15.*CS%rad_max_wind) ) then + radius_km = 10.0 * CS%rad_max_wind + if (CS%BR_Bench) radius_km = radius_km/1000. + radiusB = (10.0*US%L_to_m*CS%rad_max_wind)**CS%Holland_B + tmp = ( 0.5*radius_km*absf) * (CS%rho_a*radiusB) + U10 = (3.0 - radius/(5.0*CS%rad_max_wind)) * (CS%Holland_AxBxDp*exp(-CS%Holland_A/radiusB) ) / & + ( tmp + sqrt(CS%Holland_AxBxDp*exp(-CS%Holland_A/radiusB) * (CS%rho_a*radiusB) + tmp**2) ) + else + U10 = 0.0 endif - radiusB=radius10**CS%Holland_B - - U10 = (sqrt(CS%Holland_AxBxDp*exp(-CS%Holland_A/radiusB)/(CS%rho_A*radiusB)& - +0.25*(radius_km*absf)**2)-0.5*radius_km*absf) & - * (15.-radius/CS%rad_max_wind)/5. - else - U10 = 0. endif - Adir = atan2(YY,xx) + + Adir = atan2(YY,XX) + !\ ! Wind angle model following Zhang and Ulhorn (2012) ! ALPH is inflow angle positive outward. - RSTR = min(10.,radius / CS%rad_max_wind) - A0 = -0.9*RSTR - 0.09*CS%max_windspeed - 14.33 - A1 = -A0*(0.04*RSTR + 0.05*CS%Hurr_translation_spd + 0.14) - P1 = (6.88*RSTR - 9.60*CS%Hurr_translation_spd + 85.31) * CS%Deg2Rad + RSTR = min(10., radius / CS%rad_max_wind) + A0 = -0.9*RSTR - 0.09*US%L_T_to_m_s*CS%max_windspeed - 14.33 + A1 = -A0*(0.04*RSTR + 0.05*US%L_T_to_m_s*CS%hurr_translation_spd + 0.14) + P1 = (6.88*RSTR - 9.60*US%L_T_to_m_s*CS%hurr_translation_spd + 85.31) * CS%Deg2Rad ALPH = A0 - A1*cos(CS%hurr_translation_dir-Adir-P1) - if ( (radius/CS%rad_max_wind.gt.10.) .and.& - (radius/CS%rad_max_wind.lt.15.) ) then - ALPH = ALPH*(15.0-radius/CS%rad_max_wind)/5. - elseif (radius/CS%rad_max_wind.gt.15.) then + if ( (radius > 10.*CS%rad_max_wind) .and.& + (radius < 15.*CS%rad_max_wind) ) then + ALPH = ALPH*(15.0 - radius/CS%rad_max_wind)/5. + elseif (radius > 15.*CS%rad_max_wind) then ALPH = 0.0 endif ALPH = ALPH * CS%Deg2Rad ! Calculate translation speed components - U_TS = CS%hurr_translation_spd/2.*cos(CS%hurr_translation_dir) - V_TS = CS%hurr_translation_spd/2.*sin(CS%hurr_translation_dir) + U_TS = CS%hurr_translation_spd * 0.5*cos(CS%hurr_translation_dir) + V_TS = CS%hurr_translation_spd * 0.5*sin(CS%hurr_translation_dir) ! Set output (relative) winds - dU = U10*sin(Adir-CS%Pi-Alph) - UOCN + U_TS - dV = U10*cos(Adir-Alph) - VOCN + V_TS + dU = U10*sin(Adir-CS%Pi-Alph) - Uocn + U_TS + dV = U10*cos(Adir-Alph) - Vocn + V_TS ! Use a simple drag coefficient as a function of U10 (from Sullivan et al., 2010) du10 = sqrt(du**2+dv**2) - if (du10.lt.11.) then - Cd = 1.2e-3 - elseif (du10.lt.20.0) then - Cd = (0.49 + 0.065*U10)*1.e-3 + if (dU10 < 11.0*US%m_s_to_L_T) then + Cd = 1.2e-3 + elseif (dU10 < 20.0*US%m_s_to_L_T) then + if (CS%answers_2018) then + Cd = (0.49 + 0.065*US%L_T_to_m_s*U10)*1.e-3 + else + Cd = (0.49 + 0.065*US%L_T_to_m_s*dU10)*1.e-3 + endif else - Cd = 1.8e-3 + Cd = 1.8e-3 endif ! Compute stress vector - TX = CS%rho_A * Cd * sqrt(du**2 + dV**2) * dU - TY = CS%rho_A * Cd * sqrt(du**2 + dV**2) * dV + TX = US%L_to_Z * CS%rho_a * Cd * sqrt(dU**2 + dV**2) * dU + TY = US%L_to_Z * CS%rho_a * Cd * sqrt(dU**2 + dV**2) * dV end subroutine idealized_hurricane_wind_profile @@ -445,14 +482,22 @@ subroutine SCM_idealized_hurricane_wind_forcing(state, forces, day, G, US, CS) integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB real :: pie, Deg2Rad - real :: U10, A, B, C, r, f, du10, rkm ! For wind profile expression - real :: xx, t0 !for location - real :: dp, rB + real :: du10 ! The magnitude of the difference between the 10 m wind and the ocean flow [L T-1 ~> m s-1] + real :: U10 ! The 10 m wind speed [L T-1 ~> m s-1] + real :: A, B, C ! For wind profile expression + real :: rad ! The distance from the hurricane center [L ~> m] + real :: rkm ! The distance from the hurricane center, sometimes scaled to km [L ~> m] or [1000 L ~> km] + real :: f_local ! The local Coriolis parameter [T-1 ~> s-1] + real :: xx ! x-position [L ~> m] + real :: t0 !for location + real :: dP ! The pressure difference across the hurricane [R L2 T-2 ~> Pa] + real :: rB real :: Cd ! Air-sea drag coefficient - real :: Uocn, Vocn ! Surface ocean velocity components - real :: dU, dV ! Air-sea differential motion + real :: Uocn, Vocn ! Surface ocean velocity components [L T-1 ~> m s-1] + real :: dU, dV ! Air-sea differential motion [L T-1 ~> m s-1] !Wind angle variables - real :: Alph,Rstr, A0, A1, P1, Adir, transdir, V_TS, U_TS + real :: Alph,Rstr, A0, A1, P1, Adir, transdir + real :: V_TS, U_TS ! Components of the translation speed [L T-1 ~> m s-1] logical :: BR_Bench ! Bounds for loops and memory allocation is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -471,79 +516,85 @@ subroutine SCM_idealized_hurricane_wind_forcing(state, forces, day, G, US, CS) t0 = 129600. !TC 'eye' crosses (0,0) at 36 hours| transdir = pie !translation direction (-x) | !------------------------------------------------------| - dp = CS%pressure_ambient - CS%pressure_central - C = CS%max_windspeed / sqrt( DP ) - B = C**2 * CS%rho_a * exp(1.0) - if (BR_Bench) then - ! rho_a reset to value used in generated wind for benchmark test - B = C**2 * 1.2 * exp(1.0) + dP = CS%pressure_ambient - CS%pressure_central + if (CS%answers_2018) then + C = CS%max_windspeed / sqrt( US%R_to_kg_m3*dP ) + B = C**2 * US%R_to_kg_m3*CS%rho_a * exp(1.0) + if (BR_Bench) then ! rho_a reset to value used in generated wind for benchmark test + B = C**2 * 1.2 * exp(1.0) + endif + elseif (BR_Bench) then ! rho_a reset to value used in generated wind for benchmark test + B = (CS%max_windspeed**2 / dP ) * 1.2*US%kg_m3_to_R * exp(1.0) + else + B = (CS%max_windspeed**2 /dP ) * CS%rho_a * exp(1.0) endif - A = (CS%rad_max_wind/1000.)**B - f = US%s_to_T*G%CoriolisBu(is,js) ! f=f(x,y) but in the SCM is constant + + A = (US%L_to_m*CS%rad_max_wind / 1000.)**B + f_local = G%CoriolisBu(is,js) ! f=f(x,y) but in the SCM is constant if (BR_Bench) then - ! f reset to value used in generated wind for benchmark test - f = 5.5659e-05 !### A constant value in s-1. + ! f reset to value used in generated wind for benchmark test + f_local = 5.5659e-05*US%T_to_s endif !/ BR ! Calculate x position as a function of time. - xx = ( t0 - time_type_to_real(day)) * CS%hurr_translation_spd * cos(transdir) - r = sqrt(xx**2 + CS%DY_from_center**2) + xx = US%s_to_T*( t0 - time_type_to_real(day)) * CS%hurr_translation_spd * cos(transdir) + rad = sqrt(xx**2 + CS%dy_from_center**2) !/ BR - ! rkm - r converted to km for Holland prof. + ! rkm - rad converted to km for Holland prof. ! used in km due to error, correct implementation should ! not need rkm, but to match winds w/ experiment this must ! be maintained. Causes winds far from storm center to be a ! couple of m/s higher than the correct Holland prof. if (BR_Bench) then - rkm = r/1000. - rB = (rkm)**B + rkm = rad/1000. + rB = (US%L_to_m*rkm)**B else ! if not comparing to benchmark, then use correct Holland prof. - rkm = r - rB = r**B + rkm = rad + rB = (US%L_to_m*rad)**B endif !/ BR ! Calculate U10 in the interior (inside of 10x radius of maximum wind), ! while adjusting U10 to 0 outside of 12x radius of maximum wind. ! Note that rho_a is set to 1.2 following generated wind for experiment - if (r/CS%rad_max_wind > 0.001 .AND. r/CS%rad_max_wind < 10.) then - U10 = sqrt( A*B*dp*exp(-A/rB)/(1.2*rB) + 0.25*(rkm*f)**2 ) - 0.5*rkm*f - elseif (r/CS%rad_max_wind > 10. .AND. r/CS%rad_max_wind < 12.) then - r=CS%rad_max_wind*10. - if (BR_Bench) then - rkm = r/1000. - rB=rkm**B - else - rkm = r - rB = r**B - endif - U10 = ( sqrt( A*B*dp*exp(-A/rB)/(1.2*rB) + 0.25*(rkm*f)**2 ) - 0.5*rkm*f) & - * (12. - r/CS%rad_max_wind)/2. + if (rad > 0.001*CS%rad_max_wind .AND. rad < 10.*CS%rad_max_wind) then + U10 = sqrt( A*B*dP*exp(-A/rB)/(1.2*US%kg_m3_to_R*rB) + 0.25*(rkm*f_local)**2 ) - 0.5*rkm*f_local + elseif (rad > 10.*CS%rad_max_wind .AND. rad < 12.*CS%rad_max_wind) then + rad=(CS%rad_max_wind)*10. + if (BR_Bench) then + rkm = rad/1000. + rB = (US%L_to_m*rkm)**B + else + rkm = rad + rB = (US%L_to_m*rad)**B + endif + U10 = ( sqrt( A*B*dP*exp(-A/rB)/(1.2*US%kg_m3_to_R*rB) + 0.25*(rkm*f_local)**2 ) - 0.5*rkm*f_local) & + * (12. - rad/CS%rad_max_wind)/2. else - U10 = 0. + U10 = 0. endif - Adir = atan2(CS%DY_from_center,xx) + Adir = atan2(CS%dy_from_center,xx) !/ BR ! Wind angle model following Zhang and Ulhorn (2012) ! ALPH is inflow angle positive outward. - RSTR = min(10.,r / CS%rad_max_wind) - A0 = -0.9*RSTR -0.09*CS%max_windspeed - 14.33 - A1 = -A0 *(0.04*RSTR +0.05*CS%hurr_translation_spd+0.14) - P1 = (6.88*RSTR -9.60*CS%hurr_translation_spd+85.31)*pie/180. + RSTR = min(10., rad / CS%rad_max_wind) + A0 = -0.9*RSTR - 0.09*US%L_T_to_m_s*CS%max_windspeed - 14.33 + A1 = -A0 *(0.04*RSTR + 0.05*US%L_T_to_m_s*CS%hurr_translation_spd + 0.14) + P1 = (6.88*RSTR - 9.60*US%L_T_to_m_s*CS%hurr_translation_spd + 85.31)*pie/180. ALPH = A0 - A1*cos( (TRANSDIR - ADIR ) - P1) - if (r/CS%rad_max_wind > 10. .AND. r/CS%rad_max_wind < 12.) then - ALPH = ALPH* (12. - r/CS%rad_max_wind)/2. - elseif (r/CS%rad_max_wind > 12.) then - ALPH = 0.0 + if (rad > 10.*CS%rad_max_wind .AND. rad < 12.*CS%rad_max_wind) then + ALPH = ALPH* (12. - rad/CS%rad_max_wind)/2. + elseif (rad > 12.*CS%rad_max_wind) then + ALPH = 0.0 endif ALPH = ALPH * Deg2Rad !/BR ! Prepare for wind calculation ! X_TS is component of translation speed added to wind vector ! due to background steering wind. - U_TS = CS%hurr_translation_spd/2.*cos(transdir) - V_TS = CS%hurr_translation_spd/2.*sin(transdir) + U_TS = CS%hurr_translation_spd*0.5*cos(transdir) + V_TS = CS%hurr_translation_spd*0.5*sin(transdir) ! Set the surface wind stresses, in [Pa]. A positive taux ! accelerates the ocean to the (pseudo-)east. @@ -553,9 +604,9 @@ subroutine SCM_idealized_hurricane_wind_forcing(state, forces, day, G, US, CS) !/BR ! Turn off surface current for stress calculation to be ! consistent with test case. - Uocn = 0.!state%u(I,j) - Vocn = 0.!0.25*( (state%v(i,J) + state%v(i+1,J-1)) & - ! +(state%v(i+1,J) + state%v(i,J-1)) ) + Uocn = 0. ! state%u(I,j) + Vocn = 0. ! 0.25*( (state%v(i,J) + state%v(i+1,J-1)) & + ! +(state%v(i+1,J) + state%v(i,J-1)) ) !/BR ! Wind vector calculated from location/direction (sin/cos flipped b/c ! cyclonic wind is 90 deg. phase shifted from position angle). @@ -565,35 +616,41 @@ subroutine SCM_idealized_hurricane_wind_forcing(state, forces, day, G, US, CS) !BR ! Add a simple drag coefficient as a function of U10 | !/----------------------------------------------------| - du10=sqrt(du**2+dv**2) - if (du10 < 11.) then - Cd = 1.2e-3 - elseif (du10 < 20.) then - Cd = (0.49 + 0.065 * U10 )*0.001 + du10 = sqrt(du**2+dv**2) + if (dU10 < 11.0*US%m_s_to_L_T) then + Cd = 1.2e-3 + elseif (dU10 < 20.0*US%m_s_to_L_T) then + if (CS%answers_2018) then + Cd = (0.49 + 0.065 * US%L_T_to_m_s*U10 )*0.001 + else + Cd = (0.49 + 0.065 * US%L_T_to_m_s*dU10 )*0.001 + endif else - Cd = 0.0018 + Cd = 0.0018 endif - forces%taux(I,j) = CS%rho_a * US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z * & - G%mask2dCu(I,j) * Cd*sqrt(du**2+dV**2)*dU + forces%taux(I,j) = CS%rho_a * US%L_to_Z * G%mask2dCu(I,j) * Cd*du10*dU enddo ; enddo !/BR ! See notes above do J=js-1,Jeq ; do i=is,ie - Uocn = 0.!0.25*( (state%u(I,j) + state%u(I-1,j+1)) & - ! +(state%u(I-1,j) + state%u(I,j+1)) ) - Vocn = 0.!state%v(i,J) + Uocn = 0. ! 0.25*( (state%u(I,j) + state%u(I-1,j+1)) & + ! +(state%u(I-1,j) + state%u(I,j+1)) ) + Vocn = 0. ! state%v(i,J) dU = U10*sin(Adir-pie-Alph) - Uocn + U_TS dV = U10*cos(Adir-Alph) - Vocn + V_TS du10=sqrt(du**2+dv**2) - if (du10 < 11.) then - Cd = 1.2e-3 - elseif (du10 < 20.) then - Cd = (0.49 + 0.065 * U10 )*0.001 + if (dU10 < 11.0*US%m_s_to_L_T) then + Cd = 1.2e-3 + elseif (dU10 < 20.0*US%m_s_to_L_T) then + if (CS%answers_2018) then + Cd = (0.49 + 0.065 * US%L_T_to_m_s*U10 )*0.001 + else + Cd = (0.49 + 0.065 * US%L_T_to_m_s*dU10 )*0.001 + endif else - Cd = 0.0018 + Cd = 0.0018 endif - forces%tauy(I,j) = CS%rho_a * US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z * & - G%mask2dCv(I,j) * Cd*du10*dV + forces%tauy(I,j) = CS%rho_a * US%L_to_Z * G%mask2dCv(I,j) * Cd*dU10*dV enddo ; enddo ! Set the surface friction velocity [m s-1]. ustar is always positive. do j=js,je ; do i=is,ie diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index a048d85491..e7361bf13c 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -653,7 +653,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar) ! uniform cases. ! call DHH85_mid(GV, US, Midpoint, UStokes) ! Putting into x-direction, so setting y direction to 0 - CS%US_y(ii,JJ,kk) = 0.0 !### Note that =0 should be =US - RWH + CS%US_y(ii,JJ,kk) = 0.0 + ! For rotational symmetry there should be the option for this to become = UStokes ! bgr - see note above, but this is true ! if this is used for anything ! other than simple LES comparison @@ -1081,7 +1082,7 @@ end subroutine Get_StokesSL_LiFoxKemper subroutine Get_SL_Average_Prof( GV, AvgDepth, H, Profile, Average ) type(verticalGrid_type), & intent(in) :: GV !< Ocean vertical grid structure - real, intent(in) :: AvgDepth !< Depth to average over [Z ~> m]. + real, intent(in) :: AvgDepth !< Depth to average over (negative) [Z ~> m]. real, dimension(SZK_(GV)), & intent(in) :: H !< Grid thickness [H ~> m or kg m-2] real, dimension(SZK_(GV)), & @@ -1090,7 +1091,7 @@ subroutine Get_SL_Average_Prof( GV, AvgDepth, H, Profile, Average ) real, intent(out) :: Average !< Output quantity averaged over depth AvgDepth [arbitrary] !! (used here for Stokes drift) !Local variables - real :: top, midpoint, bottom ! Depths [Z ~> m]. + real :: top, midpoint, bottom ! Depths, negative downward [Z ~> m]. real :: Sum integer :: kk @@ -1103,17 +1104,25 @@ subroutine Get_SL_Average_Prof( GV, AvgDepth, H, Profile, Average ) Top = Bottom MidPoint = Bottom - GV%H_to_Z * 0.5*h(kk) Bottom = Bottom - GV%H_to_Z * h(kk) - if (AvgDepth < Bottom) then !Whole cell within H_LA + if (AvgDepth < Bottom) then ! The whole cell is within H_LA Sum = Sum + Profile(kk) * (GV%H_to_Z * H(kk)) - elseif (AvgDepth < Top) then !partial cell within H_LA + elseif (AvgDepth < Top) then ! A partial cell is within H_LA Sum = Sum + Profile(kk) * (Top-AvgDepth) + exit + else + exit endif enddo - ! Divide by AvgDepth !### Consider dividing by the depth in the column if that is smaller. -RWH - Average = Sum / abs(AvgDepth) + ! Divide by AvgDepth or the depth in the column, whichever is smaller. + if (abs(AvgDepth) <= abs(Bottom)) then + Average = Sum / abs(AvgDepth) + elseif (abs(Bottom) > 0.0) then + Average = Sum / abs(Bottom) + else + Average = 0.0 + endif - return end subroutine Get_SL_Average_Prof !> Get SL averaged Stokes drift from the banded Spectrum method @@ -1153,7 +1162,7 @@ end subroutine Get_SL_Average_Band subroutine DHH85_mid(GV, US, zpt, UStokes) type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, intent(in) :: ZPT !< Depth to get Stokes drift [Z ~> m]. !### THIS IS NOT USED YET. + real, intent(in) :: zpt !< Depth to get Stokes drift [Z ~> m]. real, intent(out) :: UStokes !< Stokes drift [m s-1] ! real :: ann, Bnn, Snn, Cnn, Dnn @@ -1195,7 +1204,7 @@ subroutine DHH85_mid(GV, US, zpt, UStokes) exp(-bnn*(omega_peak/omega)**4)*Cnn**Dnn ! Stokes units m (multiply by frequency range for units of m/s) Stokes = 2.0 * wavespec * omega**3 * & - exp( 2.0 * omega**2 * zpt / GV%mks_g_Earth) / GV%mks_g_Earth + exp( 2.0 * omega**2 * US%Z_to_m*zpt / GV%mks_g_Earth) / GV%mks_g_Earth UStokes = UStokes + Stokes*domega omega = omega + domega enddo