Commit fe6d52dd authored by Knut's avatar Knut

Merge branch 'postinit' into waves

Conflicts:
	src/2d/bdy_2d.F90
	src/2d/variables_2d.F90
parents 9d151262 9e917582
Pipeline #143 failed with stages
......@@ -390,16 +390,18 @@
select case (bdy_2d_type(l))
case (ZERO_GRADIENT,CLAMPED_VEL,FLATHER_VEL)
do j = wfj(n),wlj(n)
z(i,j) = z(i+1,j)
if ( az(i+1,j) .ne. 0 ) z(i,j) = z(i+1,j)
end do
case (SOMMERFELD)
do j = wfj(n),wlj(n)
if ( az(i+1,j) .ne. 0 ) then
cfl = sqrt(g*_HALF_*(D(i,j)+D(i+1,j)))*dtm/DXU
z(i,j) = ( &
(_ONE_ - _TWO_*cfl*(_ONE_-theta))*z (i ,j) &
+(_ONE_ + _TWO_*cfl*(_ONE_-theta))*zo(i+1,j) &
-(_ONE_ - _TWO_*cfl*theta )*z (i+1,j) &
)/(_ONE_ + _TWO_*cfl*theta )
end if
end do
case (CLAMPED_ELEV,CLAMPED)
do j = wfj(n),wlj(n)
......@@ -409,13 +411,16 @@
end do
case (FLATHER_ELEV)
do j = wfj(n),wlj(n)
if ( az(i+1,j) .ne. 0 ) then
! Note (KK): approximate interface depths at vel-time stage
! by spatial mean at last sse-time stage
depth = _HALF_*(D(i,j)+D(i+1,j))
! Note (KK): note approximation of sse at vel-time stage
a = ramp*bdy_data(kl) &
- _TWO_/sqrt(g*depth)*(U(i,j)-ramp*bdy_data_u(kl)*depth)
z(i,j) = max(a,-H(i,j)+min_depth)
a = _TWO_/sqrt(g*depth)*(U(i,j)-ramp*bdy_data_u(kl)*depth)
else
a = _ZERO_
end if
z(i,j) = max(ramp*bdy_data(kl)-a,-H(i,j)+min_depth)
k = k+1
kl = kl + 1
end do
......@@ -429,16 +434,18 @@
select case (bdy_2d_type(l))
case (ZERO_GRADIENT,CLAMPED_VEL,FLATHER_VEL)
do i = nfi(n),nli(n)
z(i,j) = z(i,j-1)
if ( az(i,j-1) .ne. 0 ) z(i,j) = z(i,j-1)
end do
case (SOMMERFELD)
do i = nfi(n),nli(n)
if ( az(i,j-1) .ne. 0 ) then
cfl = sqrt(g*_HALF_*(D(i,j-1)+D(i,j)))*dtm/DYVJM1
z(i,j) = ( &
(_ONE_ - _TWO_*cfl*(_ONE_-theta))*z (i,j ) &
+(_ONE_ + _TWO_*cfl*(_ONE_-theta))*zo(i,j-1) &
-(_ONE_ - _TWO_*cfl*theta )*z (i,j-1) &
)/(_ONE_ + _TWO_*cfl*theta )
end if
end do
case (CLAMPED_ELEV,CLAMPED)
do i = nfi(n),nli(n)
......@@ -448,13 +455,16 @@
end do
case (FLATHER_ELEV)
do i = nfi(n),nli(n)
if ( az(i,j-1) .ne. 0 ) then
! Note (KK): approximate interface depths at vel-time stage
! by spatial mean at last sse-time stage
depth = _HALF_*(D(i,j-1)+D(i,j))
! Note (KK): note approximation of sse at vel-time stage
a = ramp*bdy_data(kl) &
+ _TWO_/sqrt(g*depth)*(V(i,j-1)-ramp*bdy_data_v(kl)*depth)
z(i,j) = max(a,-H(i,j)+min_depth)
a = _TWO_/sqrt(g*depth)*(V(i,j-1)-ramp*bdy_data_v(kl)*depth)
else
a = _ZERO_
end if
z(i,j) = max(ramp*bdy_data(kl)+a,-H(i,j)+min_depth)
k = k+1
kl = kl + 1
end do
......@@ -468,16 +478,18 @@
select case (bdy_2d_type(l))
case (ZERO_GRADIENT,CLAMPED_VEL,FLATHER_VEL)
do j = efj(n),elj(n)
z(i,j) = z(i-1,j)
if ( az(i-1,j) .ne. 0 ) z(i,j) = z(i-1,j)
end do
case (SOMMERFELD)
do j = efj(n),elj(n)
if ( az(i-1,j) .ne. 0 ) then
cfl = sqrt(g*_HALF_*(D(i-1,j)+D(i,j)))*dtm/DXUIM1
z(i,j) = ( &
(_ONE_ - _TWO_*cfl*(_ONE_-theta))*z (i ,j) &
+(_ONE_ + _TWO_*cfl*(_ONE_-theta))*zo(i-1,j) &
-(_ONE_ - _TWO_*cfl*theta )*z (i-1,j) &
)/(_ONE_ + _TWO_*cfl*theta )
end if
end do
case (CLAMPED_ELEV,CLAMPED)
do j = efj(n),elj(n)
......@@ -487,13 +499,16 @@
end do
case (FLATHER_ELEV)
do j = efj(n),elj(n)
if ( az(i-1,j) .ne. 0 ) then
! Note (KK): approximate interface depths at vel-time stage
! by spatial mean at last sse-time stage
depth = _HALF_*(D(i-1,j)+D(i,j))
! Note (KK): note approximation of sse at vel-time stage
a = ramp*bdy_data(kl) &
+ _TWO_/sqrt(g*depth)*(U(i-1,j)-ramp*bdy_data_u(kl)*depth)
z(i,j) = max(a,-H(i,j)+min_depth)
a = _TWO_/sqrt(g*depth)*(U(i-1,j)-ramp*bdy_data_u(kl)*depth)
else
a = _ZERO_
end if
z(i,j) = max(ramp*bdy_data(kl)+a,-H(i,j)+min_depth)
k = k+1
kl = kl + 1
end do
......@@ -507,16 +522,18 @@
select case (bdy_2d_type(l))
case (ZERO_GRADIENT,CLAMPED_VEL,FLATHER_VEL)
do i = sfi(n),sli(n)
z(i,j) = z(i,j+1)
if ( az(i,j+1) .ne. 0 ) z(i,j) = z(i,j+1)
end do
case (SOMMERFELD)
do i = sfi(n),sli(n)
if ( az(i,j+1) .ne. 0 ) then
cfl = sqrt(g*_HALF_*(D(i,j)+D(i,j+1)))*dtm/DYV
z(i,j) = ( &
(_ONE_ - _TWO_*cfl*(_ONE_-theta))*z (i,j ) &
+(_ONE_ + _TWO_*cfl*(_ONE_-theta))*zo(i,j+1) &
-(_ONE_ - _TWO_*cfl*theta )*z (i,j+1) &
)/(_ONE_ + _TWO_*cfl*theta )
end if
end do
case (CLAMPED_ELEV,CLAMPED)
do i = sfi(n),sli(n)
......@@ -526,13 +543,16 @@
end do
case (FLATHER_ELEV)
do i = sfi(n),sli(n)
if ( az(i,j+1) .ne. 0 ) then
! Note (KK): approximate interface depths at vel-time stage
! by spatial mean at last sse-time stage
depth = _HALF_*(D(i,j)+D(i,j+1))
! Note (KK): note approximation of sse at vel-time stage
a = ramp*bdy_data(kl) &
- _TWO_/sqrt(g*depth)*(V(i,j)-ramp*bdy_data_v(kl)*depth)
z(i,j) = max(a,-H(i,j)+min_depth)
a = _TWO_/sqrt(g*depth)*(V(i,j)-ramp*bdy_data_v(kl)*depth)
else
a = _ZERO_
end if
z(i,j) = max(ramp*bdy_data(kl)-a,-H(i,j)+min_depth)
k = k+1
kl = kl + 1
end do
......@@ -627,6 +647,7 @@
select case (bdy_2d_type(l))
case (FLATHER_VEL)
do j = wfj(n),wlj(n)
if ( az(i+1,j) .ne. 0 ) then
! Note (KK): approximate interface depths at vel-time stage
! by spatial mean at last sse-time stage
depth = _HALF_*(D(i,j)+D(i+1,j))
......@@ -636,11 +657,13 @@
if ( waveforcing_method .ne. NO_WAVES ) then
UEuler(i,j) = U(i,j) - UStokes(i,j)
end if
end if
k = k+1
kl = kl + 1
end do
case (CLAMPED_VEL,CLAMPED)
do j = wfj(n),wlj(n)
if ( az(i+1,j) .ne. 0 ) then
! Note (KK): approximate interface depths at vel-time stage
! by spatial mean at last sse-time stage
depth = _HALF_*(D(i,j)+D(i+1,j))
......@@ -648,6 +671,7 @@
if ( waveforcing_method .ne. NO_WAVES ) then
UEuler(i,j) = U(i,j) - UStokes(i,j)
end if
end if
k = k+1
kl = kl + 1
end do
......@@ -662,6 +686,7 @@
select case (bdy_2d_type(l))
case (FLATHER_VEL)
do j = efj(n),elj(n)
if ( az(i-1,j) .ne. 0 ) then
! Note (KK): approximate interface depths at vel-time stage
! by spatial mean at last sse-time stage
depth = _HALF_*(D(i-1,j)+D(i,j))
......@@ -671,11 +696,13 @@
if ( waveforcing_method .ne. NO_WAVES ) then
UEuler(i-1,j) = U(i-1,j) - UStokes(i-1,j)
end if
end if
k = k+1
kl = kl + 1
end do
case (CLAMPED_VEL,CLAMPED)
do j = efj(n),elj(n)
if ( az(i-1,j) .ne. 0 ) then
! Note (KK): approximate interface depths at vel-time stage
! by spatial mean at last sse-time stage
depth = _HALF_*(D(i-1,j)+D(i,j))
......@@ -683,6 +710,7 @@
if ( waveforcing_method .ne. NO_WAVES ) then
UEuler(i-1,j) = U(i-1,j) - UStokes(i-1,j)
end if
end if
k = k+1
kl = kl + 1
end do
......@@ -700,6 +728,7 @@
select case (bdy_2d_type(l))
case (FLATHER_VEL)
do i = nfi(n),nli(n)
if ( az(i,j-1) .ne. 0 ) then
! Note (KK): approximate interface depths at vel-time stage
! by spatial mean at last sse-time stage
depth = _HALF_*(D(i,j-1)+D(i,j))
......@@ -709,11 +738,13 @@
if ( waveforcing_method .ne. NO_WAVES ) then
VEuler(i,j-1) = V(i,j-1) - VStokes(i,j-1)
end if
end if
k = k+1
kl = kl + 1
end do
case (CLAMPED_VEL,CLAMPED)
do i = nfi(n),nli(n)
if ( az(i,j-1) .ne. 0 ) then
! Note (KK): approximate interface depths at vel-time stage
! by spatial mean at last sse-time stage
depth = _HALF_*(D(i,j-1)+D(i,j))
......@@ -721,6 +752,7 @@
if ( waveforcing_method .ne. NO_WAVES ) then
VEuler(i,j-1) = V(i,j-1) - VStokes(i,j-1)
end if
end if
k = k+1
kl = kl + 1
end do
......@@ -735,6 +767,7 @@
select case (bdy_2d_type(l))
case (FLATHER_VEL)
do i = sfi(n),sli(n)
if ( az(i,j+1) .ne. 0 ) then
! Note (KK): approximate interface depths at vel-time stage
! by spatial mean at last sse-time stage
depth = _HALF_*(D(i,j)+D(i,j+1))
......@@ -744,11 +777,13 @@
if ( waveforcing_method .ne. NO_WAVES ) then
VEuler(i,j) = V(i,j) - VStokes(i,j)
end if
end if
k = k+1
kl = kl + 1
end do
case (CLAMPED_VEL,CLAMPED)
do i = sfi(n),sli(n)
if ( az(i,j+1) .ne. 0 ) then
! Note (KK): approximate interface depths at vel-time stage
! by spatial mean at last sse-time stage
depth = _HALF_*(D(i,j)+D(i,j+1))
......@@ -756,6 +791,7 @@
if ( waveforcing_method .ne. NO_WAVES ) then
VEuler(i,j) = V(i,j) - VStokes(i,j)
end if
end if
k = k+1
kl = kl + 1
end do
......
......@@ -168,7 +168,6 @@
integer :: i,j
integer :: elev_method=1
REALTYPE :: elev_const=_ZERO_
integer,parameter :: rk = kind(_ONE_)
character(LEN = PATH_MAX) :: elev_file='elev.nc'
namelist /m2d/ &
elev_method,elev_const,elev_file, &
......@@ -434,8 +433,8 @@
! This is only needed for proper flexible output
where (az .eq. 0)
z = -9999.0d0
zo = -9999.0d0
z = -9999._rk
zo = -9999._rk
end where
return
......
......@@ -23,6 +23,7 @@
IMPLICIT NONE
!
! !PUBLIC DATA MEMBERS:
integer, parameter :: rk = kind(_ONE_)
REALTYPE :: dtm
REALTYPE,dimension(:,:),pointer :: zo,z
logical :: deformC=.false.
......@@ -121,7 +122,6 @@
!
! !LOCAL VARIABLES:
integer :: rc
integer,parameter :: rk = kind(_ONE_)
!EOP
!-------------------------------------------------------------------------
!BOC
......@@ -174,12 +174,12 @@
end if
z = -9999*_ONE_; zo =_ZERO_
z = -9999._rk ; zo =_ZERO_
U = _ZERO_; DU = _ZERO_; Uint = _ZERO_; UEx = _ZERO_
V = _ZERO_; DV = _ZERO_; Vint = _ZERO_; VEx = _ZERO_
velx = -9999.0 ; vely = -9999.0
velx = -9999._rk ; vely = -9999._rk
where (az .gt. 0)
velx = _ZERO_ ; vely = _ZERO_
end where
......@@ -336,7 +336,6 @@
!
! !LOCAL VARIABLES:
logical :: used
integer,parameter :: rk = kind(_ONE_)
!EOP
!-----------------------------------------------------------------------
!BOC
......
......@@ -431,8 +431,8 @@
num(i,j,:) = 1.e-15
nuh(i,j,:) = 1.e-15
#ifndef NO_BAROCLINIC
S(i,j,:) = -9999.0
T(i,j,:) = -9999.0
S(i,j,:) = -9999._rk
T(i,j,:) = -9999._rk
#endif
end if
end do
......
......@@ -19,7 +19,7 @@
use domain, only: imin,jmin,imax,jmax,kmax,H,az,dry_z
!KB use get_field, only: get_3d_field
use variables_2d, only: fwf_int
use variables_3d, only: S,hn,kmin
use variables_3d, only: rk,S,hn,kmin
use meteo, only: metforcing,met_method,nudge_sss,sss
use meteo, only: METEO_CONST,METEO_FROMFILE,METEO_FROMEXT
use halo_zones, only: update_3d_halo,wait_halo,D_TAG,H_TAG
......@@ -285,8 +285,8 @@
stop 'init_salinity'
end select
S(:,:,0) = -9999*_ONE_
forall(i=imin:imax,j=jmin:jmax, az(i,j).eq.0) S(i,j,:) = -9999*_ONE_
S(:,:,0) = -9999._rk
forall(i=imin:imax,j=jmin:jmax, az(i,j).eq.0) S(i,j,:) = -9999._rk
call update_3d_halo(S,S,az,imin,jmin,imax,jmax,kmax,D_TAG)
call wait_halo(D_TAG)
......
......@@ -20,7 +20,7 @@
use domain, only: imin,jmin,imax,kmax,jmax,H,az,dry_z
use domain, only: ill,ihl,jll,jhl
use domain, only: ilg,ihg,jlg,jhg
use variables_3d, only: T,rad,hn,kmin,A,g1,g2,heatflux_net
use variables_3d, only: rk,T,rad,hn,kmin,A,g1,g2,heatflux_net
use meteo, only: metforcing,met_method,nudge_sst,sst
use meteo, only: METEO_CONST,METEO_FROMFILE,METEO_FROMEXT
!KB use get_field, only: get_3d_field
......@@ -367,8 +367,8 @@ end interface
stop 'init_temperature'
end select
T(:,:,0) = -9999*_ONE_
forall(i=imin:imax,j=jmin:jmax, az(i,j).eq.0) T(i,j,:) = -9999*_ONE_
T(:,:,0) = -9999._rk
forall(i=imin:imax,j=jmin:jmax, az(i,j).eq.0) T(i,j,:) = -9999._rk
call update_3d_halo(T,T,az,imin,jmin,imax,jmax,kmax,D_TAG)
call wait_halo(D_TAG)
......
......@@ -117,6 +117,7 @@
IMPLICIT NONE
!
! !PUBLIC DATA MEMBERS:
integer, parameter :: rk = kind(_ONE_)
REALTYPE :: dt,cnpar=0.9
REALTYPE :: avmback=_ZERO_,avhback=_ZERO_
logical :: do_numerical_analyses_3d=.false.
......@@ -213,7 +214,6 @@
!
! !LOCAL VARIABLES:
integer :: i,j, rc
integer,parameter :: rk = kind(_ONE_)
!EOP
!-------------------------------------------------------------------------
!BOC
......@@ -314,8 +314,8 @@
#endif
! must be nonzero for gotm_fabm in case of calc_temp=F
g1 = -9999*_ONE_
g2 = -9999*_ONE_
g1 = -9999._rk
g2 = -9999._rk
#ifdef STRUCTURE_FRICTION
sf = _ZERO_
......@@ -512,7 +512,7 @@
! Original author(s): Karsten Bolding & Jorn Bruggeman
!
! !LOCAL VARIABLES:
integer,parameter :: rk = kind(_ONE_)
!
!EOP
!-----------------------------------------------------------------------
!BOC
......
......@@ -183,9 +183,9 @@ stop
! calculate the z-coordinate of the cell centers
! references to mean sea level
zc(:,:,0)=-H(:,:)
zc(:,:,1)=-H(:,:) + 0.5*hn(:,:,1)
zc(:,:,1)=-H(:,:) + _HALF_*hn(:,:,1)
do k=2,kmax
zc(:,:,k)=zc(:,:,k-1)+0.5*(hn(:,:,k-1)+hn(:,:,k))
zc(:,:,k)=zc(:,:,k-1)+_HALF_*(hn(:,:,k-1)+hn(:,:,k))
end do
#ifdef SLICE_MODEL
......
......@@ -276,10 +276,11 @@
#endif
end if
#endif
call init_output_processing()
call init_les(runtype)
call init_output_processing()
call init_register_all_variables(runtype)
allocate(type_getm_host::output_manager_host)
......
......@@ -11,6 +11,12 @@
!
! !USES:
use field_manager
use variables_2d, only: register_2d_variables
use variables_3d, only: register_3d_variables
#ifdef _FABM_
use getm_fabm, only: register_fabm_variables
#endif
use output_processing, only: register_processed_variables, finalize_register_processed_variables
IMPLICIT NONE
!
! default: all is private.
......@@ -229,6 +235,7 @@
#ifdef _FABM_
call finalize_register_fabm_variables(fm)
#endif
call finalize_register_processed_variables(fm)
return
end subroutine finalize_register_all_variables
......
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: les_smagorinsky_sym - \label{les_smagorinsky_sym}
!
! !INTERFACE:
subroutine les_smagorinsky_sym(dudxC,dudxV, &
#ifndef SLICE_MODEL
dvdyC,dvdyU, &
#endif
shearX,shearU, &
AmC,AmX,AmU,AmV)
! Note (KK): keep in sync with interface in les.F90
!
! !DESCRIPTION:
!
!
! !USES:
use variables_les, only: SmagC2_2d,SmagX2_2d,SmagU2_2d,SmagV2_2d
use domain, only: imin,imax,jmin,jmax,az,ax,au,av
use domain, only: dxc,dyc,dxx,dyx,dxu,dyu,dxv,dyv
use getm_timers, only: tic,toc,TIM_SMAG2D
!$ use omp_lib
IMPLICIT NONE
!
! !INPUT PARAMETERS:
REALTYPE,dimension(E2DFIELD),intent(in) :: dudxC,dudxV
#ifndef SLICE_MODEL
REALTYPE,dimension(E2DFIELD),intent(in) :: dvdyC,dvdyU
#endif
REALTYPE,dimension(E2DFIELD),intent(in) :: shearX,shearU
!
! !OUTPUT PARAMETERS:
REALTYPE,dimension(E2DFIELD),intent(out),optional :: AmC,AmX,AmU,AmV
!
! !REVISION HISTORY:
! Original author(s): Knut Klingbeil
!
! !LOCAL VARIABLES:
REALTYPE :: dudx,dvdy
integer :: i,j
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
integer, save :: Ncall = 0
Ncall = Ncall+1
write(debug,*) 'les_smagorinsky_sym() # ',Ncall
#endif
#ifdef SLICE_MODEL
j = jmax/2 ! this MUST NOT be changed!!!
#endif
call tic(TIM_SMAG2D)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP FIRSTPRIVATE(j) &
!$OMP PRIVATE(i,dudx,dvdy)
if (present(AmC)) then
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin-1,jmax+1
#endif
do i=imin-1,imax+1
! Note (KK): AmC(az=1) needed in uv_diffusion
if (az(i,j) .eq. 1) then
! interpolate shearC
! Note (KK): in W/E open boundary cells shearC(az=2) would
! require shearU outside open boundary
! in N/S open boundary cells shearC(az=2) would
! require shearU(au=3)
! however shearC(az=2) not needed
AmC(i,j) = ( &
dudxC(i,j) &
#ifndef SLICE_MODEL
- dvdyC(i,j) &
#endif
)**2 &
+ (_HALF_*(shearU(i-1,j) + shearU(i,j)))**2
AmC(i,j) = SmagC2_2d(i,j)*DXC*DYC*sqrt(AmC(i,j))
end if
end do
#ifndef SLICE_MODEL
end do
#endif
!$OMP END DO NOWAIT
#ifdef SLICE_MODEL
!$OMP BARRIER
!$OMP SINGLE
AmC(imin-1:imax+1,j+1) = AmC(imin-1:imax+1,j)
!$OMP END SINGLE
#endif
end if
if (present(AmX)) then
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin-1,jmax
#endif
do i=imin-1,imax
if (ax(i,j) .eq. 1) then
! interpolate dudxX and dvdyX
if (av(i,j).eq.3 .or. av(i+1,j).eq.3) then
! Note (KK): western/eastern open bdy
dudx = _ZERO_
else
dudx = _HALF_*(dudxV(i,j) + dudxV(i+1,j ))