From b758406e33f1a7ef4959e72ef1896e37cd7a2dcd Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 22 Oct 2021 16:20:10 -0400 Subject: [PATCH] +Argument cleanup in vertical diffusivity code Cleaned up 26 falsely optional or unused arguments in the vertical diffusivity code, and related changes. Several descriptive comments were also corrected, including the correction of the units of 10 variables related to CVMix_KPP. This commit includes: - Made the Kd_int arguments to set_diffusivity() and 3 subsidiary routines mandatory and reordered the arguments so that the non-optional arguments come before the grid types - Made the halo_TS and double_diffuse arguments to set_diffusivity_init() mandatory. - Made the Time argument to ALE_sponge() mandatory. - Made the Kd and Kv arguments to calculate_CVMIX_conv() mandatory. - Removed the unused halo argument to adjust_salt(). - Removed the unused Kddt_convect argument to full_convection(). - Made the halo arguments to full_convection()and smoothed_dRdT_dRdS() mandatory. - Made the useALEalgorithm argument to geothermal_init() mandatory. - Removed the unused initialize_all arguments to Calculate_kappa_shear() and Calc_kappa_shear_vertex(). - Removed the unused I_Ld2_1d and dz_Int_1d arguments to kappa_shear_column(). - Made 3 arguments to calculate_projected_state() mandatory and reordered the arguments accordingly. - Eliminating the unused skip_diags arguments to calculateBuoyancyFlux() and extractFluxes(), which are now effectively always false. All answers are bitwise identical, and no output is changed. --- src/core/MOM_forcing_type.F90 | 43 ++-- .../vertical/MOM_ALE_sponge.F90 | 7 +- .../vertical/MOM_CVMix_KPP.F90 | 8 +- .../vertical/MOM_CVMix_conv.F90 | 28 +-- .../vertical/MOM_diabatic_aux.F90 | 8 +- .../vertical/MOM_diabatic_driver.F90 | 71 +++--- .../vertical/MOM_full_convection.F90 | 59 ++--- .../vertical/MOM_geothermal.F90 | 6 +- .../vertical/MOM_kappa_shear.F90 | 181 ++++++--------- .../vertical/MOM_set_diffusivity.F90 | 217 +++++++----------- 10 files changed, 248 insertions(+), 380 deletions(-) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 2429ce9d2d..a67d440dfe 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -390,7 +390,7 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & FluxRescaleDepth, useRiverHeatContent, useCalvingHeatContent, & h, T, netMassInOut, netMassOut, net_heat, net_salt, pen_SW_bnd, tv, & aggregate_FW, nonpenSW, netmassInOut_rate, net_Heat_Rate, & - net_salt_rate, pen_sw_bnd_Rate, skip_diags) + net_salt_rate, pen_sw_bnd_Rate) type(ocean_grid_type), intent(in) :: G !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure @@ -452,7 +452,6 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & real, dimension(max(1,nsw),G%isd:G%ied), & optional, intent(out) :: pen_sw_bnd_rate !< Rate of penetrative shortwave heating !! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]. - logical, optional, intent(in) :: skip_diags !< If present and true, skip calculating diagnostics ! local real :: htot(SZI_(G)) ! total ocean depth [H ~> m or kg m-2] @@ -492,7 +491,6 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & is = G%isc ; ie = G%iec ; nz = GV%ke calculate_diags = .true. - if (present(skip_diags)) calculate_diags = .not. skip_diags ! error checking @@ -914,7 +912,7 @@ end subroutine extractFluxes2d !! extractFluxes routine allows us to get "stuf per time" rather than the time integrated !! fluxes needed in other routines that call extractFluxes. subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt, tv, j, & - buoyancyFlux, netHeatMinusSW, netSalt, skip_diags) + buoyancyFlux, netHeatMinusSW, netSalt) type(ocean_grid_type), intent(in) :: G !< ocean grid type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -928,12 +926,11 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt type(thermo_var_ptrs), intent(inout) :: tv !< thermodynamics type integer, intent(in) :: j !< j-row to work on real, dimension(SZI_(G),SZK_(GV)+1), intent(inout) :: buoyancyFlux !< buoyancy fluxes [L2 T-3 ~> m2 s-3] - real, dimension(SZI_(G)), intent(inout) :: netHeatMinusSW !< surf Heat flux + real, dimension(SZI_(G)), intent(inout) :: netHeatMinusSW !< Surface heat flux excluding shortwave !! [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1] - real, dimension(SZI_(G)), intent(inout) :: netSalt !< surf salt flux + real, dimension(SZI_(G)), intent(inout) :: netSalt !< surface salt flux !! [ppt H s-1 ~> ppt m s-1 or ppt kg m-2 s-1] - logical, optional, intent(in) :: skip_diags !< If present and true, skip calculating - !! diagnostics inside extractFluxes1d() + ! local variables integer :: k real, parameter :: dt = 1. ! to return a rate from extractFluxes1d @@ -942,12 +939,12 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt ! [H s-1 ~> m s-1 or kg m-2 s-1] real, dimension(SZI_(G)) :: netHeat ! net temp flux [degC H s-1 ~> degC m s-2 or degC kg m-2 s-1] real, dimension(max(nsw,1), SZI_(G)) :: penSWbnd ! penetrating SW radiation by band - ! [degC H ~> degC m or degC kg m-2] + ! [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1] real, dimension(SZI_(G)) :: pressure ! pressure at the surface [R L2 T-2 ~> Pa] real, dimension(SZI_(G)) :: dRhodT ! density partial derivative wrt temp [R degC-1 ~> kg m-3 degC-1] real, dimension(SZI_(G)) :: dRhodS ! density partial derivative wrt saln [R ppt-1 ~> kg m-3 ppt-1] real, dimension(SZI_(G),SZK_(GV)+1) :: netPen ! The net penetrating shortwave radiation at each level - ! [degC H ~> degC m or degC kg m-2] + ! [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1] logical :: useRiverHeatContent logical :: useCalvingHeatContent @@ -978,7 +975,7 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt*US%s_to_T, & depthBeforeScalingFluxes, useRiverHeatContent, useCalvingHeatContent, & h(:,j,:), Temp(:,j,:), netH, netEvap, netHeatMinusSW, & - netSalt, penSWbnd, tv, .false., skip_diags=skip_diags) + netSalt, penSWbnd, tv, .false.) ! Sum over bands and attenuate as a function of depth ! netPen is the netSW as a function of depth @@ -1011,7 +1008,7 @@ end subroutine calculateBuoyancyFlux1d !> Calculates surface buoyancy flux by adding up the heat, FW and salt fluxes, !! for 2d arrays. This is a wrapper for calculateBuoyancyFlux1d. subroutine calculateBuoyancyFlux2d(G, GV, US, fluxes, optics, h, Temp, Salt, tv, & - buoyancyFlux, netHeatMinusSW, netSalt, skip_diags) + buoyancyFlux, netHeatMinusSW, netSalt) type(ocean_grid_type), intent(in) :: G !< ocean grid type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -1021,26 +1018,18 @@ subroutine calculateBuoyancyFlux2d(G, GV, US, fluxes, optics, h, Temp, Salt, tv, real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: Temp !< temperature [degC] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: Salt !< salinity [ppt] type(thermo_var_ptrs), intent(inout) :: tv !< thermodynamics type - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: buoyancyFlux !< buoyancy fluxes [L2 T-3 ~> m2 s-3] - real, dimension(SZI_(G),SZJ_(G)), optional, intent(inout) :: netHeatMinusSW !< surf temp flux - !! [degC H ~> degC m or degC kg m-2] - real, dimension(SZI_(G),SZJ_(G)), optional, intent(inout) :: netSalt !< surf salt flux - !! [ppt H ~> ppt m or ppt kg m-2] - logical, optional, intent(in) :: skip_diags !< If present and true, skip calculating - !! diagnostics inside extractFluxes1d() + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: buoyancyFlux !< buoyancy fluxes [L2 T-3 ~> m2 s-3] + real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: netHeatMinusSW !< surface heat flux excluding shortwave + !! [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1] + real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: netSalt !< Net surface salt flux + !! [ppt H s-1 ~> ppt m s-1 or ppt kg m-2 s-1] ! local variables - real, dimension( SZI_(G) ) :: netT ! net temperature flux [degC H s-1 ~> degC m s-2 or degC kg m-2 s-1] - real, dimension( SZI_(G) ) :: netS ! net saln flux !! [ppt H s-1 ~> ppt m s-1 or ppt kg m-2 s-1] integer :: j - netT(G%isc:G%iec) = 0. ; netS(G%isc:G%iec) = 0. - - !$OMP parallel do default(shared) firstprivate(netT,netS) + !$OMP parallel do default(shared) do j=G%jsc,G%jec call calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, optics_nbands(optics), h, Temp, Salt, & - tv, j, buoyancyFlux(:,j,:), netT, netS, skip_diags=skip_diags) - if (present(netHeatMinusSW)) netHeatMinusSW(G%isc:G%iec,j) = netT(G%isc:G%iec) - if (present(netSalt)) netSalt(G%isc:G%iec,j) = netS(G%isc:G%iec) + tv, j, buoyancyFlux(:,j,:), netHeatMinusSW(:,j), netSalt(:,j)) enddo end subroutine calculateBuoyancyFlux2d diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 1225487eaf..4d179e2bfb 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -840,7 +840,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) real, intent(in) :: dt !< The amount of time covered by this call [T ~> s]. type(ALE_sponge_CS), pointer :: CS !< A pointer to the control structure for this module !! that is set by a previous call to initialize_ALE_sponge (in). - type(time_type), optional, intent(in) :: Time !< The current model date + type(time_type), intent(in) :: Time !< The current model date real :: damp ! The timestep times the local damping coefficient [nondim]. real :: I1pdamp ! I1pdamp is 1/(1 + damp). [nondim]. @@ -885,8 +885,6 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) endif if (CS%time_varying_sponges) then - if (.not. present(Time)) & - call MOM_error(FATAL,"apply_ALE_sponge: No time information provided") do m=1,CS%fldno nz_data = CS%Ref_val(m)%nz_data call horiz_interp_and_extrap_tracer(CS%Ref_val(m)%id, Time, 1.0, G, sp_val, mask_z, z_in, & @@ -971,9 +969,6 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) if (CS%sponge_uv) then if (CS%time_varying_sponges) then - if (.not. present(Time)) & - call MOM_error(FATAL,"apply_ALE_sponge: No time information provided") - nz_data = CS%Ref_val_u%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, & diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index 53dddcf168..b9ceb85cc5 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -1389,8 +1389,8 @@ subroutine KPP_NonLocalTransport_temp(CS, G, GV, h, nonLocalTrans, surfFlux, & type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer/level thickness [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: nonLocalTrans !< Non-local transport [nondim] - real, dimension(SZI_(G),SZJ_(G)), intent(in) :: surfFlux !< Surface flux of scalar - !! [conc H s-1 ~> conc m s-1 or conc kg m-2 s-1] + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: surfFlux !< Surface flux of temperature + !! [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1] real, intent(in) :: dt !< Time-step [s] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: scalar !< temperature real, intent(in) :: C_p !< Seawater specific heat capacity [J kg-1 degC-1] @@ -1451,8 +1451,8 @@ subroutine KPP_NonLocalTransport_saln(CS, G, GV, h, nonLocalTrans, surfFlux, dt, type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer/level thickness [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: nonLocalTrans !< Non-local transport [nondim] - real, dimension(SZI_(G),SZJ_(G)), intent(in) :: surfFlux !< Surface flux of scalar - !! [conc H s-1 ~> conc m s-1 or conc kg m-2 s-1] + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: surfFlux !< Surface flux of salt + !! [ppt H s-1 ~> ppt m s-1 or ppt kg m-2 s-1] real, intent(in) :: dt !< Time-step [s] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: scalar !< Scalar (scalar units [conc]) diff --git a/src/parameterizations/vertical/MOM_CVMix_conv.F90 b/src/parameterizations/vertical/MOM_CVMix_conv.F90 index a615c9f40b..138e932c22 100644 --- a/src/parameterizations/vertical/MOM_CVMix_conv.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_conv.F90 @@ -150,10 +150,10 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl, Kd, Kv, Kd_aux) !! by a previous call to CVMix_conv_init. real, dimension(SZI_(G),SZJ_(G)), intent(in) :: hbl !< Depth of ocean boundary layer [Z ~> m] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & - optional, intent(inout) :: Kd !< Diapycnal diffusivity at each interface that + intent(inout) :: Kd !< Diapycnal diffusivity at each interface that !! will be incremented here [Z2 T-1 ~> m2 s-1]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & - optional, intent(inout) :: KV !< Viscosity at each interface that will be + intent(inout) :: KV !< Viscosity at each interface that will be !! incremented here [Z2 T-1 ~> m2 s-1]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & optional, intent(inout) :: Kd_aux !< A second diapycnal diffusivity at each @@ -243,12 +243,10 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl, Kd, Kv, Kd_aux) max_nlev=GV%ke, & OBL_ind=kOBL) - if (present(Kd)) then - ! Increment the diffusivity outside of the boundary layer. - do K=max(1,kOBL+1),GV%ke+1 - Kd(i,j,K) = Kd(i,j,K) + US%m2_s_to_Z2_T * kd_col(K) - enddo - endif + ! Increment the diffusivity outside of the boundary layer. + do K=max(1,kOBL+1),GV%ke+1 + Kd(i,j,K) = Kd(i,j,K) + US%m2_s_to_Z2_T * kd_col(K) + enddo if (present(Kd_aux)) then ! Increment the other diffusivity outside of the boundary layer. do K=max(1,kOBL+1),GV%ke+1 @@ -256,12 +254,10 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl, Kd, Kv, Kd_aux) enddo endif - if (present(Kv)) then - ! Increment the viscosity outside of the boundary layer. - do K=max(1,kOBL+1),GV%ke+1 - Kv(i,j,K) = Kv(i,j,K) + US%m2_s_to_Z2_T * kv_col(K) - enddo - endif + ! Increment the viscosity outside of the boundary layer. + do K=max(1,kOBL+1),GV%ke+1 + Kv(i,j,K) = Kv(i,j,K) + US%m2_s_to_Z2_T * kv_col(K) + enddo ! Store 3-d arrays for diagnostics. if (CS%id_kv_conv > 0) then @@ -288,8 +284,8 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl, Kd, Kv, Kd_aux) ! call hchksum(Kd_conv, "MOM_CVMix_conv: Kd_conv", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s) ! if (CS%id_kv_conv > 0) & ! call hchksum(Kv_conv, "MOM_CVMix_conv: Kv_conv", G%HI, haloshift=0, scale=US%m2_s_to_Z2_T) - if (present(Kd)) call hchksum(Kd, "MOM_CVMix_conv: Kd", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s) - if (present(Kv)) call hchksum(Kv, "MOM_CVMix_conv: Kv", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s) + call hchksum(Kd, "MOM_CVMix_conv: Kd", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s) + call hchksum(Kv, "MOM_CVMix_conv: Kv", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s) endif ! send diagnostics to post_data diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 8d53594ebb..4c822309d0 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -318,7 +318,7 @@ end subroutine differential_diffuse_T_S !> This subroutine keeps salinity from falling below a small but positive threshold. !! This usually occurs when the ice model attempts to extract more salt then !! is actually available to it from the ocean. -subroutine adjust_salt(h, tv, G, GV, CS, halo) +subroutine adjust_salt(h, tv, G, GV, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & @@ -327,7 +327,6 @@ subroutine adjust_salt(h, tv, G, GV, CS, halo) !! available thermodynamic fields. type(diabatic_aux_CS), intent(in) :: CS !< The control structure returned by a previous !! call to diabatic_aux_init. - integer, optional, intent(in) :: halo !< Halo width over which to work ! local variables real :: salt_add_col(SZI_(G),SZJ_(G)) !< The accumulated salt requirement [ppt R Z ~> gSalt m-2] @@ -336,9 +335,6 @@ subroutine adjust_salt(h, tv, G, GV, CS, halo) integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - if (present(halo)) then - is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo - endif ! call cpu_clock_begin(id_clock_adjust_salt) @@ -1024,7 +1020,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t optional, intent(out) :: dSV_dS !< Partial derivative of specific volume with !! salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1]. real, dimension(SZI_(G),SZJ_(G)), & - optional, intent(out) :: SkinBuoyFlux !< Buoyancy flux at surface [Z2 T-3 ~> m2 s-3]. + optional, intent(out) :: SkinBuoyFlux !< Buoyancy flux at surface [Z2 T-3 ~> m2 s-3]. ! Local variables integer, parameter :: maxGroundings = 5 diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index ff8c270a89..d3f92e99cc 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -87,9 +87,9 @@ module MOM_diabatic_driver ! vary with the Boussinesq approximation, the Boussinesq variant is given first. !> Control structure for this module -type, public:: diabatic_CS; private +type, public :: diabatic_CS ; private - logical :: use_legacy_diabatic !< If true (default), use the a legacy version of the diabatic + logical :: use_legacy_diabatic !< If true (default), use a legacy version of the diabatic !! algorithm. This is temporary and is needed to avoid change !! in answers. logical :: bulkmixedlayer !< If true, a refined bulk mixed layer is used with @@ -242,11 +242,14 @@ module MOM_diabatic_driver type(group_pass_type) :: pass_Kv !< For group halo pass type(diag_grid_storage) :: diag_grids_prev!< Stores diagnostic grids at some previous point in the algorithm ! Data arrays for communicating between components - real, allocatable, dimension(:,:,:) :: KPP_NLTheat !< KPP non-local transport for heat [m s-1] - real, allocatable, dimension(:,:,:) :: KPP_NLTscalar !< KPP non-local transport for scalars [m s-1] + !### Why are these arrays in this control structure, and not local variables in the various routines? + real, allocatable, dimension(:,:,:) :: KPP_NLTheat !< KPP non-local transport for heat [nondim] + real, allocatable, dimension(:,:,:) :: KPP_NLTscalar !< KPP non-local transport for scalars [nondim] real, allocatable, dimension(:,:,:) :: KPP_buoy_flux !< KPP forcing buoyancy flux [L2 T-3 ~> m2 s-3] - real, allocatable, dimension(:,:) :: KPP_temp_flux !< KPP effective temperature flux [degC m s-1] - real, allocatable, dimension(:,:) :: KPP_salt_flux !< KPP effective salt flux [ppt m s-1] + real, allocatable, dimension(:,:) :: KPP_temp_flux !< KPP effective temperature flux + !! [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1] + real, allocatable, dimension(:,:) :: KPP_salt_flux !< KPP effective salt flux + !! [ppt H s-1 ~> ppt m s-1 or ppt kg m-2 s-1] type(time_type), pointer :: Time !< Pointer to model time (needed for sponges) end type diabatic_CS @@ -274,8 +277,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m] type(forcing), intent(inout) :: fluxes !< points to forcing fields !! unused fields have NULL ptrs - type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, BBL properies, and - type(accel_diag_ptrs), intent(inout) :: ADp !< related points to accelerations in momentum + type(vertvisc_type), intent(inout) :: visc !< Structure with vertical viscosities, + !! BBL properties and related fields + type(accel_diag_ptrs), intent(inout) :: ADp !< Points to accelerations in momentum !! equations, to enable the later derived !! diagnostics, like energy budgets type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations @@ -465,8 +469,9 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m] type(forcing), intent(inout) :: fluxes !< points to forcing fields !! unused fields have NULL ptrs - type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, BBL properies, and - type(accel_diag_ptrs), intent(inout) :: ADp !< related points to accelerations in momentum + type(vertvisc_type), intent(inout) :: visc !< Structure with vertical viscosities, + !! BBL properties and related fields + type(accel_diag_ptrs), intent(inout) :: ADp !< Points to accelerations in momentum !! equations, to enable the later derived !! diagnostics, like energy budgets type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations @@ -577,11 +582,11 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim if (CS%debug) & call MOM_state_chksum("before set_diffusivity", u, v, h, G, GV, US, haloshift=CS%halo_TS_diff) if (CS%double_diffuse) then - call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, G, GV, US, CS%set_diff_CSp, & - Kd_int=Kd_int, Kd_extra_T=Kd_extra_T, Kd_extra_S=Kd_extra_S) + call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, Kd_int, G, GV, US, & + CS%set_diff_CSp, Kd_extra_T=Kd_extra_T, Kd_extra_S=Kd_extra_S) else - call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, G, GV, US, & - CS%set_diff_CSp, Kd_int=Kd_int) + call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, Kd_int, G, GV, US, & + CS%set_diff_CSp) endif call cpu_clock_end(id_clock_set_diffusivity) if (showCallTree) call callTree_waypoint("done with set_diffusivity (diabatic)") @@ -718,7 +723,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim ! Calculate vertical mixing due to convection (computed via CVMix) if (CS%use_CVMix_conv) then ! Increment vertical diffusion and viscosity due to convection - call calculate_CVMix_conv(h, tv, G, GV, US, CS%CVMix_conv_CSp, Hml, Kd=Kd_int, Kv=visc%Kv_slow) + call calculate_CVMix_conv(h, tv, G, GV, US, CS%CVMix_conv_CSp, Hml, Kd_int, visc%Kv_slow) endif ! This block sets ent_t and ent_s from h and Kd_int. @@ -1049,8 +1054,9 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m] type(forcing), intent(inout) :: fluxes !< points to forcing fields !! unused fields have NULL ptrs - type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, BBL properies, and - type(accel_diag_ptrs), intent(inout) :: ADp !< related points to accelerations in momentum + type(vertvisc_type), intent(inout) :: visc !< Structure with vertical viscosities, + !! BBL properties and related fields + type(accel_diag_ptrs), intent(inout) :: ADp !< Points to accelerations in momentum !! equations, to enable the later derived !! diagnostics, like energy budgets type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations @@ -1161,11 +1167,11 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, if (CS%debug) & call MOM_state_chksum("before set_diffusivity", u, v, h, G, GV, US, haloshift=CS%halo_TS_diff) if (CS%double_diffuse) then - call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, G, GV, US, CS%set_diff_CSp, & - Kd_int=Kd_heat, Kd_extra_T=Kd_extra_T, Kd_extra_S=Kd_extra_S) + call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, Kd_heat, G, GV, US, & + CS%set_diff_CSp, Kd_extra_T=Kd_extra_T, Kd_extra_S=Kd_extra_S) else - call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, G, GV, US, & - CS%set_diff_CSp, Kd_int=Kd_heat) + call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, Kd_heat, G, GV, US, & + CS%set_diff_CSp) endif call cpu_clock_end(id_clock_set_diffusivity) if (showCallTree) call callTree_waypoint("done with set_diffusivity (diabatic)") @@ -1270,9 +1276,9 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, if (CS%use_CVMix_conv) then ! Increment vertical diffusion and viscosity due to convection if (CS%useKPP) then - call calculate_CVMix_conv(h, tv, G, GV, US, CS%CVMix_conv_CSp, Hml, Kd=Kd_heat, Kv=visc%Kv_shear, Kd_aux=Kd_salt) + call calculate_CVMix_conv(h, tv, G, GV, US, CS%CVMix_conv_CSp, Hml, Kd_heat, visc%Kv_shear, Kd_aux=Kd_salt) else - call calculate_CVMix_conv(h, tv, G, GV, US, CS%CVMix_conv_CSp, Hml, Kd=Kd_heat, Kv=visc%Kv_slow, Kd_aux=Kd_salt) + call calculate_CVMix_conv(h, tv, G, GV, US, CS%CVMix_conv_CSp, Hml, Kd_heat, visc%Kv_slow, Kd_aux=Kd_salt) endif endif @@ -1558,8 +1564,9 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m] type(forcing), intent(inout) :: fluxes !< points to forcing fields !! unused fields have NULL ptrs - type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, BBL properies, and - type(accel_diag_ptrs), intent(inout) :: ADp !< related points to accelerations in momentum + type(vertvisc_type), intent(inout) :: visc !< Structure with vertical viscosities, + !! BBL properties and related fields + type(accel_diag_ptrs), intent(inout) :: ADp !< Points to accelerations in momentum !! equations, to enable the later derived !! diagnostics, like energy budgets type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations @@ -1764,11 +1771,11 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e if (CS%debug) & call MOM_state_chksum("before set_diffusivity", u, v, h, G, GV, US, haloshift=CS%halo_TS_diff) if (CS%double_diffuse) then - call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, G, GV, US, CS%set_diff_CSp, & - Kd_lay=Kd_lay, Kd_int=Kd_int, Kd_extra_T=Kd_extra_T, Kd_extra_S=Kd_extra_S) + call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, Kd_int, G, GV, US, & + CS%set_diff_CSp, Kd_lay=Kd_lay, Kd_extra_T=Kd_extra_T, Kd_extra_S=Kd_extra_S) else - call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, G, GV, US, & - CS%set_diff_CSp, Kd_lay=Kd_lay, Kd_int=Kd_int) + call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, Kd_int, G, GV, US, & + CS%set_diff_CSp, Kd_lay=Kd_lay) endif call cpu_clock_end(id_clock_set_diffusivity) if (showCallTree) call callTree_waypoint("done with set_diffusivity (diabatic)") @@ -1859,7 +1866,7 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e ! Add vertical diff./visc. due to convection (computed via CVMix) if (CS%use_CVMix_conv) then - call calculate_CVMix_conv(h, tv, G, GV, US, CS%CVMix_conv_CSp, Hml, Kd=Kd_int, Kv=visc%Kv_slow) + call calculate_CVMix_conv(h, tv, G, GV, US, CS%CVMix_conv_CSp, Hml, Kd_int, visc%Kv_slow) endif if (CS%useKPP) then @@ -3194,7 +3201,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di cmor_standard_name='ocean_vertical_heat_diffusivity', & cmor_long_name='Ocean vertical heat diffusivity') CS%id_Kd_salt = register_diag_field('ocean_model', 'Kd_salt', diag%axesTi, Time, & - 'Total diapycnal diffusivity for salt at interfaces', 'm2 s-1', conversion=US%Z2_T_to_m2_s, & + 'Total diapycnal diffusivity for salt at interfaces', 'm2 s-1', conversion=US%Z2_T_to_m2_s, & cmor_field_name='difvso', & cmor_standard_name='ocean_vertical_salt_diffusivity', & cmor_long_name='Ocean vertical salt diffusivity') @@ -3205,8 +3212,6 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di if (CS%useKPP) then allocate(CS%KPP_NLTheat(isd:ied,jsd:jed,nz+1), source=0.0) allocate(CS%KPP_NLTscalar(isd:ied,jsd:jed,nz+1), source=0.0) - endif - if (CS%useKPP) then allocate(CS%KPP_buoy_flux(isd:ied,jsd:jed,nz+1), source=0.0) allocate(CS%KPP_temp_flux(isd:ied,jsd:jed), source=0.0) allocate(CS%KPP_salt_flux(isd:ied,jsd:jed), source=0.0) diff --git a/src/parameterizations/vertical/MOM_full_convection.F90 b/src/parameterizations/vertical/MOM_full_convection.F90 index ceb77b52b8..aa1dfbf809 100644 --- a/src/parameterizations/vertical/MOM_full_convection.F90 +++ b/src/parameterizations/vertical/MOM_full_convection.F90 @@ -18,8 +18,7 @@ module MOM_full_convection contains !> Calculate new temperatures and salinities that have been subject to full convective mixing. -subroutine full_convection(G, GV, US, h, tv, T_adj, S_adj, p_surf, Kddt_smooth, & - Kddt_convect, halo) +subroutine full_convection(G, GV, US, h, tv, T_adj, S_adj, p_surf, Kddt_smooth, halo) type(ocean_grid_type), intent(in) :: 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 @@ -34,9 +33,7 @@ subroutine full_convection(G, GV, US, h, tv, T_adj, S_adj, p_surf, Kddt_smooth, real, dimension(:,:), pointer :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa] (or NULL). real, intent(in) :: Kddt_smooth !< A smoothing vertical !! diffusivity times a timestep [H2 ~> m2 or kg2 m-4]. - real, optional, intent(in) :: Kddt_convect !< A large convecting vertical - !! diffusivity times a timestep [H2 ~> m2 or kg2 m-4]. - integer, optional, intent(in) :: halo !< Halo width over which to compute + integer, intent(in) :: halo !< Halo width over which to compute ! Local variables real, dimension(SZI_(G),SZK_(GV)+1) :: & @@ -46,61 +43,53 @@ subroutine full_convection(G, GV, US, h, tv, T_adj, S_adj, p_surf, Kddt_smooth, ! in roundoff and can be neglected [H ~> m or kg m-2]. ! logical :: use_EOS ! If true, density is calculated from T & S using an equation of state. real, dimension(SZI_(G),SZK0_(G)) :: & - Te_a, & ! A partially updated temperature estimate including the influnce from + Te_a, & ! A partially updated temperature estimate including the influence from ! mixing with layers above rescaled by a factor of d_a [degC]. - ! This array is discreted on tracer cells, but contains an extra + ! This array is discretized on tracer cells, but contains an extra ! layer at the top for algorithmic convenience. - Se_a ! A partially updated salinity estimate including the influnce from + Se_a ! A partially updated salinity estimate including the influence from ! mixing with layers above rescaled by a factor of d_a [ppt]. - ! This array is discreted on tracer cells, but contains an extra + ! This array is discretized on tracer cells, but contains an extra ! layer at the top for algorithmic convenience. real, dimension(SZI_(G),SZK_(GV)+1) :: & - Te_b, & ! A partially updated temperature estimate including the influnce from + Te_b, & ! A partially updated temperature estimate including the influence from ! mixing with layers below rescaled by a factor of d_b [degC]. - ! This array is discreted on tracer cells, but contains an extra + ! This array is discretized on tracer cells, but contains an extra ! layer at the bottom for algorithmic convenience. - Se_b ! A partially updated salinity estimate including the influnce from + Se_b ! A partially updated salinity estimate including the influence from ! mixing with layers below rescaled by a factor of d_b [ppt]. - ! This array is discreted on tracer cells, but contains an extra + ! This array is discretized on tracer cells, but contains an extra ! layer at the bottom for algorithmic convenience. real, dimension(SZI_(G),SZK_(GV)+1) :: & c_a, & ! The fractional influence of the properties of the layer below - ! in the final properies with a downward-first solver, nondim. + ! in the final properties with a downward-first solver [nondim] d_a, & ! The fractional influence of the properties of the layer in question - ! and layers above in the final properies with a downward-first solver, nondim. + ! and layers above in the final properties with a downward-first solver [nondim] ! d_a = 1.0 - c_a c_b, & ! The fractional influence of the properties of the layer above - ! in the final properies with a upward-first solver, nondim. + ! in the final properties with a upward-first solver [nondim] d_b ! The fractional influence of the properties of the layer in question - ! and layers below in the final properies with a upward-first solver, nondim. + ! and layers below in the final properties with a upward-first solver [nondim] ! d_b = 1.0 - c_b real, dimension(SZI_(G),SZK_(GV)+1) :: & mix !< The amount of mixing across the interface between layers [H ~> m or kg m-2]. real :: mix_len ! The length-scale of mixing, when it is active [H ~> m or kg m-2] - real :: h_b, h_a ! The thicknessses of the layers above and below an interface [H ~> m or kg m-2] + real :: h_b, h_a ! The thicknesses of the layers above and below an interface [H ~> m or kg m-2] real :: b_b, b_a ! Inverse pivots used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1]. - real :: kap_dt_x2 ! The product of 2*kappa*dt [H2 ~> m2 or kg2 m-4]. - logical, dimension(SZI_(G)) :: do_i ! Do more work on this column. logical, dimension(SZI_(G)) :: last_down ! The last setup pass was downward. integer, dimension(SZI_(G)) :: change_ct ! The number of interfaces where the ! mixing has changed this iteration. - integer :: changed_col ! The number of colums whose mixing changed. + integer :: changed_col ! The number of columns whose mixing changed. integer :: i, j, k, is, ie, js, je, nz, itt - if (present(halo)) then - is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo - else - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec - endif + is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo nz = GV%ke if (.not.associated(tv%eqn_of_state)) return h_neglect = GV%H_subroundoff - kap_dt_x2 = 0.0 - if (present(Kddt_convect)) kap_dt_x2 = 2.0*Kddt_convect mix_len = (1.0e20 * nz) * (G%max_depth * GV%Z_to_H) h0 = 1.0e-16*sqrt(Kddt_smooth) + h_neglect @@ -135,7 +124,6 @@ subroutine full_convection(G, GV, US, h, tv, T_adj, S_adj, p_surf, Kddt_smooth, Te_a(i,k-2), Te_b(i,k+1), Se_a(i,k-2), Se_b(i,k+1), & d_a(i,K-1), d_b(i,K+1))) then mix(i,K) = mix_len - if (kap_dt_x2 > 0.0) mix(i,K) = kap_dt_x2 / ((h(i,j,k-1)+h(i,j,k)) + h0) change_ct(i) = change_ct(i) + 1 endif endif @@ -178,7 +166,6 @@ subroutine full_convection(G, GV, US, h, tv, T_adj, S_adj, p_surf, Kddt_smooth, Te_a(i,k-2), Te_b(i,k+1), Se_a(i,k-2), Se_b(i,k+1), & d_a(i,K-1), d_b(i,K+1))) then mix(i,K) = mix_len - if (kap_dt_x2 > 0.0) mix(i,K) = kap_dt_x2 / ((h(i,j,k-1)+h(i,j,k)) + h0) change_ct(i) = change_ct(i) + 1 endif endif @@ -260,7 +247,7 @@ subroutine full_convection(G, GV, US, h, tv, T_adj, S_adj, p_surf, Kddt_smooth, k = 1 ! A hook for debugging. - ! The following set of expressions for the final values are derived from the the partial + ! The following set of expressions for the final values are derived from the partial ! updates for the estimated temperatures and salinities around an interface, then directly ! solving for the final temperatures and salinities. They are here for later reference ! and to document an intermediate step in the stability calculation. @@ -336,7 +323,7 @@ subroutine smoothed_dRdT_dRdS(h, tv, Kddt, dR_dT, dR_dS, G, GV, US, j, p_surf, h type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type integer, intent(in) :: j !< The j-point to work on. real, dimension(:,:), pointer :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa]. - integer, optional, intent(in) :: halo !< Halo width over which to compute + integer, intent(in) :: halo !< Halo width over which to compute ! Local variables real :: mix(SZI_(G),SZK_(GV)+1) ! The diffusive mixing length (kappa*dt)/dz @@ -352,14 +339,10 @@ subroutine smoothed_dRdT_dRdS(h, tv, Kddt, dR_dT, dR_dS, G, GV, US, j, p_surf, h real :: h_neglect, h0 ! Negligible thicknesses to allow for zero thicknesses, ! [H ~> m or kg m-2]. real :: h_tr ! The thickness at tracer points, plus h_neglect [H ~> m or kg m-2]. - integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state + integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer :: i, k, is, ie, nz - if (present(halo)) then - is = G%isc-halo ; ie = G%iec+halo - else - is = G%isc ; ie = G%iec - endif + is = G%isc-halo ; ie = G%iec+halo nz = GV%ke h_neglect = GV%H_subroundoff diff --git a/src/parameterizations/vertical/MOM_geothermal.F90 b/src/parameterizations/vertical/MOM_geothermal.F90 index 7944d4b89f..90da24f170 100644 --- a/src/parameterizations/vertical/MOM_geothermal.F90 +++ b/src/parameterizations/vertical/MOM_geothermal.F90 @@ -499,7 +499,7 @@ subroutine geothermal_init(Time, G, GV, US, param_file, diag, CS, useALEalgorith type(diag_ctrl), target, intent(inout) :: diag !< Structure used to regulate diagnostic output. type(geothermal_CS), pointer :: CS !< Pointer pointing to the module control !! structure. - logical, optional, intent(in) :: useALEalgorithm !< logical for whether to use ALE remapping + logical, intent(in) :: useALEalgorithm !< logical for whether to use ALE remapping ! This include declares and sets the variable "version". #include "version_variable.h" @@ -587,13 +587,13 @@ subroutine geothermal_init(Time, G, GV, US, param_file, diag, CS, useALEalgorith 'internal_heat_temp_tendency', diag%axesTL, Time, & 'Temperature tendency (in 3D) due to internal (geothermal) sources', & 'degC s-1', conversion=US%s_to_T, v_extensive=.true.) - if (present(useALEalgorithm)) then ; if (.not.useALEalgorithm) then + if (.not.useALEalgorithm) then ! Do not offer this diagnostic if heating will be in place. CS%id_internal_heat_h_tendency=register_diag_field('ocean_model', & 'internal_heat_h_tendency', diag%axesTL, Time, & 'Thickness tendency (in 3D) due to internal (geothermal) sources', & trim(thickness_units)//' s-1', conversion=GV%H_to_MKS*US%s_to_T, v_extensive=.true.) - endif ; endif + endif end subroutine geothermal_init diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90 index 033f717091..a44a7aee95 100644 --- a/src/parameterizations/vertical/MOM_kappa_shear.F90 +++ b/src/parameterizations/vertical/MOM_kappa_shear.F90 @@ -107,7 +107,7 @@ module MOM_kappa_shear !> Subroutine for calculating shear-driven diffusivity and TKE in tracer columns subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & - kv_io, dt, G, GV, US, CS, initialize_all) + kv_io, dt, G, GV, US, CS) type(ocean_grid_type), intent(in) :: 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 @@ -137,8 +137,6 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & real, intent(in) :: dt !< Time increment [T ~> s]. type(Kappa_shear_CS), pointer :: CS !< The control structure returned by a previous !! call to kappa_shear_init. - logical, optional, intent(in) :: initialize_all !< If present and false, the previous - !! value of kappa is used to start the iterations ! Local variables real, dimension(SZI_(G),SZK_(GV)) :: & @@ -168,8 +166,6 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & real :: dz_massless ! A layer thickness that is considered massless [Z ~> m]. logical :: use_temperature ! If true, temperature and salinity have been ! allocated and are being used as state variables. - logical :: new_kappa = .true. ! If true, ignore the value of kappa from the - ! last call to this subroutine. integer, dimension(SZK_(GV)+1) :: kc ! The index map between the original ! interfaces and the interfaces with massless layers @@ -180,14 +176,13 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & is = G%isc ; ie = G%iec; js = G%jsc ; je = G%jec ; nz = GV%ke - use_temperature = .false. ; if (associated(tv%T)) use_temperature = .true. - new_kappa = .true. ; if (present(initialize_all)) new_kappa = initialize_all + use_temperature = associated(tv%T) k0dt = dt*CS%kappa_0 dz_massless = 0.1*sqrt(k0dt) - !$OMP parallel do default(private) shared(js,je,is,ie,nz,h,u_in,v_in,use_temperature,new_kappa, & - !$OMP tv,G,GV,US,CS,kappa_io,dz_massless,k0dt,p_surf,dt,tke_io,kv_io) + !$OMP parallel do default(private) shared(js,je,is,ie,nz,h,u_in,v_in,use_temperature,tv,G,GV,US, & + !$OMP CS,kappa_io,dz_massless,k0dt,p_surf,dt,tke_io,kv_io) do j=js,je do k=1,nz ; do i=is,ie h_2d(i,k) = h(i,j,k)*GV%H_to_Z @@ -198,9 +193,6 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & enddo ; enddo ; else ; do k=1,nz ; do i=is,ie rho_2d(i,k) = GV%Rlay(k) ! Could be tv%Rho(i,j,k) ? enddo ; enddo ; endif - if (.not.new_kappa) then ; do K=1,nz+1 ; do i=is,ie - kappa_2d(i,K) = kappa_io(i,j,K) - enddo ; enddo ; endif !--------------------------------------- ! Work on each column. @@ -278,11 +270,7 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & ! Set the initial guess for kappa, here defined at interfaces. ! ---------------------------------------------------- - if (new_kappa) then - do K=1,nzc+1 ; kappa(K) = US%m2_s_to_Z2_T*1.0 ; enddo - else - do K=1,nzc+1 ; kappa(K) = kappa_2d(i,K) ; enddo - endif + do K=1,nzc+1 ; kappa(K) = US%m2_s_to_Z2_T*1.0 ; enddo call kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & dz, u0xdz, v0xdz, T0xdz, S0xdz, kappa_avg, & @@ -340,7 +328,7 @@ end subroutine Calculate_kappa_shear !> Subroutine for calculating shear-driven diffusivity and TKE in corner columns subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_io, tke_io, & - kv_io, dt, G, GV, US, CS, initialize_all) + kv_io, dt, G, GV, US, CS) type(ocean_grid_type), intent(in) :: 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 @@ -373,8 +361,6 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ real, intent(in) :: dt !< Time increment [T ~> s]. type(Kappa_shear_CS), pointer :: CS !< The control structure returned by a previous !! call to kappa_shear_init. - logical, optional, intent(in) :: initialize_all !< If present and false, the previous - !! value of kappa is used to start the iterations ! Local variables real, dimension(SZIB_(G),SZK_(GV)) :: & @@ -397,7 +383,7 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ tke, & ! The Turbulent Kinetic Energy per unit mass at an interface [Z2 T-2 ~> m2 s-2]. kappa_avg, & ! The time-weighted average of kappa [Z2 T-1 ~> m2 s-1]. tke_avg ! The time-weighted average of TKE [Z2 T-2 ~> m2 s-2]. - real :: f2 ! The squared Coriolis parameter of each column [T-2 ~> s-2]. + real :: f2 ! The squared Coriolis parameter of each column [T-2 ~> s-2]. real :: surface_pres ! The top surface pressure [R L2 T-2 ~> Pa]. real :: dz_in_lay ! The running sum of the thickness in a layer [Z ~> m]. @@ -407,8 +393,6 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ real :: I_Prandtl ! The inverse of the turbulent Prandtl number [nondim]. logical :: use_temperature ! If true, temperature and salinity have been ! allocated and are being used as state variables. - logical :: new_kappa = .true. ! If true, ignore the value of kappa from the - ! last call to this subroutine. logical :: do_i ! If true, work on this column. integer, dimension(SZK_(GV)+1) :: kc ! The index map between the original @@ -421,15 +405,14 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ ! Diagnostics that should be deleted? isB = G%isc-1 ; ieB = G%iecB ; jsB = G%jsc-1 ; jeB = G%jecB ; nz = GV%ke - use_temperature = .false. ; if (associated(tv%T)) use_temperature = .true. - new_kappa = .true. ; if (present(initialize_all)) new_kappa = initialize_all + use_temperature = associated(tv%T) k0dt = dt*CS%kappa_0 dz_massless = 0.1*sqrt(k0dt) I_Prandtl = 0.0 ; if (CS%Prandtl_turb > 0.0) I_Prandtl = 1.0 / CS%Prandtl_turb - !$OMP parallel do default(private) shared(jsB,jeB,isB,ieB,nz,h,u_in,v_in,use_temperature,new_kappa, & - !$OMP tv,G,GV,US,CS,kappa_io,dz_massless,k0dt,p_surf,dt,tke_io,kv_io,I_Prandtl) + !$OMP parallel do default(private) shared(jsB,jeB,isB,ieB,nz,h,u_in,v_in,use_temperature,tv,G,GV, & + !$OMP US,CS,kappa_io,dz_massless,k0dt,p_surf,dt,tke_io,kv_io,I_Prandtl) do J=JsB,JeB J2 = mod(J,2)+1 ; J2m1 = 3-J2 ! = mod(J-1,2)+1 @@ -467,9 +450,6 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ if (.not.use_temperature) then ; do k=1,nz ; do I=IsB,IeB rho_2d(I,k) = GV%Rlay(k) enddo ; enddo ; endif - if (.not.new_kappa) then ; do K=1,nz+1 ; do I=IsB,IeB - kappa_2d(I,K,J2) = kv_io(I,J,K) * I_Prandtl - enddo ; enddo ; endif !--------------------------------------- ! Work on each column. @@ -558,11 +538,7 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ ! ---------------------------------------------------- ! Set the initial guess for kappa, here defined at interfaces. ! ---------------------------------------------------- - if (new_kappa) then - do K=1,nzc+1 ; kappa(K) = US%m2_s_to_Z2_T*1.0 ; enddo - else - do K=1,nzc+1 ; kappa(K) = kappa_2d(I,K,J2) ; enddo - endif + do K=1,nzc+1 ; kappa(K) = US%m2_s_to_Z2_T*1.0 ; enddo call kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & dz, u0xdz, v0xdz, T0xdz, S0xdz, kappa_avg, & @@ -621,9 +597,8 @@ end subroutine Calc_kappa_shear_vertex !> This subroutine calculates shear-driven diffusivity and TKE in a single column -subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & - dz, u0xdz, v0xdz, T0xdz, S0xdz, kappa_avg, & - tke_avg, tv, CS, GV, US, I_Ld2_1d, dz_Int_1d) +subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, dz, & + u0xdz, v0xdz, T0xdz, S0xdz, kappa_avg, tke_avg, tv, CS, GV, US) type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZK_(GV)+1), & intent(inout) :: kappa !< The time-weighted average of kappa [Z2 T-1 ~> m2 s-1]. @@ -654,11 +629,6 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & type(Kappa_shear_CS), pointer :: CS !< The control structure returned by a previous !! call to kappa_shear_init. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(SZK_(GV)+1), & - optional, intent(out) :: I_Ld2_1d !< The inverse of the squared mixing length [Z-2 ~> m-2]. - real, dimension(SZK_(GV)+1), & - optional, intent(out) :: dz_Int_1d !< The extent of a finite-volume space surrounding an interface, - !! as used in calculating kappa and TKE [Z ~> m]. ! Local variables real, dimension(nzc) :: & @@ -858,9 +828,8 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & ! enddo ! This call just calculates N2 and S2. - call calculate_projected_state(kappa, u, v, T, Sal, 0.0, nzc, dz, I_dz_int, & - dbuoy_dT, dbuoy_dS, u, v, T, Sal, GV, US, & - N2=N2, S2=S2, vel_underflow=CS%vel_underflow) + call calculate_projected_state(kappa, u, v, T, Sal, 0.0, nzc, dz, I_dz_int, dbuoy_dT, dbuoy_dS, & + CS%vel_underflow, u, v, T, Sal, N2, S2, GV, US) ! ---------------------------------------------------- ! Iterate ! ---------------------------------------------------- @@ -923,9 +892,8 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & ! value of max_KS_it may be unimportant, especially if it is large ! enough. call calculate_projected_state(kappa_out, u, v, T, Sal, 0.5*dt_test, nzc, dz, I_dz_int, & - dbuoy_dT, dbuoy_dS, u_test, v_test, T_test, S_test, & - GV, US, N2, S2, ks_int=ks_kappa, ke_int=ke_kappa, & - vel_underflow=CS%vel_underflow) + dbuoy_dT, dbuoy_dS, CS%vel_underflow, u_test, v_test, & + T_test, S_test, N2, S2, GV, US, ks_int=ks_kappa, ke_int=ke_kappa) valid_dt = .true. Idtt = 1.0 / dt_test do K=max(ks_kappa-1,2),min(ke_kappa+1,nzc) @@ -956,9 +924,9 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & if ((dt_test < dt_rem) .and. valid_dt) then dt_inc = 0.5*dt_test do itt_dt=1,dt_refinements - call calculate_projected_state(kappa_out, u, v, T, Sal, 0.5*(dt_test+dt_inc), & - nzc, dz, I_dz_int, dbuoy_dT, dbuoy_dS, u_test, v_test, T_test, S_test, & - GV, US, N2, S2, ks_int=ks_kappa, ke_int=ke_kappa, vel_underflow=CS%vel_underflow) + call calculate_projected_state(kappa_out, u, v, T, Sal, 0.5*(dt_test+dt_inc), nzc, dz, & + I_dz_int, dbuoy_dT, dbuoy_dS, CS%vel_underflow, u_test, v_test, T_test, S_test, & + N2, S2, GV, US, ks_int=ks_kappa, ke_int=ke_kappa) valid_dt = .true. Idtt = 1.0 / (dt_test+dt_inc) do K=max(ks_kappa-1,2),min(ke_kappa+1,nzc) @@ -1006,9 +974,8 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & else ! call cpu_clock_begin(id_clock_project) call calculate_projected_state(kappa_out, u, v, T, Sal, dt_now, nzc, dz, I_dz_int, & - dbuoy_dT, dbuoy_dS, u_test, v_test, T_test, S_test, & - GV, US, N2=N2, S2=S2, ks_int=ks_kappa, ke_int=ke_kappa, & - vel_underflow=CS%vel_underflow) + dbuoy_dT, dbuoy_dS, CS%vel_underflow, u_test, v_test, & + T_test, S_test, N2, S2, GV, US, ks_int=ks_kappa, ke_int=ke_kappa) ! call cpu_clock_end(id_clock_project) ! call cpu_clock_begin(id_clock_KQ) @@ -1026,9 +993,8 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & ! call cpu_clock_begin(id_clock_project) call calculate_projected_state(kappa_mid, u, v, T, Sal, dt_now, nzc, dz, I_dz_int, & - dbuoy_dT, dbuoy_dS, u_test, v_test, T_test, S_test, & - GV, US, N2=N2, S2=S2, ks_int=ks_kappa, ke_int=ke_kappa, & - vel_underflow=CS%vel_underflow) + dbuoy_dT, dbuoy_dS, CS%vel_underflow, u_test, v_test, & + T_test, S_test, N2, S2, GV, US, ks_int=ks_kappa, ke_int=ke_kappa) ! call cpu_clock_end(id_clock_project) ! call cpu_clock_begin(id_clock_KQ) @@ -1050,9 +1016,9 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & if (dt_rem > 0.0) then ! Update the values of u, v, T, Sal, N2, and S2 for the next iteration. ! call cpu_clock_begin(id_clock_project) - call calculate_projected_state(kappa_mid, u, v, T, Sal, dt_now, nzc, & - dz, I_dz_int, dbuoy_dT, dbuoy_dS, u, v, T, Sal, & - GV, US, N2, S2, vel_underflow=CS%vel_underflow) + call calculate_projected_state(kappa_mid, u, v, T, Sal, dt_now, nzc, dz, I_dz_int, & + dbuoy_dT, dbuoy_dS, CS%vel_underflow, u, v, T, Sal, N2, S2, & + GV, US) ! call cpu_clock_end(id_clock_project) endif @@ -1060,25 +1026,13 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & enddo ! end itt loop - if (present(I_Ld2_1d)) then - do K=1,GV%ke+1 ; I_Ld2_1d(K) = 0.0 ; enddo - do K=2,nzc ; if (TKE(K) > 0.0) & - I_Ld2_1d(K) = I_L2_bdry(K) + (N2(K) / CS%lambda**2 + f2) / TKE(K) - enddo - endif - if (present(dz_Int_1d)) then - do K=1,nzc+1 ; dz_Int_1d(K) = dz_Int(K) ; enddo - do K=nzc+2,GV%ke ; dz_Int_1d(K) = 0.0 ; enddo - endif - end subroutine kappa_shear_column !> This subroutine calculates the velocities, temperature and salinity that !! the water column will have after mixing for dt with diffusivities kappa. It !! may also calculate the projected buoyancy frequency and shear. -subroutine calculate_projected_state(kappa, u0, v0, T0, S0, dt, nz, & - dz, I_dz_int, dbuoy_dT, dbuoy_dS, & - u, v, T, Sal, GV, US, N2, S2, ks_int, ke_int, vel_underflow) +subroutine calculate_projected_state(kappa, u0, v0, T0, S0, dt, nz, dz, I_dz_int, dbuoy_dT, dbuoy_dS, & + vel_under, u, v, T, Sal, N2, S2, GV, US, ks_int, ke_int) integer, intent(in) :: nz !< The number of layers (after eliminating massless !! layers?). real, dimension(nz+1), intent(in) :: kappa !< The diapycnal diffusivity at interfaces, @@ -1087,6 +1041,7 @@ subroutine calculate_projected_state(kappa, u0, v0, T0, S0, dt, nz, & real, dimension(nz), intent(in) :: v0 !< The initial meridional velocity [L T-1 ~> m s-1]. real, dimension(nz), intent(in) :: T0 !< The initial temperature [degC]. real, dimension(nz), intent(in) :: S0 !< The initial salinity [ppt]. + real, intent(in) :: dt !< The time step [T ~> s]. real, dimension(nz), intent(in) :: dz !< The grid spacing of layers [Z ~> m]. real, dimension(nz+1), intent(in) :: I_dz_int !< The inverse of the layer's thicknesses !! [Z-1 ~> m-1]. @@ -1094,36 +1049,30 @@ subroutine calculate_projected_state(kappa, u0, v0, T0, S0, dt, nz, & !! temperature [Z T-2 degC-1 ~> m s-2 degC-1]. real, dimension(nz+1), intent(in) :: dbuoy_dS !< The partial derivative of buoyancy with !! salinity [Z T-2 ppt-1 ~> m s-2 ppt-1]. - real, intent(in) :: dt !< The time step [T ~> s]. + real, intent(in) :: vel_under !< Any velocities that are smaller in magnitude + !! than this value are set to 0 [L T-1 ~> m s-1]. real, dimension(nz), intent(inout) :: u !< The zonal velocity after dt [L T-1 ~> m s-1]. real, dimension(nz), intent(inout) :: v !< The meridional velocity after dt [L T-1 ~> m s-1]. real, dimension(nz), intent(inout) :: T !< The temperature after dt [degC]. real, dimension(nz), intent(inout) :: Sal !< The salinity after dt [ppt]. + real, dimension(nz+1), intent(inout) :: N2 !< The buoyancy frequency squared at interfaces [T-2 ~> s-2]. + real, dimension(nz+1), intent(inout) :: S2 !< The squared shear at interfaces [T-2 ~> s-2]. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(nz+1), optional, & - intent(inout) :: N2 !< The buoyancy frequency squared at interfaces [T-2 ~> s-2]. - real, dimension(nz+1), optional, & - intent(inout) :: S2 !< The squared shear at interfaces [T-2 ~> s-2]. integer, optional, intent(in) :: ks_int !< The topmost k-index with a non-zero diffusivity. integer, optional, intent(in) :: ke_int !< The bottommost k-index with a non-zero !! diffusivity. - real, optional, intent(in) :: vel_underflow !< If present and true, any velocities that - !! are smaller in magnitude than this value are - !! set to 0 [L T-1 ~> m s-1]. ! Local variables real, dimension(nz+1) :: c1 real :: L2_to_Z2 ! A conversion factor from horizontal length units to vertical depth ! units squared [Z2 s2 T-2 m-2 ~> 1]. - real :: underflow_vel ! Velocities smaller in magnitude than underflow_vel are set to 0 [L T-1 ~> m s-1]. real :: a_a, a_b, b1, d1, bd1, b1nz_0 integer :: k, ks, ke ks = 1 ; ke = nz if (present(ks_int)) ks = max(ks_int-1,1) if (present(ke_int)) ke = min(ke_int,nz) - underflow_vel = 0.0 ; if (present(vel_underflow)) underflow_vel = vel_underflow if (ks > ke) return @@ -1166,51 +1115,49 @@ subroutine calculate_projected_state(kappa, u0, v0, T0, S0, dt, nz, & endif u(ke) = b1nz_0 * (dz(ke)*u0(ke) + a_a*u(ke-1)) v(ke) = b1nz_0 * (dz(ke)*v0(ke) + a_a*v(ke-1)) - if (abs(u(ke)) < underflow_vel) u(ke) = 0.0 - if (abs(v(ke)) < underflow_vel) v(ke) = 0.0 + if (abs(u(ke)) < vel_under) u(ke) = 0.0 + if (abs(v(ke)) < vel_under) v(ke) = 0.0 do k=ke-1,ks,-1 u(k) = u(k) + c1(k+1)*u(k+1) v(k) = v(k) + c1(k+1)*v(k+1) - if (abs(u(k)) < underflow_vel) u(k) = 0.0 - if (abs(v(k)) < underflow_vel) v(k) = 0.0 + if (abs(u(k)) < vel_under) u(k) = 0.0 + if (abs(v(k)) < vel_under) v(k) = 0.0 T(k) = T(k) + c1(k+1)*T(k+1) Sal(k) = Sal(k) + c1(k+1)*Sal(k+1) enddo else ! dt <= 0.0 do k=1,nz u(k) = u0(k) ; v(k) = v0(k) ; T(k) = T0(k) ; Sal(k) = S0(k) - if (abs(u(k)) < underflow_vel) u(k) = 0.0 - if (abs(v(k)) < underflow_vel) v(k) = 0.0 + if (abs(u(k)) < vel_under) u(k) = 0.0 + if (abs(v(k)) < vel_under) v(k) = 0.0 enddo endif - if (present(S2)) then - ! L2_to_Z2 = US%m_to_Z**2 * US%T_to_s**2 - L2_to_Z2 = US%L_to_Z**2 - S2(1) = 0.0 ; S2(nz+1) = 0.0 - if (ks > 1) & - S2(ks) = ((u(ks)-u0(ks-1))**2 + (v(ks)-v0(ks-1))**2) * (L2_to_Z2*I_dz_int(ks)**2) - do K=ks+1,ke - S2(K) = ((u(k)-u(k-1))**2 + (v(k)-v(k-1))**2) * (L2_to_Z2*I_dz_int(K)**2) - enddo - if (ke 1) & - N2(ks) = max(0.0, I_dz_int(ks) * & - (dbuoy_dT(ks) * (T0(ks-1)-T(ks)) + dbuoy_dS(ks) * (S0(ks-1)-Sal(ks)))) - do K=ks+1,ke - N2(K) = max(0.0, I_dz_int(K) * & - (dbuoy_dT(K) * (T(k-1)-T(k)) + dbuoy_dS(K) * (Sal(k-1)-Sal(k)))) - enddo - if (ke 1) & + S2(ks) = ((u(ks)-u0(ks-1))**2 + (v(ks)-v0(ks-1))**2) * (L2_to_Z2*I_dz_int(ks)**2) + do K=ks+1,ke + S2(K) = ((u(k)-u(k-1))**2 + (v(k)-v(k-1))**2) * (L2_to_Z2*I_dz_int(K)**2) + enddo + if (ke 1) & + N2(ks) = max(0.0, I_dz_int(ks) * & + (dbuoy_dT(ks) * (T0(ks-1)-T(ks)) + dbuoy_dS(ks) * (S0(ks-1)-Sal(ks)))) + do K=ks+1,ke + N2(K) = max(0.0, I_dz_int(K) * & + (dbuoy_dT(K) * (T(k-1)-T(k)) + dbuoy_dS(K) * (Sal(k-1)-Sal(k)))) + enddo + if (ke s]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & + intent(out) :: Kd_int !< Diapycnal diffusivity at each interface [Z2 T-1 ~> m2 s-1]. type(set_diffusivity_CS), pointer :: CS !< Module control structure. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & optional, intent(out) :: Kd_lay !< Diapycnal diffusivity of each layer [Z2 T-1 ~> m2 s-1]. - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & - optional, intent(out) :: Kd_int !< Diapycnal diffusivity at each interface [Z2 T-1 ~> m2 s-1]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & optional, intent(out) :: Kd_extra_T !< The extra diffusivity at interfaces of !! temperature due to double diffusion relative to @@ -302,7 +302,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & ! Set Kd_lay, Kd_int and Kv_slow to constant values, mostly to fill the halos. if (present(Kd_lay)) Kd_lay(:,:,:) = CS%Kd - if (present(Kd_int)) Kd_int(:,:,:) = CS%Kd + Kd_int(:,:,:) = CS%Kd if (present(Kd_extra_T)) Kd_extra_T(:,:,:) = 0.0 if (present(Kd_extra_S)) Kd_extra_S(:,:,:) = 0.0 if (associated(visc%Kv_slow)) visc%Kv_slow(:,:,:) = CS%Kv @@ -468,98 +468,69 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & ! Add the input turbulent diffusivity. if (CS%useKappaShear .or. CS%use_CVMix_shear) then - if (present(Kd_int)) then - do K=2,nz ; do i=is,ie - Kd_int_2d(i,K) = visc%Kd_shear(i,j,K) + 0.5 * (Kd_lay_2d(i,k-1) + Kd_lay_2d(i,k)) - enddo ; enddo - do i=is,ie - Kd_int_2d(i,1) = visc%Kd_shear(i,j,1) ! This isn't actually used. It could be 0. - Kd_int_2d(i,nz+1) = 0.0 - enddo - endif + do K=2,nz ; do i=is,ie + Kd_int_2d(i,K) = visc%Kd_shear(i,j,K) + 0.5 * (Kd_lay_2d(i,k-1) + Kd_lay_2d(i,k)) + enddo ; enddo + do i=is,ie + Kd_int_2d(i,1) = visc%Kd_shear(i,j,1) ! This isn't actually used. It could be 0. + Kd_int_2d(i,nz+1) = 0.0 + enddo do k=1,nz ; do i=is,ie Kd_lay_2d(i,k) = Kd_lay_2d(i,k) + 0.5 * (visc%Kd_shear(i,j,K) + visc%Kd_shear(i,j,K+1)) enddo ; enddo else - if (present(Kd_int)) then - do i=is,ie - Kd_int_2d(i,1) = Kd_lay_2d(i,1) ; Kd_int_2d(i,nz+1) = 0.0 - enddo - do K=2,nz ; do i=is,ie - Kd_int_2d(i,K) = 0.5 * (Kd_lay_2d(i,k-1) + Kd_lay_2d(i,k)) - enddo ; enddo - endif + do i=is,ie + Kd_int_2d(i,1) = Kd_lay_2d(i,1) ; Kd_int_2d(i,nz+1) = 0.0 + enddo + do K=2,nz ; do i=is,ie + Kd_int_2d(i,K) = 0.5 * (Kd_lay_2d(i,k-1) + Kd_lay_2d(i,k)) + enddo ; enddo endif - if (present(Kd_int)) then - ! Add the ML_Rad diffusivity. - if (CS%ML_radiation) & - call add_MLrad_diffusivity(h, fluxes, j, G, GV, US, CS, TKE_to_Kd, Kd_lay_2d, Kd_int_2d) - - ! Add the Nikurashin and / or tidal bottom-driven mixing - if (CS%use_tidal_mixing) & - call calculate_tidal_mixing(h, N2_bot, j, TKE_to_Kd, maxTKE, G, GV, US, CS%tidal_mixing_CSp, & - N2_lay, N2_int, Kd_lay_2d, Kd_int_2d, CS%Kd_max, visc%Kv_slow) - - ! This adds the diffusion sustained by the energy extracted from the flow by the bottom drag. - if (CS%bottomdraglaw .and. (CS%BBL_effic>0.0)) then - if (CS%use_LOTW_BBL_diffusivity) then - call add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, G, GV, US, CS, & - dd%Kd_BBL, Kd_lay_2d, Kd_int_2d) - else - call add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & - maxTKE, kb, G, GV, US, CS, Kd_lay_2d, Kd_int_2d, dd%Kd_BBL) - endif - endif + ! Add the ML_Rad diffusivity. + if (CS%ML_radiation) & + call add_MLrad_diffusivity(h, fluxes, j, Kd_int_2d, G, GV, US, CS, TKE_to_Kd, Kd_lay_2d) - if (CS%limit_dissipation) then - ! This calculates the dissipation ONLY from Kd calculated in this routine - ! dissip has units of W/m3 (= kg/m3 * m2/s * 1/s2) - ! 1) a global constant, - ! 2) a dissipation proportional to N (aka Gargett) and - ! 3) dissipation corresponding to a (nearly) constant diffusivity. - do K=2,nz ; do i=is,ie - dissip = max( CS%dissip_min, & ! Const. floor on dissip. - CS%dissip_N0 + CS%dissip_N1 * sqrt(N2_int(i,K)), & ! Floor aka Gargett - CS%dissip_N2 * N2_int(i,K)) ! Floor of Kd_min*rho0/F_Ri - Kd_int_2d(i,K) = max(Kd_int_2d(i,K) , & ! Apply floor to Kd - dissip * (CS%FluxRi_max / (GV%Rho0 * (N2_int(i,K) + Omega2)))) - enddo ; enddo - endif + ! Add the Nikurashin and / or tidal bottom-driven mixing + if (CS%use_tidal_mixing) & + call calculate_tidal_mixing(h, N2_bot, j, TKE_to_Kd, maxTKE, G, GV, US, CS%tidal_mixing_CSp, & + N2_lay, N2_int, Kd_lay_2d, Kd_int_2d, CS%Kd_max, visc%Kv_slow) - ! Optionally add a uniform diffusivity at the interfaces. - if (CS%Kd_add > 0.0) then ; do K=1,nz+1 ; do i=is,ie - Kd_int_2d(i,K) = Kd_int_2d(i,K) + CS%Kd_add - enddo ; enddo ; endif + ! This adds the diffusion sustained by the energy extracted from the flow by the bottom drag. + if (CS%bottomdraglaw .and. (CS%BBL_effic>0.0)) then + if (CS%use_LOTW_BBL_diffusivity) then + call add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Kd_int_2d, G, GV, US, CS, & + dd%Kd_BBL, Kd_lay_2d) + else + call add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & + maxTKE, kb, G, GV, US, CS, Kd_lay_2d, Kd_int_2d, dd%Kd_BBL) + endif + endif - ! Copy the 2-d slices into the 3-d array that is exported. - do K=1,nz+1 ; do i=is,ie - Kd_int(i,j,K) = Kd_int_2d(i,K) + if (CS%limit_dissipation) then + ! This calculates the dissipation ONLY from Kd calculated in this routine + ! dissip has units of W/m3 (= kg/m3 * m2/s * 1/s2) + ! 1) a global constant, + ! 2) a dissipation proportional to N (aka Gargett) and + ! 3) dissipation corresponding to a (nearly) constant diffusivity. + do K=2,nz ; do i=is,ie + dissip = max( CS%dissip_min, & ! Const. floor on dissip. + CS%dissip_N0 + CS%dissip_N1 * sqrt(N2_int(i,K)), & ! Floor aka Gargett + CS%dissip_N2 * N2_int(i,K)) ! Floor of Kd_min*rho0/F_Ri + Kd_int_2d(i,K) = max(Kd_int_2d(i,K) , & ! Apply floor to Kd + dissip * (CS%FluxRi_max / (GV%Rho0 * (N2_int(i,K) + Omega2)))) enddo ; enddo + endif - else ! Kd_int is not present. - - ! Add the ML_Rad diffusivity. - if (CS%ML_radiation) & - call add_MLrad_diffusivity(h, fluxes, j, G, GV, US, CS, TKE_to_Kd, Kd_lay_2d) - - ! Add the Nikurashin and / or tidal bottom-driven mixing - if (CS%use_tidal_mixing) & - call calculate_tidal_mixing(h, N2_bot, j, TKE_to_Kd, maxTKE, G, GV, US, CS%tidal_mixing_CSp, & - N2_lay, N2_int, Kd_lay_2d, Kd_max=CS%Kd_max, Kv=visc%Kv_slow) - - ! This adds the diffusion sustained by the energy extracted from the flow by the bottom drag. - if (CS%bottomdraglaw .and. (CS%BBL_effic>0.0)) then - if (CS%use_LOTW_BBL_diffusivity) then - call add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, G, GV, US, CS, & - dd%Kd_BBL, Kd_lay_2d) - else - call add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & - maxTKE, kb, G, GV, US, CS, Kd_lay_2d, Kd_BBL=dd%Kd_BBL) - endif - endif + ! Optionally add a uniform diffusivity at the interfaces. + if (CS%Kd_add > 0.0) then ; do K=1,nz+1 ; do i=is,ie + Kd_int_2d(i,K) = Kd_int_2d(i,K) + CS%Kd_add + enddo ; enddo ; endif - endif + ! Copy the 2-d slices into the 3-d array that is exported. + do K=1,nz+1 ; do i=is,ie + Kd_int(i,j,K) = Kd_int_2d(i,K) + enddo ; enddo if (CS%limit_dissipation) then ! This calculates the layer dissipation ONLY from Kd calculated in this routine @@ -1163,7 +1134,7 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & !! thermodynamic fields. type(forcing), intent(in) :: fluxes !< A structure of thermodynamic surface fluxes type(vertvisc_type), intent(in) :: visc !< Structure containing vertical viscosities, bottom - !! boundary layer properies, and related fields + !! boundary layer properties and related fields integer, intent(in) :: j !< j-index of row to work on real, dimension(SZI_(G),SZK_(GV)), intent(in) :: TKE_to_Kd !< The conversion rate between the TKE !! TKE dissipated within a layer and the @@ -1177,8 +1148,7 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & type(set_diffusivity_CS), pointer :: CS !< Diffusivity control structure real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: Kd_lay !< The diapycnal diffusivity in layers, !! [Z2 T-1 ~> m2 s-1]. - real, dimension(SZI_(G),SZK_(GV)+1), & - optional, intent(inout) :: Kd_int !< The diapycnal diffusivity at interfaces, + real, dimension(SZI_(G),SZK_(GV)+1), intent(inout) :: Kd_int !< The diapycnal diffusivity at interfaces, !! [Z2 T-1 ~> m2 s-1]. real, dimension(:,:,:), pointer :: Kd_BBL !< Interface BBL diffusivity [Z2 T-1 ~> m2 s-1]. @@ -1330,10 +1300,8 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & else Kd_lay(i,k) = (TKE_to_layer + TKE_Ray) * TKE_to_Kd(i,k) endif - if (present(Kd_int)) then - Kd_int(i,K) = Kd_int(i,K) + 0.5 * delta_Kd - Kd_int(i,K+1) = Kd_int(i,K+1) + 0.5 * delta_Kd - endif + Kd_int(i,K) = Kd_int(i,K) + 0.5 * delta_Kd + Kd_int(i,K+1) = Kd_int(i,K+1) + 0.5 * delta_Kd if (do_diag_Kd_BBL) then Kd_BBL(i,j,K) = Kd_BBL(i,j,K) + 0.5 * delta_Kd Kd_BBL(i,j,K+1) = Kd_BBL(i,j,K+1) + 0.5 * delta_Kd @@ -1357,10 +1325,8 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & delta_Kd = TKE_here * TKE_to_Kd(i,k) if (CS%Kd_max >= 0.0) delta_Kd = min(delta_Kd, CS%Kd_max) Kd_lay(i,k) = Kd_lay(i,k) + delta_Kd - if (present(Kd_int)) then - Kd_int(i,K) = Kd_int(i,K) + 0.5 * delta_Kd - Kd_int(i,K+1) = Kd_int(i,K+1) + 0.5 * delta_Kd - endif + Kd_int(i,K) = Kd_int(i,K) + 0.5 * delta_Kd + Kd_int(i,K+1) = Kd_int(i,K+1) + 0.5 * delta_Kd if (do_diag_Kd_BBL) then Kd_BBL(i,j,K) = Kd_BBL(i,j,K) + 0.5 * delta_Kd Kd_BBL(i,j,K+1) = Kd_BBL(i,j,K+1) + 0.5 * delta_Kd @@ -1386,8 +1352,8 @@ end subroutine add_drag_diffusivity !> Calculates a BBL diffusivity use a Prandtl number 1 diffusivity with a law of the !! wall turbulent viscosity, up to a BBL height where the energy used for mixing has !! consumed the mechanical TKE input. -subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, & - G, GV, US, CS, Kd_BBL, Kd_lay, Kd_int) +subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Kd_int, & + G, GV, US, CS, Kd_BBL, Kd_lay) type(ocean_grid_type), intent(in) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -1401,16 +1367,16 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, & !! thermodynamic fields. type(forcing), intent(in) :: fluxes !< Surface fluxes structure type(vertvisc_type), intent(in) :: visc !< Structure containing vertical viscosities, bottom - !! boundary layer properies, and related fields. + !! boundary layer properties and related fields. integer, intent(in) :: j !< j-index of row to work on real, dimension(SZI_(G),SZK_(GV)+1), & intent(in) :: N2_int !< Square of Brunt-Vaisala at interfaces [T-2 ~> s-2] + real, dimension(SZI_(G),SZK_(GV)+1), & + intent(inout) :: Kd_int !< Interface net diffusivity [Z2 T-1 ~> m2 s-1] type(set_diffusivity_CS), pointer :: CS !< Diffusivity control structure real, dimension(:,:,:), pointer :: Kd_BBL !< Interface BBL diffusivity [Z2 T-1 ~> m2 s-1] real, dimension(SZI_(G),SZK_(GV)), & optional, intent(inout) :: Kd_lay !< Layer net diffusivity [Z2 T-1 ~> m2 s-1] - real, dimension(SZI_(G),SZK_(GV)+1), & - optional, intent(inout) :: Kd_int !< Interface net diffusivity [Z2 T-1 ~> m2 s-1] ! Local variables real :: TKE_column ! net TKE input into the column [Z3 T-3 ~> m3 s-3] @@ -1537,7 +1503,7 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, & TKE_remaining = TKE_remaining - TKE_consumed ! Note this will be non-negative ! Add this BBL diffusivity to the model net diffusivity. - if (present(Kd_int)) Kd_int(i,K) = Kd_int(i,K) + Kd_wall + Kd_int(i,K) = Kd_int(i,K) + Kd_wall if (present(Kd_lay)) Kd_lay(i,k) = Kd_lay(i,k) + 0.5 * (Kd_wall + Kd_lower) Kd_lower = Kd_wall ! Store for next layer up. if (do_diag_Kd_BBL) Kd_BBL(i,j,K) = Kd_wall @@ -1547,7 +1513,7 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, & end subroutine add_LOTW_BBL_diffusivity !> This routine adds effects of mixed layer radiation to the layer diffusivities. -subroutine add_MLrad_diffusivity(h, fluxes, j, G, GV, US, CS, TKE_to_Kd, Kd_lay, Kd_int) +subroutine add_MLrad_diffusivity(h, fluxes, j, Kd_int, G, GV, US, CS, TKE_to_Kd, Kd_lay) type(ocean_grid_type), intent(in) :: 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 @@ -1555,6 +1521,8 @@ subroutine add_MLrad_diffusivity(h, fluxes, j, G, GV, US, CS, TKE_to_Kd, Kd_lay, intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] type(forcing), intent(in) :: fluxes !< Surface fluxes structure integer, intent(in) :: j !< The j-index to work on + real, dimension(SZI_(G),SZK_(GV)+1), intent(inout) :: Kd_int !< The diapycnal diffusivity at interfaces + !! [Z2 T-1 ~> m2 s-1]. type(set_diffusivity_CS), pointer :: CS !< Diffusivity control structure real, dimension(SZI_(G),SZK_(GV)), intent(in) :: TKE_to_Kd !< The conversion rate between the TKE !! TKE dissipated within a layer and the @@ -1563,9 +1531,6 @@ subroutine add_MLrad_diffusivity(h, fluxes, j, G, GV, US, CS, TKE_to_Kd, Kd_lay, !! [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1] real, dimension(SZI_(G),SZK_(GV)), & optional, intent(inout) :: Kd_lay !< The diapycnal diffusivity in layers [Z2 T-1 ~> m2 s-1]. - real, dimension(SZI_(G),SZK_(GV)+1), & - optional, intent(inout) :: Kd_int !< The diapycnal diffusivity at interfaces - !! [Z2 T-1 ~> m2 s-1]. ! This routine adds effects of mixed layer radiation to the layer diffusivities. @@ -1639,14 +1604,12 @@ subroutine add_MLrad_diffusivity(h, fluxes, j, G, GV, US, CS, TKE_to_Kd, Kd_lay, Kd_lay(i,k) = Kd_lay(i,k) + Kd_mlr_ml(i) endif ; enddo ; enddo endif - if (present(Kd_int)) then - do K=2,kml+1 ; do i=is,ie ; if (do_i(i)) then - Kd_int(i,K) = Kd_int(i,K) + Kd_mlr_ml(i) - endif ; enddo ; enddo - if (kml<=nz-1) then ; do i=is,ie ; if (do_i(i)) then - Kd_int(i,Kml+2) = Kd_int(i,Kml+2) + 0.5 * Kd_mlr_ml(i) - endif ; enddo ; endif - endif + do K=2,kml+1 ; do i=is,ie ; if (do_i(i)) then + Kd_int(i,K) = Kd_int(i,K) + Kd_mlr_ml(i) + endif ; enddo ; enddo + if (kml<=nz-1) then ; do i=is,ie ; if (do_i(i)) then + Kd_int(i,Kml+2) = Kd_int(i,Kml+2) + 0.5 * Kd_mlr_ml(i) + endif ; enddo ; endif do k=kml+2,nz-1 do_any = .false. @@ -1674,10 +1637,8 @@ subroutine add_MLrad_diffusivity(h, fluxes, j, G, GV, US, CS, TKE_to_Kd, Kd_lay, if (present(Kd_lay)) then Kd_lay(i,k) = Kd_lay(i,k) + Kd_mlr endif - if (present(Kd_int)) then - Kd_int(i,K) = Kd_int(i,K) + 0.5 * Kd_mlr - Kd_int(i,K+1) = Kd_int(i,K+1) + 0.5 * Kd_mlr - endif + Kd_int(i,K) = Kd_int(i,K) + 0.5 * Kd_mlr + Kd_int(i,K+1) = Kd_int(i,K+1) + 0.5 * Kd_mlr TKE_ml_flux(i) = TKE_ml_flux(i) * exp(-z1) if (TKE_ml_flux(i) * I_decay(i) < 0.1 * CS%Kd_min * Omega2) then @@ -1703,7 +1664,7 @@ subroutine set_BBL_TKE(u, v, h, fluxes, visc, G, GV, US, CS, OBC) intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] type(forcing), intent(in) :: fluxes !< A structure of thermodynamic surface fluxes type(vertvisc_type), intent(in) :: visc !< Structure containing vertical viscosities, bottom - !! boundary layer properies, and related fields. + !! boundary layer properties and related fields. type(set_diffusivity_CS), pointer :: CS !< Diffusivity control structure type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. @@ -1995,10 +1956,10 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ !! structure. type(int_tide_CS), pointer :: int_tide_CSp !< A pointer to the internal tides control !! structure - integer, optional, intent(out) :: halo_TS !< The halo size of tracer points that must be + integer, intent(out) :: halo_TS !< The halo size of tracer points that must be !! valid for the calculations in set_diffusivity. - logical, optional, intent(out) :: double_diffuse !< If present, this indicates whether - !! some version of double diffusion is being used. + logical, intent(out) :: double_diffuse !< This indicates whether some version + !! of double diffusion is being used. ! Local variables real :: decay_length @@ -2323,14 +2284,10 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ 'Double-diffusion density ratio', 'nondim') endif - if (present(halo_TS)) then - halo_TS = 0 - if (CS%Vertex_Shear) halo_TS = 1 - endif + halo_TS = 0 + if (CS%Vertex_Shear) halo_TS = 1 - if (present(double_diffuse)) then - double_diffuse = (CS%double_diffusion .or. CS%use_CVMix_ddiff) - endif + double_diffuse = (CS%double_diffusion .or. CS%use_CVMix_ddiff) end subroutine set_diffusivity_init