Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Visc_rem_[uv] multiplied momentum budget diagnostics #10

Merged
merged 13 commits into from
Aug 5, 2021
Merged

Visc_rem_[uv] multiplied momentum budget diagnostics #10

Show file tree
Hide file tree
Changes from 9 commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
7a65026
Implement visc remnant versions of PFu, CAu, u_accel_bt, diffu
NoraLoose Jul 24, 2021
2601b47
Add d[uv]_dt_visc_rem diagnostics
NoraLoose Jul 24, 2021
20b8bfc
Attempt to implement KE budget modified by visc rem
NoraLoose Jul 25, 2021
4208525
Fix syntax error
NoraLoose Jul 31, 2021
f4d4839
Fix some style errors
NoraLoose Jul 31, 2021
b1357bc
Fix more style errors
NoraLoose Jul 31, 2021
be09de0
Remove added KE budget diagnostics
NoraLoose Aug 4, 2021
fbdac20
Fix small bug
NoraLoose Aug 4, 2021
86741cd
Small style fix
NoraLoose Aug 4, 2021
a40a90f
Fix allocation problems
NoraLoose Aug 4, 2021
f31e5e6
Remove du_dt_visc_rem diagnostic
NoraLoose Aug 4, 2021
ca0fb89
Fix small inconsistency: CS%visc_rem_v --> CS%ADp%visc_rem_v
NoraLoose Aug 4, 2021
a874d6d
Don't need inout for ADp anymore
NoraLoose Aug 4, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 10 additions & 0 deletions src/core/MOM_barotropic.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2685,6 +2685,16 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
ADp%diag_hv(i,J,k) = BT_cont%h_v(i,J,k)
enddo ; enddo ; enddo
endif
if (present(ADp) .and. (associated(ADp%visc_rem_u))) then
do k=1,nz ; do j=js,je ; do I=is-1,ie
ADp%visc_rem_u(I,j,k) = visc_rem_u(I,j,k)
enddo ; enddo ; enddo
endif
if (present(ADp) .and. (associated(ADp%visc_rem_u))) then
do k=1,nz ; do J=js-1,je ; do i=is,ie
ADp%visc_rem_v(i,J,k) = visc_rem_v(i,J,k)
enddo ; enddo ; enddo
endif

if (G%nonblocking_updates) then
if (find_etaav) call complete_group_pass(CS%pass_etaav, G%Domain)
Expand Down
85 changes: 85 additions & 0 deletions src/core/MOM_dynamics_split_RK2.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -163,10 +163,12 @@ module MOM_dynamics_split_RK2
integer :: id_h_PFu = -1, id_h_PFv = -1
integer :: id_hf_PFu_2d = -1, id_hf_PFv_2d = -1
integer :: id_intz_PFu_2d = -1, id_intz_PFv_2d = -1
integer :: id_PFu_visc_rem = -1, id_PFv_visc_rem = -1
! integer :: id_hf_CAu = -1, id_hf_CAv = -1
integer :: id_h_CAu = -1, id_h_CAv = -1
integer :: id_hf_CAu_2d = -1, id_hf_CAv_2d = -1
integer :: id_intz_CAu_2d = -1, id_intz_CAv_2d = -1
integer :: id_CAu_visc_rem = -1, id_CAv_visc_rem = -1

! Split scheme only.
integer :: id_uav = -1, id_vav = -1
Expand All @@ -175,6 +177,7 @@ module MOM_dynamics_split_RK2
integer :: id_h_u_BT_accel = -1, id_h_v_BT_accel = -1
integer :: id_hf_u_BT_accel_2d = -1, id_hf_v_BT_accel_2d = -1
integer :: id_intz_u_BT_accel_2d = -1, id_intz_v_BT_accel_2d = -1
integer :: id_u_BT_accel_visc_rem = -1, id_v_BT_accel_visc_rem = -1
!>@}

type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
Expand Down Expand Up @@ -354,6 +357,12 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
real, dimension(SZI_(G),SZJB_(G)) :: &
intz_PFv_2d, intz_CAv_2d, intz_v_BT_accel_2d ! [H L T-2 ~> m2 s-2].

