diff --git a/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 b/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 index 08a09dbe23..9743c7fa3f 100644 --- a/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 +++ b/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 @@ -205,7 +205,7 @@ module MOM_surface_forcing_gfdl !> This subroutine translates the Ice_ocean_boundary_type into a MOM !! thermodynamic forcing type, including changes of units, sign conventions, !! and putting the fields into arrays with MOM-standard halos. -subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, sfc_state) +subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, US, CS, sfc_state) type(ice_ocean_boundary_type), & target, intent(in) :: IOB !< An ice-ocean boundary type with fluxes to drive !! the ocean in a coupled model @@ -215,6 +215,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, sfc integer, dimension(4), intent(in) :: index_bounds !< The i- and j- size of the arrays in IOB. type(time_type), intent(in) :: Time !< The time of the fluxes, used for interpolating the !! salinity to the right time, when it is being restored. + real, intent(in) :: valid_time !< The amount of time over which these fluxes + !! should be applied [s]. type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(surface_forcing_CS),pointer :: CS !< A pointer to the control structure returned by a @@ -307,7 +309,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, sfc if (CS%restore_temp) call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed) - fluxes%dt_buoy_accum = 0.0 endif ! endif for allocation and initialization @@ -324,11 +325,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, sfc ! ocean model, rather than using haloless arrays, in which case the last line ! would be: ( (/isd,is,ie,ied/), (/jsd,js,je,jed/)) - if (CS%allow_flux_adjustments) then - fluxes%heat_added(:,:)=0.0 - fluxes%salt_flux_added(:,:)=0.0 - endif - ! allocation and initialization on first call to this routine if (CS%area_surf < 0.0) then do j=js,je ; do i=is,ie @@ -337,6 +333,16 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, sfc CS%area_surf = reproducing_sum(work_sum, isr, ier, jsr, jer) endif ! endif for allocation and initialization + + ! Indicate that there are new unused fluxes. + fluxes%fluxes_used = .false. + fluxes%dt_buoy_accum = US%s_to_T*valid_time + + if (CS%allow_flux_adjustments) then + fluxes%heat_added(:,:)=0.0 + fluxes%salt_flux_added(:,:)=0.0 + endif + do j=js,je ; do i=is,ie fluxes%salt_flux(i,j) = 0.0 fluxes%vprec(i,j) = 0.0 diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 9982754053..1f01845ae4 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -513,7 +513,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda if (do_thermo) then if (OS%fluxes%fluxes_used) then - call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%fluxes, index_bnds, OS%Time, & + call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%fluxes, index_bnds, OS%Time, dt_coupling, & OS%grid, OS%US, OS%forcing_CSp, OS%sfc_state) ! Add ice shelf fluxes @@ -528,14 +528,11 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda call MOM_generic_tracer_fluxes_accumulate(OS%fluxes, 1.0) ! Here weight=1, so just store the current fluxes call disable_averaging(OS%diag) #endif - ! Indicate that there are new unused fluxes. - OS%fluxes%fluxes_used = .false. - OS%fluxes%dt_buoy_accum = dt_coupling else ! The previous fluxes have not been used yet, so translate the input fluxes ! into a temporary type and then accumulate them in about 20 lines. OS%flux_tmp%C_p = OS%fluxes%C_p - call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%flux_tmp, index_bnds, OS%Time, & + call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%flux_tmp, index_bnds, OS%Time, dt_coupling, & OS%grid, OS%US, OS%forcing_CSp, OS%sfc_state) if (OS%use_ice_shelf) & @@ -544,7 +541,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda call iceberg_fluxes(OS%grid, OS%US, OS%flux_tmp, OS%use_ice_shelf, & OS%sfc_state, dt_coupling, OS%marine_ice_CSp) - call fluxes_accumulate(OS%flux_tmp, OS%fluxes, dt_coupling, OS%grid, weight) + call fluxes_accumulate(OS%flux_tmp, OS%fluxes, OS%grid, weight) #ifdef _USE_GENERIC_TRACER ! Incorporate the current tracer fluxes into the running averages call MOM_generic_tracer_fluxes_accumulate(OS%flux_tmp, weight) @@ -646,16 +643,11 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda if (do_thermo) OS%nstep_thermo = OS%nstep_thermo + 1 if (do_dyn) then - call enable_averaging(dt_coupling, OS%Time_dyn, OS%diag) - call mech_forcing_diags(OS%forces, dt_coupling, OS%grid, OS%diag, OS%forcing_CSp%handles) - call disable_averaging(OS%diag) + call mech_forcing_diags(OS%forces, dt_coupling, OS%grid, OS%Time_dyn, OS%diag, OS%forcing_CSp%handles) endif if (OS%fluxes%fluxes_used .and. do_thermo) then - call enable_averaging(OS%fluxes%dt_buoy_accum, OS%Time, OS%diag) - call forcing_diagnostics(OS%fluxes, OS%sfc_state, OS%fluxes%dt_buoy_accum, & - OS%grid, OS%US, OS%diag, OS%forcing_CSp%handles) - call disable_averaging(OS%diag) + call forcing_diagnostics(OS%fluxes, OS%sfc_state, OS%grid, OS%US, OS%Time, OS%diag, OS%forcing_CSp%handles) endif ! Translate state into Ocean. diff --git a/config_src/mct_driver/mom_ocean_model_mct.F90 b/config_src/mct_driver/mom_ocean_model_mct.F90 index 240766a8d2..63556c2750 100644 --- a/config_src/mct_driver/mom_ocean_model_mct.F90 +++ b/config_src/mct_driver/mom_ocean_model_mct.F90 @@ -515,7 +515,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & call enable_averaging(dt_coupling, OS%Time + Ocean_coupling_time_step, OS%diag) if (do_thermo) & - call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%fluxes, index_bnds, OS%Time, & + call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%fluxes, index_bnds, OS%Time, dt_coupling, & OS%grid, OS%US, OS%forcing_CSp, OS%sfc_state, & OS%restore_salinity, OS%restore_temp) @@ -543,16 +543,12 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & call MOM_generic_tracer_fluxes_accumulate(OS%fluxes, weight) !here weight=1, just saving the current fluxes #endif - ! Indicate that there are new unused fluxes. - OS%fluxes%fluxes_used = .false. - OS%fluxes%dt_buoy_accum = dt_coupling - else OS%flux_tmp%C_p = OS%fluxes%C_p if (do_thermo) & - call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%flux_tmp, index_bnds, OS%Time, & + call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%flux_tmp, index_bnds, OS%Time, dt_coupling, & OS%grid, OS%US, OS%forcing_CSp, OS%sfc_state, OS%restore_salinity,OS%restore_temp) if (OS%use_ice_shelf) then @@ -570,7 +566,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & OS%sfc_state, dt_coupling, OS%marine_ice_CSp) endif - call forcing_accumulate(OS%flux_tmp, OS%forces, OS%fluxes, dt_coupling, OS%grid, weight) + call forcing_accumulate(OS%flux_tmp, OS%forces, OS%fluxes, OS%grid, weight) ! Some of the fields that exist in both the forcing and mech_forcing types ! (e.g., ustar) are time-averages must be copied back to the forces type. @@ -669,15 +665,10 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & OS%Time = Master_time + Ocean_coupling_time_step OS%nstep = OS%nstep + 1 - call enable_averaging(dt_coupling, OS%Time, OS%diag) - call mech_forcing_diags(OS%forces, dt_coupling, OS%grid, OS%diag, OS%forcing_CSp%handles) - call disable_averaging(OS%diag) + call mech_forcing_diags(OS%forces, dt_coupling, OS%grid, OS%Time, OS%diag, OS%forcing_CSp%handles) if (OS%fluxes%fluxes_used) then - call enable_averaging(OS%fluxes%dt_buoy_accum, OS%Time, OS%diag) - call forcing_diagnostics(OS%fluxes, OS%sfc_state, OS%fluxes%dt_buoy_accum, & - OS%grid, OS%US, OS%diag, OS%forcing_CSp%handles) - call disable_averaging(OS%diag) + call forcing_diagnostics(OS%fluxes, OS%sfc_state, OS%grid, OS%US, OS%Time, OS%diag, OS%forcing_CSp%handles) endif ! Translate state into Ocean. diff --git a/config_src/mct_driver/mom_surface_forcing_mct.F90 b/config_src/mct_driver/mom_surface_forcing_mct.F90 index b487787a2e..981202eda8 100644 --- a/config_src/mct_driver/mom_surface_forcing_mct.F90 +++ b/config_src/mct_driver/mom_surface_forcing_mct.F90 @@ -192,7 +192,7 @@ module MOM_surface_forcing_mct !> This subroutine translates the Ice_ocean_boundary_type into a MOM !! thermodynamic forcing type, including changes of units, sign conventions, !! and putting the fields into arrays with MOM-standard halos. -subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, & +subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, US, CS, & sfc_state, restore_salt, restore_temp) type(ice_ocean_boundary_type), & @@ -205,6 +205,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, & integer, dimension(4), intent(in) :: index_bounds !< The i- and j- size of the arrays in IOB. type(time_type), intent(in) :: Time !< The time of the fluxes, used for interpolating the !! salinity to the right time, when it is being restored. + real, intent(in) :: valid_time !< The amount of time over which these fluxes + !! should be applied [s]. type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(surface_forcing_CS),pointer :: CS !< A pointer to the control structure returned by a @@ -309,7 +311,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, & if (restore_temp) call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed) - fluxes%dt_buoy_accum = 0.0 endif ! endif for allocation and initialization @@ -326,11 +327,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, & ! ocean model, rather than using haloless arrays, in which case the last line ! would be: ( (/isd,is,ie,ied/), (/jsd,js,je,jed/)) - if (CS%allow_flux_adjustments) then - fluxes%heat_added(:,:)=0.0 - fluxes%salt_flux_added(:,:)=0.0 - endif - ! allocation and initialization on first call to this routine if (CS%area_surf < 0.0) then do j=js,je ; do i=is,ie @@ -339,6 +335,16 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, & CS%area_surf = reproducing_sum(work_sum, isr, ier, jsr, jer) endif ! endif for allocation and initialization + + ! Indicate that there are new unused fluxes. + fluxes%fluxes_used = .false. + fluxes%dt_buoy_accum = US%s_to_T*valid_time + + if (CS%allow_flux_adjustments) then + fluxes%heat_added(:,:)=0.0 + fluxes%salt_flux_added(:,:)=0.0 + endif + do j=js,je ; do i=is,ie fluxes%salt_flux(i,j) = 0.0 fluxes%vprec(i,j) = 0.0 diff --git a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 index 1e13b8e536..240b576669 100644 --- a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 +++ b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 @@ -517,7 +517,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & if (OS%fluxes%fluxes_used) then if (do_thermo) & - call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%fluxes, index_bnds, OS%Time, & + call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%fluxes, index_bnds, OS%Time, dt_coupling, & OS%grid, OS%US, OS%forcing_CSp, OS%sfc_state, & OS%restore_salinity, OS%restore_temp) @@ -544,13 +544,10 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & call enable_averaging(dt_coupling, OS%Time + Ocean_coupling_time_step, OS%diag) !Is this needed? call MOM_generic_tracer_fluxes_accumulate(OS%fluxes, weight) !here weight=1, just saving the current fluxes #endif - ! Indicate that there are new unused fluxes. - OS%fluxes%fluxes_used = .false. - OS%fluxes%dt_buoy_accum = dt_coupling else OS%flux_tmp%C_p = OS%fluxes%C_p if (do_thermo) & - call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%flux_tmp, index_bnds, OS%Time, & + call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%flux_tmp, index_bnds, OS%Time, dt_coupling, & OS%grid, OS%US, OS%forcing_CSp, OS%sfc_state, OS%restore_salinity,OS%restore_temp) if (OS%use_ice_shelf) then @@ -568,7 +565,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & OS%sfc_state, dt_coupling, OS%marine_ice_CSp) endif - call forcing_accumulate(OS%flux_tmp, OS%forces, OS%fluxes, dt_coupling, OS%grid, weight) + call forcing_accumulate(OS%flux_tmp, OS%forces, OS%fluxes, OS%grid, weight) ! Some of the fields that exist in both the forcing and mech_forcing types ! (e.g., ustar) are time-averages must be copied back to the forces type. call copy_back_forcing_fields(OS%fluxes, OS%forces, OS%grid) @@ -664,15 +661,10 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & OS%Time = Master_time + Ocean_coupling_time_step OS%nstep = OS%nstep + 1 - call enable_averaging(dt_coupling, OS%Time, OS%diag) - call mech_forcing_diags(OS%forces, dt_coupling, OS%grid, OS%diag, OS%forcing_CSp%handles) - call disable_averaging(OS%diag) + call mech_forcing_diags(OS%forces, dt_coupling, OS%grid, OS%Time, OS%diag, OS%forcing_CSp%handles) if (OS%fluxes%fluxes_used) then - call enable_averaging(OS%fluxes%dt_buoy_accum, OS%Time, OS%diag) - call forcing_diagnostics(OS%fluxes, OS%sfc_state, OS%fluxes%dt_buoy_accum, & - OS%grid, US%US, OS%diag, OS%forcing_CSp%handles) - call disable_averaging(OS%diag) + call forcing_diagnostics(OS%fluxes, OS%sfc_state, OS%grid, OS%US, OS%Time, OS%diag, OS%forcing_CSp%handles) endif ! Translate state into Ocean. diff --git a/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 b/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 index 955e608ac4..270d4e9f4c 100644 --- a/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 +++ b/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 @@ -199,7 +199,7 @@ module MOM_surface_forcing_nuopc !> This subroutine translates the Ice_ocean_boundary_type into a MOM !! thermodynamic forcing type, including changes of units, sign conventions, !! and putting the fields into arrays with MOM-standard halos. -subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, & +subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, US, CS, & sfc_state, restore_salt, restore_temp) type(ice_ocean_boundary_type), & target, intent(in) :: IOB !< An ice-ocean boundary type with fluxes to drive @@ -210,6 +210,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, & integer, dimension(4), intent(in) :: index_bounds !< The i- and j- size of the arrays in IOB. type(time_type), intent(in) :: Time !< The time of the fluxes, used for interpolating the !! salinity to the right time, when it is being restored. + real, intent(in) :: valid_time !< The amount of time over which these fluxes + !! should be applied [s]. type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(surface_forcing_CS),pointer :: CS !< A pointer to the control structure returned by a @@ -314,10 +316,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, & if (restore_temp) call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed) - fluxes%dt_buoy_accum = 0.0 endif ! endif for allocation and initialization - if (((associated(IOB%ustar_berg) .and. (.not.associated(fluxes%ustar_berg))) & .or. (associated(IOB%area_berg) .and. (.not.associated(fluxes%area_berg)))) & .or. (associated(IOB%mass_berg) .and. (.not.associated(fluxes%mass_berg)))) & @@ -331,12 +331,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, & ! ocean model, rather than using haloless arrays, in which case the last line ! would be: ( (/isd,is,ie,ied/), (/jsd,js,je,jed/)) - - if (CS%allow_flux_adjustments) then - fluxes%heat_added(:,:)=0.0 - fluxes%salt_flux_added(:,:)=0.0 - endif - ! allocation and initialization on first call to this routine if (CS%area_surf < 0.0) then do j=js,je ; do i=is,ie @@ -345,6 +339,16 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, & CS%area_surf = reproducing_sum(work_sum, isr, ier, jsr, jer) endif ! endif for allocation and initialization + + ! Indicate that there are new unused fluxes. + fluxes%fluxes_used = .false. + fluxes%dt_buoy_accum = US%s_to_T*valid_time + + if (CS%allow_flux_adjustments) then + fluxes%heat_added(:,:)=0.0 + fluxes%salt_flux_added(:,:)=0.0 + endif + do j=js,je ; do i=is,ie fluxes%salt_flux(i,j) = 0.0 fluxes%vprec(i,j) = 0.0 diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index a6d6597c0e..cea90b5db4 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -133,12 +133,15 @@ program MOM_main ! if Time_step_ocean is not an exact ! representation of dt_forcing. real :: dt_forcing ! The coupling time step [s]. - real :: dt ! The baroclinic dynamics time step [s]. + real :: dt ! The nominal baroclinic dynamics time step [s]. real :: dt_off ! Offline time step [s]. integer :: ntstep ! The number of baroclinic dynamics time steps ! within dt_forcing. - real :: dt_therm - real :: dt_dyn, dtdia, t_elapsed_seg + real :: dt_therm ! The thermodynamic timestep [s] + real :: dt_dyn ! The actual dynamic timestep used [s]. The value of dt_dyn is + ! chosen so that dt_forcing is an integer multiple of dt_dyn. + real :: dtdia ! The diabatic timestep [s] + real :: t_elapsed_seg ! The elapsed time in this run segment [s] integer :: n, n_max, nts, n_last_thermo logical :: diabatic_first, single_step_call type(time_type) :: Time2, time_chg @@ -491,7 +494,7 @@ program MOM_main call add_shelf_forces(grid, US, Ice_shelf_CSp, forces) endif fluxes%fluxes_used = .false. - fluxes%dt_buoy_accum = dt_forcing + fluxes%dt_buoy_accum = US%s_to_T*dt_forcing if (use_waves) then call Update_Surface_Waves(grid, GV, US, time, time_step_ocean, waves_csp) @@ -573,16 +576,12 @@ program MOM_main call write_cputime(Time, ns+ntstep-1, nmax, write_CPU_CSp) endif ; endif - call enable_averaging(dt_forcing, Time, diag) - call mech_forcing_diags(forces, dt_forcing, grid, diag, surface_forcing_CSp%handles) - call disable_averaging(diag) + call mech_forcing_diags(forces, dt_forcing, grid, Time, diag, surface_forcing_CSp%handles) if (.not. offline_tracer_mode) then if (fluxes%fluxes_used) then - call enable_averaging(fluxes%dt_buoy_accum, Time, diag) - call forcing_diagnostics(fluxes, sfc_state, fluxes%dt_buoy_accum, grid, US, & + call forcing_diagnostics(fluxes, sfc_state, grid, US, Time, & diag, surface_forcing_CSp%handles) - call disable_averaging(diag) else call MOM_error(FATAL, "The solo MOM_driver is not yet set up to handle "//& "thermodynamic time steps that are longer than the coupling timestep.") diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 4b16730fee..690e5250db 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -626,7 +626,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_int_in, CS, & elseif (thermo_does_span_coupling) then dtdia = dt_therm if ((fluxes%dt_buoy_accum > 0.0) .and. (dtdia > time_interval) .and. & - (abs(US%s_to_T*fluxes%dt_buoy_accum - dtdia) > 1e-6*dtdia)) then + (abs(fluxes%dt_buoy_accum - dtdia) > 1e-6*dtdia)) then call MOM_error(FATAL, "step_MOM: Mismatch between long thermodynamic "//& "timestep and time over which buoyancy fluxes have been accumulated.") endif diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 9794070f20..3dd3af8fbf 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -7,6 +7,7 @@ module MOM_forcing_type use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE use MOM_diag_mediator, only : post_data, register_diag_field, register_scalar_field use MOM_diag_mediator, only : time_type, diag_ctrl, safe_alloc_alloc, query_averaging_enabled +use MOM_diag_mediator, only : enable_averages, enable_averaging, disable_averaging use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_EOS, only : calculate_density_derivs use MOM_file_parser, only : get_param, log_param, log_version, param_file_type @@ -161,7 +162,7 @@ module MOM_forcing_type logical :: fluxes_used = .true. !< If true, all of the heat, salt, and mass !! fluxes have been applied to the ocean. real :: dt_buoy_accum = -1.0 !< The amount of time over which the buoyancy fluxes - !! should be applied [s]. If negative, this forcing + !! should be applied [T ~> s]. If negative, this forcing !! type variable has not yet been inialized. real :: C_p !< heat capacity of seawater [J kg-1 degC-1]. @@ -410,7 +411,8 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt_in_T, & real :: Pen_sw_tot(SZI_(G)) ! sum across all bands of Pen_SW [degC H ~> degC m or degC kg m-2]. real :: pen_sw_tot_rate(SZI_(G)) ! Summed rate of shortwave heating across bands ! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1] - real :: Ih_limit ! inverse depth at which surface fluxes start to be limited [H-1 ~> m-1 or m2 kg-1] + real :: Ih_limit ! inverse depth at which surface fluxes start to be limited + ! or 0 for no limiting [H-1 ~> m-1 or m2 kg-1] real :: scale ! scale scales away fluxes if depth < FluxRescaleDepth real :: W_m2_to_H_T ! converts W/m^2 to H degC T-1 [degC H T-1 W-2 m2 ~> degC m3 J-1 or degC kg J-1] real :: RZ_T_to_W_m2_degC ! Converts mass fluxes to heat fluxes per degree temperature @@ -438,7 +440,7 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt_in_T, & if (present(pen_sw_bnd_rate)) do_PSWBR = .true. !}BGR - Ih_limit = 1.0 / FluxRescaleDepth + Ih_limit = 0.0 ; if (FluxRescaleDepth > 0.0) Ih_limit = 1.0 / FluxRescaleDepth RZ_T_to_W_m2_degC = fluxes%C_p*US%R_to_kg_m3*US%Z_to_m*US%s_to_T I_Cp = 1.0 / fluxes%C_p W_m2_to_H_T = 1.0 / (US%s_to_T * GV%H_to_kg_m2 * fluxes%C_p) @@ -488,8 +490,7 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt_in_T, & do i=is,ie - scale = 1.0 - if (htot(i)*Ih_limit < 1.0) scale = htot(i)*Ih_limit + scale = 1.0 ; if ((Ih_limit > 0.0) .and. (htot(i)*Ih_limit < 1.0)) scale = htot(i)*Ih_limit ! Convert the penetrating shortwave forcing to (K * H) and reduce fluxes for shallow depths. ! (H=m for Bouss, H=kg/m2 for non-Bouss) @@ -1877,40 +1878,38 @@ end subroutine register_forcing_type_diags !> Accumulate the forcing over time steps, taking input from a mechanical forcing type !! and a temporary forcing-flux type. -subroutine forcing_accumulate(flux_tmp, forces, fluxes, dt, G, wt2) +subroutine forcing_accumulate(flux_tmp, forces, fluxes, G, wt2) type(forcing), intent(in) :: flux_tmp !< A temporary structure with current !!thermodynamic forcing fields type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces type(forcing), intent(inout) :: fluxes !< A structure containing time-averaged !! thermodynamic forcing fields - real, intent(in) :: dt !< The elapsed time since the last call to this subroutine [s] - type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure - real, intent(out) :: wt2 !< The relative weight of the new fluxes + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + real, intent(out) :: wt2 !< The relative weight of the new fluxes ! This subroutine copies mechancal forcing from flux_tmp to fluxes and ! stores the time-weighted averages of the various buoyancy fluxes in fluxes, ! and increments the amount of time over which the buoyancy forcing should be ! applied, all via a call to fluxes accumulate. - call fluxes_accumulate(flux_tmp, fluxes, dt, G, wt2, forces) + call fluxes_accumulate(flux_tmp, fluxes, G, wt2, forces) end subroutine forcing_accumulate !> Accumulate the thermodynamic fluxes over time steps -subroutine fluxes_accumulate(flux_tmp, fluxes, dt, G, wt2, forces) +subroutine fluxes_accumulate(flux_tmp, fluxes, G, wt2, forces) type(forcing), intent(in) :: flux_tmp !< A temporary structure with current - !! thermodynamic forcing fields + !! thermodynamic forcing fields type(forcing), intent(inout) :: fluxes !< A structure containing time-averaged - !! thermodynamic forcing fields - real, intent(in) :: dt !< The elapsed time since the last call to this subroutine [s] - type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure - real, intent(out) :: wt2 !< The relative weight of the new fluxes + !! thermodynamic forcing fields + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + real, intent(out) :: wt2 !< The relative weight of the new fluxes type(mech_forcing), optional, intent(in) :: forces !< A structure with the driving mechanical forces ! This subroutine copies mechancal forcing from flux_tmp to fluxes and ! stores the time-weighted averages of the various buoyancy fluxes in fluxes, - ! and increments the amount of time over which the buoyancy forcing should be - ! applied. + ! and increments the amount of time over which the buoyancy forcing in fluxes should be + ! applied based on the time interval stored in flux_tmp. real :: wt1 integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, i0, j0 @@ -1921,13 +1920,13 @@ subroutine fluxes_accumulate(flux_tmp, fluxes, dt, G, wt2, forces) IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB - if (fluxes%dt_buoy_accum < 0) call MOM_error(FATAL, "forcing_accumulate: "//& + if (fluxes%dt_buoy_accum < 0) call MOM_error(FATAL, "fluxes_accumulate: "//& "fluxes must be initialzed before it can be augmented.") ! wt1 is the relative weight of the previous fluxes. - wt1 = fluxes%dt_buoy_accum / (fluxes%dt_buoy_accum + dt) - wt2 = 1.0 - wt1 ! = dt / (fluxes%dt_buoy_accum + dt) - fluxes%dt_buoy_accum = fluxes%dt_buoy_accum + dt + wt1 = fluxes%dt_buoy_accum / (fluxes%dt_buoy_accum + flux_tmp%dt_buoy_accum) + wt2 = 1.0 - wt1 ! = flux_tmp%dt_buoy_accum / (fluxes%dt_buoy_accum + flux_tmp%dt_buoy_accum) + fluxes%dt_buoy_accum = fluxes%dt_buoy_accum + flux_tmp%dt_buoy_accum ! Copy over the pressure fields and accumulate averages of ustar, either from the forcing ! type or from the temporary fluxes type. @@ -2198,11 +2197,12 @@ end subroutine copy_back_forcing_fields !> Offer mechanical forcing fields for diagnostics for those !! fields registered as part of register_forcing_type_diags. -subroutine mech_forcing_diags(forces, dt, G, diag, handles) +subroutine mech_forcing_diags(forces, dt, G, time_end, diag, handles) type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces - real, intent(in) :: dt !< time step + real, intent(in) :: dt !< time step for the forcing [s] type(ocean_grid_type), intent(in) :: G !< grid type - type(diag_ctrl), intent(in) :: diag !< diagnostic type + type(time_type), intent(in) :: time_end !< The end time of the diagnostic interval. + type(diag_ctrl), intent(inout) :: diag !< diagnostic type type(forcing_diags), intent(inout) :: handles !< diagnostic id for diag_manager integer :: i,j,is,ie,js,je @@ -2210,7 +2210,8 @@ subroutine mech_forcing_diags(forces, dt, G, diag, handles) call cpu_clock_begin(handles%id_clock_forcing) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec - if (query_averaging_enabled(diag)) then + call enable_averaging(dt, time_end, diag) + ! if (query_averaging_enabled(diag)) then if ((handles%id_taux > 0) .and. associated(forces%taux)) & call post_data(handles%id_taux, forces%taux, diag) @@ -2224,22 +2225,23 @@ subroutine mech_forcing_diags(forces, dt, G, diag, handles) if ((handles%id_area_berg > 0) .and. associated(forces%area_berg)) & call post_data(handles%id_area_berg, forces%area_berg, diag) - endif + ! endif + call disable_averaging(diag) call cpu_clock_end(handles%id_clock_forcing) end subroutine mech_forcing_diags !> Offer buoyancy forcing fields for diagnostics for those !! fields registered as part of register_forcing_type_diags. -subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, US, diag, handles) +subroutine forcing_diagnostics(fluxes, sfc_state, G, US, time_end, diag, handles) type(forcing), intent(in) :: fluxes !< A structure containing thermodynamic forcing fields type(surface), intent(in) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. - real, intent(in) :: dt !< time step type(ocean_grid_type), intent(in) :: G !< grid type - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - type(diag_ctrl), intent(in) :: diag !< diagnostic regulator + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(time_type), intent(in) :: time_end !< The end time of the diagnostic interval. + type(diag_ctrl), intent(inout) :: diag !< diagnostic regulator type(forcing_diags), intent(inout) :: handles !< diagnostic ids ! local @@ -2248,7 +2250,7 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, US, diag, handles) real :: ave_flux ! for diagnosing averaged boundary flux real :: C_p ! seawater heat capacity (J/(deg K * kg)) real :: RZ_T_conversion ! A combination of scaling factors for mass fluxes [kg T m-2 s-1 R-1 Z-1 ~> 1] - real :: I_dt ! inverse time step + real :: I_dt ! inverse time step [s-1] real :: ppt2mks ! conversion between ppt and mks integer :: i,j,is,ie,js,je @@ -2256,11 +2258,12 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, US, diag, handles) C_p = fluxes%C_p RZ_T_conversion = US%R_to_kg_m3*US%Z_to_m*US%s_to_T - I_dt = 1.0/dt + I_dt = 1.0 / (US%T_to_s*fluxes%dt_buoy_accum) ppt2mks = 1e-3 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec - if (query_averaging_enabled(diag)) then + call enable_averages(fluxes%dt_buoy_accum, time_end, diag) + ! if (query_averaging_enabled(diag)) then ! post the diagnostics for surface mass fluxes ================================== @@ -2796,7 +2799,8 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, US, diag, handles) if ((handles%id_ustar_ice_cover > 0) .and. associated(fluxes%ustar_shelf)) & call post_data(handles%id_ustar_ice_cover, fluxes%ustar_shelf, diag) - endif ! query_averaging_enabled + ! endif ! query_averaging_enabled + call disable_averaging(diag) call cpu_clock_end(handles%id_clock_forcing) end subroutine forcing_diagnostics diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index 668f185152..6affbab231 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -944,7 +944,7 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS) !! describe the surface state of the ocean. type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various !! thermodynamic variables. - real, intent(in) :: dt !< The amount of time over which to average [s]. + real, intent(in) :: dt !< The amount of time over which to average [T ~> s]. type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(Sum_output_CS), pointer :: CS !< The control structure returned by a previous call @@ -963,7 +963,6 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS) real :: heat_input ! The total heat added by boundary fluxes, integrated ! over a time step and summed over space [J]. real :: C_p ! The heat capacity of seawater [J degC-1 kg-1]. - real :: dt_in_T ! Time increment [T ~> s] real :: RZL2_to_kg ! A combination of scaling factors for mass [kg R-1 Z-1 L-2 ~> 1] type(EFP_type) :: & @@ -977,13 +976,12 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec C_p = fluxes%C_p RZL2_to_kg = US%L_to_m**2*US%R_to_kg_m3*US%Z_to_m - dt_in_T = US%s_to_T*dt FW_in(:,:) = 0.0 ; FW_input = 0.0 if (associated(fluxes%evap)) then if (associated(fluxes%lprec) .and. associated(fluxes%fprec)) then do j=js,je ; do i=is,ie - FW_in(i,j) = RZL2_to_kg * dt_in_T*G%areaT(i,j)*(fluxes%evap(i,j) + & + FW_in(i,j) = RZL2_to_kg * dt*G%areaT(i,j)*(fluxes%evap(i,j) + & (((fluxes%lprec(i,j) + fluxes%vprec(i,j)) + fluxes%lrunoff(i,j)) + & (fluxes%fprec(i,j) + fluxes%frunoff(i,j)))) enddo ; enddo @@ -994,7 +992,7 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS) endif if (associated(fluxes%seaice_melt)) then ; do j=js,je ; do i=is,ie - FW_in(i,j) = FW_in(i,j) + RZL2_to_kg*dt_in_T * & + FW_in(i,j) = FW_in(i,j) + RZL2_to_kg*dt * & G%areaT(i,j) * fluxes%seaice_melt(i,j) enddo ; enddo ; endif @@ -1002,18 +1000,18 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS) if (CS%use_temperature) then if (associated(fluxes%sw)) then ; do j=js,je ; do i=is,ie - heat_in(i,j) = heat_in(i,j) + dt*US%L_to_m**2*G%areaT(i,j) * (fluxes%sw(i,j) + & + heat_in(i,j) = heat_in(i,j) + dt*US%T_to_s*US%L_to_m**2*G%areaT(i,j) * (fluxes%sw(i,j) + & (fluxes%lw(i,j) + (fluxes%latent(i,j) + fluxes%sens(i,j)))) enddo ; enddo ; endif if (associated(fluxes%seaice_melt_heat)) then ; do j=js,je ; do i=is,ie - heat_in(i,j) = heat_in(i,j) + dt*US%L_to_m**2*G%areaT(i,j) * fluxes%seaice_melt_heat(i,j) + heat_in(i,j) = heat_in(i,j) + dt*US%T_to_s*US%L_to_m**2*G%areaT(i,j) * fluxes%seaice_melt_heat(i,j) enddo ; enddo ; endif ! smg: new code ! include heat content from water transport across ocean surface ! if (associated(fluxes%heat_content_lprec)) then ; do j=js,je ; do i=is,ie -! heat_in(i,j) = heat_in(i,j) + dt_in_T*RZL2_to_kg*G%areaT(i,j) * & +! heat_in(i,j) = heat_in(i,j) + dt*RZL2_to_kg*G%areaT(i,j) * & ! (fluxes%heat_content_lprec(i,j) + (fluxes%heat_content_fprec(i,j) & ! + (fluxes%heat_content_lrunoff(i,j) + (fluxes%heat_content_frunoff(i,j) & ! + (fluxes%heat_content_cond(i,j) + (fluxes%heat_content_vprec(i,j) & @@ -1043,7 +1041,7 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS) heat_in(i,j) = heat_in(i,j) + US%L_to_m**2*G%areaT(i,j) * tv%frazil(i,j) enddo ; enddo ; endif if (associated(fluxes%heat_added)) then ; do j=js,je ; do i=is,ie - heat_in(i,j) = heat_in(i,j) + dt*US%L_to_m**2*G%areaT(i,j)*fluxes%heat_added(i,j) + heat_in(i,j) = heat_in(i,j) + dt*US%T_to_s*US%L_to_m**2*G%areaT(i,j)*fluxes%heat_added(i,j) enddo ; enddo ; endif ! if (associated(sfc_state%sw_lost)) then ; do j=js,je ; do i=is,ie ! heat_in(i,j) = heat_in(i,j) - US%L_to_m**2*G%areaT(i,j) * sfc_state%sw_lost(i,j) @@ -1051,7 +1049,7 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS) if (associated(fluxes%salt_flux)) then ; do j=js,je ; do i=is,ie ! convert salt_flux from kg (salt)/(m^2 s) to ppt * [m s-1]. - salt_in(i,j) = RZL2_to_kg * dt_in_T * & + salt_in(i,j) = RZL2_to_kg * dt * & G%areaT(i,j)*(1000.0*fluxes%salt_flux(i,j)) enddo ; enddo ; endif endif diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90 index d315a18b16..9349cf06d7 100644 --- a/src/parameterizations/vertical/MOM_kappa_shear.F90 +++ b/src/parameterizations/vertical/MOM_kappa_shear.F90 @@ -34,10 +34,10 @@ module MOM_kappa_shear !> This control structure holds the parameters that regulate shear mixing type, public :: Kappa_shear_CS ; private real :: RiNo_crit !< The critical shear Richardson number for - !! shear-entrainment. The theoretical value is 0.25. + !! shear-entrainment [nondim]. The theoretical value is 0.25. !! The values found by Jackson et al. are 0.25-0.35. real :: Shearmix_rate !< A nondimensional rate scale for shear-driven - !! entrainment. The value given by Jackson et al. + !! entrainment [nondim]. The value given by Jackson et al. !! is 0.085-0.089. real :: FRi_curvature !< A constant giving the curvature of the function !! of the Richardson number that relates shear to @@ -50,15 +50,16 @@ module MOM_kappa_shear !! shear (i.e. proportional to |S|*tke) [nondim]. !! The values found by Jackson et al. are 0.14-0.12. real :: lambda !< The coefficient for the buoyancy length scale - !! in the kappa equation. Nondimensional. + !! in the kappa equation [nondim]. !! The values found by Jackson et al. are 0.82-0.81. real :: lambda2_N_S !< The square of the ratio of the coefficients of !! the buoyancy and shear scales in the diffusivity - !! equation, 0 to eliminate the shear scale. Nondim. + !! equation, 0 to eliminate the shear scale [nondim]. real :: TKE_bg !< The background level of TKE [Z2 T-2 ~> m2 s-2]. real :: kappa_0 !< The background diapycnal diffusivity [Z2 T-1 ~> m2 s-1]. - real :: kappa_tol_err !< The fractional error in kappa that is tolerated. - real :: Prandtl_turb !< Prandtl number used to convert Kd_shear into viscosity. + real :: kappa_trunc !< Diffusivities smaller than this are rounded to 0 [Z2 T-1 ~> m2 s-1]. + real :: kappa_tol_err !< The fractional error in kappa that is tolerated [nondim]. + real :: Prandtl_turb !< Prandtl number used to convert Kd_shear into viscosity [nondim]. integer :: nkml !< The number of layers in the mixed layer, as !! treated in this routine. If the pieces of the !! mixed layer are not to be treated collectively, @@ -67,6 +68,9 @@ module MOM_kappa_shear !! to estimate the instantaneous shear-driven mixing. integer :: max_KS_it !< The maximum number of iterations that may be used !! to estimate the time-averaged diffusivity. + logical :: dKdQ_iteration_bug !< 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. logical :: KS_at_vertex !< If true, do the calculations of the shear-driven mixing !! at the cell vertices (i.e., the vorticity points). logical :: eliminate_massless !< If true, massless layers are merged with neighboring @@ -734,7 +738,7 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & local_src_avg, & ! The time-integral of the local source [nondim]. tol_min, & ! Minimum tolerated ksrc for the corrector step [T-1 ~> s-1]. tol_max, & ! Maximum tolerated ksrc for the corrector step [T-1 ~> s-1]. - tol_chg, & ! The tolerated change integrated in time [nondim]. + tol_chg, & ! The tolerated kappa change integrated over a timestep [nondim]. dist_from_top, & ! The distance from the top surface [Z ~> m]. local_src ! The sum of all sources of kappa, including kappa_src and ! sources from the elliptic term [T-1 ~> s-1]. @@ -747,11 +751,14 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & ! [Pa Z-1 = kg m-1 s-2 Z-1 ~> kg m-2 s-2]. 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, tol2 ! ### Tolerances that need to be set better later. + real :: tol_dksrc ! Tolerance for the change in the kappa source within an iteration + ! relative to the local source [nondim]. + 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 - ! within an iteration. 0 < tol_dksrc_low < 1. + ! within an iteration [nondim]. 0 < tol_dksrc_low < 1. real :: Ri_crit ! The critical shear Richardson number for shear- - ! driven mixing. The theoretical value is 0.25. + ! driven mixing [nondim]. The theoretical value is 0.25. real :: dt_rem ! The remaining time to advance the solution [T ~> s]. real :: dt_now ! The time step used in the current iteration [T ~> s]. real :: dt_wt ! The fractional weight of the current iteration [nondim]. @@ -794,7 +801,7 @@ 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 are hard-coded for now. Perhaps these could be made dynamic later? + !### 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 dt_refinements = 5 ! Selected so that 1/2^dt_refinements < 1-tol_dksrc_low @@ -1464,7 +1471,7 @@ subroutine find_kappa_tke(N2, S2, kappa_in, Idz, dz_Int, I_L2_bdry, f2, & TKE_min = max(CS%TKE_bg, 1.0E-20*US%m_to_Z**2*US%T_to_s**2) Ri_crit = CS%Rino_crit Ilambda2 = 1.0 / CS%lambda**2 - kappa_trunc = 0.01*kappa0 ! ### CHANGE THIS HARD-WIRING LATER? + kappa_trunc = CS%kappa_trunc do_Newton = .false. ; abort_Newton = .false. tol_err = CS%kappa_tol_err Newton_err = 0.2 ! This initial value may be automatically reduced later. @@ -1691,11 +1698,13 @@ subroutine find_kappa_tke(N2, S2, kappa_in, Idz, dz_Int, I_L2_bdry, f2, & cK(K+1) = bK * Idz(k) cKcomp = bK * (Idz(k-1)*cKcomp + decay_term_k) ! = 1-cK(K+1) - !### The following expression appears to be dimensionally inconsistent in length. -RWH - dKdQ(K) = bK * (Idz(k-1)*dKdQ(K-1)*cQ(K) + & + if (CS%dKdQ_iteration_bug) then + dKdQ(K) = bK * (Idz(k-1)*dKdQ(K-1)*cQ(K) + & US%m_to_Z*(N2(K)*Ilambda2 + f2) * I_Q**2 * kappa(K) ) - ! I think that the second term needs to be multiplied by dz_Int(K): - ! dz_Int(K)*(N2(K)*Ilambda2 + f2) * I_Q**2 * kappa(K) ) + else + dKdQ(K) = bK * (Idz(k-1)*dKdQ(K-1)*cQ(K) + & + dz_Int(K)*(N2(K)*Ilambda2 + f2) * I_Q**2 * kappa(K) ) + endif dK(K) = bK * (kap_src + Idz(k-1)*dK(K-1) + Idz(k-1)*dKdQ(K-1)*dQ(K-1)) ! Truncate away negligibly small values of kappa. @@ -1822,11 +1831,9 @@ subroutine find_kappa_tke(N2, S2, kappa_in, Idz, dz_Int, I_L2_bdry, f2, & kap_src = dz_Int(K) * (K_src(K) - I_Ld2(K)*kappa_prev(K)) + & (Idz(k-1)*(kappa_prev(k-1)-kappa_prev(k)) - & Idz(k)*(kappa_prev(k)-kappa_prev(k+1))) - !### The last line of the following appears to be dimensionally inconsistent with the first two. - ! I think that the term on the last line needs to be multiplied by dz_Int(K). K_err_lin = -Idz(k-1)*(dK(K-1)-dK(K)) + Idz(k)*(dK(K)-dK(K+1)) + & dz_Int(K)*I_Ld2_debug(K)*dK(K) - kap_src - & - US%m_to_Z*(N2(K)*Ilambda2 + f2)*I_Q**2*kappa_prev(K) * dQ(K) + dz_Int(K)*(N2(K)*Ilambda2 + f2)*I_Q**2*kappa_prev(K) * dQ(K) tke_src = dz_Int(K) * ((kappa_prev(K) + kappa0)*S2(K) - & kappa_prev(K)*N2(K) - (TKE_prev(K) - q0)*TKE_decay(K)) - & @@ -1950,8 +1957,9 @@ function kappa_shear_init(Time, G, GV, US, param_file, diag, CS) ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mdl = "MOM_kappa_shear" ! This module's name. + real :: kappa_0_unscaled ! The value of kappa_0 in MKS units [m2 s-1] real :: KD_normal ! The KD of the main model, read here only as a parameter - ! for setting the default of KD_SMOOTH + ! for setting the default of KD_SMOOTH in MKS units [m2 s-1] if (associated(CS)) then call MOM_error(WARNING, "kappa_shear_init called with an associated "// & @@ -1995,7 +2003,11 @@ function kappa_shear_init(Time, G, GV, US, param_file, diag, CS) "The background diffusivity that is used to smooth the "//& "density and shear profiles before solving for the "//& "diffusivities. Defaults to value of KD.", & - units="m2 s-1", default=KD_normal, scale=US%m2_s_to_Z2_T) + units="m2 s-1", default=KD_normal, scale=US%m2_s_to_Z2_T, unscaled=kappa_0_unscaled) + call get_param(param_file, mdl, "KD_TRUNC_KAPPA_SHEAR", CS%kappa_trunc, & + "The value of shear-driven diffusivity that is considered negligible "//& + "and is rounded down to 0. The default is 1% of KD_KAPPA_SHEAR_0.", & + units="m2 s-1", default=0.01*kappa_0_unscaled, scale=US%m2_s_to_Z2_T) call get_param(param_file, mdl, "FRI_CURVATURE", CS%FRi_curvature, & "The nondimensional curvature of the function of the "//& "Richardson number in the kappa source term in the "//& @@ -2050,12 +2062,16 @@ 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, "DEBUG_KAPPA_SHEAR", CS%debug, & "If true, write debugging data for the kappa-shear code. \n"//& "Caution: this option is _very_ verbose and should only "//& "be used in single-column mode!", & default=.false., debuggingParam=.true.) - + call get_param(param_file, mdl, "KAPPA_SHEAR_ITER_BUG", CS%dKdQ_iteration_bug, & + "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.) ! 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/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index c612a1ceed..a6a23d2adf 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -100,6 +100,10 @@ module MOM_vert_friction !! calculation, perhaps based on a bulk Richardson !! number criterion, to determine the mixed layer !! thickness for viscosity. + logical :: answers_2018 !< If true, use the order of arithmetic and expressions that recover the + !! answers from the end of 2018. Otherwise, use expressions that do not + !! use an arbitary and hard-coded maximum viscous coupling coefficient + !! between layers. logical :: debug !< If true, write verbose checksums for debugging purposes. integer :: nkml !< The number of layers in the mixed layer. integer, pointer :: ntrunc !< The number of times the velocity has been @@ -363,7 +367,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & zDS = 0.0 stress = dt_Rho0 * forces%tauy(i,J) do k=1,nz - h_a = 0.5 * (h(i,J,k) + h(i,J+1,k)) + h_a = 0.5 * (h(i,J,k) + h(i,J+1,k)) + h_neglect hfr = 1.0 ; if ((zDS+h_a) > Hmix) hfr = (Hmix - zDS) / h_a v(i,J,k) = v(i,J,k) + I_Hmix * hfr * stress zDS = zDS + h_a ; if (zDS >= Hmix) exit @@ -678,8 +682,8 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) if (CS%bottomdraglaw) then ; do I=Isq,Ieq kv_bbl(I) = visc%Kv_bbl_u(I,j) - bbl_thick(I) = visc%bbl_thick_u(I,j) * GV%Z_to_H - if (do_i(I)) I_Hbbl(I) = 1.0 / (bbl_thick(I) + h_neglect) + bbl_thick(I) = visc%bbl_thick_u(I,j) * GV%Z_to_H + h_neglect + if (do_i(I)) I_Hbbl(I) = 1.0 / bbl_thick(I) enddo ; endif do k=1,nz ; do I=Isq,Ieq ; if (do_i(I)) then @@ -845,7 +849,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) if (CS%bottomdraglaw) then ; do i=is,ie kv_bbl(i) = visc%Kv_bbl_v(i,J) - bbl_thick(i) = visc%bbl_thick_v(i,J) * GV%Z_to_H + bbl_thick(i) = visc%bbl_thick_v(i,J) * GV%Z_to_H + h_neglect if (do_i(i)) I_Hbbl(i) = 1.0 / bbl_thick(i) enddo ; endif @@ -1081,13 +1085,13 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, real :: I_Hmix ! The inverse of the mixed layer thickness [H-1 ~> m-1 or m2 kg-1]. real :: a_ml ! The layer coupling coefficient across an interface in ! the mixed layer [Z T-1 ~> m s-1]. - real :: I_amax ! The inverse of the maximum coupling coefficient [T s-1 Z-1 ~> m-1].??? + real :: I_amax ! The inverse of the maximum coupling coefficient [T Z-1 ~> s m-1]. real :: temp1 ! A temporary variable [H Z ~> m2 or kg m-1] 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 :: z2 ! A copy of z_i [nondim] real :: topfn ! A function that is 1 at the top and small far from it [nondim] - real :: a_top ! A viscosity associated with the top boundary layer [Z2 T-1 ~> m2 s-1] + real :: kv_top ! A viscosity associated with the top boundary layer [Z2 T-1 ~> m2 s-1] logical :: do_shelf, do_OBCs integer :: i, k, is, ie, max_nk integer :: nz @@ -1104,7 +1108,11 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, ! The maximum coupling coefficent was originally introduced to avoid ! truncation error problems in the tridiagonal solver. Effectively, the 1e-10 ! sets the maximum coupling coefficient increment to 1e10 m per timestep. - I_amax = (1.0e-10*US%Z_to_m) * dt + if (CS%answers_2018) then + I_amax = (1.0e-10*US%Z_to_m) * dt + else + I_amax = 0.0 + endif do_shelf = .false. ; if (present(shelf)) do_shelf = shelf do_OBCs = .false. @@ -1130,12 +1138,12 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, if (CS%bottomdraglaw) then r = hvel(i,nz)*0.5 if (r < bbl_thick(i)) then - a_cpl(i,nz+1) = kv_bbl(i) / (I_amax*kv_bbl(i) + r*GV%H_to_Z) + a_cpl(i,nz+1) = kv_bbl(i) / (I_amax*kv_bbl(i) + (r+h_neglect)*GV%H_to_Z) else - a_cpl(i,nz+1) = kv_bbl(i) / (I_amax*kv_bbl(i) + bbl_thick(i)*GV%H_to_Z) + a_cpl(i,nz+1) = kv_bbl(i) / (I_amax*kv_bbl(i) + (bbl_thick(i)+h_neglect)*GV%H_to_Z) endif else - a_cpl(i,nz+1) = CS%Kvbbl / (0.5*hvel(i,nz)*GV%H_to_Z + I_amax*CS%Kvbbl) + a_cpl(i,nz+1) = CS%Kvbbl / ((0.5*hvel(i,nz)+h_neglect)*GV%H_to_Z + I_amax*CS%Kvbbl) endif endif ; enddo @@ -1198,9 +1206,9 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, Kv_tot(i,K) = Kv_tot(i,K) + (kv_bbl(i) - CS%Kv)*botfn r = 0.5*(hvel(i,k) + hvel(i,k-1)) if (r > bbl_thick(i)) then - h_shear = ((1.0 - botfn) * r + botfn*bbl_thick(i)) + h_shear = ((1.0 - botfn) * r + botfn*bbl_thick(i)) + h_neglect else - h_shear = r + h_shear = r + h_neglect endif else Kv_tot(i,K) = Kv_tot(i,K) + (CS%Kvbbl-CS%Kv)*botfn @@ -1216,10 +1224,10 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, do i=is,ie ; if (do_i(i)) then if (work_on_u) then kv_TBL(i) = visc%Kv_tbl_shelf_u(I,j) - tbl_thick(i) = visc%tbl_thick_shelf_u(I,j) * GV%Z_to_H + tbl_thick(i) = visc%tbl_thick_shelf_u(I,j) * GV%Z_to_H + h_neglect else kv_TBL(i) = visc%Kv_tbl_shelf_v(i,J) - tbl_thick(i) = visc%tbl_thick_shelf_v(i,J) * GV%Z_to_H + tbl_thick(i) = visc%tbl_thick_shelf_v(i,J) * GV%Z_to_H + h_neglect endif z_t(i) = 0.0 @@ -1227,7 +1235,7 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, if (0.5*hvel(i,1) > tbl_thick(i)) then a_cpl(i,1) = kv_TBL(i) / (tbl_thick(i)*GV%H_to_Z + I_amax*kv_TBL(i)) else - a_cpl(i,1) = kv_TBL(i) / (0.5*hvel(i,1)*GV%H_to_Z + I_amax*kv_TBL(i)) + a_cpl(i,1) = kv_TBL(i) / ((0.5*hvel(i,1)+h_neglect)*GV%H_to_Z + I_amax*kv_TBL(i)) endif endif ; enddo @@ -1237,13 +1245,13 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, r = 0.5*(hvel(i,k)+hvel(i,k-1)) if (r > tbl_thick(i)) then - h_shear = ((1.0 - topfn) * r + topfn*tbl_thick(i)) + h_shear = ((1.0 - topfn) * r + topfn*tbl_thick(i)) + h_neglect else - h_shear = r + h_shear = r + h_neglect endif - a_top = topfn * kv_TBL(i) - a_cpl(i,K) = a_cpl(i,K) + a_top / (h_shear*GV%H_to_Z + I_amax*a_top) + kv_top = topfn * kv_TBL(i) + a_cpl(i,K) = a_cpl(i,K) + kv_top / (h_shear*GV%H_to_Z + I_amax*kv_top) endif ; enddo ; enddo elseif (CS%dynamic_viscous_ML .or. (GV%nkml>0)) then max_nk = 0 @@ -1292,7 +1300,7 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, temp1 = (z_t(i)*h_ml(i) - z_t(i)*z_t(i))*GV%H_to_Z ! This viscosity is set to go to 0 at the mixed layer top and bottom (in a log-layer) ! and be further limited by rotation to give the natural Ekman length. - visc_ml = u_star(i) * 0.41 * (temp1*u_star(i)) / (absf(i)*temp1 + h_ml(i)*u_star(i)) + visc_ml = u_star(i) * 0.41 * (temp1*u_star(i)) / (absf(i)*temp1 + (h_ml(i)+h_neglect)*u_star(i)) a_ml = visc_ml / (0.25*(hvel(i,k)+hvel(i,k-1) + h_neglect) * GV%H_to_Z + 0.5*I_amax*visc_ml) ! Choose the largest estimate of a. if (a_ml > a_cpl(i,K)) a_cpl(i,K) = a_ml @@ -1339,7 +1347,7 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, US, CS maxvel = CS%maxvel truncvel = 0.9*maxvel H_report = 6.0 * GV%Angstrom_H - dt_Rho0 = (US%L_T_to_m_s*US%Z_to_m) * dt / (GV%Rho0) + dt_Rho0 = (US%L_T_to_m_s*US%Z_to_m) * dt / GV%Rho0 if (len_trim(CS%u_trunc_file) > 0) then !$OMP parallel do default(shared) private(trunc_any,CFL) @@ -1536,6 +1544,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & real :: hmix_str_dflt real :: Kv_dflt ! A default viscosity [m2 s-1]. real :: Hmix_m ! A boundary layer thickness [m]. + logical :: default_2018_answers integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz ! This include declares and sets the variable "version". #include "version_variable.h" @@ -1559,6 +1568,14 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & ! Default, read and log parameters call log_version(param_file, mdl, version, "") + 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, "VERT_FRICTION_2018_ANSWERS", CS%answers_2018, & + "If true, use the order of arithmetic and expressions that recover the answers "//& + "from the end of 2018. Otherwise, use expressions that do not use an arbitary "//& + "and hard-coded maximum viscous coupling coefficient between layers.", & + default=default_2018_answers) call get_param(param_file, mdl, "BOTTOMDRAGLAW", CS%bottomdraglaw, & "If true, the bottom stress is calculated with a drag "//& "law of the form c_drag*|u|*u. The velocity magnitude "//&