Commit 9d9484e2 authored by Knut's avatar Knut
Browse files

Dvel and Dveln available

parent e0d2f623
......@@ -5,7 +5,7 @@
! !ROUTINE: depth_update - adjust the depth to new elevations.
!
! !INTERFACE:
subroutine depth_update(zo,z,Dlast,D,DU,DV,first,from3d)
subroutine depth_update(zo,z,D,Dvel,DU,DV,from3d)
! Note (KK): keep in sync with interface in m2d.F90
!
......@@ -29,24 +29,20 @@
IMPLICIT NONE
!
! !INPUT PARAMETERS:
REALTYPE,dimension(E2DFIELD),intent(in) :: zo,z
logical,intent(in),optional :: first,from3d
!
! !INPUT/OUTPUT PARAMETERS:
REALTYPE,dimension(:,:),pointer,intent(inout) :: Dlast,D
REALTYPE,dimension(E2DFIELD),intent(in) :: zo,z
logical,intent(in),optional :: from3d
!
! !OUTPUT PARAMETERS:
REALTYPE,dimension(E2DFIELD),intent(out) :: DU,DV
REALTYPE,dimension(E2DFIELD),intent(out) :: D,Dvel,DU,DV
!
! !REVISION HISTORY:
! Original author(s): Hans Burchard & Karsten Bolding
!
! !LOCAL VARIABLES:
integer :: i,j
REALTYPE :: d1,d2i,x
REALTYPE,dimension(E2DFIELD) :: ztmp
REALTYPE,dimension(:,:),pointer :: p2d
logical :: update_dryingfactors
REALTYPE :: d1,d2i
REALTYPE,dimension(E2DFIELD) :: zvel
logical :: update_dryingfactors
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -60,31 +56,17 @@
! TODO/BJB: Why is this turned off?
#undef USE_MASK
p2d => Dlast ; Dlast => D ; D => p2d
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,d1,d2i,x)
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,d1,d2i)
! Depth in elevation points
if (present(first)) then
if (first) then
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
! TODO/BJB: Is it enough to do this on az?
Dlast(i,j) = zo(i,j) + H(i,j)
end do
end do
!$OMP END DO NOWAIT
end if
end if
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
! TODO/BJB: Is it enough to do this on az?
D(i,j) = z(i,j)+H(i,j)
ztmp(i,j) = _HALF_ * ( zo(i,j) + z(i,j) )
zvel(i,j) = _HALF_ * ( zo(i,j) + z(i,j) )
Dvel(i,j) = zvel(i,j) + H(i,j)
end do
end do
!$OMP END DO
......@@ -96,8 +78,8 @@
#ifdef USE_MASK
if(au(i,j) .gt. 0) then
#endif
x=max(_HALF_*(ztmp(i,j)+ztmp(i+1,j)),-HU(i,j)+min_depth)
DU(i,j) = x+HU(i,j)
DU(i,j) = max( min_depth , &
_HALF_*(zvel(i,j)+zvel(i+1,j)) + HU(i,j) )
#ifdef USE_MASK
end if
#endif
......@@ -112,8 +94,8 @@
#ifdef USE_MASK
if(av(i,j) .gt. 0) then
#endif
x = max(_HALF_*(ztmp(i,j)+ztmp(i,j+1)),-HV(i,j)+min_depth)
DV(i,j) = x+HV(i,j)
DV(i,j) = max( min_depth , &
_HALF_*(zvel(i,j)+zvel(i,j+1)) + HV(i,j) )
#ifdef USE_MASK
end if
#endif
......
......@@ -12,11 +12,11 @@
allocate(t_z(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (t_z)'
allocate(t_Dlast(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (t_Dlast)'
allocate(D(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (D)'
allocate(t_D(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (t_D)'
allocate(Dvel(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (Dvel)'
allocate(DU(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (DU)'
......
......@@ -2,7 +2,7 @@
integer,dimension(:,:),allocatable :: break_mask
integer,dimension(:,:),allocatable :: break_stat
#endif
REALTYPE,dimension(:,:),allocatable,target :: t_zo,t_z,t_Dlast,t_D,DU,DV
REALTYPE,dimension(:,:),allocatable,target :: t_zo,t_z,D,Dvel,DU,DV
REALTYPE,dimension(:,:),allocatable :: U,V
REALTYPE,dimension(:,:),allocatable :: UEx,VEx
REALTYPE,dimension(:,:),allocatable :: fU,fV
......
......@@ -34,20 +34,19 @@
interface
subroutine depth_update(zo,z,Dlast,D,DU,DV,first,from3d)
subroutine depth_update(zo,z,D,Dvel,DU,DV,from3d)
use domain, only: imin,imax,jmin,jmax
IMPLICIT NONE
REALTYPE,dimension(E2DFIELD),intent(in) :: zo,z
logical,intent(in),optional :: first,from3d
REALTYPE,dimension(:,:),pointer,intent(inout) :: Dlast,D
REALTYPE,dimension(E2DFIELD),intent(out) :: DU,DV
REALTYPE,dimension(E2DFIELD),intent(in) :: zo,z
logical,intent(in),optional :: from3d
REALTYPE,dimension(E2DFIELD),intent(out) :: D,Dvel,DU,DV
end subroutine depth_update
subroutine uv_advect(U,V,Dold,Dnew,DU,DV)
subroutine uv_advect(U,V,D,Dvel,DU,DV)
use domain, only: imin,imax,jmin,jmax
IMPLICIT NONE
REALTYPE,dimension(E2DFIELD),intent(in) :: U,V
REALTYPE,dimension(E2DFIELD),target,intent(in) :: Dold,Dnew,DU,DV
REALTYPE,dimension(E2DFIELD),target,intent(in) :: D,Dvel,DU,DV
end subroutine uv_advect
subroutine uv_diffusion(An_method,U,V,D,DU,DV)
......@@ -205,7 +204,7 @@
zo = z
! KK-TODO: check whether we need D[ |U|V] in init_3d
! otherwise we can move this call by default to postinit_2d
call depth_update(zo,z,Dlast,D,DU,DV,first=.true.)
call depth_update(zo,z,D,Dvel,DU,DV)
end if
if (Am .lt. _ZERO_) then
......@@ -405,7 +404,7 @@
end where
end if
call depth_update(zo,z,Dlast,D,DU,DV,first=.true.)
call depth_update(zo,z,D,Dvel,DU,DV)
end if
......@@ -471,7 +470,7 @@
end if
call tic(TIM_INTEGR2D)
call uv_advect(U,V,Dlast,D,DU,DV)
call uv_advect(U,V,D,Dvel,DU,DV)
call uv_diffusion(An_method,U,V,D,DU,DV) ! Has to be called after uv_advect.
call toc(TIM_INTEGR2D)
......@@ -484,7 +483,7 @@
end if
if (have_boundaries) call update_2d_bdy(loop,bdy2d_ramp)
call sealevel()
call depth_update(zo,z,Dlast,D,DU,DV)
call depth_update(zo,z,D,Dvel,DU,DV)
if(residual .gt. 0) then
call tic(TIM_INTEGR2D)
......
......@@ -5,7 +5,7 @@
integer break_mask(E2DFIELD)
integer break_stat(E2DFIELD)
#endif
REALTYPE,dimension(E2DFIELD),target :: t_zo,t_z,t_Dlast,t_D,DU,DV
REALTYPE,dimension(E2DFIELD),target :: t_zo,t_z,D,Dvel,DU,DV
REALTYPE U(E2DFIELD)
REALTYPE V(E2DFIELD)
REALTYPE UEx(E2DFIELD)
......
......@@ -5,7 +5,7 @@
! !ROUTINE: uv_advect - 2D advection of momentum \label{sec-uv-advect}
!
! !INTERFACE:
subroutine uv_advect(U,V,Dold,Dnew,DU,DV)
subroutine uv_advect(U,V,D,Dvel,DU,DV)
! Note (KK): keep in sync with interface in m2d.F90
!
......@@ -32,16 +32,16 @@
!
! !INPUT PARAMETERS:
REALTYPE,dimension(E2DFIELD),intent(in) :: U,V
REALTYPE,dimension(E2DFIELD),target,intent(in) :: Dold,Dnew,DU,DV
REALTYPE,dimension(E2DFIELD),target,intent(in) :: D,Dvel,DU,DV
!
! !REVISION HISTORY:
! Original author(s): Hans Burchard & Karsten Bolding
!
! !LOCAL VARIABLES:
integer :: i,j
REALTYPE,dimension(E2DFIELD) :: fadv,Uadv,Vadv,DUadv,DVadv
REALTYPE,dimension(E2DFIELD),target :: Dadv
REALTYPE,dimension(:,:),pointer :: pDadv
REALTYPE,dimension(E2DFIELD) :: fadv,Uadv,Vadv
REALTYPE,dimension(E2DFIELD),target :: Dadv,DUadv,DVadv
REALTYPE,dimension(:,:),pointer :: pDadv,pDUadv,pDVadv
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -64,6 +64,19 @@
! Here begins dimensional split advection for u-velocity
!$OMP SINGLE
! KK-TODO: _POINTER_REMAP_, but this requires that D[U|V] become
! pointers like mask_[u|v]flux in do_advection and similar
! for h[u|v] in do_advection_3d and in uv_advect_3d
!#ifdef _POINTER_REMAP_
! pDUadv(imin-HALO:,jmin-HALO:) => Dvel(1+_IRANGE_HALO_,_JRANGE_HALO_)
!#else
pDUadv => DUadv
pDUadv(_IRANGE_HALO_-1,:) = Dvel(1+_IRANGE_HALO_,:)
!#endif
pDVadv => DVadv
!$OMP END SINGLE
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin-HALO,jmax+HALO
......@@ -77,10 +90,9 @@
Uadv(i,j) = _HALF_*( U(i,j) + U(i+1,j) )
Vadv(i,j) = _HALF_*( V(i,j) + V(i+1,j) )
end if
DUadv(i,j) = _HALF_*( Dold(i+1,j) + Dnew(i+1,j) )
! Note (KK): DV only valid until jmax+1
! therefore DVadv only valid until jmax+1
DVadv(i,j) = _HALF_*( DV(i,j) + DV(i+1,j) )
! therefore pDVadv only valid until jmax+1
pDVadv(i,j) = _HALF_*( DV(i,j) + DV(i+1,j) )
end do
#ifndef SLICE_MODEL
end do
......@@ -178,14 +190,27 @@
call toc(TIM_UVADVH)
end if
call do_advection(dtm,fadv,Uadv,Vadv,DUadv,DVadv,pDadv,pDadv, &
vel2d_adv_split,vel2d_adv_hor,_ZERO_,U_TAG, &
call do_advection(dtm,fadv,Uadv,Vadv,pDUadv,pDVadv,pDadv,pDadv, &
vel2d_adv_split,vel2d_adv_hor,_ZERO_,U_TAG, &
advres=UEx)
!$OMP END SINGLE
! Here begins dimensional split advection for v-velocity
!$OMP SINGLE
! KK-TODO: _POINTER_REMAP_, but this requires that D[U|V] become
! pointers like mask_[u|v]flux in do_advection and similar
! for h[u|v] in do_advection_3d and in uv_advect_3d
pDUadv => DUadv
!#ifdef _POINTER_REMAP_
! pDVadv(imin-HALO:,jmin-HALO:) => Dvel(_IRANGE_HALO_,1+_JRANGE_HALO_)
!#else
pDVadv => DVadv
pDVadv(:,_JRANGE_HALO_-1) = Dvel(:,1+_JRANGE_HALO_)
!#endif
!$OMP END SINGLE
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin-HALO,jmax+HALO-1
......@@ -200,9 +225,8 @@
Vadv(i,j) = _HALF_*( V(i,j) + V(i,j+1) )
end if
! Note (KK): DU only valid until imax+1
! therefore DUadv only valid until imax+1
DUadv(i,j) = _HALF_*( DU(i,j) + DU(i,j+1) )
DVadv(i,j) = _HALF_*( Dold(i,j+1) + Dnew(i,j+1) )
! therefore pDUadv only valid until imax+1
pDUadv(i,j) = _HALF_*( DU(i,j) + DU(i,j+1) )
end do
#ifndef SLICE_MODEL
end do
......@@ -288,8 +312,8 @@
call toc(TIM_UVADVH)
end if
call do_advection(dtm,fadv,Uadv,Vadv,DUadv,DVadv,pDadv,pDadv, &
vel2d_adv_split,vel2d_adv_hor,_ZERO_,V_TAG, &
call do_advection(dtm,fadv,Uadv,Vadv,pDUadv,pDVadv,pDadv,pDadv, &
vel2d_adv_split,vel2d_adv_hor,_ZERO_,V_TAG, &
advres=VEx)
else ! if (vel2d_adv_hor .eq. NOADV)
......
......@@ -21,7 +21,7 @@
IMPLICIT NONE
!
! !PUBLIC DATA MEMBERS:
REALTYPE,dimension(:,:),pointer :: zo,z,Dlast,D
REALTYPE,dimension(:,:),pointer :: zo,z
#ifdef STATIC
#include "static_2d.h"
......@@ -112,7 +112,6 @@
#endif
zo => t_zo ; z => t_z
Dlast => t_Dlast ; D => t_D
#ifdef USE_BREAKS
break_stat = 0
......@@ -121,7 +120,7 @@
z = _ZERO_; zo =_ZERO_
zub=_ZERO_ ; zub0=_ZERO_
zvb=_ZERO_ ; zvb0=_ZERO_
Dlast = _ZERO_ ; D = _ZERO_
D = _ZERO_ ; Dvel = _ZERO_
U = _ZERO_; DU = _ZERO_; fU = _ZERO_; Uint = _ZERO_; UEx = _ZERO_
V = _ZERO_; DV = _ZERO_; fV = _ZERO_; Vint = _ZERO_; VEx = _ZERO_
......
......@@ -165,11 +165,11 @@
allocate(ssvn(I2DFIELD),stat=rc) ! Elevation after macro time step (v-column)
if (rc /= 0) stop 'init_3d: Error allocating memory (ssvn)'
allocate(t_Dold(I2DFIELD),stat=rc) ! depth before macro time step (z-column)
if (rc /= 0) stop 'init_3d: Error allocating memory (t_Dold)'
allocate(Dn(I2DFIELD),stat=rc) ! depth after macro time step (z-column)
if (rc /= 0) stop 'init_3d: Error allocating memory (Dn)'
allocate(t_Dn(I2DFIELD),stat=rc) ! depth after macro time step (z-column)
if (rc /= 0) stop 'init_3d: Error allocating memory (t_Dn)'
allocate(Dveln(I2DFIELD),stat=rc) ! depth during macro time step (z-column)
if (rc /= 0) stop 'init_3d: Error allocating memory (Dveln)'
allocate(Dun(I2DFIELD),stat=rc) ! depth after macro time step (u-column)
if (rc /= 0) stop 'init_3d: Error allocating memory (Dun)'
......
......@@ -51,7 +51,7 @@
REALTYPE, dimension(:,:), allocatable :: sseo,ssen
REALTYPE, dimension(:,:), allocatable :: ssuo,ssun
REALTYPE, dimension(:,:), allocatable :: ssvo,ssvn
REALTYPE,dimension(:,:),allocatable,target :: t_Dold,t_Dn,Dun,Dvn
REALTYPE,dimension(:,:),allocatable,target :: Dn,Dveln,Dun,Dvn
! 3D friction in 3D domain
REALTYPE, dimension(:,:), allocatable :: rru,rrv
......
......@@ -236,7 +236,6 @@
if (.not. hotstart) then
ssen = z
call start_macro()
Dold = Dn
call coordinates(hotstart)
call hcc_check()
end if
......@@ -387,7 +386,7 @@
end do
end do
call depth_update(sseo,ssen,Dold,Dn,Dun,Dvn,first=.true.,from3d=.true.)
call depth_update(sseo,ssen,Dn,Dveln,Dun,Dvn,from3d=.true.)
call coordinates(hotstart)
end if
......@@ -604,7 +603,7 @@
#endif
call tic(TIM_INTEGR3D)
call uv_advect(Uint,Vint,Dold,Dn,Dun,Dvn)
call uv_advect(Uint,Vint,Dn,Dveln,Dun,Dvn)
call uv_diffusion(0,Uint,Vint,Dn,Dun,Dvn) ! Has to be called after uv_advect.
call toc(TIM_INTEGR3D)
......
......@@ -29,10 +29,9 @@
!
! !USES:
use domain, only: imin,imax,jmin,jmax,H,HU,HV,min_depth
use m2d, only: depth_update
use m2d, only: z,Uint,Vint
use m3d, only: M
use variables_3d, only: sseo,ssen,ssuo,ssun,ssvo,ssvn,Dold,Dn,Dun,Dvn
use variables_3d, only: sseo,ssen,ssuo,ssun,ssvo,ssvn,Dn,Dveln,Dun,Dvn
use getm_timers, only: tic, toc, TIM_STARTMCR
IMPLICIT NONE
!
......@@ -41,8 +40,8 @@
!
! !LOCAL VARIABLES:
integer :: i,j
REALTYPE,dimension(I2DFIELD) :: ssevel
REALTYPE :: split
REALTYPE,dimension(:,:),pointer :: p2d
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -53,32 +52,35 @@
#endif
call tic(TIM_STARTMCR)
p2d => Dold ; Dold => Dn ; Dn => p2d
do j=jmin-HALO,jmax+HALO ! Defining 'old' and 'new' sea surface
do i=imin-HALO,imax+HALO ! elevation for macro time step
sseo(i,j)=ssen(i,j)
ssen(i,j)=z(i,j)
ssevel(i,j) = _HALF_ * ( sseo(i,j) + ssen(i,j) )
Dn(i,j) = ssen(i,j) + H(i,j)
Dveln(i,j) = ssevel(i,j) + H(i,j)
! KK-TODO: use of Dn & Co. in more routines (coordinates,momentum,rivers)
! and replacement of ssun+HU by Dun!
! and removement of ssun?
! calculation of Dveln,Dun,Dvn by depth_update
end do
end do
do j=jmin-HALO,jmax+HALO ! Same for U-points
do i=imin-HALO,imax+HALO-1
ssuo(i,j)=ssun(i,j)
ssun(i,j)=0.25*(sseo(i,j)+sseo(i+1,j)+ssen(i,j)+ssen(i+1,j))
ssun(i,j)=max(ssun(i,j),-HU(i,j)+min_depth)
Dun(i,j) = ssun(i,j) + HU(i,j)
Dun(i,j) = max( min_depth , &
_HALF_*(ssevel(i,j)+ssevel(i+1,j)) + HU(i,j) )
ssun(i,j) = Dun(i,j) - HU(i,j)
end do
end do
do j=jmin-HALO,jmax+HALO-1
do i=imin-HALO,imax+HALO ! Same for V-points
ssvo(i,j)=ssvn(i,j)
ssvn(i,j)=0.25*(sseo(i,j)+sseo(i,j+1)+ssen(i,j)+ssen(i,j+1))
ssvn(i,j)=max(ssvn(i,j),-HV(i,j)+min_depth)
Dvn(i,j) = ssvn(i,j) + HV(i,j)
Dvn(i,j) = max( min_depth , &
_HALF_*(ssevel(i,j)+ssevel(i,j+1)) + HV(i,j) )
ssvn(i,j) = Dvn(i,j) - HV(i,j)
end do
end do
......
......@@ -81,7 +81,7 @@
REALTYPE :: ssun(I2DFIELD)
REALTYPE :: ssvo(I2DFIELD)
REALTYPE :: ssvn(I2DFIELD)
REALTYPE,dimension(I2DFIELD),target :: t_Dold,t_Dn,Dun,Dvn
REALTYPE,dimension(I2DFIELD),target :: Dn,Dveln,Dun,Dvn
! 3D friction in 3D domain
REALTYPE :: rru(I2DFIELD)
......
......@@ -115,7 +115,6 @@
IMPLICIT NONE
!
! !PUBLIC DATA MEMBERS:
REALTYPE,dimension(:,:),pointer :: Dold,Dn
REALTYPE :: dt,cnpar=0.9
REALTYPE :: avmback=_ZERO_,avhback=_ZERO_
logical :: do_numerical_analyses=.false.
......@@ -192,8 +191,6 @@
#include "dynamic_allocations_3d.h"
#endif
Dold => t_Dold ; Dn => t_Dn
hn = _ZERO_ ; hun = _ZERO_ ; hvn = _ZERO_
uu = _ZERO_ ; vv = _ZERO_ ; ww = _ZERO_
#ifdef _MOMENTUM_TERMS_
......@@ -203,7 +200,7 @@
cor_v = _ZERO_ ; epg_v = _ZERO_ ; ipg_v = _ZERO_
#endif
ssen = _ZERO_ ; ssun = _ZERO_ ; ssvn = _ZERO_
Dold = _ZERO_ ; Dn = _ZERO_ ; Dun = _ZERO_ ; Dvn = _ZERO_
Dn = _ZERO_ ; Dveln = _ZERO_ ; Dun = _ZERO_ ; Dvn = _ZERO_
rru= _ZERO_ ; rrv= _ZERO_
uuEx= _ZERO_ ; vvEx= _ZERO_
tke=1.e-10 ; eps=1.e-10
......
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