! Diagnostics for momentum budget terms multiplied by visc_rem_[uv],
real, allocatable, dimension(:,:,:) :: &
PFu_visc_rem, PFv_visc_rem, & ! Pressure force accel. x visc_rem_[uv] [L T-2 ~> m s-2].
CAu_visc_rem, CAv_visc_rem, & ! Coriolis force accel. x visc_rem_[uv] [L T-2 ~> m s-2].
u_BT_accel_visc_rem, v_BT_accel_visc_rem ! barotropic correction accel. x visc_rem_[uv] [L T-2 ~> m s-2].

real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s].

logical :: dyn_p_surf
Expand Down Expand Up @@ -1102,6 +1111,61 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
deallocate(h_v_BT_accel)
endif

if (CS%id_PFu_visc_rem > 0) then
allocate(PFu_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
PFu_visc_rem(:,:,:) = 0.0
do k=1,nz ; do j=js,je ; do I=Isq,Ieq
PFu_visc_rem(I,j,k) = CS%PFu(I,j,k) * CS%visc_rem_u(I,j,k)
enddo ; enddo ; enddo
call post_data(CS%id_PFu_visc_rem, PFu_visc_rem, CS%diag)
deallocate(PFu_visc_rem)
endif
if (CS%id_PFv_visc_rem > 0) then
allocate(PFv_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
PFv_visc_rem(:,:,:) = 0.0
do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
PFv_visc_rem(i,J,k) = CS%PFv(i,J,k) * CS%visc_rem_v(i,J,k)
enddo ; enddo ; enddo
call post_data(CS%id_PFv_visc_rem, PFv_visc_rem, CS%diag)
deallocate(PFv_visc_rem)
endif
if (CS%id_CAu_visc_rem > 0) then
allocate(CAu_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
CAu_visc_rem(:,:,:) = 0.0
do k=1,nz ; do j=js,je ; do I=Isq,Ieq
CAu_visc_rem(I,j,k) = CS%CAu(I,j,k) * CS%visc_rem_u(I,j,k)
enddo ; enddo ; enddo
call post_data(CS%id_CAu_visc_rem, CAu_visc_rem, CS%diag)
deallocate(CAu_visc_rem)
endif
if (CS%id_CAv_visc_rem > 0) then
allocate(CAv_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
CAv_visc_rem(:,:,:) = 0.0
do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
CAv_visc_rem(i,J,k) = CS%CAv(i,J,k) * CS%visc_rem_v(i,J,k)
enddo ; enddo ; enddo
call post_data(CS%id_CAv_visc_rem, CAv_visc_rem, CS%diag)
deallocate(CAv_visc_rem)
endif
if (CS%id_u_BT_accel_visc_rem > 0) then
allocate(u_BT_accel_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
u_BT_accel_visc_rem(:,:,:) = 0.0
do k=1,nz ; do j=js,je ; do I=Isq,Ieq
u_BT_accel_visc_rem(I,j,k) = CS%u_accel_bt(I,j,k) * CS%visc_rem_u(I,j,k)
enddo ; enddo ; enddo
call post_data(CS%id_u_BT_accel_visc_rem, u_BT_accel_visc_rem, CS%diag)
deallocate(u_BT_accel_visc_rem)
endif
if (CS%id_v_BT_accel_visc_rem > 0) then
allocate(v_BT_accel_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
v_BT_accel_visc_rem(:,:,:) = 0.0
do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
v_BT_accel_visc_rem(i,J,k) = CS%v_accel_bt(i,J,k) * CS%visc_rem_v(i,J,k)
enddo ; enddo ; enddo
call post_data(CS%id_v_BT_accel_visc_rem, v_BT_accel_visc_rem, CS%diag)
deallocate(v_BT_accel_visc_rem)
endif

if (CS%debug) then
call MOM_state_chksum("Corrector ", u, v, h, uh, vh, G, GV, US, symmetric=sym)
call uvchksum("Corrector avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s)
Expand Down Expand Up @@ -1606,6 +1670,27 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
if (CS%id_intz_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)

CS%id_PFu_visc_rem = register_diag_field('ocean_model', 'PFu_visc_rem', diag%axesCuL, Time, &
'Zonal Pressure Force Acceleration multiplied by the viscous remnant', 'm2 s-2', &
conversion=GV%H_to_m*US%L_T2_to_m_s2)
CS%id_PFv_visc_rem = register_diag_field('ocean_model', 'PFv_visc_rem', diag%axesCvL, Time, &
'Meridional Pressure Force Acceleration multiplied by the viscous remnant', 'm2 s-2', &
conversion=GV%H_to_m*US%L_T2_to_m_s2)

CS%id_CAu_visc_rem = register_diag_field('ocean_model', 'CAu_visc_rem', diag%axesCuL, Time, &
'Zonal Coriolis and Advective Acceleration multiplied by the viscous remnant', 'm2 s-2', &
conversion=GV%H_to_m*US%L_T2_to_m_s2)
CS%id_CAv_visc_rem = register_diag_field('ocean_model', 'CAv_visc_rem', diag%axesCvL, Time, &
'Meridional Coriolis and Advective Acceleration multiplied by the viscous remnant', 'm2 s-2', &
conversion=GV%H_to_m*US%L_T2_to_m_s2)

CS%id_u_BT_accel_visc_rem = register_diag_field('ocean_model', 'u_BT_accel_visc_rem', diag%axesCuL, Time, &
'Barotropic Anomaly Zonal Acceleration multiplied by the viscous remnant', 'm2 s-2', &
conversion=GV%H_to_m*US%L_T2_to_m_s2)
CS%id_v_BT_accel_visc_rem = register_diag_field('ocean_model', 'v_BT_accel_visc_rem', diag%axesCvL, Time, &
'Barotropic Anomaly Meridional Acceleration multiplied by the viscous remnant', 'm2 s-2', &
conversion=GV%H_to_m*US%L_T2_to_m_s2)

id_clock_Cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE)
id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=CLOCK_MODULE)
id_clock_pres = cpu_clock_id('(Ocean pressure force)', grain=CLOCK_MODULE)
Expand Down
7 changes: 6 additions & 1 deletion src/core/MOM_variables.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -161,6 +161,8 @@ module MOM_variables

! Each of the following fields has nz layers.
real, pointer, dimension(:,:,:) :: &
du_dt => NULL(), & !< Zonal acceleration [L T-2 ~> m s-2]
dv_dt => NULL(), & !< Meridional acceleration [L T-2 ~> m s-2]
diffu => NULL(), & !< Zonal acceleration due to along isopycnal viscosity [L T-2 ~> m s-2]
diffv => NULL(), & !< Meridional acceleration due to along isopycnal viscosity [L T-2 ~> m s-2]
CAu => NULL(), & !< Zonal Coriolis and momentum advection accelerations [L T-2 ~> m s-2]
Expand All @@ -172,7 +174,7 @@ module MOM_variables
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]
v_accel_bt => NULL() !< Pointer to the meridional barotropic-solver acceleration [L T-2 ~> m s-2]
v_accel_bt => NULL() !< Pointer to the meridional barotropic-solver acceleration [L T-2 ~> m s-2]
real, pointer, dimension(:,:,:) :: du_other => NULL()
!< Zonal velocity changes due to any other processes that are
!! not due to any explicit accelerations [L T-1 ~> m s-1].
Expand All @@ -191,6 +193,9 @@ module MOM_variables
real, pointer :: diag_hu(:,:,:) => NULL() !< layer thickness at u points
real, pointer :: diag_hv(:,:,:) => NULL() !< layer thickness at v points

real, pointer :: visc_rem_u(:,:,:) => NULL() !< viscous remnant at u points
real, pointer :: visc_rem_v(:,:,:) => NULL() !< viscous remnant at v points

end type accel_diag_ptrs

!> Pointers to arrays with transports, which can later be used for derived diagnostics, like energy balances.
Expand Down
9 changes: 7 additions & 2 deletions src/diagnostics/MOM_diagnostics.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -290,6 +290,12 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
if (CS%id_dv_dt>0) call post_data(CS%id_dv_dt, CS%dv_dt, CS%diag, alt_h = diag_pre_sync%h_state)

if (CS%id_dh_dt>0) call post_data(CS%id_dh_dt, CS%dh_dt, CS%diag, alt_h = diag_pre_sync%h_state)
do k=1,nz ; do j=js,je ; do I=is-1,ie
ADp%du_dt(I,j,k) = CS%du_dt(I,j,k)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps you forgot to delete this loop?
If not, since you are modifying the ADp structure, its intent should be inout:

  type(accel_diag_ptrs),   intent(inout)    :: ADp  !< structure with pointers to
                                                 !! accelerations in momentum equation.

Also, this loop should only be exercised if CS%id_hf_du_dt_2d > 0. Otherwise, there will be a memory error.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good point about the inout! I need ADp%du_dt for the computation of du_dt_visc_rem, see here:

https://github.com/NoraLoose/MOM6/blob/86741cd1d3b049303db325d61089263eda90979b/src/parameterizations/vertical/MOM_vert_friction.F90#L535-L536

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Gotcha! Then, you should also allocate ADp%du_dt when CS%id_du_dt_visc_rem > 0.

Copy link
Member Author

@NoraLoose NoraLoose Aug 4, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you mean that I should add something like:

if (CS%id_du_dt_visc_rem > 0) then
    if (.not.associated(CS%du_dt)) then
        call safe_alloc_ptr(CS%du_dt,IsdB,IedB,jsd,jed,nz)
        call register_time_deriv(lbound(MIS%u), MIS%u, CS%du_dt, CS)
    endif
endif

in MOM_diagnostics_init?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That won't work. In MOM_diagnostics_init, CS is the diagnostics_CS type.

enddo ; enddo ; enddo
do k=1,nz ; do J=js-1,je ; do i=is,ie
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same here, it seems this should be deleted since id_hf_du_dt is commented out in MOM_diagnostics_init. This look should only be exercised if ADp%dv_dt has been allocated.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same here. I need ADp%dv_dt to compute dv_dt_visc_rem in the MOM_vert_friction module. What do you recommend to do about allocation of ADp%du_dt and ADp%dv_dt? These variables are defined to be part of the ADp structure here:

https://github.com/NoraLoose/MOM6/blob/86741cd1d3b049303db325d61089263eda90979b/src/core/MOM_variables.F90#L164-L165

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suggest allocating them in MOM_diagnostics_init.

du_dt from diagnostics_CS is only allocated when id_du_dt > 0, see here. It also needs to be allocated when id_hf_du_dt > 0. Same for dv_dt.

ADp%dv_dt(i,J,k) = CS%dv_dt(i,J,k)
enddo ; enddo ; enddo

!! Diagnostics for terms multiplied by fractional thicknesses

Expand Down Expand Up @@ -2322,7 +2328,6 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS)
call safe_alloc_ptr(ADp%gradKEu,IsdB,IedB,jsd,jed,nz)
call safe_alloc_ptr(ADp%gradKEv,isd,ied,JsdB,JedB,nz)
endif

if (associated(CS%KE_visc)) then
call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)
call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
Expand Down Expand Up @@ -2368,7 +2373,7 @@ subroutine MOM_diagnostics_end(CS, ADp)
if (associated(CS%vhGM_Rlay)) deallocate(CS%vhGM_Rlay)

if (associated(ADp%gradKEu)) deallocate(ADp%gradKEu)
if (associated(ADp%gradKEu)) deallocate(ADp%gradKEu)
if (associated(ADp%gradKEv)) deallocate(ADp%gradKEv)
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_dia)) deallocate(ADp%du_dt_dia)
Expand Down
36 changes: 36 additions & 0 deletions src/parameterizations/lateral/MOM_hor_visc.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -190,6 +190,7 @@ module MOM_hor_visc
integer :: id_h_diffu = -1, id_h_diffv = -1
integer :: id_hf_diffu_2d = -1, id_hf_diffv_2d = -1
integer :: id_intz_diffu_2d = -1, id_intz_diffv_2d = -1
integer :: id_diffu_visc_rem = -1, id_diffv_visc_rem = -1
integer :: id_Ah_h = -1, id_Ah_q = -1
integer :: id_Kh_h = -1, id_Kh_q = -1
integer :: id_GME_coeff_h = -1, id_GME_coeff_q = -1
Expand Down Expand Up @@ -280,6 +281,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,

