diff --git a/config_src/mct_driver/mom_ocean_model_mct.F90 b/config_src/mct_driver/mom_ocean_model_mct.F90 index f8a4a19532..2f94c9b7f9 100644 --- a/config_src/mct_driver/mom_ocean_model_mct.F90 +++ b/config_src/mct_driver/mom_ocean_model_mct.F90 @@ -11,56 +11,57 @@ module MOM_ocean_model_mct ! This code is a stop-gap wrapper of the MOM6 code to enable it to be called ! in the same way as MOM4. -use MOM, only : initialize_MOM, step_MOM, MOM_control_struct, MOM_end -use MOM, only : extract_surface_state, allocate_surface_state, finish_MOM_initialization -use MOM, only : get_MOM_state_elements, MOM_state_is_synchronized -use MOM, only : get_ocean_stocks, step_offline -use MOM_constants, only : CELSIUS_KELVIN_OFFSET, hlf -use MOM_diag_mediator, only : diag_ctrl, enable_averaging, disable_averaging -use MOM_diag_mediator, only : diag_mediator_close_registration, diag_mediator_end -use MOM_domains, only : pass_var, pass_vector, AGRID, BGRID_NE, CGRID_NE -use MOM_domains, only : TO_ALL, Omit_Corners -use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe -use MOM_error_handler, only : callTree_enter, callTree_leave -use MOM_file_parser, only : get_param, log_version, close_param_file, param_file_type -use MOM_forcing_type, only : allocate_forcing_type -use MOM_forcing_type, only : forcing, mech_forcing -use MOM_forcing_type, only : forcing_accumulate, copy_common_forcing_fields -use MOM_forcing_type, only : copy_back_forcing_fields, set_net_mass_forcing -use MOM_forcing_type, only : set_derived_forcing_fields -use MOM_forcing_type, only : forcing_diagnostics, mech_forcing_diags -use MOM_get_input, only : Get_MOM_Input, directories -use MOM_grid, only : ocean_grid_type -use MOM_io, only : close_file, file_exists, read_data, write_version_number -use MOM_marine_ice, only : iceberg_forces, iceberg_fluxes, marine_ice_init, marine_ice_CS -use MOM_restart, only : MOM_restart_CS, save_restart -use MOM_string_functions, only : uppercase -use MOM_surface_forcing_mct, only : surface_forcing_init, convert_IOB_to_fluxes -use MOM_surface_forcing_mct, only : convert_IOB_to_forces, ice_ocn_bnd_type_chksum -use MOM_surface_forcing_mct, only : ice_ocean_boundary_type, surface_forcing_CS -use MOM_surface_forcing_mct, only : forcing_save_restart -use MOM_time_manager, only : time_type, get_time, set_time, operator(>) -use MOM_time_manager, only : operator(+), operator(-), operator(*), operator(/) -use MOM_time_manager, only : operator(/=), operator(<=), operator(>=) -use MOM_time_manager, only : operator(<), real_to_time_type, time_type_to_real -use MOM_tracer_flow_control, only : call_tracer_register, tracer_flow_control_init -use MOM_tracer_flow_control, only : call_tracer_flux_init -use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : surface -use MOM_verticalGrid, only : verticalGrid_type -use MOM_ice_shelf, only : initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS -use MOM_ice_shelf, only : add_shelf_forces, ice_shelf_end, ice_shelf_save_restart -use coupler_types_mod, only : coupler_1d_bc_type, coupler_2d_bc_type -use coupler_types_mod, only : coupler_type_spawn, coupler_type_write_chksums -use coupler_types_mod, only : coupler_type_initialized, coupler_type_copy_data -use coupler_types_mod, only : coupler_type_set_diags, coupler_type_send_data -use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain -use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain -use fms_mod, only : stdout -use mpp_mod, only : mpp_chksum -use MOM_EOS, only : gsw_sp_from_sr, gsw_pt_from_ct -use MOM_wave_interface, only: wave_parameters_CS, MOM_wave_interface_init -use MOM_wave_interface, only: MOM_wave_interface_init_lite, Update_Surface_Waves +use MOM, only : initialize_MOM, step_MOM, MOM_control_struct, MOM_end +use MOM, only : extract_surface_state, allocate_surface_state, finish_MOM_initialization +use MOM, only : get_MOM_state_elements, MOM_state_is_synchronized +use MOM, only : get_ocean_stocks, step_offline +use MOM_constants, only : CELSIUS_KELVIN_OFFSET, hlf +use MOM_diag_mediator, only : diag_ctrl, enable_averaging, disable_averaging +use MOM_diag_mediator, only : diag_mediator_close_registration, diag_mediator_end +use MOM_domains, only : pass_var, pass_vector, AGRID, BGRID_NE, CGRID_NE +use MOM_domains, only : TO_ALL, Omit_Corners +use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe +use MOM_error_handler, only : callTree_enter, callTree_leave +use MOM_file_parser, only : get_param, log_version, close_param_file, param_file_type +use MOM_forcing_type, only : allocate_forcing_type +use MOM_forcing_type, only : forcing, mech_forcing +use MOM_forcing_type, only : forcing_accumulate, copy_common_forcing_fields +use MOM_forcing_type, only : copy_back_forcing_fields, set_net_mass_forcing +use MOM_forcing_type, only : set_derived_forcing_fields +use MOM_forcing_type, only : forcing_diagnostics, mech_forcing_diags +use MOM_get_input, only : Get_MOM_Input, directories +use MOM_grid, only : ocean_grid_type +use MOM_io, only : close_file, file_exists, read_data, write_version_number +use MOM_marine_ice, only : iceberg_forces, iceberg_fluxes, marine_ice_init, marine_ice_CS +use MOM_restart, only : MOM_restart_CS, save_restart +use MOM_string_functions, only : uppercase +use MOM_surface_forcing_mct, only : surface_forcing_init, convert_IOB_to_fluxes +use MOM_surface_forcing_mct, only : convert_IOB_to_forces, ice_ocn_bnd_type_chksum +use MOM_surface_forcing_mct, only : ice_ocean_boundary_type, surface_forcing_CS +use MOM_surface_forcing_mct, only : forcing_save_restart +use MOM_time_manager, only : time_type, get_time, set_time, operator(>) +use MOM_time_manager, only : operator(+), operator(-), operator(*), operator(/) +use MOM_time_manager, only : operator(/=), operator(<=), operator(>=) +use MOM_time_manager, only : operator(<), real_to_time_type, time_type_to_real +use MOM_tracer_flow_control, only : call_tracer_register, tracer_flow_control_init +use MOM_tracer_flow_control, only : call_tracer_flux_init +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : surface +use MOM_verticalGrid, only : verticalGrid_type +use MOM_ice_shelf, only : initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS +use MOM_ice_shelf, only : add_shelf_forces, ice_shelf_end, ice_shelf_save_restart +use coupler_types_mod, only : coupler_1d_bc_type, coupler_2d_bc_type +use coupler_types_mod, only : coupler_type_spawn, coupler_type_write_chksums +use coupler_types_mod, only : coupler_type_initialized, coupler_type_copy_data +use coupler_types_mod, only : coupler_type_set_diags, coupler_type_send_data +use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain +use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain +use fms_mod, only : stdout +use mpp_mod, only : mpp_chksum +use MOM_EOS, only : gsw_sp_from_sr, gsw_pt_from_ct +use MOM_wave_interface, only : wave_parameters_CS, MOM_wave_interface_init +use MOM_wave_interface, only : MOM_wave_interface_init_lite, Update_Surface_Waves +use time_interp_external_mod, only : time_interp_external_init ! MCT specfic routines use MOM_domains, only : MOM_infra_end @@ -265,6 +266,8 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i OS%is_ocean_pe = Ocean_sfc%is_ocean_pe if (.not.OS%is_ocean_pe) return + call time_interp_external_init + OS%Time = Time_in call initialize_MOM(OS%Time, Time_init, param_file, OS%dirs, OS%MOM_CSp, & OS%restart_CSp, Time_in, offline_tracer_mode=OS%offline_tracer_mode, & diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index b1ce9a60c0..741ce832e8 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -119,7 +119,9 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) integer :: year, month, day, hour, minute, seconds, seconds_n, seconds_d, rc character(len=240) :: runid !< Run ID character(len=32) :: runtype !< Run type - character(len=240) :: restartfile !< Path/Name of restart file + character(len=512) :: restartfile !< Path/Name of restart file + character(len=2048) :: restartfiles !< Path/Name of restart files. + !! (same as restartfile if a single restart file is to be read in) integer :: nu !< i/o unit to read pointer file character(len=240) :: restart_pointer_file !< File name for restart pointer file character(len=240) :: restartpath !< Path of the restart file @@ -164,6 +166,7 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) !logical :: lsend_precip_fact !< If T,send precip_fact to cpl for use in fw balance !! (partially-coupled option) character(len=128) :: err_msg !< Error message + integer :: iostat ! set the cdata pointers: call seq_cdata_setptrs(cdata_o, id=MOM_MCT_ID, mpicom=mpicom_ocn, & @@ -296,15 +299,27 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) nu = shr_file_getUnit() restart_pointer_file = trim(glb%pointer_filename) if (is_root_pe()) write(glb%stdout,*) 'Reading ocn pointer file: ',restart_pointer_file + restartfile = ""; restartfiles = ""; open(nu, file=restart_pointer_file, form='formatted', status='unknown') - read(nu,'(a)') restartfile + do + read(nu,'(a)', iostat=iostat) restartfile + if (len(trim(restartfiles))>1 .and. iostat<0) then + exit ! done reading restart files list. + else if (iostat/=0) then + call MOM_error(FATAL, 'Error reading rpointer.ocn') + endif + ! check if the length of restartfiles variable is sufficient: + if (len(restartfiles)-len(trim(restartfiles)) < len(trim(restartfile))) then + call MOM_error(FATAL, "Restart file name(s) too long.") + endif + restartfiles = trim(restartfiles) // " " // trim(restartfile) + enddo close(nu) - !restartfile = trim(restartpath) // trim(restartfile) if (is_root_pe()) then - write(glb%stdout,*) 'Reading restart file: ',trim(restartfile) + write(glb%stdout,*) 'Reading restart file(s): ',trim(restartfiles) end if call shr_file_freeUnit(nu) - call ocean_model_init(glb%ocn_public, glb%ocn_state, time0, time_start, input_restart_file=trim(restartfile)) + call ocean_model_init(glb%ocn_public, glb%ocn_state, time0, time_start, input_restart_file=trim(restartfiles)) endif if (is_root_pe()) then write(glb%stdout,'(/12x,a/)') '======== COMPLETED MOM INITIALIZATION ========' @@ -434,6 +449,9 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) integer :: ocn_cpl_dt !< one ocn coupling interval in seconds. (to be received from cesm) real (kind=8) :: mom_cpl_dt !< one ocn coupling interval in seconds. (internal) integer :: ncouple_per_day !< number of ocean coupled call in one day (non-dim) + integer :: num_rest_files !< number of restart files written + integer :: i + character(len=8) :: suffix ! reset shr logging to ocn log file: if (is_root_pe()) then @@ -534,7 +552,8 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) write(restartname,'(A,".mom6.r.",I4.4,"-",I2.2,"-",I2.2,"-",I5.5)') trim(runid), year, month, day, seconds call save_restart(glb%ocn_state%dirs%restart_output_dir, glb%ocn_state%Time, glb%grid, & - glb%ocn_state%restart_CSp, .false., filename=restartname, GV=glb%ocn_state%GV) + glb%ocn_state%restart_CSp, .false., filename=restartname, GV=glb%ocn_state%GV, & + num_rest_files=num_rest_files) ! write name of restart file in the rpointer file nu = shr_file_getUnit() @@ -542,6 +561,19 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) restart_pointer_file = trim(glb%pointer_filename) open(nu, file=restart_pointer_file, form='formatted', status='unknown') write(nu,'(a)') trim(restartname) //'.nc' + + if (num_rest_files > 1) then + ! append i.th restart file name to rpointer + do i=1, num_rest_files-1 + if (i < 10) then + write(suffix,'("_",I1)') i + else + write(suffix,'("_",I2)') i + endif + write(nu,'(a)') trim(restartname) // trim(suffix) // '.nc' + enddo + endif + close(nu) write(glb%stdout,*) 'ocn restart pointer file written: ',trim(restartname) endif diff --git a/config_src/nuopc_driver/mom_cap.F90 b/config_src/nuopc_driver/mom_cap.F90 index 67bae67f74..ce11cfb3f9 100644 --- a/config_src/nuopc_driver/mom_cap.F90 +++ b/config_src/nuopc_driver/mom_cap.F90 @@ -15,7 +15,6 @@ module MOM_cap_mod use mpp_mod, only: stdlog, stdout, mpp_root_pe, mpp_clock_id use mpp_mod, only: mpp_clock_begin, mpp_clock_end, MPP_CLOCK_SYNC use mpp_mod, only: MPP_CLOCK_DETAILED, CLOCK_COMPONENT, MAXPES -use time_interp_external_mod, only: time_interp_external_init use time_manager_mod, only: set_calendar_type, time_type, increment_date use time_manager_mod, only: set_time, set_date, get_time, get_date, month_name use time_manager_mod, only: GREGORIAN, JULIAN, NOLEAP, THIRTY_DAY_MONTHS, NO_CALENDAR @@ -1559,6 +1558,7 @@ subroutine ModelAdvance(gcomp, rc) return endif write(writeunit,'(a)') trim(restartname)//'.nc' + if (num_rest_files > 1) then ! append i.th restart file name to rpointer do i=1, num_rest_files-1 @@ -1737,7 +1737,6 @@ subroutine ModelSetRunClock(gcomp, rc) line=__LINE__, file=__FILE__, rcToReturn=rc) return endif - ! not used in nems call NUOPC_CompAttributeGet(gcomp, name="restart_ymd", value=cvalue, & isPresent=isPresent, isSet=isSet, rc=rc) @@ -1811,6 +1810,8 @@ subroutine ocean_model_finalize(gcomp, rc) type(TIME_TYPE) :: Time type(ESMF_Clock) :: clock type(ESMF_Time) :: currTime + type(ESMF_Alarm), allocatable :: alarmList(:) + integer :: alarmCount character(len=64) :: timestamp logical :: write_restart character(len=*),parameter :: subname='(MOM_cap:ocean_model_finalize)' diff --git a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 index 2616d99e75..493762f4bc 100644 --- a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 +++ b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 @@ -39,6 +39,7 @@ module MOM_ocean_model_nuopc use MOM_time_manager, only : operator(+), operator(-), operator(*), operator(/) use MOM_time_manager, only : operator(/=), operator(<=), operator(>=) use MOM_time_manager, only : operator(<), real_to_time_type, time_type_to_real +use time_interp_external_mod,only : time_interp_external_init use MOM_tracer_flow_control, only : call_tracer_register, tracer_flow_control_init use MOM_tracer_flow_control, only : call_tracer_flux_init use MOM_unit_scaling, only : unit_scale_type @@ -267,6 +268,8 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i OS%is_ocean_pe = Ocean_sfc%is_ocean_pe if (.not.OS%is_ocean_pe) return + call time_interp_external_init + OS%Time = Time_in call initialize_MOM(OS%Time, Time_init, param_file, OS%dirs, OS%MOM_CSp, & OS%restart_CSp, Time_in, offline_tracer_mode=OS%offline_tracer_mode, & diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 957290a04f..101269c069 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -1661,24 +1661,25 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, call find_face_areas(Datu, Datv, G, GV, US, CS, MS, eta, 1+iev-ie) endif - !GOMP parallel default(shared) + !$OMP parallel default(shared) private(vel_prev, ioff, joff) if (CS%dynamic_psurf .or. .not.project_velocity) then if (use_BT_cont) then - !GOMP do + !$OMP do do j=jsv-1,jev+1 ; do I=isv-2,iev+1 uhbt(I,j) = find_uhbt(ubt(I,j), BTCL_u(I,j), US) + uhbt0(I,j) enddo ; enddo - !GOMP do + !$OMP end do nowait + !$OMP do do J=jsv-2,jev+1 ; do i=isv-1,iev+1 vhbt(i,J) = find_vhbt(vbt(i,J), BTCL_v(i,J), US) + vhbt0(i,J) enddo ; enddo - !GOMP do + !$OMP do do j=jsv-1,jev+1 ; do i=isv-1,iev+1 eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * CS%IareaT(i,j)) * & ((uhbt(I-1,j) - uhbt(I,j)) + (vhbt(i,J-1) - vhbt(i,J))) enddo ; enddo else - !GOMP do + !$OMP do do j=jsv-1,jev+1 ; do i=isv-1,iev+1 eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * CS%IareaT(i,j)) * & (((Datu(I-1,j)*ubt(I-1,j) + uhbt0(I-1,j)) - & @@ -1689,7 +1690,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, endif if (CS%dynamic_psurf) then - !GOMP do + !$OMP do do j=jsv-1,jev+1 ; do i=isv-1,iev+1 p_surf_dyn(i,j) = dyn_coef_eta(i,j) * (eta_pred(i,j) - eta(i,j)) enddo ; enddo @@ -1700,31 +1701,31 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! eta_PF_BT => eta_pred ; if (project_velocity) eta_PF_BT => eta if (find_etaav) then - !GOMP do + !$OMP do do j=js,je ; do i=is,ie eta_sum(i,j) = eta_sum(i,j) + wt_accel2(n) * eta_PF_BT(i,j) enddo ; enddo + !$OMP end do nowait endif if (interp_eta_PF) then wt_end = n*Instep ! This could be (n-0.5)*Instep. - !GOMP do + !$OMP do do j=jsv-1,jev+1 ; do i=isv-1,iev+1 eta_PF(i,j) = eta_PF_1(i,j) + wt_end*d_eta_PF(i,j) enddo ; enddo endif if (apply_OBC_flather .or. apply_OBC_open) then - !GOMP do + !$OMP do do j=jsv,jev ; do I=isv-2,iev+1 ubt_old(I,j) = ubt(I,j) enddo ; enddo - !GOMP do + !$OMP do do J=jsv-2,jev+1 ; do i=isv,iev vbt_old(i,J) = vbt(i,J) enddo ; enddo endif - !GOMP end parallel if (apply_OBCs) then if (MOD(n+G%first_direction,2)==1) then @@ -1734,7 +1735,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, endif if (CS%BT_OBC%apply_u_OBCs) then ! save the old value of ubt and uhbt - !GOMP parallel do default(shared) + !$OMP do do J=jsv-joff,jev+joff ; do i=isv-1,iev ubt_prev(i,J) = ubt(i,J) ; uhbt_prev(i,J) = uhbt(i,J) ubt_sum_prev(i,J) = ubt_sum(i,J) ; uhbt_sum_prev(i,J) = uhbt_sum(i,J) ; ubt_wtd_prev(i,J) = ubt_wtd(i,J) @@ -1742,7 +1743,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, endif if (CS%BT_OBC%apply_v_OBCs) then ! save the old value of vbt and vhbt - !GOMP parallel do default(shared) + !$OMP do do J=jsv-1,jev ; do i=isv-ioff,iev+ioff vbt_prev(i,J) = vbt(i,J) ; vhbt_prev(i,J) = vhbt(i,J) vbt_sum_prev(i,J) = vbt_sum(i,J) ; vhbt_sum_prev(i,J) = vhbt_sum(i,J) ; vbt_wtd_prev(i,J) = vbt_wtd(i,J) @@ -1750,10 +1751,9 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, endif endif - !GOMP parallel default(shared) private(vel_prev) if (MOD(n+G%first_direction,2)==1) then ! On odd-steps, update v first. - !GOMP do + !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv-1,iev+1 Cor_v(i,J) = -1.0*((amer(I-1,j) * ubt(I-1,j) + cmer(I,j+1) * ubt(I,j+1)) + & (bmer(I,j) * ubt(I,j) + dmer(I-1,j+1) * ubt(I-1,j+1))) - Cor_ref_v(i,J) @@ -1761,20 +1761,23 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, (eta_PF_BT(i,j+1)-eta_PF(i,j+1))*gtot_S(i,j+1)) * & dgeo_de * CS%IdyCv(i,J) enddo ; enddo + !$OMP end do nowait if (CS%dynamic_psurf) then - !GOMP do + !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv-1,iev+1 PFv(i,J) = PFv(i,J) + (p_surf_dyn(i,j) - p_surf_dyn(i,j+1)) * CS%IdyCv(i,J) enddo ; enddo + !$OMP end do nowait endif if (CS%BT_OBC%apply_v_OBCs) then ! zero out PF across boundary - !GOMP do + !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv-1,iev+1 ; if (OBC%segnum_v(i,J) /= OBC_NONE) then PFv(i,J) = 0.0 endif ; enddo ; enddo + !$OMP end do nowait endif - !GOMP do + !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv-1,iev+1 vel_prev = vbt(i,J) vbt(i,J) = bt_rem_v(i,J) * (vbt(i,J) + & @@ -1790,24 +1793,26 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, enddo ; enddo if (use_BT_cont) then - !GOMP do + !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv-1,iev+1 vhbt(i,J) = find_vhbt(vbt_trans(i,J), BTCL_v(i,J), US) + vhbt0(i,J) enddo ; enddo + !$OMP end do nowait else - !GOMP do + !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv-1,iev+1 vhbt(i,J) = Datv(i,J)*vbt_trans(i,J) + vhbt0(i,J) enddo ; enddo + !$OMP end do nowait endif if (CS%BT_OBC%apply_v_OBCs) then ! copy back the value for v-points on the boundary. - !GOMP do + !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv-1,iev+1 ; if (OBC%segnum_v(i,J) /= OBC_NONE) then vbt(i,J) = vbt_prev(i,J) ; vhbt(i,J) = vhbt_prev(i,J) endif ; enddo ; enddo endif ! Now update the zonal velocity. - !GOMP do + !$OMP do schedule(static) do j=jsv,jev ; do I=isv-1,iev Cor_u(I,j) = ((azon(I,j) * vbt(i+1,J) + czon(I,j) * vbt(i,J-1)) + & (bzon(I,j) * vbt(i,J) + dzon(I,j) * vbt(i+1,J-1))) - & @@ -1816,21 +1821,24 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, (eta_PF_BT(i+1,j)-eta_PF(i+1,j))*gtot_W(i+1,j)) * & dgeo_de * CS%IdxCu(I,j) enddo ; enddo + !$OMP end do nowait if (CS%dynamic_psurf) then - !GOMP do + !$OMP do schedule(static) do j=jsv,jev ; do I=isv-1,iev PFu(I,j) = PFu(I,j) + (p_surf_dyn(i,j) - p_surf_dyn(i+1,j)) * CS%IdxCu(I,j) enddo ; enddo + !$OMP end do nowait endif if (CS%BT_OBC%apply_u_OBCs) then ! zero out pressure force across boundary - !GOMP do + !$OMP do schedule(static) do j=jsv,jev ; do I=isv-1,iev ; if (OBC%segnum_u(I,j) /= OBC_NONE) then PFu(I,j) = 0.0 endif ; enddo ; enddo + !$OMP end do nowait endif - !GOMP do + !$OMP do schedule(static) do j=jsv,jev ; do I=isv-1,iev vel_prev = ubt(I,j) ubt(I,j) = bt_rem_u(I,j) * (ubt(I,j) + & @@ -1845,27 +1853,28 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * (Cor_u(I,j) + PFu(I,j)) endif enddo ; enddo + !$OMP end do nowait if (use_BT_cont) then - !GOMP do + !$OMP do schedule(static) do j=jsv,jev ; do I=isv-1,iev uhbt(I,j) = find_uhbt(ubt_trans(I,j), BTCL_u(I,j), US) + uhbt0(I,j) enddo ; enddo else - !GOMP do + !$OMP do schedule(static) do j=jsv,jev ; do I=isv-1,iev uhbt(I,j) = Datu(I,j)*ubt_trans(I,j) + uhbt0(I,j) enddo ; enddo endif if (CS%BT_OBC%apply_u_OBCs) then ! copy back the value for u-points on the boundary. - !GOMP do + !$OMP do schedule(static) do j=jsv,jev ; do I=isv-1,iev ; if (OBC%segnum_u(I,j) /= OBC_NONE) then ubt(I,j) = ubt_prev(I,j); uhbt(I,j) = uhbt_prev(I,j) endif ; enddo ; enddo endif else ! On even steps, update u first. - !GOMP do + !$OMP do schedule(static) do j=jsv-1,jev+1 ; do I=isv-1,iev Cor_u(I,j) = ((azon(I,j) * vbt(i+1,J) + czon(I,j) * vbt(i,J-1)) + & (bzon(I,j) * vbt(i,J) + dzon(I,j) * vbt(i+1,J-1))) - & @@ -1874,22 +1883,24 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, (eta_PF_BT(i+1,j)-eta_PF(i+1,j))*gtot_W(i+1,j)) * & dgeo_de * CS%IdxCu(I,j) enddo ; enddo + !$OMP end do nowait if (CS%dynamic_psurf) then - !GOMP do + !$OMP do schedule(static) do j=jsv-1,jev+1 ; do I=isv-1,iev PFu(I,j) = PFu(I,j) + (p_surf_dyn(i,j) - p_surf_dyn(i+1,j)) * CS%IdxCu(I,j) enddo ; enddo + !$OMP end do nowait endif if (CS%BT_OBC%apply_u_OBCs) then ! zero out pressure force across boundary - !GOMP do + !$OMP do schedule(static) do j=jsv,jev ; do I=isv-1,iev ; if (OBC%segnum_u(I,j) /= OBC_NONE) then PFu(I,j) = 0.0 endif ; enddo ; enddo endif - !GOMP do + !$OMP do schedule(static) do j=jsv-1,jev+1 ; do I=isv-1,iev vel_prev = ubt(I,j) ubt(I,j) = bt_rem_u(I,j) * (ubt(I,j) + & @@ -1906,18 +1917,20 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, enddo ; enddo if (use_BT_cont) then - !GOMP do + !$OMP do schedule(static) do j=jsv-1,jev+1 ; do I=isv-1,iev uhbt(I,j) = find_uhbt(ubt_trans(I,j), BTCL_u(I,j), US) + uhbt0(I,j) enddo ; enddo + !$OMP end do nowait else - !GOMP do + !$OMP do schedule(static) do j=jsv-1,jev+1 ; do I=isv-1,iev uhbt(I,j) = Datu(I,j)*ubt_trans(I,j) + uhbt0(I,j) enddo ; enddo + !$OMP end do nowait endif if (CS%BT_OBC%apply_u_OBCs) then ! copy back the value for u-points on the boundary. - !GOMP do + !$OMP do schedule(static) do j=jsv-1,jev+1 ; do I=isv-1,iev ; if (OBC%segnum_u(I,j) /= OBC_NONE) then ubt(I,j) = ubt_prev(I,j); uhbt(I,j) = uhbt_prev(I,j) endif ; enddo ; enddo @@ -1925,7 +1938,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! Now update the meridional velocity. if (CS%use_old_coriolis_bracket_bug) then - !GOMP do + !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv,iev Cor_v(i,J) = -1.0*((amer(I-1,j) * ubt(I-1,j) + bmer(I,j) * ubt(I,j)) + & (cmer(I,j+1) * ubt(I,j+1) + dmer(I-1,j+1) * ubt(I-1,j+1))) - Cor_ref_v(i,J) @@ -1933,8 +1946,9 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, (eta_PF_BT(i,j+1)-eta_PF(i,j+1))*gtot_S(i,j+1)) * & dgeo_de * CS%IdyCv(i,J) enddo ; enddo + !$OMP end do nowait else - !GOMP do + !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv,iev Cor_v(i,J) = -1.0*((amer(I-1,j) * ubt(I-1,j) + cmer(I,j+1) * ubt(I,j+1)) + & (bmer(I,j) * ubt(I,j) + dmer(I-1,j+1) * ubt(I-1,j+1))) - Cor_ref_v(i,J) @@ -1942,23 +1956,25 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, (eta_PF_BT(i,j+1)-eta_PF(i,j+1))*gtot_S(i,j+1)) * & dgeo_de * CS%IdyCv(i,J) enddo ; enddo + !$OMP end do nowait endif if (CS%dynamic_psurf) then - !GOMP do + !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv,iev PFv(i,J) = PFv(i,J) + (p_surf_dyn(i,j) - p_surf_dyn(i,j+1)) * CS%IdyCv(i,J) enddo ; enddo + !$OMP end do nowait endif if (CS%BT_OBC%apply_v_OBCs) then ! zero out PF across boundary - !GOMP do + !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv-1,iev+1 ; if (OBC%segnum_v(i,J) /= OBC_NONE) then PFv(i,J) = 0.0 endif ; enddo ; enddo endif - !GOMP do + !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv,iev vel_prev = vbt(i,J) vbt(i,J) = bt_rem_v(i,J) * (vbt(i,J) + & @@ -1973,65 +1989,69 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, v_accel_bt(I,j) = v_accel_bt(I,j) + wt_accel(n) * (Cor_v(i,J) + PFv(i,J)) endif enddo ; enddo + !$OMP end do nowait if (use_BT_cont) then - !GOMP do + !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv,iev vhbt(i,J) = find_vhbt(vbt_trans(i,J), BTCL_v(i,J), US) + vhbt0(i,J) enddo ; enddo else - !GOMP do + !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv,iev vhbt(i,J) = Datv(i,J)*vbt_trans(i,J) + vhbt0(i,J) enddo ; enddo endif if (CS%BT_OBC%apply_v_OBCs) then ! copy back the value for v-points on the boundary. - !GOMP do + !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv,iev ; if (OBC%segnum_v(i,J) /= OBC_NONE) then vbt(i,J) = vbt_prev(i,J); vhbt(i,J) = vhbt_prev(i,J) endif ; enddo ; enddo endif endif - !GOMP end parallel - !GOMP parallel default(shared) if (find_PF) then - !GOMP do + !$OMP do do j=js,je ; do I=is-1,ie PFu_bt_sum(I,j) = PFu_bt_sum(I,j) + wt_accel2(n) * PFu(I,j) enddo ; enddo - !GOMP do + !$OMP end do nowait + !$OMP do do J=js-1,je ; do i=is,ie PFv_bt_sum(i,J) = PFv_bt_sum(i,J) + wt_accel2(n) * PFv(i,J) enddo ; enddo + !$OMP end do nowait endif if (find_Cor) then - !GOMP do + !$OMP do do j=js,je ; do I=is-1,ie Coru_bt_sum(I,j) = Coru_bt_sum(I,j) + wt_accel2(n) * Cor_u(I,j) enddo ; enddo - !GOMP do + !$OMP end do nowait + !$OMP do do J=js-1,je ; do i=is,ie Corv_bt_sum(i,J) = Corv_bt_sum(i,J) + wt_accel2(n) * Cor_v(i,J) enddo ; enddo + !$OMP end do nowait endif - !GOMP do + !$OMP do do j=js,je ; do I=is-1,ie ubt_sum(I,j) = ubt_sum(I,j) + wt_trans(n) * ubt_trans(I,j) uhbt_sum(I,j) = uhbt_sum(I,j) + wt_trans(n) * uhbt(I,j) ubt_wtd(I,j) = ubt_wtd(I,j) + wt_vel(n) * ubt(I,j) enddo ; enddo - !GOMP do + !$OMP end do nowait + !$OMP do do J=js-1,je ; do i=is,ie vbt_sum(i,J) = vbt_sum(i,J) + wt_trans(n) * vbt_trans(i,J) vhbt_sum(i,J) = vhbt_sum(i,J) + wt_trans(n) * vhbt(i,J) vbt_wtd(i,J) = vbt_wtd(i,J) + wt_vel(n) * vbt(i,J) enddo ; enddo - !GOMP end parallel + !$OMP end do nowait if (apply_OBCs) then if (CS%BT_OBC%apply_u_OBCs) then ! copy back the value for u-points on the boundary. - !GOMP parallel do default(shared) + !$OMP do do j=js,je ; do I=is-1,ie if (OBC%segnum_u(I,j) /= OBC_NONE) then ubt_sum(I,j) = ubt_sum_prev(I,j) ; uhbt_sum(I,j) = uhbt_sum_prev(I,j) @@ -2041,7 +2061,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, endif if (CS%BT_OBC%apply_v_OBCs) then ! copy back the value for v-points on the boundary. - !GOMP parallel do default(shared) + !$OMP do do J=js-1,je ; do I=is,ie if (OBC%segnum_v(i,J) /= OBC_NONE) then vbt_sum(i,J) = vbt_sum_prev(i,J) ; vhbt_sum(i,J) = vhbt_sum_prev(i,J) @@ -2050,24 +2070,32 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, enddo ; enddo endif + !$OMP single call apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, & ubt_trans, vbt_trans, eta, ubt_old, vbt_old, CS%BT_OBC, & G, MS, US, iev-ie, dtbt, bebt, use_BT_cont, Datu, Datv, BTCL_u, BTCL_v, & uhbt0, vhbt0) - if (CS%BT_OBC%apply_u_OBCs) then ; do j=js,je ; do I=is-1,ie - if (OBC%segnum_u(I,j) /= OBC_NONE) then - ubt_sum(I,j) = ubt_sum(I,j) + wt_trans(n) * ubt_trans(I,j) - uhbt_sum(I,j) = uhbt_sum(I,j) + wt_trans(n) * uhbt(I,j) - ubt_wtd(I,j) = ubt_wtd(I,j) + wt_vel(n) * ubt(I,j) - endif - enddo ; enddo ; endif - if (CS%BT_OBC%apply_v_OBCs) then ; do J=js-1,je ; do i=is,ie - if (OBC%segnum_v(i,J) /= OBC_NONE) then - vbt_sum(i,J) = vbt_sum(i,J) + wt_trans(n) * vbt_trans(i,J) - vhbt_sum(i,J) = vhbt_sum(i,J) + wt_trans(n) * vhbt(i,J) - vbt_wtd(i,J) = vbt_wtd(i,J) + wt_vel(n) * vbt(i,J) - endif - enddo ; enddo ; endif + !$OMP end single + if (CS%BT_OBC%apply_u_OBCs) then + !$OMP do + do j=js,je ; do I=is-1,ie + if (OBC%segnum_u(I,j) /= OBC_NONE) then + ubt_sum(I,j) = ubt_sum(I,j) + wt_trans(n) * ubt_trans(I,j) + uhbt_sum(I,j) = uhbt_sum(I,j) + wt_trans(n) * uhbt(I,j) + ubt_wtd(I,j) = ubt_wtd(I,j) + wt_vel(n) * ubt(I,j) + endif + enddo ; enddo + endif + if (CS%BT_OBC%apply_v_OBCs) then + !$OMP do + do J=js-1,je ; do i=is,ie + if (OBC%segnum_v(i,J) /= OBC_NONE) then + vbt_sum(i,J) = vbt_sum(i,J) + wt_trans(n) * vbt_trans(i,J) + vhbt_sum(i,J) = vhbt_sum(i,J) + wt_trans(n) * vhbt(i,J) + vbt_wtd(i,J) = vbt_wtd(i,J) + wt_vel(n) * vbt(i,J) + endif + enddo ; enddo + endif endif if (CS%debug_bt) then @@ -2075,12 +2103,13 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, scale=US%s_to_T*US%L_to_m**2*GV%H_to_m) endif - !$OMP parallel do default(shared) + !$OMP do do j=jsv,jev ; do i=isv,iev eta(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * CS%IareaT(i,j)) * & ((uhbt(I-1,j) - uhbt(I,j)) + (vhbt(i,J-1) - vhbt(i,J))) eta_wtd(i,j) = eta_wtd(i,j) + eta(i,j) * wt_eta(n) enddo ; enddo + !$OMP end parallel if (do_hifreq_output) then time_step_end = time_bt_start + real_to_time(n*US%T_to_s*dtbt) diff --git a/src/framework/MOM_get_input.F90 b/src/framework/MOM_get_input.F90 index ad48086543..b6b5b89be9 100644 --- a/src/framework/MOM_get_input.F90 +++ b/src/framework/MOM_get_input.F90 @@ -21,7 +21,8 @@ module MOM_get_input character(len=240) :: & restart_input_dir = ' ',& !< The directory to read restart and input files. restart_output_dir = ' ',&!< The directory into which to write restart files. - output_directory = ' ', & !< The directory to use to write the model output. + output_directory = ' ' !< The directory to use to write the model output. + character(len=2048) :: & input_filename = ' ' !< A string that indicates the input files or how !! the run segment should be started. end type directories @@ -46,7 +47,8 @@ subroutine get_MOM_input(param_file, dirs, check_params, default_input_filename, parameter_filename(npf), & ! List of files containing parameters. output_directory, & ! Directory to use to write the model output. restart_input_dir, & ! Directory for reading restart and input files. - restart_output_dir, & ! Directory into which to write restart files. + restart_output_dir ! Directory into which to write restart files. + character(len=2048) :: & input_filename ! A string that indicates the input files or how ! the run segment should be started. character(len=240) :: output_dir diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index ca2e37afb9..20056c15ad 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -1057,7 +1057,9 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_ num_files = num_files+1 enddo + if (present(num_rest_files)) num_rest_files = num_files + end subroutine save_restart !> restore_state reads the model state from previously generated files. All diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index cf993b8aa8..953cc6d838 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -46,6 +46,8 @@ module MOM_hor_visc !! limited to guarantee stability. logical :: better_bound_Ah !< If true, use a more careful bounding of the !! biharmonic viscosity to guarantee stability. + real :: Re_Ah !! If nonzero, the biharmonic coefficient is scaled + !< so that the biharmonic Reynolds number is equal to this. real :: bound_coef !< The nondimensional coefficient of the ratio of !! the viscosity bounds to the theoretical maximum !! for stability without considering other terms [nondim]. @@ -105,11 +107,6 @@ module MOM_hor_visc !< The background biharmonic viscosity at h points [L4 T-1 ~> m4 s-1]. !! The actual viscosity may be the larger of this !! viscosity and the Smagorinsky and Leith viscosities. -! real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: Biharm5_const2_xx - !< A constant relating the biharmonic viscosity to the - !! square of the velocity shear [L4 T ~> m4 s]. This value is - !! set to be the magnitude of the Coriolis terms once the - !! velocity differences reach a value of order 1/2 MAXVEL. real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: reduction_xx !< The amount by which stresses through h points are reduced !! due to partial barriers [nondim]. @@ -117,8 +114,9 @@ module MOM_hor_visc Kh_Max_xx, & !< The maximum permitted Laplacian viscosity [L2 T-1 ~> m2 s-1]. Ah_Max_xx, & !< The maximum permitted biharmonic viscosity [L4 T-1 ~> m4 s-1]. n1n2_h, & !< Factor n1*n2 in the anisotropic direction tensor at h-points - n1n1_m_n2n2_h !< Factor n1**2-n2**2 in the anisotropic direction tensor at h-points - + n1n1_m_n2n2_h, & !< Factor n1**2-n2**2 in the anisotropic direction tensor at h-points + grid_sp_h2, & !< Harmonic mean of the squares of the grid [L2 ~> m2] + grid_sp_h3 !< Harmonic mean of the squares of the grid^(3/2) [L3 ~> m3] real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: Kh_bg_xy !< The background Laplacian viscosity at q points [L2 T-1 ~> m2 s-1]. !! The actual viscosity may be the larger of this @@ -127,11 +125,6 @@ module MOM_hor_visc !< The background biharmonic viscosity at q points [L4 T-1 ~> m4 s-1]. !! The actual viscosity may be the larger of this !! viscosity and the Smagorinsky and Leith viscosities. -! real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: Biharm5_const2_xy - !< A constant relating the biharmonic viscosity to the - !! square of the velocity shear [L4 T ~> m4 s]. This value is - !! set to be the magnitude of the Coriolis terms once the - !! velocity differences reach a value of order 1/2 MAXVEL. real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: reduction_xy !< The amount by which stresses through q points are reduced !! due to partial barriers [nondim]. @@ -162,27 +155,31 @@ module MOM_hor_visc ! parameters and metric terms. real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: & Laplac2_const_xx, & !< Laplacian metric-dependent constants [L2 ~> m2] - Biharm5_const_xx, & !< Biharmonic metric-dependent constants [L5 ~> m5] + Biharm6_const_xx, & !< Biharmonic metric-dependent constants [L6 ~> m6] Laplac3_const_xx, & !< Laplacian metric-dependent constants [L3 ~> m3] Biharm_const_xx, & !< Biharmonic metric-dependent constants [L4 ~> m4] - Biharm_const2_xx !< Biharmonic metric-dependent constants [T L4 ~> s m4] + Biharm_const2_xx, & !< Biharmonic metric-dependent constants [T L4 ~> s m4] + Re_Ah_const_xx !< Biharmonic metric-dependent constants [L3 ~> m3] real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: & Laplac2_const_xy, & !< Laplacian metric-dependent constants [L2 ~> m2] - Biharm5_const_xy, & !< Biharmonic metric-dependent constants [L5 ~> m5] + Biharm6_const_xy, & !< Biharmonic metric-dependent constants [L6 ~> m6] Laplac3_const_xy, & !< Laplacian metric-dependent constants [L3 ~> m3] Biharm_const_xy, & !< Biharmonic metric-dependent constants [L4 ~> m4] - Biharm_const2_xy !< Biharmonic metric-dependent constants [T L4 ~> s m4] + Biharm_const2_xy, & !< Biharmonic metric-dependent constants [T L4 ~> s m4] + Re_Ah_const_xy !< Biharmonic metric-dependent constants [L3 ~> m3] type(diag_ctrl), pointer :: diag => NULL() !< structure to regulate diagnostics !>@{ !! Diagnostic id + integer :: id_grid_Re_Ah = -1, id_grid_Re_Kh = -1 integer :: id_diffu = -1, id_diffv = -1 integer :: id_Ah_h = -1, id_Ah_q = -1 integer :: id_Kh_h = -1, id_Kh_q = -1 integer :: id_GME_coeff_h = -1, id_GME_coeff_q = -1 integer :: id_vort_xy_q = -1, id_div_xx_h = -1 + integer :: id_sh_xy_q = -1, id_sh_xx_h = -1 integer :: id_FrictWork = -1, id_FrictWorkIntz = -1 integer :: id_FrictWork_GME = -1 !>@} @@ -258,6 +255,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! Leith_Ah_h, & ! Leith bi-harmonic viscosity at h-points [L4 T-1 ~> m4 s-1] grad_vort_mag_h, & ! Magnitude of vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1] grad_vort_mag_h_2d, & ! Magnitude of 2d vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1] + Del2vort_h, & ! Laplacian of vorticity at h-points [L-2 T-1 ~> m-2 s-1] grad_div_mag_h, & ! Magnitude of divergence gradient at h-points [L-1 T-1 ~> m-1 s-1] dudx, dvdy, & ! components in the horizontal tension [T-1 ~> s-1] grad_vel_mag_h, & ! Magnitude of the velocity gradient tensor squared at h-points [T-2 ~> s-2] @@ -279,6 +277,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, Leith_Ah_q, & ! Leith bi-harmonic viscosity at q-points [L4 T-1 ~> m4 s-1] grad_vort_mag_q, & ! Magnitude of vorticity gradient at q-points [L-1 T-1 ~> m-1 s-1] grad_vort_mag_q_2d, & ! Magnitude of 2d vorticity gradient at q-points [L-1 T-1 ~> m-1 s-1] + Del2vort_q, & ! Laplacian of vorticity at q-points [L-2 T-1 ~> m-2 s-1] grad_div_mag_q, & ! Magnitude of divergence gradient at q-points [L-1 T-1 ~> m-1 s-1] grad_vel_mag_q, & ! Magnitude of the velocity gradient tensor squared at q-points [T-2 ~> s-2] hq, & ! harmonic mean of the harmonic means of the u- & v point thicknesses [H ~> m or kg m-2] @@ -290,6 +289,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, Ah_q, & ! biharmonic viscosity at corner points [L4 T-1 ~> m4 s-1] Kh_q, & ! Laplacian viscosity at corner points [L2 T-1 ~> m2 s-1] vort_xy_q, & ! vertical vorticity at corner points [T-1 ~> s-1] + sh_xy_q, & ! horizontal shearing strain at corner points [T-1 ~> s-1] GME_coeff_q, & !< GME coeff. at q-points [L2 T-1 ~> m2 s-1] max_diss_rate_q ! maximum possible energy dissipated by lateral friction [L2 T-3 ~> m2 s-3] @@ -302,9 +302,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, Kh_h, & ! Laplacian viscosity at thickness points [L2 T-1 ~> m2 s-1] max_diss_rate_h, & ! maximum possible energy dissipated by lateral friction [L2 T-3 ~> m2 s-3] FrictWork, & ! work done by MKE dissipation mechanisms [R L2 T-3 ~> W m-2] - FrictWork_GME, & ! work done by GME [R L2 T-3 ~> W m-2] - div_xx_h ! horizontal divergence [T-1 ~> s-1] + FrictWork_GME, & ! work done by GME [R L2 T-3 ~> W m-2] + div_xx_h, & ! horizontal divergence [T-1 ~> s-1] + sh_xx_h ! horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1] real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & + grid_Re_Kh, & !< Grid Reynolds number for Laplacian horizontal viscosity at h points [nondim] + grid_Re_Ah, & !< Grid Reynolds number for Biharmonic horizontal viscosity at h points [nondim] GME_coeff_h !< GME coeff. at h-points [L2 T-1 ~> m2 s-1] real :: Ah ! biharmonic viscosity [L4 T-1 ~> m4 s-1] real :: Kh ! Laplacian viscosity [L2 T-1 ~> m2 s-1] @@ -337,10 +340,18 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, real :: FWfrac ! Fraction of maximum theoretical energy transfer to use when scaling GME coefficient [nondim] real :: DY_dxBu ! Ratio of meridional over zonal grid spacing at vertices [nondim] real :: DX_dyBu ! Ratio of zonal over meridiononal grid spacing at vertices [nondim] + real :: DY_dxCv ! Ratio of meridional over zonal grid spacing at faces [nondim] + real :: DX_dyCu ! Ratio of zonal over meridional grid spacing at faces [nondim] real :: Sh_F_pow ! The ratio of shear over the absolute value of f raised to some power and rescaled [nondim] real :: backscat_subround ! The ratio of f over Shear_mag that is so small that the backscatter ! calculation gives the same value as if f were 0 [nondim]. real :: H0_GME ! Depth used to scale down GME coefficient in shallow areas [Z ~> m] + real :: KE ! Local kinetic energy [L2 T-2 ~> m2 s-2] + real, parameter :: KH_min = 1.E-30 ! This is the minimun horizontal Laplacian viscosity used to estimate the + ! grid Raynolds number [L2 T-1 ~> m2 s-1] + real, parameter :: AH_min = 1.E-30 ! This is the minimun horizontal Biharmonic viscosity used to estimate the + ! grid Raynolds number [L4 T-1 ~> m4 s-1] + logical :: rescale_Kh, legacy_bound logical :: find_FrictWork logical :: apply_OBC = .false. @@ -348,7 +359,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, logical :: use_MEKE_Au integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz integer :: i, j, k, n - real :: inv_PI3, inv_PI2, inv_PI5 + real :: inv_PI3, inv_PI2, inv_PI6 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -356,7 +367,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, h_neglect3 = h_neglect**3 inv_PI3 = 1.0/((4.0*atan(1.0))**3) inv_PI2 = 1.0/((4.0*atan(1.0))**2) - inv_PI5 = inv_PI3 * inv_PI2 + inv_PI6 = inv_PI3 * inv_PI3 Ah_h(:,:,:) = 0.0 Kh_h(:,:,:) = 0.0 @@ -467,11 +478,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !$OMP apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, & !$OMP use_MEKE_Ku, use_MEKE_Au, boundary_mask_h, boundary_mask_q, & !$OMP backscat_subround, GME_coeff_limiter, & - !$OMP h_neglect, h_neglect3, FWfrac, inv_PI3, inv_PI5, H0_GME, & + !$OMP h_neglect, h_neglect3, FWfrac, inv_PI3, inv_PI6, H0_GME, & !$OMP diffu, diffv, max_diss_rate_h, max_diss_rate_q, & !$OMP Kh_h, Kh_q, Ah_h, Ah_q, FrictWork, FrictWork_GME, & - !$OMP div_xx_h, vort_xy_q, GME_coeff_h, GME_coeff_q, & - !$OMP TD, KH_u_GME, KH_v_GME & + !$OMP div_xx_h, sh_xx_h, vort_xy_q, sh_xy_q, GME_coeff_h, GME_coeff_q, & + !$OMP TD, KH_u_GME, KH_v_GME, grid_Re_Kh, grid_Re_Ah & !$OMP ) & !$OMP private( & !$OMP i, j, k, n, & @@ -485,7 +496,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !$OMP grad_vel_mag_bt_h, grad_vel_mag_bt_q, grad_d2vel_mag_h, & !$OMP meke_res_fn, Shear_mag, vert_vort_mag, hrat_min, visc_bound_rem, & !$OMP Kh, Ah, AhSm, AhLth, local_strain, Sh_F_pow, & - !$OMP dDel2vdx, dDel2udy, & + !$OMP dDel2vdx, dDel2udy, DY_dxCv, DX_dyCu, Del2vort_q, Del2vort_h, KE, & !$OMP h2uq, h2vq, hu, hv, hq, FatH, RoScl, GME_coeff & !$OMP ) do k=1,nz @@ -504,7 +515,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo ! Components for the shearing strain - do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 + do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 dvdx(I,J) = CS%DY_dxBu(I,J)*(v(i+1,J,k)*G%IdyCv(i+1,J) - v(i,J,k)*G%IdyCv(i,J)) dudy(I,J) = CS%DX_dyBu(I,J)*(u(I,j+1,k)*G%IdxCu(I,j+1) - u(I,j,k)*G%IdxCu(I,j)) enddo ; enddo @@ -625,7 +636,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, elseif (OBC%segment(n)%direction == OBC_DIRECTION_S) then if ((J >= js-1) .and. (J <= je+1)) then do I = max(Isq-1,OBC%segment(n)%HI%isd), min(Ieq+1,OBC%segment(n)%HI%ied) - h_u(I,j) = h_u(i,j+1) + h_u(I,j) = h_u(I,j+1) enddo endif elseif (OBC%segment(n)%direction == OBC_DIRECTION_E) then @@ -681,35 +692,48 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif; endif endif - if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then + ! Vorticity + if (CS%no_slip) then + do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 + vort_xy(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx(I,J) - dudy(I,J) ) + enddo ; enddo + else + do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 + vort_xy(I,J) = G%mask2dBu(I,J) * ( dvdx(I,J) - dudy(I,J) ) + enddo ; enddo + endif - ! Vorticity - if (CS%no_slip) then - do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 - vort_xy(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx(I,J) - dudy(I,J) ) - enddo ; enddo - else - do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 - vort_xy(I,J) = G%mask2dBu(I,J) * ( dvdx(I,J) - dudy(I,J) ) - enddo ; enddo - endif + ! Divergence + do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 + div_xx(i,j) = dudx(i,j) + dvdy(i,j) + enddo ; enddo + + if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then ! Vorticity gradient - do J=js-2,Jeq+1 ; do i=is-1,Ieq+1 + do J=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+2 DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J) vort_xy_dx(i,J) = DY_dxBu * (vort_xy(I,J) * G%IdyCu(I,j) - vort_xy(I-1,J) * G%IdyCu(I-1,j)) enddo ; enddo - do j=js-1,Jeq+1 ; do I=is-2,Ieq+1 + do j=Jsq-1,Jeq+2 ; do I=Isq-1,Ieq+1 DX_dyBu = G%dxBu(I,J) * G%IdyBu(I,J) vort_xy_dy(I,j) = DX_dyBu * (vort_xy(I,J) * G%IdxCv(i,J) - vort_xy(I,J-1) * G%IdxCv(i,J-1)) enddo ; enddo + ! Laplacian of vorticity + do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 + DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J) + DX_dyBu = G%dxBu(I,J) * G%IdyBu(I,J) + + Del2vort_q(I,J) = DY_dxBu * (vort_xy_dx(i+1,J) * G%IdyCv(i+1,J) - vort_xy_dx(i,J) * G%IdyCv(i,J)) + & + DX_dyBu * (vort_xy_dy(I,j+1) * G%IdyCu(I,j+1) - vort_xy_dy(I,j) * G%IdyCu(I,j)) + enddo ; enddo + do J=Jsq,Jeq+1 ; do I=Isq,Ieq+1 + Del2vort_h(i,j) = 0.25*(Del2vort_q(I,J) + Del2vort_q(I-1,J) + Del2vort_q(I,J-1) + Del2vort_q(I-1,J-1)) + enddo ; enddo + if (CS%modified_Leith) then - ! Divergence - do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 - div_xx(i,j) = dudx(i,j) + dvdy(i,j) - enddo ; enddo ! Divergence gradient do j=Jsq-1,Jeq+2 ; do I=Isq-1,Ieq+1 @@ -839,7 +863,14 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif if ((CS%id_Kh_h>0) .or. find_FrictWork .or. CS%debug) Kh_h(i,j,k) = Kh + + if (CS%id_grid_Re_Kh>0) then + KE = 0.125*((u(I,j,k)+u(I-1,j,k))**2 + (v(i,J,k)+v(i,J-1,k))**2) + grid_Re_Kh(i,j,k) = (sqrt(KE) * sqrt(CS%grid_sp_h2(i,j)))/MAX(Kh,KH_min) + endif + if (CS%id_div_xx_h>0) div_xx_h(i,j,k) = div_xx(i,j) + if (CS%id_sh_xx_h>0) sh_xx_h(i,j,k) = sh_xx(i,j) str_xx(i,j) = -Kh * sh_xx(i,j) else ! not Laplacian @@ -866,7 +897,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, AhSm = CS%Biharm_const_xx(i,j) * Shear_mag endif endif - if (CS%Leith_Ah) AhLth = CS%Biharm5_const_xx(i,j) * vert_vort_mag * inv_PI5 + if (CS%Leith_Ah) AhLth = CS%Biharm6_const_xx(i,j) * abs(Del2vort_h(i,j)) * inv_PI6 Ah = MAX(MAX(CS%Ah_bg_xx(i,j), AhSm), AhLth) if (CS%bound_Ah .and. .not.CS%better_bound_Ah) & Ah = MIN(Ah, CS%Ah_Max_xx(i,j)) @@ -876,12 +907,22 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (use_MEKE_Au) Ah = Ah + MEKE%Au(i,j) ! *Add* the MEKE contribution + if (CS%Re_Ah > 0.0) then + KE = 0.125*((u(I,j,k)+u(I-1,j,k))**2 + (v(i,J,k)+v(i,J-1,k))**2) + Ah = sqrt(KE) * CS%Re_Ah_const_xx(i,j) + endif + if (CS%better_bound_Ah) then Ah = MIN(Ah, visc_bound_rem*hrat_min*CS%Ah_Max_xx(i,j)) endif if ((CS%id_Ah_h>0) .or. find_FrictWork .or. CS%debug) Ah_h(i,j,k) = Ah + if (CS%id_grid_Re_Ah>0) then + KE = 0.125*((u(I,j,k)+u(I-1,j,k))**2 + (v(i,J,k)+v(i,J-1,k))**2) + grid_Re_Ah(i,j,k) = (sqrt(KE) * CS%grid_sp_h3(i,j))/MAX(Ah,AH_min) + endif + str_xx(i,j) = str_xx(i,j) + Ah * & (CS%DY_dxT(i,j) * (G%IdyCu(I,j)*Del2u(I,j) - G%IdyCu(I-1,j)*Del2u(I-1,j)) - & CS%DX_dyT(i,j) * (G%IdxCv(i,J)*Del2v(i,J) - G%IdxCv(i,J-1)*Del2v(i,J-1))) @@ -1010,6 +1051,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%id_Kh_q>0 .or. CS%debug) Kh_q(I,J,k) = Kh if (CS%id_vort_xy_q>0) vort_xy_q(I,J,k) = vort_xy(I,J) + if (CS%id_sh_xy_q>0) sh_xy_q(I,J,k) = sh_xy(I,J) str_xy(I,J) = -Kh * sh_xy(I,J) else ! not Laplacian @@ -1036,7 +1078,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, AhSm = CS%Biharm_const_xy(I,J) * Shear_mag endif endif - if (CS%Leith_Ah) AhLth = CS%Biharm5_const_xy(I,J) * vert_vort_mag * inv_PI5 + if (CS%Leith_Ah) AhLth = CS%Biharm6_const_xy(I,J) * abs(Del2vort_q(I,J)) * inv_PI6 Ah = MAX(MAX(CS%Ah_bg_xy(I,J), AhSm), AhLth) if (CS%bound_Ah .and. .not.CS%better_bound_Ah) & Ah = MIN(Ah, CS%Ah_Max_xy(I,J)) @@ -1045,8 +1087,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif ! Smagorinsky_Ah or Leith_Ah if (use_MEKE_Au) then ! *Add* the MEKE contribution - Ah = Ah + 0.25*( (MEKE%Au(I,J) + MEKE%Au(I+1,J+1)) + & - (MEKE%Au(I+1,J) + MEKE%Au(I,J+1)) ) + Ah = Ah + 0.25*( (MEKE%Au(i,j) + MEKE%Au(i+1,j+1)) + & + (MEKE%Au(i+1,j) + MEKE%Au(i,j+1)) ) + endif + + if (CS%Re_Ah > 0.0) then + KE = 0.125*((u(I,j,k)+u(I,j+1,k))**2 + (v(i,J,k)+v(i+1,J,k))**2) + Ah = sqrt(KE) * CS%Re_Ah_const_xy(i,j) endif if (CS%better_bound_Ah) then @@ -1147,7 +1194,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, CS%dy2h(i+1,j)*str_xx(i+1,j)) + & G%IdxCu(I,j)*(CS%dx2q(I,J-1)*str_xy(I,J-1) - & CS%dx2q(I,J) *str_xy(I,J))) * & - G%IareaCu(I,j)) / (h_u(i,j) + h_neglect) + G%IareaCu(I,j)) / (h_u(I,j) + h_neglect) enddo ; enddo if (apply_OBC) then @@ -1282,10 +1329,14 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%id_FrictWork>0) call post_data(CS%id_FrictWork, FrictWork, CS%diag) if (CS%id_FrictWork_GME>0) call post_data(CS%id_FrictWork_GME, FrictWork_GME, CS%diag) if (CS%id_Ah_h>0) call post_data(CS%id_Ah_h, Ah_h, CS%diag) + if (CS%id_grid_Re_Ah>0) call post_data(CS%id_grid_Re_Ah, grid_Re_Ah, CS%diag) if (CS%id_div_xx_h>0) call post_data(CS%id_div_xx_h, div_xx_h, CS%diag) if (CS%id_vort_xy_q>0) call post_data(CS%id_vort_xy_q, vort_xy_q, CS%diag) + if (CS%id_sh_xx_h>0) call post_data(CS%id_sh_xx_h, sh_xx_h, CS%diag) + if (CS%id_sh_xy_q>0) call post_data(CS%id_sh_xy_q, sh_xy_q, CS%diag) if (CS%id_Ah_q>0) call post_data(CS%id_Ah_q, Ah_q, CS%diag) if (CS%id_Kh_h>0) call post_data(CS%id_Kh_h, Kh_h, CS%diag) + if (CS%id_grid_Re_Kh>0) call post_data(CS%id_grid_Re_Kh, grid_Re_Kh, CS%diag) if (CS%id_Kh_q>0) call post_data(CS%id_Kh_q, Kh_q, CS%diag) if (CS%id_GME_coeff_h > 0) call post_data(CS%id_GME_coeff_h, GME_coeff_h, CS%diag) if (CS%id_GME_coeff_q > 0) call post_data(CS%id_GME_coeff_q, GME_coeff_q, CS%diag) @@ -1365,7 +1416,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) logical :: split ! If true, use the split time stepping scheme. ! If false and USE_GME = True, issue a FATAL error. logical :: default_2018_answers - character(len=64) :: inputdir, filename real :: deg2rad ! Converts degrees to radians real :: slat_fn ! sin(lat)**Kh_pwr_of_sine @@ -1374,28 +1424,22 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB integer :: i, j - ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mdl = "MOM_hor_visc" ! module name - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB - if (associated(CS)) then call MOM_error(WARNING, "hor_visc_init called with an associated "// & "control structure.") return endif allocate(CS) - CS%diag => diag - ! Read parameters and write them to the model log. call log_version(param_file, mdl, version, "") - ! It is not clear whether these initialization lines are needed for the ! cases where the corresponding parameters are not read. CS%bound_Kh = .false. ; CS%better_bound_Kh = .false. ; CS%Smagorinsky_Kh = .false. ; CS%Leith_Kh = .false. @@ -1405,13 +1449,10 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) CS%Modified_Leith = .false. CS%anisotropic = .false. CS%dynamic_aniso = .false. - Kh = 0.0 ; Ah = 0.0 - ! If GET_ALL_PARAMS is true, all parameters are read in all cases to enable ! parameter spelling checks. call get_param(param_file, mdl, "GET_ALL_PARAMS", get_all, default=.false.) - 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.) @@ -1419,9 +1460,7 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) "If true, use the order of arithmetic and expressions that recover the "//& "answers from the end of 2018. Otherwise, use updated and more robust "//& "forms of the same expressions.", default=default_2018_answers) - call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false.) - call get_param(param_file, mdl, "LAPLACIAN", CS%Laplacian, & "If true, use a Laplacian horizontal viscosity.", & default=.false.) @@ -1447,7 +1486,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) "The power used to raise SIN(LAT) when using a latitudinally "//& "dependent background viscosity.", & units = "nondim", default=4.0) - call get_param(param_file, mdl, "SMAGORINSKY_KH", CS%Smagorinsky_Kh, & "If true, use a Smagorinsky nonlinear eddy viscosity.", & default=.false.) @@ -1456,11 +1494,9 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) "The nondimensional Laplacian Smagorinsky constant, "//& "often 0.15.", units="nondim", default=0.0, & fail_if_missing = CS%Smagorinsky_Kh) - call get_param(param_file, mdl, "LEITH_KH", CS%Leith_Kh, & "If true, use a Leith nonlinear eddy viscosity.", & default=.false.) - call get_param(param_file, mdl, "MODIFIED_LEITH", CS%Modified_Leith, & "If true, add a term to Leith viscosity which is "//& "proportional to the gradient of divergence.", & @@ -1468,7 +1504,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) call get_param(param_file, mdl, "RES_SCALE_MEKE_VISC", CS%res_scale_MEKE, & "If true, the viscosity contribution from MEKE is scaled by "//& "the resolution function.", default=.false.) - if (CS%Leith_Kh .or. get_all) then call get_param(param_file, mdl, "LEITH_LAP_CONST", Leith_Lap_const, & "The nondimensional Laplacian Leith constant, "//& @@ -1527,7 +1562,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) "to the spherical coordinates.", units = "nondim", fail_if_missing=.true.) end select endif - call get_param(param_file, mdl, "BIHARMONIC", CS%biharmonic, & "If true, use a biharmonic horizontal viscosity. "//& "BIHARMONIC may be used with LAPLACIAN.", & @@ -1554,7 +1588,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) call get_param(param_file, mdl, "LEITH_AH", CS%Leith_Ah, & "If true, use a biharmonic Leith nonlinear eddy "//& "viscosity.", default=.false.) - call get_param(param_file, mdl, "BOUND_AH", CS%bound_Ah, & "If true, the biharmonic coefficient is locally limited "//& "to be stable.", default=.true.) @@ -1562,13 +1595,16 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) "If true, the biharmonic coefficient is locally limited "//& "to be stable with a better bounding than just BOUND_AH.", & default=CS%bound_Ah) + call get_param(param_file, mdl, "RE_AH", CS%Re_Ah, & + "If nonzero, the biharmonic coefficient is scaled "//& + "so that the biharmonic Reynolds number is equal to this.", & + units="nondim", default=0.0) if (CS%Smagorinsky_Ah .or. get_all) then call get_param(param_file, mdl, "SMAG_BI_CONST",Smag_bi_const, & "The nondimensional biharmonic Smagorinsky constant, "//& "typically 0.015 - 0.06.", units="nondim", default=0.0, & fail_if_missing = CS%Smagorinsky_Ah) - call get_param(param_file, mdl, "BOUND_CORIOLIS", bound_Cor_def, default=.false.) call get_param(param_file, mdl, "BOUND_CORIOLIS_BIHARM", CS%bound_Coriolis, & "If true use a viscosity that increases with the square "//& @@ -1587,29 +1623,24 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) units="m s-1", default=maxvel, scale=US%m_s_to_L_T) endif endif - if (CS%Leith_Ah .or. get_all) & call get_param(param_file, mdl, "LEITH_BI_CONST", Leith_bi_const, & "The nondimensional biharmonic Leith constant, "//& "typical values are thus far undetermined.", units="nondim", default=0.0, & fail_if_missing = CS%Leith_Ah) - endif - call get_param(param_file, mdl, "USE_LAND_MASK_FOR_HVISC", CS%use_land_mask, & "If true, use Use the land mask for the computation of thicknesses "//& "at velocity locations. This eliminates the dependence on arbitrary "//& "values over land or outside of the domain. Default is False in order to "//& "maintain answers with legacy experiments but should be changed to True "//& "for new experiments.", default=.false.) - if (CS%better_bound_Ah .or. CS%better_bound_Kh .or. get_all) & call get_param(param_file, mdl, "HORVISC_BOUND_COEF", CS%bound_coef, & "The nondimensional coefficient of the ratio of the "//& "viscosity bounds to the theoretical maximum for "//& "stability without considering other terms.", units="nondim", & default=0.8) - call get_param(param_file, mdl, "NOSLIP", CS%no_slip, & "If true, no slip boundary conditions are used; otherwise "//& "free slip boundary conditions are assumed. The "//& @@ -1617,7 +1648,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) "cleaner than the no slip BCs. The use of free slip BCs "//& "is strongly encouraged, and no slip BCs are not used with "//& "the biharmonic viscosity.", default=.false.) - call get_param(param_file, mdl, "USE_KH_BG_2D", CS%use_Kh_bg_2d, & "If true, read a file containing 2-d background harmonic "//& "viscosities. The final viscosity is the maximum of the other "//& @@ -1631,38 +1661,30 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) call get_param(param_file, mdl, "USE_GME", CS%use_GME, & "If true, use the GM+E backscatter scheme in association \n"//& "with the Gent and McWilliams parameterization.", default=.false.) - if (CS%use_GME) then call get_param(param_file, mdl, "SPLIT", split, & "Use the split time stepping if true.", default=.true., & do_not_log=.true.) if (.not. split) call MOM_error(FATAL,"ERROR: Currently, USE_GME = True "// & "cannot be used with SPLIT=False.") - call get_param(param_file, mdl, "GME_H0", CS%GME_h0, & "The strength of GME tapers quadratically to zero when the bathymetric "//& "depth is shallower than GME_H0.", units="m", scale=US%m_to_Z, & default=1000.0) - call get_param(param_file, mdl, "GME_EFFICIENCY", CS%GME_efficiency, & "The nondimensional prefactor multiplying the GME coefficient.", & units="nondim", default=1.0) - call get_param(param_file, mdl, "GME_LIMITER", CS%GME_limiter, & "The absolute maximum value the GME coefficient is allowed to take.", & units="m2 s-1", scale=US%m_to_L**2*US%T_to_s, default=1.0e7) - endif - if (CS%bound_Kh .or. CS%bound_Ah .or. CS%better_bound_Kh .or. CS%better_bound_Ah) & call get_param(param_file, mdl, "DT", dt, & "The (baroclinic) dynamics time step.", units="s", scale=US%s_to_T, & fail_if_missing=.true.) - if (CS%no_slip .and. CS%biharmonic) & call MOM_error(FATAL,"ERROR: NOSLIP and BIHARMONIC cannot be defined "// & "at the same time in MOM.") - if (.not.(CS%Laplacian .or. CS%biharmonic)) then ! Only issue inviscid warning if not in single column mode (usually 2x2 domain) if ( max(G%domain%niglobal, G%domain%njglobal)>2 ) call MOM_error(WARNING, & @@ -1670,9 +1692,7 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) "LAPLACIAN or BIHARMONIC viscosity.") return ! We are not using either Laplacian or Bi-harmonic lateral viscosity endif - deg2rad = atan(1.0) / 45. - ALLOC_(CS%dx2h(isd:ied,jsd:jed)) ; CS%dx2h(:,:) = 0.0 ALLOC_(CS%dy2h(isd:ied,jsd:jed)) ; CS%dy2h(:,:) = 0.0 ALLOC_(CS%dx2q(IsdB:IedB,JsdB:JedB)) ; CS%dx2q(:,:) = 0.0 @@ -1681,8 +1701,8 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) ALLOC_(CS%dy_dxT(isd:ied,jsd:jed)) ; CS%dy_dxT(:,:) = 0.0 ALLOC_(CS%dx_dyBu(IsdB:IedB,JsdB:JedB)) ; CS%dx_dyBu(:,:) = 0.0 ALLOC_(CS%dy_dxBu(IsdB:IedB,JsdB:JedB)) ; CS%dy_dxBu(:,:) = 0.0 - if (CS%Laplacian) then + ALLOC_(CS%grid_sp_h2(isd:ied,jsd:jed)) ; CS%grid_sp_h2(:,:) = 0.0 ALLOC_(CS%Kh_bg_xx(isd:ied,jsd:jed)) ; CS%Kh_bg_xx(:,:) = 0.0 ALLOC_(CS%Kh_bg_xy(IsdB:IedB,JsdB:JedB)) ; CS%Kh_bg_xy(:,:) = 0.0 if (CS%bound_Kh .or. CS%better_bound_Kh) then @@ -1700,7 +1720,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) endif ALLOC_(CS%reduction_xx(isd:ied,jsd:jed)) ; CS%reduction_xx(:,:) = 0.0 ALLOC_(CS%reduction_xy(IsdB:IedB,JsdB:JedB)) ; CS%reduction_xy(:,:) = 0.0 - if (CS%anisotropic) then ALLOC_(CS%n1n2_h(isd:ied,jsd:jed)) ; CS%n1n2_h(:,:) = 0.0 ALLOC_(CS%n1n1_m_n2n2_h(isd:ied,jsd:jed)) ; CS%n1n1_m_n2n2_h(:,:) = 0.0 @@ -1718,7 +1737,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) "Runtime parameter ANISOTROPIC_MODE is out of range.") end select endif - if (CS%use_Kh_bg_2d) then ALLOC_(CS%Kh_bg_2d(isd:ied,jsd:jed)) ; CS%Kh_bg_2d(:,:) = 0.0 call get_param(param_file, mdl, "KH_BG_2D_FILENAME", filename, & @@ -1730,15 +1748,14 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) G%domain, timelevel=1, scale=US%m_to_L**2*US%T_to_s) call pass_var(CS%Kh_bg_2d, G%domain) endif - if (CS%biharmonic) then ALLOC_(CS%Idx2dyCu(IsdB:IedB,jsd:jed)) ; CS%Idx2dyCu(:,:) = 0.0 ALLOC_(CS%Idx2dyCv(isd:ied,JsdB:JedB)) ; CS%Idx2dyCv(:,:) = 0.0 ALLOC_(CS%Idxdy2u(IsdB:IedB,jsd:jed)) ; CS%Idxdy2u(:,:) = 0.0 ALLOC_(CS%Idxdy2v(isd:ied,JsdB:JedB)) ; CS%Idxdy2v(:,:) = 0.0 - ALLOC_(CS%Ah_bg_xx(isd:ied,jsd:jed)) ; CS%Ah_bg_xx(:,:) = 0.0 ALLOC_(CS%Ah_bg_xy(IsdB:IedB,JsdB:JedB)) ; CS%Ah_bg_xy(:,:) = 0.0 + ALLOC_(CS%grid_sp_h3(isd:ied,jsd:jed)) ; CS%grid_sp_h3(:,:) = 0.0 if (CS%bound_Ah .or. CS%better_bound_Ah) then ALLOC_(CS%Ah_Max_xx(isd:ied,jsd:jed)) ; CS%Ah_Max_xx(:,:) = 0.0 ALLOC_(CS%Ah_Max_xy(IsdB:IedB,JsdB:JedB)) ; CS%Ah_Max_xy(:,:) = 0.0 @@ -1752,11 +1769,14 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) endif endif if (CS%Leith_Ah) then - ALLOC_(CS%biharm5_const_xx(isd:ied,jsd:jed)) ; CS%biharm5_const_xx(:,:) = 0.0 - ALLOC_(CS%biharm5_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%biharm5_const_xy(:,:) = 0.0 + ALLOC_(CS%biharm6_const_xx(isd:ied,jsd:jed)) ; CS%biharm6_const_xx(:,:) = 0.0 + ALLOC_(CS%biharm6_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%biharm6_const_xy(:,:) = 0.0 + endif + if (CS%Re_Ah > 0.0) then + ALLOC_(CS%Re_Ah_const_xx(isd:ied,jsd:jed)); CS%Re_Ah_const_xx(:,:) = 0.0 + ALLOC_(CS%Re_Ah_const_xy(IsdB:IedB,JsdB:JedB)); CS%Re_Ah_const_xy(:,:) = 0.0 endif endif - do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 CS%dx2q(I,J) = G%dxBu(I,J)*G%dxBu(I,J) ; CS%dy2q(I,J) = G%dyBu(I,J)*G%dyBu(I,J) CS%DX_dyBu(I,J) = G%dxBu(I,J)*G%IdyBu(I,J) ; CS%DY_dxBu(I,J) = G%dyBu(I,J)*G%IdxBu(I,J) @@ -1765,7 +1785,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) CS%dx2h(i,j) = G%dxT(i,j)*G%dxT(i,j) ; CS%dy2h(i,j) = G%dyT(i,j)*G%dyT(i,j) CS%DX_dyT(i,j) = G%dxT(i,j)*G%IdyT(i,j) ; CS%DY_dxT(i,j) = G%dyT(i,j)*G%IdxT(i,j) enddo ; enddo - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 CS%reduction_xx(i,j) = 1.0 if ((G%dy_Cu(I,j) > 0.0) .and. (G%dy_Cu(I,j) < G%dyCu(I,j)) .and. & @@ -1781,7 +1800,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) (G%dx_Cv(i,J-1) < G%dxCv(i,J-1) * CS%reduction_xx(i,j))) & CS%reduction_xx(i,j) = G%dx_Cv(i,J-1) / (G%dxCv(i,J-1)) enddo ; enddo - do J=js-1,Jeq ; do I=is-1,Ieq CS%reduction_xy(I,J) = 1.0 if ((G%dy_Cu(I,j) > 0.0) .and. (G%dy_Cu(I,j) < G%dyCu(I,j)) .and. & @@ -1797,38 +1815,33 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) (G%dx_Cv(i+1,J) < G%dxCv(i+1,J) * CS%reduction_xy(I,J))) & CS%reduction_xy(I,J) = G%dx_Cv(i+1,J) / (G%dxCv(i+1,J)) enddo ; enddo - if (CS%Laplacian) then ! The 0.3 below was 0.4 in MOM1.10. The change in hq requires ! this to be less than 1/3, rather than 1/2 as before. if (CS%bound_Kh .or. CS%bound_Ah) Kh_Limit = 0.3 / (dt*4.0) - ! Calculate and store the background viscosity at h-points do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 ! Static factors in the Smagorinsky and Leith schemes grid_sp_h2 = (2.0*CS%dx2h(i,j)*CS%dy2h(i,j)) / (CS%dx2h(i,j) + CS%dy2h(i,j)) + CS%grid_sp_h2(i,j) = grid_sp_h2 grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2) if (CS%Smagorinsky_Kh) CS%Laplac2_const_xx(i,j) = Smag_Lap_const * grid_sp_h2 if (CS%Leith_Kh) CS%Laplac3_const_xx(i,j) = Leith_Lap_const * grid_sp_h3 ! Maximum of constant background and MICOM viscosity CS%Kh_bg_xx(i,j) = MAX(Kh, Kh_vel_scale * sqrt(grid_sp_h2)) - ! Use the larger of the above and values read from a file if (CS%use_Kh_bg_2d) CS%Kh_bg_xx(i,j) = MAX(CS%Kh_bg_2d(i,j), CS%Kh_bg_xx(i,j)) - ! Use the larger of the above and a function of sin(latitude) if (Kh_sin_lat>0.) then slat_fn = abs( sin( deg2rad * G%geoLatT(i,j) ) ) ** Kh_pwr_of_sine CS%Kh_bg_xx(i,j) = MAX(Kh_sin_lat * slat_fn, CS%Kh_bg_xx(i,j)) endif - if (CS%bound_Kh .and. .not.CS%better_bound_Kh) then ! Limit the background viscosity to be numerically stable CS%Kh_Max_xx(i,j) = Kh_Limit * grid_sp_h2 CS%Kh_bg_xx(i,j) = MIN(CS%Kh_bg_xx(i,j), CS%Kh_Max_xx(i,j)) endif enddo ; enddo - ! Calculate and store the background viscosity at q-points do J=js-1,Jeq ; do I=is-1,Ieq ! Static factors in the Smagorinsky and Leith schemes @@ -1838,7 +1851,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) if (CS%Leith_Kh) CS%Laplac3_const_xy(I,J) = Leith_Lap_const * grid_sp_q3 ! Maximum of constant background and MICOM viscosity CS%Kh_bg_xy(I,J) = MAX(Kh, Kh_vel_scale * sqrt(grid_sp_q2)) - ! Use the larger of the above and values read from a file if (CS%use_Kh_bg_2d) then ; if (CS%Kh_bg_2d_bug) then ! This option is unambiguously wrong, and should be obsoleted as soon as possible. @@ -1848,13 +1860,11 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) 0.25*((CS%Kh_bg_2d(i,j) + CS%Kh_bg_2d(i+1,j+1)) + & (CS%Kh_bg_2d(i+1,j) + CS%Kh_bg_2d(i,j+1))) ) endif ; endif - ! Use the larger of the above and a function of sin(latitude) if (Kh_sin_lat>0.) then slat_fn = abs( sin( deg2rad * G%geoLatBu(I,J) ) ) ** Kh_pwr_of_sine CS%Kh_bg_xy(I,J) = MAX(Kh_sin_lat * slat_fn, CS%Kh_bg_xy(I,J)) endif - if (CS%bound_Kh .and. .not.CS%better_bound_Kh) then ! Limit the background viscosity to be numerically stable CS%Kh_Max_xy(I,J) = Kh_Limit * grid_sp_q2 @@ -1862,9 +1872,7 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) endif enddo ; enddo endif - if (CS%biharmonic) then - do j=js-1,Jeq+1 ; do I=Isq-1,Ieq+1 CS%Idx2dyCu(I,j) = (G%IdxCu(I,j)*G%IdxCu(I,j)) * G%IdyCu(I,j) CS%Idxdy2u(I,j) = G%IdxCu(I,j) * (G%IdyCu(I,j)*G%IdyCu(I,j)) @@ -1873,7 +1881,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) CS%Idx2dyCv(i,J) = (G%IdxCv(i,J)*G%IdxCv(i,J)) * G%IdyCv(i,J) CS%Idxdy2v(i,J) = G%IdxCv(i,J) * (G%IdyCv(i,J)*G%IdyCv(i,J)) enddo ; enddo - CS%Ah_bg_xy(:,:) = 0.0 ! The 0.3 below was 0.4 in MOM1.10. The change in hq requires ! this to be less than 1/3, rather than 1/2 as before. @@ -1883,7 +1890,7 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 grid_sp_h2 = (2.0*CS%dx2h(i,j)*CS%dy2h(i,j)) / (CS%dx2h(i,j)+CS%dy2h(i,j)) grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2) - + CS%grid_sp_h3(i,j) = grid_sp_h3 if (CS%Smagorinsky_Ah) then CS%Biharm_const_xx(i,j) = Smag_bi_const * (grid_sp_h2 * grid_sp_h2) if (CS%bound_Coriolis) then @@ -1894,9 +1901,10 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) endif endif if (CS%Leith_Ah) then - CS%biharm5_const_xx(i,j) = Leith_bi_const * (grid_sp_h3 * grid_sp_h2) + CS%biharm6_const_xx(i,j) = Leith_bi_const * (grid_sp_h3 * grid_sp_h3) endif CS%Ah_bg_xx(i,j) = MAX(Ah, Ah_vel_scale * grid_sp_h2 * sqrt(grid_sp_h2)) + if (CS%Re_Ah > 0.0) CS%Re_Ah_const_xx(i,j) = grid_sp_h3 / CS%Re_Ah if (Ah_time_scale > 0.) CS%Ah_bg_xx(i,j) = & MAX(CS%Ah_bg_xx(i,j), (grid_sp_h2 * grid_sp_h2) / Ah_time_scale) if (CS%bound_Ah .and. .not.CS%better_bound_Ah) then @@ -1907,7 +1915,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) do J=js-1,Jeq ; do I=is-1,Ieq grid_sp_q2 = (2.0*CS%dx2q(I,J)*CS%dy2q(I,J)) / (CS%dx2q(I,J)+CS%dy2q(I,J)) grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2) - if (CS%Smagorinsky_Ah) then CS%Biharm_const_xy(I,J) = Smag_bi_const * (grid_sp_q2 * grid_sp_q2) if (CS%bound_Coriolis) then @@ -1916,10 +1923,10 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) endif endif if (CS%Leith_Ah) then - CS%biharm5_const_xy(i,j) = Leith_bi_const * (grid_sp_q3 * grid_sp_q2) + CS%biharm6_const_xy(I,J) = Leith_bi_const * (grid_sp_q3 * grid_sp_q3) endif - CS%Ah_bg_xy(I,J) = MAX(Ah, Ah_vel_scale * grid_sp_q2 * sqrt(grid_sp_q2)) + if (CS%Re_Ah > 0.0) CS%Re_Ah_const_xy(i,j) = grid_sp_q3 / CS%Re_Ah if (Ah_time_scale > 0.) CS%Ah_bg_xy(i,j) = & MAX(CS%Ah_bg_xy(i,j), (grid_sp_q2 * grid_sp_q2) / Ah_time_scale) if (CS%bound_Ah .and. .not.CS%better_bound_Ah) then @@ -1928,7 +1935,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) endif enddo ; enddo endif - ! The Laplacian bounds should avoid overshoots when CS%bound_coef < 1. if (CS%Laplacian .and. CS%better_bound_Kh) then Idt = 1.0 / dt @@ -1957,7 +1963,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) call Bchksum(CS%Kh_Max_xy, "Kh_Max_xy", G%HI, haloshift=0, scale=US%L_to_m**2*US%s_to_T) endif endif - ! The biharmonic bounds should avoid overshoots when CS%bound_coef < 0.5, but ! empirically work for CS%bound_coef <~ 1.0 if (CS%biharmonic .and. CS%better_bound_Ah) then @@ -1967,7 +1972,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) CS%dy2h(i,j) * CS%DY_dxT(i,j) * (G%IdyCu(I,j) + G%IdyCu(I-1,j)) ) + & CS%Idx2dyCu(I,j)*(CS%dx2q(I,J) * CS%DX_dyBu(I,J) * (G%IdxCu(I,j+1) + G%IdxCu(I,j)) + & CS%dx2q(I,J-1)*CS%DX_dyBu(I,J-1)*(G%IdxCu(I,j) + G%IdxCu(I,j-1)) ) ) - u0v(I,j) = (CS%Idxdy2u(I,j)*(CS%dy2h(i+1,j)*CS%DX_dyT(i+1,j)*(G%IdxCv(i+1,J) + G%IdxCv(i+1,J-1)) + & CS%dy2h(i,j) * CS%DX_dyT(i,j) * (G%IdxCv(i,J) + G%IdxCv(i,J-1)) ) + & CS%Idx2dyCu(I,j)*(CS%dx2q(I,J) * CS%DY_dxBu(I,J) * (G%IdyCv(i+1,J) + G%IdyCv(i,J)) + & @@ -1978,13 +1982,11 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) CS%dy2q(I-1,J)*CS%DX_dyBu(I-1,J)*(G%IdxCu(I-1,j+1) + G%IdxCu(I-1,j)) ) + & CS%Idx2dyCv(i,J)*(CS%dx2h(i,j+1)*CS%DY_dxT(i,j+1)*(G%IdyCu(I,j+1) + G%IdyCu(I-1,j+1)) + & CS%dx2h(i,j) * CS%DY_dxT(i,j) * (G%IdyCu(I,j) + G%IdyCu(I-1,j)) ) ) - v0v(i,J) = (CS%Idxdy2v(i,J)*(CS%dy2q(I,J) * CS%DY_dxBu(I,J) * (G%IdyCv(i+1,J) + G%IdyCv(i,J)) + & CS%dy2q(I-1,J)*CS%DY_dxBu(I-1,J)*(G%IdyCv(i,J) + G%IdyCv(i-1,J)) ) + & CS%Idx2dyCv(i,J)*(CS%dx2h(i,j+1)*CS%DX_dyT(i,j+1)*(G%IdxCv(i,J+1) + G%IdxCv(i,J)) + & CS%dx2h(i,j) * CS%DX_dyT(i,j) * (G%IdxCv(i,J) + G%IdxCv(i,J-1)) ) ) enddo ; enddo - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 denom = max( & (CS%dy2h(i,j) * & @@ -1999,7 +2001,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) if (denom > 0.0) & CS%Ah_Max_xx(I,J) = CS%bound_coef * 0.5 * Idt / denom enddo ; enddo - do J=js-1,Jeq ; do I=is-1,Ieq denom = max( & (CS%dx2q(I,J) * & @@ -2019,74 +2020,62 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) call Bchksum(CS%Ah_Max_xy, "Ah_Max_xy", G%HI, haloshift=0, scale=US%L_to_m**4*US%s_to_T) endif endif - ! Register fields for output from this module. - CS%id_diffu = register_diag_field('ocean_model', 'diffu', diag%axesCuL, Time, & 'Zonal Acceleration from Horizontal Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2) - CS%id_diffv = register_diag_field('ocean_model', 'diffv', diag%axesCvL, Time, & 'Meridional Acceleration from Horizontal Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2) - if (CS%biharmonic) then CS%id_Ah_h = register_diag_field('ocean_model', 'Ahh', diag%axesTL, Time, & 'Biharmonic Horizontal Viscosity at h Points', 'm4 s-1', conversion=US%L_to_m**4*US%s_to_T, & cmor_field_name='difmxybo', & cmor_long_name='Ocean lateral biharmonic viscosity', & cmor_standard_name='ocean_momentum_xy_biharmonic_diffusivity') - CS%id_Ah_q = register_diag_field('ocean_model', 'Ahq', diag%axesBL, Time, & 'Biharmonic Horizontal Viscosity at q Points', 'm4 s-1', conversion=US%L_to_m**4*US%s_to_T) + CS%id_grid_Re_Ah = register_diag_field('ocean_model', 'grid_Re_Ah', diag%axesTL, Time, & + 'Grid Reynolds number for the Biharmonic horizontal viscosity at h points', 'nondim') endif - if (CS%Laplacian) then CS%id_Kh_h = register_diag_field('ocean_model', 'Khh', diag%axesTL, Time, & 'Laplacian Horizontal Viscosity at h Points', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T, & cmor_field_name='difmxylo', & cmor_long_name='Ocean lateral Laplacian viscosity', & cmor_standard_name='ocean_momentum_xy_laplacian_diffusivity') - CS%id_Kh_q = register_diag_field('ocean_model', 'Khq', diag%axesBL, Time, & 'Laplacian Horizontal Viscosity at q Points', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T) - - if (CS%Leith_Kh) then - CS%id_vort_xy_q = register_diag_field('ocean_model', 'vort_xy_q', diag%axesBL, Time, & - 'Vertical vorticity at q Points', 's-1', conversion=US%s_to_T) - - CS%id_div_xx_h = register_diag_field('ocean_model', 'div_xx_h', diag%axesTL, Time, & - 'Horizontal divergence at h Points', 's-1', conversion=US%s_to_T) - endif - + CS%id_grid_Re_Kh = register_diag_field('ocean_model', 'grid_Re_Kh', diag%axesTL, Time, & + 'Grid Reynolds number for the Laplacian horizontal viscosity at h points', 'nondim') + CS%id_vort_xy_q = register_diag_field('ocean_model', 'vort_xy_q', diag%axesBL, Time, & + 'Vertical vorticity at q Points', 's-1', conversion=US%s_to_T) + CS%id_div_xx_h = register_diag_field('ocean_model', 'div_xx_h', diag%axesTL, Time, & + 'Horizontal divergence at h Points', 's-1', conversion=US%s_to_T) + CS%id_sh_xy_q = register_diag_field('ocean_model', 'sh_xy_q', diag%axesBL, Time, & + 'Shearing strain at q Points', 's-1', conversion=US%s_to_T) + CS%id_sh_xx_h = register_diag_field('ocean_model', 'sh_xx_h', diag%axesTL, Time, & + 'Horizontal tension at h Points', 's-1', conversion=US%s_to_T) endif - if (CS%use_GME) then CS%id_GME_coeff_h = register_diag_field('ocean_model', 'GME_coeff_h', diag%axesTL, Time, & 'GME coefficient at h Points', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T) - CS%id_GME_coeff_q = register_diag_field('ocean_model', 'GME_coeff_q', diag%axesBL, Time, & 'GME coefficient at q Points', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T) - CS%id_FrictWork_GME = register_diag_field('ocean_model','FrictWork_GME',diag%axesTL,Time,& 'Integral work done by lateral friction terms in GME (excluding diffusion of energy)', & 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) endif - CS%id_FrictWork = register_diag_field('ocean_model','FrictWork',diag%axesTL,Time,& 'Integral work done by lateral friction terms', & 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) - CS%id_FrictWorkIntz = register_diag_field('ocean_model','FrictWorkIntz',diag%axesT1,Time, & 'Depth integrated work done by lateral friction', & 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2, & cmor_field_name='dispkexyfo', & cmor_long_name='Depth integrated ocean kinetic energy dissipation due to lateral friction',& cmor_standard_name='ocean_kinetic_energy_dissipation_per_unit_area_due_to_xy_friction') - if (CS%Laplacian .or. get_all) then endif - end subroutine hor_visc_init - !> Calculates factors in the anisotropic orientation tensor to be align with the grid. !! With n1=1 and n2=0, this recovers the approach of Large et al, 2001. subroutine align_aniso_tensor_to_grid(CS, n1, n2) @@ -2095,18 +2084,14 @@ subroutine align_aniso_tensor_to_grid(CS, n1, n2) real, intent(in) :: n2 !< j-component of direction vector [nondim] ! Local variables real :: recip_n2_norm - ! For normalizing n=(n1,n2) in case arguments are not a unit vector recip_n2_norm = n1**2 + n2**2 if (recip_n2_norm > 0.) recip_n2_norm = 1./recip_n2_norm - CS%n1n2_h(:,:) = 2. * ( n1 * n2 ) * recip_n2_norm CS%n1n2_q(:,:) = 2. * ( n1 * n2 ) * recip_n2_norm CS%n1n1_m_n2n2_h(:,:) = ( n1 * n1 - n2 * n2 ) * recip_n2_norm CS%n1n1_m_n2n2_q(:,:) = ( n1 * n1 - n2 * n2 ) * recip_n2_norm - end subroutine align_aniso_tensor_to_grid - !> Apply a 1-1-4-1-1 Laplacian filter one time on GME diffusive flux to reduce any !! horizontal two-grid-point noise subroutine smooth_GME(CS,G,GME_flux_h,GME_flux_q) @@ -2117,15 +2102,12 @@ subroutine smooth_GME(CS,G,GME_flux_h,GME_flux_q) !! at h points real, dimension(SZIB_(G),SZJB_(G)), optional, intent(inout) :: GME_flux_q!< GME diffusive flux !! at q points - ! local variables real, dimension(SZI_(G),SZJ_(G)) :: GME_flux_h_original real, dimension(SZIB_(G),SZJB_(G)) :: GME_flux_q_original real :: wc, ww, we, wn, ws ! averaging weights for smoothing integer :: i, j, k, s - do s=1,1 - ! Update halos if (present(GME_flux_h)) then call pass_var(GME_flux_h, G%Domain) @@ -2135,14 +2117,12 @@ subroutine smooth_GME(CS,G,GME_flux_h,GME_flux_q) do i = G%isc, G%iec ! skip land points if (G%mask2dT(i,j)==0.) cycle - ! compute weights ww = 0.125 * G%mask2dT(i-1,j) we = 0.125 * G%mask2dT(i+1,j) ws = 0.125 * G%mask2dT(i,j-1) wn = 0.125 * G%mask2dT(i,j+1) wc = 1.0 - (ww+we+wn+ws) - GME_flux_h(i,j) = wc * GME_flux_h_original(i,j) & + ww * GME_flux_h_original(i-1,j) & + we * GME_flux_h_original(i+1,j) & @@ -2150,7 +2130,6 @@ subroutine smooth_GME(CS,G,GME_flux_h,GME_flux_q) + wn * GME_flux_h_original(i,j+1) enddo; enddo endif - ! Update halos if (present(GME_flux_q)) then call pass_var(GME_flux_q, G%Domain, position=CORNER, complete=.true.) @@ -2160,14 +2139,12 @@ subroutine smooth_GME(CS,G,GME_flux_h,GME_flux_q) do I = G%IscB, G%IecB ! skip land points if (G%mask2dBu(I,J)==0.) cycle - ! compute weights ww = 0.125 * G%mask2dBu(I-1,J) we = 0.125 * G%mask2dBu(I+1,J) ws = 0.125 * G%mask2dBu(I,J-1) wn = 0.125 * G%mask2dBu(I,J+1) wc = 1.0 - (ww+we+wn+ws) - GME_flux_q(I,J) = wc * GME_flux_q_original(I,J) & + ww * GME_flux_q_original(I-1,J) & + we * GME_flux_q_original(I+1,J) & @@ -2175,24 +2152,20 @@ subroutine smooth_GME(CS,G,GME_flux_h,GME_flux_q) + wn * GME_flux_q_original(I,J+1) enddo; enddo endif - enddo ! s-loop - end subroutine smooth_GME - !> Deallocates any variables allocated in hor_visc_init. subroutine hor_visc_end(CS) type(hor_visc_CS), pointer :: CS !< The control structure returned by a !! previous call to hor_visc_init. - if (CS%Laplacian .or. CS%biharmonic) then DEALLOC_(CS%dx2h) ; DEALLOC_(CS%dx2q) ; DEALLOC_(CS%dy2h) ; DEALLOC_(CS%dy2q) DEALLOC_(CS%dx_dyT) ; DEALLOC_(CS%dy_dxT) ; DEALLOC_(CS%dx_dyBu) ; DEALLOC_(CS%dy_dxBu) DEALLOC_(CS%reduction_xx) ; DEALLOC_(CS%reduction_xy) endif - if (CS%Laplacian) then DEALLOC_(CS%Kh_bg_xx) ; DEALLOC_(CS%Kh_bg_xy) + DEALLOC_(CS%grid_sp_h2) if (CS%bound_Kh) then DEALLOC_(CS%Kh_Max_xx) ; DEALLOC_(CS%Kh_Max_xy) endif @@ -2203,8 +2176,8 @@ subroutine hor_visc_end(CS) DEALLOC_(CS%Laplac3_const_xx) ; DEALLOC_(CS%Laplac3_const_xy) endif endif - if (CS%biharmonic) then + DEALLOC_(CS%grid_sp_h3) DEALLOC_(CS%Idx2dyCu) ; DEALLOC_(CS%Idx2dyCv) DEALLOC_(CS%Idxdy2u) ; DEALLOC_(CS%Idxdy2v) DEALLOC_(CS%Ah_bg_xx) ; DEALLOC_(CS%Ah_bg_xy) @@ -2212,14 +2185,14 @@ subroutine hor_visc_end(CS) DEALLOC_(CS%Ah_Max_xx) ; DEALLOC_(CS%Ah_Max_xy) endif if (CS%Smagorinsky_Ah) then - DEALLOC_(CS%Biharm5_const_xx) ; DEALLOC_(CS%Biharm5_const_xy) - ! if (CS%bound_Coriolis) then - ! DEALLOC_(CS%Biharm5_const2_xx) ; DEALLOC_(CS%Biharm5_const2_xy) - ! endif + DEALLOC_(CS%Biharm6_const_xx) ; DEALLOC_(CS%Biharm6_const_xy) endif if (CS%Leith_Ah) then DEALLOC_(CS%Biharm_const_xx) ; DEALLOC_(CS%Biharm_const_xy) endif + if (CS%Re_Ah > 0.0) then + DEALLOC_(CS%Re_Ah_const_xx) ; DEALLOC_(CS%Re_Ah_const_xy) + endif endif if (CS%anisotropic) then DEALLOC_(CS%n1n2_h) @@ -2228,10 +2201,7 @@ subroutine hor_visc_end(CS) DEALLOC_(CS%n1n1_m_n2n2_q) endif deallocate(CS) - end subroutine hor_visc_end - - !> \namespace mom_hor_visc !! !! This module contains the subroutine horizontal_viscosity() that calculates the @@ -2532,5 +2502,4 @@ end subroutine hor_visc_end !! Smith, R.D., and McWilliams, J.C., 2003: Anisotropic horizontal viscosity for !! ocean models. Ocean Modelling, 5(2), 129-156. !! https://doi.org/10.1016/S1463-5003(02)00016-1 - end module MOM_hor_visc diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 37e549f3f1..a5a545e85d 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -751,8 +751,6 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo type(ocean_grid_type), intent(in) :: G !< Ocean 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 -! real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal flow [L T-1 ~> m s-1] -! real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional flow [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] integer, intent(in) :: k !< Layer for which to calculate vorticity magnitude real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: div_xx_dx !< x-derivative of horizontal divergence @@ -763,15 +761,6 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo !! (d/dx(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1] real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: vort_xy_dy !< y-derivative of vertical vorticity !! (d/dy(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1] -! real, dimension(SZI_(G),SZJ_(G)), intent(out) :: Leith_Kh_h !< Leith Laplacian viscosity - !! at h-points [L2 T-1 ~> m2 s-1] -! real, dimension(SZIB_(G),SZJB_(G)), intent(out) :: Leith_Kh_q !< Leith Laplacian viscosity - !! at q-points [L2 T-1 ~> m2 s-1] -! real, dimension(SZI_(G),SZJ_(G)), intent(out) :: Leith_Ah_h !< Leith bi-harmonic viscosity - !! at h-points [L4 T-1 ~> m4 s-1] -! real, dimension(SZIB_(G),SZJB_(G)), intent(out) :: Leith_Ah_q !< Leith bi-harmonic viscosity - !! at q-points [L4 T-1 ~> m4 s-1] - ! Local variables real, dimension(SZI_(G),SZJB_(G)) :: & dslopey_dz, & ! z-derivative of y-slope at v-points [Z-1 ~> m-1] @@ -799,16 +788,9 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo inv_PI3 = 1.0/((4.0*atan(1.0))**3) - !### I believe this halo update to be unnecessary. -RWH - call pass_var(h, G%Domain) - if ((k > 1) .and. (k < nz)) then - ! Add in stretching term for the QG Leith vsicosity -! if (CS%use_QG_Leith) then - - !### do j=js-1,je+1 ; do I=is-2,Ieq+1 - do j=js-2,Jeq+2 ; do I=is-2,Ieq+1 + do j=js-1,je+1 ; do I=is-2,Ieq+1 h_at_slope_above = 2. * ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) * h(i+1,j,k) ) / & ( ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) + h(i+1,j,k) ) & + ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k-1) + h(i+1,j,k-1) ) + GV%H_subroundoff ) @@ -820,8 +802,7 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo h_at_u(I,j) = 2. * ( h_at_slope_above * h_at_slope_below ) * Ih enddo ; enddo - !### do J=js-2,Jeq+1 ; do i=is-1,ie+1 - do J=js-2,Jeq+1 ; do i=is-2,Ieq+2 + do J=js-2,Jeq+1 ; do i=is-1,ie+1 h_at_slope_above = 2. * ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) * h(i,j+1,k) ) / & ( ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) + h(i,j+1,k) ) & + ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k-1) + h(i,j+1,k-1) ) + GV%H_subroundoff ) @@ -833,8 +814,7 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo h_at_v(i,J) = 2. * ( h_at_slope_above * h_at_slope_below ) * Ih enddo ; enddo - !### do J=js-1,je ; do i=is-1,Ieq+1 - do J=js-2,Jeq+1 ; do i=is-1,Ieq+1 + do J=js-1,je ; do i=is-1,Ieq+1 f = 0.5 * ( G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J) ) vort_xy_dx(i,J) = vort_xy_dx(i,J) - f * US%L_to_Z * & ( ( h_at_u(I,j) * dslopex_dz(I,j) + h_at_u(I-1,j+1) * dslopex_dz(I-1,j+1) ) & @@ -842,33 +822,25 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo ( ( h_at_u(I,j) + h_at_u(I-1,j+1) ) + ( h_at_u(I-1,j) + h_at_u(I,j+1) ) + GV%H_subroundoff) enddo ; enddo - !### do j=js-1,Jeq+1 ; do I=is-1,ie - do j=js-1,Jeq+1 ; do I=is-2,Ieq+1 + do j=js-1,Jeq+1 ; do I=is-1,ie f = 0.5 * ( G%CoriolisBu(I,J) + G%CoriolisBu(I,J-1) ) - !### I think that this should be vort_xy_dy(I,j) = vort_xy_dy(I,j) - f * & - vort_xy_dy(I,j) = vort_xy_dx(I,j) - f * US%L_to_Z * & + vort_xy_dy(I,j) = vort_xy_dy(I,j) - f * US%L_to_Z * & ( ( h_at_v(i,J) * dslopey_dz(i,J) + h_at_v(i+1,J-1) * dslopey_dz(i+1,J-1) ) & + ( h_at_v(i,J-1) * dslopey_dz(i,J-1) + h_at_v(i+1,J) * dslopey_dz(i+1,J) ) ) / & ( ( h_at_v(i,J) + h_at_v(i+1,J-1) ) + ( h_at_v(i,J-1) + h_at_v(i+1,J) ) + GV%H_subroundoff) enddo ; enddo endif ! k > 1 - !### I believe this halo update to be unnecessary. -RWH - call pass_vector(vort_xy_dy,vort_xy_dx,G%Domain) - if (CS%use_QG_Leith_GM) then do j=js,je ; do I=is-1,Ieq - !### These expressions are not rotationally symmetric. Add parentheses and regroup, as in: - ! grad_vort_mag_u(I,j) = SQRT(vort_xy_dy(I,j)**2 + (0.25*((vort_xy_dx(i,J) + vort_xy_dx(i+1,J-1)) + - ! (vort_xy_dx(i+1,J) + vort_xy_dx(i,J-1))))**2 ) - grad_vort_mag_u(I,j) = SQRT(vort_xy_dy(I,j)**2 + (0.25*(vort_xy_dx(i,J) + vort_xy_dx(i+1,J) & - + vort_xy_dx(i,J-1) + vort_xy_dx(i+1,J-1)))**2) - grad_div_mag_u(I,j) = SQRT(div_xx_dx(I,j)**2 + (0.25*(div_xx_dy(i,J) + div_xx_dy(i+1,J) & - + div_xx_dy(i,J-1) + div_xx_dy(i+1,J-1)))**2) + grad_vort_mag_u(I,j) = SQRT(vort_xy_dy(I,j)**2 + (0.25*((vort_xy_dx(i,J) + vort_xy_dx(i+1,J-1)) & + + (vort_xy_dx(i+1,J) + vort_xy_dx(i,J-1))))**2) + grad_div_mag_u(I,j) = SQRT(div_xx_dx(I,j)**2 + (0.25*((div_xx_dy(i,J) + div_xx_dy(i+1,J-1)) & + + (div_xx_dy(i+1,J) + div_xx_dy(i,J-1))))**2) if (CS%use_beta_in_QG_Leith) then - beta_u(I,j) = sqrt( (0.5*(G%dF_dx(i,j)+G%dF_dx(i+1,j))**2) + & - (0.5*(G%dF_dy(i,j)+G%dF_dy(i+1,j))**2) ) + beta_u(I,j) = sqrt((0.5*(G%dF_dx(i,j)+G%dF_dx(i+1,j))**2) + & + (0.5*(G%dF_dy(i,j)+G%dF_dy(i+1,j))**2)) CS%KH_u_QG(I,j,k) = MIN(grad_vort_mag_u(I,j) + grad_div_mag_u(I,j), 3.0*beta_u(I,j)) * & CS%Laplac3_const_u(I,j) * inv_PI3 else @@ -878,14 +850,13 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo enddo ; enddo do J=js-1,Jeq ; do i=is,ie - !### These expressions are not rotationally symmetric. Add parentheses and regroup. - grad_vort_mag_v(i,J) = SQRT(vort_xy_dx(i,J)**2 + (0.25*(vort_xy_dy(I,j) + vort_xy_dy(I-1,j) & - + vort_xy_dy(I,j+1) + vort_xy_dy(I-1,j+1)))**2) - grad_div_mag_v(i,J) = SQRT(div_xx_dy(i,J)**2 + (0.25*(div_xx_dx(I,j) + div_xx_dx(I-1,j) & - + div_xx_dx(I,j+1) + div_xx_dx(I-1,j+1)))**2) + grad_vort_mag_v(i,J) = SQRT(vort_xy_dx(i,J)**2 + (0.25*((vort_xy_dy(I,j) + vort_xy_dy(I-1,j+1)) & + + (vort_xy_dy(I,j+1) + vort_xy_dy(I-1,j))))**2) + grad_div_mag_v(i,J) = SQRT(div_xx_dy(i,J)**2 + (0.25*((div_xx_dx(I,j) + div_xx_dx(I-1,j+1)) & + + (div_xx_dx(I,j+1) + div_xx_dx(I-1,j))))**2) if (CS%use_beta_in_QG_Leith) then - beta_v(i,J) = sqrt( (0.5*(G%dF_dx(i,j)+G%dF_dx(i,j+1))**2) + & - (0.5*(G%dF_dy(i,j)+G%dF_dy(i,j+1))**2) ) + beta_v(i,J) = sqrt((0.5*(G%dF_dx(i,j)+G%dF_dx(i,j+1))**2) + & + (0.5*(G%dF_dy(i,j)+G%dF_dy(i,j+1))**2)) CS%KH_v_QG(i,J,k) = MIN(grad_vort_mag_v(i,J) + grad_div_mag_v(i,J), 3.0*beta_v(i,J)) * & CS%Laplac3_const_v(i,J) * inv_PI3 else @@ -1297,13 +1268,12 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) do j=Jsq,Jeq+1 ; do I=is-1,Ieq ! Static factors in the Leith schemes grid_sp_u2 = G%dyCu(I,j)*G%dxCu(I,j) - grid_sp_u3 = sqrt(grid_sp_u2) + grid_sp_u3 = grid_sp_u2*sqrt(grid_sp_u2) CS%Laplac3_const_u(I,j) = Leith_Lap_const * grid_sp_u3 enddo ; enddo do j=js-1,Jeq ; do I=Isq,Ieq+1 ! Static factors in the Leith schemes - !### The second factor here is wrong. It should be G%dxCv(i,J). - grid_sp_v2 = G%dyCv(i,J)*G%dxCu(i,J) + grid_sp_v2 = G%dyCv(i,J)*G%dxCv(i,J) grid_sp_v3 = grid_sp_v2*sqrt(grid_sp_v2) CS%Laplac3_const_v(i,J) = Leith_Lap_const * grid_sp_v3 enddo ; enddo diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index 3b7420aa54..01a39d394b 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -17,6 +17,8 @@ module MOM_CVMix_KPP use MOM_verticalGrid, only : verticalGrid_type use MOM_wave_interface, only : wave_parameters_CS, Get_Langmuir_Number use MOM_domains, only : pass_var +use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end +use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE use CVMix_kpp, only : CVMix_init_kpp, CVMix_put_kpp, CVMix_get_kpp_real use CVMix_kpp, only : CVMix_coeffs_kpp @@ -169,6 +171,10 @@ module MOM_CVMix_KPP end type KPP_CS +!>@{ CPU time clocks +integer :: id_clock_KPP_calc, id_clock_KPP_compute_BLD, id_clock_KPP_smoothing +!!@} + #define __DO_SAFETY_CHECKS__ contains @@ -226,11 +232,17 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) 'The number of times the 1-1-4-1-1 Laplacian filter is applied on '// & 'OBL depth.', & default=0) + if (CS%n_smooth > G%domain%nihalo) then + call MOM_error(FATAL,'KPP smoothing number (N_SMOOTH) cannot be greater than NIHALO.') + elseif (CS%n_smooth > G%domain%njhalo) then + call MOM_error(FATAL,'KPP smoothing number (N_SMOOTH) cannot be greater than NJHALO.') + endif if (CS%n_smooth > 0) then call get_param(paramFile, mdl, 'DEEPEN_ONLY_VIA_SMOOTHING', CS%deepen_only, & 'If true, apply OBLdepth smoothing at a cell only if the OBLdepth '// & 'gets deeper via smoothing.', & default=.false.) + id_clock_KPP_smoothing = cpu_clock_id('(Ocean KPP BLD smoothing)', grain=CLOCK_ROUTINE) endif call get_param(paramFile, mdl, 'RI_CRIT', CS%Ri_crit, & 'Critical bulk Richardson number used to define depth of the '// & @@ -579,6 +591,8 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) if (CS%id_EnhK > 0) allocate( CS%EnhK( SZI_(G), SZJ_(G), SZK_(G)+1 ) ) if (CS%id_EnhK > 0) CS%EnhK(:,:,:) = 0. + id_clock_KPP_calc = cpu_clock_id('Ocean KPP calculate)', grain=CLOCK_MODULE) + id_clock_KPP_compute_BLD = cpu_clock_id('(Ocean KPP comp BLD)', grain=CLOCK_ROUTINE) end function KPP_init @@ -640,6 +654,7 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, & if (CS%id_Kd_in > 0) call post_data(CS%id_Kd_in, Kt, CS%diag) + call cpu_clock_begin(id_clock_KPP_calc) buoy_scale = US%L_to_m**2*US%s_to_T**3 !$OMP parallel do default(none) firstprivate(nonLocalTrans) & @@ -860,6 +875,7 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, & enddo ! i enddo ! j + call cpu_clock_end(id_clock_KPP_calc) #ifdef __DO_SAFETY_CHECKS__ if (CS%debug) then @@ -960,21 +976,23 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl endif #endif + call cpu_clock_begin(id_clock_KPP_compute_BLD) + ! some constants GoRho = US%L_T_to_m_s**2*US%m_to_Z * GV%g_Earth / GV%Rho0 buoy_scale = US%L_to_m**2*US%s_to_T**3 ! loop over horizontal points on processor - !GOMP parallel do default(none) private(surfFricVel, iFaceHeight, hcorr, dh, cellHeight, & - !GOMP surfBuoyFlux, U_H, V_H, u, v, Coriolis, pRef, SLdepth_0d, & - !GOMP ksfc, surfHtemp, surfHsalt, surfHu, surfHv, surfHuS, & - !GOMP surfHvS, hTot, delH, surftemp, surfsalt, surfu, surfv, & - !GOMP surfUs, surfVs, Uk, Vk, deltaU2, km1, kk, pres_1D, & - !GOMP Temp_1D, salt_1D, surfBuoyFlux2, MLD_GUESS, LA, rho_1D, & - !GOMP deltarho, N2_1d, ws_1d, LangEnhVT2, enhvt2, wst, & - !GOMP BulkRi_1d, zBottomMinusOffset) & - !GOMP shared(G, GV, CS, US, uStar, h, buoy_scale, buoyFlux, & - !GOMP Temp, Salt, waves, tv, GoRho) + !$OMP parallel do default(none) private(surfFricVel, iFaceHeight, hcorr, dh, cellHeight, & + !$OMP surfBuoyFlux, U_H, V_H, Coriolis, pRef, SLdepth_0d, & + !$OMP ksfc, surfHtemp, surfHsalt, surfHu, surfHv, surfHuS, & + !$OMP surfHvS, hTot, delH, surftemp, surfsalt, surfu, surfv, & + !$OMP surfUs, surfVs, Uk, Vk, deltaU2, km1, kk, pres_1D, & + !$OMP Temp_1D, salt_1D, surfBuoyFlux2, MLD_GUESS, LA, rho_1D, & + !$OMP deltarho, N2_1d, ws_1d, LangEnhVT2, enhvt2, wst, & + !$OMP BulkRi_1d, zBottomMinusOffset) & + !$OMP shared(G, GV, CS, US, uStar, h, buoy_scale, buoyFlux, & + !$OMP Temp, Salt, waves, tv, GoRho, u, v) do j = G%jsc, G%jec do i = G%isc, G%iec @@ -1245,6 +1263,8 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl enddo enddo + call cpu_clock_end(id_clock_KPP_compute_BLD) + ! send diagnostics to post_data if (CS%id_BulkRi > 0) call post_data(CS%id_BulkRi, CS%BulkRi, CS%diag) if (CS%id_N > 0) call post_data(CS%id_N, CS%N, CS%diag) @@ -1274,7 +1294,7 @@ subroutine KPP_smooth_BLD(CS,G,GV,h) real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer/level thicknesses [H ~> m or kg m-2] ! local - real, dimension(SZI_(G),SZJ_(G)) :: OBLdepth_original ! Original OBL depths computed by CVMix + real, dimension(SZI_(G),SZJ_(G)) :: OBLdepth_prev ! OBLdepth before s.th smoothing iteration real, dimension( G%ke ) :: cellHeight ! Cell center heights referenced to surface [m] ! (negative in the ocean) real, dimension( G%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface [m] @@ -1284,15 +1304,20 @@ subroutine KPP_smooth_BLD(CS,G,GV,h) real :: hcorr ! A cumulative correction arising from inflation of vanished layers [m] integer :: i, j, k, s - do s=1,CS%n_smooth + call cpu_clock_begin(id_clock_KPP_smoothing) + + ! Update halos + call pass_var(CS%OBLdepth, G%Domain, halo=CS%n_smooth) - ! Update halos - call pass_var(CS%OBLdepth, G%Domain) + if (CS%id_OBLdepth_original > 0) CS%OBLdepth_original = CS%OBLdepth + + do s=1,CS%n_smooth - OBLdepth_original = CS%OBLdepth - if (CS%id_OBLdepth_original > 0) CS%OBLdepth_original = OBLdepth_original + OBLdepth_prev = CS%OBLdepth ! apply smoothing on OBL depth + !$OMP parallel do default(none) shared(G, GV, CS, h, OBLdepth_prev) & + !$OMP private(wc, ww, we, wn, ws, dh, hcorr, cellHeight, iFaceHeight) do j = G%jsc, G%jec do i = G%isc, G%iec @@ -1319,14 +1344,14 @@ subroutine KPP_smooth_BLD(CS,G,GV,h) wn = 0.125 * G%mask2dT(i,j+1) wc = 1.0 - (ww+we+wn+ws) - CS%OBLdepth(i,j) = wc * OBLdepth_original(i,j) & - + ww * OBLdepth_original(i-1,j) & - + we * OBLdepth_original(i+1,j) & - + ws * OBLdepth_original(i,j-1) & - + wn * OBLdepth_original(i,j+1) + CS%OBLdepth(i,j) = wc * OBLdepth_prev(i,j) & + + ww * OBLdepth_prev(i-1,j) & + + we * OBLdepth_prev(i+1,j) & + + ws * OBLdepth_prev(i,j-1) & + + wn * OBLdepth_prev(i,j+1) ! Apply OBLdepth smoothing at a cell only if the OBLdepth gets deeper via smoothing. - if (CS%deepen_only) CS%OBLdepth(i,j) = max(CS%OBLdepth(i,j),CS%OBLdepth_original(i,j)) + if (CS%deepen_only) CS%OBLdepth(i,j) = max(CS%OBLdepth(i,j), OBLdepth_prev(i,j)) ! prevent OBL depths deeper than the bathymetric depth CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -iFaceHeight(G%ke+1) ) ! no deeper than bottom @@ -1336,30 +1361,7 @@ subroutine KPP_smooth_BLD(CS,G,GV,h) enddo ! s-loop - ! Update kOBL for smoothed OBL depths - do j = G%jsc, G%jec - do i = G%isc, G%iec - - ! skip land points - if (G%mask2dT(i,j)==0.) cycle - - iFaceHeight(1) = 0.0 ! BBL is all relative to the surface - hcorr = 0. - do k=1,G%ke - - ! cell center and cell bottom in meters (negative values in the ocean) - dh = h(i,j,k) * GV%H_to_m ! Nominal thickness to use for increment - dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0) - hcorr = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0 - dh = max( dh, CS%min_thickness ) ! Limit increment dh>=min_thickness - cellHeight(k) = iFaceHeight(k) - 0.5 * dh - iFaceHeight(k+1) = iFaceHeight(k) - dh - enddo - - CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) ) - - enddo - enddo + call cpu_clock_end(id_clock_KPP_smoothing) end subroutine KPP_smooth_BLD @@ -1380,6 +1382,7 @@ subroutine KPP_get_BLD(CS, BLD, G, US, m_to_BLD_units) scale = US%m_to_Z ; if (present(m_to_BLD_units)) scale = m_to_BLD_units + !$OMP parallel do default(none) shared(BLD, CS, G, scale) do j = G%jsc, G%jec ; do i = G%isc, G%iec BLD(i,j) = scale * CS%OBLdepth(i,j) enddo ; enddo @@ -1418,6 +1421,7 @@ subroutine KPP_NonLocalTransport_temp(CS, G, GV, h, nonLocalTrans, surfFlux, & ! Update tracer due to non-local redistribution of surface flux if (CS%applyNonLocalTrans) then + !$OMP parallel do default(none) shared(dt, scalar, dtracer, G) do k = 1, G%ke do j = G%jsc, G%jec do i = G%isc, G%iec @@ -1432,6 +1436,7 @@ subroutine KPP_NonLocalTransport_temp(CS, G, GV, h, nonLocalTrans, surfFlux, & if (CS%id_NLT_dTdt > 0) call post_data(CS%id_NLT_dTdt, dtracer, CS%diag) if (CS%id_NLT_temp_budget > 0) then dtracer(:,:,:) = 0.0 + !$OMP parallel do default(none) shared(dtracer, nonLocalTrans, surfFlux, C_p, G, GV) do k = 1, G%ke do j = G%jsc, G%jec do i = G%isc, G%iec @@ -1477,6 +1482,7 @@ subroutine KPP_NonLocalTransport_saln(CS, G, GV, h, nonLocalTrans, surfFlux, dt, ! Update tracer due to non-local redistribution of surface flux if (CS%applyNonLocalTrans) then + !$OMP parallel do default(none) shared(G, dt, scalar, dtracer) do k = 1, G%ke do j = G%jsc, G%jec do i = G%isc, G%iec @@ -1491,6 +1497,7 @@ subroutine KPP_NonLocalTransport_saln(CS, G, GV, h, nonLocalTrans, surfFlux, dt, if (CS%id_NLT_dSdt > 0) call post_data(CS%id_NLT_dSdt, dtracer, CS%diag) if (CS%id_NLT_saln_budget > 0) then dtracer(:,:,:) = 0.0 + !$OMP parallel do default(none) shared(G, GV, dtracer, nonLocalTrans, surfFlux) do k = 1, G%ke do j = G%jsc, G%jec do i = G%isc, G%iec diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index 801e573023..708d6a7f46 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -687,7 +687,7 @@ subroutine calculate_tidal_mixing(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, US, C real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(inout) :: Kd_lay !< The diapycnal diffusvity in layers [Z2 T-1 ~> m2 s-1]. real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), & - optional, intent(inout) :: Kd_int !< The diapycnal diffusvity at interfaces, + optional, intent(inout) :: Kd_int !< The diapycnal diffusvity at interfaces, !! [Z2 T-1 ~> m2 s-1]. real, intent(in) :: Kd_max !< The maximum increment for diapycnal !! diffusivity due to TKE-based processes, @@ -698,7 +698,7 @@ subroutine calculate_tidal_mixing(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, US, C if (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation .or. CS%Lowmode_itidal_dissipation) then if (CS%use_CVMix_tidal) then - call calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kv) + call calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kd_int, Kv) else call add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, & G, GV, US, CS, N2_lay, Kd_lay, Kd_int, Kd_max) @@ -709,7 +709,7 @@ end subroutine calculate_tidal_mixing !> Calls the CVMix routines to compute tidal dissipation and to add the effect of internal-tide-driven !! mixing to the interface diffusivities. -subroutine calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kv) +subroutine calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kd_int, Kv) integer, intent(in) :: j !< The j-index to work on type(ocean_grid_type), intent(in) :: G !< Grid structure. type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure @@ -721,6 +721,9 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kv) intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(inout) :: Kd_lay!< The diapycnal diffusivities in the layers [Z2 T-1 ~> m2 s-1]. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), & + optional, intent(inout) :: Kd_int !< The diapycnal diffusvity at interfaces, + !! [Z2 T-1 ~> m2 s-1]. real, dimension(:,:,:), pointer :: Kv !< The "slow" vertical viscosity at each interface !! (not layer!) [Z2 T-1 ~> m2 s-1]. ! Local variables @@ -780,8 +783,8 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kv) ! XXX: Temporary de-scaling of N2_int(i,:) into a temporary variable - do k=1,G%ke+1 - N2_int_i(k) = US%s_to_T**2 * N2_int(i,k) + do K=1,G%ke+1 + N2_int_i(K) = US%s_to_T**2 * N2_int(i,K) enddo call CVMix_coeffs_tidal( Mdiff_out = Kv_tidal, & @@ -799,11 +802,15 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kv) do k=1,G%ke Kd_lay(i,j,k) = Kd_lay(i,j,k) + 0.5 * US%m2_s_to_Z2_T * (Kd_tidal(k) + Kd_tidal(k+1)) enddo - + if (present(Kd_int)) then + do K=1,G%ke+1 + Kd_int(i,j,K) = Kd_int(i,j,K) + (US%m2_s_to_Z2_T * Kd_tidal(K)) + enddo + endif ! Update viscosity with the proper unit conversion. if (associated(Kv)) then - do k=1,G%ke+1 - Kv(i,j,k) = Kv(i,j,k) + US%m2_s_to_Z2_T * Kv_tidal(k) ! Rescale from m2 s-1 to Z2 T-1. + do K=1,G%ke+1 + Kv(i,j,K) = Kv(i,j,K) + US%m2_s_to_Z2_T * Kv_tidal(K) ! Rescale from m2 s-1 to Z2 T-1. enddo endif @@ -895,11 +902,16 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kv) do k=1,G%ke Kd_lay(i,j,k) = Kd_lay(i,j,k) + 0.5 * US%m2_s_to_Z2_T * (Kd_tidal(k) + Kd_tidal(k+1)) enddo + if (present(Kd_int)) then + do K=1,G%ke+1 + Kd_int(i,j,K) = Kd_int(i,j,K) + (US%m2_s_to_Z2_T * Kd_tidal(K)) + enddo + endif ! Update viscosity if (associated(Kv)) then - do k=1,G%ke+1 - Kv(i,j,k) = Kv(i,j,k) + US%m2_s_to_Z2_T * Kv_tidal(k) ! Rescale from m2 s-1 to Z2 T-1. + do K=1,G%ke+1 + Kv(i,j,K) = Kv(i,j,K) + US%m2_s_to_Z2_T * Kv_tidal(K) ! Rescale from m2 s-1 to Z2 T-1. enddo endif diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index f244931376..73e4669734 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -23,8 +23,7 @@ module MOM_lateral_boundary_diffusion use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member - -use iso_fortran_env, only : stdout=>output_unit, stderr=>error_unit +use iso_fortran_env, only : stdout=>output_unit, stderr=>error_unit implicit none ; private @@ -38,15 +37,18 @@ module MOM_lateral_boundary_diffusion !> Sets parameters for lateral boundary mixing module. type, public :: lateral_boundary_diffusion_CS ; private - integer :: method !< Determine which of the three methods calculate - !! and apply near boundary layer fluxes - !! 1. Bulk-layer approach - !! 2. Along layer - integer :: deg !< Degree of polynomial reconstruction - integer :: surface_boundary_scheme !< Which boundary layer scheme to use - !! 1. ePBL; 2. KPP - logical :: limiter !< Controls wether a flux limiter is applied. - !! Only valid when method = 1. + integer :: method !< Determine which of the three methods calculate + !! and apply near boundary layer fluxes + !! 1. Along layer + !! 2. Bulk-layer approach (not recommended) + integer :: deg !< Degree of polynomial reconstruction + integer :: surface_boundary_scheme !< Which boundary layer scheme to use + !! 1. ePBL; 2. KPP + logical :: limiter !< Controls wether a flux limiter is applied. + !! Only valid when method = 2. + logical :: linear !< If True, apply a linear transition at the base/top of the boundary. + !! The flux will be fully applied at k=k_min and zero at k=k_max. + type(remapping_CS) :: remap_CS !< Control structure to hold remapping configuration type(KPP_CS), pointer :: KPP_CSp => NULL() !< KPP control structure needed to get BLD type(energetic_PBL_CS), pointer :: energetic_PBL_CSp => NULL() !< ePBL control structure needed to get BLD @@ -103,13 +105,16 @@ logical function lateral_boundary_diffusion_init(Time, G, param_file, diag, diab ! Read all relevant parameters and write them to the model log. call get_param(param_file, mdl, "LATERAL_BOUNDARY_METHOD", CS%method, & "Determine how to apply boundary lateral diffusion of tracers: \n"//& - "1. Bulk layer approach \n"//& - "2. Along layer approach", default=1) - if (CS%method == 1) then + "1. Along layer approach \n"//& + "2. Bulk layer approach (this option is not recommended)", default=1) + if (CS%method == 2) then call get_param(param_file, mdl, "APPLY_LIMITER", CS%limiter, & "If True, apply a flux limiter in the LBD. This is only available \n"//& - "when LATERAL_BOUNDARY_METHOD=1.", default=.false.) + "when LATERAL_BOUNDARY_METHOD=2.", default=.false.) endif + call get_param(param_file, mdl, "LBD_LINEAR_TRANSITION", CS%linear, & + "If True, apply a linear transition at the base/top of the boundary. \n"//& + "The flux will be fully applied at k=k_min and zero at k=k_max.", default=.false.) call get_param(param_file, mdl, "LBD_BOUNDARY_EXTRAP", boundary_extrap, & "Use boundary extrapolation in LBD code", & default=.false.) @@ -179,58 +184,63 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) call build_reconstructions_1d( CS%remap_CS, G%ke, h(i,j,:), tracer%t(i,j,:), ppoly0_coefs(i,j,:,:), & ppoly0_E(i,j,:,:), ppoly_S, remap_method, GV%H_subroundoff, GV%H_subroundoff) enddo ; enddo + ! Diffusive fluxes in the i-direction uFlx(:,:,:) = 0. vFlx(:,:,:) = 0. uFlx_bulk(:,:) = 0. vFlx_bulk(:,:) = 0. - ! Method #1 - if ( CS%method == 1 ) then + ! Method #1 (layer by layer) + if (CS%method == 1) then do j=G%jsc,G%jec do i=G%isc-1,G%iec if (G%mask2dCu(I,j)>0.) then - call fluxes_bulk_method(SURFACE, GV%ke, CS%deg, h(I,j,:), h(I+1,j,:), hbl(I,j), hbl(I+1,j), & - G%areaT(I,j), G%areaT(I+1,j), tracer%t(I,j,:), tracer%t(I+1,j,:), & - ppoly0_coefs(I,j,:,:), ppoly0_coefs(I+1,j,:,:), ppoly0_E(I,j,:,:), & - ppoly0_E(I+1,j,:,:), remap_method, Coef_x(I,j), uFlx_bulk(I,j), uFlx(I,j,:), CS%limiter) + call fluxes_layer_method(SURFACE, GV%ke, CS%deg, h(I,j,:), h(I+1,j,:), hbl(I,j), hbl(I+1,j), & + G%areaT(I,j), G%areaT(I+1,j), tracer%t(I,j,:), tracer%t(I+1,j,:), ppoly0_coefs(I,j,:,:), & + ppoly0_coefs(I+1,j,:,:), ppoly0_E(I,j,:,:), ppoly0_E(I+1,j,:,:), remap_method, Coef_x(I,j), & + uFlx(I,j,:), CS%linear) endif enddo enddo do J=G%jsc-1,G%jec do i=G%isc,G%iec if (G%mask2dCv(i,J)>0.) then - call fluxes_bulk_method(SURFACE, GV%ke, CS%deg, h(i,J,:), h(i,J+1,:), hbl(i,J), hbl(i,J+1), & - G%areaT(i,J), G%areaT(i,J+1), tracer%t(i,J,:), tracer%t(i,J+1,:), & - ppoly0_coefs(i,J,:,:), ppoly0_coefs(i,J+1,:,:), ppoly0_E(i,J,:,:), & - ppoly0_E(i,J+1,:,:), remap_method, Coef_y(i,J), vFlx_bulk(i,J), vFlx(i,J,:), CS%limiter) + call fluxes_layer_method(SURFACE, GV%ke, CS%deg, h(i,J,:), h(i,J+1,:), hbl(i,J), hbl(i,J+1), & + G%areaT(i,J), G%areaT(i,J+1), tracer%t(i,J,:), tracer%t(i,J+1,:), ppoly0_coefs(i,J,:,:), & + ppoly0_coefs(i,J+1,:,:), ppoly0_E(i,J,:,:), ppoly0_E(i,J+1,:,:), remap_method, Coef_y(i,J), & + vFlx(i,J,:), CS%linear) endif enddo enddo - ! Post tracer bulk diags - if (tracer%id_lbd_bulk_dfx>0) call post_data(tracer%id_lbd_bulk_dfx, uFlx_bulk*Idt, CS%diag) - if (tracer%id_lbd_bulk_dfy>0) call post_data(tracer%id_lbd_bulk_dfy, vFlx_bulk*Idt, CS%diag) - ! Method #2 + ! Method #2 (bulk approach) elseif (CS%method == 2) then do j=G%jsc,G%jec do i=G%isc-1,G%iec if (G%mask2dCu(I,j)>0.) then - call fluxes_layer_method(SURFACE, GV%ke, CS%deg, h(I,j,:), h(I+1,j,:), hbl(I,j), hbl(I+1,j), & - G%areaT(I,j), G%areaT(I+1,j), tracer%t(I,j,:), tracer%t(I+1,j,:), ppoly0_coefs(I,j,:,:), & - ppoly0_coefs(I+1,j,:,:), ppoly0_E(I,j,:,:), ppoly0_E(I+1,j,:,:), remap_method, Coef_x(I,j), uFlx(I,j,:)) + call fluxes_bulk_method(SURFACE, GV%ke, CS%deg, h(I,j,:), h(I+1,j,:), hbl(I,j), hbl(I+1,j), & + G%areaT(I,j), G%areaT(I+1,j), tracer%t(I,j,:), tracer%t(I+1,j,:), & + ppoly0_coefs(I,j,:,:), ppoly0_coefs(I+1,j,:,:), ppoly0_E(I,j,:,:), & + ppoly0_E(I+1,j,:,:), remap_method, Coef_x(I,j), uFlx_bulk(I,j), uFlx(I,j,:), CS%limiter, & + CS%linear) endif enddo enddo do J=G%jsc-1,G%jec do i=G%isc,G%iec if (G%mask2dCv(i,J)>0.) then - call fluxes_layer_method(SURFACE, GV%ke, CS%deg, h(i,J,:), h(i,J+1,:), hbl(i,J), hbl(i,J+1), & - G%areaT(i,J), G%areaT(i,J+1), tracer%t(i,J,:), tracer%t(i,J+1,:), ppoly0_coefs(i,J,:,:), & - ppoly0_coefs(i,J+1,:,:), ppoly0_E(i,J,:,:), ppoly0_E(i,J+1,:,:), remap_method, Coef_y(i,J), vFlx(i,J,:)) + call fluxes_bulk_method(SURFACE, GV%ke, CS%deg, h(i,J,:), h(i,J+1,:), hbl(i,J), hbl(i,J+1), & + G%areaT(i,J), G%areaT(i,J+1), tracer%t(i,J,:), tracer%t(i,J+1,:), & + ppoly0_coefs(i,J,:,:), ppoly0_coefs(i,J+1,:,:), ppoly0_E(i,J,:,:), & + ppoly0_E(i,J+1,:,:), remap_method, Coef_y(i,J), vFlx_bulk(i,J), vFlx(i,J,:), CS%limiter, & + CS%linear) endif enddo enddo + ! Post tracer bulk diags + if (tracer%id_lbd_bulk_dfx>0) call post_data(tracer%id_lbd_bulk_dfx, uFlx_bulk*Idt, CS%diag) + if (tracer%id_lbd_bulk_dfy>0) call post_data(tracer%id_lbd_bulk_dfy, vFlx_bulk*Idt, CS%diag) endif ! Update the tracer fluxes @@ -298,26 +308,26 @@ end subroutine lateral_boundary_diffusion !< Calculate bulk layer value of a scalar quantity as the thickness weighted average real function bulk_average(boundary, nk, deg, h, hBLT, phi, ppoly0_E, ppoly0_coefs, method, k_top, zeta_top, k_bot, & zeta_bot) - integer :: boundary !< SURFACE or BOTTOM [nondim] - integer :: nk !< Number of layers [nondim] - integer :: deg !< Degree of polynomial [nondim] - real, dimension(nk) :: h !< Layer thicknesses [H ~> m or kg m-2] - real :: hBLT !< Depth of the boundary layer [H ~> m or kg m-2] + integer :: boundary !< SURFACE or BOTTOM [nondim] + integer :: nk !< Number of layers [nondim] + integer :: deg !< Degree of polynomial [nondim] + real, dimension(nk) :: h !< Layer thicknesses [H ~> m or kg m-2] + real :: hBLT !< Depth of the boundary layer [H ~> m or kg m-2] real, dimension(nk) :: phi !< Scalar quantity - real, dimension(nk,2) :: ppoly0_E !< Edge value of polynomial - real, dimension(nk,deg+1) :: ppoly0_coefs !< Coefficients of polynomial - integer :: method !< Remapping scheme to use + real, dimension(nk,2) :: ppoly0_E !< Edge value of polynomial + real, dimension(nk,deg+1) :: ppoly0_coefs!< Coefficients of polynomial + integer :: method !< Remapping scheme to use integer :: k_top !< Index of the first layer within the boundary real :: zeta_top !< Fraction of the layer encompassed by the bottom boundary layer !! (0 if none, 1. if all). For the surface, this is always 0. because - !! integration starts at the surface [nondim] + !! integration starts at the surface [nondim] integer :: k_bot !< Index of the last layer within the boundary real :: zeta_bot !< Fraction of the layer encompassed by the surface boundary layer !! (0 if none, 1. if all). For the bottom boundary layer, this is always 1. - !! because integration starts at the bottom [nondim] + !! because integration starts at the bottom [nondim] ! Local variables - real :: htot !< Running sum of the thicknesses (top to bottom) [H ~> m or kg m-2] + real :: htot !< Running sum of the thicknesses (top to bottom) integer :: k !< k indice @@ -426,45 +436,50 @@ end subroutine boundary_k_range !> Calculate the lateral boundary diffusive fluxes using the layer by layer method. -!! See \ref section_method2 +!! See \ref section_method1 subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, & - ppoly0_coefs_L, ppoly0_coefs_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_layer) + ppoly0_coefs_L, ppoly0_coefs_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, & + F_layer, linear_decay) integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim] integer, intent(in ) :: nk !< Number of layers [nondim] integer, intent(in ) :: deg !< order of the polynomial reconstruction [nondim] - real, dimension(nk), intent(in ) :: h_L !< Layer thickness (left) [H ~> m or kg m-2] - real, dimension(nk), intent(in ) :: h_R !< Layer thickness (right) [H ~> m or kg m-2] + real, dimension(nk), intent(in ) :: h_L !< Layer thickness (left) [H ~> m or kg m-2] + real, dimension(nk), intent(in ) :: h_R !< Layer thickness (right) [H ~> m or kg m-2] real, intent(in ) :: hbl_L !< Thickness of the boundary boundary - !! layer (left) [H ~> m or kg m-2] + !! layer (left) [H ~> m or kg m-2] real, intent(in ) :: hbl_R !< Thickness of the boundary boundary - !! layer (right) [H ~> m or kg m-2] - real, intent(in ) :: area_L !< Area of the horizontal grid (left) [L2 ~> m2] - real, intent(in ) :: area_R !< Area of the horizontal grid (right) [L2 ~> m2] - real, dimension(nk), intent(in ) :: phi_L !< Tracer values (left) [conc] - real, dimension(nk), intent(in ) :: phi_R !< Tracer values (right) [conc] - real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_L !< Tracer reconstruction (left) [conc] - real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_R !< Tracer reconstruction (right) [conc] - real, dimension(nk,2), intent(in ) :: ppoly0_E_L !< Polynomial edge values (left) [ nondim ] - real, dimension(nk,2), intent(in ) :: ppoly0_E_R !< Polynomial edge values (right) [ nondim ] - integer, intent(in ) :: method !< Method of polynomial integration [ nondim ] + !! layer (right) [H ~> m or kg m-2] + real, intent(in ) :: area_L !< Area of the horizontal grid (left) [L2 ~> m2] + real, intent(in ) :: area_R !< Area of the horizontal grid (right) [L2 ~> m2] + real, dimension(nk), intent(in ) :: phi_L !< Tracer values (left) [conc] + real, dimension(nk), intent(in ) :: phi_R !< Tracer values (right) [conc] + real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_L !< Tracer reconstruction (left) [conc] + real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_R !< Tracer reconstruction (right) [conc] + real, dimension(nk,2), intent(in ) :: ppoly0_E_L !< Polynomial edge values (left) [nondim] + real, dimension(nk,2), intent(in ) :: ppoly0_E_R !< Polynomial edge values (right) [nondim] + integer, intent(in ) :: method !< Method of polynomial integration [nondim] real, intent(in ) :: khtr_u !< Horizontal diffusivities times delta t !! at a velocity point [L2 ~> m2] real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U- or V-point !! [H L2 conc ~> m3 conc] - + logical, optional, intent(in ) :: linear_decay !< If True, apply a linear transition at the base of + !! the boundary layer ! Local variables - real, dimension(nk) :: h_means !< Calculate the layer-wise harmonic means [H ~> m or kg m-2] - real :: khtr_avg !< Thickness-weighted diffusivity at the u-point [m^2 s^-1] + real, dimension(nk) :: h_means !< Calculate the layer-wise harmonic means [H ~> m or kg m-2] + real :: khtr_avg !< Thickness-weighted diffusivity at the u-point [m^2 s^-1] !! This is just to remind developers that khtr_avg should be !! computed once khtr is 3D. real :: heff !< Harmonic mean of layer thicknesses [H ~> m or kg m-2] real :: inv_heff !< Inverse of the harmonic mean of layer thicknesses !! [H-1 ~> m-1 or m2 kg-1] real :: phi_L_avg, phi_R_avg !< Bulk, thickness-weighted tracer averages (left and right column) - !! [conc m^-3 ] + !! [conc m^-3 ] real :: htot !< Total column thickness [H ~> m or kg m-2] - integer :: k, k_bot_min, k_top_max !< k-indices, min and max for top and bottom, respectively + real :: heff_tot !< Total effective column thickness in the transition layer [m] + integer :: k, k_bot_min, k_top_max !< k-indices, min and max for bottom and top, respectively + integer :: k_bot_max, k_top_min !< k-indices, max and min for bottom and top, respectively + integer :: k_bot_diff, k_top_diff !< different between left and right k-indices for bottom and top, respectively integer :: k_top_L, k_bot_L !< k-indices left integer :: k_top_R, k_bot_R !< k-indices right real :: zeta_top_L, zeta_top_R !< distance from the top of a layer to the boundary @@ -472,19 +487,30 @@ subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L real :: zeta_bot_L, zeta_bot_R !< distance from the bottom of a layer to the boundary !!layer depth [nondim] real :: h_work_L, h_work_R !< dummy variables - real :: hbl_min !< minimum BLD (left and right) [m] + real :: hbl_min !< minimum BLD (left and right) [m] + real :: wgt !< weight to be used in the linear transition to the interior [nondim] + real :: a !< coefficient to be used in the linear transition to the interior [nondim] + logical :: linear !< True if apply a linear transition F_layer(:) = 0.0 if (hbl_L == 0. .or. hbl_R == 0.) then return endif + linear = .false. + if (PRESENT(linear_decay)) then + linear = linear_decay + endif + ! Calculate vertical indices containing the boundary layer call boundary_k_range(boundary, nk, h_L, hbl_L, k_top_L, zeta_top_L, k_bot_L, zeta_bot_L) call boundary_k_range(boundary, nk, h_R, hbl_R, k_top_R, zeta_top_R, k_bot_R, zeta_bot_R) if (boundary == SURFACE) then k_bot_min = MIN(k_bot_L, k_bot_R) + k_bot_max = MAX(k_bot_L, k_bot_R) + k_bot_diff = (k_bot_max - k_bot_min) + ! make sure left and right k indices span same range if (k_bot_min .ne. k_bot_L) then k_bot_L = k_bot_min @@ -503,15 +529,37 @@ subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L heff = harmonic_mean(h_work_L, h_work_R) ! tracer flux where the minimum BLD intersets layer ! GMM, khtr_avg should be computed once khtr is 3D - F_layer(k_bot_min) = -(heff * khtr_u) * (phi_R_avg - phi_L_avg) + if ((linear) .and. (k_bot_diff .gt. 1)) then + ! apply linear decay at the base of hbl + do k = k_bot_min,1,-1 + heff = harmonic_mean(h_L(k), h_R(k)) + F_layer(k) = -(heff * khtr_u) * (phi_R(k) - phi_L(k)) + enddo + ! heff_total + heff_tot = 0.0 + do k = k_bot_min+1,k_bot_max, 1 + heff_tot = heff_tot + harmonic_mean(h_L(k), h_R(k)) + enddo - do k = k_bot_min-1,1,-1 - heff = harmonic_mean(h_L(k), h_R(k)) - F_layer(k) = -(heff * khtr_u) * (phi_R(k) - phi_L(k)) - enddo + a = -1.0/heff_tot + heff_tot = 0.0 + do k = k_bot_min+1,k_bot_max, 1 + heff = harmonic_mean(h_L(k), h_R(k)) + wgt = (a*(heff_tot + (heff * 0.5))) + 1.0 + F_layer(k) = -(heff * khtr_u) * (phi_R(k) - phi_L(k)) * wgt + heff_tot = heff_tot + heff + enddo + else + F_layer(k_bot_min) = -(heff * khtr_u) * (phi_R_avg - phi_L_avg) + do k = k_bot_min-1,1,-1 + heff = harmonic_mean(h_L(k), h_R(k)) + F_layer(k) = -(heff * khtr_u) * (phi_R(k) - phi_L(k)) + enddo + endif endif if (boundary == BOTTOM) then + ! TODO: GMM add option to apply linear decay k_top_max = MAX(k_top_L, k_top_R) ! make sure left and right k indices span same range if (k_top_max .ne. k_top_L) then @@ -542,28 +590,29 @@ subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L end subroutine fluxes_layer_method !> Apply the lateral boundary diffusive fluxes calculated from a 'bulk model' -!! See \ref section_method1 +!! See \ref section_method2 subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, ppoly0_coefs_L, & - ppoly0_coefs_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer, F_limit) + ppoly0_coefs_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer, F_limit, & + linear_decay) integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim] integer, intent(in ) :: nk !< Number of layers [nondim] integer, intent(in ) :: deg !< order of the polynomial reconstruction [nondim] - real, dimension(nk), intent(in ) :: h_L !< Layer thickness (left) [H ~> m or kg m-2] - real, dimension(nk), intent(in ) :: h_R !< Layer thickness (right) [H ~> m or kg m-2] + real, dimension(nk), intent(in ) :: h_L !< Layer thickness (left) [H ~> m or kg m-2] + real, dimension(nk), intent(in ) :: h_R !< Layer thickness (right) [H ~> m or kg m-2] real, intent(in ) :: hbl_L !< Thickness of the boundary boundary - !! layer (left) [H ~> m or kg m-2] + !! layer (left) [H ~> m or kg m-2] real, intent(in ) :: hbl_R !< Thickness of the boundary boundary - !! layer (left) [H ~> m or kg m-2] - real, intent(in ) :: area_L !< Area of the horizontal grid (left) [L2 ~> m2] - real, intent(in ) :: area_R !< Area of the horizontal grid (right) [L2 ~> m2] - real, dimension(nk), intent(in ) :: phi_L !< Tracer values (left) [conc] - real, dimension(nk), intent(in ) :: phi_R !< Tracer values (right) [conc] - real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_L !< Tracer reconstruction (left) [conc] - real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_R !< Tracer reconstruction (right) [conc] - real, dimension(nk,2), intent(in ) :: ppoly0_E_L !< Polynomial edge values (left) [nondim] - real, dimension(nk,2), intent(in ) :: ppoly0_E_R !< Polynomial edge values (right) [nondim] - integer, intent(in ) :: method !< Method of polynomial integration [nondim] + !! layer (left) [H ~> m or kg m-2] + real, intent(in ) :: area_L !< Area of the horizontal grid (left) [L2 ~> m2] + real, intent(in ) :: area_R !< Area of the horizontal grid (right) [L2 ~> m2] + real, dimension(nk), intent(in ) :: phi_L !< Tracer values (left) [conc] + real, dimension(nk), intent(in ) :: phi_R !< Tracer values (right) [conc] + real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_L !< Tracer reconstruction (left) [conc] + real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_R !< Tracer reconstruction (right) [conc] + real, dimension(nk,2), intent(in ) :: ppoly0_E_L !< Polynomial edge values (left) [nondim] + real, dimension(nk,2), intent(in ) :: ppoly0_E_R !< Polynomial edge values (right) [nondim] + integer, intent(in ) :: method !< Method of polynomial integration [nondim] real, intent(in ) :: khtr_u !< Horizontal diffusivities times delta t !! at a velocity point [L2 ~> m2] real, intent( out) :: F_bulk !< The bulk mixed layer lateral flux @@ -571,6 +620,8 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U- or V-point !! [H L2 conc ~> m3 conc] logical, optional, intent(in ) :: F_limit !< If True, apply a limiter + logical, optional, intent(in ) :: linear_decay !< If True, apply a linear transition at the base of + !! the boundary layer ! Local variables real, dimension(nk) :: h_means !< Calculate the layer-wise harmonic means [H ~> m or kg m-2] @@ -578,12 +629,14 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, !! This is just to remind developers that khtr_avg should be !! computed once khtr is 3D. real :: heff !< Harmonic mean of layer thicknesses [H ~> m or kg m-2] + real :: heff_tot !< Total effective column thickness in the transition layer [m] real :: inv_heff !< Inverse of the harmonic mean of layer thicknesses !! [H-1 ~> m-1 or m2 kg-1] real :: phi_L_avg, phi_R_avg !< Bulk, thickness-weighted tracer averages (left and right column) !! [conc m^-3 ] - real :: htot !< Total column thickness [H ~> m or kg m-2] + real :: htot ! Total column thickness [H ~> m or kg m-2] integer :: k, k_min, k_max !< k-indices, min and max for top and bottom, respectively + integer :: k_diff !< difference between k_max and k_min integer :: k_top_L, k_bot_L !< k-indices left integer :: k_top_R, k_bot_R !< k-indices right real :: zeta_top_L, zeta_top_R !< distance from the top of a layer to the @@ -594,12 +647,17 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, real :: F_max !< The maximum amount of flux that can leave a !! cell [m^3 conc] logical :: limiter !< True if flux limiter should be applied + logical :: linear !< True if apply a linear transition real :: hfrac !< Layer fraction wrt sum of all layers [nondim] real :: dphi !< tracer gradient [conc m^-3] + real :: wgt !< weight to be used in the linear transition to the + !! interior [nondim] + real :: a !< coefficient to be used in the linear transition to the + !! interior [nondim] + F_bulk = 0. + F_layer(:) = 0. if (hbl_L == 0. .or. hbl_R == 0.) then - F_bulk = 0. - F_layer(:) = 0. return endif @@ -607,6 +665,10 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, if (PRESENT(F_limit)) then limiter = F_limit endif + linear = .false. + if (PRESENT(linear_decay)) then + linear = linear_decay + endif ! Calculate vertical indices containing the boundary layer call boundary_k_range(boundary, nk, h_L, hbl_L, k_top_L, zeta_top_L, k_bot_L, zeta_bot_L) @@ -617,7 +679,6 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, zeta_top_L, k_bot_L, zeta_bot_L) phi_R_avg = bulk_average(boundary, nk, deg, h_R, hbl_R, phi_R, ppoly0_E_R, ppoly0_coefs_R, method, k_top_R, & zeta_top_R, k_bot_R, zeta_bot_R) - ! Calculate the 'bulk' diffusive flux from the bulk averaged quantities ! GMM, khtr_avg should be computed once khtr is 3D heff = harmonic_mean(hbl_L, hbl_R) @@ -625,31 +686,53 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, ! Calculate the layerwise sum of the vertical effective thickness. This is different than the heff calculated ! above, but is used as a way to decompose the fluxes onto the individual layers h_means(:) = 0. - if (boundary == SURFACE) then k_min = MIN(k_bot_L, k_bot_R) + k_max = MAX(k_bot_L, k_bot_R) + k_diff = (k_max - k_min) + if ((linear) .and. (k_diff .gt. 1)) then + do k=1,k_min + h_means(k) = harmonic_mean(h_L(k),h_R(k)) + enddo + ! heff_total + heff_tot = 0.0 + do k = k_min+1,k_max, 1 + heff_tot = heff_tot + harmonic_mean(h_L(k), h_R(k)) + enddo - ! left hand side - if (k_bot_L == k_min) then - h_work_L = h_L(k_min) * zeta_bot_L + a = -1.0/heff_tot + heff_tot = 0.0 + ! fluxes will decay linearly at base of hbl + do k = k_min+1,k_max, 1 + heff = harmonic_mean(h_L(k), h_R(k)) + wgt = (a*(heff_tot + (heff * 0.5))) + 1.0 + h_means(k) = harmonic_mean(h_L(k), h_R(k)) * wgt + heff_tot = heff_tot + heff + enddo else - h_work_L = h_L(k_min) - endif + ! left hand side + if (k_bot_L == k_min) then + h_work_L = h_L(k_min) * zeta_bot_L + else + h_work_L = h_L(k_min) + endif - ! right hand side - if (k_bot_R == k_min) then - h_work_R = h_R(k_min) * zeta_bot_R - else - h_work_R = h_R(k_min) - endif + ! right hand side + if (k_bot_R == k_min) then + h_work_R = h_R(k_min) * zeta_bot_R + else + h_work_R = h_R(k_min) + endif - h_means(k_min) = harmonic_mean(h_work_L,h_work_R) + h_means(k_min) = harmonic_mean(h_work_L,h_work_R) - do k=1,k_min-1 - h_means(k) = harmonic_mean(h_L(k),h_R(k)) - enddo + do k=1,k_min-1 + h_means(k) = harmonic_mean(h_L(k),h_R(k)) + enddo + endif elseif (boundary == BOTTOM) then + !TODO, GMM linear decay is not implemented here k_max = MAX(k_top_L, k_top_R) ! left hand side if (k_top_L == k_max) then @@ -672,14 +755,14 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, enddo endif - if ( SUM(h_means) == 0. ) then + if ( SUM(h_means) == 0. .or. F_bulk == 0.) then return - ! Decompose the bulk flux onto the individual layers + ! Decompose the bulk flux onto the individual layers else ! Initialize remaining thickness inv_heff = 1./SUM(h_means) do k=1,nk - if (h_means(k) > 0.) then + if ((h_means(k) > 0.) .and. (phi_L(k) /= phi_R(k))) then hfrac = h_means(k)*inv_heff F_layer(k) = F_bulk * hfrac @@ -1035,10 +1118,6 @@ logical function test_layer_fluxes(verbose, nk, test_name, F_calc, F_ans) test_layer_fluxes = .true. write(stdunit,*) "MOM_lateral_boundary_diffusion, UNIT TEST FAILED: ", test_name write(stdunit,10) k, F_calc(k), F_ans(k) - ! ### Once these unit tests are passing, and failures are caught properly, - ! we will post failure notifications to both stdout and stderr. - !write(stderr,*) "MOM_lateral_boundary_diffusion, UNIT TEST FAILED: ", test_name - !write(stderr,10) k, F_calc(k), F_ans(k) elseif (verbose) then write(stdunit,10) k, F_calc(k), F_ans(k) endif @@ -1096,12 +1175,37 @@ end function test_boundary_k_range !! !! Boundary lateral diffusion can be applied using one of the three methods: !! -!! * [Method #1: Bulk layer](@ref section_method1) (default); -!! * [Method #2: Along layer](@ref section_method2); +!! * [Method #1: Along layer](@ref section_method2) (default); +!! * [Method #2: Bulk layer](@ref section_method1); !! !! A brief summary of these methods is provided below. !! -!! \subsection section_method1 Bulk layer approach (Method #1) +!! \subsection section_method1 Along layer approach (Method #1) +!! +!! This is the recommended and more straight forward method where diffusion is +!! applied layer by layer using only information from neighboring cells. +!! +!! Step #1: compute vertical indices containing boundary layer (boundary_k_range). +!! For the TOP boundary layer, these are: +!! +!! k_top, k_bot, zeta_top, zeta_bot +!! +!! Step #2: calculate the diffusive flux at each layer: +!! +!! \f[ F_{k} = -KHTR \times h_{eff}(k) \times (\phi_R(k) - \phi_L(k)), \f] +!! where h_eff is the [harmonic mean](@ref section_harmonic_mean) of the layer thickness +!! in the left and right columns. This method does not require a limiter since KHTR +!! is already limted based on a diffusive CFL condition prior to the call of this +!! module. +!! +!! Step #3: option to linearly decay the flux from k_bot_min to k_bot_max: +!! +!! If LBD_LINEAR_TRANSITION = True and k_bot_diff > 1, the diffusive flux will decay +!! linearly between the top interface of the layer containing the minimum boundary +!! layer depth (k_bot_min) and the lower interface of the layer containing the +!! maximum layer depth (k_bot_max). +!! +!! \subsection section_method2 Bulk layer approach (Method #2) !! !! Apply the lateral boundary diffusive fluxes calculated from a 'bulk model'.This !! is a lower order representation (Kraus-Turner like approach) which assumes that @@ -1131,7 +1235,14 @@ end function test_boundary_k_range !! h_u is the [harmonic mean](@ref section_harmonic_mean) of thicknesses at each layer. !! Special care (layer reconstruction) must be taken at k_min = min(k_botL, k_bot_R). !! -!! Step #4: limit the tracer flux so that 1) only down-gradient fluxes are applied, +!! Step #4: option to linearly decay the flux from k_bot_min to k_bot_max: +!! +!! If LBD_LINEAR_TRANSITION = True and k_bot_diff > 1, the diffusive flux will decay +!! linearly between the top interface of the layer containing the minimum boundary +!! layer depth (k_bot_min) and the lower interface of the layer containing the +!! maximum layer depth (k_bot_max). +!! +!! Step #5: limit the tracer flux so that 1) only down-gradient fluxes are applied, !! and 2) the flux cannot be larger than F_max, which is defined using the tracer !! gradient: !! @@ -1142,25 +1253,6 @@ end function test_boundary_k_range !! 0 1 0 .2.2.2 !! 0 .2 !! -!! \subsection section_method2 Along layer approach (Method #2) -!! -!! This is a more straight forward method where diffusion is applied layer by layer using -!! only information from neighboring cells. -!! -!! Step #1: compute vertical indices containing boundary layer (boundary_k_range). -!! For the TOP boundary layer, these are: -!! -!! k_top, k_bot, zeta_top, zeta_bot -!! -!! Step #2: calculate the diffusive flux at each layer: -!! -!! \f[ F_{k} = -KHTR \times h_{eff}(k) \times (\phi_R(k) - \phi_L(k)), \f] -!! where h_eff is the [harmonic mean](@ref section_harmonic_mean) of the layer thickness -!! in the left and right columns. Special care (layer reconstruction) must be taken at -!! k_min = min(k_botL, k_bot_R). This method does not require a limiter since KHTR -!! is already limted based on a diffusive CFL condition prior to the call of this -!! module. -!! !! \subsection section_harmonic_mean Harmonic Mean !! !! The harmonic mean (HM) betwen h1 and h2 is defined as: diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index 30cdec3b37..d60aade72b 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -312,7 +312,7 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS, p_surf) pa_to_H = 1. / (GV%H_to_RZ * GV%g_Earth) k_top(:,:) = 1 ; k_bot(:,:) = 1 - zeta_top(:,:) = 0. ; zeta_bot(:,:) = 1. + zeta_top(:,:) = 0. ; zeta_bot(:,:) = 0. ! Check if hbl needs to be extracted if (CS%interior_only) then @@ -461,7 +461,7 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS, p_surf) CS%Pint(i,j,:), CS%Tint(i,j,:), CS%Sint(i,j,:), CS%dRdT(i,j,:), CS%dRdS(i,j,:), & CS%Pint(i+1,j,:), CS%Tint(i+1,j,:), CS%Sint(i+1,j,:), CS%dRdT(i+1,j,:), CS%dRdS(i+1,j,:), & CS%uPoL(I,j,:), CS%uPoR(I,j,:), CS%uKoL(I,j,:), CS%uKoR(I,j,:), CS%uhEff(I,j,:), & - k_bot(I,j), k_bot(I+1,j), 1.-zeta_bot(I,j), 1.-zeta_bot(I+1,j)) + k_bot(I,j), k_bot(I+1,j), zeta_bot(I,j), zeta_bot(I+1,j)) else call find_neutral_surface_positions_discontinuous(CS, G%ke, & CS%P_i(i,j,:,:), h(i,j,:), CS%T_i(i,j,:,:), CS%S_i(i,j,:,:), CS%ppoly_coeffs_T(i,j,:,:), & @@ -482,7 +482,7 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS, p_surf) CS%Pint(i,j,:), CS%Tint(i,j,:), CS%Sint(i,j,:), CS%dRdT(i,j,:), CS%dRdS(i,j,:), & CS%Pint(i,j+1,:), CS%Tint(i,j+1,:), CS%Sint(i,j+1,:), CS%dRdT(i,j+1,:), CS%dRdS(i,j+1,:), & CS%vPoL(i,J,:), CS%vPoR(i,J,:), CS%vKoL(i,J,:), CS%vKoR(i,J,:), CS%vhEff(i,J,:), & - k_bot(i,J), k_bot(i,J+1), 1.-zeta_bot(i,J), 1.-zeta_bot(i,J+1)) + k_bot(i,J), k_bot(i,J+1), zeta_bot(i,J), zeta_bot(i,J+1)) else call find_neutral_surface_positions_discontinuous(CS, G%ke, & CS%P_i(i,j,:,:), h(i,j,:), CS%T_i(i,j,:,:), CS%S_i(i,j,:,:), CS%ppoly_coeffs_T(i,j,:,:), &