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

User/bgr/kvisc fix #267

Merged
merged 2 commits into from
Apr 21, 2016
+77 −70
Merged

User/bgr/kvisc fix #267

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ee9c682
Disable diffusion of momentum by diabatic_driver in ALE mode
breichl Apr 20, 2016
5726297
Missing factor of 2 in application of KV_turb
breichl Apr 20, 2016
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
137 changes: 69 additions & 68 deletions src/parameterizations/vertical/MOM_diabatic_driver.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1240,84 +1240,85 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS)
call do_group_pass(CS%pass_hold_eb_ea,G%Domain)
call cpu_clock_end(id_clock_pass)


! Use a tridiagonal solver to determine effect of the diapycnal
! advection on velocity field. It is assumed that water leaves
! or enters the ocean with the surface velocity.
if (CS%debug) then
call MOM_state_chksum("before u/v tridiag ", u, v, h, G, GV)
call hchksum(ea, "before u/v tridiag ea",G)
call hchksum(eb, "before u/v tridiag eb",G)
call hchksum(hold, "before u/v tridiag hold",G)
endif
call cpu_clock_begin(id_clock_tridiag)
if (.not. CS%useALEalgorithm) then
! Use a tridiagonal solver to determine effect of the diapycnal
! advection on velocity field. It is assumed that water leaves
! or enters the ocean with the surface velocity.
if (CS%debug) then
call MOM_state_chksum("before u/v tridiag ", u, v, h, G, GV)
call hchksum(ea, "before u/v tridiag ea",G)
call hchksum(eb, "before u/v tridiag eb",G)
call hchksum(hold, "before u/v tridiag hold",G)
endif
call cpu_clock_begin(id_clock_tridiag)
!$OMP parallel do default(none) shared(js,je,Isq,Ieq,ADp,u,hold,ea,h_neglect,eb,nz,Idt) &
!$OMP private(hval,b1,d1,c1,eaval)
do j=js,je
do I=Isq,Ieq
if (ASSOCIATED(ADp%du_dt_dia)) ADp%du_dt_dia(I,j,1) = u(I,j,1)
hval = (hold(i,j,1) + hold(i+1,j,1)) + (ea(i,j,1) + ea(i+1,j,1)) + h_neglect
b1(I) = 1.0 / (hval + (eb(i,j,1) + eb(i+1,j,1)))
d1(I) = hval * b1(I)
u(I,j,1) = b1(I) * (hval * u(I,j,1))
enddo
do k=2,nz ; do I=Isq,Ieq
if (ASSOCIATED(ADp%du_dt_dia)) ADp%du_dt_dia(I,j,k) = u(I,j,k)
c1(I,k) = (eb(i,j,k-1)+eb(i+1,j,k-1)) * b1(I)
eaval = ea(i,j,k) + ea(i+1,j,k)
hval = hold(i,j,k) + hold(i+1,j,k) + h_neglect
b1(I) = 1.0 / ((eb(i,j,k) + eb(i+1,j,k)) + (hval + d1(I)*eaval))
d1(I) = (hval + d1(I)*eaval) * b1(I)
u(I,j,k) = (hval*u(I,j,k) + eaval*u(I,j,k-1))*b1(I)
enddo ; enddo
do k=nz-1,1,-1 ; do I=Isq,Ieq
u(I,j,k) = u(I,j,k) + c1(I,k+1)*u(I,j,k+1)
if (ASSOCIATED(ADp%du_dt_dia)) &
ADp%du_dt_dia(I,j,k) = (u(I,j,k) - ADp%du_dt_dia(I,j,k)) * Idt
enddo ; enddo
if (ASSOCIATED(ADp%du_dt_dia)) then
do j=js,je
do I=Isq,Ieq
ADp%du_dt_dia(I,j,nz) = (u(I,j,nz)-ADp%du_dt_dia(I,j,nz)) * Idt
if (ASSOCIATED(ADp%du_dt_dia)) ADp%du_dt_dia(I,j,1) = u(I,j,1)
hval = (hold(i,j,1) + hold(i+1,j,1)) + (ea(i,j,1) + ea(i+1,j,1)) + h_neglect
b1(I) = 1.0 / (hval + (eb(i,j,1) + eb(i+1,j,1)))
d1(I) = hval * b1(I)
u(I,j,1) = b1(I) * (hval * u(I,j,1))
enddo
do k=2,nz ; do I=Isq,Ieq
if (ASSOCIATED(ADp%du_dt_dia)) ADp%du_dt_dia(I,j,k) = u(I,j,k)
c1(I,k) = (eb(i,j,k-1)+eb(i+1,j,k-1)) * b1(I)
eaval = ea(i,j,k) + ea(i+1,j,k)
hval = hold(i,j,k) + hold(i+1,j,k) + h_neglect
b1(I) = 1.0 / ((eb(i,j,k) + eb(i+1,j,k)) + (hval + d1(I)*eaval))
d1(I) = (hval + d1(I)*eaval) * b1(I)
u(I,j,k) = (hval*u(I,j,k) + eaval*u(I,j,k-1))*b1(I)
enddo ; enddo
do k=nz-1,1,-1 ; do I=Isq,Ieq
u(I,j,k) = u(I,j,k) + c1(I,k+1)*u(I,j,k+1)
if (ASSOCIATED(ADp%du_dt_dia)) &
ADp%du_dt_dia(I,j,k) = (u(I,j,k) - ADp%du_dt_dia(I,j,k)) * Idt
enddo ; enddo
if (ASSOCIATED(ADp%du_dt_dia)) then
do I=Isq,Ieq
ADp%du_dt_dia(I,j,nz) = (u(I,j,nz)-ADp%du_dt_dia(I,j,nz)) * Idt
enddo
endif
enddo
if (CS%debug) then
call MOM_state_chksum("aft 1st loop tridiag ", u, v, h, G, GV)
endif
enddo
if (CS%debug) then
call MOM_state_chksum("aft 1st loop tridiag ", u, v, h, G, GV)
endif
!$OMP parallel do default(none) shared(Jsq,Jeq,is,ie,ADp,v,hold,ea,h_neglect,eb,nz,Idt) &
!$OMP private(hval,b1,d1,c1,eaval)
do J=Jsq,Jeq
do i=is,ie
if (ASSOCIATED(ADp%dv_dt_dia)) ADp%dv_dt_dia(i,J,1) = v(i,J,1)
hval = (hold(i,j,1) + hold(i,j+1,1)) + (ea(i,j,1) + ea(i,j+1,1)) + h_neglect
b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i,j+1,1)))
d1(I) = hval * b1(I)
v(i,J,1) = b1(i) * (hval * v(i,J,1))
enddo
do k=2,nz ; do i=is,ie
if (ASSOCIATED(ADp%dv_dt_dia)) ADp%dv_dt_dia(i,J,k) = v(i,J,k)
c1(i,k) = (eb(i,j,k-1)+eb(i,j+1,k-1)) * b1(i)
eaval = ea(i,j,k) + ea(i,j+1,k)
hval = hold(i,j,k) + hold(i,j+1,k) + h_neglect
b1(i) = 1.0 / ((eb(i,j,k) + eb(i,j+1,k)) + (hval + d1(i)*eaval))
d1(i) = (hval + d1(i)*eaval) * b1(i)
v(i,J,k) = (hval*v(i,J,k) + eaval*v(i,J,k-1))*b1(i)
enddo ; enddo
do k=nz-1,1,-1 ; do i=is,ie
v(i,J,k) = v(i,J,k) + c1(i,k+1)*v(i,J,k+1)
if (ASSOCIATED(ADp%dv_dt_dia)) &
ADp%dv_dt_dia(i,J,k) = (v(i,J,k) - ADp%dv_dt_dia(i,J,k)) * Idt
enddo ; enddo
if (ASSOCIATED(ADp%dv_dt_dia)) then
do J=Jsq,Jeq
do i=is,ie
ADp%dv_dt_dia(i,J,nz) = (v(i,J,nz)-ADp%dv_dt_dia(i,J,nz)) * Idt
if (ASSOCIATED(ADp%dv_dt_dia)) ADp%dv_dt_dia(i,J,1) = v(i,J,1)
hval = (hold(i,j,1) + hold(i,j+1,1)) + (ea(i,j,1) + ea(i,j+1,1)) + h_neglect
b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i,j+1,1)))
d1(I) = hval * b1(I)
v(i,J,1) = b1(i) * (hval * v(i,J,1))
enddo
do k=2,nz ; do i=is,ie
if (ASSOCIATED(ADp%dv_dt_dia)) ADp%dv_dt_dia(i,J,k) = v(i,J,k)
c1(i,k) = (eb(i,j,k-1)+eb(i,j+1,k-1)) * b1(i)
eaval = ea(i,j,k) + ea(i,j+1,k)
hval = hold(i,j,k) + hold(i,j+1,k) + h_neglect
b1(i) = 1.0 / ((eb(i,j,k) + eb(i,j+1,k)) + (hval + d1(i)*eaval))
d1(i) = (hval + d1(i)*eaval) * b1(i)
v(i,J,k) = (hval*v(i,J,k) + eaval*v(i,J,k-1))*b1(i)
enddo ; enddo
do k=nz-1,1,-1 ; do i=is,ie
v(i,J,k) = v(i,J,k) + c1(i,k+1)*v(i,J,k+1)
if (ASSOCIATED(ADp%dv_dt_dia)) &
ADp%dv_dt_dia(i,J,k) = (v(i,J,k) - ADp%dv_dt_dia(i,J,k)) * Idt
enddo ; enddo
if (ASSOCIATED(ADp%dv_dt_dia)) then
do i=is,ie
ADp%dv_dt_dia(i,J,nz) = (v(i,J,nz)-ADp%dv_dt_dia(i,J,nz)) * Idt
enddo
endif
enddo
call cpu_clock_end(id_clock_tridiag)
if (CS%debug) then
call MOM_state_chksum("after u/v tridiag ", u, v, h, G, GV)
endif
enddo
call cpu_clock_end(id_clock_tridiag)
if (CS%debug) then
call MOM_state_chksum("after u/v tridiag ", u, v, h, G, GV)
endif
endif ! useALEalgorithm

