Skip to content

GitLab

G
getm
Commit 767fb613 authored by Knut's avatar Knut
Browse files
Options

tovel

parent 798b2303
Hide whitespace changes
Inline Side-by-side
Showing with 699 additions and 354 deletions
schemas/getm-2.3.defaults
View file @ 767fb613
......@@ -475,12 +475,12 @@
<save_3d>
True
</save_3d>
<save_transports>
True
</save_transports>
<save_vel>
True
</save_vel>
<destag>
True
</destag>
<save_strho>
True
</save_strho>
......
schemas/getm-2.3.schema
View file @ 767fb613
......@@ -710,8 +710,8 @@
<element name="save_2d" type="bool" label="save 2d fields "/>
<element name="save_meteo" type="bool" label="save meteo data"/>
<element name="save_3d" type="bool" label="save 3d fields"/>
<element name="save_vel" type="bool" label="save velocities"/>
<element name="destag" type="bool" label="save velocities at T-points"/>
<element name="save_transports" type="bool" label="save grid-related transports"/>
<element name="save_vel" type="bool" label="save velocities at T-points"/>
<element name="save_strho" type="bool" label="save any of salinity, temperature, density"/>
<element name="save_s" type="bool" label="save salinity"/>
<element name="save_t" type="bool" label="save temperature"/>
......
src/2d/depth_update.F90
View file @ 767fb613
......@@ -22,7 +22,7 @@
! !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 variables_2d, only: D,Dlast,z,zo,DU,DV
use getm_timers, only: tic, toc, TIM_DPTHUPDATE
!$ use omp_lib
IMPLICIT NONE
......@@ -33,6 +33,8 @@
! !LOCAL VARIABLES:
integer :: i,j
REALTYPE :: d1,d2i,x
REALTYPE,dimension(E2DFIELD) :: ztmp
REALTYPE,dimension(:,:),pointer :: p2d
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -46,6 +48,8 @@
! 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
......@@ -55,6 +59,7 @@
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
......@@ -66,7 +71,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 +87,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
......
src/2d/dynamic_allocations_2d.h
View file @ 767fb613
......@@ -6,8 +6,11 @@
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_D(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (t_D)'
allocate(t_Dlast(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (t_Dlast)'
allocate(DU(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (DU)'
......@@ -15,11 +18,11 @@
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(t_z(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (t_z)'
allocate(zo(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (zo)'
allocate(t_zo(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (t_zo)'
allocate(U(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (U)'
......
src/2d/dynamic_declarations_2d.h
View file @ 767fb613
......@@ -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_D,t_Dlast,DU,DV,t_z,t_zo
REALTYPE,dimension(:,:),allocatable :: U,V
REALTYPE,dimension(:,:),allocatable :: UEx,VEx
REALTYPE,dimension(:,:),allocatable :: fU,fV
......
src/2d/m2d.F90
View file @ 767fb613
......@@ -194,6 +194,7 @@
end where
zo = z
call depth_update()
Dlast = D
end if
if (Am .lt. _ZERO_) then
......
src/2d/sealevel.F90
View file @ 767fb613
......@@ -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)
......
src/2d/static_2d.h
View file @ 767fb613
......@@ -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_D,t_Dlast,DU,DV,t_z,t_zo
REALTYPE U(E2DFIELD)
REALTYPE V(E2DFIELD)
REALTYPE UEx(E2DFIELD)
......
src/2d/variables_2d.F90
View file @ 767fb613
......@@ -21,6 +21,8 @@
IMPLICIT NONE
!
! !PUBLIC DATA MEMBERS:
REALTYPE,dimension(:,:),pointer :: D,Dlast,z,zo
#ifdef STATIC
#include "static_2d.h"
#else
......@@ -108,6 +110,8 @@
#ifndef STATIC
#include "dynamic_allocations_2d.h"
#endif
D => t_D ; Dlast => t_Dlast
z => t_z ; zo => t_zo
#ifdef USE_BREAKS
break_stat = 0
......@@ -116,7 +120,7 @@
z = _ZERO_; zo =_ZERO_
zub=_ZERO_ ; zub0=_ZERO_
zvb=_ZERO_ ; zvb0=_ZERO_
D = _ZERO_;
D = _ZERO_; Dlast = _ZERO_
U = _ZERO_; DU = _ZERO_; fU = _ZERO_; Uint = _ZERO_; UEx = _ZERO_
V = _ZERO_; DV = _ZERO_; fV = _ZERO_; Vint = _ZERO_; VEx = _ZERO_
......
src/3d/dynamic_allocations_3d.h
View file @ 767fb613
......@@ -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_Dn(I2DFIELD),stat=rc) ! depth after macro time step (z-column)
if (rc /= 0) stop 'init_3d: Error allocating memory (t_Dn)'
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(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)'
......
src/3d/dynamic_declarations_3d.h
View file @ 767fb613
......@@ -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_Dn,t_Dold,Dun,Dvn
! 3D friction in 3D domain
REALTYPE, dimension(:,:), allocatable :: rru,rrv
......
src/3d/m3d.F90
View file @ 767fb613
......@@ -229,6 +229,7 @@
if (.not. hotstart) then
ssen = z
call start_macro()
Dold = Dn
call coordinates(hotstart)
call hcc_check()
end if
......
src/3d/start_macro.F90
View file @ 767fb613
......@@ -31,7 +31,7 @@
use domain, only: imin,imax,jmin,jmax,H,HU,HV,min_depth
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,Dn,Dold,Dun,Dvn
use getm_timers, only: tic, toc, TIM_STARTMCR
IMPLICIT NONE
!
......@@ -41,6 +41,7 @@
! !LOCAL VARIABLES:
integer :: i,j
REALTYPE :: split
REALTYPE,dimension(:,:),pointer :: p2d
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -51,6 +52,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)
......
src/3d/static_3d.h
View file @ 767fb613
......@@ -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_Dn,t_Dold,Dun,Dvn
! 3D friction in 3D domain
REALTYPE :: rru(I2DFIELD)
......
src/3d/variables_3d.F90
View file @ 767fb613
......@@ -115,6 +115,7 @@
IMPLICIT NONE
!
! !PUBLIC DATA MEMBERS:
REALTYPE,dimension(:,:),pointer :: Dn,Dold
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
Dn => t_Dn ; Dold => t_Dold
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_
Dn = _ZERO_ ; Dold = _ZERO_ ; Dun = _ZERO_ ; Dvn = _ZERO_
rru= _ZERO_ ; rrv= _ZERO_
uuEx= _ZERO_ ; vvEx= _ZERO_
tke=1.e-10 ; eps=1.e-10
......
src/futils/Makefile
View file @ 767fb613
......@@ -13,7 +13,7 @@ LIBSRC = ver_interpol.F90 kbk_interpol.F90 tridiagonal.F90 pos.F90 \
cnv_2d.F90 cnv_3d.F90 eta_mask.F90 col_interpol.F90 \
to_2d_vel.F90 to_2d_u.F90 to_2d_v.F90 \
to_3d_vel.F90 to_3d_uu.F90 to_3d_vv.F90 \
tow.F90 c2x.F90 check_3d_fields.F90
to_velx.F90 to_vely.F90 tow.F90 c2x.F90 check_3d_fields.F90
SRC = $(MODSRC) $(LIBSRC)
......@@ -47,6 +47,8 @@ ${LIB}(to_2d_v.o) \
${LIB}(to_3d_vel.o) \
${LIB}(to_3d_uu.o) \
${LIB}(to_3d_vv.o) \
${LIB}(to_velx.o) \
${LIB}(to_vely.o) \
${LIB}(tow.o) \
${LIB}(c2x.o) \
${LIB}(check_3d_fields.o)
......
src/futils/to_velx.F90 0 → 100644
View file @ 767fb613
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: to_velx() - calculates cell centered physical vertical velocity
!
! !INTERFACE:
subroutine to_velx(imin,jmin,imax,jmax,az, &
dt,grid_type, &
#if defined(CURVILINEAR) || defined(SPHERICAL)
dxv,dyu,arcd1, &
#else
dx,dy,ard1, &
#endif
xc,xu,xv,D,Dlast,U,DU,V,DV,wwm,wwp,missing,velx)
!
! !DESCRIPTION:
!
! Here the lateral velocity in the x-y-system at cell centres
! (but velocity time stages) is calculated.
!
! !USES:
!$ use omp_lib
IMPLICIT NONE
!
! !INPUT PARAMETERS:
integer,intent(in) :: imin,jmin,imax,jmax,grid_type
integer,dimension(E2DFIELD),intent(in) :: az
REALTYPE,intent(in) :: dt
#if defined(CURVILINEAR) || defined(SPHERICAL)
REALTYPE,dimension(E2DFIELD),intent(in) :: dxv,dyu,arcd1
#else
REALTYPE,intent(in) :: dx,dy,ard1
#endif
REALTYPE,dimension(E2DFIELD),intent(in) :: xc,xu,xv,D,Dlast,U,DU,V,DV,wwm,wwp
REALTYPE,intent(in) :: missing
!
! !OUTPUT PARAMETERS:
REALTYPE,dimension(E2DFIELD),intent(out) :: velx
!
! !REVISION HISTORY:
! Original author(s): Knut Klingbeil
!
! !LOCAL VARIABLES:
integer :: i,j
REALTYPE :: dtm1
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
integer, save :: Ncall = 0
Ncall = Ncall+1
write(debug,*) 'to_velx() # ',Ncall
#endif
#ifdef SLICE_MODEL
! Note (KK): this value MUST NOT be changed !!!
j = jmax/2
#endif
dtm1 = _ONE_/dt
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP FIRSTPRIVATE(j) &
!$OMP PRIVATE(i)
select case(grid_type)
case (1,2)
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin-HALO,jmax+HALO
#endif
do i=imin-HALO+1,imax+HALO
if (az(i,j) .eq. 1) then
velx(i,j) = ( U(i-1,j) + U(i,j) ) / ( Dlast(i,j) + D(i,j) )
else
velx(i,j) = missing
end if
end do
#ifndef SLICE_MODEL
end do
#endif
!$OMP END DO
case (3)
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin-HALO+1,jmax+HALO
#endif
do i=imin-HALO+1,imax+HALO
if (az(i,j) .eq. 1) then
velx(i,j) = ( &
( &
( D(i,j) - Dlast(i,j) )*dtm1 &
+ ( wwp(i,j) - wwm(i,j) ) &
) &
*xc(i,j) &
+ ( &
U(i ,j )*xu(i ,j )*DYU &
- U(i-1,j )*xu(i-1,j )*DYUIM1 &
#ifndef SLICE_MODEL
+ V(i ,j )*xv(i ,j )*DXV &
- V(i ,j-1)*xv(i ,j-1)*DXVJM1 &
#endif
) &
*ARCD1 &
) &
/(_HALF_*(Dlast(i,j)+D(i,j)))
else
velx(i,j) = missing
end if
end do
#ifndef SLICE_MODEL
end do
#endif
!$OMP END DO
case (4)
stop 'tovel: grid_type=4 not implemented yet'
case default
stop 'tovel: invalid grid_type'
end select
!$OMP END PARALLEL
#ifdef SLICE_MODEL
velx(:,j+1) = velx(:,j)
#endif
#ifdef DEBUG
write(debug,*) 'Leaving to_velx()'
write(debug,*)
#endif
return
end subroutine to_velx
!EOC
!-----------------------------------------------------------------------
! Copyright (C) 2012 - Hans Burchard and Karsten Bolding (BBH) !
!-----------------------------------------------------------------------
src/futils/to_vely.F90 0 → 100644
View file @ 767fb613
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: to_vely() - calculates cell centered physical vertical velocity
!
! !INTERFACE:
subroutine to_vely(imin,jmin,imax,jmax,az, &
dt,grid_type, &
#if defined(CURVILINEAR) || defined(SPHERICAL)
dxv,dyu,arcd1, &
#else
dx,dy,ard1, &
#endif
yc,yu,yv,D,Dlast,U,DU,V,DV,wwm,wwp,missing,vely)
!
! !DESCRIPTION:
!
! Here the lateral velocity in the x-y-system at cell centres
! (but velocity time stages) is calculated.
!
! !USES:
!$ use omp_lib
IMPLICIT NONE
!
! !INPUT PARAMETERS:
integer,intent(in) :: imin,jmin,imax,jmax,grid_type
integer,dimension(E2DFIELD),intent(in) :: az
REALTYPE,intent(in) :: dt
#if defined(CURVILINEAR) || defined(SPHERICAL)
REALTYPE,dimension(E2DFIELD),intent(in) :: dxv,dyu,arcd1
#else
REALTYPE,intent(in) :: dx,dy,ard1
#endif
REALTYPE,dimension(E2DFIELD),intent(in) :: yc,yu,yv,D,Dlast,U,DU,V,DV,wwm,wwp
REALTYPE,intent(in) :: missing
!
! !OUTPUT PARAMETERS:
REALTYPE,dimension(E2DFIELD),intent(out) :: vely
!
! !REVISION HISTORY:
! Original author(s): Knut Klingbeil
!
! !LOCAL VARIABLES:
integer :: i,j
REALTYPE :: dtm1
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
integer, save :: Ncall = 0
Ncall = Ncall+1
write(debug,*) 'to_vely() # ',Ncall
#endif
#ifdef SLICE_MODEL
! Note (KK): this value MUST NOT be changed !!!
j = jmax/2
#endif
dtm1 = _ONE_/dt
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP FIRSTPRIVATE(j) &
!$OMP PRIVATE(i)
select case(grid_type)
case (1,2)
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin-HALO+1,jmax+HALO
#endif
do i=imin-HALO,imax+HALO
if (az(i,j) .eq. 1) then
vely(i,j) = ( V(i,j-1) + V(i,j) ) / ( Dlast(i,j) + D(i,j) )
else
vely(i,j) = missing
end if
end do
#ifndef SLICE_MODEL
end do
#endif