Skip to content

Commit 144f8a3

Browse files
JanStreffingJanStreffing
authored
Merge pull request #714 from manufrn/uke_advection
Add uke advection to dynamic backscatter
2 parents 84ab1bc + 5094433 commit 144f8a3

File tree

3 files changed

+163
-8
lines changed

3 files changed

+163
-8
lines changed

‎src/MOD_DYN.F90‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -110,6 +110,7 @@ MODULE MOD_DYN
110110
! coefficients and options for opt_visc=8 (dynamic backscatter)
111111
logical :: uke_scaling = .true.
112112
real(kind=WP) :: uke_scaling_factor = 1._WP
113+
logical :: uke_advection = .false.
113114
real(kind=WP) :: rosb_dis = 1._WP
114115
integer :: smooth_back = 2
115116
integer :: smooth_dis = 2

‎src/gen_modules_backscatter.F90‎

Lines changed: 157 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ module g_backscatter
2020
!___________________________________________________________________________
2121
! allocate backscatter arrays
2222
real(kind=WP), allocatable, dimension(:,:) :: v_back
23-
real(kind=WP), allocatable, dimension(:,:) :: uke, uke_back, uke_dis, uke_dif
23+
real(kind=WP), allocatable, dimension(:,:) :: uke, uke_back, uke_dis, uke_dif, uke_adv
2424
real(kind=WP), allocatable, dimension(:,:) :: uke_rhs, uke_rhs_old
2525
real(kind=WP), allocatable, dimension(:,:) :: UV_dis_posdef_b2, UV_dis_posdef, UV_back_posdef
2626
real(kind=WP), allocatable, dimension(:,:,:):: UV_back, UV_dis
@@ -31,9 +31,10 @@ module g_backscatter
3131
!
3232
!___________________________________________________________________________
3333
! allocate/initialise backscatter arrays
34-
subroutine init_backscatter(partit, mesh)
34+
subroutine init_backscatter(dynamics, partit, mesh)
3535
implicit none
3636
integer :: elem_size
37+
type(t_dyn) , intent(inout), target :: dynamics
3738
type(t_mesh), intent(in), target :: mesh
3839
type(t_partit), intent(inout), target :: partit
3940
#include "associate_part_def.h"
@@ -66,6 +67,11 @@ subroutine init_backscatter(partit, mesh)
6667
UV_back = 0.0_WP
6768
UV_back_tend = 0.0_WP
6869
UV_total_tend = 0.0_WP
70+
71+
if (dynamics%uke_advection) then
72+
allocate(uke_adv( nl-1, elem_size))
73+
uke_adv = 0.0_WP
74+
endif
6975

7076
end subroutine init_backscatter
7177

@@ -216,6 +222,139 @@ subroutine visc_filt_dbcksc(dynamics, partit, mesh)
216222
deallocate(uuu)
217223
end subroutine visc_filt_dbcksc
218224

