From d7dece9717c629ca92c7dcff2a2c2bb1723920cf Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 1 Jul 2021 20:37:15 -0400 Subject: [PATCH] +Add wind-stress acceleration diagnostics Added new diagnostics of the acceleration driven by the wind stress accelerations as redistributed by a timestep of the vertical viscosity and not lost to bottom drag within a timestep. This is also in the diagnostics of the accelerations due to the vertical viscosity, but the redistribution can be found from the difference of the two. Also added a diagnostic of the contribution of the wind stresses to kinetic energy, and applied an underflow limiter on both the new acceleration diagnostic and the existing viscous acceleration diagnostic. All solutions are bitwise identical, but there are new diagnostics. --- src/core/MOM_variables.F90 | 4 + src/diagnostics/MOM_diagnostics.F90 | 54 ++++++++--- .../vertical/MOM_vert_friction.F90 | 95 ++++++++++++++----- 3 files changed, 119 insertions(+), 34 deletions(-) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index f966ab2ad2..3f98a97052 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -169,6 +169,10 @@ module MOM_variables PFv => NULL(), & !< Meridional acceleration due to pressure forces [L T-2 ~> m s-2] du_dt_visc => NULL(), &!< Zonal acceleration due to vertical viscosity [L T-2 ~> m s-2] dv_dt_visc => NULL(), &!< Meridional acceleration due to vertical viscosity [L T-2 ~> m s-2] + du_dt_str => NULL(), & !< Zonal acceleration due to the surface stress (included + !! in du_dt_visc) [L T-2 ~> m s-2] + dv_dt_str => NULL(), & !< Meridional acceleration due to the surface stress (included + !! in dv_dt_visc) [L T-2 ~> m s-2] du_dt_dia => NULL(), & !< Zonal acceleration due to diapycnal mixing [L T-2 ~> m s-2] dv_dt_dia => NULL(), & !< Meridional acceleration due to diapycnal mixing [L T-2 ~> m s-2] u_accel_bt => NULL(), &!< Pointer to the zonal barotropic-solver acceleration [L T-2 ~> m s-2] diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index e6b01af33d..44b05cc081 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -107,6 +107,7 @@ module MOM_diagnostics !! of this spurious Coriolis source. KE_adv => NULL(), & !< KE source from along-layer advection [H L2 T-3 ~> m3 s-3] KE_visc => NULL(), & !< KE source from vertical viscosity [H L2 T-3 ~> m3 s-3] + KE_stress => NULL(), & !< KE source from surface stress (included in KE_visc) [H L2 T-3 ~> m3 s-3] KE_horvisc => NULL(), & !< KE source from horizontal viscosity [H L2 T-3 ~> m3 s-3] KE_dia => NULL() !< KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3] @@ -121,8 +122,8 @@ module MOM_diagnostics integer :: id_col_ht = -1, id_dh_dt = -1 integer :: id_KE = -1, id_dKEdt = -1 integer :: id_PE_to_KE = -1, id_KE_BT = -1 - integer :: id_KE_Coradv = -1 - integer :: id_KE_adv = -1, id_KE_visc = -1 + integer :: id_KE_Coradv = -1, id_KE_adv = -1 + integer :: id_KE_visc = -1, id_KE_stress = -1 integer :: id_KE_horvisc = -1, id_KE_dia = -1 integer :: id_uh_Rlay = -1, id_vh_Rlay = -1 integer :: id_uhGM_Rlay = -1, id_vhGM_Rlay = -1 @@ -1060,7 +1061,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS KE_h(i,j) = CS%KE(i,j,k) * CS%dh_dt(i,j,k) enddo ; enddo if (.not.G%symmetric) & - call do_group_pass(CS%pass_KE_uv, G%domain) + call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie CS%dKE_dt(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) & * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) @@ -1078,7 +1079,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%PFv(i,J,k) enddo ; enddo if (.not.G%symmetric) & - call do_group_pass(CS%pass_KE_uv, G%domain) + call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie CS%PE_to_KE(i,j,k) = 0.5 * G%IareaT(i,j) & * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) @@ -1096,7 +1097,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%v_accel_bt(i,J,k) enddo ; enddo if (.not.G%symmetric) & - call do_group_pass(CS%pass_KE_uv, G%domain) + call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie CS%KE_BT(i,j,k) = 0.5 * G%IareaT(i,j) & * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) @@ -1118,7 +1119,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS * (uh(I,j,k) - uh(I-1,j,k) + vh(i,J,k) - vh(i,J-1,k)) enddo ; enddo if (.not.G%symmetric) & - call do_group_pass(CS%pass_KE_uv, G%domain) + call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie CS%KE_CorAdv(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) & * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) @@ -1146,7 +1147,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS * (uh(I,j,k) - uh(I-1,j,k) + vh(i,J,k) - vh(i,J-1,k)) enddo ; enddo if (.not.G%symmetric) & - call do_group_pass(CS%pass_KE_uv, G%domain) + call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie CS%KE_adv(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) & * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) @@ -1164,7 +1165,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%dv_dt_visc(i,J,k) enddo ; enddo if (.not.G%symmetric) & - call do_group_pass(CS%pass_KE_uv, G%domain) + call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie CS%KE_visc(i,j,k) = 0.5 * G%IareaT(i,j) & * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) @@ -1173,6 +1174,24 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS if (CS%id_KE_visc > 0) call post_data(CS%id_KE_visc, CS%KE_visc, CS%diag) endif + if (associated(CS%KE_stress)) then + do k=1,nz + do j=js,je ; do I=Isq,Ieq + KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%du_dt_str(I,j,k) + enddo ; enddo + do J=Jsq,Jeq ; do i=is,ie + KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%dv_dt_str(i,J,k) + enddo ; enddo + if (.not.G%symmetric) & + call do_group_pass(CS%pass_KE_uv, G%domain) + do j=js,je ; do i=is,ie + CS%KE_stress(i,j,k) = 0.5 * G%IareaT(i,j) * & + ((KE_u(I,j) + KE_u(I-1,j)) + (KE_v(i,J) + KE_v(i,J-1))) + enddo ; enddo + enddo + if (CS%id_KE_stress > 0) call post_data(CS%id_KE_stress, CS%KE_stress, CS%diag) + endif + if (associated(CS%KE_horvisc)) then do k=1,nz do j=js,je ; do I=Isq,Ieq @@ -1182,7 +1201,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%diffv(i,J,k) enddo ; enddo if (.not.G%symmetric) & - call do_group_pass(CS%pass_KE_uv, G%domain) + call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie CS%KE_horvisc(i,j,k) = 0.5 * G%IareaT(i,j) & * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) @@ -1203,7 +1222,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS KE_h(i,j) = CS%KE(i,j,k) * (CDp%diapyc_vel(i,j,k) - CDp%diapyc_vel(i,j,k+1)) enddo ; enddo if (.not.G%symmetric) & - call do_group_pass(CS%pass_KE_uv, G%domain) + call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie CS%KE_dia(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) & * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) @@ -1894,6 +1913,11 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) if (CS%id_KE_visc>0) call safe_alloc_ptr(CS%KE_visc,isd,ied,jsd,jed,nz) + CS%id_KE_stress = register_diag_field('ocean_model', 'KE_stress', diag%axesTL, Time, & + 'Kinetic Energy Source from Surface Stresses or Body Wind Stress', & + 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) + if (CS%id_KE_stress>0) call safe_alloc_ptr(CS%KE_stress,isd,ied,jsd,jed,nz) + CS%id_KE_horvisc = register_diag_field('ocean_model', 'KE_horvisc', diag%axesTL, Time, & 'Kinetic Energy Source from Horizontal Viscosity', & 'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) @@ -2294,7 +2318,7 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS) if (associated(CS%dKE_dt) .or. associated(CS%PE_to_KE) .or. & associated(CS%KE_BT) .or. associated(CS%KE_CorAdv) .or. & - associated(CS%KE_adv) .or. associated(CS%KE_visc) .or. & + associated(CS%KE_adv) .or. associated(CS%KE_visc) .or. associated(CS%KE_stress) .or. & associated(CS%KE_horvisc) .or. associated(CS%KE_dia)) & call safe_alloc_ptr(CS%KE,isd,ied,jsd,jed,nz) @@ -2323,6 +2347,11 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS) call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz) endif + if (associated(CS%KE_stress)) then + call safe_alloc_ptr(ADp%du_dt_str,IsdB,IedB,jsd,jed,nz) + call safe_alloc_ptr(ADp%dv_dt_str,isd,ied,JsdB,JedB,nz) + endif + if (associated(CS%KE_dia)) then call safe_alloc_ptr(ADp%du_dt_dia,IsdB,IedB,jsd,jed,nz) call safe_alloc_ptr(ADp%dv_dt_dia,isd,ied,JsdB,JedB,nz) @@ -2353,6 +2382,7 @@ subroutine MOM_diagnostics_end(CS, ADp, CDp) if (associated(CS%KE_Coradv)) deallocate(CS%KE_Coradv) if (associated(CS%KE_adv)) deallocate(CS%KE_adv) if (associated(CS%KE_visc)) deallocate(CS%KE_visc) + if (associated(CS%KE_stress)) deallocate(CS%KE_stress) if (associated(CS%KE_horvisc)) deallocate(CS%KE_horvisc) if (associated(CS%KE_dia)) deallocate(CS%KE_dia) if (associated(CS%dv_dt)) deallocate(CS%dv_dt) @@ -2368,6 +2398,8 @@ subroutine MOM_diagnostics_end(CS, ADp, CDp) if (associated(ADp%gradKEu)) deallocate(ADp%gradKEu) if (associated(ADp%du_dt_visc)) deallocate(ADp%du_dt_visc) if (associated(ADp%dv_dt_visc)) deallocate(ADp%dv_dt_visc) + if (associated(ADp%du_dt_str)) deallocate(ADp%du_dt_str) + if (associated(ADp%dv_dt_str)) deallocate(ADp%dv_dt_str) if (associated(ADp%du_dt_dia)) deallocate(ADp%du_dt_dia) if (associated(ADp%dv_dt_dia)) deallocate(ADp%dv_dt_dia) if (associated(ADp%du_other)) deallocate(ADp%du_other) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 5b85c5f5f6..1d46f9aee3 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -121,6 +121,7 @@ module MOM_vert_friction !>@{ Diagnostic identifiers integer :: id_du_dt_visc = -1, id_dv_dt_visc = -1, id_au_vv = -1, id_av_vv = -1 + integer :: id_du_dt_str = -1, id_dv_dt_str = -1 integer :: id_h_u = -1, id_h_v = -1, id_hML_u = -1 , id_hML_v = -1 integer :: id_taux_bot = -1, id_tauy_bot = -1 integer :: id_Kv_slow = -1, id_Kv_u = -1, id_Kv_v = -1 @@ -207,6 +208,8 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & real :: stress ! The surface stress times the time step, divided ! by the density [H L T-1 ~> m2 s-1 or kg m-1 s-1]. + real :: accel_underflow ! An acceleration magnitude that is so small that values that are less + ! than this are diagnosed as 0 [L T-2 ~> m s-2]. real :: zDS, hfr, h_a ! Temporary variables used with direct_stress. real :: surface_stress(SZIB_(G))! The same as stress, unless the wind stress ! stress is applied as a body force [H L T-1 ~> m2 s-1 or kg m-1 s-1]. @@ -236,6 +239,8 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & h_neglect = GV%H_subroundoff Idt = 1.0 / dt + accel_underflow = CS%vel_underflow * Idt + !Check if Stokes mixing allowed if requested (present and associated) DoStokesMixing=.false. if (CS%StokesMixing) then @@ -265,9 +270,13 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & ADp%du_dt_visc(I,j,k) = u(I,j,k) enddo ; enddo ; endif -! One option is to have the wind stress applied as a body force -! over the topmost Hmix fluid. If DIRECT_STRESS is not defined, -! the wind stress is applied as a stress boundary condition. + if (associated(ADp%du_dt_str)) then ; do k=1,nz ; do I=Isq,Ieq + ADp%du_dt_str(I,j,k) = 0.0 + enddo ; enddo ; endif + + ! One option is to have the wind stress applied as a body force + ! over the topmost Hmix fluid. If DIRECT_STRESS is not defined, + ! the wind stress is applied as a stress boundary condition. if (CS%direct_stress) then do I=Isq,Ieq ; if (do_i(I)) then surface_stress(I) = 0.0 @@ -277,6 +286,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & 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 zDS = zDS + h_a ; if (zDS >= Hmix) exit enddo endif ; enddo ! end of i loop @@ -316,6 +326,8 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & b1(I) = 1.0 / (b_denom_1 + dt_Z_to_H*CS%a_u(I,j,2)) d1(I) = b_denom_1 * b1(I) u(I,j,1) = b1(I) * (CS%h_u(I,j,1) * u(I,j,1) + surface_stress(I)) + if (associated(ADp%du_dt_str)) & + ADp%du_dt_str(I,j,1) = b1(I) * (CS%h_u(I,j,1) * ADp%du_dt_str(I,j,1) + surface_stress(I)*Idt) endif ; enddo do k=2,nz ; do I=Isq,Ieq ; if (do_i(I)) then c1(I,k) = dt_Z_to_H * CS%a_u(I,j,K) * b1(I) @@ -324,6 +336,9 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & d1(I) = b_denom_1 * b1(I) u(I,j,k) = (CS%h_u(I,j,k) * u(I,j,k) + & dt_Z_to_H * CS%a_u(I,j,K) * u(I,j,k-1)) * b1(I) + if (associated(ADp%du_dt_str)) & + ADp%du_dt_str(I,j,k) = (CS%h_u(I,j,k) * ADp%du_dt_str(I,j,k) + & + dt_Z_to_H * CS%a_u(I,j,K) * ADp%du_dt_str(I,j,k-1)) * b1(I) endif ; enddo ; enddo ! back substitute to solve for the new velocities @@ -332,8 +347,17 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & u(I,j,k) = u(I,j,k) + c1(I,k+1) * u(I,j,k+1) endif ; enddo ; enddo ! i and k loops + if (associated(ADp%du_dt_str)) then + do i=is,ie ; if (abs(ADp%du_dt_str(I,j,nz)) < accel_underflow) ADp%du_dt_str(I,j,nz) = 0.0 ; enddo + do k=nz-1,1,-1 ; do I=Isq,Ieq ; if (do_i(I)) then + ADp%du_dt_str(I,j,k) = ADp%du_dt_str(I,j,k) + c1(I,k+1) * ADp%du_dt_str(I,j,k+1) + if (abs(ADp%du_dt_str(I,j,k)) < accel_underflow) ADp%du_dt_str(I,j,k) = 0.0 + endif ; enddo ; enddo + endif + if (associated(ADp%du_dt_visc)) then ; do k=1,nz ; do I=Isq,Ieq ADp%du_dt_visc(I,j,k) = (u(I,j,k) - ADp%du_dt_visc(I,j,k))*Idt + 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 @@ -373,9 +397,13 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & ADp%dv_dt_visc(i,J,k) = v(i,J,k) enddo ; enddo ; endif -! One option is to have the wind stress applied as a body force -! over the topmost Hmix fluid. If DIRECT_STRESS is not defined, -! the wind stress is applied as a stress boundary condition. + if (associated(ADp%dv_dt_str)) then ; do k=1,nz ; do i=is,ie + ADp%dv_dt_str(i,J,k) = 0.0 + enddo ; enddo ; endif + + ! One option is to have the wind stress applied as a body force + ! over the topmost Hmix fluid. If DIRECT_STRESS is not defined, + ! the wind stress is applied as a stress boundary condition. if (CS%direct_stress) then do i=is,ie ; if (do_i(i)) then surface_stress(i) = 0.0 @@ -385,6 +413,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & h_a = 0.5 * (h(i,J,k) + h(i,J+1,k)) + h_neglect hfr = 1.0 ; if ((zDS+h_a) > Hmix) hfr = (Hmix - zDS) / h_a v(i,J,k) = v(i,J,k) + I_Hmix * hfr * stress + if (associated(ADp%dv_dt_str)) ADp%dv_dt_str(i,J,k) = (I_Hmix * hfr * stress) * Idt zDS = zDS + h_a ; if (zDS >= Hmix) exit enddo endif ; enddo ! end of i loop @@ -401,6 +430,8 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & b1(i) = 1.0 / (b_denom_1 + dt_Z_to_H*CS%a_v(i,J,2)) d1(i) = b_denom_1 * b1(i) v(i,J,1) = b1(i) * (CS%h_v(i,J,1) * v(i,J,1) + surface_stress(i)) + if (associated(ADp%dv_dt_str)) & + ADp%dv_dt_str(i,J,1) = b1(i) * (CS%h_v(i,J,1) * ADp%dv_dt_str(i,J,1) + surface_stress(i)*Idt) endif ; enddo do k=2,nz ; do i=is,ie ; if (do_i(i)) then c1(i,k) = dt_Z_to_H * CS%a_v(i,J,K) * b1(i) @@ -408,13 +439,25 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & b1(i) = 1.0 / (b_denom_1 + dt_Z_to_H * CS%a_v(i,J,K+1)) d1(i) = b_denom_1 * b1(i) v(i,J,k) = (CS%h_v(i,J,k) * v(i,J,k) + dt_Z_to_H * CS%a_v(i,J,K) * v(i,J,k-1)) * b1(i) + if (associated(ADp%dv_dt_str)) & + ADp%dv_dt_str(i,J,k) = (CS%h_v(i,J,k) * ADp%dv_dt_str(i,J,k) + & + dt_Z_to_H * CS%a_v(i,J,K) * ADp%dv_dt_str(i,J,k-1)) * b1(i) endif ; enddo ; enddo do k=nz-1,1,-1 ; do i=is,ie ; if (do_i(i)) then v(i,J,k) = v(i,J,k) + c1(i,k+1) * v(i,J,k+1) endif ; enddo ; enddo ! i and k loops + if (associated(ADp%dv_dt_str)) then + do i=is,ie ; if (abs(ADp%dv_dt_str(i,J,nz)) < accel_underflow) ADp%dv_dt_str(i,J,nz) = 0.0 ; enddo + do k=nz-1,1,-1 ; do i=is,ie ; if (do_i(i)) then + ADp%dv_dt_str(i,J,k) = ADp%dv_dt_str(i,J,k) + c1(i,k+1) * ADp%dv_dt_str(i,J,k+1) + if (abs(ADp%dv_dt_str(i,J,k)) < accel_underflow) ADp%dv_dt_str(i,J,k) = 0.0 + endif ; enddo ; enddo + endif + if (associated(ADp%dv_dt_visc)) then ; do k=1,nz ; do i=is,ie ADp%dv_dt_visc(i,J,k) = (v(i,J,k) - ADp%dv_dt_visc(i,J,k))*Idt + 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 @@ -458,7 +501,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & enddo endif -! Offer diagnostic fields for averaging. + ! Offer diagnostic fields for averaging. if (CS%id_du_dt_visc > 0) & call post_data(CS%id_du_dt_visc, ADp%du_dt_visc, CS%diag) if (CS%id_dv_dt_visc > 0) & @@ -467,6 +510,10 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & call post_data(CS%id_taux_bot, taux_bot, CS%diag) if (present(tauy_bot) .and. (CS%id_tauy_bot > 0)) & call post_data(CS%id_tauy_bot, tauy_bot, CS%diag) + if (CS%id_du_dt_str > 0) & + call post_data(CS%id_du_dt_str, ADp%du_dt_str, CS%diag) + if (CS%id_dv_dt_str > 0) & + call post_data(CS%id_dv_dt_str, ADp%dv_dt_str, CS%diag) ! Diagnostics for terms multiplied by fractional thicknesses @@ -524,10 +571,9 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & end subroutine vertvisc -!> Calculate the fraction of momentum originally in a layer that remains -!! after a time-step of viscosity, and the fraction of a time-step's -!! worth of barotropic acceleration that a layer experiences after -!! viscosity is applied. +!> Calculate the fraction of momentum originally in a layer that remains in the water column +!! after a time-step of viscosity, equivalently the fraction of a time-step's worth of +!! barotropic acceleration that a layer experiences after viscosity is applied. subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, US, CS) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure @@ -566,10 +612,8 @@ subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, US, CS) do k=1,nz ; do i=Isq,Ieq ; Ray(i,k) = 0.0 ; enddo ; enddo - ! Find the zonal viscous using a modification of a standard tridagonal solver. -!$OMP parallel do default(none) shared(G,Isq,Ieq,CS,nz,visc,dt_Z_to_H,visc_rem_u) & -!$OMP firstprivate(Ray) & -!$OMP private(do_i,b_denom_1,b1,d1,c1) + ! Find the zonal viscous remnant using a modification of a standard tridagonal solver. + !$OMP parallel do default(shared) firstprivate(Ray) private(do_i,b_denom_1,b1,d1,c1) do j=G%jsc,G%jec do I=Isq,Ieq ; do_i(I) = (G%mask2dCu(I,j) > 0) ; enddo @@ -597,10 +641,8 @@ subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, US, CS) enddo ! end u-component j loop - ! Now find the meridional viscous using a modification. -!$OMP parallel do default(none) shared(Jsq,Jeq,is,ie,G,CS,visc,dt_Z_to_H,visc_rem_v,nz) & -!$OMP firstprivate(Ray) & -!$OMP private(do_i,b_denom_1,b1,d1,c1) + ! Now find the meridional viscous remnant using the robust tridiagonal solver. + !$OMP parallel do default(shared) firstprivate(Ray) private(do_i,b_denom_1,b1,d1,c1) do J=Jsq,Jeq do i=is,ie ; do_i(i) = (G%mask2dCv(i,J) > 0) ; enddo @@ -1813,13 +1855,20 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & 'Mixed Layer Thickness at Meridional Velocity Points for Viscosity', & thickness_units, conversion=GV%H_to_MKS) - CS%id_du_dt_visc = register_diag_field('ocean_model', 'du_dt_visc', diag%axesCuL, & - Time, 'Zonal Acceleration from Vertical Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2) + CS%id_du_dt_visc = register_diag_field('ocean_model', 'du_dt_visc', diag%axesCuL, Time, & + 'Zonal Acceleration from Vertical Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2) if (CS%id_du_dt_visc > 0) call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz) - CS%id_dv_dt_visc = register_diag_field('ocean_model', 'dv_dt_visc', diag%axesCvL, & - Time, 'Meridional Acceleration from Vertical Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2) + CS%id_dv_dt_visc = register_diag_field('ocean_model', 'dv_dt_visc', diag%axesCvL, Time, & + 'Meridional Acceleration from Vertical Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2) if (CS%id_dv_dt_visc > 0) call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz) + CS%id_du_dt_str = register_diag_field('ocean_model', 'du_dt_str', diag%axesCuL, Time, & + 'Zonal Acceleration from Surface Wind Stresses', 'm s-2', conversion=US%L_T2_to_m_s2) + if (CS%id_du_dt_str > 0) call safe_alloc_ptr(ADp%du_dt_str,IsdB,IedB,jsd,jed,nz) + CS%id_dv_dt_str = register_diag_field('ocean_model', 'dv_dt_str', diag%axesCvL, Time, & + 'Meridional Acceleration from Surface Wind Stresses', 'm s-2', conversion=US%L_T2_to_m_s2) + if (CS%id_dv_dt_str > 0) call safe_alloc_ptr(ADp%dv_dt_str,isd,ied,JsdB,JedB,nz) + CS%id_taux_bot = register_diag_field('ocean_model', 'taux_bot', diag%axesCu1, & Time, 'Zonal Bottom Stress from Ocean to Earth', & 'Pa', conversion=US%RZ_to_kg_m2*US%L_T2_to_m_s2)