Commit 45d5bce7 authored by kb's avatar kb
Browse files

added new 2D open boundaries

parent db1ffae0
......@@ -51,6 +51,13 @@
#define ARVD1 ard1
#endif
! For 2D boundary conditions
#define ZERO_GRADIENT 1
#define SOMMERFELDT 2
#define CLAMPED 3
#define FLATHER_ELEV 4
#define FLATHER_VEL 5
! Reserved Fortran units
#define stdin 5
#define stdout 6
......@@ -64,6 +71,7 @@
#define BDYDATA 22
! Data/file formats
#define NO_DATA -1
#define ANALYTICAL 0
#define ASCII 1
#define NETCDF 2
......
!$Id: have_bdy.F90,v 1.8 2006-03-01 15:54:07 kbk Exp $
!$Id: have_bdy.F90,v 1.9 2008-12-09 00:31:57 kb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -29,7 +29,7 @@
! Original author(s): Karsten Bolding & Hans Burchard
!
! !LOCAL VARIABLES:
integer :: i,j,k,n
integer :: i,j,k,m,n
integer :: nbdy
integer :: f,l
!EOP
......@@ -43,8 +43,10 @@
nbdy = 0
i = 0
m = 0
if (NWB .ge. 1) then
do n = 1,NWB
m = m+1
if (wi(n) .ge. imin+ioff .and. wi(n) .le. imax+ioff) then
f = max(jmin+joff,wfj(n)) - joff
l = min(jmax+joff,wlj(n)) - joff
......@@ -54,6 +56,7 @@
wfj(i) = f
wlj(i) = l
nbdy = nbdy+1
bdy_2d_type(nbdy) = bdy_2d_type(m)
do k=1,nsbv
if (bdy_map(k,1) .eq. wi(i)+ioff .and. &
bdy_map(k,2) .eq. f+joff) then
......@@ -69,6 +72,7 @@
i = 0
if (NNB .ge. 1) then
do n = 1,NNB
m = m+1
if (nj(n) .ge. jmin+joff .and. nj(n) .le. jmax+joff) then
f = max(imin+ioff,nfi(n)) - ioff
l = min(imax+ioff,nli(n)) - ioff
......@@ -78,6 +82,7 @@
nli(i) = l
nj(i) = nj(n) - joff
nbdy = nbdy+1
bdy_2d_type(nbdy) = bdy_2d_type(m)
do k=1,nsbv
if (bdy_map(k,1) .eq. f+ioff .and. &
bdy_map(k,2) .eq. nj(i)+joff) then
......@@ -93,6 +98,7 @@
i = 0
if (NEB .ge. 1) then
do n = 1,NEB
m = m+1
if (ei(n) .ge. imin+ioff .and. ei(n) .le. imax+ioff) then
f = max(jmin+joff,efj(n)) - joff
l = min(jmax+joff,elj(n)) - joff
......@@ -102,6 +108,7 @@
efj(i) = f
elj(i) = l
nbdy = nbdy+1
bdy_2d_type(nbdy) = bdy_2d_type(m)
do k=1,nsbv
if (bdy_map(k,1) .eq. ei(i)+ioff .and. &
bdy_map(k,2) .eq. f+joff) then
......@@ -117,6 +124,7 @@
i = 0
if (NSB .ge. 1) then
do n = 1,NSB
m = m+1
if (sj(n) .ge. jmin+joff .and. sj(n) .le. jmax+joff) then
f = max(imin+ioff,sfi(n)) - ioff
l = min(imax+ioff,sli(n)) - ioff
......@@ -126,6 +134,7 @@
sli(i) = l
sj(i) = sj(n) - joff
nbdy = nbdy+1
bdy_2d_type(nbdy) = bdy_2d_type(m)
do k=1,nsbv
if (bdy_map(k,1) .eq. f+ioff .and. &
bdy_map(k,2) .eq. sj(i)+joff) then
......
!$Id: m2d.F90,v 1.24 2008-04-14 11:25:05 kb Exp $
!$Id: m2d.F90,v 1.25 2008-12-09 00:31:57 kb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -36,9 +36,9 @@
logical :: bdy2d=.false.
integer :: bdyfmt_2d,bdytype,bdyramp_2d=-1
character(len=PATH_MAX) :: bdyfile_2d
REAL_4B :: bdy_old(1500)
REAL_4B :: bdy_new(1500)
REAL_4B :: bdy_data(1500)
REAL_4B :: bdy_data_u(1500)
REAL_4B :: bdy_data_v(1500)
REAL_4B, allocatable :: bdy_times(:)
integer, parameter :: comm_method=-1
!
......
!$Id: update_2d_bdy.F90,v 1.8 2006-03-01 15:54:07 kbk Exp $
!$Id: update_2d_bdy.F90,v 1.9 2008-12-09 00:31:57 kb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -23,8 +23,15 @@
use domain, only: NWB,NNB,NEB,NSB,H,min_depth,imin,imax,jmin,jmax,az
use domain, only: wi,wfj,wlj,nj,nfi,nli,ei,efj,elj,sj,sfi,sli
use domain, only: bdy_index,nsbv
use m2d, only: dtm,bdyfmt_2d,bdy_data
use variables_2d, only: z
use domain, only: bdy_2d_type
use m2d, only: dtm,bdyfmt_2d,bdy_data,bdy_data_u,bdy_data_v
use variables_2d, only: z,D,U,DU,V,DV
#if defined(SPHERICAL) || defined(CURVILINEAR)
use domain, only: dxc,dyc
#else
use domain, only: dx,dy
#endif
IMPLICIT NONE
!
! !INPUT PARAMETERS:
......@@ -60,6 +67,7 @@
t = loop*dtm
select case (bdyfmt_2d)
case (NO_DATA)
case (ANALYTICAL)
#define OMEGA 2.*3.141592654/43200.
#ifdef COAST_TEST
......@@ -145,13 +153,33 @@
fac = _ONE_
if(bdyramp .gt. 1) fac=min( _ONE_ ,FOUR*loop/float(bdyramp))
! implicit Sommerfeldt
#if 0
eta = 1/(1+Cx)(eta_old + Cx*eta_new)
Cx=dtm/dx*srt(g*D)
z(i,j) = _ONE_/(_ONE_+Cx)*(z(i,j)+Cx*z(i,j-1))
#endif
l = 0
do n = 1,NWB
l = l+1
k = bdy_index(l)
i = wi(n)
do j = wfj(n),wlj(n)
z(i,j) = max(fac*bdy_data(k),-H(i,j)+min_depth)
select case (bdy_2d_type(l))
case (ZERO_GRADIENT)
z(i,j) = z(i+1,j)
case (SOMMERFELDT)
z(i,j) = z(i,j) + dtm/DXC*sqrt(9.81*D(i,j))*(z(i+1,j)-z(i,j))
case (CLAMPED)
z(i,j) = max(fac*bdy_data(k),-H(i,j)+min_depth)
case (FLATHER_ELEV)
a = sqrt(DU(i,j)/9.81)*(U(i,j)/DU(i,j)-bdy_data_u(k))
z(i,j) = max(fac*(bdy_data(k) - a),-H(i,j)+min_depth)
case default
FATAL 'Illegal NWB 2D boundary type selection'
stop 'update_2d_bdy()'
end select
k = k+1
end do
end do
......@@ -161,7 +189,20 @@
k = bdy_index(l)
j = nj(n)
do i = nfi(n),nli(n)
z(i,j) = max(fac*bdy_data(k),-H(i,j)+min_depth)
select case (bdy_2d_type(l))
case (ZERO_GRADIENT)
z(i,j) = z(i,j-1)
case (SOMMERFELDT)
z(i,j) = z(i,j) - dtm/DXC*sqrt(9.81*D(i,j))*(z(i,j)-z(i,j-1))
case (CLAMPED)
z(i,j) = max(fac*bdy_data(k),-H(i,j)+min_depth)
case (FLATHER_ELEV)
a = sqrt(DV(i,j)/9.81)*(V(i,j-1)/DV(i,j-1)+bdy_data_v(k))
z(i,j) = max(fac*(bdy_data(k) + a),-H(i,j)+min_depth)
case default
FATAL 'Illegal NNB 2D boundary type selection'
stop 'update_2d_bdy()'
end select
k = k+1
end do
end do
......@@ -171,7 +212,20 @@
k = bdy_index(l)
i = ei(n)
do j = efj(n),elj(n)
z(i,j) = max(fac*bdy_data(k),-H(i,j)+min_depth)
select case (bdy_2d_type(l))
case (ZERO_GRADIENT)
z(i,j) = z(i-1,j)
case (SOMMERFELDT)
z(i,j) = z(i,j) - dtm/DXC*sqrt(9.81*D(i,j))*(z(i,j)-z(i-1,j))
case (CLAMPED)
z(i,j) = max(fac*bdy_data(k),-H(i,j)+min_depth)
case (FLATHER_ELEV)
a = sqrt(DU(i,j)/9.81)*(U(i-1,j)/DU(i-1,j)-bdy_data_u(k))
z(i,j) = max(fac*(bdy_data(k) + a),-H(i,j)+min_depth)
case default
FATAL 'Illegal NEB 2D boundary type selection'
stop 'update_2d_bdy()'
end select
k = k+1
end do
end do
......@@ -181,7 +235,20 @@
k = bdy_index(l)
j = sj(n)
do i = sfi(n),sli(n)
z(i,j) = max(fac*bdy_data(k),-H(i,j)+min_depth)
select case (bdy_2d_type(l))
case (ZERO_GRADIENT)
z(i,j) = z(i,j+1)
case (SOMMERFELDT)
z(i,j) = z(i,j) + dtm/DXC*sqrt(9.81*D(i,j))*(z(i,j+1)-z(i,j))
case (CLAMPED)
z(i,j) = max(fac*bdy_data(k),-H(i,j)+min_depth)
case (FLATHER_ELEV)
a = sqrt(DV(i,j)/9.81)*(V(i,j)/DV(i,j)-bdy_data_v(k))
z(i,j) = max(fac*(bdy_data(k) - a),-H(i,j)+min_depth)
case default
FATAL 'Illegal NSB 2D boundary type selection'
stop 'update_2d_bdy()'
end select
k = k+1
end do
end do
......
!$Id: bdy_spec.F90,v 1.5 2005-04-29 12:55:31 kbk Exp $
!$Id: bdy_spec.F90,v 1.6 2008-12-09 00:31:57 kb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -14,6 +14,8 @@
use domain, only: NWB,NNB,NEB,NSB,NOB
use domain, only: wi,wfj,wlj,nj,nfi,nli,ei,efj,elj,sj,sfi,sli
use domain, only: bdy_index,bdy_map,nsbv
use domain, only: bdy_2d_type,bdy_3d_type
use domain, only: need_2d_bdy_elev,need_2d_bdy_u,need_2d_bdy_v
IMPLICIT NONE
!
! !INPUT PARAMETERS:
......@@ -27,6 +29,9 @@
! Original author(s): Karsten Bolding & Hans Burchard
!
! $Log: bdy_spec.F90,v $
! Revision 1.6 2008-12-09 00:31:57 kb
! added new 2D open boundaries
!
! Revision 1.5 2005-04-29 12:55:31 kbk
! removing print statement
!
......@@ -50,6 +55,7 @@
!
! !LOCAL VARIABLES:
integer :: i,j,k,l,n,rc
integer :: type_2d(4,10),type_3d(4,10)
!
!EOP
!-----------------------------------------------------------------------
......@@ -75,9 +81,17 @@
allocate(wlj(NWB),stat=rc)
if (rc /= 0) stop 'bdy_spec: Error allocating memory (wlj)'
do n = 1,NWB
read(BDYINFO,*,END=91,ERR=92) wi(n),wfj(n),wlj(n)
read(BDYINFO,*,END=91,ERR=92) wi(n),wfj(n),wlj(n), &
type_2d(1,n),type_3d(1,n)
nsbv = nsbv + (wlj(n)-wfj(n)+1)
end do
do n = 1,NWB
if (type_2d(1,n) .eq. CLAMPED) need_2d_bdy_elev = .true.
if (type_2d(1,n) .eq. FLATHER_ELEV) then
need_2d_bdy_elev = .true.
need_2d_bdy_u = .true.
end if
end do
end if
! Northern boundary info
......@@ -90,9 +104,17 @@
allocate(nli(NNB),stat=rc)
if (rc /= 0) stop 'bdy_spec: Error allocating memory (nli)'
do n = 1,NNB
read(BDYINFO,*,END=91,ERR=92) nj(n),nfi(n),nli(n)
read(BDYINFO,*,END=91,ERR=92) nj(n),nfi(n),nli(n), &
type_2d(2,n),type_3d(2,n)
nsbv = nsbv + (nli(n)-nfi(n)+1)
end do
do n = 1,NNB
if (type_2d(2,n) .eq. CLAMPED) need_2d_bdy_elev = .true.
if (type_2d(2,n) .eq. FLATHER_ELEV) then
need_2d_bdy_elev = .true.
need_2d_bdy_v = .true.
end if
end do
end if
! Eastern boundary info
......@@ -105,9 +127,17 @@
allocate(elj(NEB),stat=rc)
if (rc /= 0) stop 'bdy_spec: Error allocating memory (elj)'
do n = 1,NEB
read(BDYINFO,*,END=91,ERR=92) ei(n),efj(n),elj(n)
read(BDYINFO,*,END=91,ERR=92) ei(n),efj(n),elj(n), &
type_2d(3,n),type_3d(3,n)
nsbv = nsbv + (elj(n)-efj(n)+1)
end do
do n = 1,NEB
if (type_2d(3,n) .eq. CLAMPED) need_2d_bdy_elev = .true.
if (type_2d(3,n) .eq. FLATHER_ELEV) then
need_2d_bdy_elev = .true.
need_2d_bdy_u = .true.
end if
end do
end if
! Southern boundary info
......@@ -120,15 +150,51 @@
allocate(sli(NSB),stat=rc)
if (rc /= 0) stop 'bdy_spec: Error allocating memory (sli)'
do n = 1,NSB
read(BDYINFO,*,END=91,ERR=92) sj(n),sfi(n),sli(n)
read(BDYINFO,*,END=91,ERR=92) sj(n),sfi(n),sli(n), &
type_2d(4,n),type_3d(4,n)
nsbv = nsbv + (sli(n)-sfi(n)+1)
end do
do n=1,NSB
if (type_2d(4,n) .eq. CLAMPED) need_2d_bdy_elev = .true.
if (type_2d(4,n) .eq. FLATHER_ELEV) then
need_2d_bdy_elev = .true.
need_2d_bdy_v = .true.
end if
end do
end if
close(BDYINFO)
NOB = NWB+NNB+NEB+NSB
allocate(bdy_2d_type(NOB),stat=rc)
if (rc /= 0) stop 'bdy_spec: Error allocating memory (bdy_2d_type)'
allocate(bdy_3d_type(NOB),stat=rc)
if (rc /= 0) stop 'bdy_spec: Error allocating memory (bdy_3d_type)'
l=1
do n = 1,NWB
bdy_2d_type(l) = type_2d(1,n)
bdy_3d_type(l) = type_3d(1,n)
l=l+1
end do
do n = 1,NNB
bdy_2d_type(l) = type_2d(2,n)
bdy_3d_type(l) = type_3d(2,n)
l=l+1
end do
do n = 1,NEB
bdy_2d_type(l) = type_2d(3,n)
bdy_3d_type(l) = type_3d(3,n)
l=l+1
end do
do n = 1,NSB
bdy_2d_type(l) = type_2d(4,n)
bdy_3d_type(l) = type_3d(4,n)
l=l+1
end do
LEVEL2 '# of open boundaries ',NOB
LEVEL2 '# of open bdy points ',nsbv
......
!$Id: domain.F90,v 1.28 2007-10-16 06:22:56 kbk Exp $
!$Id: domain.F90,v 1.29 2008-12-09 00:31:57 kb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -62,12 +62,19 @@
integer :: nsbv
integer :: ioff=0,joff=0
integer, dimension(:), allocatable :: bdy_2d_type
integer, dimension(:), allocatable :: bdy_3d_type
integer, dimension(:), allocatable :: wi,wfj,wlj
integer, dimension(:), allocatable :: nj,nfi,nli
integer, dimension(:), allocatable :: ei,efj,elj
integer, dimension(:), allocatable :: sj,sfi,sli
integer, allocatable :: bdy_index(:),bdy_map(:,:)
character(len=64) :: bdy_2d_desc(5)
logical :: need_2d_bdy_elev = .false.
logical :: need_2d_bdy_u = .false.
logical :: need_2d_bdy_v = .false.
REALTYPE :: cori= _ZERO_
! method for specifying bottom roughness (0=const, 1=from topo.nc)
......@@ -84,6 +91,9 @@
! Original author(s): Karsten Bolding & Hans Burchard
!
! $Log: domain.F90,v $
! Revision 1.29 2008-12-09 00:31:57 kb
! added new 2D open boundaries
!
! Revision 1.28 2007-10-16 06:22:56 kbk
! curvi-linear now runs in parallel
!
......@@ -294,6 +304,12 @@
write(debug,*) 'init_domain()'
#endif
bdy_2d_desc(ZERO_GRADIENT) = "Zero gradient"
bdy_2d_desc(SOMMERFELDT) = "Sommerfeldt rad."
bdy_2d_desc(CLAMPED) = "Clamped"
bdy_2d_desc(FLATHER_ELEV) = "Flather (elev)"
bdy_2d_desc(FLATHER_VEL) = "Flather (vel)"
LEVEL1 'init_domain'
! Read domain specific things from the namelist.
......
!$Id: print_bdy.F90,v 1.2 2003-04-23 11:59:39 kbk Exp $
!$Id: print_bdy.F90,v 1.3 2008-12-09 00:31:57 kb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -13,6 +13,8 @@
! !USES:
use domain, only: NWB,NNB,NEB,NSB
use domain, only: wi,wfj,wlj,nj,nfi,nli,ei,efj,elj,sj,sfi,sli
use domain, only: bdy_2d_type,bdy_3d_type
use domain, only: bdy_2d_desc
IMPLICIT NONE
!
! !INPUT PARAMETERS:
......@@ -26,6 +28,9 @@
! Original author(s): Karsten Bolding & Hans Burchard
!
! $Log: print_bdy.F90,v $
! Revision 1.3 2008-12-09 00:31:57 kb
! added new 2D open boundaries
!
! Revision 1.2 2003-04-23 11:59:39 kbk
! update_2d_halo on spherical variables + TABS to spaces
!
......@@ -37,7 +42,7 @@
!
!
! !LOCAL VARIABLES:
integer :: n
integer :: m,n
!
!EOP
!-----------------------------------------------------------------------
......@@ -48,30 +53,35 @@
write(debug,*) 'print_bdy() # ',Ncall
#endif
m=0
LEVEL2 TRIM(header)
LEVEL2 'There are ',NWB+NNB+NEB+NSB,' open boundaries.'
if (NWB .ge. 1) then
LEVEL3 'Western Boundary'
do n = 1,NWB
LEVEL3 wi(n),wfj(n),wlj(n)
m=m+1
LEVEL3 wi(n),wfj(n),wlj(n),trim(bdy_2d_desc(bdy_2d_type(m)))
end do
end if
if (NNB .ge. 1) then
LEVEL3 'Northern Boundary'
do n = 1,NNB
LEVEL3 nj(n),nfi(n),nli(n)
m=m+1
LEVEL3 nj(n),nfi(n),nli(n),trim(bdy_2d_desc(bdy_2d_type(m)))
end do
end if
if (NEB .ge. 1) then
LEVEL3 'Eastern Boundary'
do n = 1,NEB
LEVEL3 ei(n),efj(n),elj(n)
m=m+1
LEVEL3 ei(n),efj(n),elj(n),trim(bdy_2d_desc(bdy_2d_type(m)))
end do
end if
if (NSB .ge. 1) then
LEVEL3 'Southern Boundary'
do n = 1,NSB
LEVEL3 sj(n),sfi(n),sli(n)
m=m+1
LEVEL3 sj(n),sfi(n),sli(n),trim(bdy_2d_desc(bdy_2d_type(m)))
end do
end if
......
!$Id: get_2d_bdy.F90,v 1.2 2003-04-23 12:04:08 kbk Exp $
!$Id: get_2d_bdy.F90,v 1.3 2008-12-09 00:31:58 kb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -26,6 +26,9 @@
! Original author(s): Karsten Bolding & Hans Burchard
!
! $Log: get_2d_bdy.F90,v $
! Revision 1.3 2008-12-09 00:31:58 kb
! added new 2D open boundaries
!
! Revision 1.2 2003-04-23 12:04:08 kbk
! cleaned code + TABS to spaces
!
......@@ -52,6 +55,7 @@
write(debug,*) 'get_2d_bdy() # ',ncall
#endif
select case (fmt)
case (NO_DATA)
case (ANALYTICAL)
case (ASCII)
STDERR 'should get ASCII boundary data'
......
!$Id: init_2d_bdy.F90,v 1.2 2003-04-23 12:04:08 kbk Exp $
!$Id: init_2d_bdy.F90,v 1.3 2008-12-09 00:31:58 kb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -27,6 +27,9 @@
! Original author(s): Karsten Bolding & Hans Burchard
!
! $Log: init_2d_bdy.F90,v $
! Revision 1.3 2008-12-09 00:31:58 kb
! added new 2D open boundaries
!
! Revision 1.2 2003-04-23 12:04:08 kbk
! cleaned code + TABS to spaces
!
......@@ -61,6 +64,7 @@
LEVEL2 'init_2d_bdy'
select case (fmt)
case (NO_DATA)
case (ANALYTICAL)
LEVEL3 'Analytical boundary formulations'
stop 'init_2d_bdy'
......
!$Id: ncdf_2d_bdy.F90,v 1.6 2007-09-30 13:00:43 kbk Exp $
!$Id: ncdf_2d_bdy.F90,v 1.7 2008-12-09 00:31:58 kb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -11,9 +11,11 @@
! !DESCRIPTION:
!
! !USES:
use m2d, only: dtm,bdy_times,bdy_old,bdy_new,bdy_data
!KB use m2d, only: dtm,bdy_times,bdy_old,bdy_new,bdy_data
use m2d, only: dtm,bdy_times,bdy_data,bdy_data_u,bdy_data_v
use time, only: string_to_julsecs,time_diff,julianday,secondsofday
use time, only: write_time_string,timestr
use domain, only: need_2d_bdy_elev,need_2d_bdy_u,need_2d_bdy_v
IMPLICIT NONE
!
private
......@@ -22,14 +24,26 @@
!
! !PRIVATE DATA MEMBERS:
integer :: ncid
integer :: time_id,elev_id,nsets,bdy_len
integer :: time_id,elev_id=-1,nsets,bdy_len
integer :: u_id=-1, v_id=-1
integer :: start(2),edges(2)
REALTYPE :: offset
REAL_4B :: bdy_old(1500)
REAL_4B :: bdy_new(1500)
REAL_4B :: bdy_old_u(1500)
REAL_4B :: bdy_new_u(1500)
REAL_4B :: bdy_old_v(1500)
REAL_4B :: bdy_new_v(1500)
!
!
! !REVISION HISTORY:
! Original author(s): Karsten Bolding & Hans Burchard