225+
226+
subroutine backscatter_adv(dynamics, partit, mesh)
227+
IMPLICIT NONE
228+
type(t_dyn) , intent(inout), target :: dynamics
229+
type(t_mesh), intent(in), target :: mesh
230+
type(t_partit), intent(inout), target :: partit
231+
integer :: n, nz, el1, el2, ul1
232+
integer :: nl1, nl2, nod(2), el, ed, k, nle
233+
real(kind=WP), dimension(:,:,:), pointer :: UV
234+
real(kind=WP), dimension(:,:), pointer :: Wvel_e
235+
real(kind=WP), allocatable, dimension(:) :: contr1, contr2, wuke
236+
real(kind=WP), allocatable, dimension(:,:) :: advnode
237+
#include "associate_part_def.h"
238+
#include "associate_mesh_def.h"
239+
#include "associate_part_ass.h"
240+
#include "associate_mesh_ass.h"
241+
Wvel_e =>dynamics%w_e(:,:)
242+
UV =>dynamics%UV(:,:,:)
243+
244+
allocate(contr1(1:nl-1))
245+
allocate(contr2(1:nl-1))
246+
allocate(wuke(1:nl-1))
247+
allocate(advnode(nl-1, myDim_nod2d+eDim_nod2D))
248+
249+
advnode=0.0_WP
250+
contr1=0.0_WP
251+
contr2=0.0_WP
252+
wuke=0.0_WP
253+
254+
!vertical advection
255+
do n=1,myDim_nod2d
256+
nl1 = nlevels_nod2D(n)-1
257+
258+
!find average of uke between vertical layers
259+
do k=1,nod_in_elem2D_num(n)
260+
el = nod_in_elem2D(k,n)
261+
nle = nlevels(el)-1
262+
263+
wuke(1) = wuke(1) + uke(1,el)*elem_area(el)
264+
wuke(2:nle) = wuke(2:nle) + 0.5_WP*(uke(2:nle,el)+uke(1:nle-1,el))*elem_area(el)
265+
266+
end do
267+
268+
wuke(1:nl) = wuke(1:nl1)*Wvel_e(1:nl1,n)
269+
270+
do nz=1,nl1
271+
! Here 1/3 because 1/3 of the area is related to the node
272+
advnode(nz,n) = - (wuke(nz) - wuke(nz+1) ) / (3._WP*hnode(nz,n))
273+
end do
274+
275+
advnode(nl1+1:nl-1,n) = 0._WP
276+
277+
end do
278+
279+
!horizontal advection
280+
DO ed=1, myDim_edge2D
281+
nod = edges(:,ed)
282+
el1 = edge_tri(1,ed)
283+
el2 = edge_tri(2,ed)
284+
nl1 = nlevels(el1)-1
285+
286+
!calculate fluxes throught the each edge
287+
!for the edges that are on the boundary and have only one element
288+
contr1(1:nl1) = UV(2,1:nl1,el1)*edge_cross_dxdy(1,ed) - UV(1,1:nl1,el1)*edge_cross_dxdy(2,ed)
289+
290+
! if we have the other element as well
291+
if (el2>0) then
292+
nl2 = nlevels(el2)-1
293+
294+
contr2(1:nl2) = - UV(2,1:nl2,el2)*edge_cross_dxdy(3,ed) + UV(1,1:nl2,el2)*edge_cross_dxdy(4,ed)
295+
296+
contr1(nl1+1:max(nl1,nl2)) = 0._WP
297+
contr2(nl2+1:max(nl1,nl2)) = 0._WP
298+
299+
if (nod(1) <= myDim_nod2d) then
300+
do nz=1, max(nl1,nl2)
301+
advnode(nz,nod(1)) = advnode(nz,nod(1)) + contr1(nz)*uke(nz,el1) + contr2(nz)*uke(nz,el2)
302+
end do
303+
endif
304+
305+
if (nod(2) <= myDim_nod2d) then
306+
do nz=1, max(nl1,nl2)
307+
advnode(nz,nod(2)) = advnode(nz,nod(2)) - contr1(nz)*uke(nz,el1) - contr2(nz)*uke(nz,el2)
308+
end do
309+
endif
310+
311+
else ! ed is a boundary edge, there is only the contribution from el1
312+
if (nod(1) <= myDim_nod2d) then
313+
do nz=1, nl1
314+
advnode(nz,nod(1)) = advnode(nz,nod(1)) + contr1(nz)*uke(nz,el1)
315+
end do
316+
endif
317+
! second edge node
318+
if (nod(2) <= myDim_nod2d) then
319+
do nz=1, nl1
320+
advnode(nz,nod(2)) = advnode(nz,nod(2)) - contr1(nz)*uke(nz,el1)
321+
end do
322+
endif
323+
endif
324+
end do
325+
326+
! Multiply by the segment area
327+
do n=1,myDim_nod2d
328+
nl1 = nlevels_nod2D(n)-1
329+
advnode(1:nl1,n) = advnode(1:nl1,n) * area_inv(1:nl1,n)
330+
end do
331+
332+
call exchange_nod(advnode, partit)
333+
334+
!Calculate advection as an average across the nodes
335+
do n=1,myDim_nod2d
336+
nl1 = nlevels_nod2D(n)-1
337+
advnode(1:nl1,n) = advnode(1:nl1,n) *areasvol_inv(1:nl1,n)
338+
end do
339+
340+
do el=1, myDim_elem2D
341+
nl1 = nlevels(el)-1
342+
uke_adv(1:nl1,el) = uke_adv(1:nl1,el) &
343+
+ elem_area(el)*(advnode(1:nl1,elem2D_nodes(1,el)) &
344+
+ advnode(1:nl1,elem2D_nodes(2,el)) &
345+
+ advnode(1:nl1,elem2D_nodes(3,el)))/ 3.0_WP
346+
347+
do nz=1,nlevels(el)-1
348+
uke_adv(nz,el)=dt*uke_adv(nz,el)/elem_area(el)
349+
end do
350+
end do
351+
352+
deallocate(contr1)
353+
deallocate(contr2)
354+
deallocate(wuke)
355+
deallocate(advnode)
356+
end subroutine backscatter_adv
357+
219358
!
220359
!
221360
!_______________________________________________________________________________
@@ -374,13 +513,26 @@ subroutine uke_update(dynamics, partit, mesh)
374513
uke_dis(nz,:)=uuu
375514
END DO
376515