real, allocatable, dimension(:,:,:) :: h_diffu ! h x diffu [H L T-2 ~> m2 s-2]
real, allocatable, dimension(:,:,:) :: h_diffv ! h x diffv [H L T-2 ~> m2 s-2]
real, allocatable, dimension(:,:,:) :: diffu_visc_rem ! diffu x visc_rem_u [L T-2 ~> m2 s-2]
real, allocatable, dimension(:,:,:) :: diffv_visc_rem ! diffv x visc_rem_v [L T-2 ~> m2 s-2]

real, dimension(SZIB_(G),SZJB_(G)) :: &
dvdx, dudy, & ! components in the shearing strain [T-1 ~> s-1]
Expand Down Expand Up @@ -1704,6 +1707,25 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
deallocate(h_diffv)
endif

if (present(ADp) .and. (CS%id_diffu_visc_rem > 0)) then
allocate(diffu_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
diffu_visc_rem(:,:,:) = 0.0
do k=1,nz ; do j=js,je ; do I=Isq,Ieq
diffu_visc_rem(I,j,k) = diffu(I,j,k) * ADp%visc_rem_u(I,j,k)
enddo ; enddo ; enddo
call post_data(CS%id_diffu_visc_rem, diffu_visc_rem, CS%diag)
deallocate(diffu_visc_rem)
endif
if (present(ADp) .and. (CS%id_diffv_visc_rem > 0)) then
allocate(diffv_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
diffv_visc_rem(:,:,:) = 0.0
do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
diffv_visc_rem(i,J,k) = diffv(i,J,k) * ADp%visc_rem_v(i,J,k)
enddo ; enddo ; enddo
call post_data(CS%id_diffv_visc_rem, diffv_visc_rem, CS%diag)
deallocate(diffv_visc_rem)
endif

end subroutine horizontal_viscosity

!> Allocates space for and calculates static variables used by horizontal_viscosity().
Expand Down Expand Up @@ -2448,6 +2470,20 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp)
call safe_alloc_ptr(ADp%diag_hv,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)
endif

CS%id_diffu_visc_rem = register_diag_field('ocean_model', 'diffu_visc_rem', diag%axesCuL, Time, &
'Zonal Acceleration from Horizontal Viscosity multiplied by viscous remnant', 'm2 s-2', &
conversion=GV%H_to_m*US%L_T2_to_m_s2)
if ((CS%id_diffu_visc_rem > 0) .and. (present(ADp))) then
call safe_alloc_ptr(ADp%visc_rem_u,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke)
endif

CS%id_diffv_visc_rem = register_diag_field('ocean_model', 'diffv_visc_rem', diag%axesCvL, Time, &
'Meridional Acceleration from Horizontal Viscosity multiplied by viscous remnant', 'm2 s-2', &
conversion=GV%H_to_m*US%L_T2_to_m_s2)
if ((CS%id_diffv_visc_rem > 0) .and. (present(ADp))) then
call safe_alloc_ptr(ADp%visc_rem_v,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)
endif

if (CS%biharmonic) then
CS%id_Ah_h = register_diag_field('ocean_model', 'Ahh', diag%axesTL, Time, &
'Biharmonic Horizontal Viscosity at h Points', 'm4 s-1', conversion=US%L_to_m**4*US%s_to_T, &
Expand Down
44 changes: 44 additions & 0 deletions src/parameterizations/vertical/MOM_vert_friction.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -127,6 +127,7 @@ module MOM_vert_friction
! integer :: id_hf_du_dt_visc = -1, id_hf_dv_dt_visc = -1
integer :: id_h_du_dt_visc = -1, id_h_dv_dt_visc = -1
integer :: id_hf_du_dt_visc_2d = -1, id_hf_dv_dt_visc_2d = -1
integer :: id_du_dt_visc_rem = -1, id_dv_dt_visc_rem = -1
!>@}

type(PointAccel_CS), pointer :: PointAccel_CSp => NULL() !< A pointer to the control structure
Expand Down Expand Up @@ -217,6 +218,11 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
real, allocatable, dimension(:,:,:) :: h_du_dt_visc ! h x du_dt_visc [H L T-2 ~> m2 s-2]
real, allocatable, dimension(:,:,:) :: h_dv_dt_visc ! h x dv_dt_visc [H L T-2 ~> m2 s-2]

real, allocatable, dimension(:,:,:) :: du_dt_visc_rem ! du_dt_visc x visc_rem_u + du_dt x (1-visc_rem_u)
! [L T-2 ~> m2 s-2]
real, allocatable, dimension(:,:,:) :: dv_dt_visc_rem ! dv_dt_visc x visc_rem_v + dv_dt x (1-visc_rem_v)
! [L T-2 ~> m2 s-2]

logical :: do_i(SZIB_(G))
logical :: DoStokesMixing

Expand Down Expand Up @@ -522,6 +528,27 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
deallocate(h_dv_dt_visc)
endif

if (CS%id_du_dt_visc_rem > 0) then
allocate(du_dt_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
du_dt_visc_rem(:,:,:) = 0.0
do k=1,nz ; do j=js,je ; do I=Isq,Ieq
du_dt_visc_rem(I,j,k) = ADp%du_dt_visc(I,j,k) * ADp%visc_rem_u(I,j,k) + &
(1-ADp%visc_rem_u(I,j,k)) * ADp%du_dt(I,j,k)
enddo ; enddo ; enddo
call post_data(CS%id_du_dt_visc_rem, du_dt_visc_rem, CS%diag)
deallocate(du_dt_visc_rem)
endif
if (CS%id_dv_dt_visc_rem > 0) then
allocate(dv_dt_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
dv_dt_visc_rem(:,:,:) = 0.0
do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
dv_dt_visc_rem(i,J,k) = ADp%dv_dt_visc(i,J,k) * ADp%visc_rem_v(i,J,k) + &
(1-ADp%visc_rem_v(i,J,k)) * ADp%dv_dt(i,J,k)
enddo ; enddo ; enddo
call post_data(CS%id_dv_dt_visc_rem, dv_dt_visc_rem, CS%diag)
deallocate(dv_dt_visc_rem)
endif

end subroutine vertvisc

!> Calculate the fraction of momentum originally in a layer that remains
Expand Down Expand Up @@ -1877,6 +1904,23 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
call safe_alloc_ptr(ADp%diag_hv,isd,ied,JsdB,JedB,nz)
endif

CS%id_du_dt_visc_rem = register_diag_field('ocean_model', 'du_dt_visc_rem', diag%axesCuL, Time, &
'Zonal Acceleration from Horizontal Viscosity multiplied by viscous remnant fraction', 'm2 s-2', &
conversion=GV%H_to_m*US%L_T2_to_m_s2)
if (CS%id_du_dt_visc_rem > 0) then
call safe_alloc_ptr(ADp%du_dt,IsdB,IedB,jsd,jed,nz)
call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)
call safe_alloc_ptr(ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz)
endif
CS%id_dv_dt_visc_rem = register_diag_field('ocean_model', 'dv_dt_visc_rem', diag%axesCvL, Time, &
'Meridional Acceleration from Horizontal Viscosity multiplied by viscous remnant fraction', 'm2 s-2', &
conversion=GV%H_to_m*US%L_T2_to_m_s2)
if (CS%id_dv_dt_visc_rem > 0) then
call safe_alloc_ptr(ADp%dv_dt,isd,ied,JsdB,JedB,nz)
call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
call safe_alloc_ptr(ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
endif

if ((len_trim(CS%u_trunc_file) > 0) .or. (len_trim(CS%v_trunc_file) > 0)) &
call PointAccel_init(MIS, Time, G, param_file, diag, dirs, CS%PointAccel_CSp)

Expand Down