diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 6daa1e9c49..ab9cf0c444 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -1447,6 +1447,44 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, enddo ; enddo endif + ! Compute the instantaneous M2 and K1 velocities + ! Note here the filtering is applied to the velocity fields at the previous time step. + if (CS%use_filter_m2) then + call Filt_accum(ubt, um2, CS%Time, US, CS%Filt_CS_um2) + call Filt_accum(vbt, vm2, CS%Time, US, CS%Filt_CS_vm2) + endif + if (CS%use_filter_k1) then + call Filt_accum(ubt, uk1, CS%Time, US, CS%Filt_CS_uk1) + call Filt_accum(vbt, vk1, CS%Time, US, CS%Filt_CS_vk1) + endif + + ! Apply frequency-dependent linear wave drag + if (CS%linear_freq_drag) then + Drag_u(:,:) = 0.0 ; Drag_v(:,:) = 0.0 + !$OMP do + do j=js,je ; do I=is-1,ie + if (CS%lin_drag_um2(I,j) > 0.0 .or. CS%lin_drag_uk1(I,j) > 0.0) then + Htot = 0.5 * (eta(i,j) + eta(i+1,j)) + if (GV%Boussinesq) & + Htot = Htot + 0.5*GV%Z_to_H * (CS%bathyT(i,j) + CS%bathyT(i+1,j)) + Drag_u(I,j) = (um2(I,j) * CS%lin_drag_um2(I,j) + & + uk1(I,j) * CS%lin_drag_uk1(I,j)) / Htot + BT_force_u(I,j) = BT_force_u(I,j) + Drag_u(I,j) + endif + enddo ; enddo + !$OMP do + do J=js-1,je ; do i=is,ie + if (CS%lin_drag_vm2(i,J) > 0.0 .or. CS%lin_drag_vk1(i,J) > 0.0) then + Htot = 0.5 * (eta(i,j) + eta(i,j+1)) + if (GV%Boussinesq) & + Htot = Htot + 0.5*GV%Z_to_H * (CS%bathyT(i,j) + CS%bathyT(i,j+1)) + Drag_v(i,J) = (vm2(i,J) * CS%lin_drag_vm2(i,J) + & + vk1(i,J) * CS%lin_drag_vk1(i,J)) / Htot + BT_force_v(i,J) = BT_force_v(i,J) + Drag_v(i,J) + endif + enddo ; enddo + endif + if ((Isq > is-1) .or. (Jsq > js-1)) then ! Non-symmetric memory is being used, so the edge values need to be ! filled in with a halo update of a non-symmetric array. @@ -1618,46 +1656,6 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, endif ; enddo ; enddo endif - ! Compute the instantaneous M2 and K1 velocities - if (CS%use_filter_m2) then - call Filt_accum(ubt, um2, CS%Time, US, CS%Filt_CS_um2) - call Filt_accum(vbt, vm2, CS%Time, US, CS%Filt_CS_vm2) - endif - if (CS%use_filter_k1) then - call Filt_accum(ubt, uk1, CS%Time, US, CS%Filt_CS_uk1) - call Filt_accum(vbt, vk1, CS%Time, US, CS%Filt_CS_vk1) - endif - - ! Apply frequency-dependent linear wave drag - Drag_u(:,:) = 0.0 ; Drag_v(:,:) = 0.0 - - if (CS%linear_freq_drag) then - !$OMP do - do j=js,je ; do I=is-1,ie - if (CS%lin_drag_um2(I,j) > 0.0 .or. CS%lin_drag_uk1(I,j) > 0.0) then - Htot = 0.5 * (eta(i,j) + eta(i+1,j)) - if (GV%Boussinesq) & - Htot = Htot + 0.5*GV%Z_to_H * (CS%bathyT(i,j) + CS%bathyT(i+1,j)) - if (CS%lin_drag_um2(I,j) > 0.0) & - Drag_u(I,j) = Drag_u(I,j) + um2(I,j) * CS%lin_drag_um2(I,j) / Htot - if (CS%lin_drag_uk1(I,j) > 0.0) & - Drag_u(I,j) = Drag_u(I,j) + uk1(I,j) * CS%lin_drag_uk1(I,j) / Htot - endif - enddo ; enddo - !$OMP do - do J=js-1,je ; do i=is,ie - if (CS%lin_drag_vm2(i,J) > 0.0 .or. CS%lin_drag_vk1(i,J) > 0.0) then - Htot = 0.5 * (eta(i,j) + eta(i,j+1)) - if (GV%Boussinesq) & - Htot = Htot + 0.5*GV%Z_to_H * (CS%bathyT(i,j) + CS%bathyT(i,j+1)) - if (CS%lin_drag_vm2(i,J) > 0.0) & - Drag_v(i,J) = Drag_v(i,J) + vm2(i,J) * CS%lin_drag_vm2(i,J) / Htot - if (CS%lin_drag_vk1(i,J) > 0.0) & - Drag_v(i,J) = Drag_v(i,J) + vk1(i,J) * CS%lin_drag_vk1(i,J) / Htot - endif - enddo ; enddo - endif - ! Zero out the arrays for various time-averaged quantities. if (find_etaav) then !$OMP do @@ -2103,21 +2101,33 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, 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) + & - dtbt * ((BT_force_v(i,J) + Cor_v(i,J)) + PFv(i,J) - Drag_v(i,J))) + dtbt * ((BT_force_v(i,J) + Cor_v(i,J)) + PFv(i,J))) if (abs(vbt(i,J)) < CS%vel_underflow) vbt(i,J) = 0.0 vbt_trans(i,J) = trans_wt1*vbt(i,J) + trans_wt2*vel_prev enddo ; enddo - if (CS%linear_wave_drag) then + if (CS%linear_wave_drag .and. CS%linear_freq_drag) then !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv-1,iev+1 v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * & - ((Cor_v(i,J) + PFv(i,J)) - vbt(i,J)*Rayleigh_v(i,J) - Drag_v(i,J)) + ((Cor_v(i,J) + PFv(i,J)) - (vbt(i,J)*Rayleigh_v(i,J) + Drag_v(i,J))) + enddo ; enddo + elseif (CS%linear_wave_drag) then + !$OMP do schedule(static) + do J=jsv-1,jev ; do i=isv-1,iev+1 + v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * & + ((Cor_v(i,J) + PFv(i,J)) - vbt(i,J)*Rayleigh_v(i,J)) + enddo ; enddo + elseif (CS%linear_freq_drag) then + !$OMP do schedule(static) + do J=jsv-1,jev ; do i=isv-1,iev+1 + v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * & + ((Cor_v(i,J) + PFv(i,J)) - Drag_v(i,J)) enddo ; enddo else !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv-1,iev+1 - v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * (Cor_v(i,J) + PFv(i,J) - Drag_v(i,J)) + v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * (Cor_v(i,J) + PFv(i,J)) enddo ; enddo endif @@ -2180,23 +2190,37 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, 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) + & - dtbt * ((BT_force_u(I,j) + Cor_u(I,j)) + PFu(I,j) - Drag_u(I,j))) + dtbt * ((BT_force_u(I,j) + Cor_u(I,j)) + PFu(I,j))) if (abs(ubt(I,j)) < CS%vel_underflow) ubt(I,j) = 0.0 ubt_trans(I,j) = trans_wt1*ubt(I,j) + trans_wt2*vel_prev enddo ; enddo !$OMP end do nowait - if (CS%linear_wave_drag) then + if (CS%linear_wave_drag .and. CS%linear_freq_drag) then !$OMP do schedule(static) do j=jsv,jev ; do I=isv-1,iev u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * & - ((Cor_u(I,j) + PFu(I,j)) - ubt(I,j)*Rayleigh_u(I,j) - Drag_u(I,j)) + ((Cor_u(I,j) + PFu(I,j)) - (ubt(I,j)*Rayleigh_u(I,j) + Drag_u(I,j))) + enddo ; enddo + !$OMP end do nowait + elseif (CS%linear_wave_drag) then + !$OMP do schedule(static) + do j=jsv,jev ; do I=isv-1,iev + u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * & + ((Cor_u(I,j) + PFu(I,j)) - ubt(I,j)*Rayleigh_u(I,j)) + enddo ; enddo + !$OMP end do nowait + elseif (CS%linear_freq_drag) then + !$OMP do schedule(static) + do j=jsv,jev ; do I=isv-1,iev + u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * & + ((Cor_u(I,j) + PFu(I,j)) - Drag_u(I,j)) enddo ; enddo !$OMP end do nowait else !$OMP do schedule(static) do j=jsv,jev ; do I=isv-1,iev - u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * (Cor_u(I,j) + PFu(I,j) - Drag_u(I,j)) + u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * (Cor_u(I,j) + PFu(I,j)) enddo ; enddo !$OMP end do nowait endif @@ -2258,21 +2282,33 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, 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) + & - dtbt * ((BT_force_u(I,j) + Cor_u(I,j)) + PFu(I,j) - Drag_u(I,j))) + dtbt * ((BT_force_u(I,j) + Cor_u(I,j)) + PFu(I,j))) if (abs(ubt(I,j)) < CS%vel_underflow) ubt(I,j) = 0.0 ubt_trans(I,j) = trans_wt1*ubt(I,j) + trans_wt2*vel_prev enddo ; enddo - if (CS%linear_wave_drag) then + if (CS%linear_wave_drag .and. CS%linear_freq_drag) then !$OMP do schedule(static) do j=jsv-1,jev+1 ; do I=isv-1,iev u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * & - ((Cor_u(I,j) + PFu(I,j)) - ubt(I,j)*Rayleigh_u(I,j) - Drag_u(I,j)) + ((Cor_u(I,j) + PFu(I,j)) - (ubt(I,j)*Rayleigh_u(I,j) + Drag_u(I,j))) + enddo ; enddo + elseif (CS%linear_wave_drag) then + !$OMP do schedule(static) + do j=jsv,jev ; do I=isv-1,iev + u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * & + ((Cor_u(I,j) + PFu(I,j)) - ubt(I,j)*Rayleigh_u(I,j)) + enddo ; enddo + elseif (CS%linear_freq_drag) then + !$OMP do schedule(static) + do j=jsv,jev ; do I=isv-1,iev + u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * & + ((Cor_u(I,j) + PFu(I,j)) - Drag_u(I,j)) enddo ; enddo else !$OMP do schedule(static) do j=jsv-1,jev+1 ; do I=isv-1,iev - u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * (Cor_u(I,j) + PFu(I,j) - Drag_u(I,j)) + u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * (Cor_u(I,j) + PFu(I,j)) enddo ; enddo endif @@ -2346,23 +2382,37 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, 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) + & - dtbt * ((BT_force_v(i,J) + Cor_v(i,J)) + PFv(i,J) - Drag_v(i,J))) + dtbt * ((BT_force_v(i,J) + Cor_v(i,J)) + PFv(i,J))) if (abs(vbt(i,J)) < CS%vel_underflow) vbt(i,J) = 0.0 vbt_trans(i,J) = trans_wt1*vbt(i,J) + trans_wt2*vel_prev enddo ; enddo !$OMP end do nowait - if (CS%linear_wave_drag) then + if (CS%linear_wave_drag .and. CS%linear_freq_drag) then !$OMP do schedule(static) - do J=jsv-1,jev ; do i=isv,iev + do J=jsv-1,jev ; do i=isv-1,iev+1 + v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * & + ((Cor_v(i,J) + PFv(i,J)) - (vbt(i,J)*Rayleigh_v(i,J) + Drag_v(i,J))) + enddo ; enddo + !$OMP end do nowait + elseif (CS%linear_wave_drag) then + !$OMP do schedule(static) + do J=jsv-1,jev ; do i=isv-1,iev+1 + v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * & + ((Cor_v(i,J) + PFv(i,J)) - vbt(i,J)*Rayleigh_v(i,J)) + enddo ; enddo + !$OMP end do nowait + elseif (CS%linear_freq_drag) then + !$OMP do schedule(static) + do J=jsv-1,jev ; do i=isv-1,iev+1 v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * & - ((Cor_v(i,J) + PFv(i,J)) - vbt(i,J)*Rayleigh_v(i,J) - Drag_v(i,J)) + ((Cor_v(i,J) + PFv(i,J)) - Drag_v(i,J)) enddo ; enddo !$OMP end do nowait else !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv,iev - v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * (Cor_v(i,J) + PFv(i,J) - Drag_v(i,J)) + v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * (Cor_v(i,J) + PFv(i,J)) enddo ; enddo !$OMP end do nowait endif @@ -5068,25 +5118,24 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, if (len_trim(wave_drag_u) > 0 .and. len_trim(wave_drag_v) > 0) then call MOM_read_data(wave_drag_file, wave_drag_u, CS%lin_drag_u, G%Domain, & - position=EAST_FACE, scale=GV%m_to_H*US%T_to_s) + position=EAST_FACE, scale=wave_drag_scale*GV%m_to_H*US%T_to_s) call pass_var(CS%lin_drag_u, G%Domain) - CS%lin_drag_u(:,:) = wave_drag_scale * CS%lin_drag_u(:,:) call MOM_read_data(wave_drag_file, wave_drag_v, CS%lin_drag_v, G%Domain, & - position=NORTH_FACE, scale=GV%m_to_H*US%T_to_s) + position=NORTH_FACE, scale=wave_drag_scale*GV%m_to_H*US%T_to_s) call pass_var(CS%lin_drag_v, G%Domain) - CS%lin_drag_v(:,:) = wave_drag_scale * CS%lin_drag_v(:,:) elseif (len_trim(wave_drag_var) > 0) then allocate(lin_drag_h(isd:ied,jsd:jed), source=0.0) - call MOM_read_data(wave_drag_file, wave_drag_var, lin_drag_h, G%Domain, scale=GV%m_to_H*US%T_to_s) + call MOM_read_data(wave_drag_file, wave_drag_var, lin_drag_h, G%Domain, & + scale=wave_drag_scale*GV%m_to_H*US%T_to_s) call pass_var(lin_drag_h, G%Domain) do j=js,je ; do I=is-1,ie - CS%lin_drag_u(I,j) = wave_drag_scale * 0.5 * (lin_drag_h(i,j) + lin_drag_h(i+1,j)) + CS%lin_drag_u(I,j) = 0.5 * (lin_drag_h(i,j) + lin_drag_h(i+1,j)) enddo ; enddo do J=js-1,je ; do i=is,ie - CS%lin_drag_v(i,J) = wave_drag_scale * 0.5 * (lin_drag_h(i,j) + lin_drag_h(i,j+1)) + CS%lin_drag_v(i,J) = 0.5 * (lin_drag_h(i,j) + lin_drag_h(i,j+1)) enddo ; enddo deallocate(lin_drag_h) endif ! (len_trim(wave_drag_u) > 0 .and. len_trim(wave_drag_v) > 0) @@ -5109,25 +5158,24 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, if (CS%use_filter_m2 .and. m2_drag_scale > 0.0) then if (len_trim(m2_drag_u) > 0 .and. len_trim(m2_drag_v) > 0) then call MOM_read_data(wave_drag_file, m2_drag_u, CS%lin_drag_um2, G%Domain, & - position=EAST_FACE, scale=GV%m_to_H*US%T_to_s) + position=EAST_FACE, scale=m2_drag_scale*GV%m_to_H*US%T_to_s) call pass_var(CS%lin_drag_um2, G%Domain) - CS%lin_drag_um2(:,:) = m2_drag_scale * CS%lin_drag_um2(:,:) call MOM_read_data(wave_drag_file, m2_drag_v, CS%lin_drag_vm2, G%Domain, & - position=NORTH_FACE, scale=GV%m_to_H*US%T_to_s) + position=NORTH_FACE, scale=m2_drag_scale*GV%m_to_H*US%T_to_s) call pass_var(CS%lin_drag_vm2, G%Domain) - CS%lin_drag_vm2(:,:) = m2_drag_scale * CS%lin_drag_vm2(:,:) elseif (len_trim(m2_drag_var) > 0) then allocate(lin_drag_h(isd:ied,jsd:jed), source=0.0) - call MOM_read_data(wave_drag_file, m2_drag_var, lin_drag_h, G%Domain, scale=GV%m_to_H*US%T_to_s) + call MOM_read_data(wave_drag_file, m2_drag_var, lin_drag_h, G%Domain, & + scale=m2_drag_scale*GV%m_to_H*US%T_to_s) call pass_var(lin_drag_h, G%Domain) do j=js,je ; do I=is-1,ie - CS%lin_drag_um2(I,j) = m2_drag_scale * 0.5 * (lin_drag_h(i,j) + lin_drag_h(i+1,j)) + CS%lin_drag_um2(I,j) = 0.5 * (lin_drag_h(i,j) + lin_drag_h(i+1,j)) enddo ; enddo do J=js-1,je ; do i=is,ie - CS%lin_drag_vm2(i,J) = m2_drag_scale * 0.5 * (lin_drag_h(i,j) + lin_drag_h(i,j+1)) + CS%lin_drag_vm2(i,J) = 0.5 * (lin_drag_h(i,j) + lin_drag_h(i,j+1)) enddo ; enddo deallocate(lin_drag_h) endif ! (len_trim(m2_drag_u) > 0 .and. len_trim(m2_drag_v) > 0) @@ -5136,25 +5184,24 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, if (CS%use_filter_k1 .and. k1_drag_scale > 0.0) then if (len_trim(k1_drag_u) > 0 .and. len_trim(k1_drag_v) > 0) then call MOM_read_data(wave_drag_file, k1_drag_u, CS%lin_drag_uk1, G%Domain, & - position=EAST_FACE, scale=GV%m_to_H*US%T_to_s) + position=EAST_FACE, scale=k1_drag_scale*GV%m_to_H*US%T_to_s) call pass_var(CS%lin_drag_uk1, G%Domain) - CS%lin_drag_uk1(:,:) = k1_drag_scale * CS%lin_drag_uk1(:,:) call MOM_read_data(wave_drag_file, k1_drag_v, CS%lin_drag_vk1, G%Domain, & - position=NORTH_FACE, scale=GV%m_to_H*US%T_to_s) + position=NORTH_FACE, scale=k1_drag_scale*GV%m_to_H*US%T_to_s) call pass_var(CS%lin_drag_vk1, G%Domain) - CS%lin_drag_vk1(:,:) = k1_drag_scale * CS%lin_drag_vk1(:,:) elseif (len_trim(k1_drag_var) > 0) then allocate(lin_drag_h(isd:ied,jsd:jed), source=0.0) - call MOM_read_data(wave_drag_file, k1_drag_var, lin_drag_h, G%Domain, scale=GV%m_to_H*US%T_to_s) + call MOM_read_data(wave_drag_file, k1_drag_var, lin_drag_h, G%Domain, & + scale=k1_drag_scale*GV%m_to_H*US%T_to_s) call pass_var(lin_drag_h, G%Domain) do j=js,je ; do I=is-1,ie - CS%lin_drag_uk1(I,j) = k1_drag_scale * 0.5 * (lin_drag_h(i,j) + lin_drag_h(i+1,j)) + CS%lin_drag_uk1(I,j) = 0.5 * (lin_drag_h(i,j) + lin_drag_h(i+1,j)) enddo ; enddo do J=js-1,je ; do i=is,ie - CS%lin_drag_vk1(i,J) = k1_drag_scale * 0.5 * (lin_drag_h(i,j) + lin_drag_h(i,j+1)) + CS%lin_drag_vk1(i,J) = 0.5 * (lin_drag_h(i,j) + lin_drag_h(i,j+1)) enddo ; enddo deallocate(lin_drag_h) endif ! (len_trim(k1_drag_u) > 0 .and. len_trim(k1_drag_v) > 0)