Commit 69425337 authored by Knut's avatar Knut
Browse files

central depth_update routine + pointer swapping

parent f959deb6
......@@ -5,7 +5,9 @@
! !ROUTINE: depth_update - adjust the depth to new elevations.
!
! !INTERFACE:
subroutine depth_update
subroutine depth_update(zo,z,Dlast,D,DU,DV,first,from3d)
! Note (KK): keep in sync with interface in m2d.F90
!
! !DESCRIPTION:
!
......@@ -22,17 +24,29 @@
! !USES:
use domain, only: imin,imax,jmin,jmax,H,HU,HV,min_depth,crit_depth
use domain, only: az,au,av,dry_z,dry_u,dry_v
use variables_2d, only: D,z,zo,DU,DV
use getm_timers, only: tic, toc, TIM_DPTHUPDATE
!$ use omp_lib
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
!
! !OUTPUT PARAMETERS:
REALTYPE,dimension(E2DFIELD),intent(out) :: 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
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -46,18 +60,34 @@
! 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)
! 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) )
end do
end do
!$OMP END DO NOWAIT
!$OMP END DO
! U-points
!$OMP DO SCHEDULE(RUNTIME)
......@@ -66,7 +96,7 @@
#ifdef USE_MASK
if(au(i,j) .gt. 0) then
#endif
x=max(_QUART_*(zo(i,j)+zo(i+1,j)+z(i,j)+z(i+1,j)),-HU(i,j)+min_depth)
x=max(_HALF_*(ztmp(i,j)+ztmp(i+1,j)),-HU(i,j)+min_depth)
DU(i,j) = x+HU(i,j)
#ifdef USE_MASK
end if
......@@ -82,7 +112,7 @@
#ifdef USE_MASK
if(av(i,j) .gt. 0) then
#endif
x = max(_QUART_*(zo(i,j)+zo(i,j+1)+z(i,j)+z(i,j+1)),-HV(i,j)+min_depth)
x = max(_HALF_*(ztmp(i,j)+ztmp(i,j+1)),-HV(i,j)+min_depth)
DV(i,j) = x+HV(i,j)
#ifdef USE_MASK
end if
......@@ -91,7 +121,15 @@
end do
!$OMP END DO
update_dryingfactors = .true.
if (present(from3d)) then
if (from3d) then
update_dryingfactors = .false.
end if
end if
! KK-TODO: why not set the dry masks in any case?
! however, consider that last call is from postinit_3d
if (update_dryingfactors) then
d1 = 2*min_depth
d2i = _ONE_/(crit_depth-2*min_depth)
!$OMP DO SCHEDULE(RUNTIME)
......@@ -109,6 +147,7 @@
end do
end do
!$OMP END DO
end if
!$OMP END PARALLEL
......
......@@ -6,8 +6,17 @@
if (rc /= 0) stop 'init_2d: Error allocating memory (iteration_stat)'
#endif
allocate(D(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (D)'
allocate(t_zo(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (t_zo)'
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(t_D(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (t_D)'
allocate(DU(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (DU)'
......@@ -15,12 +24,6 @@
allocate(DV(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (DV)'
allocate(z(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (z)'
allocate(zo(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (zo)'
allocate(U(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (U)'
......
......@@ -2,9 +2,7 @@
integer,dimension(:,:),allocatable :: break_mask
integer,dimension(:,:),allocatable :: break_stat
#endif
REALTYPE,dimension(:,:),allocatable :: D
REALTYPE,dimension(:,:),allocatable,target :: DU,DV
REALTYPE,dimension(:,:),allocatable :: z,zo
REALTYPE,dimension(:,:),allocatable,target :: t_zo,t_z,t_Dlast,t_D,DU,DV
REALTYPE,dimension(:,:),allocatable :: U,V
REALTYPE,dimension(:,:),allocatable :: UEx,VEx
REALTYPE,dimension(:,:),allocatable :: fU,fV
......
......@@ -34,6 +34,15 @@
interface
subroutine depth_update(zo,z,Dlast,D,DU,DV,first,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
end subroutine depth_update
subroutine uv_advect(U,V,DU,DV)
use domain, only: imin,imax,jmin,jmax
IMPLICIT NONE
......@@ -193,7 +202,9 @@
z = -H+min_depth
end where
zo = z
call depth_update()
! 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.)
end if
if (Am .lt. _ZERO_) then
......@@ -393,7 +404,7 @@
end where
end if
call depth_update()
call depth_update(zo,z,Dlast,D,DU,DV,first=.true.)
end if
......@@ -472,7 +483,7 @@
end if
if (have_boundaries) call update_2d_bdy(loop,bdy2d_ramp)
call sealevel()
call depth_update()
call depth_update(zo,z,Dlast,D,DU,DV)
if(residual .gt. 0 .and. loop .ge. residual) then
call tic(TIM_INTEGR2D)
......
......@@ -43,6 +43,7 @@
!
! !LOCAL VARIABLES:
integer :: i,j
REALTYPE,dimension(:,:),pointer :: p2d
#ifdef USE_BREAKS
integer :: n,break_flag,break_flags(nprocs)
#endif
......@@ -59,14 +60,7 @@
#endif
call tic(TIM_SEALEVEL)
! OMP-NOTE: This loop does not really improve from threading.
! It is bound by memory access, and everything is nicely
! lined up anyway, so we keep it in serial.
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
zo(i,j) = z(i,j)
end do
end do
p2d => zo ; zo => z ; z => p2d
#ifdef USE_BREAKS
break_flag=1
......@@ -93,9 +87,9 @@
do j=jmin,jmax
do i=imin,imax
if (az(i,j) .eq. 1) then
z(i,j)=z(i,j)-dtm*((U(i,j)*DYU-U(i-1,j )*DYUIM1) &
+(V(i,j)*DXV-V(i ,j-1)*DXVJM1))*ARCD1 &
+dtm*fwf(i,j)
z(i,j)=zo(i,j)-dtm*((U(i,j)*DYU-U(i-1,j )*DYUIM1) &
+(V(i,j)*DXV-V(i ,j-1)*DXVJM1))*ARCD1 &
+dtm*fwf(i,j)
#if 0
#ifdef FRESHWATER_LENSE_TEST
kk=1.0
......@@ -152,7 +146,7 @@
end do
if(break_flag .ne. 0) then
z=zo
p2d => z ; z => zo ; zo => p2d
call tic(TIM_SEALEVELH)
call update_2d_halo(U,U,au,imin,jmin,imax,jmax,u_TAG)
call wait_halo(u_TAG)
......
......@@ -5,10 +5,7 @@
integer break_mask(E2DFIELD)
integer break_stat(E2DFIELD)
#endif
REALTYPE D(E2DFIELD)
REALTYPE,dimension(E2DFIELD),target :: DU,DV
REALTYPE z(E2DFIELD)
REALTYPE zo(E2DFIELD)
REALTYPE,dimension(E2DFIELD),target :: t_zo,t_z,t_Dlast,t_D,DU,DV
REALTYPE U(E2DFIELD)
REALTYPE V(E2DFIELD)
REALTYPE UEx(E2DFIELD)
......
......@@ -21,6 +21,8 @@
IMPLICIT NONE
!
! !PUBLIC DATA MEMBERS:
REALTYPE,dimension(:,:),pointer :: zo,z,Dlast,D
#ifdef STATIC
#include "static_2d.h"
#else
......@@ -109,6 +111,9 @@
#include "dynamic_allocations_2d.h"
#endif
zo => t_zo ; z => t_z
Dlast => t_Dlast ; D => t_D
#ifdef USE_BREAKS
break_stat = 0
#endif
......@@ -116,7 +121,7 @@
z = _ZERO_; zo =_ZERO_
zub=_ZERO_ ; zub0=_ZERO_
zvb=_ZERO_ ; zvb0=_ZERO_
D = _ZERO_;
Dlast = _ZERO_ ; D = _ZERO_
U = _ZERO_; DU = _ZERO_; fU = _ZERO_; Uint = _ZERO_; UEx = _ZERO_
V = _ZERO_; DV = _ZERO_; fV = _ZERO_; Vint = _ZERO_; VEx = _ZERO_
......
......@@ -153,24 +153,27 @@
allocate(ssen(I2DFIELD),stat=rc) ! Elevation after macro time step (z-column)
if (rc /= 0) stop 'init_3d: Error allocating memory (ssen)'
allocate(Dn(I2DFIELD),stat=rc) ! depth after macro time step (z-column)
if (rc /= 0) stop 'init_3d: Error allocating memory (Dn)'
allocate(ssuo(I2DFIELD),stat=rc) ! Elevation before macro time step (u-column)
if (rc /= 0) stop 'init_3d: Error allocating memory (ssuo)'
allocate(ssun(I2DFIELD),stat=rc) ! Elevation after macro time step (u-column)
if (rc /= 0) stop 'init_3d: Error allocating memory (ssun)'
allocate(Dun(I2DFIELD),stat=rc) ! depth after macro time step (u-column)
if (rc /= 0) stop 'init_3d: Error allocating memory (Dun)'
allocate(ssvo(I2DFIELD),stat=rc) ! Elevation before macro time step (v-column)
if (rc /= 0) stop 'init_3d: Error allocating memory (ssvo)'
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(t_Dn(I2DFIELD),stat=rc) ! depth after macro time step (z-column)
if (rc /= 0) stop 'init_3d: Error allocating memory (t_Dn)'
allocate(Dun(I2DFIELD),stat=rc) ! depth after macro time step (u-column)
if (rc /= 0) stop 'init_3d: Error allocating memory (Dun)'
allocate(Dvn(I2DFIELD),stat=rc) ! depth after macro time step (v-column)
if (rc /= 0) stop 'init_3d: Error allocating memory (Dvn)'
......
......@@ -48,10 +48,10 @@
#endif
! 2D fields in 3D domain
REALTYPE, dimension(:,:), allocatable :: sseo,ssen,Dn
REALTYPE, dimension(:,:), allocatable :: sseo,ssen
REALTYPE, dimension(:,:), allocatable :: ssuo,ssun
REALTYPE, dimension(:,:), allocatable :: ssvo,ssvn
REALTYPE,dimension(:,:),allocatable,target :: Dun,Dvn
REALTYPE,dimension(:,:),allocatable,target :: t_Dold,t_Dn,Dun,Dvn
! 3D friction in 3D domain
REALTYPE, dimension(:,:), allocatable :: rru,rrv
......
......@@ -26,7 +26,7 @@
use exceptions
use parameters, only: avmmol
use domain, only: openbdy,maxdepth,vert_cord,az
use m2d, only: uv_advect,uv_diffusion
use m2d, only: depth_update,uv_advect,uv_diffusion
use variables_2d, only: z,Uint,Vint,UEx,VEx
#ifndef NO_BAROCLINIC
use temperature,only: init_temperature, do_temperature, &
......@@ -229,6 +229,7 @@
if (.not. hotstart) then
ssen = z
call start_macro()
Dold = Dn
call coordinates(hotstart)
call hcc_check()
end if
......@@ -379,6 +380,7 @@
end do
end do
call depth_update(sseo,ssen,Dold,Dn,Dun,Dvn,first=.true.,from3d=.true.)
call coordinates(hotstart)
end if
......
......@@ -29,9 +29,10 @@
!
! !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,Dn,Dun,Dvn
use variables_3d, only: sseo,ssen,ssuo,ssun,ssvo,ssvn,Dold,Dn,Dun,Dvn
use getm_timers, only: tic, toc, TIM_STARTMCR
IMPLICIT NONE
!
......@@ -41,6 +42,7 @@
! !LOCAL VARIABLES:
integer :: i,j
REALTYPE :: split
REALTYPE,dimension(:,:),pointer :: p2d
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -51,6 +53,8 @@
#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)
......
......@@ -77,12 +77,11 @@
! 2D fields in 3D domain
REALTYPE :: sseo(I2DFIELD)
REALTYPE :: ssen(I2DFIELD)
REALTYPE :: Dn(I2DFIELD)
REALTYPE :: ssuo(I2DFIELD)
REALTYPE :: ssun(I2DFIELD)
REALTYPE :: ssvo(I2DFIELD)
REALTYPE :: ssvn(I2DFIELD)
REALTYPE,dimension(I2DFIELD),target :: Dun,Dvn
REALTYPE,dimension(I2DFIELD),target :: t_Dold,t_Dn,Dun,Dvn
! 3D friction in 3D domain
REALTYPE :: rru(I2DFIELD)
......
......@@ -115,6 +115,7 @@
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.
......@@ -191,6 +192,8 @@
#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_
......@@ -200,7 +203,7 @@
cor_v = _ZERO_ ; epg_v = _ZERO_ ; ipg_v = _ZERO_
#endif
ssen = _ZERO_ ; ssun = _ZERO_ ; ssvn = _ZERO_
Dn = _ZERO_ ; Dun = _ZERO_ ; Dvn = _ZERO_
Dold = _ZERO_ ; Dn = _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