! Frazil formation keeps temperature above the freezing point.
! make_frazil is deliberately called at both the beginning and at
Expand Down
10 changes: 8 additions & 2 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 @@ -995,6 +995,8 @@ subroutine find_coupling_coef(a, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_m
integer :: nz
real :: botfn

a(:,:) = 0.0

if (work_on_u) then ; is = G%IscB ; ie = G%IecB
else ; is = G%isc ; ie = G%iec ; endif
nz = G%ke
Expand Down Expand Up @@ -1034,13 +1036,17 @@ subroutine find_coupling_coef(a, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_m
endif ; enddo

if (associated(visc%Kv_turb)) then
! BGR/ Add factor of 2. * the averaged Kv_turb.
! this is needed to reproduce the analytical solution to
! a simple diffusion problem, likely due to h_shear being
! equal to 2 x \delta z
if (work_on_u) then
do K=nz,2,-1 ; do i=is,ie ; if (do_i(i)) then
a(i,K) = a(i,K) + 0.5*(visc%Kv_turb(i,j,k) + visc%Kv_turb(i+1,j,k))
a(i,K) = a(i,K) + (2.*0.5)*(visc%Kv_turb(i,j,k) + visc%Kv_turb(i+1,j,k))
endif ; enddo ; enddo
else
do K=nz,2,-1 ; do i=is,ie ; if (do_i(i)) then
a(i,K) = a(i,K) + 0.5*(visc%Kv_turb(i,j,k) + visc%Kv_turb(i,j+1,k))
a(i,K) = a(i,K) + (2.*0.5)*(visc%Kv_turb(i,j,k) + visc%Kv_turb(i,j+1,k))
endif ; enddo ; enddo
endif
endif
Expand Down