Commit d51b4a55 authored by bjb's avatar bjb
Browse files

updated names for temporal relax bdy scheme

parent 7cafd2d3
!$Id: bdy_3d.F90,v 1.13 2009-09-04 07:58:19 bjb Exp $
!$Id: bdy_3d.F90,v 1.14 2009-09-04 11:50:34 bjb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -26,13 +26,14 @@
! !PUBLIC DATA MEMBERS:
public init_bdy_3d, do_bdy_3d
REALTYPE, public, allocatable :: S_bdy(:,:),T_bdy(:,:)
REALTYPE, public :: bdy3d_rlxucut=_ONE_/50
REALTYPE, public :: bdy3d_rlxmax=_ONE_/4
REALTYPE, public :: bdy3d_rlxmin=_ZERO_
REALTYPE, allocatable :: bdyvertS(:), bdyvertT(:)
REALTYPE, allocatable :: rlxcoef(:)
logical, public :: bdy3d_tmrlx=.false.
REALTYPE, public :: bdy3d_tmrlx_ucut=_ONE_/50
REALTYPE, public :: bdy3d_tmrlx_max=_ONE_/4
REALTYPE, public :: bdy3d_tmrlx_min=_ZERO_
!
! !PRIVATE DATA MEMBERS:
REALTYPE, allocatable :: bdyvertS(:), bdyvertT(:)
REALTYPE, allocatable :: rlxcoef(:)
!
! !REVISION HISTORY:
! Original author(s): Karsten Bolding & Hans Burchard
......@@ -139,7 +140,7 @@
! !LOCAL VARIABLES:
integer :: i,j,k,l,n,ii,jj,kk
REALTYPE :: sp(1:4),rat
REALTYPE :: bdy3d_rlxumin
REALTYPE :: bdy3d_tmrlx_umin
REALTYPE :: wsum
!EOP
!-----------------------------------------------------------------------
......@@ -153,7 +154,7 @@
#if 0
select case (tag)
case (1)
! Lateral zero-gradient boundary condition (north & south)
! Lateral zero-gradient boundary condition (north & south)
do k=1,kmax
do i=imin,imax
if (au(i,jmin) .eq. 3) field(i,jmin,k)=field(i,jmin+1,k)
......@@ -161,7 +162,7 @@
end do
end do
case (2)
! Lateral zero-gradient boundary conditions (west & east)
! Lateral zero-gradient boundary conditions (west & east)
do k=1,kmax
do j=jmin,jmax
if (av(imin,j) .eq. 3) field(imin,j,k)=field(imin+1,j,k)
......@@ -181,8 +182,8 @@
sp(3)=0.25
sp(4)=0.0625
! Hardcoding of lower limit of velocity cut-off for temporal relaxation.
! Linear variation between bdy3d_rlxumin and bdy3d_rlxucut.
bdy3d_rlxumin = -0.25*bdy3d_rlxucut
! Linear variation between bdy3d_tmrlx_umin and bdy3d_tmrlx_ucut.
bdy3d_tmrlx_umin = -_QUART_*bdy3d_tmrlx_ucut
l = 0
k = 0
......@@ -191,44 +192,58 @@
k = bdy_index(l)
i = wi(n)
do j=wfj(n),wlj(n)
if (au(i,j).gt.0) then
! Local temporal relaxation coeficient depends on
! local current just *inside* domain:
do kk=1,kmax
if (uu(i,j,kk).ge.bdy3d_rlxucut) then
rlxcoef(kk) = bdy3d_rlxmax
else if (uu(i,j,kk).le.bdy3d_rlxumin) then
rlxcoef(kk) = bdy3d_rlxmin
else
rlxcoef(kk) = (bdy3d_rlxmax-bdy3d_rlxmin)*(uu(i,j,kk)-bdy3d_rlxumin) &
/(bdy3d_rlxucut-bdy3d_rlxumin)+ bdy3d_rlxmin
end if
end do
else
rlxcoef(:)=bdy3d_rlxmax
end if
! Temporal relaxation: Weight inner (actual) solution near boundary
! with boundary condition (outer solution.)
wsum= MIN(az(i-1+2,j),1)*sp(2)+MIN(az(i-1+3,j),1)*sp(3)+MIN(az(i-1+4,j),1)*sp(4)
if (wsum>_ZERO_) then
! Get (weighted avr of) inner near-bdy solution (sponge) cells:
bdyvertS(:) = (_ONE_/wsum)*( MIN(az(i-1+2,j),1) * sp(2) * S(i-1+2,j,:) &
+MIN(az(i-1+3,j),1) * sp(3) * S(i-1+3,j,:) &
+MIN(az(i-1+4,j),1) * sp(4) * S(i-1+4,j,:) )
bdyvertT(:) = (_ONE_/wsum)*( MIN(az(i-1+2,j),1) * sp(2) * T(i-1+2,j,:) &
+MIN(az(i-1+3,j),1) * sp(3) * T(i-1+3,j,:) &
+MIN(az(i-1+4,j),1) * sp(4) * T(i-1+4,j,:) )
! Weight inner and outer (bc) solutions for use
! in spatial relaxation/sponge
bdyvertS(:) = (_ONE_-rlxcoef(:))*bdyvertS(:) + rlxcoef(:)*S_bdy(:,k)
bdyvertT(:) = (_ONE_-rlxcoef(:))*bdyvertT(:) + rlxcoef(:)*T_bdy(:,k)
if (bdy3d_tmrlx) then
if (au(i,j).gt.0) then
! Local temporal relaxation coeficient depends on
! local current just *inside* domain:
do kk=1,kmax
if (uu(i,j,kk).ge.bdy3d_tmrlx_ucut) then
rlxcoef(kk) = bdy3d_tmrlx_max
else if (uu(i,j,kk).le.bdy3d_tmrlx_umin) then
rlxcoef(kk) = bdy3d_tmrlx_min
else
rlxcoef(kk) = (bdy3d_tmrlx_max-bdy3d_tmrlx_min) &
*(uu(i,j,kk)-bdy3d_tmrlx_umin) &
/(bdy3d_tmrlx_ucut-bdy3d_tmrlx_umin) &
+ bdy3d_tmrlx_min
end if
end do
else
rlxcoef(:)=bdy3d_tmrlx_max
end if
! Temporal relaxation: Weight inner (actual) solution near boundary
! with boundary condition (outer solution.)
wsum= MIN(az(i-1+2,j),1)*sp(2) &
+MIN(az(i-1+3,j),1)*sp(3) &
+MIN(az(i-1+4,j),1)*sp(4)
if (wsum>_ZERO_) then
! Get (weighted avr of) inner near-bdy solution (sponge) cells:
bdyvertS(:) = (_ONE_/wsum)*( &
MIN(az(i-1+2,j),1) * sp(2) * S(i-1+2,j,:) &
+MIN(az(i-1+3,j),1) * sp(3) * S(i-1+3,j,:) &
+MIN(az(i-1+4,j),1) * sp(4) * S(i-1+4,j,:) &
)
bdyvertT(:) = (_ONE_/wsum)*( &
MIN(az(i-1+2,j),1) * sp(2) * T(i-1+2,j,:) &
+MIN(az(i-1+3,j),1) * sp(3) * T(i-1+3,j,:) &
+MIN(az(i-1+4,j),1) * sp(4) * T(i-1+4,j,:) &
)
! Weight inner and outer (bc) solutions for use
! in spatial relaxation/sponge
bdyvertS(:) = (_ONE_-rlxcoef(:))*bdyvertS(:) + rlxcoef(:)*S_bdy(:,k)
bdyvertT(:) = (_ONE_-rlxcoef(:))*bdyvertT(:) + rlxcoef(:)*T_bdy(:,k)
else
! No near-bdy points. Just clamp bdy temporally:
bdyvertS(:) = S_bdy(:,k)
bdyvertT(:) = T_bdy(:,k)
end if
else
! No near-bdy points. Just clamp bdy temporally:
! No time-relaxation. Just clamp at bondary points.
bdyvertS(:) = S_bdy(:,k)
bdyvertT(:) = T_bdy(:,k)
end if
do ii=1,4
! Spatial relaxation / sponge:
! Spatial relaxation / sponge:
if (az(i-1+ii,j).gt.0) then
S(i-1+ii,j,:) = sp(ii)*bdyvertS(:)+(_ONE_-sp(ii))*S(i-1+ii,j,:)
T(i-1+ii,j,:) = sp(ii)*bdyvertT(:)+(_ONE_-sp(ii))*T(i-1+ii,j,:)
......@@ -248,30 +263,43 @@
k = bdy_index(l)
j = nj(n)
do i = nfi(n),nli(n)
if (av(i,j-1).gt.0) then
do kk=1,kmax
if (vv(i,j-1,kk).le.-bdy3d_rlxucut) then
rlxcoef(kk) = bdy3d_rlxmax
else if (vv(i,j-1,kk).ge.-bdy3d_rlxumin) then
rlxcoef(kk) = bdy3d_rlxmin
else
rlxcoef(kk) = -(bdy3d_rlxmax-bdy3d_rlxmin)*(vv(i,j-1,kk)+bdy3d_rlxumin) &
/(bdy3d_rlxucut-bdy3d_rlxumin)+ bdy3d_rlxmin
end if
end do
else
rlxcoef(:)=bdy3d_rlxmax
end if
wsum= MIN(az(i,j+1-2),1)*sp(2)+MIN(az(i,j+1-3),1)*sp(3)+MIN(az(i,j+1-4),1)*sp(4)
if (wsum>_ZERO_) then
bdyvertS(:) = (_ONE_/wsum)*( MIN(az(i,j+1-2),1) * sp(2) * S(i,j+1-2,:) &
+MIN(az(i,j+1-3),1) * sp(3) * S(i,j+1-3,:) &
+MIN(az(i,j+1-4),1) * sp(4) * S(i,j+1-4,:) )
bdyvertT(:) = (_ONE_/wsum)*( MIN(az(i,j+1-2),1) * sp(2) * T(i,j+1-2,:) &
+MIN(az(i,j+1-3),1) * sp(3) * T(i,j+1-3,:) &
+MIN(az(i,j+1-4),1) * sp(4) * T(i,j+1-4,:) )
bdyvertS(:) = (_ONE_-rlxcoef(:))*bdyvertS(:) + rlxcoef(:)*S_bdy(:,k)
bdyvertT(:) = (_ONE_-rlxcoef(:))*bdyvertT(:) + rlxcoef(:)*T_bdy(:,k)
if (bdy3d_tmrlx) then
if (av(i,j-1).gt.0) then
do kk=1,kmax
if (vv(i,j-1,kk).le.-bdy3d_tmrlx_ucut) then
rlxcoef(kk) = bdy3d_tmrlx_max
else if (vv(i,j-1,kk).ge.-bdy3d_tmrlx_umin) then
rlxcoef(kk) = bdy3d_tmrlx_min
else
rlxcoef(kk) = -(bdy3d_tmrlx_max-bdy3d_tmrlx_min) &
*(vv(i,j-1,kk)+bdy3d_tmrlx_umin) &
/(bdy3d_tmrlx_ucut-bdy3d_tmrlx_umin) &
+ bdy3d_tmrlx_min
end if
end do
else
rlxcoef(:)=bdy3d_tmrlx_max
end if
wsum= MIN(az(i,j+1-2),1)*sp(2) &
+MIN(az(i,j+1-3),1)*sp(3) &
+MIN(az(i,j+1-4),1)*sp(4)
if (wsum>_ZERO_) then
bdyvertS(:) = (_ONE_/wsum)*( &
MIN(az(i,j+1-2),1) * sp(2) * S(i,j+1-2,:) &
+MIN(az(i,j+1-3),1) * sp(3) * S(i,j+1-3,:) &
+MIN(az(i,j+1-4),1) * sp(4) * S(i,j+1-4,:) &
)
bdyvertT(:) = (_ONE_/wsum)*( &
MIN(az(i,j+1-2),1) * sp(2) * T(i,j+1-2,:) &
+MIN(az(i,j+1-3),1) * sp(3) * T(i,j+1-3,:) &
+MIN(az(i,j+1-4),1) * sp(4) * T(i,j+1-4,:) &
)
bdyvertS(:) = (_ONE_-rlxcoef(:))*bdyvertS(:) + rlxcoef(:)*S_bdy(:,k)
bdyvertT(:) = (_ONE_-rlxcoef(:))*bdyvertT(:) + rlxcoef(:)*T_bdy(:,k)
else
bdyvertS(:) = S_bdy(:,k)
bdyvertT(:) = T_bdy(:,k)
end if
else
bdyvertS(:) = S_bdy(:,k)
bdyvertT(:) = T_bdy(:,k)
......@@ -296,30 +324,43 @@
k = bdy_index(l)
i = ei(n)
do j=efj(n),elj(n)
if (au(i-1,j).gt.0) then
do kk=1,kmax
if (uu(i-1,j,kk).le.-bdy3d_rlxucut) then
rlxcoef(kk) = bdy3d_rlxmax
else if (uu(i-1,j,kk).ge.-bdy3d_rlxumin) then
rlxcoef(kk) = bdy3d_rlxmin
else
rlxcoef(kk) = -(bdy3d_rlxmax-bdy3d_rlxmin)*(uu(i-1,j,kk)+bdy3d_rlxumin) &
/(bdy3d_rlxucut-bdy3d_rlxumin)+ bdy3d_rlxmin
end if
end do
else
rlxcoef(:)=bdy3d_rlxmax
end if
wsum= MIN(az(i+1-2,j),1)*sp(2)+MIN(az(i+1-3,j),1)*sp(3)+MIN(az(i+1-4,j),1)*sp(4)
if (wsum>_ZERO_) then
bdyvertS(:) = (_ONE_/wsum)*( MIN(az(i+1-2,j),1) * sp(2) * S(i+1-2,j,:) &
+MIN(az(i+1-3,j),1) * sp(3) * S(i+1-3,j,:) &
+MIN(az(i+1-4,j),1) * sp(4) * S(i+1-4,j,:) )
bdyvertT(:) = (_ONE_/wsum)*( MIN(az(i+1-2,j),1) * sp(2) * T(i+1-2,j,:) &
+MIN(az(i+1-3,j),1) * sp(3) * T(i+1-3,j,:) &
+MIN(az(i+1-4,j),1) * sp(4) * T(i+1-4,j,:) )
bdyvertS(:) = (_ONE_-rlxcoef(:))*bdyvertS(:) + rlxcoef(:)*S_bdy(:,k)
bdyvertT(:) = (_ONE_-rlxcoef(:))*bdyvertT(:) + rlxcoef(:)*T_bdy(:,k)
if (bdy3d_tmrlx) then
if (au(i-1,j).gt.0) then
do kk=1,kmax
if (uu(i-1,j,kk).le.-bdy3d_tmrlx_ucut) then
rlxcoef(kk) = bdy3d_tmrlx_max
else if (uu(i-1,j,kk).ge.-bdy3d_tmrlx_umin) then
rlxcoef(kk) = bdy3d_tmrlx_min
else
rlxcoef(kk) = -(bdy3d_tmrlx_max-bdy3d_tmrlx_min) &
*(uu(i-1,j,kk)+bdy3d_tmrlx_umin) &
/(bdy3d_tmrlx_ucut-bdy3d_tmrlx_umin) &
+ bdy3d_tmrlx_min
end if
end do
else
rlxcoef(:)=bdy3d_tmrlx_max
end if
wsum= MIN(az(i+1-2,j),1)*sp(2) &
+MIN(az(i+1-3,j),1)*sp(3) &
+MIN(az(i+1-4,j),1)*sp(4)
if (wsum>_ZERO_) then
bdyvertS(:) = (_ONE_/wsum)*( &
MIN(az(i+1-2,j),1) * sp(2) * S(i+1-2,j,:) &
+MIN(az(i+1-3,j),1) * sp(3) * S(i+1-3,j,:) &
+MIN(az(i+1-4,j),1) * sp(4) * S(i+1-4,j,:) &
)
bdyvertT(:) = (_ONE_/wsum)*( &
MIN(az(i+1-2,j),1) * sp(2) * T(i+1-2,j,:) &
+MIN(az(i+1-3,j),1) * sp(3) * T(i+1-3,j,:) &
+MIN(az(i+1-4,j),1) * sp(4) * T(i+1-4,j,:) &
)
bdyvertS(:) = (_ONE_-rlxcoef(:))*bdyvertS(:) + rlxcoef(:)*S_bdy(:,k)
bdyvertT(:) = (_ONE_-rlxcoef(:))*bdyvertT(:) + rlxcoef(:)*T_bdy(:,k)
else
bdyvertS(:) = S_bdy(:,k)
bdyvertT(:) = T_bdy(:,k)
end if
else
bdyvertS(:) = S_bdy(:,k)
bdyvertT(:) = T_bdy(:,k)
......@@ -345,29 +386,42 @@
j = sj(n)
do i = sfi(n),sli(n)
if (av(i,j).gt.0) then
do kk=1,kmax
if (vv(i,j,kk).ge.bdy3d_rlxucut) then
rlxcoef(kk) = bdy3d_rlxmax
else if (vv(i,j,kk).le.bdy3d_rlxumin) then
rlxcoef(kk) = bdy3d_rlxmin
else
rlxcoef(kk) = (bdy3d_rlxmax-bdy3d_rlxmin)*(vv(i,j,kk)-bdy3d_rlxumin) &
/(bdy3d_rlxucut-bdy3d_rlxumin)+ bdy3d_rlxmin
end if
end do
else
rlxcoef(:)=bdy3d_rlxmax
end if
wsum= MIN(az(i,j-1+2),1)*sp(2)+MIN(az(i,j-1+3),1)*sp(3)+MIN(az(i,j-1+4),1)*sp(4)
if (wsum>_ZERO_) then
bdyvertS(:) = (_ONE_/wsum)*( MIN(az(i,j-1+2),1) * sp(2) * S(i,j-1+2,:) &
+MIN(az(i,j-1+3),1) * sp(3) * S(i,j-1+3,:) &
+MIN(az(i,j-1+4),1) * sp(4) * S(i,j-1+4,:) )
bdyvertT(:) = (_ONE_/wsum)*( MIN(az(i,j-1+2),1) * sp(2) * T(i,j-1+2,:) &
+MIN(az(i,j-1+3),1) * sp(3) * T(i,j-1+3,:) &
+MIN(az(i,j-1+4),1) * sp(4) * T(i,j-1+4,:) )
bdyvertS(:) = (_ONE_-rlxcoef(:))*bdyvertS(:) + rlxcoef(:)*S_bdy(:,k)
bdyvertT(:) = (_ONE_-rlxcoef(:))*bdyvertT(:) + rlxcoef(:)*T_bdy(:,k)
if (bdy3d_tmrlx) then
do kk=1,kmax
if (vv(i,j,kk).ge.bdy3d_tmrlx_ucut) then
rlxcoef(kk) = bdy3d_tmrlx_max
else if (vv(i,j,kk).le.bdy3d_tmrlx_umin) then
rlxcoef(kk) = bdy3d_tmrlx_min
else
rlxcoef(kk) = (bdy3d_tmrlx_max-bdy3d_tmrlx_min) &
*(vv(i,j,kk)-bdy3d_tmrlx_umin) &
/(bdy3d_tmrlx_ucut-bdy3d_tmrlx_umin) &
+ bdy3d_tmrlx_min
end if
end do
else
rlxcoef(:)=bdy3d_tmrlx_max
end if
wsum= MIN(az(i,j-1+2),1)*sp(2) &
+MIN(az(i,j-1+3),1)*sp(3) &
+MIN(az(i,j-1+4),1)*sp(4)
if (wsum>_ZERO_) then
bdyvertS(:) = (_ONE_/wsum)*( &
MIN(az(i,j-1+2),1) * sp(2) * S(i,j-1+2,:) &
+MIN(az(i,j-1+3),1) * sp(3) * S(i,j-1+3,:) &
+MIN(az(i,j-1+4),1) * sp(4) * S(i,j-1+4,:) &
)
bdyvertT(:) = (_ONE_/wsum)*( &
MIN(az(i,j-1+2),1) * sp(2) * T(i,j-1+2,:) &
+MIN(az(i,j-1+3),1) * sp(3) * T(i,j-1+3,:) &
+MIN(az(i,j-1+4),1) * sp(4) * T(i,j-1+4,:) &
)
bdyvertS(:) = (_ONE_-rlxcoef(:))*bdyvertS(:) + rlxcoef(:)*S_bdy(:,k)
bdyvertT(:) = (_ONE_-rlxcoef(:))*bdyvertT(:) + rlxcoef(:)*T_bdy(:,k)
else
bdyvertS(:) = S_bdy(:,k)
bdyvertT(:) = T_bdy(:,k)
end if
else
bdyvertS(:) = S_bdy(:,k)
bdyvertT(:) = T_bdy(:,k)
......
!$Id: m3d.F90,v 1.45 2009-09-04 07:58:19 bjb Exp $
!$Id: m3d.F90,v 1.46 2009-09-04 11:50:33 bjb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -43,7 +43,7 @@
use variables_3d
use advection_3d, only: init_advection_3d
use bdy_3d, only: init_bdy_3d, do_bdy_3d
use bdy_3d, only: bdy3d_rlxucut, bdy3d_rlxmax, bdy3d_rlxmin
use bdy_3d, only: bdy3d_tmrlx, bdy3d_tmrlx_ucut, bdy3d_tmrlx_max, bdy3d_tmrlx_min
! Necessary to use halo_zones because update_3d_halos() have been moved out
! temperature.F90 and salinity.F90 - should be changed at a later stage
use halo_zones, only: update_3d_halo,wait_halo,D_TAG
......@@ -119,7 +119,8 @@
NAMELIST /m3d/ &
M,cnpar,cord_relax, &
bdy3d,bdyfmt_3d,bdyramp_3d,bdyfile_3d, &
bdy3d_rlxucut,bdy3d_rlxmax,bdy3d_rlxmin, &
bdy3d_tmrlx,bdy3d_tmrlx_ucut, &
bdy3d_tmrlx_max,bdy3d_tmrlx_min, &
vel_hor_adv,vel_ver_adv,vel_adv_split, &
calc_temp,calc_salt, &
avmback,avhback,ip_method,ip_ramp, &
......@@ -219,20 +220,20 @@
end select
! Sanity checks for bdy 3d
if (bdy3d) then
LEVEL2 'bdy3d=.true.'
LEVEL2 'bdy3d_rlxmax= ',bdy3d_rlxmax
LEVEL2 'bdy3d_rlxmin= ',bdy3d_rlxmin
LEVEL2 'bdy3d_rlxucut= ',bdy3d_rlxucut
if (bdy3d_rlxmin<_ZERO_ .or. bdy3d_rlxmin>_ONE_) &
call getm_error("init_3d()", &
"bdy3d_rlxmin is out of valid range [0:1]")
if (bdy3d_rlxmax<bdy3d_rlxmin .or. bdy3d_rlxmin>_ONE_) &
call getm_error("init_3d()", &
"bdy3d_rlxmax is out of valid range [bdy3d_rlxmax:1]")
if (bdy3d_rlxucut<_ZERO_) &
call getm_error("init_3d()", &
"bdy3d_rlxmax is out of valid range [0:inf[")
if (bdy3d .and. bdy3d_tmrlx) then
LEVEL2 'bdy3d_tmrlx=.true.'
LEVEL2 'bdy3d_tmrlx_max= ',bdy3d_tmrlx_max
LEVEL2 'bdy3d_tmrlx_min= ',bdy3d_tmrlx_min
LEVEL2 'bdy3d_tmrlx_ucut= ',bdy3d_tmrlx_ucut
if (bdy3d_tmrlx_min<_ZERO_ .or. bdy3d_tmrlx_min>_ONE_) &
call getm_error("init_3d()", &
"bdy3d_tmrlx_min is out of valid range [0:1]")
if (bdy3d_tmrlx_max<bdy3d_tmrlx_min .or. bdy3d_tmrlx_min>_ONE_) &
call getm_error("init_3d()", &
"bdy3d_tmrlx_max is out of valid range [bdy3d_tmrlx_max:1]")
if (bdy3d_tmrlx_ucut<_ZERO_) &
call getm_error("init_3d()", &
"bdy3d_tmrlx_max is out of valid range [0:inf[")
end if
LEVEL2 'ip_ramp=',ip_ramp
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment