diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index f4a39eeabb..8c39faf198 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -3515,13 +3515,13 @@ subroutine extract_surface_state(CS, sfc_state_in) enddo ! end of j loop endif ! melt_potential - if (allocated(sfc_state%taux_shelf) .and. associated(CS%visc%taux_shelf)) then + if (allocated(sfc_state%taux_shelf) .and. allocated(CS%visc%taux_shelf)) then !$OMP parallel do default(shared) do j=js,je ; do I=is-1,ie sfc_state%taux_shelf(I,j) = CS%visc%taux_shelf(I,j) enddo ; enddo endif - if (allocated(sfc_state%tauy_shelf) .and. associated(CS%visc%tauy_shelf)) then + if (allocated(sfc_state%tauy_shelf) .and. allocated(CS%visc%tauy_shelf)) then !$OMP parallel do default(shared) do J=js-1,je ; do i=is,ie sfc_state%tauy_shelf(i,J) = CS%visc%tauy_shelf(i,J) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 4adfa42673..c8fcfc52eb 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -220,39 +220,39 @@ module MOM_variables type, public :: vertvisc_type real :: Prandtl_turb !< The Prandtl number for the turbulent diffusion !! that is captured in Kd_shear [nondim]. - real, pointer, dimension(:,:) :: & - bbl_thick_u => NULL(), & !< The bottom boundary layer thickness at the u-points [Z ~> m]. - bbl_thick_v => NULL(), & !< The bottom boundary layer thickness at the v-points [Z ~> m]. - kv_bbl_u => NULL(), & !< The bottom boundary layer viscosity at the u-points [Z2 T-1 ~> m2 s-1]. - kv_bbl_v => NULL(), & !< The bottom boundary layer viscosity at the v-points [Z2 T-1 ~> m2 s-1]. - ustar_BBL => NULL() !< The turbulence velocity in the bottom boundary layer at h points [Z T-1 ~> m s-1]. - real, pointer, dimension(:,:) :: TKE_BBL => NULL() - !< A term related to the bottom boundary layer source of turbulent kinetic - !! energy, currently in [Z3 T-3 ~> m3 s-3], but may at some time be changed - !! to [R Z3 T-3 ~> W m-2]. - real, pointer, dimension(:,:) :: & - taux_shelf => NULL(), & !< The zonal stresses on the ocean under shelves [R Z L T-2 ~> Pa]. - tauy_shelf => NULL() !< The meridional stresses on the ocean under shelves [R Z L T-2 ~> Pa]. - real, pointer, dimension(:,:) :: tbl_thick_shelf_u => NULL() + real, allocatable, dimension(:,:) :: & + bbl_thick_u, & !< The bottom boundary layer thickness at the u-points [Z ~> m]. + bbl_thick_v, & !< The bottom boundary layer thickness at the v-points [Z ~> m]. + kv_bbl_u, & !< The bottom boundary layer viscosity at the u-points [Z2 T-1 ~> m2 s-1]. + kv_bbl_v, & !< The bottom boundary layer viscosity at the v-points [Z2 T-1 ~> m2 s-1]. + ustar_BBL, & !< The turbulence velocity in the bottom boundary layer at h points [Z T-1 ~> m s-1]. + TKE_BBL, & !< A term related to the bottom boundary layer source of turbulent kinetic + !! energy, currently in [Z3 T-3 ~> m3 s-3], but may at some time be changed + !! to [R Z3 T-3 ~> W m-2]. + taux_shelf, & !< The zonal stresses on the ocean under shelves [R Z L T-2 ~> Pa]. + tauy_shelf !< The meridional stresses on the ocean under shelves [R Z L T-2 ~> Pa]. + real, allocatable, dimension(:,:) :: tbl_thick_shelf_u !< Thickness of the viscous top boundary layer under ice shelves at u-points [Z ~> m]. - real, pointer, dimension(:,:) :: tbl_thick_shelf_v => NULL() + real, allocatable, dimension(:,:) :: tbl_thick_shelf_v !< Thickness of the viscous top boundary layer under ice shelves at v-points [Z ~> m]. - real, pointer, dimension(:,:) :: kv_tbl_shelf_u => NULL() + real, allocatable, dimension(:,:) :: kv_tbl_shelf_u !< Viscosity in the viscous top boundary layer under ice shelves at u-points [Z2 T-1 ~> m2 s-1]. - real, pointer, dimension(:,:) :: kv_tbl_shelf_v => NULL() + real, allocatable, dimension(:,:) :: kv_tbl_shelf_v !< Viscosity in the viscous top boundary layer under ice shelves at v-points [Z2 T-1 ~> m2 s-1]. - real, pointer, dimension(:,:) :: nkml_visc_u => NULL() + real, allocatable, dimension(:,:) :: nkml_visc_u !< The number of layers in the viscous surface mixed layer at u-points [nondim]. !! This is not an integer because there may be fractional layers, and it is stored in !! terms of layers, not depth, to facilitate the movement of the viscous boundary layer !! with the flow. - real, pointer, dimension(:,:) :: nkml_visc_v => NULL() + real, allocatable, dimension(:,:) :: nkml_visc_v !< The number of layers in the viscous surface mixed layer at v-points [nondim]. + real, allocatable, dimension(:,:,:) :: & + Ray_u, & !< The Rayleigh drag velocity to be applied to each layer at u-points [Z T-1 ~> m s-1]. + Ray_v !< The Rayleigh drag velocity to be applied to each layer at v-points [Z T-1 ~> m s-1]. + + ! The following elements are pointers so they can be used as targets for pointers in the restart registry. real, pointer, dimension(:,:) :: & - MLD => NULL() !< Instantaneous active mixing layer depth [Z ~> m]. - real, pointer, dimension(:,:,:) :: & - Ray_u => NULL(), & !< The Rayleigh drag velocity to be applied to each layer at u-points [Z T-1 ~> m s-1]. - Ray_v => NULL() !< The Rayleigh drag velocity to be applied to each layer at v-points [Z T-1 ~> m s-1]. + MLD => NULL() !< Instantaneous active mixing layer depth [Z ~> m]. real, pointer, dimension(:,:,:) :: Kd_shear => NULL() !< The shear-driven turbulent diapycnal diffusivity at the interfaces between layers !! in tracer columns [Z2 T-1 ~> m2 s-1]. diff --git a/src/diagnostics/MOM_PointAccel.F90 b/src/diagnostics/MOM_PointAccel.F90 index 423ef2b4f9..d53b2e6636 100644 --- a/src/diagnostics/MOM_PointAccel.F90 +++ b/src/diagnostics/MOM_PointAccel.F90 @@ -39,6 +39,8 @@ module MOM_PointAccel !! written by this PE during the current run. integer :: max_writes !< The maximum number of times any PE can write out !! a column's worth of accelerations during a run. + logical :: full_column !< If true, write out the accelerations in all massive layers, + !! otherwise just document the ones with large velocities. type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock. type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to !! regulate the timing of diagnostic output. @@ -80,11 +82,12 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st !! call to PointAccel_init. real, intent(in) :: vel_rpt !< The velocity magnitude that triggers a report [L T-1 ~> m s-1] real, optional, intent(in) :: str !< The surface wind stress [R L Z T-2 ~> Pa] - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), & optional, intent(in) :: a !< The layer coupling coefficients from vertvisc [Z T-1 ~> m s-1]. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & optional, intent(in) :: hv !< The layer thicknesses at velocity grid points, !! from vertvisc [H ~> m or kg m-2]. + ! Local variables real :: CFL ! The local velocity-based CFL number [nondim] real :: Angstrom ! A negligibly small thickness [H ~> m or kg m-2] @@ -145,6 +148,9 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st if (ke < ks) then ks = 1; ke = nz; write(file,'("U: Unable to set ks & ke.")') endif + if (CS%full_column) then + ks = 1 ; ke = nz + endif call get_date(CS%Time, yr, mo, day, hr, minute, sec) call get_time((CS%Time - set_date(yr, 1, 1, 0, 0, 0)), sec, yearday) @@ -217,21 +223,23 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st (vel_scale*ADp%du_other(I,j,k)) ; enddo endif if (present(a)) then - write(file,'(/,"a: ")', advance='no') - do k=ks,ke+1 ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (US%Z_to_m*a(I,j,k)*dt) ; enddo + write(file,'(/,"a: ",ES10.3," ")', advance='no') US%Z_to_m*a(I,j,ks)*dt + do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') (US%Z_to_m*a(I,j,K)*dt) ; enddo endif if (present(hv)) then write(file,'(/,"hvel: ")', advance='no') do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') h_scale*hv(I,j,k) ; enddo endif - write(file,'(/,"Stress: ",ES10.3)') vel_scale*US%Z_to_m * (str*dt / GV%Rho0) + if (present(str)) then + write(file,'(/,"Stress: ",ES10.3)', advance='no') vel_scale*US%Z_to_m * (str*dt / GV%Rho0) + endif if (associated(CS%u_accel_bt)) then - write(file,'("dubt: ")', advance='no') + write(file,'(/,"dubt: ")', advance='no') do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (vel_scale*dt*CS%u_accel_bt(I,j,k)) ; enddo - write(file,'(/)') endif + write(file,'(/)') write(file,'(/,"h--: ")', advance='no') do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (h_scale*hin(i,j-1,k)) ; enddo @@ -249,14 +257,12 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st e(nz+1) = -US%Z_to_m*(G%bathyT(i,j) + G%Z_ref) do k=nz,1,-1 ; e(K) = e(K+1) + h_scale*hin(i,j,k) ; enddo - write(file,'(/,"e-: ")', advance='no') - write(file,'(ES10.3," ")', advance='no') e(ks) + write(file,'(/,"e-: ",ES10.3," ")', advance='no') e(ks) do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') e(K) ; enddo e(nz+1) = -US%Z_to_m*(G%bathyT(i+1,j) + G%Z_ref) do k=nz,1,-1 ; e(K) = e(K+1) + h_scale*hin(i+1,j,k) ; enddo - write(file,'(/,"e+: ")', advance='no') - write(file,'(ES10.3," ")', advance='no') e(ks) + write(file,'(/,"e+: ",ES10.3," ")', advance='no') e(ks) do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') e(K) ; enddo if (associated(CS%T)) then write(file,'(/,"T-: ")', advance='no') @@ -415,11 +421,12 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st !! call to PointAccel_init. real, intent(in) :: vel_rpt !< The velocity magnitude that triggers a report [L T-1 ~> m s-1] real, optional, intent(in) :: str !< The surface wind stress [R L Z T-2 ~> Pa] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), & optional, intent(in) :: a !< The layer coupling coefficients from vertvisc [Z T-1 ~> m s-1]. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & optional, intent(in) :: hv !< The layer thicknesses at velocity grid points, !! from vertvisc [H ~> m or kg m-2]. + ! Local variables real :: CFL ! The local velocity-based CFL number [nondim] real :: Angstrom ! A negligibly small thickness [H ~> m or kg m-2] @@ -479,6 +486,9 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st if (ke < ks) then ks = 1; ke = nz; write(file,'("V: Unable to set ks & ke.")') endif + if (CS%full_column) then + ks = 1 ; ke = nz + endif call get_date(CS%Time, yr, mo, day, hr, minute, sec) call get_time((CS%Time - set_date(yr, 1, 1, 0, 0, 0)), sec, yearday) @@ -556,21 +566,23 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st (vel_scale*ADp%dv_other(i,J,k)) ; enddo endif if (present(a)) then - write(file,'(/,"a: ")', advance='no') - do k=ks,ke+1 ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (US%Z_to_m*a(i,j,k)*dt) ; enddo + write(file,'(/,"a: ",ES10.3," ")', advance='no') US%Z_to_m*a(i,J,ks)*dt + do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') (US%Z_to_m*a(i,J,K)*dt) ; enddo endif if (present(hv)) then write(file,'(/,"hvel: ")', advance='no') do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') h_scale*hv(i,J,k) ; enddo endif - write(file,'(/,"Stress: ",ES10.3)') vel_scale*US%Z_to_m * (str*dt / GV%Rho0) + if (present(str)) then + write(file,'(/,"Stress: ",ES10.3)', advance='no') vel_scale*US%Z_to_m * (str*dt / GV%Rho0) + endif if (associated(CS%v_accel_bt)) then write(file,'("dvbt: ")', advance='no') do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (vel_scale*dt*CS%v_accel_bt(i,J,k)) ; enddo - write(file,'(/)') endif + write(file,'(/)') write(file,'("h--: ")', advance='no') do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') h_scale*hin(i-1,j,k) ; enddo @@ -587,14 +599,12 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st e(nz+1) = -US%Z_to_m*(G%bathyT(i,j) + G%Z_ref) do k=nz,1,-1 ; e(K) = e(K+1) + h_scale*hin(i,j,k) ; enddo - write(file,'(/,"e-: ")', advance='no') - write(file,'(ES10.3," ")', advance='no') e(ks) + write(file,'(/,"e-: ",ES10.3," ")', advance='no') e(ks) do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') e(K) ; enddo e(nz+1) = -US%Z_to_m*(G%bathyT(i,j+1) + G%Z_ref) do k=nz,1,-1 ; e(K) = e(K+1) + h_scale*hin(i,j+1,k) ; enddo - write(file,'(/,"e+: ")', advance='no') - write(file,'(ES10.3," ")', advance='no') e(ks) + write(file,'(/,"e+: ",ES10.3," ")', advance='no') e(ks) do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') e(K) ; enddo if (associated(CS%T)) then write(file,'(/,"T-: ")', advance='no') @@ -773,6 +783,10 @@ subroutine PointAccel_init(MIS, Time, G, param_file, diag, dirs, CS) call get_param(param_file, mdl, "MAX_TRUNC_FILE_SIZE_PER_PE", CS%max_writes, & "The maximum number of columns of truncations that any PE "//& "will write out during a run.", default=50, debuggingParam=.true.) + call get_param(param_file, mdl, "DEBUG_FULL_COLUMN", CS%full_column, & + "If true, write out the accelerations in all massive layers; otherwise "//& + "just document the ones with large velocities.", & + default=.false., debuggingParam=.true.) if (len_trim(dirs%output_directory) > 0) then if (len_trim(CS%u_trunc_file) > 0) & diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 1cd20d3c96..2661251766 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -248,7 +248,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h ! Calculate drag_rate_visc(i,j) which accounts for the model bottom mean flow - if (CS%visc_drag) then + if (CS%visc_drag .and. allocated(visc%Kv_bbl_u) .and. allocated(visc%Kv_bbl_v)) then !$OMP parallel do default(shared) do j=js,je ; do I=is-1,ie drag_vel_u(I,j) = 0.0 diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 20f92af86b..eff9d7ff72 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -588,19 +588,19 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i call hchksum(Kd_extra_S, "MOM_set_diffusivity: Kd_extra_S", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s) endif - if (associated(visc%kv_bbl_u) .and. associated(visc%kv_bbl_v)) then + if (allocated(visc%kv_bbl_u) .and. allocated(visc%kv_bbl_v)) then call uvchksum("BBL Kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, G%HI, & haloshift=0, symmetric=.true., scale=US%Z2_T_to_m2_s, & scalar_pair=.true.) endif - if (associated(visc%bbl_thick_u) .and. associated(visc%bbl_thick_v)) then + if (allocated(visc%bbl_thick_u) .and. allocated(visc%bbl_thick_v)) then call uvchksum("BBL bbl_thick_[uv]", visc%bbl_thick_u, visc%bbl_thick_v, & G%HI, haloshift=0, symmetric=.true., scale=US%Z_to_m, & scalar_pair=.true.) endif - if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) then + if (allocated(visc%Ray_u) .and. allocated(visc%Ray_v)) then call uvchksum("Ray_[uv]", visc%Ray_u, visc%Ray_v, G%HI, 0, symmetric=.true., scale=US%Z_to_m*US%s_to_T) endif @@ -1198,7 +1198,7 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & cdrag_sqrt = sqrt(CS%cdrag) TKE_Ray = 0.0 ; Rayleigh_drag = .false. - if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) Rayleigh_drag = .true. + if (allocated(visc%Ray_u) .and. allocated(visc%Ray_v)) Rayleigh_drag = .true. I_Rho0 = 1.0 / (GV%Rho0) R0_g = GV%Rho0 / (US%L_to_Z**2 * GV%g_Earth) @@ -1416,7 +1416,7 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Kd_int ! Determine whether to add Rayleigh drag contribution to TKE Rayleigh_drag = .false. - if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) Rayleigh_drag = .true. + if (allocated(visc%Ray_u) .and. allocated(visc%Ray_v)) Rayleigh_drag = .true. I_Rho0 = 1.0 / (GV%Rho0) cdrag_sqrt = sqrt(CS%cdrag) @@ -1668,7 +1668,7 @@ subroutine set_BBL_TKE(u, v, h, fluxes, visc, G, GV, US, CS, OBC) real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] type(forcing), intent(in) :: fluxes !< A structure of thermodynamic surface fluxes - type(vertvisc_type), intent(in) :: visc !< Structure containing vertical viscosities, bottom + type(vertvisc_type), intent(inout) :: visc !< Structure containing vertical viscosities, bottom !! boundary layer properties and related fields. type(set_diffusivity_CS), pointer :: CS !< Diffusivity control structure type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. @@ -1692,6 +1692,7 @@ subroutine set_BBL_TKE(u, v, h, fluxes, visc, G, GV, US, CS, OBC) v2_bbl ! square of average meridional velocity in BBL [L2 T-2 ~> m2 s-2] real :: cdrag_sqrt ! square root of the drag coefficient [nondim] + real :: I_cdrag_sqrt ! The inverse of the square root of the drag coefficient [nondim] real :: hvel ! thickness at velocity points [Z ~> m]. logical :: domore, do_i(SZI_(G)) @@ -1716,29 +1717,34 @@ subroutine set_BBL_TKE(u, v, h, fluxes, visc, G, GV, US, CS, OBC) "Module must be initialized before it is used.") if (.not.CS%bottomdraglaw .or. (CS%BBL_effic<=0.0)) then - if (associated(visc%ustar_BBL)) then + if (allocated(visc%ustar_BBL)) then do j=js,je ; do i=is,ie ; visc%ustar_BBL(i,j) = 0.0 ; enddo ; enddo endif - if (associated(visc%TKE_BBL)) then + if (allocated(visc%TKE_BBL)) then do j=js,je ; do i=is,ie ; visc%TKE_BBL(i,j) = 0.0 ; enddo ; enddo endif return endif cdrag_sqrt = sqrt(CS%cdrag) + I_cdrag_sqrt = 0.0 ; if (cdrag_sqrt > 0.0) I_cdrag_sqrt = 1.0 / cdrag_sqrt !$OMP parallel default(shared) private(do_i,vhtot,htot,domore,hvel,uhtot,ustar,u2_bbl) !$OMP do do J=js-1,je - ! Determine ustar and the square magnitude of the velocity in the - ! bottom boundary layer. Together these give the TKE source and - ! vertical decay scale. - do i=is,ie ; if ((G%mask2dCv(i,J) > 0.0) .and. (cdrag_sqrt*visc%bbl_thick_v(i,J) > 0.0)) then - do_i(i) = .true. ; vhtot(i) = 0.0 ; htot(i) = 0.0 - vstar(i,J) = visc%Kv_bbl_v(i,J) / (cdrag_sqrt*visc%bbl_thick_v(i,J)) - else - do_i(i) = .false. ; vstar(i,J) = 0.0 ; htot(i) = 0.0 - endif ; enddo + ! Determine ustar and the square magnitude of the velocity in the bottom boundary layer. + ! Together these give the TKE source and vertical decay scale. + do i=is,ie + do_i(i) = .false. ; vstar(i,J) = 0.0 ; vhtot(i) = 0.0 ; htot(i) = 0.0 + enddo + if (allocated(visc%Kv_bbl_v)) then + do i=is,ie ; if ((G%mask2dCv(i,J) > 0.0) .and. (cdrag_sqrt*visc%bbl_thick_v(i,J) > 0.0)) then + do_i(i) = .true. + vstar(i,J) = visc%Kv_bbl_v(i,J) / (cdrag_sqrt*visc%bbl_thick_v(i,J)) + endif ; enddo + endif + !### What about terms from visc%Ray? + do k=nz,1,-1 domore = .false. do i=is,ie ; if (do_i(i)) then @@ -1782,12 +1788,16 @@ subroutine set_BBL_TKE(u, v, h, fluxes, visc, G, GV, US, CS, OBC) enddo !$OMP do do j=js,je - do I=is-1,ie ; if ((G%mask2dCu(I,j) > 0.0) .and. (cdrag_sqrt*visc%bbl_thick_u(I,j) > 0.0)) then - do_i(I) = .true. ; uhtot(I) = 0.0 ; htot(I) = 0.0 - ustar(I) = visc%Kv_bbl_u(I,j) / (cdrag_sqrt*visc%bbl_thick_u(I,j)) - else - do_i(I) = .false. ; ustar(I) = 0.0 ; htot(I) = 0.0 - endif ; enddo + do I=is-1,ie + do_i(I) = .false. ; ustar(I) = 0.0 ; uhtot(I) = 0.0 ; htot(I) = 0.0 + enddo + if (allocated(visc%bbl_thick_u)) then + do I=is-1,ie ; if ((G%mask2dCu(I,j) > 0.0) .and. (cdrag_sqrt*visc%bbl_thick_u(I,j) > 0.0)) then + do_i(I) = .true. + ustar(I) = visc%Kv_bbl_u(I,j) / (cdrag_sqrt*visc%bbl_thick_u(I,j)) + endif ; enddo + endif + do k=nz,1,-1 ; domore = .false. do I=is-1,ie ; if (do_i(I)) then ! Determine if grid point is an OBC diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 7934b6b019..22d65110be 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -67,6 +67,9 @@ module MOM_set_visc !! actual velocity in the bottommost `HBBL`, depending !! on whether linear_drag is true. !! Runtime parameter `BOTTOMDRAGLAW`. + logical :: body_force_drag !< If true, the bottom stress is imposed as an explicit body force + !! applied over a fixed distance from the bottom, rather than as an + !! implicit calculation based on an enhanced near-bottom viscosity. logical :: BBL_use_EOS !< If true, use the equation of state in determining !! the properties of the bottom boundary layer. logical :: linear_drag !< If true, the drag law is cdrag*`DRAG_BG_VEL`*u. @@ -128,7 +131,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any !! available thermodynamic fields. Absent fields - !! have NULL ptrs.. + !! have NULL ptrs. type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and !! related fields. type(set_visc_CS), intent(inout) :: CS !< The control structure returned by a previous @@ -146,7 +149,9 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) ! layer with temperature [R C-1 ~> kg m-3 degC-1]. dR_dS, & ! Partial derivative of the density in the bottom boundary ! layer with salinity [R S-1 ~> kg m-3 ppt-1]. - press ! The pressure at which dR_dT and dR_dS are evaluated [R L2 T-2 ~> Pa]. + press, & ! The pressure at which dR_dT and dR_dS are evaluated [R L2 T-2 ~> Pa]. + umag_avg, & ! The average magnitude of velocities in the bottom boundary layer [L T-1 ~> m s-1]. + h_bbl_drag ! The thickness over which to apply drag as a body force [H ~> m or kg m-2]. real :: htot ! Sum of the layer thicknesses up to some point [H ~> m or kg m-2]. real :: htot_vel ! Sum of the layer thicknesses up to some point [H ~> m or kg m-2]. @@ -199,8 +204,9 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) ! quadratic bottom drag [L2 T-2 ~> m2 s-2]. real :: hwtot ! Sum of the thicknesses used to calculate ! the near-bottom velocity magnitude [H ~> m or kg m-2]. - real :: hutot ! Running sum of thicknesses times the - ! velocity magnitudes [H L T-1 ~> m2 s-1 or kg m-1 s-1]. + real :: I_hwtot ! The Adcroft reciprocal of hwtot [H-1 ~> m-1 or m2 kg-1]. + real :: hutot ! Running sum of thicknesses times the velocity + ! magnitudes [H L T-1 ~> m2 s-1 or kg m-1 s-1]. real :: Thtot ! Running sum of thickness times temperature [C H ~> degC m or degC kg m-2]. real :: Shtot ! Running sum of thickness times salinity [S H ~> ppt m or ppt kg m-2]. real :: hweight ! The thickness of a layer that is within Hbbl @@ -265,6 +271,9 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) ! viscous bottom boundary layer [nondim]. real :: BBL_visc_frac ! The fraction of all the drag that is expressed as ! a viscous bottom boundary layer [nondim]. + real :: h_bbl_fr ! The fraction of the bottom boundary layer in a layer [nondim]. + real :: h_sum ! The sum of the thicknesses of the layers below the one being + ! worked on [H ~> m or kg m-2]. real, parameter :: C1_3 = 1.0/3.0, C1_6 = 1.0/6.0, C1_12 = 1.0/12.0 real :: C2pi_3 ! An irrational constant, 2/3 pi. real :: tmp ! A temporary variable. @@ -373,6 +382,9 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) if (.not.use_BBL_EOS) Rml_vel(:,:) = 0.0 + if (allocated(visc%Ray_u)) visc%Ray_u(:,:,:) = 0.0 + if (allocated(visc%Ray_v)) visc%Ray_v(:,:,:) = 0.0 + !$OMP parallel do default(private) shared(u,v,h,tv,visc,G,GV,US,CS,Rml,nz,nkmb, & !$OMP nkml,Isq,Ieq,Jsq,Jeq,h_neglect,Rho0x400_G,C2pi_3, & !$OMP U_bg_sq,cdrag_sqrt_Z,cdrag_sqrt,K2,use_BBL_EOS, & @@ -505,7 +517,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) endif endif ; endif - if (use_BBL_EOS .or. .not.CS%linear_drag) then + if (use_BBL_EOS .or. CS%body_force_drag .or. .not.CS%linear_drag) then ! Calculate the mean velocity magnitude over the bottommost CS%Hbbl of ! the water column for determining the quadratic bottom drag. ! Used in ustar(i) @@ -528,16 +540,14 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) U_bg_sq = 0.5*( G%mask2dT(i,j)*(CS%tideamp(i,j)*CS%tideamp(i,j))+ & G%mask2dT(i+1,j)*(CS%tideamp(i+1,j)*CS%tideamp(i+1,j)) ) endif - hutot = hutot + hweight * sqrt(u(I,j,k)*u(I,j,k) + & - v_at_u*v_at_u + U_bg_sq) + hutot = hutot + hweight * sqrt(u(I,j,k)*u(I,j,k) + v_at_u*v_at_u + U_bg_sq) else u_at_v = set_u_at_v(u, h, G, GV, i, j, k, mask_u, OBC) if (CS%BBL_use_tidal_bg) then U_bg_sq = 0.5*( G%mask2dT(i,j)*(CS%tideamp(i,j)*CS%tideamp(i,j))+ & G%mask2dT(i,j+1)*(CS%tideamp(i,j+1)*CS%tideamp(i,j+1)) ) endif - hutot = hutot + hweight * sqrt(v(i,J,k)*v(i,J,k) + & - u_at_v*u_at_v + U_bg_sq) + hutot = hutot + hweight * sqrt(v(i,J,k)*v(i,J,k) + u_at_v*u_at_v + U_bg_sq) endif ; endif if (use_BBL_EOS .and. (hweight >= 0.0)) then @@ -548,11 +558,17 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) ! Set u* based on u*^2 = Cdrag u_bbl^2 if (.not.CS%linear_drag .and. (hwtot > 0.0)) then - ustar(i) = cdrag_sqrt_Z*hutot/hwtot + ustar(i) = cdrag_sqrt_Z*hutot / hwtot else ustar(i) = cdrag_sqrt_Z*CS%drag_bg_vel endif + ! Find the Adcroft reciprocal of the total thickness weights + I_hwtot = 0.0 ; if (hwtot > 0.0) I_hwtot = 1.0 / hwtot + + umag_avg(i) = hutot * I_hwtot + h_bbl_drag(i) = hwtot + if (use_BBL_EOS) then ; if (hwtot > 0.0) then T_EOS(i) = Thtot/hwtot ; S_EOS(i) = Shtot/hwtot else @@ -931,14 +947,12 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) if (m==1) then if (Rayleigh > 0.0) then v_at_u = set_v_at_u(v, h, G, GV, i, j, k, mask_v, OBC) - visc%Ray_u(I,j,k) = Rayleigh*sqrt(u(I,j,k)*u(I,j,k) + & - v_at_u*v_at_u + U_bg_sq) + visc%Ray_u(I,j,k) = Rayleigh * sqrt(u(I,j,k)*u(I,j,k) + v_at_u*v_at_u + U_bg_sq) else ; visc%Ray_u(I,j,k) = 0.0 ; endif else if (Rayleigh > 0.0) then u_at_v = set_u_at_v(u, h, G, GV, i, j, k, mask_u, OBC) - visc%Ray_v(i,J,k) = Rayleigh*sqrt(v(i,J,k)*v(i,J,k) + & - u_at_v*u_at_v + U_bg_sq) + visc%Ray_v(i,J,k) = Rayleigh * sqrt(v(i,J,k)*v(i,J,k) + u_at_v*u_at_v + U_bg_sq) else ; visc%Ray_v(i,J,k) = 0.0 ; endif endif @@ -995,13 +1009,32 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) kv_bbl = cdrag_sqrt*ustar(i)*bbl_thick_Z endif endif + + if (CS%body_force_drag .and. (h_bbl_drag(i) > 0.0)) then + ! Increment the Rayleigh drag as a way introduce the bottom drag as a body force. + h_sum = 0.0 + I_hwtot = 1.0 / h_bbl_drag(i) + do k=nz,1,-1 + h_bbl_fr = min(h_bbl_drag(i) - h_sum, h_at_vel(i,k)) * I_hwtot + if (m==1) then + visc%Ray_u(I,j,k) = visc%Ray_u(I,j,k) + (CS%cdrag*US%L_to_Z*umag_avg(I)) * h_bbl_fr + else + visc%Ray_v(i,J,k) = visc%Ray_v(i,J,k) + (CS%cdrag*US%L_to_Z*umag_avg(i)) * h_bbl_fr + endif + h_sum = h_sum + h_at_vel(i,k) + if (h_sum >= bbl_thick) exit ! The top of this layer is above the drag zone. + enddo + ! Do not enhance the near-bottom viscosity in this case. + Kv_bbl = CS%Kv_BBL_min + endif + kv_bbl = max(CS%Kv_BBL_min, kv_bbl) if (m==1) then - visc%Kv_bbl_u(I,j) = kv_bbl visc%bbl_thick_u(I,j) = bbl_thick_Z + if (allocated(visc%Kv_bbl_u)) visc%Kv_bbl_u(I,j) = kv_bbl else - visc%Kv_bbl_v(i,J) = kv_bbl visc%bbl_thick_v(i,J) = bbl_thick_Z + if (allocated(visc%Kv_bbl_v)) visc%Kv_bbl_v(i,J) = kv_bbl endif endif ; enddo ! end of i loop enddo ; enddo ! end of m & j loops @@ -1025,12 +1058,12 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) call post_data(CS%id_Ray_v, visc%Ray_v, CS%diag) if (CS%debug) then - if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) & + if (allocated(visc%Ray_u) .and. allocated(visc%Ray_v)) & call uvchksum("Ray [uv]", visc%Ray_u, visc%Ray_v, G%HI, haloshift=0, scale=US%Z_to_m*US%s_to_T) - if (associated(visc%kv_bbl_u) .and. associated(visc%kv_bbl_v)) & + if (allocated(visc%kv_bbl_u) .and. allocated(visc%kv_bbl_v)) & call uvchksum("kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, G%HI, & haloshift=0, scale=US%Z2_T_to_m2_s, scalar_pair=.true.) - if (associated(visc%bbl_thick_u) .and. associated(visc%bbl_thick_v)) & + if (allocated(visc%bbl_thick_u) .and. allocated(visc%bbl_thick_v)) & call uvchksum("bbl_thick_[uv]", visc%bbl_thick_u, visc%bbl_thick_v, & G%HI, haloshift=0, scale=US%Z_to_m, scalar_pair=.true.) endif @@ -1193,8 +1226,8 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS) real :: htot_vel ! Sum of the layer thicknesses up to some point [H ~> m or kg m-2]. real :: hwtot ! Sum of the thicknesses used to calculate ! the near-bottom velocity magnitude [H ~> m or kg m-2]. - real :: hutot ! Running sum of thicknesses times the - ! velocity magnitudes [H L T-1 ~> m2 s-1 or kg m-1 s-1]. + real :: hutot ! Running sum of thicknesses times the velocity + ! magnitudes [H L T-1 ~> m2 s-1 or kg m-1 s-1]. real :: hweight ! The thickness of a layer that is within Hbbl ! of the bottom [H ~> m or kg m-2]. real :: tbl_thick_Z ! The thickness of the top boundary layer [Z ~> m]. @@ -1277,13 +1310,18 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS) if (associated(forces%frac_shelf_u)) then ! This configuration has ice shelves, and the appropriate variables need to be ! allocated. If the arrays have already been allocated, these calls do nothing. - call safe_alloc_ptr(visc%tauy_shelf, G%isd, G%ied, G%JsdB, G%JedB) - call safe_alloc_ptr(visc%tbl_thick_shelf_u, G%IsdB, G%IedB, G%jsd, G%jed) - call safe_alloc_ptr(visc%tbl_thick_shelf_v, G%isd, G%ied, G%JsdB, G%JedB) - call safe_alloc_ptr(visc%kv_tbl_shelf_u, G%IsdB, G%IedB, G%jsd, G%jed) - call safe_alloc_ptr(visc%kv_tbl_shelf_v, G%isd, G%ied, G%JsdB, G%JedB) - call safe_alloc_ptr(visc%taux_shelf, G%IsdB, G%IedB, G%jsd, G%jed) - call safe_alloc_ptr(visc%tauy_shelf, G%isd, G%ied, G%JsdB, G%JedB) + if (.not.allocated(visc%taux_shelf)) & + allocate(visc%taux_shelf(G%IsdB:G%IedB, G%jsd:G%jed), source=0.0) + if (.not.allocated(visc%tauy_shelf)) & + allocate(visc%tauy_shelf(G%isd:G%ied, G%JsdB:G%JedB), source=0.0) + if (.not.allocated(visc%tbl_thick_shelf_u)) & + allocate(visc%tbl_thick_shelf_u(G%IsdB:G%IedB, G%jsd:G%jed), source=0.0) + if (.not.allocated(visc%tbl_thick_shelf_v)) & + allocate(visc%tbl_thick_shelf_v(G%isd:G%ied, G%JsdB:G%JedB), source=0.0) + if (.not.allocated(visc%kv_tbl_shelf_u)) & + allocate(visc%kv_tbl_shelf_u(G%IsdB:G%IedB, G%jsd:G%jed), source=0.0) + if (.not.allocated(visc%kv_tbl_shelf_v)) & + allocate(visc%kv_tbl_shelf_v(G%isd:G%ied, G%JsdB:G%JedB), source=0.0) ! With a linear drag law under shelves, the friction velocity is already known. ! if (CS%linear_drag) ustar(:) = cdrag_sqrt_Z*CS%drag_bg_vel @@ -1456,8 +1494,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS) if (.not.CS%linear_drag) then v_at_u = set_v_at_u(v, h, G, GV, i, j, k, mask_v, OBC) - hutot = hutot + hweight * sqrt(u(I,j,k)**2 + & - v_at_u**2 + U_bg_sq) + hutot = hutot + hweight * sqrt(u(I,j,k)**2 + v_at_u**2 + U_bg_sq) endif if (use_EOS) then Thtot(I) = Thtot(I) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k)) @@ -1466,7 +1503,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS) enddo ; endif if ((.not.CS%linear_drag) .and. (hwtot > 0.0)) then - ustar(I) = cdrag_sqrt_Z * hutot/hwtot + ustar(I) = cdrag_sqrt_Z * hutot / hwtot else ustar(I) = cdrag_sqrt_Z * CS%drag_bg_vel endif @@ -1694,8 +1731,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS) if (.not.CS%linear_drag) then u_at_v = set_u_at_v(u, h, G, GV, i, J, k, mask_u, OBC) - hutot = hutot + hweight * sqrt(v(i,J,k)**2 + & - u_at_v**2 + U_bg_sq) + hutot = hutot + hweight * sqrt(v(i,J,k)**2 + u_at_v**2 + U_bg_sq) endif if (use_EOS) then Thtot(i) = Thtot(i) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k)) @@ -1704,7 +1740,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS) enddo ; endif if (.not.CS%linear_drag) then ; if (hwtot > 0.0) then - ustar(i) = cdrag_sqrt_Z * hutot/hwtot + ustar(i) = cdrag_sqrt_Z * hutot / hwtot else ustar(i) = cdrag_sqrt_Z * CS%drag_bg_vel endif ; endif @@ -1793,14 +1829,12 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS) enddo ! J-loop at v-points if (CS%debug) then - if (associated(visc%nkml_visc_u) .and. associated(visc%nkml_visc_v)) & + if (allocated(visc%nkml_visc_u) .and. allocated(visc%nkml_visc_v)) & call uvchksum("nkml_visc_[uv]", visc%nkml_visc_u, visc%nkml_visc_v, & G%HI, haloshift=0, scalar_pair=.true.) endif - if (CS%id_nkml_visc_u > 0) & - call post_data(CS%id_nkml_visc_u, visc%nkml_visc_u, CS%diag) - if (CS%id_nkml_visc_v > 0) & - call post_data(CS%id_nkml_visc_v, visc%nkml_visc_v, CS%diag) + if (CS%id_nkml_visc_u > 0) call post_data(CS%id_nkml_visc_u, visc%nkml_visc_u, CS%diag) + if (CS%id_nkml_visc_v > 0) call post_data(CS%id_nkml_visc_v, visc%nkml_visc_v, CS%diag) end subroutine set_viscous_ML @@ -1966,6 +2000,11 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS "may be an assumed value or it may be based on the "//& "actual velocity in the bottommost HBBL, depending on "//& "LINEAR_DRAG.", default=.true.) + call get_param(param_file, mdl, "DRAG_AS_BODY_FORCE", CS%body_force_drag, & + "If true, the bottom stress is imposed as an explicit body force "//& + "applied over a fixed distance from the bottom, rather than as an "//& + "implicit calculation based on an enhanced near-bottom viscosity", & + default=.false., do_not_log=.not.CS%bottomdraglaw) call get_param(param_file, mdl, "CHANNEL_DRAG", CS%Channel_drag, & "If true, the bottom drag is exerted directly on each "//& "layer proportional to the fraction of the bottom it "//& @@ -2145,8 +2184,8 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS if (CS%bottomdraglaw) then allocate(visc%bbl_thick_u(IsdB:IedB,jsd:jed), source=0.0) - allocate(visc%kv_bbl_u(IsdB:IedB,jsd:jed), source=0.0) allocate(visc%bbl_thick_v(isd:ied,JsdB:JedB), source=0.0) + allocate(visc%kv_bbl_u(IsdB:IedB,jsd:jed), source=0.0) allocate(visc%kv_bbl_v(isd:ied,JsdB:JedB), source=0.0) allocate(visc%ustar_bbl(isd:ied,jsd:jed), source=0.0) allocate(visc%TKE_bbl(isd:ied,jsd:jed), source=0.0) @@ -2177,7 +2216,7 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS call pass_var(CS%tideamp,G%domain) endif endif - if (CS%Channel_drag) then + if (CS%Channel_drag .or. CS%body_force_drag) then allocate(visc%Ray_u(IsdB:IedB,jsd:jed,nz), source=0.0) allocate(visc%Ray_v(isd:ied,JsdB:JedB,nz), source=0.0) CS%id_Ray_u = register_diag_field('ocean_model', 'Rayleigh_u', diag%axesCuL, & @@ -2243,31 +2282,30 @@ subroutine set_visc_end(visc, CS) !! related fields. Elements are deallocated here. type(set_visc_CS), intent(inout) :: CS !< The control structure returned by a previous !! call to set_visc_init. - if (CS%bottomdraglaw) then - deallocate(visc%bbl_thick_u) ; deallocate(visc%bbl_thick_v) - deallocate(visc%kv_bbl_u) ; deallocate(visc%kv_bbl_v) - if (allocated(CS%bbl_u)) deallocate(CS%bbl_u) - if (allocated(CS%bbl_v)) deallocate(CS%bbl_v) - endif - if (CS%Channel_drag) then - deallocate(visc%Ray_u) ; deallocate(visc%Ray_v) - endif - if (CS%dynamic_viscous_ML) then - deallocate(visc%nkml_visc_u) ; deallocate(visc%nkml_visc_v) - endif + + if (allocated(visc%bbl_thick_u)) deallocate(visc%bbl_thick_u) + if (allocated(visc%bbl_thick_v)) deallocate(visc%bbl_thick_v) + if (allocated(visc%kv_bbl_u)) deallocate(visc%kv_bbl_u) + if (allocated(visc%kv_bbl_v)) deallocate(visc%kv_bbl_v) + if (allocated(CS%bbl_u)) deallocate(CS%bbl_u) + if (allocated(CS%bbl_v)) deallocate(CS%bbl_v) + if (allocated(visc%Ray_u)) deallocate(visc%Ray_u) + if (allocated(visc%Ray_v)) deallocate(visc%Ray_v) + if (allocated(visc%nkml_visc_u)) deallocate(visc%nkml_visc_u) + if (allocated(visc%nkml_visc_v)) deallocate(visc%nkml_visc_v) if (associated(visc%Kd_shear)) deallocate(visc%Kd_shear) if (associated(visc%Kv_slow)) deallocate(visc%Kv_slow) if (associated(visc%TKE_turb)) deallocate(visc%TKE_turb) if (associated(visc%Kv_shear)) deallocate(visc%Kv_shear) if (associated(visc%Kv_shear_Bu)) deallocate(visc%Kv_shear_Bu) - if (associated(visc%ustar_bbl)) deallocate(visc%ustar_bbl) - if (associated(visc%TKE_bbl)) deallocate(visc%TKE_bbl) - if (associated(visc%taux_shelf)) deallocate(visc%taux_shelf) - if (associated(visc%tauy_shelf)) deallocate(visc%tauy_shelf) - if (associated(visc%tbl_thick_shelf_u)) deallocate(visc%tbl_thick_shelf_u) - if (associated(visc%tbl_thick_shelf_v)) deallocate(visc%tbl_thick_shelf_v) - if (associated(visc%kv_tbl_shelf_u)) deallocate(visc%kv_tbl_shelf_u) - if (associated(visc%kv_tbl_shelf_v)) deallocate(visc%kv_tbl_shelf_v) + if (allocated(visc%ustar_bbl)) deallocate(visc%ustar_bbl) + if (allocated(visc%TKE_bbl)) deallocate(visc%TKE_bbl) + if (allocated(visc%taux_shelf)) deallocate(visc%taux_shelf) + if (allocated(visc%tauy_shelf)) deallocate(visc%tauy_shelf) + if (allocated(visc%tbl_thick_shelf_u)) deallocate(visc%tbl_thick_shelf_u) + if (allocated(visc%tbl_thick_shelf_v)) deallocate(visc%tbl_thick_shelf_v) + if (allocated(visc%kv_tbl_shelf_u)) deallocate(visc%kv_tbl_shelf_u) + if (allocated(visc%kv_tbl_shelf_v)) deallocate(visc%kv_tbl_shelf_v) end subroutine set_visc_end !> \namespace mom_set_visc diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 4067f13757..855d563efc 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -86,9 +86,6 @@ module MOM_vert_friction !! may be an assumed value or it may be based on the !! actual velocity in the bottommost HBBL, depending !! on whether linear_drag is true. - logical :: Channel_drag !< If true, the drag is exerted directly on each - !! layer according to what fraction of the bottom - !! they overlie. logical :: harmonic_visc !< If true, the harmonic mean thicknesses are used !! to calculate the viscous coupling between layers !! except near the bottom. Otherwise the arithmetic @@ -280,7 +277,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & zDS = 0.0 stress = dt_Rho0 * forces%taux(I,j) do k=1,nz - h_a = 0.5 * (h(I,j,k) + h(I+1,j,k)) + h_neglect + h_a = 0.5 * (h(i,j,k) + h(i+1,j,k)) + h_neglect hfr = 1.0 ; if ((zDS+h_a) > Hmix) hfr = (Hmix - zDS) / h_a u(I,j,k) = u(I,j,k) + I_Hmix * hfr * stress if (associated(ADp%du_dt_str)) ADp%du_dt_str(i,J,k) = (I_Hmix * hfr * stress) * Idt @@ -291,7 +288,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & surface_stress(I) = dt_Rho0 * (G%mask2dCu(I,j)*forces%taux(I,j)) enddo ; endif ! direct_stress - if (CS%Channel_drag) then ; do k=1,nz ; do I=Isq,Ieq + if (allocated(visc%Ray_u)) then ; do k=1,nz ; do I=Isq,Ieq Ray(I,k) = visc%Ray_u(I,j,k) enddo ; enddo ; endif @@ -309,7 +306,10 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & ! calculate c'_k = - c_k / (b_k + a_k c'_(k-1)) ! and d'_k = (d_k - a_k d'_(k-1)) / (b_k + a_k c'_(k-1)) ! where c'_1 = c_1/b_1 and d'_1 = d_1/b_1 - ! (see Thomas' tridiagonal matrix algorithm) + ! + ! This form is mathematically equivalent to Thomas' tridiagonal matrix algorithm, but it + ! does not suffer from the acute sensitivity to truncation errors of the Thomas algorithm + ! because it involves no subtraction, as discussed by Schopf & Loughe, MWR, 1995. ! ! b1 is the denominator term 1 / (b_k + a_k c'_(k-1)) ! b_denom_1 is (b_k + a_k + c_k) - a_k(1 - c'_(k-1)) @@ -357,7 +357,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & if (abs(ADp%du_dt_visc(I,j,k)) < accel_underflow) ADp%du_dt_visc(I,j,k) = 0.0 enddo ; enddo ; endif - if (associated(visc%taux_shelf)) then ; do I=Isq,Ieq + if (allocated(visc%taux_shelf)) then ; do I=Isq,Ieq visc%taux_shelf(I,j) = -GV%Rho0*CS%a1_shelf_u(I,j)*u(I,j,1) ! - u_shelf? enddo ; endif @@ -365,7 +365,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & do I=Isq,Ieq taux_bot(I,j) = GV%Rho0 * (u(I,j,nz)*CS%a_u(I,j,nz+1)) enddo - if (CS%Channel_drag) then ; do k=1,nz ; do I=Isq,Ieq + if (allocated(visc%Ray_u)) then ; do k=1,nz ; do I=Isq,Ieq taux_bot(I,j) = taux_bot(I,j) + GV%Rho0 * (Ray(I,k)*u(I,j,k)) enddo ; enddo ; endif endif @@ -418,7 +418,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & surface_stress(i) = dt_Rho0 * (G%mask2dCv(i,J)*forces%tauy(i,J)) enddo ; endif ! direct_stress - if (CS%Channel_drag) then ; do k=1,nz ; do i=is,ie + if (allocated(visc%Ray_v)) then ; do k=1,nz ; do i=is,ie Ray(i,k) = visc%Ray_v(i,J,k) enddo ; enddo ; endif @@ -457,7 +457,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & if (abs(ADp%dv_dt_visc(i,J,k)) < accel_underflow) ADp%dv_dt_visc(i,J,k) = 0.0 enddo ; enddo ; endif - if (associated(visc%tauy_shelf)) then ; do i=is,ie + if (allocated(visc%tauy_shelf)) then ; do i=is,ie visc%tauy_shelf(i,J) = -GV%Rho0*CS%a1_shelf_v(i,J)*v(i,J,1) ! - v_shelf? enddo ; endif @@ -465,7 +465,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & do i=is,ie tauy_bot(i,J) = GV%Rho0 * (v(i,J,nz)*CS%a_v(i,J,nz+1)) enddo - if (CS%Channel_drag) then ; do k=1,nz ; do i=is,ie + if (allocated(visc%Ray_v)) then ; do k=1,nz ; do i=is,ie tauy_bot(i,J) = tauy_bot(i,J) + GV%Rho0 * (Ray(i,k)*v(i,J,k)) enddo ; enddo ; endif endif @@ -595,7 +595,7 @@ subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, US, CS) do j=G%jsc,G%jec do I=Isq,Ieq ; do_i(I) = (G%mask2dCu(I,j) > 0.0) ; enddo - if (CS%Channel_drag) then ; do k=1,nz ; do I=Isq,Ieq + if (allocated(visc%Ray_u)) then ; do k=1,nz ; do I=Isq,Ieq Ray(I,k) = visc%Ray_u(I,j,k) enddo ; enddo ; endif @@ -624,7 +624,7 @@ subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, US, CS) do J=Jsq,Jeq do i=is,ie ; do_i(i) = (G%mask2dCv(i,J) > 0.0) ; enddo - if (CS%Channel_drag) then ; do k=1,nz ; do i=is,ie + if (allocated(visc%Ray_v)) then ; do k=1,nz ; do i=is,ie Ray(i,k) = visc%Ray_v(i,J,k) enddo ; enddo ; endif @@ -753,11 +753,11 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) if (CS%debug .or. (CS%id_hML_u > 0)) allocate(hML_u(G%IsdB:G%IedB,G%jsd:G%jed), source=0.0) if (CS%debug .or. (CS%id_hML_v > 0)) allocate(hML_v(G%isd:G%ied,G%JsdB:G%JedB), source=0.0) - if ((associated(visc%taux_shelf) .or. associated(forces%frac_shelf_u)) .and. & + if ((allocated(visc%taux_shelf) .or. associated(forces%frac_shelf_u)) .and. & .not.associated(CS%a1_shelf_u)) then allocate(CS%a1_shelf_u(G%IsdB:G%IedB,G%jsd:G%jed), source=0.0) endif - if ((associated(visc%tauy_shelf) .or. associated(forces%frac_shelf_v)) .and. & + if ((allocated(visc%tauy_shelf) .or. associated(forces%frac_shelf_v)) .and. & .not.associated(CS%a1_shelf_v)) then allocate(CS%a1_shelf_v(G%isd:G%ied,G%JsdB:G%JedB), source=0.0) endif @@ -1420,8 +1420,8 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, US, CS real :: maxvel ! Velocities components greater than maxvel real :: truncvel ! are truncated to truncvel, both [L T-1 ~> m s-1]. - real :: CFL ! The local CFL number. - real :: H_report ! A thickness below which not to report truncations. + real :: CFL ! The local CFL number [nondim] + real :: H_report ! A thickness below which not to report truncations [H ~> m or kg m-2] real :: vel_report(SZIB_(G),SZJB_(G)) ! The velocity to report [L T-1 ~> m s-1] real :: u_old(SZIB_(G),SZJ_(G),SZK_(GV)) ! The previous u-velocity [L T-1 ~> m s-1] real :: v_old(SZI_(G),SZJB_(G),SZK_(GV)) ! The previous v-velocity [L T-1 ~> m s-1] @@ -1512,10 +1512,9 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, US, CS if (len_trim(CS%u_trunc_file) > 0) then do j=js,je ; do I=Isq,Ieq ; if (dowrite(I,j)) then -! Here the diagnostic reporting subroutines are called if -! unphysically large values were found. + ! Call a diagnostic reporting subroutines are called if unphysically large values are found. call write_u_accel(I, j, u_old, h, ADp, CDp, dt, G, GV, US, CS%PointAccel_CSp, & - vel_report(I,j), forces%taux(I,j), a=CS%a_u, hv=CS%h_u) + vel_report(I,j), forces%taux(I,j), a=CS%a_u, hv=CS%h_u) endif ; enddo ; enddo endif @@ -1597,10 +1596,9 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, US, CS if (len_trim(CS%v_trunc_file) > 0) then do J=Jsq,Jeq ; do i=is,ie ; if (dowrite(i,J)) then -! Here the diagnostic reporting subroutines are called if -! unphysically large values were found. + ! Call a diagnostic reporting subroutines are called if unphysically large values are found. call write_v_accel(i, J, v_old, h, ADp, CDp, dt, G, GV, US, CS%PointAccel_CSp, & - vel_report(i,J), forces%tauy(i,J), a=CS%a_v, hv=CS%h_v) + vel_report(i,J), forces%tauy(i,J), a=CS%a_v, hv=CS%h_v) endif ; enddo ; enddo endif @@ -1668,10 +1666,6 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & "may be an assumed value or it may be based on the "//& "actual velocity in the bottommost HBBL, depending on "//& "LINEAR_DRAG.", default=.true.) - call get_param(param_file, mdl, "CHANNEL_DRAG", CS%Channel_drag, & - "If true, the bottom drag is exerted directly on each "//& - "layer proportional to the fraction of the bottom it "//& - "overlies.", default=.false.) call get_param(param_file, mdl, "DIRECT_STRESS", CS%direct_stress, & "If true, the wind stress is distributed over the "//& "topmost HMIX_STRESS of fluid (like in HYCOM), and KVML "//&