Commit d776175e authored by Karsten Bolding's avatar Karsten Bolding

adaptive coordinates - Hofmeister

parent ce686b07
......@@ -184,3 +184,10 @@
! that is the points that are independent of neighbours.
#define DO_EILOOP DO i=imin,imax
#define DO_EJLOOP DO j=jmin,jmax
! Definition of vertical coordinate identifiers
#define _SIGMA_COORDS_ 1
#define _Z_COORDS_ 2
#define _GENERAL_COORDS_ 3
#define _HYBRID_COORDS_ 4
#define _ADAPTIVE_COORDS_ 5
......@@ -520,4 +520,21 @@
</meanout>
</io_spec>
</getm>
<adaptcoord>
<adapt_coord>
<faclag>0.0</faclag>
<faciso>0.0</faciso>
<fachor>0.1</fachor>
<facdif>0.3</facdif>
<depthmin>0.2</depthmin>
<cNN>0.2</cNN>
<cSS>0.0</cSS>
<cdd>0.2</cdd>
<d_vel>0.1</d_vel>
<d_dens>0.5</d_dens>
<dsurf>20.</dsurf>
<tgrid>21600.</tgrid>
<preadapt>0</preadapt>
</adapt_coord>
</adaptcoord>
</scenario>
......@@ -45,6 +45,7 @@
<option value="1" label="sigma"/>
<option value="2" label="z-level"/>
<option value="3" label="general vertical coordinates (gvc)"/>
<option value="5" label="adaptive coordinates"/>
</options>
</element>
<element name="maxdepth" type="float" label="maximum depth in active calculation domain" unit="m">
......@@ -618,4 +619,23 @@
</element>
</element>
<element name="adaptcoord">
<condition type="eq" variable="getm/domain/vert_cord" value="5"/>
<element name="adapt_coord" label="adaptive coordinates parameters">
<element name="faclag" type="float" label="Lagrangian tendency" minInclusive="0" maxInclusive="1"/>
<!-- <condition type="eq" variable="../../getm/domain/vert_cord" value="5"/></element>-->
<element name="faciso" type="float" label="Isopycnal tendency" minInclusive="0" maxInclusive="1"/>
<element name="fachor" type="float" label="interf. postition filter" minInclusive="0" maxInclusive="1"/>
<element name="facdif" type="float" label="layer height filter" minInclusive="0" maxInclusive="1"/>
<element name="depthmin" type="float" label="minimum layer height" unit="m"/>
<element name="cNN" type="float" label="zoom. tendency stratification" minInclusive="0" maxInclusive="1"/>
<element name="cSS" type="float" label="zoom. tendency shear" minInclusive="0" maxInclusive="1"/>
<element name="cNN" type="float" label="zoom. tendency boundaries" minInclusive="0" maxInclusive="1"/>
<element name="d_vel" type="float" unit="m/s" label="vel. difference to resolve"/>
<element name="d_dens" type="float" unit="kg/m³" label="dens. difference to resolve"/>
<element name="dsurf" type="float" label="reference distance to boundaries" unit="m"/>
<element name="tgrid" type="float" label="adaptation timescale" unit="s"/>
<element name="preadapt" type="int" label="pre-adaptation iterations"/>
</element>
</element>
</element>
......@@ -13,7 +13,7 @@ MODSRC = m3d.F90 variables_3d.F90 advection_3d.F90 eqstate.F90 \
LIBSRC = start_macro.F90 hcc_check.F90 bdy_3d.F90 coordinates.F90 sigma_coordinates.F90 general_coordinates.F90 hybrid_coordinates.F90 check_h.F90 bottom_friction_3d.F90 internal_pressure.F90 uu_momentum_3d.F90 vv_momentum_3d.F90 ww_momentum_3d.F90 structure_friction_3d.F90 uv_advect_3d.F90 uv_diffusion_3d.F90 tke_eps_advect_3d.F90 stresses_3d.F90 ss_nn.F90 gotm.F90 slow_bottom_friction.F90 slow_advection.F90 slow_diffusion.F90 slow_terms.F90 divergence.F90 stop_macro.F90
TEXSRC = m3d.F90 variables_3d.F90 coordinates.F90 sigma_coordinates.F90 general_coordinates.F90 hybrid_coordinates.F90 check_h.F90 hcc_check.F90 uu_momentum_3d.F90 vv_momentum_3d.F90 ww_momentum_3d.F90 structure_friction_3d.F90 uv_advect_3d.F90 uv_diffusion_3d.F90 tke_eps_advect_3d.F90 bottom_friction_3d.F90 internal_pressure.F90 ip_blumberg_mellor.F90 ip_blumberg_mellor_lin.F90 ip_z_interpol.F90 ip_song_wright.F90 ip_chu_fan.F90 ip_shchepetkin_mcwilliams.F90 ip_stelling_vankester.F90 numerical_mixing.F90 physical_mixing.F90 advection_3d.F90 upstream_adv.F90 upstream_2dh_adv.F90 u_split_adv.F90 v_split_adv.F90 w_split_adv.F90 w_split_it_adv.F90 fct_2dh_adv.F90 temperature.F90 salinity.F90 spm.F90 eqstate.F90 ss_nn.F90 stresses_3d.F90 gotm.F90 rivers.F90 bdy_3d.F90 slow_bottom_friction.F90 slow_advection.F90 slow_diffusion.F90 slow_terms.F90 start_macro.F90 stop_macro.F90 divergence.F90
TEXSRC = m3d.F90 variables_3d.F90 coordinates.F90 sigma_coordinates.F90 general_coordinates.F90 hybrid_coordinates.F90 adaptive_coordinates.F90 check_h.F90 hcc_check.F90 uu_momentum_3d.F90 vv_momentum_3d.F90 ww_momentum_3d.F90 structure_friction_3d.F90 uv_advect_3d.F90 uv_diffusion_3d.F90 tke_eps_advect_3d.F90 bottom_friction_3d.F90 internal_pressure.F90 ip_blumberg_mellor.F90 ip_blumberg_mellor_lin.F90 ip_z_interpol.F90 ip_song_wright.F90 ip_chu_fan.F90 ip_shchepetkin_mcwilliams.F90 ip_stelling_vankester.F90 numerical_mixing.F90 physical_mixing.F90 advection_3d.F90 upstream_adv.F90 upstream_2dh_adv.F90 u_split_adv.F90 v_split_adv.F90 w_split_adv.F90 w_split_it_adv.F90 fct_2dh_adv.F90 temperature.F90 salinity.F90 spm.F90 eqstate.F90 ss_nn.F90 stresses_3d.F90 gotm.F90 rivers.F90 bdy_3d.F90 slow_bottom_friction.F90 slow_advection.F90 slow_diffusion.F90 slow_terms.F90 start_macro.F90 stop_macro.F90 divergence.F90
SRC = $(MODSRC) $(LIBSRC)
......@@ -60,6 +60,8 @@ $(LIB)(uv_advect_3d.o) \
$(LIB)(uv_diffusion_3d.o)
ifneq ($(GETM_NO_BAROCLINIC),true)
OBJ += \
$(LIB)(adaptive_coordinates.o) \
$(LIB)(preadapt_coordinates.o) \
$(LIB)(ip_blumberg_mellor.o) \
$(LIB)(ip_blumberg_mellor_lin.o) \
$(LIB)(ip_z_interpol.o) \
......@@ -102,6 +104,13 @@ $(MOD): $(INCS)
$(OBJ): $(MOD)
tests: test_eqstate
test_eqstate: modules objects test_eqstate.o
$(FC) -o $@ $@.o $(LDFLAGS) $(LIB)
$(RM) $@.o
./$@
doc: $(TEXSRC)
$(PROTEX) $(TEXSRC) > $(DOCDIR)/3d.tex
touch doc
......
......@@ -7,7 +7,7 @@
! \label{sec-coordinates}
!
! !INTERFACE:
subroutine coordinates(cord_type,cord_relax,maxdepth)
subroutine coordinates(hotstart)
!
! !DESCRIPTION:
!
......@@ -74,18 +74,23 @@
use variables_3d, only: kvmin,hvo,hvn
#endif
use getm_timers, only: tic, toc,TIM_COORDS
use m3d
use domain, only: vert_cord
IMPLICIT NONE
!
! !INPUT PARAMETERS:
integer, intent(in) :: cord_type
REALTYPE, intent(in) :: cord_relax
REALTYPE, intent(in) :: maxdepth
! integer, intent(in) :: cord_type
! REALTYPE, intent(in) :: cord_relax
! REALTYPE, intent(in) :: maxdepth
logical, intent(in) :: hotstart
!
! !REVISION HISTORY:
! Original author(s): Hans Burchard & Karsten Bolding
!
! !LOCAL VARIABLES:
logical, save :: first=.true.
integer :: ii
! integer :: preadapt=0
#ifdef SLICE_MODEL
integer :: i,j,k
#endif
......@@ -100,37 +105,40 @@
call tic(TIM_COORDS)
if (first) then
select case (cord_type)
case (1) ! sigma coordinates
select case (vert_cord)
case (_SIGMA_COORDS_) ! sigma coordinates
LEVEL2 'using sigma vertical coordinates'
call sigma_coordinates(.true.)
case (2) ! z-level
case (3) ! general vertical coordinates
case (_Z_COORDS_) ! z-level
case (_GENERAL_COORDS_) ! general vertical coordinates
LEVEL2 'using general vertical coordinates'
call general_coordinates(.true.,cord_relax,maxdepth)
case (4) ! hybrid vertical coordinates
case (_HYBRID_COORDS_) ! hybrid vertical coordinates
LEVEL2 'using hybrid vertical coordinates'
call hybrid_coordinates(.true.)
STDERR 'coordinates(): hybrid_coordinates not coded yet'
stop
case (5) ! adaptive vertical coordinates
case (_ADAPTIVE_COORDS_) ! adaptive vertical coordinates
LEVEL2 'using adaptive vertical coordinates'
!KB call adaptive_coordinates(.true.)
STDERR 'coordinates(): adaptive_coordinates not coded yet'
#ifndef NO_BAROCLINIC
call adaptive_coordinates(.true.,hotstart)
#endif
case default
end select
first = .false.
else
select case (cord_type)
case (1) ! sigma coordinates
select case (vert_cord)
case (_SIGMA_COORDS_) ! sigma coordinates
call sigma_coordinates(.false.)
case (2) ! z-level
case (3) ! general vertical coordinates
case (_Z_COORDS_) ! z-level
case (_GENERAL_COORDS_) ! general vertical coordinates
call general_coordinates(.false.,cord_relax,maxdepth)
case (4) ! hybrid vertical coordinates
case (_HYBRID_COORDS_) ! hybrid vertical coordinates
call hybrid_coordinates(.false.)
case (5) ! adaptive vertical coordinates
!KB call adaptive_coordinates(.false.)
case (_ADAPTIVE_COORDS_) ! adaptive vertical coordinates
#ifndef NO_BAROCLINIC
call adaptive_coordinates(.false.,hotstart)
#endif
case default
end select
end if ! first
......
......@@ -30,8 +30,9 @@
use m2d, only: Am
use variables_2d, only: D,z,UEx,VEx
#ifndef NO_BAROCLINIC
use temperature,only: init_temperature, do_temperature
use salinity, only: init_salinity, do_salinity
use temperature,only: init_temperature, do_temperature, &
init_temperature_field
use salinity, only: init_salinity, do_salinity, init_salinity_field
use eqstate, only: init_eqstate, do_eqstate
use internal_pressure, only: init_internal_pressure, do_internal_pressure
use internal_pressure, only: ip_method
......@@ -260,7 +261,7 @@
! Needed for interpolation of temperature and salinity
if (.not. hotstart) then
call start_macro()
call coordinates(vert_cord,cord_relax,maxdepth)
call coordinates(hotstart)
call hcc_check()
end if
......@@ -277,17 +278,20 @@
call init_eqstate()
#ifndef PECS
call do_eqstate()
call ss_nn()
#endif
if (.not. openbdy) bdy3d=.false.
if (bdy3d) call init_bdy_3d()
if (runtype .ge. 3) call init_internal_pressure()
if (runtype .eq. 3) call do_internal_pressure()
if (runtype .eq. 4) then
call init_advection_3d(2)
end if
end if
#endif
if (.not. openbdy) bdy3d=.false.
if (bdy3d) call init_bdy_3d()
if (vert_cord .eq. _ADAPTIVE_COORDS_) call preadapt_coordinates(preadapt)
#endif
#ifdef DEBUG
write(debug,*) 'Leaving init_3d()'
......@@ -401,7 +405,7 @@
#ifndef NO_BAROCLINIC
#endif
#ifdef MUDFLAT
call coordinates(vert_cord,cord_relax,maxdepth)
call coordinates(.false.)
#endif
#ifndef NO_BOTTFRIC
if (kmax .gt. 1) then
......@@ -435,8 +439,12 @@
call uu_momentum_3d(n,bdy3d)
ufirst=.true.
end if
#ifndef MUDFLAT
call coordinates(vert_cord,cord_relax,maxdepth)
if (kmax .gt. 1) then
if (vert_cord .eq. _ADAPTIVE_COORDS_) call ss_nn()
end if
call coordinates(.false.)
#endif
if (kmax .gt. 1) then
call ww_momentum_3d()
......@@ -451,14 +459,13 @@
#else
STDERR 'NO_ADVECT 3D'
#endif
if (kmax .gt. 1) then
#ifndef NO_BOTTFRIC
call stresses_3d()
#endif
#ifndef CONSTANT_VISCOSITY
#ifndef PARABOLIC_VISCOSITY
call ss_nn()
if (vert_cord .ne. _ADAPTIVE_COORDS_) call ss_nn()
#endif
call gotm()
#ifdef TURB_ADV
......
......@@ -27,11 +27,12 @@
private
!
! !PUBLIC DATA MEMBERS:
public init_salinity, do_salinity
public init_salinity, do_salinity, init_salinity_field
!
! !PRIVATE DATA MEMBERS:
integer :: salt_method=1,salt_format=2
character(len=PATH_MAX) :: salt_file="t_and_s.nc"
integer :: salt_field_no=1
character(len=32) :: salt_name='salt'
REALTYPE :: salt_const=35.
integer :: salt_hor_adv=1,salt_ver_adv=1
......@@ -86,12 +87,6 @@
!
! !LOCAL VARIABLES:
integer :: i,j,k,n
#ifdef PECS_TEST
integer :: cc(1:30)
#endif
integer, parameter :: nmax=100
REALTYPE :: zlev(nmax),prof(nmax)
integer :: salt_field_no=1
integer :: status
NAMELIST /salt/ &
salt_method,salt_const,salt_file, &
......@@ -105,49 +100,13 @@
! integer, save :: Ncall = 0
! Ncall = Ncall+1
! write(debug,*) 'init_salinity() # ',Ncall
#endif
#ifdef NS_FRESHWATER_LENSE_TEST
salt_field_no=1
#endif
#ifdef MED_15X15MINS_TEST
salt_field_no=1
#endif
LEVEL2 'init_salinity()'
read(NAMLST,salt)
select case (salt_method)
case(0)
LEVEL3 'getting initial fields from hotstart'
case(1)
LEVEL3 'setting to constant value'
forall(i=imin:imax,j=jmin:jmax, az(i,j) .ne. 0) &
S(i,j,:) = salt_const
case(2)
LEVEL3 'using profile'
call read_profile(salt_file,nmax,zlev,prof,n)
call ver_interpol(n,zlev,prof,imin,jmin,imax,jmax,kmax, &
az,H,hn,S)
case(3)
LEVEL3 'interpolating from 3D field'
call get_field(salt_file,salt_name,salt_field_no,S)
#ifdef SALTWEDGE_TEST
case(4)
!I need this here for hotstart of salinity!!!
LEVEL3 'initializing with #ifdef SALTWEDGE'
S = _ZERO_
do i=1,100
do k=1,kmax
S(i,2,k)=30.*(1.- tanh(float(i-1)*0.05))
end do
end do
#endif
case default
FATAL 'Not valid salt_method specified'
stop 'init_salinity'
end select
call init_salinity_field()
! Sanity checks for advection specifications
LEVEL3 'salt_hor_adv= ',salt_hor_adv
LEVEL3 'salt_ver_adv= ',salt_ver_adv
......@@ -209,6 +168,114 @@ salt_field_no=1
case default
end select
LEVEL3 'salt_check=',salt_check
if (salt_check .ne. 0) then
LEVEL4 'doing sanity check on salinity'
LEVEL4 'min_salt=',min_salt
LEVEL4 'max_salt=',max_salt
if (salt_check .gt. 0) then
LEVEL4 'out-of-bound values result in termination of program'
end if
if (salt_check .lt. 0) then
LEVEL4 'out-of-bound values result in warnings only'
end if
call check_3d_fields(imin,jmin,imax,jmax,kmin,kmax,az, &
S,min_salt,max_salt,status)
if (status .gt. 0) then
if (salt_check .gt. 0) then
call getm_error("do_salinity()", &
"out-of-bound values encountered")
end if
if (salt_check .lt. 0) then
LEVEL1 'do_salinity(): ',status, &
' out-of-bound values encountered'
end if
end if
end if
#ifdef DEBUG
write(debug,*) 'Leaving init_salinity()'
write(debug,*)
#endif
return
end subroutine init_salinity
!EOC
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: init_salinity_field - initialisation of the salinity field
! \label{sec-init-salinity}
!
! !INTERFACE:
subroutine init_salinity_field()
!
! !DESCRIPTION:
! Initialisation of the salinity field as specified by salt_method
! and exchange of the HALO zones
!
! !USES:
IMPLICIT NONE
!
! !INPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
integer :: i,j,k,n
#ifdef PECS_TEST
integer :: cc(1:30)
#endif
integer, parameter :: nmax=100
REALTYPE :: zlev(nmax),prof(nmax)
integer :: salt_field_no=1
integer :: status
!EOP
!-------------------------------------------------------------------------
!BOC
#ifdef NS_FRESHWATER_LENSE_TEST
salt_field_no=1
#endif
#ifdef MED_15X15MINS_TEST
salt_field_no=1
#endif
select case (salt_method)
case(0)
LEVEL3 'getting initial fields from hotstart'
case(1)
LEVEL3 'setting to constant value'
forall(i=imin:imax,j=jmin:jmax, az(i,j) .ne. 0) &
S(i,j,:) = salt_const
case(2)
LEVEL3 'using profile'
call read_profile(salt_file,nmax,zlev,prof,n)
call ver_interpol(n,zlev,prof,imin,jmin,imax,jmax,kmax, &
az,H,hn,S)
case(3)
LEVEL3 'interpolating from 3D field'
call get_field(salt_file,salt_name,salt_field_no,S)
#ifdef SALTWEDGE_TEST
case(4)
!I need this here for hotstart of salinity!!!
LEVEL3 'initializing with #ifdef SALTWEDGE'
S = _ZERO_
do i=1,100
do k=1,kmax
S(i,2,k)=30.*(1.- tanh(float(i-1)*0.05))
end do
end do
#endif
case default
FATAL 'Not valid salt_method specified'
stop 'init_salinity'
end select
#ifdef ARKONA_TEST
do i=100,135
do j=256,257
......@@ -267,40 +334,12 @@ salt_field_no=1
call wait_halo(D_TAG)
call mirror_bdy_3d(S,D_TAG)
LEVEL3 'salt_check=',salt_check
if (salt_check .ne. 0) then
LEVEL4 'doing sanity check on salinity'
LEVEL4 'min_salt=',min_salt
LEVEL4 'max_salt=',max_salt
if (salt_check .gt. 0) then
LEVEL4 'out-of-bound values result in termination of program'
end if
if (salt_check .lt. 0) then
LEVEL4 'out-of-bound values result in warnings only'
end if
call check_3d_fields(imin,jmin,imax,jmax,kmin,kmax,az, &
S,min_salt,max_salt,status)
if (status .gt. 0) then
if (salt_check .gt. 0) then
call getm_error("do_salinity()", &
"out-of-bound values encountered")
end if
if (salt_check .lt. 0) then
LEVEL1 'do_salinity(): ',status, &
' out-of-bound values encountered'
end if
end if
end if
#ifdef DEBUG
write(debug,*) 'Leaving init_salinity()'
write(debug,*) 'Leaving init_salinity_field()'
write(debug,*)
#endif
return
end subroutine init_salinity
end subroutine init_salinity_field
!EOC
!-----------------------------------------------------------------------
......
......@@ -259,6 +259,18 @@
!$OMP END DO
#endif
#ifdef SLICE_MODEL
do k=1,kmax-1
do i = imin,imax
if (az(i,2) .ge. 1 ) then
SS(i,3,k)=SS(i,2,k)
NN(i,3,k)=NN(i,2,k)
end if
end do
end do
#endif
!$OMP END PARALLEL
call toc(TIM_SSNN)
#ifdef DEBUG
......
......@@ -26,11 +26,12 @@
private
!
! !PUBLIC DATA MEMBERS:
public init_temperature, do_temperature
public init_temperature, do_temperature, init_temperature_field
!
! !PRIVATE DATA MEMBERS:
integer :: temp_method=1,temp_format=2
character(len=PATH_MAX) :: temp_file="t_and_s.nc"
integer :: temp_field_no=1
character(len=32) :: temp_name='temp'
REALTYPE :: temp_const=20.
integer :: temp_hor_adv=1,temp_ver_adv=1
......@@ -87,9 +88,6 @@
!
! !LOCAL VARIABLES:
integer :: k,i,j,n
integer, parameter :: nmax=10000
REALTYPE :: zlev(nmax),prof(nmax)
integer :: temp_field_no=1
integer :: status
namelist /temp/ &
temp_method,temp_const,temp_file, &
......@@ -108,35 +106,10 @@
write(debug,*) 'init_temperature() # ',Ncall
#endif
#ifdef NS_FRESHWATER_LENSE_TEST
temp_field_no=1
#endif
#ifdef MED_15X15MINS_TEST
temp_field_no=1
#endif
LEVEL2 'init_temperature()'
read(NAMLST,temp)
select case (temp_method)
case(0)
LEVEL3 'getting initial fields from hotstart'
case(1)
LEVEL3 'setting to constant value'
forall(i=imin:imax,j=jmin:jmax, az(i,j) .ne. 0) &
T(i,j,:) = temp_const
case(2)
LEVEL3 'using profile'
call read_profile(temp_file,nmax,zlev,prof,n)
call ver_interpol(n,zlev,prof,imin,jmin,imax,jmax,kmax, &
az,H,hn,T)
case(3)
LEVEL3 'interpolating from 3D field'
call get_field(temp_file,temp_name,temp_field_no,T)
case default
FATAL 'Not valid temp_method specified'
stop 'init_temperature'
end select
call init_temperature_field()
! Sanity checks for advection specifications
LEVEL3 'temp_hor_adv= ',temp_hor_adv
......@@ -200,10 +173,6 @@ temp_field_no=1
case default
end select
call update_3d_halo(T,T,az,imin,jmin,imax,jmax,kmax,D_TAG)
call wait_halo(D_TAG)
call mirror_bdy_3d(T,D_TAG)
select case (attenuation_method)
case (0)
LEVEL3 'setting attenuation coefficients to constant values:'
......@@ -287,6 +256,76 @@ temp_field_no=1
end subroutine init_temperature
!EOC
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: init_temperature_field - initialisation of temperature field
! \label{sec-init-temperature-field}
!
! !INTERFACE:
subroutine init_temperature_field()
!
! !DESCRIPTION:
! Initialise the temperature field as specified with temp_method
! and exchange the HALO zones
!
! !USES:
IMPLICIT NONE
!
! !INPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
integer :: k,i,j,n
integer, parameter :: nmax=10000
REALTYPE :: zlev(nmax),prof(nmax)
integer :: status
!EOP