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)