377-
DO ed=1, myDim_elem2D
378-
DO nz=1,nlevels(ed)-1
516+
if (dynamics%uke_advection) then
517+
call backscatter_adv(dynamics, partit, mesh)
518+
call exchange_elem(uke_adv, partit)
519+
520+
DO ed=1, myDim_elem2D
521+
DO nz=1,nlevels(ed)-1
522+
uke_rhs_old(nz,ed)=uke_rhs(nz,ed)
523+
uke_rhs(nz,ed)=-uke_dis(nz,ed)-uke_back(nz,ed)+uke_dif(nz,ed)-uke_adv(nz, ed)
524+
uke(nz,ed)=uke(nz,ed)+1.5_WP*uke_rhs(nz,ed)-0.5_WP*uke_rhs_old(nz,ed)
525+
END DO
526+
END DO
527+
else
528+
DO ed=1, myDim_elem2D
529+
DO nz=1,nlevels(ed)-1
379530
uke_rhs_old(nz,ed)=uke_rhs(nz,ed)
380531
uke_rhs(nz,ed)=-uke_dis(nz,ed)-uke_back(nz,ed)+uke_dif(nz,ed)
381532
uke(nz,ed)=uke(nz,ed)+1.5_WP*uke_rhs(nz,ed)-0.5_WP*uke_rhs_old(nz,ed)
533+
END DO
382534
END DO
383-
END DO
535+
end if
384536

385537
call exchange_elem(uke, partit)
386538
deallocate(uuu)

‎src/oce_setup_step.F90‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -255,7 +255,7 @@ subroutine ocean_setup(dynamics, tracers, partit, mesh)
255255

256256
!___________________________________________________________________________
257257
! initialise arrays that are needed for backscatter_coef
258-
if(dynamics%opt_visc==8) call init_backscatter(partit, mesh)
258+
if(dynamics%opt_visc==8) call init_backscatter(dynamics, partit, mesh)
259259

260260

261261
!___________________________________________________________________________
@@ -495,6 +495,7 @@ SUBROUTINE dynamics_init(dynamics, partit, mesh)
495495
logical :: use_ivertvisc=.true.
496496
logical :: uke_scaling=.true.
497497
real(kind=WP) :: uke_scaling_factor=1._WP
498+
logical :: uke_advection=.false.
498499
real(kind=WP) :: rosb_dis=1._WP
499500
integer :: smooth_back=2
500501
integer :: smooth_dis=2
@@ -516,8 +517,8 @@ SUBROUTINE dynamics_init(dynamics, partit, mesh)
516517

517518
namelist /dynamics_visc / opt_visc, check_opt_visc, visc_gamma0, visc_gamma1, visc_gamma2, &
518519
use_ivertvisc, visc_easybsreturn, &
519-
uke_scaling, uke_scaling_factor, rosb_dis, smooth_back, &
520-
smooth_dis, smooth_back_tend, K_back, c_back
520+
uke_scaling, uke_scaling_factor, uke_advection, &
521+
rosb_dis, smooth_back, smooth_dis, smooth_back_tend, K_back, c_back
521522

522523
namelist /dynamics_general/ momadv_opt, use_freeslip, use_wsplit, wsplit_maxcfl, &
523524
ldiag_KE, AB_order, &
@@ -557,6 +558,7 @@ SUBROUTINE dynamics_init(dynamics, partit, mesh)
557558
dynamics%visc_easybsreturn = visc_easybsreturn
558559
dynamics%uke_scaling = uke_scaling
559560
dynamics%uke_scaling_factor = uke_scaling_factor
561+
dynamics%uke_advection = uke_advection
560562
dynamics%rosb_dis = rosb_dis
561563
dynamics%smooth_back = smooth_back
562564
dynamics%smooth_dis = smooth_dis

0 commit comments

Comments
 (0)