adaptive vertical coordinates

INTERFACE:

    subroutine adaptive_coordinates(first,hotstart)
DESCRIPTION:

The vertical grid adaptivity is partially given by a vertical diffusion equation for the vertical layer positions, with diffusivities being proportional to shear, stratification and distance from the boundaries. In the horizontal, the grid can be smoothed with respect to $z$-levels, grid layer slope and density. Lagrangian tendency of the grid movement is supported. The adaptive terrain-following grid can be set to be an Eulerian-Lagrangian grid, a hybrid $\sigma $-$\rho$ or $\sigma $-$z$ grid and combinations of these with great flexibility. With this, internal flow structures such as thermoclines can be well resolved and followed by the grid. A set of idealised examples is presented in Hofmeister et al. (2009), which show that the introduced adaptive grid strategy reduces pressure gradient errors and numerical mixing significantly.

For the configuration of parameters, a seperate namelist file adaptcoord.inp has to be given with parameters as following:
faclag - Factor on Lagrangian coords., 0.le.faclag.le.1
facdif - Factor on thickness filter, 0.le.faclag.le.1
fachor - Factor on position filter, 0.le.faclag.le.1
cNN - dependence on stratification
cSS - dependence on shear
cdd - dep. on distance from surface and bottom
d_vel - Typical velocity difference for scaling cNN adaption
d_dens - Typical density difference for scaling cSS adaption
dsurf - reference value for surface/bottom distance [m]
tgrid - Time scale of grid adaptation [s]
preadapt - number of iterations for pre-adaptation

The parameters cNN,cSS,cdd are used for the vertical adaption and have to be less or equal 1 in sum. The difference to 1 is describing a background value which forces the coordinates back to a sigma distribution. The values ddu and ddl from the domain namelist are used for weighting the zooming to surface and bottom if cdd>0. The option preadapt allows for a pre-adaption of coordinates to the initial density field and bathymetry. The number defines the number of iterations (change coordinates, vertically advect tracer, calculate vertical gradients) used for the preadaption. The initial temperature and salinity fields are re-interpolated onto the adapted grid afterwards. USES:

    use domain, only: ga,imin,imax,jmin,jmax,kmax,H,HU,HV,az,au,av
    use variables_3d, only: dt,kmin,kumin,kvmin,ho,hn,huo,hvo,hun,hvn
    use variables_3d, only: sseo,ssen,ssuo,ssun,ssvo,ssvn
    use variables_3d, only: kmin_pmz,kumin_pmz,kvmin_pmz
    use variables_3d, only: preadapt
 
   ADAPTIVE-BEGIN
    use  parameters,  only: g,rho_0
    use variables_3d, only: uu,vv,SS
 #ifndef NO_BAROCLINIC
    use variables_3d, only: NN
    use variables_3d, only: rho
 #endif
    use domain,       only: ddu,ddl
    use halo_zones, only: update_3d_halo,wait_halo
    use halo_zones, only: H_TAG,U_TAG,V_TAG
 #if defined CURVILINEAR || defined SPHERICAL
    use domain,       only: dxv,dyu,arcd1
 #else
    use domain,       only: dx,dy,ard1
 #endif
 
  ADAPTIVE-END
 
    IMPLICIT NONE
INPUT PARAMETERS:
    logical, intent(in)                 :: first
    logical, intent(in)                 :: hotstart
OUTPUT PARAMETERS:
     integer, intent(out)                :: preadapt
REVISION HISTORY:
    Original author(s): Richard Hofmeister and Hans Burchard
LOCAL VARIABLES:
    integer         :: i,j,k,rc
    REALTYPE        :: kmaxm1
    REALTYPE        :: deltaiso
    REALTYPE, save, dimension(:),     allocatable  :: be
    REALTYPE, save, dimension(:),     allocatable  :: NNloc ! local NN vector
    REALTYPE, save, dimension(:),     allocatable  :: SSloc ! local SS vector
    REALTYPE, save, dimension(:),     allocatable  :: gaa   ! new relative coord.
    REALTYPE, save, dimension(:),     allocatable  :: gaaold! old relative coord.
    REALTYPE, save, dimension(:),     allocatable  :: aav ! total grid diffus.
    REALTYPE, save, dimension(:),     allocatable  :: avn ! NN-rel. grid diffus.
    REALTYPE, save, dimension(:),     allocatable  :: avs ! SS-rel. grid diffus.
    REALTYPE, save, dimension(:),     allocatable  :: avd ! dist.-rel. grid diff.
    REALTYPE, save, dimension(:,:,:), allocatable  :: zpos ! new pos. of z-levels
    REALTYPE, save, dimension(:,:,:), allocatable  :: zposo! old pos. of z-levels
    REALTYPE, save, dimension(:,:,:), allocatable  :: work2,work3
    REALTYPE     :: faclag=_ZERO_    ! Factor on Lagrangian coords., 0.le.faclag.le.1
    REALTYPE     :: facdif=3*_TENTH_ ! Factor on thickness filter,   0.le.faclag.le.1
    REALTYPE     :: fachor=_TENTH_   ! Factor on position filter,  0.le.faclag.le.1
    REALTYPE     :: faciso=_ZERO_    ! Factor for isopycnal tendency
    REALTYPE     :: depthmin=_ONE_/5
    REALTYPE     :: Ncrit=_ONE_/1000000
    integer      :: mhor=1    ! this number is experimental - it has to be 1 for now-
    integer      :: iw=2      ! stencil for isopycnal tendency
    REALTYPE     :: rm
    INTEGER      :: im,iii,jjj,ii
    integer      :: split=1           ! splits the vertical adaption into #split steps
    REALTYPE     :: c1ad=_ONE_/5      ! dependence on NN
    REALTYPE     :: c2ad=_ZERO_       ! dependence on SS
    REALTYPE     :: c3ad=_ONE_/5      ! distance from surface and bottom
    REALTYPE     :: c4ad=6*_TENTH_    ! background value
    REALTYPE     :: d_vel=_TENTH_     ! Typical value of absolute shear
    REALTYPE     :: d_dens=_HALF_     ! Typical value of BVF squared
    REALTYPE     :: dsurf=20*_ONE_    ! reference value for surface/bottom distance
    REALTYPE     :: tgrid=21600*_ONE_ ! Time scale of grid adaptation
    REALTYPE     :: dtgrid
    REALTYPE     :: aau(0:kmax),bu(0:kmax)
    REALTYPE     :: cu(0:kmax),du(0:kmax)
    REALTYPE     :: facupper=_ONE_
    REALTYPE     :: faclower=_ONE_
    REALTYPE     :: cNN,cSS,cdd,csum
    REALTYPE     :: cbg=6*_TENTH_
    REALTYPE     :: tfac_hor=3600*_ONE_ ! factor introducing a hor. adaption timescale = dt/tgrid_hor where tgrid_hor=tgrid/N
    integer      :: iip
 
    integer,save :: count=0
       namelist /adapt_coord/   faclag,facdif,fachor,faciso,  &
                         depthmin,Ncrit, &
                         cNN,cSS,cdd,cbg,d_vel,d_dens, &
                         dsurf,tgrid,split,preadapt
 #if (defined GETM_PARALLEL && defined INPUT_DIR)
    character(len=PATH_MAX)   :: input_dir=INPUT_DIR
 #else
    character(len=PATH_MAX)   :: input_dir='./'
 #endif