diff --git a/.gitignore b/.gitignore
index c3b074a7a2bf1db27187ab2b651ebcd2f4126599..f531e9199adb33270686c22a61c239834ce0b555 100644
--- a/.gitignore
+++ b/.gitignore
@@ -5,6 +5,8 @@
*.so
*.dylib
*.vtu
+*.nc
lib
include
vis
+data
diff --git a/amatos2d/compile/Makefiles/Makefile.common b/amatos2d/compile/Makefiles/Makefile.common
index 4f106850b25625ec177cb2dac61ac9e0147d91b9..01f1cc20542256180b97aa9867fd1ed34e5823f0 100644
--- a/amatos2d/compile/Makefiles/Makefile.common
+++ b/amatos2d/compile/Makefiles/Makefile.common
@@ -46,6 +46,12 @@ ifeq ($(strip $(SYSTEM)), gfortran)
SEROBJ+= MISC_wrapper.o
endif
+AMATLIBCALL = -lamatos
+ifeq ($(strip $(MAKETHING)), SAMATOS)
+ AMATLIBCALL = -lsamatos
+endif
+
+
# -------- include additional Macros ----------------------------#
FFLAGS += $(MACROS:%=$(MF)%)
CFLAGS += $(MACROS:%=$(MF)%)
@@ -103,7 +109,7 @@ libncugrid.a: IO_ncugrid.o
$(NCUSLIB): $(NCULIB)
@echo "make: Creating dynamic shared object (dso) library"
- $(LOADER) $(SHFLAGS) -o $(NCUSLIB) IO_ncugrid.o -L$(LIBDIR) -lamatos $(NETCDFLIB)
+ $(LOADER) $(SHFLAGS) -o $(NCUSLIB) IO_ncugrid.o -L$(LIBDIR) $(AMATLIBCALL) $(NETCDFLIB)
testncugrid: test_io_ncugrid.o
@echo "make: Linking object modules and libraries"
diff --git a/amatos2d/compile/Makefiles/Makefile_macosxGfortran_amatos b/amatos2d/compile/Makefiles/Makefile_macosxGfortran_amatos
index b21840771cf935892d1867e78a4e14ffd76289a7..7ed8a01524a728fbd925d84a4e3636b689609ffa 100644
--- a/amatos2d/compile/Makefiles/Makefile_macosxGfortran_amatos
+++ b/amatos2d/compile/Makefiles/Makefile_macosxGfortran_amatos
@@ -71,8 +71,8 @@ endif
BITMAPLIB =
BITMAPINC =
ifeq ($(strip $(IO_PGM)),yes)
- BITMAPLIB = -L/sw2/lib -lnetpbm
- BITMAPINC = -I/sw2/include
+ BITMAPLIB = -L/opt/local/lib -lnetpbm
+ BITMAPINC = -I/opt/local/include
endif
# SET NETCDF PATHS
@@ -81,8 +81,8 @@ NETCDFLIB =
NETCDFINC =
NCUGRIDLIB =
ifneq ($(strip $(NO_NETCDF)), yes)
- NETCDFLIB = -L/sw2/lib -lnetcdff -lnetcdf
- NETCDFINC = -I/sw2/include
+ NETCDFLIB = -L/opt/local/lib -lnetcdff -lnetcdf
+ NETCDFINC = -I/opt/local/include
NCUGRIDLIB = -lncugrid
endif
diff --git a/amatos2d/compile/Makefiles/Makefile_macosxGfortran_amatos.alternative b/amatos2d/compile/Makefiles/Makefile_macosxGfortran_amatos.alternative
new file mode 100644
index 0000000000000000000000000000000000000000..a54b7e75438ce322af9c1533151eb9f14604ccf5
--- /dev/null
+++ b/amatos2d/compile/Makefiles/Makefile_macosxGfortran_amatos.alternative
@@ -0,0 +1,230 @@
+##################################################################
+# AMATOS #
+# Adaptive Mesh generator for #
+# ATmospheric and Oceanic Simulation #
+##################################################################
+# makefile to build AMATOS / library objects of amatos API #
+# j. behrens 2/95, 3/96, 5/97, 12/2003 #
+# l. mentrup 1/2005 #
+# --- this is for Mac OS X on Intel architecture --- #
+# --- using the intel Fortran compiler --- #
+##################################################################
+
+# MACHINE
+MACHINE = macosx_gfortran
+
+# SYSTEM
+SYSTEM = gfortran
+
+# SET MAKETHING CORRESPONDING TO MACHINE:
+MAKETHING= AMATOS
+AMATFLAG = $(MF)$(MAKETHING)
+
+SAVETHING= SAVETEST
+
+# OPTIMIZATION SETTINGS [debug|opt|norm|opt]
+MODE := debug
+
+# CREATE SHARED LIBRARIES
+CREATE_SHARED := no #yes
+
+# FURTHER OPTIONS
+DBL_DBL := yes
+
+# PARALLELISM SETTINGS
+OMP := no
+MPI := no
+
+# LIBRARY SETTINGS [yes|no]
+# Usage: $(MAKE) NO_NETCDF=yes
+NO_NETCDF := no
+
+# LIBRARY SETTINGS [yes|no]
+# Usage: make IO_EMIT=yes IO_PGM=yes
+IO_EMIT := no
+IO_PGM := yes
+
+# ADDITIONAL MACRO DEFS
+MACROS =
+
+# SET MAIN DIRECTORY PATH
+# !! This has to be alterd by user !!
+ROOTDIR = $(HOME)/Documents/Teaching/WS2021/Transport-20/tracertransportsoftware
+
+# SET atlas/blas DIRECTORY PATH
+# !! This has to be alterd by user !!
+BLASLIB = -L/opt/local/lib/lapack -lblas
+
+# SET LAPACK DIRECTORY PATH
+# !! This has to be alterd by user !!
+LAPACKLIB = -L/opt/local/lib/lapack -llapack \
+ -L/$(ROOTDIR)/amatos2d/3rdparty/LAPACK95 -llapack95
+
+# SET IO_EMIT LIBRARY PATH
+# !! This has to be alterd by user !!
+EMITLIB =
+ifeq ($(strip $(IO_EMIT)),yes)
+ EMITLIB = -lstdc++ #-lresolv -lsocket -lnsl -lCstd
+endif
+
+# SET BITMAP LIBRARY PATH
+# !! This has to be alterd by user !!
+BITMAPLIB =
+BITMAPINC =
+ifeq ($(strip $(IO_PGM)),yes)
+ BITMAPLIB = -L/opt/local/lib -lnetpbm
+ BITMAPINC = -I/opt/local/include
+endif
+
+# SET NETCDF PATHS
+# !! This has to be alterd by user !!
+NETCDFLIB =
+NETCDFINC =
+NCUGRIDLIB =
+ifneq ($(strip $(NO_NETCDF)), yes)
+ NETCDFLIB = -L/opt/local/lib -lnetcdff -lnetcdf
+ NETCDFINC = -I/opt/local/include
+ NCUGRIDLIB = -lncugrid
+endif
+
+# SET MORE DIRECTORY PATHS
+MAINDIR = $(ROOTDIR)/amatos2d
+
+BUILDIR = $(MAINDIR)/compile/$(MACHINE)
+LIBDIR = $(ROOTDIR)/lib/$(MACHINE)
+INCDIR = $(ROOTDIR)/include/$(MACHINE)
+DATDIR = $(MAINDIR)/data
+SRCDIR = $(MAINDIR)/src
+
+TSTDIR = $(SRCDIR)/test
+GRIDDIR = $(SRCDIR)/gridgen
+SYSDIR = $(SRCDIR)/system/$(SYSTEM)
+TIMDIR = $(SRCDIR)/timing
+
+GRIDLIB = libamatos.a
+GRDSLIB = libamatos.dylib
+NCULIB = libncugrid.a
+NCUSLIB = libncugrid.dylib
+MODEND = mod
+GRDMOD = *.$(MODEND)
+
+#----------------------------------------------------------------#
+# library and include paths #
+#----------------------------------------------------------------#
+LIB_SH = $(LAPACKLIB) $(EMITLIB)
+LIBS = -L$(LIBDIR) -lamatos $(NCUGRIDLIB) $(NETCDFLIB) \
+ $(LAPACKLIB) $(BLASLIB) $(BITMAPLIB)
+INCPATH = $(NETCDFINC) -I$(INCDIR) -I.
+INCPATH_CC = $(BITMAPINC)
+
+#----------------------------------------------------------------#
+# FLAGS FOR LINUX / gnu Fortran Compiler #
+#----------------------------------------------------------------#
+F90 = gfortran
+cc = gcc
+CC = g++
+LOADER = gfortran
+AR = ar
+CP = cp
+CPFLAGS =
+ARFLAGS = vru
+# Compiler flag for macros
+MF = -D
+
+# --------------------- next are for debugging ------------------#
+ifeq ($(strip $(MODE)),debug)
+ FFLAGS = -fbounds-check -ggdb -fpic -funderscoring # -C
+ CFLAGS = -ggdb -fpic -I/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include
+ LDFLAGS = -ggdb -mmacosx-version-min=10.11
+ SHFLAGS = -fpic -dynamiclib -mmacosx-version-min=10.11
+endif
+# --------------------- next are for optimized debugging --------#
+ifeq ($(strip $(MODE)),optdebug)
+ FFLAGS = -fPIC -funderscoring # -C
+ CFLAGS = -g -fpic
+ LDFLAGS = -g
+ SHFLAGS = -fPIC -shared
+endif
+# --------------------- next are for normal compilation ---------#
+ifeq ($(strip $(MODE)),norm)
+ FFLAGS = -fPIC -funderscoring
+ CFLAGS = -fpic
+ LDFLAGS =
+ SHFLAGS = -fPIC -shared
+endif
+# --------------------- next are for optimization ---------------#
+ifeq ($(strip $(MODE)),opt)
+ FFLAGS = -O3 -fPIC -funderscoring
+ CFLAGS = -O3 -fpic
+ LDFLAGS = -O3
+ SHFLAGS = -fPIC -shared
+endif
+
+# ----- next flags identify the fortran compiler ----------------#
+FFLAGS += $(MF)$(F90)
+CFLAGS += $(MF)$(F90) -D$(MACHINE)
+
+# --------------------- next are for Module IO_EMIT -------------#
+ifeq ($(strip $(IO_EMIT)),yes)
+ LIB_SH += $(EMITLIB)
+ FFLAGS += $(MF)IO_EMIT
+ CFLAGS += $(MF)IO_EMIT
+endif
+# --------------------- next are for Module IO_PGM --------------#
+ifeq ($(strip $(IO_PGM)),yes)
+ LIBS += $(BITMAPLIB)
+ FFLAGS += $(MF)IO_PGM
+ CFLAGS += $(MF)IO_PGM
+endif
+# --------------------- next are for NetCDF support -------------#
+ifeq ($(strip $(NO_NETCDF)),yes)
+ FFLAGS += $(MF)NO_NETCDF
+endif
+# --------------------- next are for OpenMP ---------------------#
+ifeq ($(strip $(OMP)),yes)
+ #LIBS +=
+ #FFLAGS +=
+ #CFLAGS +=
+endif
+# --------------------- next are for MPI ------------------------#
+ifeq ($(strip $(MPI)),yes)
+ #LIBS +=
+ #FFLAGS +=
+ #CFLAGS +=
+endif
+# --------------------- next are for double precision -----------#
+ifeq ($(strip $(DBL_DBL)), yes)
+ FFLAGS += -DDBL_DBL
+endif
+
+#----------------------------------------------------------------#
+# additional clear/clean-Items #
+#----------------------------------------------------------------#
+CLEARSRC_EXT =
+CLEAREX_EXT = $(SAVETHING) *.dylib
+CLEARDAT_EXT = $(SAVETHING:%=%_*)
+CLEARLIB_EXT =
+CLEAN_EXT = *.vtu
+TIDY_EXT = *.i90 *.s
+
+#----------------------------------------------------------------#
+# common stuff #
+#----------------------------------------------------------------#
+
+include $(MAINDIR)/compile/Makefiles/Makefile.common
+
+#----------------------------------------------------------------#
+# copy source files #
+#----------------------------------------------------------------#
+
+include $(MAINDIR)/compile/Makefiles/Makefile.cpsrc
+
+#----------------------------------------------------------------#
+# test programs #
+#----------------------------------------------------------------#
+
+include $(MAINDIR)/compile/Makefiles/Makefile.TestsAmatos
+
+#----------------------------------------------------------------#
+# DEPENDENCIES ON INCLUDE FILES #
+#----------------------------------------------------------------#
diff --git a/amatos2d/compile/macosx_gfortran/Makefile b/amatos2d/compile/macosx_gfortran/Makefile
index a54b7e75438ce322af9c1533151eb9f14604ccf5..acddd0255fe782ee1d741637daf9b4dbae1cbbfa 100644
--- a/amatos2d/compile/macosx_gfortran/Makefile
+++ b/amatos2d/compile/macosx_gfortran/Makefile
@@ -17,7 +17,7 @@ MACHINE = macosx_gfortran
SYSTEM = gfortran
# SET MAKETHING CORRESPONDING TO MACHINE:
-MAKETHING= AMATOS
+MAKETHING= SAMATOS
AMATFLAG = $(MF)$(MAKETHING)
SAVETHING= SAVETEST
@@ -26,7 +26,7 @@ SAVETHING= SAVETEST
MODE := debug
# CREATE SHARED LIBRARIES
-CREATE_SHARED := no #yes
+CREATE_SHARED := yes
# FURTHER OPTIONS
DBL_DBL := yes
@@ -49,16 +49,15 @@ MACROS =
# SET MAIN DIRECTORY PATH
# !! This has to be alterd by user !!
-ROOTDIR = $(HOME)/Documents/Teaching/WS2021/Transport-20/tracertransportsoftware
+ROOTDIR = $(HOME)/Documents/Teaching/WS2021/Transport-20/TracerTransportSoftware
# SET atlas/blas DIRECTORY PATH
# !! This has to be alterd by user !!
-BLASLIB = -L/opt/local/lib/lapack -lblas
+BLASLIB =
# SET LAPACK DIRECTORY PATH
# !! This has to be alterd by user !!
-LAPACKLIB = -L/opt/local/lib/lapack -llapack \
- -L/$(ROOTDIR)/amatos2d/3rdparty/LAPACK95 -llapack95
+LAPACKLIB = -framework Accelerate
# SET IO_EMIT LIBRARY PATH
# !! This has to be alterd by user !!
@@ -101,8 +100,8 @@ GRIDDIR = $(SRCDIR)/gridgen
SYSDIR = $(SRCDIR)/system/$(SYSTEM)
TIMDIR = $(SRCDIR)/timing
-GRIDLIB = libamatos.a
-GRDSLIB = libamatos.dylib
+GRIDLIB = libsamatos.a
+GRDSLIB = libsamatos.dylib
NCULIB = libncugrid.a
NCUSLIB = libncugrid.dylib
MODEND = mod
@@ -112,7 +111,7 @@ GRDMOD = *.$(MODEND)
# library and include paths #
#----------------------------------------------------------------#
LIB_SH = $(LAPACKLIB) $(EMITLIB)
-LIBS = -L$(LIBDIR) -lamatos $(NCUGRIDLIB) $(NETCDFLIB) \
+LIBS = -L$(LIBDIR) -lsamatos $(NCUGRIDLIB) $(NETCDFLIB) \
$(LAPACKLIB) $(BLASLIB) $(BITMAPLIB)
INCPATH = $(NETCDFINC) -I$(INCDIR) -I.
INCPATH_CC = $(BITMAPINC)
@@ -134,7 +133,7 @@ MF = -D
# --------------------- next are for debugging ------------------#
ifeq ($(strip $(MODE)),debug)
FFLAGS = -fbounds-check -ggdb -fpic -funderscoring # -C
- CFLAGS = -ggdb -fpic -I/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include
+ CFLAGS = -ggdb -fpic
LDFLAGS = -ggdb -mmacosx-version-min=10.11
SHFLAGS = -fpic -dynamiclib -mmacosx-version-min=10.11
endif
@@ -223,7 +222,7 @@ include $(MAINDIR)/compile/Makefiles/Makefile.cpsrc
# test programs #
#----------------------------------------------------------------#
-include $(MAINDIR)/compile/Makefiles/Makefile.TestsAmatos
+include $(MAINDIR)/compile/Makefiles/Makefile.TestsSamatos
#----------------------------------------------------------------#
# DEPENDENCIES ON INCLUDE FILES #
diff --git a/amatos2d/compile/macosx_intel/Makefile b/amatos2d/compile/macosx_intel/Makefile
index 2d5716749b86ffa6af8a4ca029dda6785650b409..ce87c760220e6632960491c7e85b89e9da7742f1 100644
--- a/amatos2d/compile/macosx_intel/Makefile
+++ b/amatos2d/compile/macosx_intel/Makefile
@@ -49,7 +49,7 @@ MACROS =
# SET MAIN DIRECTORY PATH
# !! This has to be alterd by user !!
-ROOTDIR = $(HOME)/Documents/Development/amatos
+ROOTDIR = $(HOME)/Documents/Teaching/WS2021/Transport-20/TracerTransportSoftware
# SET atlas/blas DIRECTORY PATH
# !! This has to be alterd by user !!
@@ -72,8 +72,8 @@ endif
BITMAPLIB =
BITMAPINC =
ifeq ($(strip $(IO_PGM)),yes)
- BITMAPLIB = -L/sw2/lib -lnetpbm
- BITMAPINC = -I/sw2/include
+ BITMAPLIB = -L/opt/local/lib -lnetpbm
+ BITMAPINC = -I/opt/local/include
endif
# SET NETCDF PATHS
@@ -88,11 +88,11 @@ ifneq ($(strip $(NO_NETCDF)), yes)
endif
# SET MORE DIRECTORY PATHS
-MAINDIR = $(ROOTDIR)/amatos2d/trunk
+MAINDIR = $(ROOTDIR)/amatos2d
BUILDIR = $(MAINDIR)/compile/$(MACHINE)
-LIBDIR = $(MAINDIR)/lib/$(MACHINE)
-INCDIR = $(MAINDIR)/include/$(MACHINE)
+LIBDIR = $(ROOTDIR)/lib/$(MACHINE)
+INCDIR = $(ROOTDIR)/include/$(MACHINE)
DATDIR = $(MAINDIR)/data
SRCDIR = $(MAINDIR)/src
#SRCDIR = $(ROOTDIR)/amatos2d/branches/FixPeriodicBC/src/
@@ -113,8 +113,8 @@ GRDMOD = *.$(MODEND)
# library and include paths #
#----------------------------------------------------------------#
LIB_SH = $(LAPACKLIB) $(EMITLIB)
-LIBS = -L$(LIBDIR) -lamatos $(NCUGRIDLIB) $(NETCDFLIB) \
- $(LAPACKLIB) $(BLASLIB) $(BITMAPLIB)
+LIBS = $(NCUGRIDLIB) $(NETCDFLIB) \
+ $(LAPACKLIB) $(BLASLIB) $(BITMAPLIB) -L$(LIBDIR) -lamatos
INCPATH = $(NETCDFINC) -I$(INCDIR) -I.
INCPATH_CC = $(BITMAPINC)
diff --git a/flash2d/compile/Makefiles/Makefile.common b/flash2d/compile/Makefiles/Makefile.common
index c81b7d589445c6bfbadfd015f44ebe264960afb1..49b85db9c8d1580cc068323ecb240e5d1367ccd3 100644
--- a/flash2d/compile/Makefiles/Makefile.common
+++ b/flash2d/compile/Makefiles/Makefile.common
@@ -8,7 +8,6 @@
# OBJECTS #
#----------------------------------------------------------------#
MAINOBJ = \
-FLASH_metadata.o \
FLASH_parameters.o \
MISC_timing.o \
MISC_system.o \
@@ -100,6 +99,10 @@ ShallowWaterCopy::
@cp $(DATDIR)/Triang.SWCosine.dat Triang.dat
@cp $(DATDIR)/Parameters.SW.dat Parameters.dat
+StackedCopy:
+ @make datasphere
+ @cp $(DATDIR)/Parameters_stacked.dat Parameters.dat
+
#----------------------------------------------------------------#
# THIS CREATES PREDEFINED OPTIONS #
#----------------------------------------------------------------#
@@ -240,6 +243,16 @@ ShallowWater:
@cp $(OPTDIR)/SLM_initial.SWCosine.f90 SLM_initial.f90
@cp $(OPTDIR)/IO_vtuplot.SW.F90 IO_vtuplot.F90
+Stacked:
+ @make maincopy
+ @cp $(OPTDIR)/FLASH_parameters_stacked.f90 FLASH_parameters.f90
+ @cp $(OPTDIR)/ADV_rhs_stacked.f90 ADV_rhs.f90
+ @cp $(OPTDIR)/ADV_semilagrange_stacked.F90 ADV_semilagrange.F90
+ @cp $(OPTDIR)/ADV_wind_stacked.f90 ADV_wind.f90
+ @cp $(OPTDIR)/SLM_advanced_stacked.f90 SLM_advanced.f90
+ @cp $(OPTDIR)/SLM_errorestimate_stacked.f90 SLM_errorestimate.f90
+ @cp $(OPTDIR)/SLM_simple_stacked.f90 SLM_simple.f90
+
#----------------------------------------------------------------#
# THIS COMPILES THE MAIN PROGRAM #
#----------------------------------------------------------------#
diff --git a/flash2d/compile/macosx_gfortran/Makefile b/flash2d/compile/macosx_gfortran/Makefile
index 7b58fd417218f1e963dcf9e936e5b07a145a3fca..a813df455f350119492527ecf8c8943be1d44bab 100644
--- a/flash2d/compile/macosx_gfortran/Makefile
+++ b/flash2d/compile/macosx_gfortran/Makefile
@@ -31,7 +31,7 @@ NO_MPSLM := yes
# SET MAIN DIRECTORY PATH
# !! This has to be alterd by user !!
-ROOTDIR = $(HOME)/Documents/Development/amatos
+ROOTDIR = $(HOME)/Documents/Teaching/WS2021/Transport-20/TracerTransportSoftware
# SET atlas/blas DIRECTORY PATH
# !! This has to be alterd by user !!
@@ -49,10 +49,10 @@ CCLIBDIR = -L/usr/lib/
CCLIB = -lstdc++
# SET MORE DIRECTORY PATHS
-LIBDIR = $(ROOTDIR)/amatos2d/trunk/lib/$(MACHINE)
-INCDIR = $(ROOTDIR)/amatos2d/trunk/include/$(MACHINE)
-MODDIR = $(ROOTDIR)/amatos2d/trunk/include/$(MACHINE)
-MAINDIR = $(ROOTDIR)/flash2d/trunk
+LIBDIR = $(ROOTDIR)/lib/$(MACHINE)
+INCDIR = $(ROOTDIR)/include/$(MACHINE)
+MODDIR = $(ROOTDIR)/include/$(MACHINE)
+MAINDIR = $(ROOTDIR)/flash2d/
SRCDIR = $(MAINDIR)/src/flash-sphere
OPTDIR = $(MAINDIR)/src/options-sphere
@@ -65,8 +65,8 @@ BUILDIR = $(MAINDIR)/compile/$(MACHINE)
LIBNETCDF =
INCNETCDF =
ifneq ($(strip $(NO_NETCDF)), yes)
-LIBNETCDF = -L$(LIBDIR) -lncugrid -L/sw2/lib -lnetcdff -lnetcdf
-INCNETCDF = -I/sw2/include
+LIBNETCDF = -L$(LIBDIR) -lncugrid -L/opt/local/lib -lnetcdff -lnetcdf
+INCNETCDF = -I/opt/local/include
endif
# SET VISNET DIRECTORY PATH
@@ -93,19 +93,19 @@ MODEND = mod
F90 = gfortran
cc = gcc
CC = g++
-LOADER = gfortran
+LOADER = /usr/bin/gcc
# --------------------- next are for debugging ------------------#
ifeq ($(strip $(MODE)),debug)
- FFLAGS = -fbounds-check -ggdb -fpic -funderscoring # -C
+ FFLAGS = -fbounds-check -ggdb #-fpic -funderscoring -C
CFLAGS = -ggdb -fpic
- LDFLAGS = -ggdb
+ LDFLAGS = -v -L/opt/local/lib/libgcc -lgfortran
endif
# --------------------- next are for normal compilation ---------#
ifeq ($(strip $(MODE)),norm)
FFLAGS = -fPIC -funderscoring
CFLAGS = -fpic
- LDFLAGS =
+ LDFLAGS = -v
endif
# --------------------- next are for optimization ---------------#
ifeq ($(strip $(MODE)),opt)
diff --git a/flash2d/src/options-sphere/ADV_rhs_stacked.f90 b/flash2d/src/options-sphere/ADV_rhs_stacked.f90
new file mode 100644
index 0000000000000000000000000000000000000000..ea76ba6b98c8c89743049f3c442270cab82f39a1
--- /dev/null
+++ b/flash2d/src/options-sphere/ADV_rhs_stacked.f90
@@ -0,0 +1,105 @@
+!*****************************************************************
+!
+! MODULE NAME:
+! ADV_rhs
+! FUNCTION:
+! calculate the (nonhomogeneous) right hand side
+! CONTAINS:
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_righthand
+! FUNCTION:
+! calculate the rhs of the advection equation
+! SYNTAX:
+! real= slm_righthand(real.arr, real)
+! ON INPUT:
+! r_coord: coordinates of point real
+! r_time: time coordinate (optional) real
+! ON OUTPUT:
+! r_rhs: right hand side value real
+! CALLS:
+!
+! COMMENTS:
+!
+!-----------------------------------------------------------------
+!
+! PUBLIC:
+!
+! COMMENTS:
+! this is the homogeneous case!
+! USES:
+! MISC_globalparam, MISC_error
+! LIBRARIES:
+!
+! REFERENCES:
+!
+! VERSION(S):
+! 1. original version j. behrens 2/98
+! 2. compliant to amatos 1.0 j. behrens 12/2000
+! 3. compliant to amatos 1.2 j. behrens 3/2002
+!
+!*****************************************************************
+ MODULE ADV_rhs
+ USE GRID_api
+ PRIVATE
+ PUBLIC :: slm_righthand, slm_verticalrhs
+ CONTAINS
+!*****************************************************************
+ FUNCTION slm_righthand(r_coord, i_layer, r_time) RESULT (r_rhs)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension), INTENT(in) :: r_coord
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_layer
+ REAL (KIND = GRID_SR), INTENT(in), OPTIONAL :: r_time
+ REAL (KIND = GRID_SR) :: r_rhs
+ REAL (KIND = GRID_SR) :: r_tim
+
+!---------- set time
+
+ IF(present(r_time)) THEN
+ r_tim= r_time
+ ELSE
+ r_tim= 0.0
+ END IF
+
+!---------- calculate the advection at (x,y) (velocity increasing)
+
+ r_rhs= 0.0
+
+ RETURN
+ END FUNCTION slm_righthand
+
+!*****************************************************************
+ FUNCTION slm_verticalrhs(r_coord, i_layer, r_time) RESULT (r_rhs)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension), INTENT(in) :: r_coord
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_layer
+ REAL (KIND = GRID_SR), INTENT(in), OPTIONAL :: r_time
+ REAL (KIND = GRID_SR) :: r_rhs
+ REAL (KIND = GRID_SR) :: r_tim
+
+!---------- set time
+
+ IF(present(r_time)) THEN
+ r_tim= r_time
+ ELSE
+ r_tim= 0.0
+ END IF
+
+!---------- calculate the advection at (x,y) (velocity increasing)
+
+ r_rhs= 0.0
+
+ RETURN
+ END FUNCTION slm_verticalrhs
+
+!*****************************************************************
+ END MODULE ADV_rhs
diff --git a/flash2d/src/options-sphere/ADV_semilagrange_stacked.F90 b/flash2d/src/options-sphere/ADV_semilagrange_stacked.F90
new file mode 100644
index 0000000000000000000000000000000000000000..ea1c1ed5307be0288efe1f040f5289fe3b22029e
--- /dev/null
+++ b/flash2d/src/options-sphere/ADV_semilagrange_stacked.F90
@@ -0,0 +1,874 @@
+!*****************************************************************
+!
+! MODULE NAME:
+! ADV_semilagrange
+! FUNCTION:
+! perform semi-Lagrangian advection
+! CONTAINS:
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_adapt
+! FUNCTION:
+! adapt the grid according to an error estimate
+! SYNTAX:
+! CALL slm_adapt(grid, param, logical)
+! ON INPUT:
+! p_ghand: handle for the grid TYPE (grid_handle)
+! p_param: global parameter structure TYPE (global_param)
+! ON OUTPUT:
+! l_changed: flag for changed grid LOGICAL
+! CALLS:
+!
+! COMMENTS:
+!
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_initialize
+! FUNCTION:
+! initialize the advection problem
+! SYNTAX:
+! CALL slm_initialize(grid, param)
+! ON INPUT:
+! p_param: parameter data structure TYPE (global_param)
+! ON OUTPUT:
+! p_ghand: grid handling data structure TYPE (grid_handle)
+! CALLS:
+!
+! COMMENTS:
+!
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_finish
+! FUNCTION:
+! terminate slm (free dynamically alloc. memory, ...)
+! SYNTAX:
+! CALL slm_finish(grid, param)
+! ON INPUT:
+! p_ghand: grid handling data structure TYPE (grid_handle)
+! p_param: parameter data structure TYPE (global_param)
+! ON OUTPUT:
+!
+! CALLS:
+!
+! COMMENTS:
+!
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_timestepping
+! FUNCTION:
+! perform the timestepping in the slm
+! SYNTAX:
+! CALL slm_timestepping(grid, param, cmd)
+! ON INPUT:
+! p_ghand: grid handling data structure TYPE (grid_handle)
+! p_param: parameter data structure TYPE (global_param)
+! p_cmdln: command line argument struct. TYPE (cmdline)
+! ON OUTPUT:
+! p_ghand: grid handling data structure TYPE (grid_handle)
+! CALLS:
+!
+! COMMENTS:
+!
+!-----------------------------------------------------------------
+!
+! PUBLIC:
+! slm_displace, slm_update, slm_upstream
+! slm_initialize, slm_finish, slm_timestepping
+! COMMENTS:
+!
+! USES:
+! MISC_globalparam, MISC_error, FEM_handle
+! FEM_errorestimate, FEM_param
+! LIBRARIES:
+!
+! REFERENCES:
+!
+! VERSION(S):
+! 1. original version j. behrens 10/96
+! 2. several improvements/fixes j. behrens 11/96-1/97
+! 3. nodal values time depend. j. behrens 1/97
+! 4. stop_watch removed, plot j. behrens 1/97
+! (position) changed, inner
+! iteration counter added
+! 5. slm_adapt changed j. behrens 2/97
+! 6. slm_adapt changed to hide
+! grid data structures j. behrens 7/97
+! 7. control data structure j. behrens 12/97
+! 8. non-homog. rhs added j. behrens 2/98
+! 9. compliant to amatos 1.0 j. behrens 12/2000
+! 10. compliant to amatos 1.2 j. behrens 3/2002
+! 11. compliant to amatos 2.0 j. behrens 7/2003
+! 12. added visnetplot f. klaschka 12/2003
+!
+!*****************************************************************
+ MODULE ADV_semilagrange
+ USE FLASH_parameters
+ USE MISC_timing
+ USE IO_vtuplot
+#ifndef NO_NETCDF
+ USE IO_netcdfplot
+#endif
+ USE IO_utils
+ USE GRID_api
+ USE SLM_errorestimate
+ USE SLM_initial
+ USE SLM_simple
+ USE ADV_wind
+ USE ADV_rhs
+ PRIVATE
+ PUBLIC :: slm_initialize, slm_finish, slm_timestepping
+
+! USE SLM_advanced
+
+ CONTAINS
+!*****************************************************************
+ SUBROUTINE slm_adapt(p_ghand, p_param, l_changed, l_water)
+
+
+!---------- local declarations
+
+ IMPLICIT NONE
+ TYPE (grid_handle), INTENT(inout) :: p_ghand
+ TYPE (control_struct), INTENT(in) :: p_param
+ LOGICAL, INTENT(out) :: l_changed
+ LOGICAL, OPTIONAL, INTENT(in) :: l_water
+ LOGICAL :: l_switch
+ REAL (KIND = GRID_SR) :: r_errmx, &
+ r_refcrit, r_crscrit, r_fac
+ INTEGER :: i_size, &
+ i_manyc, i_manyr, i_alct, i_cnt
+ REAL (KIND = GRID_SR), DIMENSION(:), ALLOCATABLE :: r_aux1
+ INTEGER, DIMENSION(:), ALLOCATABLE :: i_aux1, i_aux2, i_aux3
+ LOGICAL :: l_ref, l_crs
+
+!---------- initialize refinement flag
+
+ l_changed= .FALSE.
+
+!---------- handle watermark switch
+
+ wat_present: IF(present(l_water)) THEN
+ l_switch= l_water
+ ELSE wat_present
+ l_switch= .TRUE.
+ END IF wat_present
+
+!---------- allocate work arrays
+
+ i_size= p_ghand%i_enumfine
+ allocate(r_aux1(i_size), i_aux1(i_size), i_aux2(i_size), &
+ i_aux3(i_size), stat=i_alct)
+ not_alloc: IF(i_alct /= 0) THEN
+ CALL grid_error(35)
+ END IF not_alloc
+
+!---------- estimate the local error
+
+ CALL slm_errorest(p_ghand, i_size, r_aux1)
+
+!---------- set coarsening/refinement criterion
+
+ r_errmx= maxval(r_aux1(1:i_size))
+ r_crscrit= r_errmx* p_param%num%r_crstolerance
+ r_refcrit= r_errmx* p_param%num%r_reftolerance
+
+!---------- get level information and set up flags for refinement/coarsening
+
+ CALL grid_getinfo(p_ghand, l_finelevel= .TRUE., i_elementlevel= i_aux1, &
+ i_elementstatus= i_aux3)
+ DO i_cnt=1,i_size
+ i_aux2(i_cnt)= 0
+ IF((i_aux1(i_cnt) > p_param%num%i_crslevel) .AND. &
+ (r_aux1(i_cnt) < r_crscrit)) i_aux2(i_cnt)= GRID_pleasecoarse
+ IF((i_aux1(i_cnt) < p_param%num%i_reflevel) .AND. &
+ (r_aux1(i_cnt) > r_refcrit)) i_aux2(i_cnt)= GRID_pleaserefine
+ END DO
+
+!---------- determine if there is enough to be done (this can be
+! switched off by l_water=.FALSE.)
+
+ IF(l_switch) THEN
+ i_manyr= count(i_aux2 == GRID_pleaserefine)
+ r_fac= real(i_manyr,GRID_SR)/ real(i_size,GRID_SR)
+ enough_ref: IF(r_fac > p_param%num%r_refwatermark) THEN
+ l_ref= .TRUE.
+ ELSE
+ l_ref= .FALSE.
+ END IF enough_ref
+
+ i_manyc= count(i_aux2 == GRID_pleasecoarse)
+ r_fac= real(i_manyc,GRID_SR)/ real(i_size,GRID_SR)
+ enough_crs: IF(r_fac > p_param%num%r_crswatermark) THEN
+ l_crs= .TRUE.
+ ELSE
+ l_crs= .FALSE.
+ END IF enough_crs
+ ELSE
+ l_ref= .TRUE.
+ l_crs= .TRUE.
+ END IF
+
+!---------- update grid flags
+
+ update: IF(l_ref .OR. l_crs) THEN
+ IF(l_ref) i_aux3= merge(i_aux2, i_aux3, i_aux2==GRID_pleaserefine)
+ IF(l_crs) i_aux3= merge(i_aux2, i_aux3, i_aux2==GRID_pleasecoarse)
+ CALL grid_putinfo(p_ghand, l_finelevel= .TRUE., i_elementstatus= i_aux3)
+ END IF update
+
+!---------- deallocate work arrays
+
+ deallocate(r_aux1, i_aux1, i_aux2, i_aux3)
+
+!---------- adapt the grid
+
+ CALL grid_adapt(p_ghand, l_changed)
+
+ RETURN
+ END SUBROUTINE slm_adapt
+
+!*****************************************************************
+ SUBROUTINE slm_diagnostics(p_ghand, p_param, p_tinfo, c_action)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+ TYPE (grid_handle), INTENT(in) :: p_ghand
+ TYPE (control_struct), INTENT(inout) :: p_param
+ TYPE (rt_info), INTENT(in) :: p_tinfo
+ CHARACTER (len=4), INTENT(in), OPTIONAL :: c_action
+ INTEGER, SAVE :: i_iodiag
+ CHARACTER (len=32) :: c_file
+ CHARACTER (len=28) :: c_tmp
+ REAL (KIND = GRID_SR), PARAMETER :: r_1o3= (1./3.)
+ INTEGER :: i_fst, i_tmp, &
+ i_size, i_alct, i_1, i_2, i_3, i_4, i_5, i_6
+ REAL (KIND = GRID_SR), SAVE :: r_rfm0, r_rsm0
+ REAL (KIND = GRID_SR) :: r_dispn, r_rfm, r_rsm, &
+ r_ts, r_calci, r_calcs, r_mxnrm, r_l2nrm, r_max, r_min, r_diffn, &
+ r_medln
+ REAL (KIND = GRID_SR), DIMENSION(:), ALLOCATABLE :: r_aux1, r_aux2, &
+ r_aux3, r_aux4
+ REAL (KIND = GRID_SR), DIMENSION(:,:), ALLOCATABLE :: r_auxx
+ INTEGER, DIMENSION(1) :: i_valind
+
+!---------- action init
+
+ present_act: IF(present(c_action)) THEN
+ action_type: IF(c_action == 'init') THEN
+
+!---------- open file for diagnostic output
+
+ i_iodiag= 9
+ i_tmp = p_param%num%i_experiment
+ write(c_tmp,*) trim(GRID_parameters%program_name), '_diag.'
+ write(c_file,1010) trim(c_tmp), i_tmp
+ c_file= adjustl(c_file)
+ open(i_iodiag, file= c_file, action= 'write', form= 'formatted', &
+ iostat= i_fst)
+ not_opened: IF(i_fst /= 0) THEN
+ CALL grid_error(36)
+ END IF not_opened
+ IF(GRID_parameters%iolog > 0) &
+ write(GRID_parameters%iolog,*) 'INFO: Filename: ', c_file, ' opened on unit: ', i_iodiag
+
+!---------- allocate workspace
+
+ i_size= p_ghand%i_nnumber
+ allocate(r_aux1(i_size), r_aux2(i_size), r_aux3(i_size), &
+ r_aux4(i_size), r_auxx(1,i_size), stat=i_alct)
+ not_alloc: IF(i_alct /= 0) THEN
+ CALL grid_error(37)
+ END IF not_alloc
+ r_aux1= 0.0; r_aux2= 0.0; r_aux3= 0.0; r_aux4= 0.0
+
+!---------- get minimum edge length
+
+ CALL grid_edgelength(p_ghand, r_min=r_medln)
+
+!---------- calculate reference values, ... extract actual calculated concentration
+
+ i_valind= (/ GRID_tracer /)
+ CALL grid_getinfo(p_ghand, i_arraypoint=i_valind, &
+ r_nodevalues= r_auxx)
+ r_aux1(:)= r_auxx(1,:)
+ DEALLOCATE(r_auxx)
+
+!---------- calculate area pieces for each node
+
+ CALL grid_nodearea(p_ghand, i_size, r_aux2)
+
+!---------- calculate analytical solution
+
+ r_ts= p_param%num%r_deltatime* float(p_tinfo%i_step)
+ CALL slm_analyticsolution(p_ghand, r_ts, i_size, r_aux3)
+
+!---------- now the integral of the concentration (mass) is
+
+ r_calci= dot_product(r_aux1, r_aux2)
+ r_rfm0 = r_calci
+
+!---------- the integral of the squared concentration ("entropy"(?)) is
+
+ r_aux4 = r_aux1* r_aux1
+ r_calcs= dot_product(r_aux4, r_aux2)
+ r_rsm0 = r_calcs
+
+!---------- the maximum-norm of the error is
+
+ r_aux4 = abs(r_aux1- r_aux3)
+ r_mxnrm= maxval(r_aux4)
+
+!---------- the l2-norm of the error is
+
+ r_aux4 = r_aux4* r_aux4
+ r_l2nrm= dot_product(r_aux4, r_aux2)
+
+!---------- maximum and minimum
+
+ r_max = maxval(r_aux1)
+ r_min = minval(r_aux1)
+
+!---------- diffusion and dispersion (not yet implemented)
+
+ r_diffn= 0.0
+ r_dispn= 0.0
+
+!---------- print it
+
+ r_rfm= r_calci/r_rfm0
+ r_rsm= r_calcs/r_rsm0
+ write(i_iodiag,1100) GRID_parameters%program_name, GRID_parameters%version, &
+ GRID_parameters%subversion, GRID_parameters%patchversion
+ i_1= p_tinfo%i_step
+ i_2= p_ghand%i_enumber
+ i_3= p_ghand%i_enumfine
+ i_4= p_ghand%i_gnumber
+ i_5= p_ghand%i_gnumfine
+ i_6= p_ghand%i_nnumber
+ write(i_iodiag,1000) i_1, i_2, i_3, i_4, i_5, i_6, r_min, r_max, &
+ r_rfm, r_rsm, r_mxnrm, r_l2nrm, r_diffn, &
+ r_dispn, r_medln
+
+!---------- deallocate workspace
+
+ deallocate(r_aux1, r_aux2, r_aux3, r_aux4)
+
+!---------- initialization done
+
+ RETURN
+
+!---------- action quit
+
+ ELSE IF(c_action == 'quit') THEN action_type
+
+!---------- close diagnostic output file
+
+ close(i_iodiag)
+ IF(GRID_parameters%iolog > 0) &
+ write(GRID_parameters%iolog,*) 'INFO: Closed file on unit: ', i_iodiag
+
+!---------- action quit done
+
+ RETURN
+ END IF action_type
+ END IF present_act
+
+!---------- action diag (default): allocate workspace
+
+ i_size= p_ghand%i_nnumber
+ allocate(r_aux1(i_size), r_aux2(i_size), r_aux3(i_size), &
+ r_aux4(i_size), r_auxx(1,i_size), stat=i_alct)
+ not_allc: IF(i_alct /= 0) THEN
+ CALL grid_error(37)
+ END IF not_allc
+ r_aux1= 0.0; r_aux2= 0.0; r_aux3= 0.0; r_aux4= 0.0
+
+!---------- get minimum edge length
+
+ CALL grid_edgelength(p_ghand, r_min=r_medln)
+
+!---------- calculate reference values, ... extract actual calculated concentration
+
+ i_valind= (/ GRID_tracer /)
+ CALL grid_getinfo(p_ghand, i_arraypoint=i_valind, &
+ r_nodevalues= r_auxx)
+ r_aux1(:)= r_auxx(1,:)
+ DEALLOCATE(r_auxx)
+
+!---------- calculate area pieces for each node
+
+ CALL grid_nodearea(p_ghand, i_size, r_aux2)
+
+!---------- calculate analytical solution
+
+ r_ts= p_param%num%r_deltatime* float(p_tinfo%i_step)
+ CALL slm_analyticsolution(p_ghand, r_ts, i_size, r_aux3)
+
+!---------- now the integral of the concentration (mass) is
+
+ r_calci= dot_product(r_aux1, r_aux2)
+
+!---------- the integral of the squared concentration ("entropy"(?)) is
+
+ r_aux4 = r_aux1* r_aux1
+ r_calcs= dot_product(r_aux4, r_aux2)
+
+!---------- the maximum-norm of the error is
+
+ r_aux4 = abs(r_aux1- r_aux3)
+ r_mxnrm= maxval(r_aux4)
+
+!---------- the l2-norm of the error is
+
+ r_aux4 = r_aux4* r_aux4
+ r_l2nrm= dot_product(r_aux4, r_aux2)
+
+!---------- maximum and minimum
+
+ r_max = maxval(r_aux1)
+ r_min = minval(r_aux1)
+
+!---------- diffusion and dispersion (not yet implemented)
+
+ r_diffn= 0.0
+ r_dispn= 0.0
+
+!---------- print it
+
+ r_rfm= r_calci/r_rfm0
+ r_rsm= r_calcs/r_rsm0
+ i_1= p_tinfo%i_step
+ i_2= p_ghand%i_enumber
+ i_3= p_ghand%i_enumfine
+ i_4= p_ghand%i_gnumber
+ i_5= p_ghand%i_gnumfine
+ i_6= p_ghand%i_nnumber
+ write(i_iodiag,1000) i_1, i_2, i_3, i_4, i_5, i_6, r_min, r_max, &
+ r_rfm, r_rsm, r_mxnrm, r_l2nrm, r_diffn, &
+ r_dispn, r_medln
+
+!---------- deallocate workspace
+
+ deallocate(r_aux1, r_aux2, r_aux3, r_aux4)
+
+ RETURN
+ 1000 FORMAT(1x, i10, 1x, i10, 1x, i10, 1x, i10, 1x, i10, 1x, i10, &
+ 1x, e15.8, 1x, e15.8, 1x, e15.8, 1x, e15.8,&
+ 1x, e15.8, 1x, e15.8, 1x, e15.8, 1x, e15.8,&
+ 1x, e15.8)
+ 1010 FORMAT(a28,i4.4)
+ 1100 FORMAT(1x,'*******************************************', &
+ '*******************************************', &
+ '*******************************************', &
+ '*******************************************', &
+ '*************************************',/ &
+ 1x,'***** PROGRAM: ',a15,174x,'*****',/ &
+ 1x,'***** VERSION: ',i2.2,'.',i2.2,'.',i2.2,181x,'*****',/ &
+ 1x,'***** Diagnostic output ',180x,'*****',/ &
+ 1x,'*******************************************', &
+ '*******************************************', &
+ '*******************************************', &
+ '*******************************************', &
+ '*************************************',/ &
+ 1x,'* timestep ',' elements ',' fine el. ',' edges ', &
+ ' fine ed. ',' nodes ',' minimum ', &
+ ' maximum ',' RFM ',' RSM ', &
+ ' max-norm ',' l2-norm ',' diffusion ', &
+ ' dispersion ',' min.edge len.*',/ &
+ 1x,'*******************************************', &
+ '*******************************************', &
+ '*******************************************', &
+ '******************************************', &
+ '**************************************')
+ END SUBROUTINE slm_diagnostics
+
+!*****************************************************************
+ SUBROUTINE slm_initialize(p_ghand, p_param)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ TYPE (grid_handle), DIMENSION(GRID_timesteps), INTENT(out) :: p_ghand
+ TYPE (control_struct), INTENT(inout) :: p_param
+
+ INTEGER :: i_steps
+ CHARACTER (len=32) :: c_file
+ CHARACTER (len=28) :: c_tmp
+ INTEGER :: i_tmp, i_cnt
+ LOGICAL :: l_refined
+ INTEGER :: i_vertnum, i_alct
+ REAL (KIND = GRID_SR), DIMENSION(:,:), POINTER :: r_vertinit
+
+!---------- decide whether a new experiment is startet or an old one is continued
+
+ new_experiment: IF(p_param%num%i_experiment <= 0) THEN
+
+!---------- reset timesteps (start with 1 in any case)
+
+ time_one: IF(p_param%num%i_frsttimestep /= 1) THEN
+ IF(GRID_parameters%iolog > 0) &
+ write(GRID_parameters%iolog,*) 'WARNING : Timestep counters reset due to new experiment'
+ i_steps= p_param%num%i_lasttimestep- p_param%num%i_frsttimestep
+ p_param%num%i_frsttimestep= 1
+ p_param%num%i_lasttimestep= p_param%num%i_frsttimestep+ i_steps
+ END IF time_one
+
+!---------- initialize variables for stacked layers
+
+ i_stcnumlayers = p_param%tst%tst_int(1,1)
+ i_stctracer = grid_registerfemvar(1,i_length=i_stcnumlayers)
+ allocate(r_stclayerpressure(i_stcnumlayers), stat=i_alct)
+ not_allc: IF(i_alct /= 0) THEN
+ CALL grid_error(c_error='[slm_initialize]: Could not allocate layers pressure array')
+ END IF not_allc
+
+
+!---------- initialize grid parameters
+
+ CALL grid_setparameter(p_ghand, i_coarselevel= p_param%num%i_crslevel, &
+ i_finelevel= p_param%num%i_reflevel)
+
+!---------- define domain, first read data from file (compiled here)
+
+ CALL grid_readdomain(i_vertnum, r_vertinit, c_readfile=p_param%io%c_domainfile)
+ CALL grid_definegeometry(i_vertnum, r_vertexarr= r_vertinit)
+
+!---------- create initial triangulation
+
+ CALL grid_createinitial(p_ghand, c_filename=p_param%io%c_triangfile)
+
+!---------- initialize grid and adapt at steep gradients
+
+ i_cnt= 0
+ l_refined= .TRUE.
+ refine_loop: DO WHILE (l_refined)
+ CALL slm_initialvalues(p_ghand(i_timeplus))
+ CALL slm_adapt(p_ghand(i_timeplus), p_param, l_refined, &
+ l_water=.FALSE.)
+ END DO refine_loop
+
+!---------- duplicate grid (old time)
+
+ CALL grid_timeduplicate(p_ghand(i_timeplus), p_ghand(i_time))
+
+!---------- initialize wind field calculation
+
+ CALL slm_windinit(p_param)
+
+!---------- if an old experiment is to be continued from stored data:
+
+ ELSE new_experiment
+
+!---------- create grid from saveset, first compile filename
+
+ i_tmp= p_param%num%i_experiment- 1
+ write(c_tmp,*) trim(GRID_parameters%program_name), '_save.'
+ write(c_file,1010) trim(c_tmp), i_tmp
+ c_file= adjustl(c_file)
+
+ CALL grid_readinitial(p_ghand, c_file)
+
+!---------- initialize wind field calculation
+
+ CALL slm_windinit(p_param)
+
+ END IF new_experiment
+
+ RETURN
+ 1010 FORMAT(a28,i4.4)
+ END SUBROUTINE slm_initialize
+
+!*****************************************************************
+ SUBROUTINE slm_finish(p_ghand, p_param)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ TYPE (grid_handle), DIMENSION(GRID_timesteps), INTENT(in) :: p_ghand
+ TYPE (control_struct), INTENT(in) :: p_param
+ CHARACTER (len=32) :: c_file
+ CHARACTER (len=28) :: c_tmp
+ INTEGER :: i_tmp
+
+!---------- open and write saveset, if required
+
+ save_req: IF(p_param%io%i_savelast /= 0) THEN
+
+ i_tmp= p_param%num%i_experiment
+ write(c_tmp,*) trim(GRID_parameters%program_name), '_save.'
+ write(c_file,1010) trim(c_tmp), i_tmp
+ c_file= adjustl(c_file)
+ CALL grid_writesaveset(c_file, p_ghand)
+
+!---------- write parameter file for next experiment
+
+ CALL io_putinputfile(p_param)
+ END IF save_req
+
+!---------- gracefully terminate wind field calculations
+
+ CALL slm_windquit
+
+ RETURN
+ 1010 FORMAT(a28,i4.4)
+ END SUBROUTINE slm_finish
+
+!*****************************************************************
+ SUBROUTINE slm_timestepping(p_ghand, p_param)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ INTEGER (KIND = GRID_SI), PARAMETER :: i_innermax=15
+ TYPE (grid_handle), DIMENSION(GRID_timesteps), &
+ INTENT(inout) :: p_ghand
+ TYPE (control_struct), INTENT(inout) :: p_param
+ INTEGER (KIND = GRID_SI) :: i_timecount
+ TYPE (sw_info) :: p_time, p_timeaux
+ LOGICAL :: l_refined
+ REAL (KIND = GRID_SR), DIMENSION(:), ALLOCATABLE :: r_tracer
+ REAL (KIND = GRID_SR), DIMENSION(:,:), ALLOCATABLE :: r_coord, r_aux
+ REAL (KIND = GRID_SR), DIMENSION(:,:), ALLOCATABLE :: r_3dtracer
+ CHARACTER (len=32) :: c_file, c_matfile
+ CHARACTER (len=28) :: c_tmp
+ INTEGER (KIND = GRID_SI) :: i_tmp, &
+ i_size, i_alct, i_tst, i_fst, i_lay, i_start
+ INTEGER (KIND = GRID_SI) :: i_iomatl=21
+ REAL (KIND = GRID_SR) :: r_modtime
+ INTEGER (KIND = GRID_SI) :: i_newlen, i_ndep
+ INTEGER (KIND = GRID_SI) :: i_layerid
+ INTEGER (KIND = GRID_SI), DIMENSION(i_stcnumlayers) :: i_valind
+
+!---------- initialize constant
+
+ i_ndep = 1_GRID_SI
+
+!---------- initialize timestep info structure
+
+ p_timestepinfo%i_step = 0
+ p_timestepinfo%i_adapit = 0
+ p_timestepinfo%l_ploted = .FALSE.
+ p_timestepinfo%l_saved = .FALSE.
+ p_timestepinfo%r_modeltime = 0.0
+
+!---------- initialize timing structure
+
+ p_time%p_tim%r_tim = 0.0
+ p_time%p_tim%r_lap = 0.0
+ p_time%p_tim%c_tim = ' '
+ p_timeaux%p_tim%r_tim= 0.0
+ p_timeaux%p_tim%r_lap= 0.0
+ p_timeaux%p_tim%c_tim= ' '
+
+!---------- initialize stop watches
+
+ CALL stop_watch_init(1,(/'total time '/),p_timeaux)
+ CALL stop_watch_init(8,(/'plotting ', 'grid duplication', &
+ 'trajectory calc.', 'right hand side ', &
+ 'grid value updt.', 'grid adaption ', &
+ 'diagnostics ', 'whole timestep '/), p_time)
+
+!---------- if diagnostics are demanded, initialize diagnostical output
+
+ IF(p_param%io%l_diagnostics) THEN
+ p_timestepinfo%i_step= 0
+ CALL slm_diagnostics(p_grid(i_timeplus), p_param, p_timestepinfo, c_action='init')
+ END IF
+
+!---------- plot initial data
+
+ i_timecount= 0
+#ifndef NO_NETCDF
+ IF(p_param%io%l_netcdf) THEN
+ CALL plot_netcdf(p_ghand(i_timeplus), i_time=i_timecount)
+ END IF
+#endif
+ IF(p_param%io%l_vtu) THEN
+ CALL generate_vtu(p_ghand(i_timeplus), i_time=i_timecount)
+ END IF
+!---------- put out initial information
+
+ CALL io_putruntimeinfo(p_ghand(i_timeplus), p_timestepinfo, p_time)
+
+!---------- timestep loop
+
+ CALL stop_watch('start',1,p_timeaux)
+ CALL stop_watch('start',8,p_time)
+ i_timecount = 0_GRID_SI
+ p_timestepinfo%r_modeltime = p_param%num%r_starttime
+ time_loop: DO WHILE (p_timestepinfo%r_modeltime < p_param%num%r_finaltime - p_param%num%r_deltatime)
+ i_timecount = i_timecount+ 1_GRID_SI
+ p_timestepinfo%i_step = i_timecount
+ p_timestepinfo%r_modeltime = p_timestepinfo%r_modeltime + p_param%num%r_deltatime
+ p_timestepinfo%i_adapit = 0_GRID_SI
+
+!---------- duplicate old grid, use it as first guess for new grid
+
+ CALL stop_watch('start',2,p_time)
+ CALL grid_timeduplicate(p_ghand(i_time), p_ghand(i_timeplus))
+ CALL stop_watch('stop ',2,p_time)
+
+!---------- adaptive (inner) loop
+
+ l_refined= .TRUE.
+ adap_loop: DO WHILE(l_refined .AND. p_timestepinfo%i_adapit < i_innermax)
+ p_timestepinfo%i_adapit= p_timestepinfo%i_adapit+ 1
+
+!---------- allocate and extract working arrays
+!---------- use amatos 1.2 functionality to calculate only new nodes
+
+ i_size= p_ghand(i_timeplus)%i_nnumber
+ allocate(r_aux(GRID_dimension,i_size), stat=i_alct)
+ not_alloc: IF(i_alct /= 0) THEN
+ CALL grid_error(38)
+ END IF not_alloc
+
+!-SLM--------- do the following SLM calculations in arrays (grid-point-wise)
+
+ CALL grid_getinfo(p_ghand(i_timeplus), r_nodecoordinates=r_aux, &
+ i_newsdepth= 1, i_nlength= i_newlen)
+ allocate(r_tracer(i_newlen), r_coord(GRID_dimension,i_newlen), stat=i_alct)
+ not_alloc0: IF(i_alct /= 0) THEN
+ CALL grid_error(38)
+ END IF not_alloc0
+ r_coord(:,1:i_newlen)= r_aux(:,1:i_newlen)
+ deallocate(r_aux)
+
+!-SLM--------- allocate 3d tracer field
+
+ allocate(r_3dtracer(i_stcnumlayers,i_newlen), stat=i_alct)
+ not_alloc1: IF(i_alct /= 0) THEN
+ CALL grid_error(38)
+ END IF not_alloc1
+
+!-SLM--------- loop over stacked layers
+
+ layer_loop: DO i_lay=1,i_stcnumlayers
+
+!-SLM--------- call the SLM step
+
+ r_modtime= p_timestepinfo%r_modeltime- p_param%num%r_deltatime
+ CALL slm_step(p_ghand, p_param, p_time, r_modtime, i_newlen, &
+ r_coord, r_tracer, i_ndep, i_lay)
+
+!-SLM--------- store layer data in 3D array
+
+ r_3dtracer(i_lay,:)= r_tracer(:)
+
+ END DO layer_loop
+
+!-SLM--------- call vertical advection step
+
+ CALL slm_verticalstep(p_ghand, p_param, p_time, r_modtime, i_newlen, &
+ r_coord, r_3dtracer, i_ndep)
+
+!-SLM--------- update grid data structure and deallocate work arrays
+!-SLM--------- change back from (grid-point)arrays to grid data structure
+
+ i_valind= (/ ((i_stctracer-1+i_lay), i_lay=1,i_stcnumlayers) /)
+ CALL grid_putinfo(p_ghand(i_timeplus), i_arraypoint= i_valind, &
+ i_newsdepth= 1, r_nodevalues= r_3dtracer)
+ deallocate(r_coord, r_tracer, r_3dtracer)
+
+!-SLM--------- adapt the grid corresponding to an error estimate
+
+ CALL stop_watch('start',6,p_time)
+ CALL slm_adapt(p_ghand(i_timeplus), p_param, l_refined)
+ CALL stop_watch('stop ',6,p_time)
+
+ END DO adap_loop
+
+!---------- diagnostics, if requested
+
+ IF(p_param%io%l_diagnostics) THEN
+ CALL stop_watch('start',7,p_time)
+ CALL slm_diagnostics(p_grid(i_timeplus), p_param, p_timestepinfo, c_action='diag')
+ CALL stop_watch('stop ',7,p_time)
+ END IF
+
+!---------- plot data (every [i_plotoffset]th timestep)
+
+ CALL stop_watch('start',1,p_time)
+ plot_step: IF(mod(i_timecount, p_param%io%i_plotoffset) == 0) THEN
+#ifndef NO_NETCDF
+ IF(p_param%io%l_netcdf) THEN
+ CALL plot_netcdf(p_ghand(i_timeplus), i_time=i_timecount)
+ p_timestepinfo%l_ploted= .TRUE.
+ END IF
+#endif
+ IF(p_param%io%l_vtu) THEN
+ CALL generate_vtu(p_ghand(i_timeplus), i_time=i_timecount)
+ p_timestepinfo%l_ploted= .TRUE.
+ END IF
+ END IF plot_step
+ CALL stop_watch('stop ',1,p_time)
+
+!---------- put a saveset to disc every ... timesteps
+
+ save_step: IF((mod(i_timecount, p_param%io%i_saveoffset) == 0) .AND. &
+ (i_timecount > 1)) THEN
+ i_tmp= p_param%num%i_experiment
+ write(c_tmp,*) trim(GRID_parameters%program_name), '_save.'
+ write(c_file,1010) trim(c_tmp), i_tmp
+ c_file= adjustl(c_file)
+ CALL grid_writesaveset(c_file,p_ghand)
+ p_timestepinfo%l_saved= .TRUE.
+ END IF save_step
+
+!---------- runtime information output
+
+ CALL stop_watch('stop ',8,p_time)
+ CALL io_putruntimeinfo(p_ghand(i_timeplus), p_timestepinfo, p_time)
+ CALL stop_watch_init(8,(/'plotting ', 'grid duplication', &
+ 'trajectory calc.', 'right hand side ', &
+ 'grid value updt.', 'grid adaption ', &
+ 'diagnostics ', 'whole timestep '/), p_time)
+ CALL stop_watch('start',8,p_time)
+
+!---------- remove obsolecent grid items
+
+ CALL grid_sweep
+
+!---------- toggle time handles for next step if gfx-proces has not exited
+
+ CALL grid_timetoggle
+
+ END DO time_loop
+ CALL stop_watch('stop ',1,p_timeaux)
+
+!---------- print total time
+
+ write(GRID_parameters%ioout,1005)
+ write(GRID_parameters%ioout,1003) p_timeaux%p_tim(1)%r_tim
+ write(GRID_parameters%ioout,1004)
+ IF(GRID_parameters%iolog > 0) THEN
+ write(GRID_parameters%iolog,1003) p_timeaux%p_tim(1)%r_tim
+ END IF
+
+!---------- terminate diagnostics
+
+ IF(p_param%io%l_diagnostics) THEN
+ CALL slm_diagnostics(p_grid(i_timeplus), p_param, p_timestepinfo, c_action='quit')
+ END IF
+
+ RETURN
+ 1000 FORMAT(1x,'***** ***** ***** ***** ***** ***** ***** ***** ***** *****',/ &
+ 1x,'***** Runtime Information Output *****',/ &
+ 1x,'***** ----- ----- ----- ----- ----- ----- ----- ----- *****')
+ 1003 FORMAT(1x,'***** Total time for timesteps ',10x,e12.4,' *****')
+ 1004 FORMAT(1x,'***** ***** ***** ***** ***** ***** ***** ***** ***** *****',/)
+ 1005 FORMAT(1x,'***** ***** ***** ***** ***** ***** ***** ***** ***** *****',/ &
+ 1x,'***** Final Information Output *****',/ &
+ 1x,'***** ----- ----- ----- ----- ----- ----- ----- ----- *****')
+ 1010 FORMAT(a28,i4.4)
+ END SUBROUTINE slm_timestepping
+
+ END MODULE ADV_semilagrange
diff --git a/flash2d/src/options-sphere/ADV_wind_stacked.f90 b/flash2d/src/options-sphere/ADV_wind_stacked.f90
new file mode 100644
index 0000000000000000000000000000000000000000..b1267a9ce7b744a5a1b7d6b86211530d9bc3bef3
--- /dev/null
+++ b/flash2d/src/options-sphere/ADV_wind_stacked.f90
@@ -0,0 +1,864 @@
+!*****************************************************************
+!
+! MODULE NAME:
+! ADV_wind
+! FUNCTION:
+! calculate the windfield for the advection problem
+! CONTAINS:
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_windfield
+! FUNCTION:
+! calculate the advecting force for simple advection
+! SYNTAX:
+! real.arr= slm_windfield(real.arr, real)
+! ON INPUT:
+! r_coord: coordinates of point real
+! r_time: time coordinate (optional) real
+! ON OUTPUT:
+! r_field: windfield real
+! CALLS:
+!
+! COMMENTS:
+!
+!-----------------------------------------------------------------
+!
+! PUBLIC:
+!
+! COMMENTS:
+!
+! USES:
+! MISC_globalparam, MISC_error
+! LIBRARIES:
+!
+! REFERENCES:
+!
+! VERSION(S):
+! 1. original version j. behrens 12/97
+! 2. bug fix concerning interp. j. behrens 2/98
+! 3. compliant to amatos 1.0 j. behrens 12/2000
+! 4. compliant to amatos 1.2 j. behrens 3/2002
+! 5. new version for course tracer transport j. behrens 01/2012
+!
+!*****************************************************************
+ MODULE ADV_wind
+ USE GRID_api
+ USE FLASH_parameters
+ PRIVATE
+ INTEGER, PARAMETER :: i_ioerr=0
+ REAL (KIND = GRID_SR) :: r_intervallen, &
+ r_scalfacx, r_scalfacy, r_scalfacz, &
+ r_readlast, r_scalinvx, r_scalinvy, r_scalinvz
+! REAL (KIND = GRID_SR), DIMENSION(:,:,:,:), ALLOCATABLE :: r_flowx
+! REAL (KIND = GRID_SR), DIMENSION(:,:,:,:), ALLOCATABLE :: r_flowy
+! REAL (KIND = GRID_SR), DIMENSION(:,:,:,:), ALLOCATABLE :: r_flowz
+! REAL (KIND = GRID_SR), DIMENSION(:), ALLOCATABLE :: r_lat
+! REAL (KIND = GRID_SR), DIMENSION(:), ALLOCATABLE :: r_lon
+! REAL (KIND = GRID_SR), DIMENSION(:), ALLOCATABLE :: r_alt
+! REAL (KIND = GRID_SR), DIMENSION(:), ALLOCATABLE :: r_times
+ REAL*4, DIMENSION(:,:,:,:), ALLOCATABLE :: r_flowx
+ REAL*4, DIMENSION(:,:,:,:), ALLOCATABLE :: r_flowy
+ REAL*4, DIMENSION(:,:,:,:), ALLOCATABLE :: r_flowz
+ REAL*4, DIMENSION(:), ALLOCATABLE :: r_lat
+ REAL*4, DIMENSION(:), ALLOCATABLE :: r_lon
+ REAL*4, DIMENSION(:), ALLOCATABLE :: r_alt
+ REAL*8, DIMENSION(:), ALLOCATABLE :: r_times
+ INTEGER (KIND = GRID_SI) :: i_lon, i_lat, i_layer, i_timesteps
+ INTEGER (KIND = GRID_SI) :: i_pointsx, i_pointsy, i_pointsz, &
+ i_timeinterval, i_intervalnum
+
+ PUBLIC :: slm_windfield, slm_windinit, slm_windquit
+ PUBLIC :: slm_vertwind
+
+ CONTAINS
+!*****************************************************************
+ FUNCTION slm_windfield(r_coord, i_layerid, r_time) RESULT (r_field)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension), INTENT(in) :: r_coord
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_layerid
+ REAL (KIND = GRID_SR), INTENT(in), OPTIONAL :: r_time
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_field
+
+ REAL (KIND = GRID_SR) :: r_tim
+
+!---------- set time
+
+ IF(present(r_time)) THEN
+ r_tim= r_time
+ ELSE
+ r_tim= 0.0_GRID_SR
+ END IF
+
+!---------- decide, if data has to be read
+
+ data_read: IF(r_readlast <= r_tim) THEN
+
+!---------- update values for next open
+
+ r_readlast = r_readlast+ r_intervallen
+ i_timeinterval= i_timeinterval+ 1
+ END IF data_read
+
+!---------- interpolate to coordinate
+
+ r_field= data_interpol(r_coord, i_layerid, i_timeinterval)
+
+ RETURN
+ END FUNCTION slm_windfield
+
+!*****************************************************************
+ FUNCTION slm_vertwind(r_coord, i_layerid) RESULT (r_wwind)
+
+!---------- declarations
+
+ IMPLICIT NONE
+
+!---------- i/o parameters
+
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension), INTENT(in) :: r_coord
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_layerid
+ REAL (KIND = GRID_SR) :: r_wwind
+
+!---------- interpolate to coordinate
+! Note that the time interval will be globally adjusted in the
+! horizontal wind field computation!
+
+ r_wwind= vertical_interpol(r_coord, i_layerid, i_timeinterval)
+
+ RETURN
+ END FUNCTION slm_vertwind
+!*****************************************************************
+ SUBROUTINE slm_windinit(p_control)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ TYPE (control_struct) :: p_control
+ INTEGER (KIND = GRID_SI) :: i_iofil= 19
+ CHARACTER (len=80) :: a_filrow
+ INTEGER (KIND = GRID_SI) :: i_iost, i_ioend, i_alct, i_cnt
+ INTEGER (KIND = GRID_SI) :: i_mon, i_day, i_year
+ REAL (KIND = GRID_SR) :: r_tim, r_start, r_t
+
+
+!---------- read current data from NetCDF file
+
+ CALL read_netcdf_currents(p_control)
+
+!---------- initialize some values
+
+ r_readlast = 0.0_GRID_SR
+ i_timeinterval= 0_GRID_SI
+ r_intervallen = 21600._GRID_SR ! we get a new wind value every 6 hours
+ r_start = 1937256._GRID_SR ! this is the hours from 1800-01-01 for 2021...
+
+!---------- compute starting time frame
+
+ i_year= p_control%tst%tst_int(1,2)
+ i_mon = p_control%tst%tst_int(1,3)
+ i_day = p_control%tst%tst_int(1,4)
+ r_tim = time_from_datecomponents(i_year,i_mon,i_day,0,0,0,r_start)
+ time_loop: DO i_cnt=1,i_timesteps-1
+ r_t = r_times(i_cnt)
+ IF(r_times(i_cnt) .LT. r_tim) THEN
+ IF(r_times(i_cnt+1) .GE. r_tim) THEN
+ i_timeinterval = i_cnt
+ exit time_loop
+ END IF
+ END IF
+ END DO time_loop
+
+! RETURN
+ END SUBROUTINE slm_windinit
+
+!*****************************************************************
+ SUBROUTINE slm_windquit
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+!---------- deallocate wind data arrays
+
+ DEALLOCATE(r_flowx, r_flowy, r_flowz)
+ DEALLOCATE(r_lon, r_lat, r_alt, r_times)
+
+! RETURN
+ END SUBROUTINE slm_windquit
+
+!*****************************************************************
+ SUBROUTINE read_netcdf_currents(p_control)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+ INCLUDE "netcdf.inc"
+
+!---------- input parameters
+ TYPE (control_struct), INTENT(in) :: p_control
+
+!---------- local variables
+ CHARACTER (len=io_fillen) :: c_uname, c_vname, c_wname
+ INTEGER (KIND = GRID_SI) :: i_alct, i_ncstat
+ INTEGER (KIND = GRID_SI) :: i_filexid, i_fileyid, i_filezid
+ INTEGER (KIND = GRID_SI) :: i_dimid, i_varid
+ INTEGER (KIND = GRID_SI) :: i_tm, i_ln, i_lt, i_ly
+! REAL (KIND = GRID_SR), DIMENSION(:,:,:,:), ALLOCATABLE :: r_auxx
+! REAL (KIND = GRID_SR), DIMENSION(:,:,:,:), ALLOCATABLE :: r_auxy
+! REAL (KIND = GRID_SR), DIMENSION(:,:,:,:), ALLOCATABLE :: r_auxz
+ REAL*4, DIMENSION(:,:,:,:), ALLOCATABLE :: r_auxx
+ REAL*4, DIMENSION(:,:,:,:), ALLOCATABLE :: r_auxy
+ REAL*4, DIMENSION(:,:,:,:), ALLOCATABLE :: r_auxz
+
+!---------- initialize file names
+ WRITE(GRID_parameters%ioout,*) &
+ '***** ***** ***** ***** INITIALIZE ***** ***** ***** *****'
+ WRITE(GRID_parameters%ioout,*) &
+ '***** reading NetCDF file header ... *****'
+
+ WRITE(c_uname,*) trim(p_control%tst%tst_char(1))
+ c_uname= adjustl(c_uname)
+ WRITE(c_vname,*) trim(p_control%tst%tst_char(2))
+ c_vname= adjustl(c_vname)
+ WRITE(c_wname,*) trim(p_control%tst%tst_char(3))
+ c_wname= adjustl(c_wname)
+
+!---------- open current files
+
+ i_ncstat= nf_open(c_uname,NF_NOWRITE,i_filexid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not open data file for u')
+ i_ncstat= nf_open(c_vname,NF_NOWRITE,i_fileyid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not open data file for v')
+ i_ncstat= nf_open(c_wname,NF_NOWRITE,i_filezid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not open data file for w')
+
+!---------- determine lon/lat/time/levels dimension sizes (grid size of currents field)
+
+ i_ncstat= nf_inq_dimid(i_filexid, 'lon', i_dimid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not identify lon dimension')
+ i_ncstat= nf_inq_dimlen(i_filexid, i_dimid, i_lon)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not read lon dimension')
+
+ i_ncstat= nf_inq_dimid(i_filexid, 'lat', i_dimid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not identify lat dimension')
+ i_ncstat= nf_inq_dimlen(i_filexid, i_dimid, i_lat)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not read lat dimension')
+
+ i_ncstat= nf_inq_dimid(i_filexid, 'time', i_dimid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not identify time dimension')
+ i_ncstat= nf_inq_dimlen(i_filexid, i_dimid, i_timesteps)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not read time dimension')
+
+ i_ncstat= nf_inq_dimid(i_filexid, 'level', i_dimid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not identify time dimension')
+ i_ncstat= nf_inq_dimlen(i_filexid, i_dimid, i_layer)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not read time dimension')
+
+!---------- make sure layer number in file is larger than required from parameters
+
+ IF(i_layer .LT. p_control%tst%tst_int(1,1)) &
+ CALL grid_error(c_error='[read_netcdf_currents]: too few layers in NetCDF file')
+
+!---------- check dimension match in v- and w-component files
+
+ i_ncstat= nf_inq_dimid(i_fileyid, 'lon', i_dimid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not identify lon dimension')
+ i_ncstat= nf_inq_dimlen(i_fileyid, i_dimid, i_ln)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not read lon dimension')
+ IF(i_lon /= i_ln) &
+ CALL grid_error(c_error='[read_netcdf_currents]: lon dimension mismatch')
+
+ i_ncstat= nf_inq_dimid(i_fileyid, 'lat', i_dimid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not identify lat dimension')
+ i_ncstat= nf_inq_dimlen(i_fileyid, i_dimid, i_lt)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not read lat dimension')
+ IF(i_lat /= i_lt) &
+ CALL grid_error(c_error='[read_netcdf_currents]: lat dimension mismatch')
+
+ i_ncstat= nf_inq_dimid(i_fileyid, 'time', i_dimid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not identify time dimension')
+ i_ncstat= nf_inq_dimlen(i_fileyid, i_dimid, i_tm)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not read time dimension')
+ IF(i_timesteps /= i_tm) &
+ CALL grid_error(c_error='[read_netcdf_currents]: time dimension mismatch')
+
+ i_ncstat= nf_inq_dimid(i_fileyid, 'level', i_dimid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not identify time dimension')
+ i_ncstat= nf_inq_dimlen(i_fileyid, i_dimid, i_ly)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not read time dimension')
+ IF(i_layer /= i_ly) &
+ CALL grid_error(c_error='[read_netcdf_currents]: layer dimension mismatch')
+
+ i_ncstat= nf_inq_dimid(i_filezid, 'lon', i_dimid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not identify lon dimension')
+ i_ncstat= nf_inq_dimlen(i_filezid, i_dimid, i_ln)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not read lon dimension')
+ IF(i_lon /= i_ln) &
+ CALL grid_error(c_error='[read_netcdf_currents]: lon dimension mismatch')
+
+ i_ncstat= nf_inq_dimid(i_filezid, 'lat', i_dimid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not identify lat dimension')
+ i_ncstat= nf_inq_dimlen(i_filezid, i_dimid, i_lt)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not read lat dimension')
+ IF(i_lat /= i_lt) &
+ CALL grid_error(c_error='[read_netcdf_currents]: lat dimension mismatch')
+
+ i_ncstat= nf_inq_dimid(i_filezid, 'time', i_dimid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not identify time dimension')
+ i_ncstat= nf_inq_dimlen(i_filezid, i_dimid, i_tm)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not read time dimension')
+ IF(i_timesteps /= i_tm) &
+ CALL grid_error(c_error='[read_netcdf_currents]: time dimension mismatch')
+
+ i_ncstat= nf_inq_dimid(i_filezid, 'level', i_dimid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not identify time dimension')
+ i_ncstat= nf_inq_dimlen(i_filezid, i_dimid, i_ly)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not read time dimension')
+ IF(i_ly .LT. p_control%tst%tst_int(1,1)) &
+ CALL grid_error(c_error='[read_netcdf_currents]: w-layer dimension mismatch')
+
+!---------- allocate latitude and longitude coordinate arrays
+
+ ALLOCATE(r_lat(i_lat), r_lon(i_lon), r_alt(i_layer), r_times(i_timesteps), stat= i_alct)
+ IF(i_alct /= 0) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not allocate lat/lon/alt/time field')
+
+!---------- read latitude and longitude coordinate values
+
+ i_ncstat= nf_inq_varid(i_filexid, 'lon', i_varid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not determine lon varid')
+ i_ncstat= nf_get_var_real(i_filexid, i_varid, r_lon)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not read lon data')
+
+ i_ncstat= nf_inq_varid(i_filexid, 'lat', i_varid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not determine lat varid')
+ i_ncstat= nf_get_var_real(i_filexid, i_varid, r_lat)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not read lat data')
+
+ i_ncstat= nf_inq_varid(i_filexid, 'level', i_varid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not determine level varid')
+ i_ncstat= nf_get_var_real(i_filexid, i_varid, r_alt)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not read level data')
+
+ i_ncstat= nf_inq_varid(i_filexid, 'time', i_varid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not determine time varid')
+ i_ncstat= nf_get_var_double(i_filexid, i_varid, r_times)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not read time data')
+
+!---------- copy globally available pressure level information
+
+ r_stclayerpressure(1:i_stcnumlayers) = r_alt(1:i_stcnumlayers)
+
+!---------- allocate current data arrays
+ WRITE(GRID_parameters%ioout,*) &
+ '***** reading NetCDF file core data ... *****'
+
+ ALLOCATE(r_flowx(i_lon, i_lat, i_layer, i_timesteps), &
+ r_flowy(i_lon, i_lat, i_layer, i_timesteps), &
+ r_flowz(i_lon, i_lat, i_layer, i_timesteps), &
+ r_auxx(i_timesteps, i_layer, i_lat, i_lon), &
+ r_auxy(i_timesteps, i_layer, i_lat, i_lon), &
+ r_auxz(i_timesteps, i_layer, i_lat, i_lon), stat= i_alct)
+ IF(i_alct /= 0) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not allocate currents fields')
+
+!---------- read x-/y-direction data of currents
+
+ i_ncstat= nf_inq_varid(i_filexid, 'uwnd', i_varid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not determine varid of uwnd')
+ i_ncstat= nf_get_var_real(i_filexid, i_varid, r_auxx)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not read uwnd data')
+
+ i_ncstat= nf_inq_varid(i_fileyid, 'vwnd', i_varid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not determine varid of vwnd')
+ i_ncstat= nf_get_var_real(i_fileyid, i_varid, r_auxy)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not read vwnd data')
+
+ i_ncstat= nf_inq_varid(i_filezid, 'omega', i_varid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not determine varid of omega')
+ i_ncstat= nf_get_var_real(i_filezid, i_varid, r_auxz)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not read omega data')
+
+!---------- reshape fields for better accessibility
+ WRITE(GRID_parameters%ioout,*) &
+ '***** reshaping core data arrays ... *****'
+
+ r_flowx = reshape(r_auxx, shape(r_flowx), order=(/4,3,2,1/))
+ r_flowy = reshape(r_auxy, shape(r_flowy), order=(/4,3,2,1/))
+ r_flowz = reshape(r_auxz, shape(r_flowz), order=(/4,3,2,1/))
+ DEALLOCATE(r_auxx, r_auxy, r_auxz)
+
+!---------- Fix mask values
+ WRITE(GRID_parameters%ioout,*) &
+ '***** cleaning core data arrays ... *****'
+
+ DO i_tm=1,i_timesteps
+ DO i_ly= 1,i_layer
+ DO i_lt= 1,i_lat
+ DO i_ln= 1,i_lon
+ IF(r_flowx(i_ln,i_lt,i_ly,i_tm) .LE. -8.99E+30) &
+ r_flowx(i_ln,i_lt,i_ly,i_tm)= 0.0
+ IF(r_flowy(i_ln,i_lt,i_ly,i_tm) .LE. -8.99E+30) &
+ r_flowy(i_ln,i_lt,i_ly,i_tm)= 0.0
+ IF(r_flowz(i_ln,i_lt,i_ly,i_tm) .LE. -8.99E+30) &
+ r_flowz(i_ln,i_lt,i_ly,i_tm)= 0.0
+ END DO
+ END DO
+ END DO
+ END DO
+
+!---------- close currents files
+
+ i_ncstat= nf_close(i_filexid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not close uwnd data file')
+ i_ncstat= nf_close(i_fileyid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not close vwnd data file')
+ i_ncstat= nf_close(i_filezid)
+ IF(i_ncstat /= NF_NOERR) &
+ CALL grid_error(c_error='[read_netcdf_currents]: could not close omega data file')
+
+ WRITE(GRID_parameters%ioout,*) &
+ '***** finished ingesting wind data *****'
+ WRITE(GRID_parameters%ioout,*) &
+ '***** ***** ***** ***** ***** ***** ***** ***** ***** *****'
+
+! RETURN
+ END SUBROUTINE read_netcdf_currents
+
+!*****************************************************************
+
+ FUNCTION data_interpol(r_cart, i_layer, i_tim) RESULT (r_field)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_cart
+ INTEGER (KIND = GRID_SI) :: i_layer
+ INTEGER (KIND = GRID_SI) :: i_tim
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_field
+
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimspherical) :: r_inter
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimspherical) :: r_coord
+ INTEGER (KIND = GRID_SI) :: i_lox, i_hix, i_loy, i_hiy, i_cnt
+ REAL (KIND = GRID_SR) :: r_dx, r_dy, r_lx, r_hx, &
+ r_ly, r_hy, r_l1, r_h1, r_l2, r_h2, r_dxi, r_dyi
+ REAL (KIND = GRID_SR) :: r_scalx, r_scaly
+ REAL (KIND = GRID_SR) :: r_deg
+ REAL (KIND = GRID_SR) :: r_earth
+
+!---------- initialize radians
+
+ r_deg= 360._GRID_SR/(r_earthrad*2.*GRID_PI)
+
+!---------- compute lat/lon coordinates from Cartesian (given) coordinates
+
+ r_coord = grid_kartgeo(r_cart)
+
+!---------- find wind box corresponding to coordinate
+
+ i_lox=1_GRID_SI
+ i_hix=2_GRID_SI
+ i_loy=1_GRID_SI
+ i_hiy=2_GRID_SI
+
+ determine_lon: DO i_cnt=1, i_lon-1
+ IF(r_lon(i_cnt) .LT. r_coord(1)) THEN
+ IF(r_lon(i_cnt +1) .GE. r_coord(1)) THEN
+ i_lox= i_cnt
+ i_hix= i_cnt+1
+ exit determine_lon
+ END IF
+ END IF
+ END DO determine_lon
+ IF(r_lon(i_lon) .LT. r_coord(1)) THEN
+ i_lox = i_lon-1
+ i_hix = i_lon
+ END IF
+
+ determine_lat: DO i_cnt=1, i_lat-1
+ IF(r_lat(i_cnt) .LT. r_coord(2)) THEN
+ IF(r_lat(i_cnt +1) .GE. r_coord(2)) THEN
+ i_loy= i_cnt
+ i_hiy= i_cnt+1
+ exit determine_lat
+ END IF
+ END IF
+ END DO determine_lat
+ IF(r_lat(i_lat) .LT. r_coord(2)) THEN
+ i_loy = i_lat-1
+ i_hiy = i_lat
+ END IF
+
+
+!---------- calculate weights for bilinear interpolation
+! NOTE: We do not consider the distance stretching with height!
+
+ r_dx= r_lon(i_hix)- r_lon(i_lox)
+ r_dy= r_lat(i_hiy)- r_lat(i_loy)
+ r_lx= r_coord(1) - r_lon(i_lox)
+ r_ly= r_coord(2) - r_lat(i_loy)
+ r_hx= r_lon(i_hix) - r_coord(1)
+ r_hy= r_lat(i_hiy) - r_coord(2)
+ r_scalx=r_deg
+ r_scaly=r_deg*cos(r_coord(2))
+
+!---------- linear interpolation in x-direction
+
+ IF(r_dx /= 0.0_GRID_SR) THEN
+ r_dxi= 1._GRID_SR/r_dx
+ r_l1= (r_hx* r_flowx(i_lox, i_loy, i_layer, i_tim)+ &
+ r_lx* r_flowx(i_hix, i_loy, i_layer, i_tim))* r_dxi
+ r_h1= (r_hx* r_flowx(i_lox, i_hiy, i_layer, i_tim)+ &
+ r_lx* r_flowx(i_hix, i_hiy, i_layer, i_tim))* r_dxi
+ r_l2= (r_hx* r_flowy(i_lox, i_loy, i_layer, i_tim)+ &
+ r_lx* r_flowy(i_hix, i_loy, i_layer, i_tim))* r_dxi
+ r_h2= (r_hx* r_flowy(i_lox, i_hiy, i_layer, i_tim)+ &
+ r_lx* r_flowy(i_hix, i_hiy, i_layer, i_tim))* r_dxi
+ ELSE
+ r_l1= r_flowx(i_lox, i_loy, i_layer, i_tim)
+ r_h1= r_flowx(i_lox, i_hiy, i_layer, i_tim)
+ r_l2= r_flowy(i_lox, i_loy, i_layer, i_tim)
+ r_h2= r_flowy(i_lox, i_hiy, i_layer, i_tim)
+ END IF
+
+!---------- linear interpolation in y-direction
+
+ IF(r_dy /= 0.0_GRID_SR) THEN
+ r_dyi= 1._GRID_SR/r_dy
+ r_inter(1)= (r_hy* r_l1+ r_ly* r_h1)* r_dyi* r_scalx
+ r_inter(2)= (r_hy* r_l2+ r_ly* r_h2)* r_dyi* r_scaly
+ ELSE
+ r_inter(1)= r_l1* r_scalx
+ r_inter(2)= r_l2* r_scaly
+ END IF
+
+!---------- calculate wind vector in kartesian coordinates
+
+ r_field= kartwind(r_inter,r_coord)
+
+ RETURN
+ END FUNCTION data_interpol
+
+!*****************************************************************
+ FUNCTION vertical_interpol(r_cart, i_layer, i_tim) RESULT (r_inter)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_cart
+ INTEGER (KIND = GRID_SI) :: i_layer
+ INTEGER (KIND = GRID_SI) :: i_tim
+ REAL (KIND = GRID_SR) :: r_inter
+
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimspherical) :: r_coord
+ INTEGER (KIND = GRID_SI) :: i_lox, i_hix, i_loy, i_hiy, i_cnt
+ REAL (KIND = GRID_SR) :: r_dx, r_dy, r_lx, r_hx, &
+ r_ly, r_hy, r_l1, r_h1, r_l2, r_h2, r_dxi, r_dyi
+ REAL (KIND = GRID_SR) :: r_scalx, r_scaly
+ REAL (KIND = GRID_SR) :: r_deg
+ REAL (KIND = GRID_SR) :: r_earth
+
+!---------- compute lat/lon coordinates from Cartesian (given) coordinates
+
+ r_coord = grid_kartgeo(r_cart)
+
+!---------- find wind box corresponding to coordinate
+
+ i_lox=1_GRID_SI
+ i_hix=2_GRID_SI
+ i_loy=1_GRID_SI
+ i_hiy=2_GRID_SI
+
+ determine_lon: DO i_cnt=1, i_lon-1
+ IF(r_lon(i_cnt) .LT. r_coord(1)) THEN
+ IF(r_lon(i_cnt +1) .GE. r_coord(1)) THEN
+ i_lox= i_cnt
+ i_hix= i_cnt+1
+ exit determine_lon
+ END IF
+ END IF
+ END DO determine_lon
+ IF(r_lon(i_lon) .LT. r_coord(1)) THEN
+ i_lox = i_lon-1
+ i_hix = i_lon
+ END IF
+
+ determine_lat: DO i_cnt=1, i_lat-1
+ IF(r_lat(i_cnt) .LT. r_coord(2)) THEN
+ IF(r_lat(i_cnt +1) .GE. r_coord(2)) THEN
+ i_loy= i_cnt
+ i_hiy= i_cnt+1
+ exit determine_lat
+ END IF
+ END IF
+ END DO determine_lat
+ IF(r_lat(i_lat) .LT. r_coord(2)) THEN
+ i_loy = i_lat-1
+ i_hiy = i_lat
+ END IF
+
+
+!---------- calculate weights for bilinear interpolation
+! NOTE: We do not consider the distance stretching with height!
+
+ r_dx= r_lon(i_hix)- r_lon(i_lox)
+ r_dy= r_lat(i_hiy)- r_lat(i_loy)
+ r_lx= r_coord(1) - r_lon(i_lox)
+ r_ly= r_coord(2) - r_lat(i_loy)
+ r_hx= r_lon(i_hix) - r_coord(1)
+ r_hy= r_lat(i_hiy) - r_coord(2)
+
+!---------- linear interpolation in x-direction
+
+ IF(r_dx /= 0.0_GRID_SR) THEN
+ r_dxi= 1._GRID_SR/r_dx
+ r_l1= (r_hx* r_flowz(i_lox, i_loy, i_layer, i_tim)+ &
+ r_lx* r_flowz(i_hix, i_loy, i_layer, i_tim))* r_dxi
+ r_h1= (r_hx* r_flowz(i_lox, i_hiy, i_layer, i_tim)+ &
+ r_lx* r_flowz(i_hix, i_hiy, i_layer, i_tim))* r_dxi
+ ELSE
+ r_l1= r_flowz(i_lox, i_loy, i_layer, i_tim)
+ r_h1= r_flowz(i_lox, i_hiy, i_layer, i_tim)
+ END IF
+
+!---------- linear interpolation in y-direction
+
+ IF(r_dy /= 0.0_GRID_SR) THEN
+ r_dyi= 1._GRID_SR/r_dy
+ r_inter= (r_hy* r_l1+ r_ly* r_h1)* r_dyi
+ ELSE
+ r_inter= r_l1
+ END IF
+
+ RETURN
+ END FUNCTION vertical_interpol
+
+!*****************************************************************
+ FUNCTION time_from_datecomponents(i_year, i_month, i_day, i_hour, i_minute, &
+ i_second, r_start, l_leap) RESULT (r_longtime)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_year
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_month
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_day
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_hour
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_minute
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_second
+ REAL (KIND = GRID_SR), INTENT(in) :: r_start
+ LOGICAL, INTENT(in), OPTIONAL :: l_leap
+ REAL (KIND = GRID_SR) :: r_longtime
+
+ INTEGER (KIND = GRID_SI), DIMENSION(12) :: i_montharray
+ LOGICAL :: l_leapyear
+ REAL (KIND = GRID_SR) :: r_accumhours
+ INTEGER (KIND = GRID_SI) :: i_cnt, i_current
+
+!---------- conversion function -- DOES NOT WORK (2022-01-25, JB)
+! reference: https://wiki.usask.ca/display/MESH/How+to+manually+calculate+NetCDF+times
+
+! r_longtime = real((1461*(i_year + 4800 + (i_month - 14)/12))/4 + &
+! (367*(i_month - 2 - 12*((i_month - 14)/12)))/12 - &
+! (3*((i_year + 4900 + (i_month - 14)/12)/100))/4 + i_day - 32075) + &
+! i_hour/24.0 + i_minute/24.0/60.0 + i_second/24.0/60.0/60.0
+
+!---------- handle optional leap year variable
+
+ IF(present(l_leap)) THEN
+ l_leapyear= l_leap
+ ELSE
+ l_leapyear= .FALSE.
+ END IF
+
+!---------- define month's days
+
+ IF(l_leapyear) THEN
+ i_montharray = (/ 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
+ ELSE
+ i_montharray = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
+ END IF
+
+!---------- compute accumulated hours from start for the months
+
+ r_accumhours = 0.
+ i_current = 0
+ DO i_cnt=1,i_month
+ r_accumhours = r_accumhours+ 24. * i_current
+ i_current = i_montharray(i_cnt)
+ END DO
+
+!---------- compute accumulated hours from days, hours, minutes and seconds
+
+ r_accumhours = r_accumhours+ 24. * (i_day-1) + i_hour + i_minute/60 + i_second/3600
+
+!---------- add offset
+
+ r_longtime = r_accumhours + r_start
+
+ RETURN
+ END FUNCTION time_from_datecomponents
+
+!*****************************************************************
+ SUBROUTINE datecomponents_from_time(r_longtime, r_start, i_year, i_month, i_day, i_hour, &
+ i_minute, i_second, l_leap)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+ REAL (KIND = GRID_SR), INTENT(in) :: r_longtime
+ REAL (KIND = GRID_SR), INTENT(in) :: r_start
+ INTEGER (KIND = GRID_SI), INTENT(inout) :: i_year
+ INTEGER (KIND = GRID_SI), INTENT(out) :: i_month
+ INTEGER (KIND = GRID_SI), INTENT(out) :: i_day
+ INTEGER (KIND = GRID_SI), INTENT(out) :: i_hour
+ INTEGER (KIND = GRID_SI), INTENT(out) :: i_minute
+ INTEGER (KIND = GRID_SI), INTENT(out) :: i_second
+ LOGICAL, OPTIONAL :: l_leap
+
+ INTEGER (KIND = GRID_SI) :: i_auxa, i_auxb, i_auxc
+
+ INTEGER (KIND = GRID_SI), DIMENSION(12) :: i_montharray
+ LOGICAL :: l_leapyear
+ REAL (KIND = GRID_SR) :: r_accumhours, r_time, r_remhours, r_hours
+ INTEGER (KIND = GRID_SI) :: i_cnt, i_current
+
+!---------- conversion function -- DOES NOT WORK (2022-01-25, JB)
+! reference: https://wiki.usask.ca/display/MESH/How+to+manually+calculate+NetCDF+times
+
+! i_auxa = 4*(floor(r_longtime) + 1401 + (((4*floor(r_longtime) + &
+! 274277)/146097)*3)/4 - 38) + 3
+! i_auxb = mod(i_auxa, 1461)/4
+! i_auxc = 5*i_auxb + 2
+!
+! i_day = (mod(i_auxc, 153))/5 + 1
+! i_month = mod((i_auxc/153 + 2), 12) + 1
+! i_year = i_auxa/1461 - 4716 + (12 + 2 - i_month)/12
+! i_second = mod(nint((r_longtime - floor(r_longtime))*24.0*60.0*60.0), 60)
+! i_minute = mod(nint((r_longtime - floor(r_longtime))*24.0*60.0 - i_second/60.0), 60)
+! i_hour = nint((r_longtime - floor(r_longtime))*24.0 - i_minute/60.0 - i_second/60.0/60.0)
+
+!---------- handle optional leap year variable
+
+ IF(present(l_leap)) THEN
+ l_leapyear= l_leap
+ ELSE
+ l_leapyear= .FALSE.
+ END IF
+
+!---------- define month's days
+
+ IF(l_leapyear) THEN
+ i_montharray = (/ 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
+ ELSE
+ i_montharray = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
+ END IF
+
+!---------- compute month
+
+ r_time = r_longtime - r_start
+ r_accumhours = 0.
+ i_current = 0
+ i_cnt = 1
+ DO WHILE ((r_accumhours .LT. r_time) .AND. (i_cnt .LE. 12))
+ i_current = i_montharray(i_cnt)
+ r_accumhours = r_accumhours+ 24. * i_current
+ i_cnt = i_cnt + 1
+ END DO
+ i_month = i_cnt
+
+!---------- from remaining hours compute days
+
+ r_remhours = r_time - (24. * sum(i_montharray(1:i_month)))
+ i_day = floor(r_remhours/24) + 1
+ r_hours = modulo(r_remhours,24._GRID_SR)
+
+!---------- now hours, minutes, seconds
+
+ i_hour = floor(r_hours)
+ r_hours = r_hours - REAL(i_hour,GRID_SR)
+
+ i_minute = floor(r_hours*60.)
+ r_hours = r_hours*60. - REAL(i_minute,GRID_SR)
+
+ i_second = floor(r_hours*60.)
+
+ END SUBROUTINE datecomponents_from_time
+
+!*****************************************************************
+ FUNCTION kartwind(r_lamphi, r_lpcoor) RESULT (r_kart)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimspherical) :: r_lamphi, r_lpcoor
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_kart
+ REAL (KIND = GRID_SR) :: r_sl, r_cl, r_sp, r_cp
+
+!---------- calculate constants
+
+ r_sl= sin(r_lpcoor(1))
+ r_cl= cos(r_lpcoor(1))
+ r_sp= sin(r_lpcoor(2))
+ r_cp= cos(r_lpcoor(2))
+
+!---------- calculate 3D kartesian components
+
+ r_kart(1)= -r_sl* r_lamphi(1)- r_sp* r_cl* r_lamphi(2)
+ r_kart(2)= r_cl* r_lamphi(1)- r_sp* r_sl* r_lamphi(2)
+ r_kart(3)= r_cp* r_lamphi(2)
+
+ RETURN
+ END FUNCTION kartwind
+
+!*****************************************************************
+ END MODULE ADV_wind
diff --git a/flash2d/src/options-sphere/FLASH_metadata_stacked.f90 b/flash2d/src/options-sphere/FLASH_metadata_stacked.f90
new file mode 100644
index 0000000000000000000000000000000000000000..daaecc2e29298d24cfbade188ef176fbfc16a998
--- /dev/null
+++ b/flash2d/src/options-sphere/FLASH_metadata_stacked.f90
@@ -0,0 +1,75 @@
+!*******************************************************************************
+!
+!> @file FLASH_metadata.f90
+!> @brief contains metadata definitions for FLASH_parameters
+!
+!*******************************************************************************
+! MODULE DESCRIPTION:
+!> @brief This is a meta data structure that contains info for the
+!> FLASH_parameters physical parameters, related to different
+!> Test cases.
+!> @author J. Behrens
+!> @version 1.0
+!> @date 5/2016
+!*****************************************************************
+ MODULE FLASH_metadata
+
+ IMPLICIT NONE
+
+!---------- character length for comparison
+
+ INTEGER, PARAMETER :: i_comparlen=12
+
+!---------- meta data for logical parameters
+
+ LOGICAL, PARAMETER :: l_logused= .FALSE.
+ INTEGER, PARAMETER :: i_lognum=1
+ CHARACTER (len=i_comparlen), DIMENSION(i_lognum) :: c_logkeywds= &
+ (/ ' ' /)
+
+!---------- meta data for integer parameters
+
+ LOGICAL, PARAMETER :: l_intused= .TRUE.
+ INTEGER, PARAMETER :: i_intnum=4
+ INTEGER, DIMENSION(i_intnum) :: i_intsizes= (/ 1, 1, 1, 1 /)
+ CHARACTER (len=i_comparlen), DIMENSION(i_intnum) :: c_intkeywds= &
+ (/ 'NUM_VERT_LAY', &
+ 'START_TIME_Y', &
+ 'START_TIME_M', &
+ 'START_TIME_D' /)
+
+!---------- meta data for character parameters
+
+ LOGICAL, PARAMETER :: l_charused= .TRUE.
+ INTEGER, PARAMETER :: i_charnum=3
+ INTEGER, PARAMETER :: i_charlength= 64
+ CHARACTER (len=i_comparlen), DIMENSION(i_charnum) :: c_charkeywds= &
+ (/ 'UWIND_FILE_N', &
+ 'VWIND_FILE_N', &
+ 'OMEGA_FILE_N' /)
+
+!---------- meta data for real parameters
+
+ LOGICAL, PARAMETER :: l_realused= .TRUE.
+ INTEGER, PARAMETER :: i_realnum=5
+ INTEGER, DIMENSION(i_realnum) :: i_realsizes= (/ 1, 1, 1, 1, 1 /)
+ CHARACTER (len=i_comparlen), DIMENSION(i_realnum) :: c_realkeywds= &
+ (/ 'SEA_LEVEL_PR', &
+ 'SEA_LEVEL_TE', &
+ 'TEMP_LAPSE_R', &
+ 'IDEAL_GAS_CO', &
+ 'DRY_AIR_MOLA' /)
+
+!---------- convenience pointers
+
+ INTEGER, PARAMETER :: PARAM_SL_PRESS= 1
+ INTEGER, PARAMETER :: PARAM_SL_TMPTR= 2
+ INTEGER, PARAMETER :: PARAM_TEMPLAPS= 3
+ INTEGER, PARAMETER :: PARAM_IDEALGAS= 4
+ INTEGER, PARAMETER :: PARAM_DRYAIRML= 5
+ INTEGER, PARAMETER :: PARAM_UWINDFIL= 1
+ INTEGER, PARAMETER :: PARAM_VWINDFIL= 2
+ INTEGER, PARAMETER :: PARAM_OWINDFIL= 3
+ INTEGER, PARAMETER :: PARAM_NUMLAYER= 1
+
+ END MODULE FLASH_metadata
diff --git a/flash2d/src/options-sphere/FLASH_parameters_stacked.f90 b/flash2d/src/options-sphere/FLASH_parameters_stacked.f90
new file mode 100644
index 0000000000000000000000000000000000000000..9dd8c8ccad383f2486c5317674287b7f023c0e71
--- /dev/null
+++ b/flash2d/src/options-sphere/FLASH_parameters_stacked.f90
@@ -0,0 +1,126 @@
+!*****************************************************************
+!
+! MODULE NAME:
+! FLASH_parameters
+! FUNCTION:
+! defines global control structure
+! CONTAINS:
+!
+! PUBLIC:
+! all
+! COMMENTS:
+!
+! LIBRARIES:
+!
+! REFERENCES:
+!
+! VERSION(S):
+! 1. original version j. behrens 12/2000
+!
+!*****************************************************************
+ MODULE FLASH_parameters
+ USE GRID_api
+
+ IMPLICIT NONE
+ INTEGER, PARAMETER :: io_fillen=128
+ INTEGER, PARAMETER :: i_redirout=8
+ INTEGER, PARAMETER :: i_redirlog=7
+ INTEGER, PARAMETER :: i_comparlen=14
+
+!---------- structure for the command line
+
+ TYPE cmdline_param
+ SEQUENCE
+ CHARACTER (len=io_fillen) :: c_infile ! input file name
+ CHARACTER (len=io_fillen) :: c_directory ! input file name
+ LOGICAL :: l_output ! redirect std output
+ LOGICAL :: l_logging ! enable logging (verbose)
+ END TYPE cmdline_param
+
+!---------- structure for the i/o behaviour
+
+ TYPE io_param
+ SEQUENCE
+ LOGICAL :: l_diagnostics ! switch on diagnostics
+ LOGICAL :: l_vtu ! switch on vtu output
+ LOGICAL :: l_netcdf ! switch on NetCDF output
+ INTEGER :: i_plotoffset ! timesteps between plots
+ INTEGER :: i_saveoffset ! timesteps between savesets
+ INTEGER :: i_savelast ! indicator for last step saving
+ CHARACTER (len=io_fillen) :: c_domainfile ! file with definitions for domain
+ CHARACTER (len=io_fillen) :: c_triangfile ! file with initial triangulation
+ END TYPE io_param
+
+!---------- structure for global physical and steering parameters
+
+ TYPE num_param
+ SEQUENCE
+ REAL :: r_deltatime ! timestep length [s]
+ REAL :: r_reftolerance ! tolerance for refinement
+ REAL :: r_crstolerance ! tolerance for coarsening
+ REAL :: r_refwatermark ! watermark for refinement
+ REAL :: r_crswatermark ! watermark for coarsening
+ REAL :: r_starttime ! first time (overwrites i_frsttimestep)
+ REAL :: r_finaltime ! last time (overwrites i_lasttimestep)
+ INTEGER :: i_experiment ! current experiment identification
+ INTEGER :: i_crslevel ! coarsest requested level
+ INTEGER :: i_reflevel ! finest requested level
+ INTEGER :: i_frsttimestep ! first timestep of experiment
+ INTEGER :: i_lasttimestep ! last timestep of experiment
+ INTEGER :: i_adviterations ! iterations in trajectory estimation
+ END TYPE num_param
+
+!---------- structure for global physical and steering parameters
+
+ TYPE test_param
+ INTEGER :: i_lognum
+ INTEGER :: i_intnum
+ INTEGER :: i_realnum
+ INTEGER :: i_charnum
+ INTEGER, DIMENSION(:), POINTER :: i_intsizes
+ INTEGER, DIMENSION(:), POINTER :: i_realsizes
+ CHARACTER (len=i_comparlen), DIMENSION(:), POINTER :: c_logkeywds
+ CHARACTER (len=i_comparlen), DIMENSION(:), POINTER :: c_intkeywds
+ CHARACTER (len=i_comparlen), DIMENSION(:), POINTER :: c_realkeywds
+ CHARACTER (len=i_comparlen), DIMENSION(:), POINTER :: c_charkeywds
+ LOGICAL, DIMENSION(:), POINTER :: tst_log
+ CHARACTER (len=io_fillen), DIMENSION(:), POINTER :: tst_char
+ INTEGER, DIMENSION(:,:), POINTER :: tst_int
+ REAL, DIMENSION(:,:), POINTER :: tst_real
+ END TYPE test_param
+
+!---------- global control structure
+
+ TYPE control_struct
+ TYPE (num_param) :: num
+ TYPE (cmdline_param) :: cmd
+ TYPE (io_param) :: io
+ TYPE (test_param) :: tst
+ END TYPE control_struct
+ TYPE (control_struct) :: p_contr
+
+!---------- structure for runtime information
+
+ TYPE rt_info
+ REAL :: r_modeltime
+ INTEGER :: i_step
+ INTEGER :: i_adapit
+ LOGICAL :: l_saved
+ LOGICAL :: l_ploted
+ END TYPE rt_info
+ TYPE (rt_info) :: p_timestepinfo
+
+!---------- convenience pointers to stacked layers
+
+ INTEGER (KIND = GRID_SI) :: i_stctracer
+ INTEGER (KIND = GRID_SI) :: i_stcnumlayers
+
+!---------- array for storing stacked layer pressure levels
+
+ REAL (KIND = GRID_SR), DIMENSION(:), POINTER :: r_stclayerpressure
+
+!---------- array for storing stacked layer pressure levels
+
+ REAL (KIND = GRID_SR), PARAMETER :: r_earthrad=6378134._GRID_SR
+
+ END MODULE FLASH_parameters
diff --git a/flash2d/src/options-sphere/IO_vtu_stacked.F90 b/flash2d/src/options-sphere/IO_vtu_stacked.F90
new file mode 100644
index 0000000000000000000000000000000000000000..f2d0481aac47e8267a7dee53a174feaebea8a231
--- /dev/null
+++ b/flash2d/src/options-sphere/IO_vtu_stacked.F90
@@ -0,0 +1,585 @@
+!-----------------------------------------------------------------------
+!> \file IO_vtu.F90
+!> \brief Contains the module IO_vtu
+!-----------------------------------------------------------------------
+!> The MODULE IO_vtu create a vtu file, unstructered grid, readable by
+!! paraview
+!
+!> Actually, amatos' structure denies a unified handling for
+!! 2d, 3d and spherical case. Therefore the corresponding definition
+!! has to be uncommented
+!-----------------------------------------------------------------------
+
+!#define VTU_OUTPUT2D
+!#define VTU_OUTPUT3D
+!#define VTU_OUTPUTSPH
+#define VTU_OUTPUTSTC
+
+MODULE IO_vtu
+
+USE GRID_api
+
+PUBLIC :: plot_vtu, t_vtu_data
+
+#if defined VTU_OUTPUT2D || defined VTU_OUTPUTSPH
+ INTEGER(KIND=GRID_SI), PARAMETER :: i_nodespercell = GRID_elementnodes
+ INTEGER(KIND=GRID_SI), PARAMETER :: i_nodesperface = GRID_edgenodes
+
+ ! This defines a triangle in VTK format
+ INTEGER(KIND=GRID_SI), PARAMETER :: i_vtucelltype = 5
+#elif defined VTU_OUTPUT3D
+ INTEGER(KIND=GRID_SI), PARAMETER :: i_nodespercell = GRID_tetranodes
+ INTEGER(KIND=GRID_SI), PARAMETER :: i_nodesperface = GRID_elementnodes
+
+ ! This defines a tetrahedron in VTK format
+ INTEGER(KIND=GRID_SI), PARAMETER :: i_vtucelltype = 10
+#elif defined VTU_OUTPUTSTC
+ INTEGER(KIND=GRID_SI), PARAMETER :: i_nodespercell = GRID_elementnodes*2
+ INTEGER(KIND=GRID_SI), PARAMETER :: i_nodesperface = GRID_edgenodes
+
+ ! This defines a wedge/prism in VTK format
+ INTEGER(KIND=GRID_SI), PARAMETER :: i_vtucelltype = 13
+#endif
+
+
+!> Structure for variable information
+TYPE t_vtu_data
+ CHARACTER(LEN=32) :: c_name !< VTU variable name
+ INTEGER(KIND=GRID_SI) :: i_size !< data dimension
+ REAL(KIND=GRID_SR),POINTER, DIMENSION(:,:) :: p_vdata !< pointer to data
+END TYPE t_vtu_data
+
+!> Very small values can cause errors when read in paraview.
+!! This is the threshold for filtering.
+REAL(KIND = GRID_SR), PARAMETER :: r_vtueps = 1e-12
+CONTAINS
+
+!---------------------------------------------------------------------
+!> \brief Writes a VTU file
+!> Use this routine to plot continuous or discontinuous data.
+!> \param[in] p_mesh The mesh
+!> \param[in] c_filename The name of the file to write to
+!> \param[in] i_nodedata Size of the array p_nodedata
+!> \param[in] p_nodedata Array of type t_vtu_data for node
+!! data. For each node one entry has to exist in
+!! p_vdata.
+!> \param[in] i_celldata Size of the array p_celldata
+!> \param[in] p_celldata Array of type t_vtu_data for cell
+!! data. For each tetrahedron one entry has to
+!! exist in p_vdata.
+#ifdef VTU_OUTPUT2D
+!> \param[in] r_zcoordinate The z-coordinate values
+!> \param[in] i_zcoordinate Use the i_zcoordinate component of p_nodedata
+!! also as z-coordinate
+#endif
+!> \param[in] l_continuous Write continuous node data.
+!! This needs less space than discontinuous.
+!! default: .true.
+!> \param[in] l_grid_info Write grid numbering information
+!! default: .false.
+!---------------------------------------------------------------------
+SUBROUTINE plot_vtu(p_mesh, c_filename, &
+ i_nodedata, p_nodedata, i_celldata, p_celldata, &
+#ifdef VTU_OUTPUT2D
+ r_zcoordinate, i_zcoordinate, &
+#endif
+#ifdef VTU_OUTPUTSTC
+ i_layers, r_height, &
+#endif
+ l_continuous, l_grid_info)
+ IMPLICIT NONE
+
+ TYPE (grid_handle), INTENT(in) :: p_mesh
+ CHARACTER (len=*), INTENT(IN) :: c_filename
+ INTEGER(KIND=GRID_SI), INTENT(IN), OPTIONAL :: i_nodedata
+ TYPE(t_vtu_data), INTENT(IN), DIMENSION(:), OPTIONAL &
+ :: p_nodedata
+ INTEGER(KIND=GRID_SI), INTENT(IN), OPTIONAL :: i_celldata
+ TYPE(t_vtu_data), INTENT(IN), DIMENSION(:), OPTIONAL &
+ :: p_celldata
+#ifdef VTU_OUTPUT2D
+ REAL(KIND=GRID_SR), POINTER, DIMENSION(:,:), OPTIONAL &
+ :: r_zcoordinate
+ INTEGER(KIND=GRID_SI), OPTIONAL, INTENT(IN) :: i_zcoordinate
+#endif
+#ifdef VTU_OUTPUTSTC
+ REAL(KIND=GRID_SR), POINTER, DIMENSION(:), OPTIONAL &
+ :: r_height
+ INTEGER(KIND=GRID_SI), OPTIONAL, INTENT(IN) :: i_layers
+#endif
+ LOGICAL, INTENT(IN), OPTIONAL :: l_grid_info
+ LOGICAL, INTENT(IN), OPTIONAL :: l_continuous
+
+ INTEGER( KIND = GRID_SI) :: i_alct
+ INTEGER( KIND = GRID_SI) :: i_cnt
+ INTEGER( KIND = GRID_SI) :: i_fst
+ INTEGER( KIND = GRID_SI) :: i_ncnt
+ INTEGER( KIND = GRID_SI) :: i_lay
+ INTEGER( KIND = GRID_SI), PARAMETER :: i_fhandle = 89
+
+ ! this variables result from different naming conventions in two and
+ ! three dimensions
+ INTEGER( KIND = GRID_SI) :: i_numberofcells
+ INTEGER( KIND = GRID_SI) :: i_numberoffaces
+ INTEGER( KIND = GRID_SI) :: i_numberofpoints
+
+ INTEGER (KIND = GRID_SI), DIMENSION(:,:), ALLOCATABLE :: i_cellnodes
+ REAL (KIND = GRID_SR), DIMENSION(:,:), ALLOCATABLE :: r_nodecoor
+ INTEGER (KIND = GRID_SI), DIMENSION(:), ALLOCATABLE :: i_nodes
+
+ LOGICAL :: l_write_grid_info
+ LOGICAL :: l_continuous_data
+
+
+#ifdef VTU_OUTPUT2D
+ IF(PRESENT(i_zcoordinate) .AND. PRESENT(r_zcoordinate)) THEN
+ CALL GRID_error(c_error='[plot_vtu] z-coordinate is defined twice')
+ END IF
+#endif
+#ifdef VTU_OUTPUTSTC
+ IF(PRESENT(i_layers)) THEN
+ i_lay = i_layers
+ ELSE
+ i_lay = 2
+ END IF
+ IF (i_lay .LT. 2) THEN
+ CALL GRID_error(c_error='[plot_vtu] in stacked mode, you cannot plot less than one layer')
+ END IF
+ IF(.NOT. PRESENT(r_height)) THEN
+ CALL GRID_error(c_error='[plot_vtu] in stacked mode, a height array is needed')
+ END IF
+#endif
+
+ ! set the value of l_continuous. The default is .true.
+ IF(PRESENT(l_continuous)) THEN
+ l_continuous_data = l_continuous
+ ELSE
+ l_continuous_data = .TRUE.
+ END IF
+
+ ! set the value of l_grid_info. The default is .false..
+ IF(PRESENT(l_grid_info)) THEN
+ l_write_grid_info = l_grid_info
+ ELSE
+ l_write_grid_info = .FALSE.
+ END IF
+
+ ! Set dimension depending types
+#if defined VTU_OUTPUT2D || defined VTU_OUTPUTSPH
+ i_numberofcells = p_mesh%i_enumfine
+ i_numberoffaces = p_mesh%i_gnumfine
+#elif defined VTU_OUTPUT3D
+ i_numberofcells = p_mesh%i_tnumfine
+ i_numberoffaces = p_mesh%i_enumfine
+#elif defined VTU_OUTPUTSTC
+ i_numberofcells = p_mesh%i_enumfine* (i_lay-1)
+ i_numberoffaces = p_mesh%i_gnumfine* i_lay + p_mesh%i_nnumber* (i_lay-1)
+#endif
+
+ IF(l_continuous_data) THEN
+ i_numberofpoints = p_mesh%i_nnumber
+ ELSE
+ i_numberofpoints = i_numberofcells * i_nodespercell
+ END IF
+
+ ! We only need the node coordinates and the nodes per tetra
+ ! as additional information.
+ ALLOCATE(r_nodecoor(GRID_dimension, p_mesh%i_nnumber), &
+ i_cellnodes(i_nodespercell, i_numberofcells), &
+ stat = i_alct)
+
+ IF(i_alct /= 0) THEN
+ CALL grid_error(c_error='[write_vtu]: could not allocate data arrays')
+ END IF
+
+ CALL GRID_getinfo(p_mesh, r_nodecoordinates = r_nodecoor, &
+#if defined VTU_OUTPUT2D || defined VTU_OUTPUTSPH
+ i_elementnodes = i_cellnodes)
+#elif defined VTU_OUTPUT3D
+ i_tetranodes = i_cellnodes)
+#endif
+
+ ! open the file and write header information
+ OPEN(i_fhandle, file = c_filename, iostat= i_fst)
+ IF(i_fst /= 0) THEN
+ RETURN
+ END IF
+
+ WRITE(i_fhandle, "(A)") '<?xml version="1.0"?>'
+ WRITE(i_fhandle, "(A)") '<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">'
+ WRITE(i_fhandle, *) '<UnstructuredGrid>'
+ WRITE(i_fhandle, "(A,I8,A,I8,A)") '<Piece NumberOfPoints="', &
+ i_numberofpoints, &
+ '" NumberOfCells="', &
+ i_numberofcells, &
+ '">'
+
+ ! write the node coordinates to VTU file
+ WRITE(i_fhandle, *) '<Points>'
+ WRITE(i_fhandle, "(A,I1,A)") '<DataArray type="Float32" NumberOfComponents="', &
+ 3, '" format="ascii">'
+
+ ! Distunguish between continuous and discontinuous case
+ IF(l_continuous_data) THEN
+ !> The continuous case is represented by amatos' grid topology.
+ DO i_cnt = 1, p_mesh%i_nnumber
+#if defined VTU_OUTPUT3D || defined VTU_OUTPUTSPH
+ WRITE(i_fhandle, *) r_nodecoor(:, i_cnt)
+#elif defined VTU_OUTPUT2D
+ IF(PRESENT(r_zcoordinate)) THEN
+ WRITE(i_fhandle, *) r_nodecoor(:, i_cnt), r_zcoordinate(1, i_cnt)
+ ELSE IF(PRESENT(i_zcoordinate)) THEN
+ WRITE(i_fhandle, *) r_nodecoor(:, i_cnt), &
+ p_nodedata(i_zcoordinate)%p_vdata(1, i_cnt)
+ ELSE
+ WRITE(i_fhandle, *) r_nodecoor(:, i_cnt), 0.0
+ END IF
+#endif
+ END DO
+
+ ELSE
+ ! Write the node coordinates. In this discontinuous case one node in
+ ! paraview belongs only to one cell.
+ DO i_cnt = 1, i_numberofcells
+ DO i_ncnt = 1, i_nodespercell
+#if defined VTU_OUTPUT2D || defined VTU_OUTPUTSPH
+ WRITE(i_fhandle, *) r_nodecoor(:, i_cellnodes(i_ncnt, i_cnt))
+#elif defined VTU_OUTPUT2D
+ IF(PRESENT(r_zcoordinate)) THEN
+ WRITE(i_fhandle, *) r_nodecoor(:, i_cellnodes(i_ncnt, i_cnt)), &
+ r_zcoordinate(1, i_cellnodes(i_ncnt, i_cnt))
+ ELSE IF(PRESENT(i_zcoordinate)) THEN
+ WRITE(i_fhandle, *) r_nodecoor(:, i_cellnodes(i_ncnt, i_cnt)), &
+ p_nodedata(i_zcoordinate)%p_vdata(1, i_cellnodes(i_ncnt, i_cnt))
+ ELSE
+ WRITE(i_fhandle, *) r_nodecoor(:, i_cellnodes(i_ncnt, i_cnt)), 0.0
+ END IF
+#endif
+ END DO
+ END DO
+ END IF
+
+ WRITE(i_fhandle, *) '</DataArray>'
+ WRITE(i_fhandle, *) '</Points>'
+
+ ! Write the node connectivity
+ WRITE(i_fhandle, *) '<Cells>'
+ WRITE(i_fhandle, *) '<DataArray type="Int32" Name="connectivity" format="ascii">'
+
+ ! check wether continuous data is written
+ IF(l_continuous_data) THEN
+ ! The indexing in paraview starts with 0. Therefore we subtract one from
+ ! the node connectivity information.
+#if defined VTU_OUTPUT2D || defined VTU_OUTPUTSPH
+ DO i_cnt = 1, p_mesh%i_enumfine
+#elif defined VTU_OUTPUT3D
+ DO i_cnt = 1, p_mesh%i_tnumfine
+#endif
+ WRITE(i_fhandle, *) i_cellnodes(:, i_cnt) - 1
+ END DO
+ ELSE
+ ! Write element nodes (indexing starts by 0!)
+ ! In the discontinuous case each tetrahedron has its own nodes
+ ALLOCATE(i_nodes(i_nodespercell), stat = i_alct)
+
+ IF(i_alct /= 0) THEN
+ CALL grid_error(c_error='[write_vtu]: could not allocate node array for VTU output')
+ END IF
+
+ DO i_cnt = 1, i_numberofcells * i_nodespercell, i_nodespercell
+
+ DO i_ncnt = 1, i_nodespercell
+ i_nodes(i_ncnt) = i_ncnt + i_cnt - 2
+ END DO
+
+ WRITE(i_fhandle, *) i_nodes
+ END DO
+
+ DEALLOCATE(i_nodes)
+ END IF
+
+ WRITE(i_fhandle, *) '</DataArray>'
+
+ ! write the cell type. Tetrahedra are represented by 10 and triangles by 5
+ WRITE(i_fhandle, *) '<DataArray type="UInt8" Name="types" format="ascii">'
+ DO i_cnt = 1, i_numberofcells
+ WRITE(i_fhandle, *) i_vtucelltype
+ END DO
+ WRITE(i_fhandle, *) '</DataArray>'
+ WRITE(i_fhandle, *) '<DataArray type="Int32" Name="offsets" format="ascii">'
+ DO i_cnt = 1, i_numberofcells
+ WRITE(i_fhandle, *) i_cnt * i_nodespercell
+ END DO
+ WRITE(i_fhandle, *) '</DataArray>'
+ WRITE(i_fhandle, *) '</Cells>'
+
+ ! Write the point data / node data
+ ! If the grid information is requested save the node numbering of amatos.
+ ! Since every node number in paraview is decreased by one, this option is
+ ! in continuous plotting not too useful.
+ WRITE(i_fhandle, *) '<PointData>'
+
+ IF(l_write_grid_info) THEN
+ WRITE(i_fhandle, *) '<DataArray type="Int32" Name="Nodenumber" format="ascii">'
+
+ DO i_cnt = 1, p_mesh%i_nnumber
+ WRITE(i_fhandle, *) i_cnt
+ END DO
+
+ WRITE(i_fhandle, *) '</DataArray>'
+ END IF
+
+ ! When variables for the node data are present, save them.
+ IF(PRESENT(i_nodedata) .AND. PRESENT(p_nodedata)) THEN
+ DO i_cnt = 1, i_nodedata
+ CALL write_vtu_data(i_fhandle, p_nodedata(i_cnt))
+ END DO
+ END IF
+
+ WRITE(i_fhandle, *) '</PointData>'
+
+
+
+ ! write the cell data.
+ WRITE(i_fhandle, *) '<CellData>'
+
+ ! Write the element numbers when requested.
+ IF(l_write_grid_info) THEN
+ WRITE(i_fhandle, *) '<DataArray type="Int32" Name="Elementnumber" format="ascii">'
+
+ DO i_cnt = 1, i_numberofcells
+ WRITE(i_fhandle, *) i_cnt
+ END DO
+
+ WRITE(i_fhandle, *) '</DataArray>'
+ END IF
+
+ ! Write the variable belonging to the tetradra.
+ IF(PRESENT(i_celldata) .AND. PRESENT(p_celldata)) THEN
+ DO i_cnt = 1, i_celldata
+ CALL write_vtu_data(i_fhandle, p_celldata(i_cnt))
+ END DO
+ END IF
+
+ ! write the end of cell data and the footer.
+ WRITE(i_fhandle, *) '</CellData>'
+ WRITE(i_fhandle, *) '</Piece>'
+ WRITE(i_fhandle, *) '</UnstructuredGrid>'
+ WRITE(i_fhandle, *) '</VTKFile>'
+
+ ! tidy up
+ CLOSE(i_fhandle)
+ DEALLOCATE(r_nodecoor, i_cellnodes)
+END SUBROUTINE plot_vtu
+
+
+#ifdef VTU_OUTPUT3D
+!---------------------------------------------------------------------
+!> Use this routine to plot the mesh consisting of face and
+!! continuous data on it.
+!> \param p_mesh - the mesh
+!> \param c_filename - the name of the file to write to
+!> \param i_nodedata - (optional) size of the array p_nodedata
+!> \param p_nodedata - (optional) array of type t_vtu_data for node
+!! data. For each node one entry has to exist in
+!! p_vdata.
+!> \param i_celldata - (optional) size of the array p_celldata
+!> \param p_nodedata - (optional) array of type t_vtu_data for cell
+!! data. For each face one entry has to
+!! exist in p_vdata.
+!> \param l_grid_info - (optional) Write grid numbering information
+!! default: .false.
+!---------------------------------------------------------------------
+SUBROUTINE plot_vtu_elements(p_mesh, c_filename, i_nodedata, &
+ p_nodedata, i_celldata, p_celldata, &
+ l_grid_info)
+ IMPLICIT NONE
+
+ TYPE (grid_handle), INTENT(in) :: p_mesh
+ CHARACTER (len=*), INTENT(IN) :: c_filename
+ INTEGER(KIND=GRID_SI), INTENT(IN), OPTIONAL :: i_nodedata
+ TYPE(t_vtu_data), INTENT(IN), DIMENSION(:), OPTIONAL &
+ :: p_nodedata
+ INTEGER(KIND=GRID_SI), INTENT(IN), OPTIONAL :: i_celldata
+ TYPE(t_vtu_data), INTENT(IN), DIMENSION(:), OPTIONAL &
+ :: p_celldata
+
+ LOGICAL, INTENT(IN), OPTIONAL :: l_grid_info
+
+ LOGICAL :: l_write_grid_info
+
+ INTEGER( KIND = GRID_SI) :: i_alct, i_cnt, i_fst
+ INTEGER( KIND = GRID_SI) :: i_fhandle = 89
+
+ INTEGER (KIND = GRID_SI), DIMENSION(:,:), ALLOCATABLE :: i_enodes
+ REAL (KIND = GRID_SR), DIMENSION(:,:), ALLOCATABLE :: r_nodecoor
+
+ !> set the optional data
+ IF(PRESENT(l_grid_info)) THEN
+ l_write_grid_info = l_grid_info
+ ELSE
+ l_write_grid_info = .FALSE.
+ END IF
+
+ ALLOCATE(r_nodecoor(GRID_dimension, p_mesh%i_nnumber), &
+ i_enodes(GRID_elementnodes, p_mesh%i_enumfine), &
+ stat = i_alct)
+
+ IF(i_alct /= 0) THEN
+ CALL grid_error(c_error='[write_vtu]: could not allocate data arrays')
+ END IF
+
+ CALL GRID_getinfo(p_mesh, r_nodecoordinates = r_nodecoor, &
+ i_elementnodes = i_enodes)
+
+
+ OPEN(i_fhandle, file = c_filename, iostat= i_fst)
+ IF(i_fst /= 0) THEN
+ RETURN
+ END IF
+
+ !> Write the VTK header
+ WRITE(i_fhandle, "(A)") '<?xml version="1.0"?>'
+ WRITE(i_fhandle, "(A)") '<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">'
+ WRITE(i_fhandle, *) '<UnstructuredGrid>'
+ WRITE(i_fhandle, "(A,I8,A,I8,A)") '<Piece NumberOfPoints="', &
+ p_mesh%i_nnumber, &
+ '" NumberOfCells="', &
+ p_mesh%i_enumfine, &
+ '">'
+
+ !> Save the coordinates of all nodes
+ WRITE(i_fhandle, *) '<Points>'
+ WRITE(i_fhandle, "(A,I1,A)") '<DataArray type="Float32" NumberOfComponents="', &
+ 3 , '" format="ascii">'
+ DO i_cnt = 1, p_mesh%i_nnumber
+ WRITE(i_fhandle, *) r_nodecoor(:, i_cnt)
+ END DO
+ WRITE(i_fhandle, *) '</DataArray>'
+ WRITE(i_fhandle, *) '</Points>'
+
+ !> Write the cell (element) information
+ WRITE(i_fhandle, *) '<Cells>'
+
+ !> Save the node connectivity
+ WRITE(i_fhandle, *) '<DataArray type="Int32" Name="connectivity" format="ascii">'
+
+
+ !> write element nodes. Write connectivity (indexing starts by 0!)
+ DO i_cnt = 1, p_mesh%i_enumfine
+ WRITE(i_fhandle, *) i_enodes(:, i_cnt) - 1
+ END DO
+ WRITE(i_fhandle, *) '</DataArray>'
+
+ !> Write the cell type (triangle = 5, tetrahedron = 10).
+ WRITE(i_fhandle, *) '<DataArray type="UInt8" Name="types" format="ascii">'
+ DO i_cnt = 1, p_mesh%i_enumfine
+ WRITE(i_fhandle, *) '5'
+ END DO
+ WRITE(i_fhandle, *) '</DataArray>'
+ WRITE(i_fhandle, *) '<DataArray type="Int32" Name="offsets" format="ascii">'
+ DO i_cnt = 1, p_mesh%i_enumfine
+ WRITE(i_fhandle, *) i_cnt * GRID_elementnodes
+ END DO
+ WRITE(i_fhandle, *) '</DataArray>'
+ WRITE(i_fhandle, *) '</Cells>'
+
+ !> write point data
+ WRITE(i_fhandle, *) '<PointData>'
+
+ !> Save the node numbers, when requested.
+ IF(l_write_grid_info) THEN
+ WRITE(i_fhandle, *) '<DataArray type="Int32" Name="Nodenumber" format="ascii">'
+
+ DO i_cnt = 1, p_mesh%i_nnumber
+ WRITE(i_fhandle, *) i_cnt
+ END DO
+
+ WRITE(i_fhandle, *) '</DataArray>'
+ END IF
+
+ !> Write the nodal variables, when present.
+ IF(PRESENT(i_nodedata) .AND. PRESENT(p_nodedata)) THEN
+ DO i_cnt = 1, i_nodedata
+ CALL write_vtu_data(i_fhandle, p_nodedata(i_cnt))
+ END DO
+ END IF
+
+ WRITE(i_fhandle, *) '</PointData>'
+
+
+
+ !> write cell data
+ WRITE(i_fhandle, *) '<CellData>'
+
+ !> save the element numbers, when requested.
+ IF(l_write_grid_info) THEN
+ WRITE(i_fhandle, *) '<DataArray type="Int32" Name="Elementnumber" format="ascii">'
+
+ DO i_cnt = 1, p_mesh%i_enumfine
+ WRITE(i_fhandle, *) i_cnt
+ END DO
+
+ WRITE(i_fhandle, *) '</DataArray>'
+ END IF
+
+ !> Write the cell variables (when present)
+ IF(PRESENT(i_celldata) .AND. PRESENT(p_celldata)) THEN
+ DO i_cnt = 1, i_celldata
+ CALL write_vtu_data(i_fhandle, p_celldata(i_cnt))
+ END DO
+ END IF
+
+ !> write the footer.
+ WRITE(i_fhandle, *) '</CellData>'
+ WRITE(i_fhandle, *) '</Piece>'
+ WRITE(i_fhandle, *) '</UnstructuredGrid>'
+ WRITE(i_fhandle, *) '</VTKFile>'
+
+ !> finally close the file and deallocate the data
+ CLOSE(i_fhandle)
+ DEALLOCATE(r_nodecoor, i_enodes)
+END SUBROUTINE plot_vtu_elements
+#endif !----3D output only
+
+!---------------------------------------------------------------------
+!> write_vtu_data writes a single variable (node or cell) to the
+!! VTK file.
+!> \param i_fhandle - the file handle
+!> \param p_data - the variable
+!---------------------------------------------------------------------
+SUBROUTINE write_vtu_data(i_fhandle, p_data)
+ IMPLICIT NONE
+
+ INTEGER(KIND=GRID_SI), INTENT(IN) :: i_fhandle
+ TYPE(t_vtu_data) :: p_data
+ INTEGER(KIND=GRID_SI) :: i_size
+ INTEGER(KIND=GRID_SI) :: i_cnt
+ INTEGER(KIND=GRID_SI) :: i_cnt2
+
+ !> Special treatment for vector valued data
+ IF(p_data%i_size > 1) THEN
+ WRITE(i_fhandle, "(A,A,A,I1, A)") &
+ '<DataArray type="Float32" Name="', &
+ TRIM(p_data%c_name), '" NumberOfComponents="', &
+ p_data%i_size, '" format="ascii">'
+ ELSE
+ WRITE(i_fhandle, "(A,A,A)") '<DataArray type="Float32" Name="', &
+ TRIM(p_data%c_name), '" format="ascii">'
+ END IF
+ i_size = SIZE(p_data%p_vdata, 2)
+
+ !> Write the data to file with respect to the threshold value
+ DO i_cnt = 1, i_size
+ DO i_cnt2 = 1, p_data%i_size
+ IF(ABS(p_data%p_vdata(i_cnt2, i_cnt)) < r_vtueps) &
+ p_data%p_vdata(i_cnt2, i_cnt) = 0
+ END DO
+
+ WRITE(i_fhandle, "(10e15.7)") p_data%p_vdata(1:p_data%i_size, i_cnt)
+ END DO
+ WRITE(i_fhandle, *) '</DataArray>'
+
+END SUBROUTINE write_vtu_data
+
+END MODULE IO_vtu
diff --git a/flash2d/src/options-sphere/IO_vtuplot_stacked.F90 b/flash2d/src/options-sphere/IO_vtuplot_stacked.F90
new file mode 100644
index 0000000000000000000000000000000000000000..aebd0d4b13537459da479355ab5f19b8f479e541
--- /dev/null
+++ b/flash2d/src/options-sphere/IO_vtuplot_stacked.F90
@@ -0,0 +1,185 @@
+!*****************************************************************
+!
+!> @file IO_vtuplot.F90
+!> @brief includes module IO_vtuplot
+!
+!*****************************************************************
+!
+! VERSION(S):
+! 1. original version j. behrens 04/2000
+! 2. amatos-1.0 and 2D compliant j. behrens 11/2000
+! 3. adapted to flash2d j. behrens 12/2014
+!
+!*****************************************************************
+! MODULE DESCRIPTION:
+!> prints data in VTU compatible file format
+!
+MODULE IO_vtuplot
+ USE FLASH_parameters
+ USE MISC_outputdata
+ USE GRID_api
+ USE IO_vtu
+ PRIVATE
+ INTEGER, SAVE :: i_timecounter= 0
+ PUBLIC :: generate_vtu
+ CONTAINS
+
+!*****************************************************************
+! DESCRIPTION of [SUBROUTINE generate_vtu]:
+!> @brief creates output in GMV compatible file format
+!>
+!> @param[in] p_handle grid handle for the linked lists
+!> @param[in] i_time time step number
+!> @param[in] i_newslevel the time level
+!>
+!> @note ON OUTPUT: a number of VTU compatible ascii files for each plotted time step
+!
+ SUBROUTINE generate_vtu(p_handle, i_time, i_newslevel)
+
+ IMPLICIT NONE
+
+!---------- local declarations
+
+ TYPE (grid_handle), INTENT(in) :: p_handle
+ INTEGER (KIND = GRID_SI), OPTIONAL, INTENT(in) :: i_time
+ INTEGER (KIND = GRID_SI), OPTIONAL, INTENT(in) :: i_newslevel
+ INTEGER (KIND = GRID_SI) :: i_alct, i_nnum, i_tnum
+ INTEGER (KIND = GRID_SI) :: i_nwl, i_elen, i_tcnt
+ REAL (KIND = GRID_SR), DIMENSION(:,:), POINTER :: r_val
+ REAL (KIND = GRID_SR), DIMENSION(:,:), POINTER :: r_velo
+ REAL (KIND = GRID_SR), DIMENSION(:,:), POINTER :: r_sta, r_lvl
+ INTEGER (KIND = GRID_SI), DIMENSION(:), ALLOCATABLE :: i_ids, i_sta
+ CHARACTER (len=32) :: c_mfile
+ CHARACTER (len=26) :: c_tmp
+ LOGICAL :: l_news
+ INTEGER, PARAMETER :: i_vallen=5
+ INTEGER (KIND = GRID_SI), DIMENSION(:), ALLOCATABLE :: i_valind
+ INTEGER (KIND = GRID_SI) :: i_trcvals, i_lay, i_nodnum
+
+ TYPE(t_vtu_data), DIMENSION(2) :: celldata
+ ! CAUTION: This is a fixed dimension, but will only work for max 16 layers...
+ TYPE(t_vtu_data), DIMENSION(20) :: nodedata
+
+!---------- file handling (open)
+
+ IF(present(i_time)) THEN
+ i_tcnt= i_time
+ i_timecounter= i_timecounter+1
+ ELSE
+ i_tcnt= i_timecounter
+ i_timecounter= i_timecounter+1
+ END IF
+ write(c_tmp,*) trim(GRID_parameters%program_name(1:20))
+ write(c_mfile,10101) trim(c_tmp), '_', i_tcnt
+ c_tmp = adjustl(c_mfile)
+ write(c_mfile,*) TRIM(c_tmp), '.vtu'
+ c_mfile= adjustl(c_mfile)
+
+!---------- handle news level
+
+ IF(present(i_newslevel)) THEN
+ l_news=.TRUE.
+ i_nwl= i_newslevel
+ ELSE
+ l_news=.FALSE.
+ END IF
+
+!---------- the nodes
+
+ i_nnum = p_handle%i_nnumber
+
+!---------- determine number of tracer values from layers and allocate accordingly
+
+ i_trcvals = i_vallen + i_stcnumlayers
+
+!---------- extract nodal grid data and nodal values
+
+ ALLOCATE(r_val(i_trcvals, i_nnum), i_valind(i_trcvals), stat=i_alct)
+ IF(i_alct /= 0) THEN
+ CALL grid_error(c_error='[generate_vtu]: could not allocate values array')
+ END IF
+ i_valind= (/GRID_ucomp, GRID_vcomp, GRID_phi, GRID_zeta, GRID_tracer, &
+ ((i_stctracer-1+i_lay), i_lay=1,i_stcnumlayers) /)
+ CALL grid_getinfo(p_handle, r_nodevalues= r_val, i_arraypoint=i_valind)
+
+!---------- assign nodal values, but before allocate aux. arrays
+
+ ALLOCATE(r_velo(2,i_nnum), stat=i_alct)
+ IF(i_alct /= 0) THEN
+ CALL grid_error(c_error='[generate_vtu]: could not allocate aux. velo arrays')
+ END IF
+
+!---------- generate 2D velocity field
+
+ nodedata(1)%c_name = 'u-comp'
+ nodedata(1)%i_size = 1
+ nodedata(1)%p_vdata => r_val(1:1,:)
+
+ nodedata(2)%c_name = 'v-comp'
+ nodedata(2)%i_size = 1
+ nodedata(2)%p_vdata => r_val(2:2,:)
+
+ nodedata(3)%c_name = 'height'
+ nodedata(3)%i_size = 1
+ nodedata(3)%p_vdata => r_val(3:3,:)
+
+ nodedata(4)%c_name = 'vorticity'
+ nodedata(4)%i_size = 1
+ nodedata(4)%p_vdata => r_val(4:4,:)
+
+ nodedata(5)%c_name = 'tracer'
+ nodedata(5)%i_size = 1
+ nodedata(5)%p_vdata => r_val(5:5,:)
+
+!---------- loop over all layers and store them in separate tracer layers...
+
+ DO i_lay=1,i_stcnumlayers
+ write(c_tmp,*) 'layer',i_lay
+ nodedata(i_lay+5)%c_name = c_tmp
+ nodedata(i_lay+5)%i_size = 1
+ nodedata(i_lay+5)%p_vdata => r_val(i_lay+5:i_lay+5,:)
+ END DO
+ i_nodnum= i_stcnumlayers+ 4
+
+!---------- extract cells grid data
+
+ i_tnum= p_handle%i_enumfine
+ ALLOCATE(i_ids(i_tnum), i_sta(i_tnum), r_sta(1,i_tnum), r_lvl(1,i_tnum), stat=i_alct)
+ IF(i_alct /= 0) THEN
+ CALL grid_error(c_error='[generate_vtu]: could not allocate cell value arrays')
+ END IF
+ IF(l_news) THEN
+ CALL grid_getinfo(p_handle, i_newsdepth= i_nwl, i_elength= i_elen, &
+ i_elementlevel= i_ids, i_elementstatus= i_sta)
+ i_tnum= i_elen
+ ELSE
+ CALL grid_getinfo(p_handle, i_elementlevel= i_ids, i_elementstatus= i_sta)
+ END IF
+
+ r_sta(1, :) = REAL(i_sta, GRID_SR)
+ r_lvl(1, :) = REAL(i_ids, GRID_SR)
+
+ ! Assign all values to driver
+ celldata(1)%c_name = 'status'
+ celldata(1)%i_size = 1
+ celldata(1)%p_vdata => r_sta
+
+ celldata(2)%c_name = 'level'
+ celldata(2)%i_size = 1
+ celldata(2)%p_vdata => r_lvl
+
+ ! plot the vtu data
+ CALL plot_vtu(p_handle, c_mfile, i_nodedata = i_nodnum, p_nodedata = nodedata, &
+ i_celldata = 2, p_celldata = celldata) !, i_zcoordinate = 4)
+
+ DEALLOCATE(i_ids, i_sta, r_sta, r_lvl, &
+ r_val, r_velo)
+
+ RETURN
+
+!---------- FORMAT
+ 10101 FORMAT(a19,a1,i6.6)
+
+ END SUBROUTINE generate_vtu
+
+END MODULE IO_vtuplot
diff --git a/flash2d/src/options-sphere/SLM_advanced_stacked.f90 b/flash2d/src/options-sphere/SLM_advanced_stacked.f90
new file mode 100644
index 0000000000000000000000000000000000000000..b0e7ae57923e9950c7d85c5f5027c75a780cf20d
--- /dev/null
+++ b/flash2d/src/options-sphere/SLM_advanced_stacked.f90
@@ -0,0 +1,525 @@
+!*****************************************************************
+!
+! MODULE NAME:
+! SLM_advanced
+! FUNCTION:
+! provide simple semi-Lagrangian routines
+! CONTAINS:
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_step
+! FUNCTION:
+! one step of the basic SLM algorithm
+! SYNTAX:
+! CALL slm_step(int, real.arr, real.arr)
+! ON INPUT:
+! ...
+! ON OUTPUT:
+! r_tracer: array with tracer values real
+! CALLS:
+!
+! COMMENTS:
+!
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_displace
+! FUNCTION:
+! extrapolate the alpha, values for the displacements of the upstream
+! points from the gridpoints
+! SYNTAX:
+! CALL slm_displace(int, real.arr, real.arr)
+! ON INPUT:
+! i_arlen: array length for the real arrays integer
+! r_coord: real array of xy-coordinates real
+! ON OUTPUT:
+! r_alpha: displacement vectors to each point real
+! CALLS:
+! wind_field
+! COMMENTS:
+!
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_update
+! FUNCTION:
+! calculate the update to the velocity
+! SYNTAX:
+! CALL slm_update(int, real.arr, real.arr)
+! ON INPUT:
+! i_arlen: array length for the real arrays integer
+! r_rside: array with right hand side values real
+! ON OUTPUT:
+! r_udate: array with new (updated) gid values real
+! CALLS:
+!
+! COMMENTS:
+! this routine is trivial for linear advection
+!
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_upstream
+! FUNCTION:
+! calculate right hand side of the equation (upstream values)
+! SYNTAX:
+! CALL slm_upstream(int, real.arr, real.arr)
+! ON INPUT:
+! i_arlen: array length for the real arrays integer
+! r_alpha: displacement vectors to each point real
+! ON OUTPUT:
+! r_rside: array with right hand side values real
+! CALLS:
+!
+! COMMENTS:
+! this routine is just interpolation for linear advection
+!
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_interpolate
+! FUNCTION:
+! do the interpolation
+! SYNTAX:
+! CALL slm_interpolate(grid, int, real, real.arr, real.arr, real.arr)
+! ON INPUT:
+! p_ogrid: grid handle to old grid (with data) TYPE (grid_handle)
+! r_fac: factor at which point to interpolate REAL
+! i_arlen: array length for the following arrays INTEGER
+! r_coord: coordinate array (new grid) REAL
+! r_alpha: displacement array (corr. to r_coord) REAL
+! r_value: values on the old grid (array) REAL
+! ON OUTPUT:
+! r_rside: right hand side (interpolated) REAL
+! CALLS:
+!
+! COMMENTS:
+! this one is plain bi-cubic spline interpolation
+!
+!-----------------------------------------------------------------
+!
+! PUBLIC:
+! slm_step
+! COMMENTS:
+!
+! USES:
+! FLASH_parameters, GRID_api, ADV_wind, ADV_rhs
+! LIBRARIES:
+!
+! REFERENCES:
+!
+! VERSION(S):
+! 1. original version j. behrens 4/2002
+! 2. compliant to amatos 2.0 j. behrens 7/2003
+!
+!*****************************************************************
+ MODULE SLM_advanced
+ USE FLASH_parameters
+ USE MISC_timing
+ USE GRID_api
+ USE ADV_wind
+ USE ADV_rhs
+ PRIVATE
+ PUBLIC :: slm_astep, slm_averticalstep
+ CONTAINS
+
+!*****************************************************************
+ SUBROUTINE slm_astep(p_ghand, p_param, p_time, r_modtime, i_size, &
+ r_coord, r_tracer, i_newsdepth, i_layer)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ TYPE (grid_handle), DIMENSION(GRID_timesteps), INTENT(in) :: p_ghand
+ TYPE (control_struct), INTENT(in) :: p_param
+ TYPE (sw_info), INTENT(inout) :: p_time
+ REAL (KIND = GRID_SR), INTENT(in) :: r_modtime
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_size
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension,i_size), INTENT(in) :: r_coord
+ REAL (KIND = GRID_SR), DIMENSION(i_size), INTENT(out) :: r_tracer
+ INTEGER (KIND = GRID_SI), OPTIONAL, INTENT(in) :: i_newsdepth
+ INTEGER (KIND = GRID_SI), OPTIONAL, INTENT(in) :: i_layer
+
+ REAL (KIND = GRID_SR), DIMENSION(:), ALLOCATABLE :: r_newvl
+ REAL (KIND = GRID_SR), DIMENSION(:,:), ALLOCATABLE :: r_alpha
+ INTEGER (KIND = GRID_SI) :: i_alct, i_layerid
+
+!---------- evaluate layer information, if given
+
+ layer_info: IF(present(i_layer)) THEN
+ i_layerid = i_layer
+ ELSE layer_info
+ i_layerid = 1 ! Default is the lower layer (layer 1)
+ END IF layer_info
+
+!---------- check size!
+
+ IF(i_size <= 0) THEN
+ IF(GRID_parameters%iolog > 0) &
+ write(GRID_parameters%iolog,*) 'INFO [slm_step]: Zero step size, returning to calling routine'
+ RETURN
+ END IF
+
+!---------- allocate auxiliary arrays
+
+ allocate(r_newvl(i_size), r_alpha(GRID_dimension,i_size), stat=i_alct)
+ not_alloc: IF(i_alct /= 0) THEN
+ CALL grid_error(38)
+ END IF not_alloc
+
+!-SLM--------- calculate trajectory pieces (displacements)
+
+ CALL stop_watch('start',3,p_time)
+ CALL slm_displace(p_param, i_size, i_layerid, r_coord, r_alpha, r_time=r_modtime)
+ CALL stop_watch('stop ',3,p_time)
+
+!-SLM--------- calculate right hand side
+
+ CALL stop_watch('start',4,p_time)
+ CALL slm_upstream(p_ghand, i_size, i_layerid, r_coord, r_alpha, r_newvl)
+ CALL stop_watch('stop ',4,p_time)
+
+!-SLM--------- calculate new grid values
+
+ CALL stop_watch('start',5,p_time)
+ CALL slm_update(p_param, i_size, i_layerid, r_coord, r_newvl, r_tracer, r_time=r_modtime)
+ CALL stop_watch('stop ',5,p_time)
+
+!-SLM--------- put alpha values to u and v field entries
+
+! r_alpha= -r_alpha
+! IF(present(i_newsdepth)) THEN
+! CALL grid_putinfo(p_ghand(i_timeplus), r_nodevalues=r_alpha, &
+! i_newsdepth=i_newsdepth, i_arraypoint=(/ GRID_ucomp, GRID_vcomp /))
+! ELSE
+! CALL grid_putinfo(p_ghand(i_timeplus), r_nodevalues=r_alpha, &
+! i_arraypoint=(/ GRID_ucomp, GRID_vcomp /))
+! END IF
+
+!-SLM--------- deallocate work arrays
+
+ deallocate(r_alpha, r_newvl)
+
+ RETURN
+ END SUBROUTINE slm_astep
+
+!*****************************************************************
+ SUBROUTINE slm_averticalstep(p_ghand, p_param, p_time, r_modtime, i_size, &
+ r_coord, r_tracer, i_newsdepth)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ TYPE (grid_handle), DIMENSION(GRID_timesteps), INTENT(in) :: p_ghand
+ TYPE (control_struct), INTENT(in) :: p_param
+ TYPE (sw_info), INTENT(inout) :: p_time
+ REAL (KIND = GRID_SR), INTENT(in) :: r_modtime
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_size
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension,i_size), INTENT(in) :: r_coord
+ REAL (KIND = GRID_SR), DIMENSION(:,:), INTENT(inout) :: r_tracer
+ INTEGER (KIND = GRID_SI), OPTIONAL, INTENT(in) :: i_newsdepth
+
+ REAL (KIND = GRID_SR), DIMENSION(:), ALLOCATABLE :: r_newvl
+ REAL (KIND = GRID_SR), DIMENSION(:,:), ALLOCATABLE :: r_alpha
+ INTEGER (KIND = GRID_SI) :: i_alct, i_layerid
+
+ REAL (KIND = GRID_SR) :: r_upstr, r_dd
+ INTEGER (KIND = GRID_SI) :: i_cnt, i_lyr, &
+ i_lcnt
+
+!---------- allocate intermediate storage for the new tracer values
+
+ allocate(r_newvl(i_stcnumlayers), stat=i_alct)
+ not_alloc: IF(i_alct /= 0) THEN
+ CALL grid_error(c_error='[slm_verticalstep]: could not allocate r_newvl')
+ END IF not_alloc
+
+!---------- outer loop over all coordinates
+ coord_loop: DO i_cnt=1,i_size
+
+!---------- inner loop over height levels
+ layer_loop: DO i_lyr=1,i_stcnumlayers
+
+!---------- an explicit Euler step for the upwind position
+ r_upstr = r_stclayerpressure(i_lyr) - p_param%num%r_deltatime* &
+ slm_vertwind(r_coord(:,i_cnt),i_lyr)
+
+!---------- upstream value interpolation (linear)
+! Step 1: check for boundary
+! upstream position is higher than highest layer
+ bound_check: IF(r_upstr <= r_stclayerpressure(i_stcnumlayers)) THEN
+ r_newvl(i_lyr) = 0._GRID_SR
+! upstream position is lower than ground layer
+ ELSE IF(r_upstr > r_stclayerpressure(1)) THEN bound_check
+ r_newvl(i_lyr) = 0._GRID_SR
+! upstream position is within vertical domain
+ ELSE bound_check
+! Step 2: interpolate upstream position - use linear interpolation
+! find vertical box of upstream position
+ layer_search: DO i_lcnt=1,i_stcnumlayers-1
+ IF((r_upstr <= r_stclayerpressure(i_lcnt)) .AND. &
+ (r_upstr > r_stclayerpressure(i_lcnt+1))) THEN
+ i_layerid = i_lcnt
+ EXIT layer_search
+ END IF
+ END DO layer_search
+! now the linear interpolation
+ r_dd = (r_tracer(i_layerid+1, i_cnt)- r_tracer(i_layerid, i_cnt))/ &
+ (r_stclayerpressure(i_layerid+1)- r_stclayerpressure(i_layerid))
+ r_newvl(i_lyr) = r_tracer(i_layerid,i_cnt) + r_dd* &
+ (r_upstr - r_stclayerpressure(i_layerid))
+ END IF bound_check
+
+!---------- finally, the update with a possible right hand side
+ r_newvl(i_lyr) = r_newvl(i_lyr) + p_param%num%r_deltatime* &
+ slm_verticalrhs(r_coord(:,i_cnt),i_lyr,r_modtime)
+ END DO layer_loop
+
+!---------- overwrite the tracer values with the newly computed advected values in this column
+ r_tracer(:,i_cnt) = r_newvl(:)
+
+ END DO coord_loop
+
+!---------- deallocate intermediate storage
+
+ deallocate(r_newvl)
+
+ RETURN
+ END SUBROUTINE slm_averticalstep
+
+!*****************************************************************
+ SUBROUTINE slm_displace(p_param, i_arlen, i_layerid, r_coord, r_alpha, r_time)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ TYPE (control_struct), INTENT(in) :: p_param
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_arlen
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_layerid
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension,i_arlen), INTENT(in) :: r_coord
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension,i_arlen), INTENT(out) :: r_alpha
+ REAL (KIND = GRID_SR), INTENT(in), OPTIONAL :: r_time
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_fac, r_caf, &
+ r_axy, r_xyc, r_sxy
+ REAL (KIND = GRID_SR) :: r_dt0, r_dt1, &
+ r_dt2, r_tim
+ INTEGER (KIND = GRID_SI) :: i_cnt1, i_cnt2
+
+!---------- set constants
+
+ r_dt0= p_param%num%r_deltatime
+ r_dt1= 0.5* p_param%num%r_deltatime
+ r_dt2= 1.5* p_param%num%r_deltatime
+ r_fac= 0.5
+ r_caf= 2.0
+ IF(present(r_time)) THEN
+ r_tim= r_time
+ ELSE
+ r_tim= 0.0
+ END IF
+
+!---------- calculate in an iteration process the displacements
+
+ unknown_loop: DO i_cnt1=1,i_arlen
+ r_axy= 0.0
+
+ iter_loop: DO i_cnt2=1, p_param%num%i_adviterations
+ r_xyc= r_coord(:,i_cnt1)- r_fac* r_axy
+ r_sxy= r_dt0* slm_windfield(r_xyc, i_layerid, r_time=r_tim)
+ r_axy= sphere_correct(r_coord(:,i_cnt1),r_sxy)
+ END DO iter_loop
+
+ r_alpha(:,i_cnt1)= r_axy
+ END DO unknown_loop
+
+ RETURN
+ END SUBROUTINE slm_displace
+
+!*****************************************************************
+ SUBROUTINE slm_update(p_param, i_arlen, i_layer, r_coord, r_rside, r_udate, r_time)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ TYPE (control_struct), INTENT(in) :: p_param
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_arlen
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_layer
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension,i_arlen), INTENT(in) :: r_coord
+ REAL (KIND = GRID_SR), DIMENSION(i_arlen), INTENT(in) :: r_rside
+ REAL (KIND = GRID_SR), DIMENSION(i_arlen), INTENT(out) :: r_udate
+ REAL (KIND = GRID_SR), INTENT(in), OPTIONAL :: r_time
+
+ INTEGER (KIND = GRID_SI) :: i_cnt
+ REAL (KIND = GRID_SR) :: r_dt, r_tim
+
+!---------- in the linear advection case and with f90 this is just
+
+! r_udate= r_rside
+
+!---------- including a non-zero right hand side, we have
+
+ r_dt= p_param%num%r_deltatime
+ IF(present(r_time)) THEN
+ r_tim= r_time
+ ELSE
+ r_tim= 0.0
+ END IF
+
+ main_loop: DO i_cnt=1, i_arlen
+ r_udate(i_cnt)= r_rside(i_cnt)+ r_dt* slm_righthand(r_coord(:,i_cnt), i_layer)
+ END DO main_loop
+
+ RETURN
+ END SUBROUTINE slm_update
+
+!*****************************************************************
+ SUBROUTINE slm_upstream(p_mesh, i_arlen, i_layer, r_coord, &
+ r_alpha, r_rside)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ TYPE (grid_handle), DIMENSION(GRID_timesteps) :: p_mesh
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_arlen
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_layer
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension,i_arlen), INTENT(in) :: r_coord
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension,i_arlen), INTENT(in) :: r_alpha
+ REAL (KIND = GRID_SR), DIMENSION(i_arlen), INTENT(out) :: r_rside
+
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_fac
+
+!---------- set factor (at which point of trajectory shall i interpolate)
+
+ r_fac= 1.0
+
+!---------- in the linear advection case this is just interpolation
+
+ CALL slm_interpolate(p_mesh, r_fac, i_arlen, i_layer, r_coord, &
+ r_alpha, r_rside)
+
+ RETURN
+ END SUBROUTINE slm_upstream
+!*****************************************************************
+ SUBROUTINE slm_interpolate(p_mesh, r_fac, i_arlen, i_layer, &
+ r_coord, r_alpha, r_rside)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ TYPE (grid_handle), DIMENSION(GRID_timesteps) :: p_mesh
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension), INTENT(in) :: r_fac
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_arlen
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_layer
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension,i_arlen), INTENT(in) :: r_coord
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension,i_arlen), INTENT(in) :: r_alpha
+ REAL (KIND = GRID_SR), DIMENSION(i_arlen), INTENT(out) :: r_rside
+
+ REAL (KIND = GRID_SR), DIMENSION(:,:), ALLOCATABLE :: r_upstr
+ REAL (KIND = GRID_SR) :: r_eps
+ INTEGER (KIND = GRID_SI) :: i_cnt, &
+ i_alct, i_val, i_out, i_stat
+
+!---------- initialize constant
+
+ i_val= i_stctracer-1+i_layer
+ r_eps= GRID_EPS
+
+!---------- allocate work array
+
+ ALLOCATE(r_upstr(GRID_dimension,i_arlen), stat=i_alct)
+ not_allocated: IF(i_alct /= 0) THEN
+ CALL grid_error(60)
+ END IF not_allocated
+
+!---------- calculate upstream coordinates
+
+ dim_loop: DO i_cnt=1, GRID_dimension
+ r_upstr(i_cnt,:) = r_coord(i_cnt,:)- r_fac(i_cnt)* r_alpha(i_cnt,:)
+ END DO dim_loop
+
+!---------- loop over nodes: find element containing upstream point
+
+ node_loop: DO i_cnt=1, i_arlen
+
+!---------- check if upstream value is outside of the domain
+
+ i_out= grid_domaincheck(p_mesh(i_time), r_upstr(:,i_cnt))
+
+!---------- take the intersection of the trajectory with the boundary as new upstream point
+
+ out_domain: IF(i_out /= 0) then
+ r_upstr(:,i_cnt)= grid_boundintersect(p_mesh(i_time), &
+ r_coord(:,i_cnt), r_upstr(:,i_cnt), i_info=i_stat)
+ no_intersect: IF(i_stat /= 0) THEN
+ r_rside(i_cnt)= 0.0
+ CYCLE node_loop
+ END IF no_intersect
+ END IF out_domain
+
+!---------- interpolate
+
+ r_rside(i_cnt)= grid_coordvalue(p_mesh(i_time), r_upstr(:,i_cnt), &
+ i_interpolorder=GRID_loworder, i_valpoint=i_val)
+ small_val: IF(abs(r_rside(i_cnt)) < r_eps) THEN
+ r_rside(i_cnt)= 0.0
+ END IF small_val
+
+ END DO node_loop
+
+!---------- deallocate work array
+
+ DEALLOCATE(r_upstr)
+
+ RETURN
+ END SUBROUTINE slm_interpolate
+
+!*****************************************************************
+ FUNCTION sphere_correct(r_coord, r_displ) RESULT (r_corct)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_coord, r_displ
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_corct
+
+ REAL (KIND = GRID_SR) :: r_e, r_rat
+ REAL (KIND = GRID_SR) :: r_c
+
+!---------- calculate Euklidean norm
+
+ r_e = eukl_norm(r_displ)
+ r_rat= r_e/ GRID_RADIUS
+
+!---------- calculate correction
+
+ r_c = tan(r_rat)/r_rat
+ r_corct= r_displ* r_c
+
+ RETURN
+ END FUNCTION sphere_correct
+
+!*****************************************************************
+ FUNCTION eukl_norm(r_vec) RESULT(r_rst)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_vec
+ REAL (KIND = GRID_SR) :: r_rst, r_tmp
+
+!---------- calculate vector crossproduct, 3D only
+
+ r_tmp= dot_product(r_vec, r_vec)
+ r_rst= sqrt(r_tmp)
+
+ END FUNCTION eukl_norm
+
+ END MODULE SLM_advanced
diff --git a/flash2d/src/options-sphere/SLM_errorestimate_stacked.f90 b/flash2d/src/options-sphere/SLM_errorestimate_stacked.f90
new file mode 100644
index 0000000000000000000000000000000000000000..49740b135d819ca59f34ee42e35f903202969537
--- /dev/null
+++ b/flash2d/src/options-sphere/SLM_errorestimate_stacked.f90
@@ -0,0 +1,154 @@
+!*****************************************************************
+!
+! MODULE NAME:
+! SLM_errorestimate
+! FUNCTION:
+! provide an error estimator for the adaptive SLM scheme
+! CONTAINS:
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_error
+! FUNCTION:
+! estimate the error elementwise
+! SYNTAX:
+! real= slm_error(real.arr)
+! ON INPUT:
+! r_v: values at element nodes REAL
+! ON OUTPUT:
+! r_esterr: estimated error REAL
+! CALLS:
+!
+! COMMENTS:
+! this is only a simple estimator based on local gradients
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_errorest
+! FUNCTION:
+! this hides the loop structure from the advection part
+! SYNTAX:
+! CALL slm_errorest(grid, int, real.arr)
+! ON INPUT:
+! p_ghand: handle for the grid TYPE (grid_handle)
+! i_siz: size of array for local errors INTEGER
+! ON OUTPUT:
+! i_siz: size of array for local errors INTEGER
+! r_arr: estimated error REAL
+! CALLS:
+!
+! COMMENTS:
+!
+!-----------------------------------------------------------------
+!
+! PUBLIC:
+! slm_errorest
+! COMMENTS:
+!
+! USES:
+! MISC_globalparam, MISC_error, GRID_api
+! LIBRARIES:
+!
+! REFERENCES:
+!
+! VERSION(S):
+! 1. original version j. behrens 11/96
+! 2. nodal values time depend. j. behrens 1/97
+! 3. FEM_handle added j. behrens 7/97
+! 4. version with hashing j. behrens 7/97
+! 5. changed to use GRID_api j. behrens 11/97
+! 6. compliant to amatos 1.0 j. behrens 12/2000
+! 7. compliant to amatos 1.2 j. behrens 3/2002
+! 8. compliant to amatos 2.0 f. klaschka 8/2006
+!
+!*****************************************************************
+ MODULE SLM_errorestimate
+ USE FLASH_parameters
+ USE GRID_api
+ PRIVATE
+ PUBLIC :: slm_errorest
+ CONTAINS
+!*****************************************************************
+ FUNCTION slm_error(r_v) RESULT (r_esterr)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+ REAL (KIND = GRID_SR), DIMENSION(GRID_elementnodes) :: r_v
+ REAL (KIND = GRID_SR) :: r_esterr
+ REAL (KIND = GRID_SR), PARAMETER :: r_1o3= (1.0/3.0)
+ REAL (KIND = GRID_SR) :: r_d1, r_d2, r_d3
+
+!---------- calculate differences
+
+ r_d1 = abs(r_v(1)- r_v(2))
+ r_d2 = abs(r_v(2)- r_v(3))
+ r_d3 = abs(r_v(3)- r_v(1))
+
+!---------- this is the estimated error
+
+ r_esterr= (r_d1+ r_d2+ r_d3)* r_1o3
+
+ RETURN
+ END FUNCTION slm_error
+
+!*****************************************************************
+ SUBROUTINE slm_errorest(p_ghand, i_siz, r_arr)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+ TYPE (grid_handle), INTENT(in) :: p_ghand
+ INTEGER, INTENT(in) :: i_siz
+ REAL (KIND = GRID_SR), DIMENSION(:), INTENT(out) :: r_arr
+ REAL (KIND = GRID_SR), DIMENSION(:,:), ALLOCATABLE :: r_aux
+ INTEGER, DIMENSION(:,:), ALLOCATABLE :: i_enods
+ INTEGER :: i_cnt, i_size, i_alct, i_nnum, i_lay, i_start
+ INTEGER, DIMENSION(1) :: i_valind
+ REAL (KIND = GRID_SR), DIMENSION(GRID_elementnodes) :: r_tmp
+
+!---------- allocate auxilliary array
+
+ i_size= p_ghand%i_enumfine
+ i_nnum= p_ghand%i_nnumber
+ IF(i_siz /= i_size) THEN
+ CALL grid_error(48)
+ END IF
+
+ ALLOCATE(r_aux(1,i_nnum), &
+ i_enods(GRID_elementnodes,i_size), stat= i_alct)
+ IF(i_alct /= 0) THEN
+ CALL grid_error(49)
+ END IF
+
+!---------- initialize r_arr, since we want to sum over all layers
+
+ r_arr= 0._GRID_SR
+
+!---------- loop over stacked layers
+
+ i_start = i_stctracer - 1
+ layer_loop: DO i_lay=1,i_stcnumlayers
+
+!---------- get values
+
+ i_valind= (/ i_start + i_lay /)
+ CALL grid_getinfo(p_ghand, l_finelevel= .TRUE., i_arraypoint=i_valind, &
+ r_nodevalues= r_aux, i_elementnodes=i_enods)
+
+!---------- loop through all elements of finest loop
+
+ elmt_loop: DO i_cnt=1,i_siz
+ r_tmp(1:GRID_elementnodes)= r_aux(1,i_enods(1:GRID_elementnodes,i_cnt))
+ r_arr(i_cnt)= r_arr(i_cnt) + slm_error(r_tmp)
+ END DO elmt_loop
+ END DO layer_loop
+
+!---------- deallocate and exit
+
+ DEALLOCATE(r_aux, i_enods)
+
+ RETURN
+ END SUBROUTINE slm_errorest
+
+ END MODULE SLM_errorestimate
diff --git a/flash2d/src/options-sphere/SLM_initial_stacked.f90 b/flash2d/src/options-sphere/SLM_initial_stacked.f90
new file mode 100644
index 0000000000000000000000000000000000000000..566891c478d0c376e57a509f012daaca217486bb
--- /dev/null
+++ b/flash2d/src/options-sphere/SLM_initial_stacked.f90
@@ -0,0 +1,360 @@
+!*****************************************************************
+!
+! MODULE NAME:
+! SLM_initial
+! FUNCTION:
+! initialize slotted cylinder for semi-Lagrangian advection
+! CONTAINS:
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_initialvalues
+! FUNCTION:
+! initialize a grid with values
+! SYNTAX:
+! CALL slm_initialvalues(grid)
+! ON INPUT:
+! p_ghand: grid handle TYPE (grid_handle)
+! ON OUTPUT:
+! p_ghand: grid handle TYPE (grid_handle)
+! CALLS:
+!
+! COMMENTS:
+! the routine is made for two dimensions only
+!
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_initslot
+! FUNCTION:
+! initialize a grid with values of slotted cylinder test case
+! SYNTAX:
+! CALL slm_initslot(grid)
+! ON INPUT:
+! p_ghand: grid handle TYPE (grid_handle)
+! ON OUTPUT:
+! p_ghand: grid handle TYPE (grid_handle)
+! CALLS:
+!
+! COMMENTS:
+! the routine is made for two dimensions only
+!
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_analyticsolution
+! FUNCTION:
+! calculates the 'analytic solution' to compare with in diagnostics
+! SYNTAX:
+! CALL slm_analyticsolution(grid, real, int, real.arr)
+! ON INPUT:
+! p_ghand: grid handle TYPE (grid_handle)
+! r_time: model time REAL
+! i_arlen: array length for values array INTEGER
+! ON OUTPUT:
+! r_array: values at gridpoints REAL
+! CALLS:
+!
+! COMMENTS:
+! the routine is made for two dimensions only
+!
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_initcylndr
+! FUNCTION:
+! initialize a grid with values of a cylinder test case
+! SYNTAX:
+! CALL slm_initcylndr(grid)
+! ON INPUT:
+! p_ghand: grid handle TYPE (grid_handle)
+! r_centr: coordinates of cyl. centre REAL
+! ON OUTPUT:
+! p_ghand: grid handle TYPE (grid_handle)
+! CALLS:
+!
+! COMMENTS:
+! the routine is made for two dimensions only
+!
+!-----------------------------------------------------------------
+!
+! NAME:
+! in_side
+! FUNCTION:
+! checks, if a given point (x,y) lies within a given triangle
+! SYNTAX:
+! logical= in_side(real.arr, real.arr, real.arr, real.arr)
+! ON INPUT:
+! r_coord: coordinate array REAL
+! r_node1: node1 of triangle REAL
+! r_node2: node2 of triangle REAL
+! r_node3: node3 of triangle REAL
+! ON OUTPUT:
+! l_in: .true. if r_coord \in p_elem logical
+! CALLS:
+!
+! COMMENTS:
+! this routine decides whether a point lies in or out of a
+! triangle. this is done with the following approach:
+! calculate the area for the given triangle and for the
+! three triangles resulting from drawing lines from the given
+! point to all three triangle nodes.
+! if the sum of the areas of those three trianlges exceeds
+! the area of the given trianlge, the the point lies outside
+! of it.
+! for calculation of the areas following formula is used:
+! (Herons formula(?))
+!
+! A = sqrt(s* (s- a)* (s- b)* (s- c)),
+!
+! where s= 1/2* (a+ b+ c)
+! a, b, c sides.
+!
+! in order to calculate the sidelengths |x-y|= sqrt(x**2-y**2)
+! the complex absolute value intrinsic function is used.
+! hopefully this fuction is faster than using sqrt and
+! power-of-two instead.
+!
+! internally double precision is used in this function, because
+! machine precision is crucial when a coordinate pair appears
+! near the edge or corner of an element.
+!
+!-----------------------------------------------------------------
+!
+! NAME:
+! dc_area
+! FUNCTION:
+! calculate area of a triangle (in a plain) in double precision
+! SYNTAX:
+! real= dc_area(real.arr, real.arr, real.arr)
+! ON INPUT:
+! r_coord1: node1 of triangle REAL
+! r_coord2: node2 of triangle REAL
+! r_coord3: node3 of triangle REAL
+! ON OUTPUT:
+! r_area: area of triangle REAL
+! CALLS:
+!
+! COMMENTS:
+!
+!-----------------------------------------------------------------
+!
+! PUBLIC:
+! slm_initialvalues, slm_analyticsolution
+! COMMENTS:
+!
+! USES:
+! MISC_globalparam, MISC_error, GRID_api
+! LIBRARIES:
+!
+! REFERENCES:
+!
+! VERSION(S):
+! 1. original version j. behrens 7/97
+! 2. names changed j. behrens 7/97
+! 3. changed to use GRID_api j. behrens 11/97
+! 4. changed interfaces j. behrens 12/97
+! 5. compliant to amatos 1.0 j. behrens 12/2000
+! 6. compliant to amatos 1.2 j. behrens 3/2002
+! 7. compliant to amatos 2.0 f. klaschka 8/2006
+!
+!*****************************************************************
+ MODULE SLM_initial
+ USE FLASH_parameters
+ USE GRID_api
+ PRIVATE
+ PUBLIC slm_initialvalues, slm_analyticsolution
+! REAL, DIMENSION(GRID_dimension) :: r_cntr=(/ -0.25, 0.0 /)
+! REAL :: r_hgt=4.0
+! REAL :: r_srd=0.15
+! REAL :: r_sln=0.22
+! REAL :: r_swd=0.06
+ CONTAINS
+!*****************************************************************
+ SUBROUTINE slm_initialvalues(p_ghand)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ TYPE (grid_handle), INTENT(in) :: p_ghand
+ INTEGER :: i_lev= 6
+
+!---------- initialize some constant for the slotted cylinder
+
+ CALL slm_initcoshill(p_ghand)
+
+ RETURN
+ END SUBROUTINE slm_initialvalues
+!*****************************************************************
+ SUBROUTINE slm_analyticsolution(p_ghand, r_time, i_arlen, r_array)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ TYPE (grid_handle), INTENT(in) :: p_ghand
+ REAL (KIND = GRID_SR), INTENT(in) :: r_time
+ INTEGER, INTENT(in) :: i_arlen
+ REAL (KIND = GRID_SR), DIMENSION(i_arlen), INTENT(out) :: r_array
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_centr
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimspherical) :: r_sentr=(/ 0.0, 0.0 /)
+ REAL (KIND = GRID_SR) :: r_rds, r_dpt
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_tmp
+ INTEGER :: i_count, i_num, i_alct
+ REAL (KIND = GRID_SR), DIMENSION(:,:), ALLOCATABLE :: r_coo
+
+!---------- calculate center ! CAUTION: THIS IS NOT IMPLEMENTED YET !!!
+
+ r_centr= grid_geokart(r_sentr)
+
+!---------- allocate workspace
+
+ i_num= p_ghand%i_nnumber
+ ALLOCATE(r_coo(GRID_dimension,i_num), stat= i_alct)
+ IF(i_alct /= 0) THEN
+ CALL grid_error(55)
+ END IF
+
+!---------- get information
+
+ CALL grid_getinfo(p_ghand, r_nodecoordinates= r_coo)
+
+!---------- check array length
+
+ IF(i_arlen < i_num) THEN
+ IF(GRID_parameters%iolog > 0) THEN
+ write(GRID_parameters%iolog,*) '[slm_analyticsolution]: Array sizes did not match, adjusting to the smaller size!'
+ END IF
+ i_num= i_arlen
+ END IF
+
+!---------- loop over the nodes
+
+ node_loop: DO i_count= 1, i_num
+ r_array(i_count)= coshill(r_coo(:,i_count),r_centr)
+ END DO node_loop
+
+!---------- deallocate workspace
+
+ DEALLOCATE(r_coo)
+
+ RETURN
+ END SUBROUTINE slm_analyticsolution
+
+!*****************************************************************
+ SUBROUTINE slm_initcoshill(p_ghand)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ TYPE (grid_handle), INTENT(in) :: p_ghand
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_centr
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimspherical) :: r_sentr=(/ 0.0, 0.0 /)
+ REAL (KIND = GRID_SR) :: r_rds, r_dpt
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_tmp
+ INTEGER (KIND = GRID_SI) :: i_count, i_num, i_alct, &
+ i_lay, i_ind
+ REAL (KIND = GRID_SR), DIMENSION(:,:), ALLOCATABLE :: r_aux
+ REAL (KIND = GRID_SR), DIMENSION(:,:), ALLOCATABLE :: r_coo
+ INTEGER (KIND = GRID_SI), DIMENSION(1) :: i_valind
+
+!---------- get center
+
+ r_centr= grid_geokart(r_sentr)
+
+!---------- allocate workspace
+
+ i_num= p_ghand%i_nnumber
+
+ ALLOCATE(r_aux(1,i_num), r_coo(GRID_dimension,i_num), stat= i_alct)
+ IF(i_alct /= 0) THEN
+ CALL grid_error(55)
+ END IF
+
+!---------- get information
+
+ CALL grid_getinfo(p_ghand, r_nodecoordinates= r_coo)
+
+!---------- loop over layers
+
+ layer_loop: DO i_lay= 1, i_stcnumlayers
+
+!---------- compute radii of initial ball-type cosine hill density...
+
+! i_ind = max(4 - abs(floor(0.5*i_stcnumlayers)-i_lay),0)
+! r_rds = (GRID_RADIUS/ 6.0)*i_ind/4
+
+ r_aux = 0.
+ IF ((i_lay .GE. 2) .AND. (i_lay .LE. i_stcnumlayers -2)) then
+!---------- loop over the nodes
+
+ node_loop: DO i_count= 1, i_num
+ r_aux(1,i_count)= coshill(r_coo(:,i_count),r_centr, r_radius=r_rds)
+ END DO node_loop
+ END IF
+
+!---------- update grid information
+
+ i_valind= (/ i_stctracer-1+i_lay /)
+ CALL grid_putinfo(p_ghand, i_arraypoint=i_valind, r_nodevalues= r_aux)
+ END DO layer_loop
+
+!---------- deallocate workspace
+
+ DEALLOCATE(r_aux, r_coo)
+
+ RETURN
+ END SUBROUTINE slm_initcoshill
+
+!*****************************************************************
+ FUNCTION coshill(r_coor, r_centr, r_radius) RESULT (r_hill)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_coor
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_centr
+ REAL (KIND = GRID_SR), OPTIONAL :: r_radius
+ REAL (KIND = GRID_SR) :: r_hill
+ REAL (KIND = GRID_SR) :: r_maxheight=1.0
+ REAL (KIND = GRID_SR) :: r_maxrad
+ REAL (KIND = GRID_SR) :: r_dist, r_tmp
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_xy
+
+!---------- watch out for optional argument
+
+ opt_arg: IF(present(r_radius)) THEN
+ r_maxrad = r_radius
+ ELSE opt_arg
+!---------- initialize maximal radius with default values
+ r_maxrad= GRID_RADIUS/ 6.0
+ END IF opt_arg
+
+!---------- initialize r_hill
+
+ r_hill= 0.0
+
+!---------- calculate distance to center
+
+ r_xy= r_centr- r_coor
+ r_dist= DOT_PRODUCT(r_xy,r_xy)
+
+!---------- calculate inner values
+
+ IF(r_dist < r_maxrad) THEN
+ r_tmp = (GRID_PI* r_dist)/r_maxrad
+ r_hill= r_maxheight *(1.0+ cos(r_tmp))
+ END IF
+
+ RETURN
+ END FUNCTION coshill
+
+!*****************************************************************
+
+ END MODULE SLM_initial
+
+
+
+
diff --git a/flash2d/src/options-sphere/SLM_simple_stacked.f90 b/flash2d/src/options-sphere/SLM_simple_stacked.f90
new file mode 100644
index 0000000000000000000000000000000000000000..4c42f69631ed3fa20badd5039a04db16040c2413
--- /dev/null
+++ b/flash2d/src/options-sphere/SLM_simple_stacked.f90
@@ -0,0 +1,539 @@
+!*****************************************************************
+!
+! MODULE NAME:
+! SLM_simple
+! FUNCTION:
+! provide simple semi-Lagrangian routines
+! CONTAINS:
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_step
+! FUNCTION:
+! one step of the basic SLM algorithm
+! SYNTAX:
+! CALL slm_step(int, real.arr, real.arr)
+! ON INPUT:
+! ...
+! ON OUTPUT:
+! r_tracer: array with tracer values real
+! CALLS:
+!
+! COMMENTS:
+!
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_displace
+! FUNCTION:
+! extrapolate the alpha, values for the displacements of the upstream
+! points from the gridpoints
+! SYNTAX:
+! CALL slm_displace(int, real.arr, real.arr)
+! ON INPUT:
+! i_arlen: array length for the real arrays integer
+! r_coord: real array of xy-coordinates real
+! ON OUTPUT:
+! r_alpha: displacement vectors to each point real
+! CALLS:
+! wind_field
+! COMMENTS:
+!
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_update
+! FUNCTION:
+! calculate the update to the velocity
+! SYNTAX:
+! CALL slm_update(int, real.arr, real.arr)
+! ON INPUT:
+! i_arlen: array length for the real arrays integer
+! r_rside: array with right hand side values real
+! ON OUTPUT:
+! r_udate: array with new (updated) gid values real
+! CALLS:
+!
+! COMMENTS:
+! this routine is trivial for linear advection
+!
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_upstream
+! FUNCTION:
+! calculate right hand side of the equation (upstream values)
+! SYNTAX:
+! CALL slm_upstream(int, real.arr, real.arr)
+! ON INPUT:
+! i_arlen: array length for the real arrays integer
+! r_alpha: displacement vectors to each point real
+! ON OUTPUT:
+! r_rside: array with right hand side values real
+! CALLS:
+!
+! COMMENTS:
+! this routine is just interpolation for linear advection
+!
+!-----------------------------------------------------------------
+!
+! NAME:
+! slm_interpolate
+! FUNCTION:
+! do the interpolation
+! SYNTAX:
+! CALL slm_interpolate(grid, int, real, real.arr, real.arr, real.arr)
+! ON INPUT:
+! p_ogrid: grid handle to old grid (with data) TYPE (grid_handle)
+! r_fac: factor at which point to interpolate REAL
+! i_arlen: array length for the following arrays INTEGER
+! r_coord: coordinate array (new grid) REAL
+! r_alpha: displacement array (corr. to r_coord) REAL
+! r_value: values on the old grid (array) REAL
+! ON OUTPUT:
+! r_rside: right hand side (interpolated) REAL
+! CALLS:
+!
+! COMMENTS:
+! this one is plain bi-cubic spline interpolation
+!
+!-----------------------------------------------------------------
+!
+! PUBLIC:
+! slm_step
+! COMMENTS:
+!
+! USES:
+! FLASH_parameters, GRID_api, ADV_wind, ADV_rhs
+! LIBRARIES:
+!
+! REFERENCES:
+!
+! VERSION(S):
+! 1. original version j. behrens 4/2002
+! 2. compliant to amatos 2.0 j. behrens 7/2003
+!
+!*****************************************************************
+ MODULE SLM_simple
+ USE FLASH_parameters
+ USE MISC_timing
+ USE GRID_api
+ USE ADV_wind
+ USE ADV_rhs
+ PUBLIC :: slm_step, slm_verticalstep
+ PRIVATE
+ CONTAINS
+
+!*****************************************************************
+ SUBROUTINE slm_step(p_ghand, p_param, p_time, r_modtime, i_size, &
+ r_coord, r_tracer, i_newsdepth, i_layer)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ TYPE (grid_handle), DIMENSION(GRID_timesteps), INTENT(in) :: p_ghand
+ TYPE (control_struct), INTENT(in) :: p_param
+ TYPE (sw_info), INTENT(inout) :: p_time
+ REAL (KIND = GRID_SR), INTENT(in) :: r_modtime
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_size
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension,i_size), INTENT(in) :: r_coord
+ REAL (KIND = GRID_SR), DIMENSION(i_size), INTENT(out) :: r_tracer
+ INTEGER (KIND = GRID_SI), OPTIONAL, INTENT(in) :: i_newsdepth
+ INTEGER (KIND = GRID_SI), OPTIONAL, INTENT(in) :: i_layer
+
+ REAL (KIND = GRID_SR), DIMENSION(:), ALLOCATABLE :: r_newvl
+ REAL (KIND = GRID_SR), DIMENSION(:,:), ALLOCATABLE :: r_alpha
+ INTEGER (KIND = GRID_SI) :: i_alct, i_layerid, i_nd
+
+!---------- evaluate layer information, if given
+
+ layer_info: IF(present(i_layer)) THEN
+ i_layerid = i_layer
+ ELSE layer_info
+ i_layerid = 1 ! Default is the lower layer (layer 1)
+ END IF layer_info
+
+ IF(present(i_newsdepth)) THEN
+ i_nd= i_newsdepth
+ ELSE
+ i_nd= 1_GRID_SI
+ END IF
+
+!---------- check size!
+
+ IF(i_size <= 0) THEN
+ IF(GRID_parameters%iolog > 0) &
+ write(GRID_parameters%iolog,*) 'INFO [slm_step]: Zero step size, returning to calling routine'
+ RETURN
+ END IF
+
+!---------- allocate auxiliary arrays
+
+ allocate(r_newvl(i_size), r_alpha(GRID_dimension,i_size), stat=i_alct)
+ not_alloc: IF(i_alct /= 0) THEN
+ CALL grid_error(38)
+ END IF not_alloc
+
+!-SLM--------- calculate trajectory pieces (displacements)
+
+ CALL stop_watch('start',3,p_time)
+ CALL slm_displace(p_param, i_size, i_layerid, r_coord, r_alpha, r_time=r_modtime)
+ CALL stop_watch('stop ',3,p_time)
+
+!-SLM--------- calculate right hand side
+
+ CALL stop_watch('start',4,p_time)
+ CALL slm_upstream(p_ghand, i_size, i_layerid, r_coord, r_alpha, r_newvl)
+ CALL stop_watch('stop ',4,p_time)
+
+!-SLM--------- calculate new grid values
+
+ CALL stop_watch('start',5,p_time)
+ CALL slm_update(p_param, i_size, i_layerid, r_coord, r_newvl, r_tracer, r_time=r_modtime)
+ CALL stop_watch('stop ',5,p_time)
+
+!-SLM--------- put alpha values to u and v field entries
+
+ IF(i_layerid == 1) THEN
+ r_alpha= -r_alpha
+ IF(present(i_newsdepth)) THEN
+ CALL grid_putinfo(p_ghand(i_timeplus), r_nodevalues=r_alpha(1:2,:), &
+ i_newsdepth=i_newsdepth, i_arraypoint=(/ GRID_ucomp, GRID_vcomp /))
+ ELSE
+ CALL grid_putinfo(p_ghand(i_timeplus), r_nodevalues=r_alpha(1:2,:), &
+ i_arraypoint=(/ GRID_ucomp, GRID_vcomp /))
+ END IF
+ END IF
+
+!-SLM--------- deallocate work arrays
+
+ deallocate(r_alpha, r_newvl)
+
+ RETURN
+ END SUBROUTINE slm_step
+
+!*****************************************************************
+ SUBROUTINE slm_verticalstep(p_ghand, p_param, p_time, r_modtime, i_size, &
+ r_coord, r_tracer, i_newsdepth)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ TYPE (grid_handle), DIMENSION(GRID_timesteps), INTENT(in) :: p_ghand
+ TYPE (control_struct), INTENT(in) :: p_param
+ TYPE (sw_info), INTENT(inout) :: p_time
+ REAL (KIND = GRID_SR), INTENT(in) :: r_modtime
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_size
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension,i_size), INTENT(in) :: r_coord
+ REAL (KIND = GRID_SR), DIMENSION(:,:), INTENT(inout) :: r_tracer
+ INTEGER (KIND = GRID_SI), OPTIONAL, INTENT(in) :: i_newsdepth
+
+ REAL (KIND = GRID_SR), DIMENSION(:), ALLOCATABLE :: r_newvl
+ REAL (KIND = GRID_SR), DIMENSION(:,:), ALLOCATABLE :: r_alpha
+ INTEGER (KIND = GRID_SI) :: i_alct, i_layerid
+
+ REAL (KIND = GRID_SR) :: r_upstr, r_dd
+ INTEGER (KIND = GRID_SI) :: i_cnt, i_lyr, &
+ i_lcnt, i_nd
+
+!---------- allocate intermediate storage for the new tracer values
+
+ IF(present(i_newsdepth)) THEN
+ i_nd= i_newsdepth
+ ELSE
+ i_nd= 1_GRID_SI
+ END IF
+
+ allocate(r_newvl(i_stcnumlayers), stat=i_alct)
+ not_alloc: IF(i_alct /= 0) THEN
+ CALL grid_error(c_error='[slm_verticalstep]: could not allocate r_newvl')
+ END IF not_alloc
+
+!---------- outer loop over all coordinates
+ coord_loop: DO i_cnt=1,i_size
+
+!---------- inner loop over height levels
+ layer_loop: DO i_lyr=1,i_stcnumlayers
+
+!---------- an explicit Euler step for the upwind position
+ r_upstr = r_stclayerpressure(i_lyr) - p_param%num%r_deltatime* &
+ slm_vertwind(r_coord(:,i_cnt),i_lyr)
+
+!---------- upstream value interpolation (linear)
+! Step 1: check for boundary
+! upstream position is higher than highest layer
+ bound_check: IF(r_upstr <= r_stclayerpressure(i_stcnumlayers)) THEN
+ r_newvl(i_lyr) = 0._GRID_SR
+! upstream position is lower than ground layer
+ ELSE IF(r_upstr > r_stclayerpressure(1)) THEN bound_check
+ r_newvl(i_lyr) = 0._GRID_SR
+! upstream position is within vertical domain
+ ELSE bound_check
+! Step 2: interpolate upstream position - use linear interpolation
+! find vertical box of upstream position
+ layer_search: DO i_lcnt=1,i_stcnumlayers-1
+ IF((r_upstr <= r_stclayerpressure(i_lcnt)) .AND. &
+ (r_upstr > r_stclayerpressure(i_lcnt+1))) THEN
+ i_layerid = i_lcnt
+ EXIT layer_search
+ END IF
+ END DO layer_search
+! now the linear interpolation
+ r_dd = (r_tracer(i_layerid+1, i_cnt)- r_tracer(i_layerid, i_cnt))/ &
+ (r_stclayerpressure(i_layerid+1)- r_stclayerpressure(i_layerid))
+ r_newvl(i_lyr) = r_tracer(i_layerid,i_cnt) + r_dd* &
+ (r_upstr - r_stclayerpressure(i_layerid))
+ END IF bound_check
+
+!---------- finally, the update with a possible right hand side
+ r_newvl(i_lyr) = r_newvl(i_lyr) + p_param%num%r_deltatime* &
+ slm_verticalrhs(r_coord(:,i_cnt),i_lyr,r_modtime)
+ END DO layer_loop
+
+!---------- overwrite the tracer values with the newly computed advected values in this column
+ r_tracer(:,i_cnt) = r_newvl(:)
+
+ END DO coord_loop
+
+!---------- deallocate intermediate storage
+
+ deallocate(r_newvl)
+
+ RETURN
+ END SUBROUTINE slm_verticalstep
+
+!*****************************************************************
+ SUBROUTINE slm_displace(p_param, i_arlen, i_layerid, r_coord, r_alpha, r_time)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ TYPE (control_struct), INTENT(in) :: p_param
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_arlen
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_layerid
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension,i_arlen), INTENT(in) :: r_coord
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension,i_arlen), INTENT(out) :: r_alpha
+ REAL (KIND = GRID_SR), INTENT(in), OPTIONAL :: r_time
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_fac, r_caf, &
+ r_axy, r_xyc, r_sxy
+ REAL (KIND = GRID_SR) :: r_dt0, r_dt1, &
+ r_dt2, r_tim
+ INTEGER (KIND = GRID_SI) :: i_cnt1, i_cnt2
+
+!---------- set constants
+
+ r_dt0= p_param%num%r_deltatime
+ r_dt1= 0.5* p_param%num%r_deltatime
+ r_dt2= 1.5* p_param%num%r_deltatime
+ r_fac= 0.5
+ r_caf= 2.0
+ IF(present(r_time)) THEN
+ r_tim= r_time
+ ELSE
+ r_tim= 0.0
+ END IF
+
+!---------- calculate in an iteration process the displacements
+
+ unknown_loop: DO i_cnt1=1,i_arlen
+ r_axy= 0.0
+
+ iter_loop: DO i_cnt2=1, p_param%num%i_adviterations
+ r_xyc= r_coord(:,i_cnt1)- r_fac* r_axy
+ r_sxy= r_dt0* slm_windfield(r_xyc, i_layerid, r_time=r_tim)
+ r_axy= sphere_correct(r_coord(:,i_cnt1),r_sxy)
+ END DO iter_loop
+
+ r_alpha(:,i_cnt1)= r_axy
+ END DO unknown_loop
+
+ RETURN
+ END SUBROUTINE slm_displace
+
+!*****************************************************************
+ SUBROUTINE slm_update(p_param, i_arlen, i_layer, r_coord, r_rside, r_udate, r_time)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ TYPE (control_struct), INTENT(in) :: p_param
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_arlen
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_layer
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension,i_arlen), INTENT(in) :: r_coord
+ REAL (KIND = GRID_SR), DIMENSION(i_arlen), INTENT(in) :: r_rside
+ REAL (KIND = GRID_SR), DIMENSION(i_arlen), INTENT(out) :: r_udate
+ REAL (KIND = GRID_SR), INTENT(in), OPTIONAL :: r_time
+
+ INTEGER (KIND = GRID_SI) :: i_cnt
+ REAL (KIND = GRID_SR) :: r_dt, r_tim
+
+!---------- in the linear advection case and with f90 this is just
+
+! r_udate= r_rside
+
+!---------- including a non-zero right hand side, we have
+
+ r_dt= p_param%num%r_deltatime
+ IF(present(r_time)) THEN
+ r_tim= r_time
+ ELSE
+ r_tim= 0.0
+ END IF
+
+ main_loop: DO i_cnt=1, i_arlen
+ r_udate(i_cnt)= r_rside(i_cnt)+ r_dt* slm_righthand(r_coord(:,i_cnt), i_layer)
+ END DO main_loop
+
+ RETURN
+ END SUBROUTINE slm_update
+
+!*****************************************************************
+ SUBROUTINE slm_upstream(p_mesh, i_arlen, i_layer, r_coord, &
+ r_alpha, r_rside)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ TYPE (grid_handle), DIMENSION(GRID_timesteps) :: p_mesh
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_arlen
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_layer
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension,i_arlen), INTENT(in) :: r_coord
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension,i_arlen), INTENT(in) :: r_alpha
+ REAL (KIND = GRID_SR), DIMENSION(i_arlen), INTENT(out) :: r_rside
+
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_fac
+
+!---------- set factor (at which point of trajectory shall i interpolate)
+
+ r_fac= 1.0
+
+!---------- in the linear advection case this is just interpolation
+
+ CALL slm_interpolate(p_mesh, r_fac, i_arlen, i_layer, r_coord, &
+ r_alpha, r_rside)
+
+ RETURN
+ END SUBROUTINE slm_upstream
+!*****************************************************************
+ SUBROUTINE slm_interpolate(p_mesh, r_fac, i_arlen, i_layer, &
+ r_coord, r_alpha, r_rside)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ TYPE (grid_handle), DIMENSION(GRID_timesteps) :: p_mesh
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension), INTENT(in) :: r_fac
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_arlen
+ INTEGER (KIND = GRID_SI), INTENT(in) :: i_layer
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension,i_arlen), INTENT(in) :: r_coord
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension,i_arlen), INTENT(in) :: r_alpha
+ REAL (KIND = GRID_SR), DIMENSION(i_arlen), INTENT(out) :: r_rside
+
+ REAL (KIND = GRID_SR), DIMENSION(:,:), ALLOCATABLE :: r_upstr
+ REAL (KIND = GRID_SR) :: r_eps
+ INTEGER (KIND = GRID_SI) :: i_cnt, &
+ i_alct, i_val, i_out, i_stat
+
+!---------- initialize constant
+
+ i_val= i_stctracer-1+i_layer
+ r_eps= GRID_EPS
+
+!---------- allocate work array
+
+ ALLOCATE(r_upstr(GRID_dimension,i_arlen), stat=i_alct)
+ not_allocated: IF(i_alct /= 0) THEN
+ CALL grid_error(60)
+ END IF not_allocated
+
+!---------- calculate upstream coordinates
+
+ dim_loop: DO i_cnt=1, GRID_dimension
+ r_upstr(i_cnt,:) = r_coord(i_cnt,:)- r_fac(i_cnt)* r_alpha(i_cnt,:)
+ END DO dim_loop
+
+!---------- loop over nodes: find element containing upstream point
+
+ node_loop: DO i_cnt=1, i_arlen
+
+!---------- check if upstream value is outside of the domain
+
+ i_out= grid_domaincheck(p_mesh(i_time), r_upstr(:,i_cnt))
+
+!---------- take the intersection of the trajectory with the boundary as new upstream point
+
+ out_domain: IF(i_out /= 0) then
+ r_upstr(:,i_cnt)= grid_boundintersect(p_mesh(i_time), &
+ r_coord(:,i_cnt), r_upstr(:,i_cnt), i_info=i_stat)
+ no_intersect: IF(i_stat /= 0) THEN
+ r_rside(i_cnt)= 0.0
+ CYCLE node_loop
+ END IF no_intersect
+ END IF out_domain
+
+!---------- interpolate
+
+ r_rside(i_cnt)= grid_coordvalue(p_mesh(i_time), r_upstr(:,i_cnt), &
+ i_interpolorder=GRID_loworder, i_valpoint=i_val)
+ small_val: IF(abs(r_rside(i_cnt)) < r_eps) THEN
+ r_rside(i_cnt)= 0.0
+ END IF small_val
+
+ END DO node_loop
+
+!---------- deallocate work array
+
+ DEALLOCATE(r_upstr)
+
+ RETURN
+ END SUBROUTINE slm_interpolate
+
+!*****************************************************************
+ FUNCTION sphere_correct(r_coord, r_displ) RESULT (r_corct)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_coord, r_displ
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_corct
+
+ REAL (KIND = GRID_SR) :: r_e, r_rat
+ REAL (KIND = GRID_SR) :: r_c
+
+!---------- calculate Euklidean norm
+
+ r_e = eukl_norm(r_displ)
+ r_rat= r_e/ GRID_RADIUS
+
+!---------- calculate correction
+
+ r_c = tan(r_rat)/r_rat
+ r_corct= r_displ* r_c
+
+ RETURN
+ END FUNCTION sphere_correct
+
+!*****************************************************************
+ FUNCTION eukl_norm(r_vec) RESULT(r_rst)
+
+!---------- local declarations
+
+ IMPLICIT NONE
+ REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_vec
+ REAL (KIND = GRID_SR) :: r_rst, r_tmp
+
+!---------- calculate vector crossproduct, 3D only
+
+ r_tmp= dot_product(r_vec, r_vec)
+ r_rst= sqrt(r_tmp)
+
+ END FUNCTION eukl_norm
+
+ END MODULE SLM_simple