diff --git a/src/2d/variables_2d.F90 b/src/2d/variables_2d.F90
index 238078c72aa829cce16da01462752498ce150d7a..c0ab8e64d1a6b4319d70f90207d5bd654fa7856f 100644
--- a/src/2d/variables_2d.F90
+++ b/src/2d/variables_2d.F90
@@ -18,6 +18,7 @@
 !
 ! !USES:
    use domain, only: imin,imax,jmin,jmax
+   use field_manager
    IMPLICIT NONE
 !
 ! !PUBLIC DATA MEMBERS:
@@ -147,6 +148,81 @@
    end subroutine init_variables_2d
 !EOC
 
+!-----------------------------------------------------------------------
+!BOP
+!
+! !ROUTINE: register_2d_variables() - register GETM variables.
+!
+! !INTERFACE:
+   subroutine register_2d_variables(fm)
+!
+! !DESCRIPTION:
+!
+! !USES:
+!   use variables_2d
+   IMPLICIT NONE
+!
+! !INPUT PARAMETERS:
+   type (type_field_manager) :: fm
+!
+! !REVISION HISTORY:
+!  Original author(s): Karsten Bolding & Jorn Bruggeman
+!
+! !LOCAL VARIABLES:
+   logical :: used
+   integer,parameter :: rk = kind(_ONE_)
+!EOP
+!-----------------------------------------------------------------------
+!BOC
+   LEVEL2 'register_2d_variables()'
+
+!D(E2DFIELD)
+!DU,DV
+!z(E2DFIELD)
+!zo(E2DFIELD)
+!U(E2DFIELD)
+!V(E2DFIELD)
+!UEx(E2DFIELD)
+!VEx(E2DFIELD)
+!fU(E2DFIELD)
+!fV(E2DFIELD)
+!ru(E2DFIELD)
+!rv(E2DFIELD)
+!Uint(E2DFIELD)
+!Vint(E2DFIELD)
+!Uinto(E2DFIELD)
+!Vinto(E2DFIELD)
+!res_du(E2DFIELD)
+!res_u(E2DFIELD)
+!res_dv(E2DFIELD)
+!res_v(E2DFIELD)
+!kbk
+!SlUx(E2DFIELD)
+!SlVx(E2DFIELD)
+!Slru(E2DFIELD)
+!Slrv(E2DFIELD)
+!zub(E2DFIELD)
+!zvb(E2DFIELD)
+!zub0(E2DFIELD)
+!zvb0(E2DFIELD)
+!An(E2DFIELD)
+!AnX(E2DFIELD)
+!fwf(E2DFIELD)
+!fwf_int(E2DFIELD)
+!EWbdy(jmax),ENbdy(imax),EEbdy(jmax),ESbdy(imax)
+
+
+!  category - 2d
+   call fm%register('z', 'm', 'sea surface elevation', standard_name='sea surface elevation', fill_value=-9999.0_rk, data2d=z(_2D_W_), category="2d", part_of_state=.true.)
+   call fm%register('zo', 'm', 'sea surface elevation', standard_name='sea surface elevation', fill_value=-9999.0_rk, data2d=zo(_2D_W_), category="2d", output_level=output_level_debug)
+   call fm%register('D', 'm', 'water depth', standard_name='water depth', fill_value=-10009.0_rk, data2d=D(_2D_W_), category="2d", part_of_state=.true.)
+   call fm%register('U', 'm2/s', 'transport in local x-direction', standard_name='', data2d=U(_2D_W_), category='2d', output_level=output_level_debug, part_of_state=.true.)
+   call fm%register('V', 'm2/s', 'transport in local y-direction', standard_name='', data2d=V(_2D_W_), category='2d', output_level=output_level_debug, part_of_state=.true.)
+
+   return
+   end subroutine register_2d_variables
+!EOC
+
 !-----------------------------------------------------------------------
 !BOP
 !
diff --git a/src/3d/coordinates.F90 b/src/3d/coordinates.F90
index 31f7d249d1ee44d3d1e4ae31383e3ad6c988ffc1..3b189183c6d479ecb0f9fa2249a0a30592a4ba05 100644
--- a/src/3d/coordinates.F90
+++ b/src/3d/coordinates.F90
@@ -30,8 +30,8 @@
 !
 !
 ! !USES:
+   use domain, only: imin,imax,jmin,jmax,kmax,H
 #ifdef SLICE_MODEL
-   use domain, only: imin,imax,jmin,jmax,kmax
    use variables_3d, only: kvmin,hvo,hvn
 #endif
    use getm_timers, only: tic, toc,TIM_COORDS
@@ -52,9 +52,7 @@
    logical, save   :: first=.true.
    integer         :: ii
 !   integer         :: preadapt=0
-#ifdef SLICE_MODEL
    integer          :: i,j,k
-#endif
 !EOP
 !-----------------------------------------------------------------------
 !BOC
@@ -100,6 +98,14 @@ stop
       end select
    end if ! first
 
+   ! calculate the z-coordinate of the cell centers
+   ! references to mean sea level
+   zc(:,:,0)=-H(:,:)
+   zc(:,:,1)=-H(:,:) + 0.5*hn(:,:,1)
+   do k=2,kmax
+      zc(:,:,k)=zc(:,:,k-1)+0.5*(hn(:,:,k-1)+hn(:,:,k))
+   end do
+
 #ifdef SLICE_MODEL
    do i=imin,imax
       do k=kvmin(i,2),kmax
diff --git a/src/3d/dynamic_allocations_3d.h b/src/3d/dynamic_allocations_3d.h
index 7dee1aee100d9d0f726628150820c14f7410025f..4d7f07f8dcb009a139d0982f27ed3bee9215c057 100644
--- a/src/3d/dynamic_allocations_3d.h
+++ b/src/3d/dynamic_allocations_3d.h
@@ -58,6 +58,9 @@
   if (rc /= 0) stop 'init_3d: Error allocating memory (sf)'
 #endif
 
+  allocate(zc(I3DFIELD),stat=rc)    ! 3D field for cell center coordinates
+  if (rc /= 0) stop 'init_3d: Error allocating memory (zc)'
+
   allocate(ho(I3DFIELD),stat=rc)    ! 3D field for old box height (z-column)
   if (rc /= 0) stop 'init_3d: Error allocating memory (ho)'
 
diff --git a/src/3d/dynamic_declarations_3d.h b/src/3d/dynamic_declarations_3d.h
index eaca9cfaa15a27836bd69e3ba9d858ffd000f036..d24e00ebf28e9a68924385b3fc375ee2b67ad8be 100644
--- a/src/3d/dynamic_declarations_3d.h
+++ b/src/3d/dynamic_declarations_3d.h
@@ -21,6 +21,7 @@
 #ifdef STRUCTURE_FRICTION
   REALTYPE, dimension(:,:,:), allocatable   :: sf
 #endif
+  REALTYPE, dimension(:,:,:), allocatable, target :: zc
   REALTYPE, dimension(:,:,:), allocatable, target :: hn,hun,hvn
   REALTYPE, dimension(:,:,:), allocatable   :: ho,huo,hvo
   REALTYPE, dimension(:,:,:), allocatable   :: hcc
diff --git a/src/3d/getm_fabm.F90 b/src/3d/getm_fabm.F90
index 3fcf68529930b77a383ab8147bab429583c10e53..cde81bfecb24c067417581c0e278ae480fc79f4e 100644
--- a/src/3d/getm_fabm.F90
+++ b/src/3d/getm_fabm.F90
@@ -39,6 +39,7 @@
 ! !PUBLIC DATA MEMBERS:
    public init_getm_fabm_fields
    public init_getm_fabm, do_getm_fabm, model, output_none
+   public register_fabm_variables
    integer, public :: fabm_init_method=0
    character(len=PATH_MAX)   :: fabm_init_file
    integer                   :: fabm_init_format, fabm_field_no
@@ -318,6 +319,73 @@ end interface
    end subroutine init_getm_fabm_fields
 !EOC
 
+#ifdef _FABM_
+!-----------------------------------------------------------------------
+!BOP
+!
+! !ROUTINE: register_fabm_variables() - register FABM variables.
+!
+! !INTERFACE:
+   subroutine register_fabm_variables(fm)
+!
+! !DESCRIPTION:
+!
+! !USES:
+   use field_manager
+   IMPLICIT NONE
+!
+! !INPUT PARAMETERS:
+   type (type_field_manager) :: fm
+!
+! !REVISION HISTORY:
+!  Original author(s): Karsten Bolding & Jorn Bruggeman
+!
+! !LOCAL VARIABLES:
+  integer :: i,output_level
+  logical :: in_output
+!EOP
+!-----------------------------------------------------------------------
+!BOC
+   if (.not. fabm_calc) return
+   LEVEL2 'register_fabm_variables()'
+
+   do i=1,size(model%state_variables)
+      output_level = output_level_default
+      if (model%state_variables(i)%output==output_none) output_level = output_level_debug
+      call fm%register(model%state_variables(i)%name, model%state_variables(i)%units, &
+         model%state_variables(i)%long_name, minimum=model%state_variables(i)%minimum, maximum=model%state_variables(i)%maximum, &
+         fill_value=model%state_variables(i)%missing_value, dimensions=(/id_dim_z/), data3d=fabm_pel(_3D_W_,i), category='fabm'//model%state_variables(i)%target%owner%get_path(), output_level=output_level)
+   end do
+   do i=1,size(model%bottom_state_variables)
+      output_level = output_level_default
+      if (model%bottom_state_variables(i)%output==output_none) output_level = output_level_debug
+      call fm%register(model%bottom_state_variables(i)%name, model%bottom_state_variables(i)%units, &
+         model%bottom_state_variables(i)%long_name, minimum=model%bottom_state_variables(i)%minimum, &
+         maximum=model%bottom_state_variables(i)%maximum, fill_value=model%bottom_state_variables(i)%missing_value, &
+         data2d=fabm_ben(_2D_W_,i), category='fabm'//model%bottom_state_variables(i)%target%owner%get_path(), output_level=output_level)
+   end do
+   do i=1,size(model%diagnostic_variables)
+      output_level = output_level_default
+      if (model%diagnostic_variables(i)%output==output_none) output_level = output_level_debug
+      call fm%register(model%diagnostic_variables(i)%name, model%diagnostic_variables(i)%units, &
+         model%diagnostic_variables(i)%long_name, minimum=model%diagnostic_variables(i)%minimum, maximum=model%diagnostic_variables(i)%maximum, &
+         fill_value=model%diagnostic_variables(i)%missing_value, dimensions=(/id_dim_z/), data3d=fabm_diag(_3D_W_,i), category='fabm'//model%diagnostic_variables(i)%target%owner%get_path(), output_level=output_level, used=in_output)
+      if (in_output) model%diagnostic_variables(i)%save = .true.
+   end do
+   do i=1,size(model%horizontal_diagnostic_variables)
+      output_level = output_level_default
+      if (model%horizontal_diagnostic_variables(i)%output==output_none) output_level = output_level_debug
+      call fm%register(model%horizontal_diagnostic_variables(i)%name, model%horizontal_diagnostic_variables(i)%units, &
+         model%horizontal_diagnostic_variables(i)%long_name, minimum=model%horizontal_diagnostic_variables(i)%minimum, maximum=model%horizontal_diagnostic_variables(i)%maximum, &
+         fill_value=model%horizontal_diagnostic_variables(i)%missing_value, data2d=fabm_diag_hz(_2D_W_,i), category='fabm'//model%horizontal_diagnostic_variables(i)%target%owner%get_path(), output_level=output_level, used=in_output)
+      if (in_output) model%horizontal_diagnostic_variables(i)%save = .true.
+   end do
+
+   return
+   end subroutine register_fabm_variables
+!EOC
+#endif
+
 !-----------------------------------------------------------------------
 !BOP
 !
diff --git a/src/3d/static_3d.h b/src/3d/static_3d.h
index 42bf6c0459459b1b46ef3803636fbed1d60346ed..88fd8f7ef835ef328e24b39c459ead6d66498491 100644
--- a/src/3d/static_3d.h
+++ b/src/3d/static_3d.h
@@ -35,6 +35,7 @@
 #ifdef STRUCTURE_FRICTION
   REALTYPE                             :: sf(I3DFIELD)
 #endif
+  REALTYPE                             :: zc(I3DFIELD)
   REALTYPE                             :: ho(I3DFIELD)
   REALTYPE, target                     :: hn(I3DFIELD)
   REALTYPE                             :: huo(I3DFIELD)
diff --git a/src/3d/variables_3d.F90 b/src/3d/variables_3d.F90
index af077c33a7345c8027fdb6098c35051e4a6e188d..7c4c83f1d84ff9bac155d9eb7d19ef3f884588b0 100644
--- a/src/3d/variables_3d.F90
+++ b/src/3d/variables_3d.F90
@@ -112,6 +112,7 @@
 !
 ! !USES:
    use domain,     only: imin,imax,jmin,jmax,kmax
+   use field_manager
    IMPLICIT NONE
 !
 ! !PUBLIC DATA MEMBERS:
@@ -229,6 +230,172 @@
    end subroutine init_variables_3d
 !EOC
 
+#ifndef NO_3D
+!-----------------------------------------------------------------------
+!BOP
+!
+! !ROUTINE: register_3d_variables() - register GETM variables.
+!
+! !INTERFACE:
+   subroutine register_3d_variables(fm,runtype)
+!
+! !DESCRIPTION:
+!
+! !USES:
+!KB   use variables_3d
+   IMPLICIT NONE
+!
+! !INPUT PARAMETERS:
+   type (type_field_manager) :: fm
+   integer, intent(in)       :: runtype
+!
+! !REVISION HISTORY:
+!  Original author(s): Karsten Bolding & Jorn Bruggeman
+!
+! !LOCAL VARIABLES:
+   integer,parameter :: rk = kind(_ONE_)
+!EOP
+!-----------------------------------------------------------------------
+!BOC
+   LEVEL2 'register_3d_variables()'
+
+!:: kmin(I2DFIELD)
+!:: kumin(I2DFIELD)
+!:: kvmin(I2DFIELD)
+!:: kmin_pmz(I2DFIELD)
+!:: kumin_pmz(I2DFIELD)
+!:: kvmin_pmz(I2DFIELD)
+
+!:: uu(I3DFIELD)
+!:: vv(I3DFIELD)
+!:: ww(I3DFIELD)
+#ifdef _MOMENTUM_TERMS_
+!:: tdv_u(I3DFIELD)
+!:: adv_u(I3DFIELD)
+!:: vsd_u(I3DFIELD)
+!:: hsd_u(I3DFIELD)
+!:: cor_u(I3DFIELD)
+!:: epg_u(I3DFIELD)
+!:: ipg_u(I3DFIELD)
+
+!:: tdv_v(I3DFIELD)
+!:: adv_v(I3DFIELD)
+!:: vsd_v(I3DFIELD)
+!:: hsd_v(I3DFIELD)
+!:: cor_v(I3DFIELD)
+!:: epg_v(I3DFIELD)
+!:: ipg_v(I3DFIELD)
+#endif
+#ifdef STRUCTURE_FRICTION
+!:: sf(I3DFIELD)
+#endif
+!:: ho(I3DFIELD)
+!:: hn(I3DFIELD)
+!:: huo(I3DFIELD)
+!:: hun(I3DFIELD)
+!:: hvo(I3DFIELD)
+!:: hvn(I3DFIELD)
+!:: hcc(I3DFIELD)
+!:: uuEx(I3DFIELD)
+!:: vvEx(I3DFIELD)
+!:: num(I3DFIELD)
+!:: nuh(I3DFIELD)
+
+! 3D turbulent fields
+!:: tke(I3DFIELD)
+!:: eps(I3DFIELD)
+!:: SS(I3DFIELD)
+#ifndef NO_BAROCLINIC
+! 3D baroclinic fields
+!:: NN(I3DFIELD)
+!:: S(I3DFIELD)
+!:: T(I3DFIELD)
+!:: rho(I3DFIELD)
+!:: rad(I3DFIELD)
+!:: buoy(I3DFIELD)
+!:: alpha(I3DFIELD)
+!:: beta(I3DFIELD)
+!:: idpdx(I3DFIELD)
+!:: idpdy(I3DFIELD)
+!:: light(I3DFIELD)
+#endif
+
+#ifdef SPM
+! suspended matter
+!:: spm(I3DFIELD)
+!:: spm_ws(I3DFIELD)
+!:: spm_pool(I2DFIELD)
+#endif
+
+! 2D fields in 3D domain
+!:: sseo(I2DFIELD)
+!:: ssen(I2DFIELD)
+!:: Dn(I2DFIELD)
+!:: ssuo(I2DFIELD)
+!:: ssun(I2DFIELD)
+!:: ssvo(I2DFIELD)
+!:: ssvn(I2DFIELD)
+!:: Dun,Dvn
+
+! 3D friction in 3D domain
+!:: rru(I2DFIELD)
+!:: rrv(I2DFIELD)
+!:: taus(I2DFIELD)
+!:: taubx(I2DFIELD)
+!:: tauby(I2DFIELD)
+!:: taub(I2DFIELD)
+
+! light attenuation
+!:: A(I2DFIELD)
+!:: g1(I2DFIELD)
+!:: g2(I2DFIELD)
+
+!  category - 3d
+   if (runtype .ge. 2) then
+      call fm%register('zc', 'm', 'center coordinate', standard_name='', dimensions=(/id_dim_z/),data3d=zc(_3D_W_), category='grid', part_of_state=.false.)
+      call fm%register('hn', 'm', 'layer thickness', standard_name='cell_thickness', dimensions=(/id_dim_z/),data3d=hn(_3D_W_), category='grid', part_of_state=.true.)
+      call fm%register('hun', 'm', 'layer thickness - U-points', standard_name='cell_thickness', dimensions=(/id_dim_z/),data3d=hun(_3D_W_), category='grid', output_level=output_level_debug)
+      call fm%register('hvn', 'm', 'layer thickness - V-points', standard_name='cell_thickness', dimensions=(/id_dim_z/),data3d=hvn(_3D_W_), category='grid', output_level=output_level_debug)
+      call fm%register('ho', 'm', 'old layer thickness', standard_name='cell_thickness', dimensions=(/id_dim_z/),data3d=ho(_3D_W_), category='grid', output_level=output_level_debug)
+      call fm%register('ssen', 'm', 'elevation at T-points (3D)', standard_name='', data2d=ssen(_2D_W_), category='3d', fill_value=-9999.0_rk, output_level=output_level_debug, part_of_state=.true.)
+      call fm%register('ssun', 'm', 'elevation at U-points (3D)', standard_name='', data2d=ssun(_2D_W_), category='3d', output_level=output_level_debug, part_of_state=.true.)
+      call fm%register('ssvn', 'm', 'elevation at V-points (3D)', standard_name='', data2d=ssvn(_2D_W_), category='3d', output_level=output_level_debug, part_of_state=.true.)
+      call fm%register('sseo', 'm', 'old elevation at T-points (3D)', standard_name='', data2d=sseo(_2D_W_), category='3d', fill_value=-9999.0_rk, output_level=output_level_debug, part_of_state=.true.)
+      call fm%register('uu', 'm2/s', 'transport in local x-direction (3D)', standard_name='', dimensions=(/id_dim_z/), data3d=uu(_3D_W_), category='3d', output_level=output_level_debug, part_of_state=.true.)
+      call fm%register('vv', 'm2/s', 'transport in local y-direction (3D)', standard_name='', dimensions=(/id_dim_z/), data3d=vv(_3D_W_), category='3d', output_level=output_level_debug, part_of_state=.true.)
+      call fm%register('ww', 'm/s', 'grid-related vertical velocity', standard_name='', dimensions=(/id_dim_z/), data3d=ww(_3D_W_), category='3d', output_level=output_level_debug)
+      call fm%register('SS', 's-2', 'shear frequency squared', standard_name='', dimensions=(/id_dim_z/), data3d=SS(_3D_W_), category='3d', output_level=output_level_debug)
+
+   end if
+
+!  category - turbulence
+   if (runtype .ge. 2) then
+      call fm%register('tke' , 'm2/s2', 'TKE'        , standard_name='', dimensions=(/id_dim_z/), data3d=tke(_3D_W_), category='turbulence', output_level=output_level_debug)
+      call fm%register('diss', 'm2/s3', 'dissipation', standard_name='', dimensions=(/id_dim_z/), data3d=eps(_3D_W_), category='turbulence', output_level=output_level_debug)
+      call fm%register('num' , 'm2/s' , 'viscosity'  , standard_name='', dimensions=(/id_dim_z/), data3d=num(_3D_W_), category='turbulence', output_level=output_level_debug)
+      call fm%register('nuh' , 'm2/s' , 'diffusivity', standard_name='', dimensions=(/id_dim_z/), data3d=nuh(_3D_W_), category='turbulence', output_level=output_level_debug)
+   end if
+
+#ifndef NO_BAROCLINIC
+!  category - baroclinic
+   if (runtype .ge. 3) then
+      call fm%register('temp', 'Celsius', 'temperature', standard_name='', dimensions=(/id_dim_z/), fill_value=-9999.0_rk, data3d=T  (_3D_W_), category='baroclinic', part_of_state=.true.)
+      call fm%register('salt', '1e-3'   , 'salinity'   , standard_name='', dimensions=(/id_dim_z/), fill_value=-9999.0_rk, data3d=S  (_3D_W_), category='baroclinic', part_of_state=.true.)
+      call fm%register('rho' , 'kg/m3'  , 'density'    , standard_name='', dimensions=(/id_dim_z/), fill_value=-9999.0_rk, data3d=rho(_3D_W_), category='baroclinic', output_level=output_level_debug)
+      call fm%register('NN', 's-2', 'buoyancy frequency squared', standard_name='', dimensions=(/id_dim_z/), data3d=NN(_3D_W_), category='baroclinic', output_level=output_level_debug)
+      call fm%register('idpdx', 'm2/s2', 'baroclinic pressure gradient - x', standard_name='', dimensions=(/id_dim_z/),data3d=idpdx(_3D_W_), category='baroclinic', output_level=output_level_debug)
+#ifndef SLICE_MODEL
+      call fm%register('idpdy', 'm2/s2', 'baroclinic pressure gradient - y', standard_name='', dimensions=(/id_dim_z/),data3d=idpdy(_3D_W_), category='baroclinic', output_level=output_level_debug)
+#endif
+   end if
+#endif
+
+   return
+   end subroutine register_3d_variables
+!EOC
+#endif
+
+
 !-----------------------------------------------------------------------
 !BOP
 !
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 2db832fabcbe29f8410081a66844e78bda9f306a..bac83881d31e550ef14f18b76363824d2e5fe0a3 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -79,11 +79,6 @@ if(GETM_USE_FABM)
 
 endif(GETM_USE_FABM)
 
-option(GETM_USE_FLEXIBLE_OUTPUT "Use new output manager" ON)
-if(GETM_USE_FLEXIBLE_OUTPUT)
-  add_definitions(-D_FLEXIBLE_OUTPUT_)
-endif()
-
 option(GETM_USE_STRUCTURE_FRICTION "Include structure friction" OFF)
 if(GETM_USE_STRUCTURE_FRICTION)
    add_definitions(-DSTRUCTURE_FRICTION)
@@ -358,6 +353,7 @@ add_library(output OBJECT
             output/diagnostic_variables.F90
             output/ncdf_out.F90
             output/nesting.F90
+            output/output_processing.F90
             output/output.F90
            )
 
@@ -369,12 +365,12 @@ add_library(output OBJECT
 
 add_dependencies(domain futils)
 add_dependencies(meteo futils domain)
-add_dependencies(2d domain)
+add_dependencies(2d domain util)
 add_dependencies(3d meteo domain 2d)
 add_dependencies(output futils domain 2d 3d meteo)
 add_dependencies(ncdf futils domain output)
 add_dependencies(input futils 2d 3d ncdf)
-add_dependencies(getm input)
+add_dependencies(getm input output)
 
 set_property(TARGET futils ncdf input APPEND PROPERTY INCLUDE_DIRECTORIES "${NetCDF_INCLUDE_DIRS}")
 add_definitions(-DNETCDF_FMT -DREAL_4B=real\(4\))
@@ -385,13 +381,10 @@ find_path(GOTM_BASE src/gotm/gotm.F90 DOC "Path to GOTM source directory.")
 if(GOTM_BASE)
   set(GOTM_BUILD_LIBRARIES_ONLY ON)
   set(GOTM_USE_FABM ${GETM_USE_FABM} CACHE BOOL "Include support for Framework for Aquatic Biogeochemical Models (fabm.net)" FORCE)
-  set(GOTM_USE_FLEXIBLE_OUTPUT ${GETM_USE_FLEXIBLE_OUTPUT} CACHE BOOL "Use new output manager" FORCE)
   add_subdirectory(${GOTM_BASE}/src gotm)
   set(GOTM_INCLUDE_DIRS ${CMAKE_BINARY_DIR}/gotm/modules)
-  set(GOTM_LIBRARIES turbulence util)
-  if(GETM_USE_FLEXIBLE_OUTPUT)
-    set(GOTM_LIBRARIES ${GOTM_LIBRARIES} output_manager)
-  endif()
+  set(GOTM_LIBRARIES turbulence util output_manager)
+#KB  set(GOTM_LIBRARIES ${GOTM_LIBRARIES} output_manager)
   if(GETM_USE_FABM)
     set(GOTM_LIBRARIES ${GOTM_LIBRARIES} gotm_fabm)
     set_property(TARGET 3d getm APPEND PROPERTY INCLUDE_DIRECTORIES $<TARGET_PROPERTY:gotm_fabm,INTERFACE_INCLUDE_DIRECTORIES>)
diff --git a/src/Rules.make b/src/Rules.make
index c0cce165784f2523e88351f9dbbdc1791cd29a8c..fb35bde7f978957ec9983c8611ddd93682c3d6db 100644
--- a/src/Rules.make
+++ b/src/Rules.make
@@ -113,10 +113,7 @@ ifdef GOTM_PREFIX
 
 GOTMLIBDIR	= $(GOTM_PREFIX)/lib
 LINKDIRS	+= -L$(GOTMLIBDIR)
-ifeq ($(GETM_FLEXIBLE_OUTPUT),true)
-DEFINES += -D_FLEXIBLE_OUTPUT_
 EXTRA_LIBS	+= -loutput_manager
-endif
 ifeq ($(FABM),true)
 EXTRA_LIBS	+= -lgotm_fabm
 endif
diff --git a/src/cmake/Modules/FindGOTM.cmake b/src/cmake/Modules/FindGOTM.cmake
index 4f00f93b228dc6a8ce1f0050327c26d98fcd428b..7d71fb5f0948475c4603b85c6f500b424ed31eb6 100644
--- a/src/cmake/Modules/FindGOTM.cmake
+++ b/src/cmake/Modules/FindGOTM.cmake
@@ -6,13 +6,6 @@ find_path(GOTM_PREFIX
   DOC "Installation prefix for General Ocean Turbulence Models - gotm.net"
 )
 
-# Find GOTM output_manager library
-if(GETM_USE_FLEXIBLE_OUTPUT)
-  find_library(GOTM_OUTMAN NAMES output_manager
-               HINTS ${GOTM_PREFIX}/lib
-               DOC "GOTM output_manager library")
-endif(GETM_USE_FLEXIBLE_OUTPUT)
-
 # Find GOTM/FABM coupling library if USE_FABM
 if(GETM_USE_FABM)
   find_library(GOTM_FABM NAMES gotm_fabm
@@ -25,12 +18,17 @@ find_library(GOTM_TURBULENCE NAMES turbulence
              HINTS ${GOTM_PREFIX}/lib
              DOC "GOTM turbulence library")
 
+# Find GOTM output_manager library
+find_library(GOTM_OUTMAN NAMES output_manager
+             HINTS ${GOTM_PREFIX}/lib
+             DOC "GOTM output_manager library")
+
 # Find GOTM utility library
 find_library(GOTM_UTIL NAMES util
              HINTS ${GOTM_PREFIX}/lib
              DOC "GOTM utility library")
 
-set(GOTM_LIBRARIES ${GOTM_OUTMAN} ${GOTM_FABM} ${GOTM_TURBULENCE} ${GOTM_UTIL})
+set(GOTM_LIBRARIES ${GOTM_FABM} ${GOTM_TURBULENCE} ${GOTM_OUTMAN} ${GOTM_UTIL})
 
 # Store configurable path of GOTM include directory
 find_path(GOTM_INCLUDE_DIRS
@@ -39,7 +37,7 @@ find_path(GOTM_INCLUDE_DIRS
           DOC "GOTM include directories"
 )
 
-mark_as_advanced(GOTM_LIBRARIES GOTM_INCLUDE_DIRS GOTM_TURBULENCE GOTM_UTIL GOTM_FABM GOTM_OUTMAN)
+mark_as_advanced(GOTM_LIBRARIES GOTM_INCLUDE_DIRS GOTM_TURBULENCE GOTM_OUTMAN GOTM_UTIL GOTM_FABM)
 
 # Process default arguments (QUIET, REQUIRED)
 include(FindPackageHandleStandardArgs) 
diff --git a/src/domain/domain.F90 b/src/domain/domain.F90
index fe3a1e97a52f82fa095af61bee97c785783f0953..c3fd703ccdf321d9fcde856ecd606dcf67c85813 100644
--- a/src/domain/domain.F90
+++ b/src/domain/domain.F90
@@ -215,25 +215,25 @@
    call read_topo_file(bathy_format,bathymetry)
 
    if ( runtype .ge. 2 ) then
-   select case (vert_cord)
-      case(_SIGMA_COORDS_)
-         LEVEL2 'Using sigma coordinates'
-      case(_Z_COORDS_)
-         LEVEL2 'Using z-level coordinates'
-      case(_GENERAL_COORDS_)
-         LEVEL2 'Using general vertical coordinates'
-      case (_HYBRID_COORDS_) ! hybrid vertical coordinates
-         LEVEL2 'using hybrid vertical coordinates'
-         STDERR 'domain: hybrid_coordinates not coded yet'
-         stop
-      case (_ADAPTIVE_COORDS_) ! adaptive vertical coordinates
-         LEVEL2 'using adaptive vertical coordinates'
-      case default
-         call getm_error("init_domain()", &
-                         "A non valid vertical coordinate system has been chosen");
-   end select
-   allocate(ga(0:kmax),stat=rc)
-   if (rc /= 0) stop 'init_domain: Error allocating memory (ga)'
+      select case (vert_cord)
+         case(_SIGMA_COORDS_)
+            LEVEL2 'Using sigma coordinates'
+         case(_Z_COORDS_)
+            LEVEL2 'Using z-level coordinates'
+         case(_GENERAL_COORDS_)
+            LEVEL2 'Using general vertical coordinates'
+         case (_HYBRID_COORDS_) ! hybrid vertical coordinates
+            LEVEL2 'using hybrid vertical coordinates'
+            STDERR 'domain: hybrid_coordinates not coded yet'
+            stop
+         case (_ADAPTIVE_COORDS_) ! adaptive vertical coordinates
+            LEVEL2 'using adaptive vertical coordinates'
+         case default
+            call getm_error("init_domain()", &
+                            "A non valid vertical coordinate system has been chosen");
+      end select
+      allocate(ga(0:kmax),stat=rc)
+      if (rc /= 0) stop 'init_domain: Error allocating memory (ga)'
    end if
 
 !  Calculation masks
diff --git a/src/futils/getm_timers.F90 b/src/futils/getm_timers.F90
index 319c4219430e9bda8b996d698dc285dd8ca02ed9..476288d8b305b60d13accf5ff40f4fcf769472a0 100644
--- a/src/futils/getm_timers.F90
+++ b/src/futils/getm_timers.F90
@@ -73,6 +73,8 @@
    ! This is for do_input and do_output
    integer, parameter :: TIM_INPUT       = 90   ! input
    integer, parameter :: TIM_OUTPUT      = 92   ! output
+   integer, parameter :: TIM_FLEX_OUTPUT = 95   ! output_manager
+   integer, parameter :: TIM_OUTPUT_PROC = 96   ! output_processing
    ! These catch stuff that are *also* measured somewhere else:
    integer, parameter :: TIM_ADV         = 100  ! 2d advection
    integer, parameter :: TIM_ADVH        = 101  ! 2d advection halo parts
@@ -178,6 +180,8 @@
    timernames(TIM_METEO)       = 'do_meteo'
    timernames(TIM_INPUT)       = 'do_input'
    timernames(TIM_OUTPUT)      = 'do_output'
+   timernames(TIM_FLEX_OUTPUT) = 'output_manager_save'
+   timernames(TIM_OUTPUT_PROC) = 'do_output_processing'
 
 #ifdef GETM_PARALLEL
    timernames(TIM_MOMENTUMH)   = ' momentum-halo'
diff --git a/src/getm/cleanup.F90 b/src/getm/cleanup.F90
index 25c22e8f7107f0e66db0f826fc79c45f08b05c54..a505e9314645b83376739c747668c91ca8d3314c 100644
--- a/src/getm/cleanup.F90
+++ b/src/getm/cleanup.F90
@@ -22,9 +22,7 @@
 #endif
    use register_all_variables, only: fm
    use output, only: clean_output
-#if defined(_FLEXIBLE_OUTPUT_)
    use output_manager
-#endif
    use kurt_parallel, only: clean_parallel
    IMPLICIT NONE
 !
@@ -63,9 +61,7 @@
    if( .not. dryrun ) then
       call clean_output(runtype,loop)
    end if
-#if defined(_FLEXIBLE_OUTPUT_)
    call output_manager_clean()
-#endif
 
    call fm%finalize()
 
diff --git a/src/getm/compilation_options.F90 b/src/getm/compilation_options.F90
index bcfd0329cb080b2a336dffffb3b464e418609f4b..47d8a3e41dcf92d3d948d57dd57418f425194a55 100644
--- a/src/getm/compilation_options.F90
+++ b/src/getm/compilation_options.F90
@@ -159,9 +159,6 @@
 #ifdef _NCDF_SAVE_DOUBLE_
    LEVEL1 '_NCDF_SAVE_DOUBLE_'
 #endif
-#ifdef _FLEXIBLE_OUTPUT_
-   LEVEL1 '_FLEXIBLE_OUTPUT_'
-#endif
 
    STDERR LINE
 
diff --git a/src/getm/initialise.F90 b/src/getm/initialise.F90
index fb2a6ece65555c7409c786fafac74180a84010e7..4b4484e3760fc0c84dbb06c973ce79d6ba4fbd65 100644
--- a/src/getm/initialise.F90
+++ b/src/getm/initialise.F90
@@ -11,28 +11,25 @@
 !
 ! !USES:
    use register_all_variables
-#ifdef _FLEXIBLE_OUTPUT_
    use output_manager_core, only:output_manager_host=>host, type_output_manager_host=>type_host
    use time, only: CalDat,JulDay
    use output_manager
-#endif
    IMPLICIT NONE
 !
 ! !PUBLIC DATA MEMBERS:
    public                              :: init_model
    integer                             :: runtype=1
    logical                             :: dryrun=.false.
+   logical                             :: list_variables=.false.
 !
 ! !REVISION HISTORY:
 !  Original author(s): Karsten Bolding & Hans Burchard
 
-#ifdef _FLEXIBLE_OUTPUT_
    type,extends(type_output_manager_host) :: type_getm_host
    contains
       procedure :: julian_day => getm_host_julian_day
       procedure :: calendar_date => getm_host_calendar_date
    end type
-#endif
 !
 !EOP
 !-----------------------------------------------------------------------
@@ -53,6 +50,7 @@
    use halo_mpi, only: init_mpi,print_MPI_info
 #endif
    use output, only: init_output,do_output,restart_file,out_dir
+   use output_processing
    use input,  only: init_input
    use domain, only: init_domain
    use domain, only: H
@@ -269,10 +267,10 @@
 #endif
    end if
 #endif
+   call init_output_processing()
 
    call init_register_all_variables(runtype)
 
-#ifdef _FLEXIBLE_OUTPUT_
    allocate(type_getm_host::output_manager_host)
    if (myid .ge. 0) then
       write(postfix,'(A,I4.4)') '.',myid
@@ -280,9 +278,9 @@
    else
       call output_manager_init(fm,title)
    end if
-#endif
 
    call do_register_all_variables(runtype)
+   if (list_variables) call fm%list()
 
 !   call init_output(runid,title,start,runtype,dryrun,myid)
    call init_output(runid,title,start,runtype,dryrun,myid,MinN,MaxN,save_initial)
@@ -353,10 +351,12 @@
    end if
 
    if (.not. dryrun) then
+      if (save_initial) then
+      call output_manager_prepare_save(julianday, int(secondsofday), 0, int(MinN-1))
+      call do_output_processing()
+      call output_manager_save(julianday,secondsofday,MinN-1)
+      end if
       call do_output(runtype,MinN-1,timestep)
-#ifdef _FLEXIBLE_OUTPUT_
-      if (save_initial) call output_manager_save(julianday,secondsofday,MinN-1)
-#endif
    end if
 
 #ifdef DEBUG
@@ -369,7 +369,6 @@
 
 !-----------------------------------------------------------------------
 
-#ifdef _FLEXIBLE_OUTPUT_
    subroutine getm_host_julian_day(self,yyyy,mm,dd,julian)
       class (type_getm_host), intent(in) :: self
       integer, intent(in)  :: yyyy,mm,dd
@@ -383,7 +382,6 @@
       integer, intent(out) :: yyyy,mm,dd
       call CalDat(julian,yyyy,mm,dd)
    end subroutine
-#endif
 
    end module initialise
 
diff --git a/src/getm/integration.F90 b/src/getm/integration.F90
index 101a99a7321f9b3d9b49466f200e40795fbf8df8..fac1c83de262029c7fc13adaa095eef65cf48b8c 100644
--- a/src/getm/integration.F90
+++ b/src/getm/integration.F90
@@ -62,13 +62,13 @@
    use suspended_matter, only: spm_calc,do_spm
 #endif
    use input,    only: do_input
+   use output_processing, only: do_output_processing
    use output,   only: do_output
 #ifdef TEST_NESTING
    use nesting,   only: nesting_file
 #endif
-#ifdef _FLEXIBLE_OUTPUT_
    use output_manager
-#endif
+   use getm_timers, only: tic, toc, TIM_FLEX_OUTPUT, TIM_OUTPUT_PROC
    IMPLICIT NONE
 !
 ! !INPUT PARAMETERS:
@@ -153,12 +153,19 @@
          call nesting_file(WRITING)
       end if
 #endif
+
       call update_time(n)
 
-      call do_output(runtype,n,timestep)
-#ifdef _FLEXIBLE_OUTPUT_
+      call tic(TIM_FLEX_OUTPUT)
+      call output_manager_prepare_save(julianday, int(secondsofday), 0, int(n))
+      call toc(TIM_FLEX_OUTPUT)
+      call tic(TIM_OUTPUT_PROC)
+      call do_output_processing()
+      call toc(TIM_OUTPUT_PROC)
+      call tic(TIM_FLEX_OUTPUT)
       call output_manager_save(julianday,secondsofday,n)
-#endif
+      call toc(TIM_FLEX_OUTPUT)
+      call do_output(runtype,n,timestep)
 #ifdef DIAGNOSE
       call diagnose(n,MaxN,runtype)
 #endif
diff --git a/src/getm/main.F90 b/src/getm/main.F90
index c469ca61c663567f49e50ec9a73a408241c53099..48869cb398ef06deb812cfd9d4269c44a2670af6 100644
--- a/src/getm/main.F90
+++ b/src/getm/main.F90
@@ -116,7 +116,7 @@
 
 !-----------------------------------------------------------------------
    subroutine cmdline
-   use initialise, only: dryrun
+   use initialise, only: dryrun,list_variables
    IMPLICIT NONE
    character(len=64)    :: arg
    integer              :: i
@@ -137,6 +137,8 @@
          call compilation_options()
          LEVEL0
          stop
+      case ('-l', '--list_variables')
+         list_variables=.true.
       case ('-h', '--help')
          call print_help()
          stop
@@ -165,9 +167,10 @@
      print '(a)', ''
      print '(a)', 'cmdline options:'
      print '(a)', ''
-     print '(a)', '  -v, --version     print version information and exit'
-     print '(a)', '  -c, --compile     print compilation options'
-     print '(a)', '  -h, --help        print usage information and exit'
+     print '(a)', '  -v, --version        print version information and exit'
+     print '(a)', '  -c, --compile        print compilation options'
+     print '(a)', '  -l, --list_variables list possible output variables'
+     print '(a)', '  -h, --help           print usage information and exit'
      print '(a)', ''
      print '(a)', 'visit getm.eu for further info'
      print '(a)', 'consider subscribing to getm-users@googlegroups.com'
diff --git a/src/getm/register_all_variables.F90 b/src/getm/register_all_variables.F90
index 90368732c9f32cdbe0b0f4007dc9f9d537b40935..a4b65501292c7689f26849b1fff32fc5c2ba7fb0 100644
--- a/src/getm/register_all_variables.F90
+++ b/src/getm/register_all_variables.F90
@@ -11,6 +11,12 @@
 !
 ! !USES:
    use field_manager
+   use variables_2d, only: register_2d_variables
+   use variables_3d, only: register_3d_variables
+#ifdef _FABM_
+   use getm_fabm, only: register_fabm_variables
+#endif
+   use output_processing, only: register_processed_variables
    IMPLICIT NONE
 !
 !  default: all is private.
@@ -19,6 +25,7 @@
 ! !PUBLIC MEMBER FUNCTIONS:
    public :: init_register_all_variables
    public :: do_register_all_variables
+   public :: finalize_register_all_variables
 !
 ! !PUBLIC DATA MEMBERS:
    type (type_field_manager), public, target :: fm
@@ -37,6 +44,12 @@
    character(len=16)         :: zname='sigma'
    character(len=64)         :: zlongname='sigma'
    character(len=16)         :: zunits='sigma'
+   character(len=16),parameter :: lonname     = 'lonc'
+   character(len=16),parameter :: lonlongname = 'longitude'
+   character(len=16),parameter :: lonunits    = 'degrees_east'
+   character(len=16),parameter :: latname     = 'latc'
+   character(len=16),parameter :: latlongname = 'latitude'
+   character(len=16),parameter :: latunits    = 'degrees_north'
 !
 !-----------------------------------------------------------------------
 
@@ -77,13 +90,13 @@
          ylongname = 'y'
          yunits    = 'm'
       case (2)
-         xname     = 'lonc'
-         xlongname = 'longitude'
-         xunits    = 'degrees_east'
-         yname     = 'latc'
-         ylongname = 'latitude'
-         yunits    = 'degrees_north'
-      case (3)
+         xname     = lonname
+         xlongname = lonlongname
+         xunits    = lonunits
+         yname     = latname
+         ylongname = latlongname
+         yunits    = latunits
+      case (3,4)
          xname     = 'xic'
          xlongname = 'xic'
          yname     = 'etac'
@@ -153,16 +166,17 @@
    LEVEL1 'do_register_all_variables()'
    call register_domain_variables(runtype)
    call register_meteo_variables()
-   call register_2d_variables()
+   call register_2d_variables(fm)
 #ifndef NO_3D
-   call register_3d_variables(runtype)
+   call register_3d_variables(fm,runtype)
 #endif
 #ifdef _FABM_
-   call register_fabm_variables()
+   call register_fabm_variables(fm)
 #endif
 #if 0
    call register_diagnostic_variables()
 #endif
+   call register_processed_variables(fm)
 
    return
    end subroutine do_register_all_variables
@@ -203,6 +217,11 @@
    end if
 #endif
 
+   if (have_lonlat .and. grid_type.ne.2) then
+      call fm%register(trim(lonname),trim(lonunits),trim(lonlongname),dimensions=(/id_dim_lon,id_dim_lat/),no_default_dimensions=.true.,data2d=lonc(_2D_W_), category='domain',output_level=output_level_required)
+      call fm%register(trim(latname),trim(latunits),trim(latlongname),dimensions=(/id_dim_lon,id_dim_lat/),no_default_dimensions=.true.,data2d=latc(_2D_W_), category='domain',output_level=output_level_required)
+   end if
+
    call fm%register('bathymetry', 'm', 'bathymetry', standard_name='bathymetry', dimensions=(/id_dim_lon,id_dim_lat/), no_default_dimensions=.true., fill_value=-10._rk, data2d=H(_2D_W_), category='domain',output_level=output_level_required)
 
 !  register -  metric
@@ -279,239 +298,6 @@
    end subroutine register_meteo_variables
 !EOC
 
-!-----------------------------------------------------------------------
-!BOP
-!
-! !ROUTINE: register_2d_variables() - register GETM variables.
-!
-! !INTERFACE:
-   subroutine register_2d_variables()
-!
-! !DESCRIPTION:
-!
-! !USES:
-   use variables_2d
-   IMPLICIT NONE
-!
-! !REVISION HISTORY:
-!  Original author(s): Karsten Bolding & Jorn Bruggeman
-!
-! !LOCAL VARIABLES:
-!EOP
-!-----------------------------------------------------------------------
-!BOC
-   LEVEL2 'register_2d_variables()'
-
-!D(E2DFIELD)
-!DU,DV
-!z(E2DFIELD)
-!zo(E2DFIELD)
-!U(E2DFIELD)
-!V(E2DFIELD)
-!UEx(E2DFIELD)
-!VEx(E2DFIELD)
-!fU(E2DFIELD)
-!fV(E2DFIELD)
-!ru(E2DFIELD)
-!rv(E2DFIELD)
-!Uint(E2DFIELD)
-!Vint(E2DFIELD)
-!Uinto(E2DFIELD)
-!Vinto(E2DFIELD)
-!res_du(E2DFIELD)
-!res_u(E2DFIELD)
-!res_dv(E2DFIELD)
-!res_v(E2DFIELD)
-!kbk
-!ruu(E2DFIELD)
-!rvv(E2DFIELD)
-!kbk
-!SlUx(E2DFIELD)
-!SlVx(E2DFIELD)
-!Slru(E2DFIELD)
-!Slrv(E2DFIELD)
-!zub(E2DFIELD)
-!zvb(E2DFIELD)
-!zub0(E2DFIELD)
-!zvb0(E2DFIELD)
-!An(E2DFIELD)
-!AnX(E2DFIELD)
-!fwf(E2DFIELD)
-!fwf_int(E2DFIELD)
-!EWbdy(jmax),ENbdy(imax),EEbdy(jmax),ESbdy(imax)
-
-
-!  category - 2d
-   call fm%register('z', 'm', 'sea surface elevation', standard_name='sea surface elevation', fill_value=-9999.0_rk, data2d=z(_2D_W_), category="2d")
-
-   call fm%register('zo', 'm', 'sea surface elevation', standard_name='sea surface elevation', fill_value=-9999.0_rk, data2d=zo(_2D_W_), category="2d", output_level=output_level_debug)
-   call fm%register('D', 'm', 'water depth', standard_name='water depth', fill_value=-10009.0_rk, data2d=D(_2D_W_), category="2d")
-   call fm%register('U', 'm2/s', 'transport in local x-direction', standard_name='', data2d=U(_2D_W_), category='2d', output_level=output_level_debug)
-   call fm%register('V', 'm2/s', 'transport in local y-direction', standard_name='', data2d=V(_2D_W_), category='2d', output_level=output_level_debug)
-
-   return
-   end subroutine register_2d_variables
-!EOC
-
-#ifndef NO_3D
-!-----------------------------------------------------------------------
-!BOP
-!
-! !ROUTINE: register_3d_variables() - register GETM variables.
-!
-! !INTERFACE:
-   subroutine register_3d_variables(runtype)
-!
-! !DESCRIPTION:
-!
-! !USES:
-   use variables_3d
-   IMPLICIT NONE
-!
-! !INPUT PARAMETERS:
-   integer, intent(in)               :: runtype
-!
-! !REVISION HISTORY:
-!  Original author(s): Karsten Bolding & Jorn Bruggeman
-!
-! !LOCAL VARIABLES:
-!EOP
-!-----------------------------------------------------------------------
-!BOC
-   LEVEL2 'register_3d_variables()'
-
-!:: kmin(I2DFIELD)
-!:: kumin(I2DFIELD)
-!:: kvmin(I2DFIELD)
-!:: kmin_pmz(I2DFIELD)
-!:: kumin_pmz(I2DFIELD)
-!:: kvmin_pmz(I2DFIELD)
-
-!:: uu(I3DFIELD)
-!:: vv(I3DFIELD)
-!:: ww(I3DFIELD)
-#ifdef _MOMENTUM_TERMS_
-!:: tdv_u(I3DFIELD)
-!:: adv_u(I3DFIELD)
-!:: vsd_u(I3DFIELD)
-!:: hsd_u(I3DFIELD)
-!:: cor_u(I3DFIELD)
-!:: epg_u(I3DFIELD)
-!:: ipg_u(I3DFIELD)
-
-!:: tdv_v(I3DFIELD)
-!:: adv_v(I3DFIELD)
-!:: vsd_v(I3DFIELD)
-!:: hsd_v(I3DFIELD)
-!:: cor_v(I3DFIELD)
-!:: epg_v(I3DFIELD)
-!:: ipg_v(I3DFIELD)
-#endif
-#ifdef STRUCTURE_FRICTION
-!:: sf(I3DFIELD)
-#endif
-!:: ho(I3DFIELD)
-!:: hn(I3DFIELD)
-!:: huo(I3DFIELD)
-!:: hun(I3DFIELD)
-!:: hvo(I3DFIELD)
-!:: hvn(I3DFIELD)
-!:: hcc(I3DFIELD)
-!:: uuEx(I3DFIELD)
-!:: vvEx(I3DFIELD)
-!:: num(I3DFIELD)
-!:: nuh(I3DFIELD)
-
-! 3D turbulent fields
-!:: tke(I3DFIELD)
-!:: eps(I3DFIELD)
-!:: SS(I3DFIELD)
-#ifndef NO_BAROCLINIC
-! 3D baroclinic fields
-!:: NN(I3DFIELD)
-!:: S(I3DFIELD)
-!:: T(I3DFIELD)
-!:: rho(I3DFIELD)
-!:: rad(I3DFIELD)
-!:: buoy(I3DFIELD)
-!:: alpha(I3DFIELD)
-!:: beta(I3DFIELD)
-!:: idpdx(I3DFIELD)
-!:: idpdy(I3DFIELD)
-!:: light(I3DFIELD)
-#endif
-
-#ifdef SPM
-! suspended matter
-!:: spm(I3DFIELD)
-!:: spm_ws(I3DFIELD)
-!:: spm_pool(I2DFIELD)
-#endif
-
-! 2D fields in 3D domain
-!:: sseo(I2DFIELD)
-!:: ssen(I2DFIELD)
-!:: Dn(I2DFIELD)
-!:: ssuo(I2DFIELD)
-!:: ssun(I2DFIELD)
-!:: ssvo(I2DFIELD)
-!:: ssvn(I2DFIELD)
-!:: Dun,Dvn
-
-! 3D friction in 3D domain
-!:: rru(I2DFIELD)
-!:: rrv(I2DFIELD)
-!:: taus(I2DFIELD)
-!:: taubx(I2DFIELD)
-!:: tauby(I2DFIELD)
-!:: taub(I2DFIELD)
-
-! light attenuation
-!:: A(I2DFIELD)
-!:: g1(I2DFIELD)
-!:: g2(I2DFIELD)
-
-!  category - 3d
-   if (runtype .ge. 2) then
-      call fm%register('hn', 'm', 'layer thickness', standard_name='cell_thickness', dimensions=(/id_dim_z/),data3d=hn(_3D_W_), category='grid')
-      call fm%register('hun', 'm', 'layer thickness - U-points', standard_name='cell_thickness', dimensions=(/id_dim_z/),data3d=hun(_3D_W_), category='grid', output_level=output_level_debug)
-      call fm%register('hvn', 'm', 'layer thickness - V-points', standard_name='cell_thickness', dimensions=(/id_dim_z/),data3d=hvn(_3D_W_), category='grid', output_level=output_level_debug)
-      call fm%register('ho', 'm', 'old layer thickness', standard_name='cell_thickness', dimensions=(/id_dim_z/),data3d=ho(_3D_W_), category='grid', output_level=output_level_debug)
-      call fm%register('ssen', 'm', 'elevation at T-points (3D)', standard_name='', data2d=ssen(_2D_W_), category='3d', fill_value=-9999.0_rk, output_level=output_level_debug)
-      call fm%register('ssun', 'm', 'elevation at U-points (3D)', standard_name='', data2d=ssun(_2D_W_), category='3d', output_level=output_level_debug)
-      call fm%register('ssvn', 'm', 'elevation at V-points (3D)', standard_name='', data2d=ssvn(_2D_W_), category='3d', output_level=output_level_debug)
-      call fm%register('sseo', 'm', 'old elevation at T-points (3D)', standard_name='', data2d=sseo(_2D_W_), category='3d', fill_value=-9999.0_rk, output_level=output_level_debug)
-      call fm%register('uu', 'm2/s', 'transport in local x-direction (3D)', standard_name='', dimensions=(/id_dim_z/), data3d=uu(_3D_W_), category='3d', output_level=output_level_debug)
-      call fm%register('vv', 'm2/s', 'transport in local y-direction (3D)', standard_name='', dimensions=(/id_dim_z/), data3d=vv(_3D_W_), category='3d', output_level=output_level_debug)
-      call fm%register('ww', 'm/s', 'grid-related vertical velocity', standard_name='', dimensions=(/id_dim_z/), data3d=ww(_3D_W_), category='3d', output_level=output_level_debug)
-      call fm%register('SS', 's-2', 'shear frequency squared', standard_name='', dimensions=(/id_dim_z/), data3d=SS(_3D_W_), category='3d', output_level=output_level_debug)
-   end if
-
-!  category - turbulence
-   if (runtype .ge. 2) then
-      call fm%register('tke' , 'm2/s2', 'TKE'        , standard_name='', dimensions=(/id_dim_z/), data3d=tke(_3D_W_), category='turbulence', output_level=output_level_debug)
-      call fm%register('diss', 'm2/s3', 'dissipation', standard_name='', dimensions=(/id_dim_z/), data3d=eps(_3D_W_), category='turbulence', output_level=output_level_debug)
-      call fm%register('num' , 'm2/s' , 'viscosity'  , standard_name='', dimensions=(/id_dim_z/), data3d=num(_3D_W_), category='turbulence', output_level=output_level_debug)
-      call fm%register('nuh' , 'm2/s' , 'diffusivity', standard_name='', dimensions=(/id_dim_z/), data3d=nuh(_3D_W_), category='turbulence', output_level=output_level_debug)
-   end if
-
-#ifndef NO_BAROCLINIC
-!  category - baroclinic
-   if (runtype .ge. 3) then
-      call fm%register('temp', 'Celsius', 'temperature', standard_name='', dimensions=(/id_dim_z/), fill_value=-9999.0_rk, data3d=T  (_3D_W_), category='baroclinic')
-      call fm%register('salt', '1e-3'   , 'salinity'   , standard_name='', dimensions=(/id_dim_z/), fill_value=-9999.0_rk, data3d=S  (_3D_W_), category='baroclinic')
-      call fm%register('rho' , 'kg/m3'  , 'density'    , standard_name='', dimensions=(/id_dim_z/), fill_value=-9999.0_rk, data3d=rho(_3D_W_), category='baroclinic', output_level=output_level_debug)
-      call fm%register('NN', 's-2', 'buoyancy frequency squared', standard_name='', dimensions=(/id_dim_z/), data3d=NN(_3D_W_), category='baroclinic', output_level=output_level_debug)
-      call fm%register('idpdx', 'm2/s2', 'baroclinic pressure gradient - x', standard_name='', dimensions=(/id_dim_z/),data3d=idpdx(_3D_W_), category='baroclinic', output_level=output_level_debug)
-      call fm%register('idpdy', 'm2/s2', 'baroclinic pressure gradient - y', standard_name='', dimensions=(/id_dim_z/),data3d=idpdy(_3D_W_), category='baroclinic', output_level=output_level_debug)
-   end if
-#endif
-
-   return
-   end subroutine register_3d_variables
-!EOC
-#endif
-
 !-----------------------------------------------------------------------
 !BOP
 !
@@ -541,71 +327,41 @@
    end subroutine register_diagnostic_variables
 !EOC
 
-#ifdef _FABM_
 !-----------------------------------------------------------------------
 !BOP
 !
-! !ROUTINE: register_fabm_variables() - register FABM variables.
+! !ROUTINE: finalize_register_all_variables() - send optional variables.
 !
 ! !INTERFACE:
-   subroutine register_fabm_variables()
+   subroutine finalize_register_all_variables(runtype)
 !
 ! !DESCRIPTION:
 !
 ! !USES:
+   use variables_2d
+#ifndef NO_3D
+   use variables_3d
+#endif
+#ifdef _FABM_
    use getm_fabm
+#endif
    IMPLICIT NONE
 !
 ! !INPUT PARAMETERS:
+   integer, intent(in)               :: runtype
 !
 ! !REVISION HISTORY:
 !  Original author(s): Karsten Bolding & Jorn Bruggeman
 !
 ! !LOCAL VARIABLES:
-  integer :: i,output_level
-  logical :: in_output
 !EOP
 !-----------------------------------------------------------------------
 !BOC
-   if (.not. fabm_calc) return
-   LEVEL2 'register_fabm_variables()'
-
-   do i=1,size(model%state_variables)
-      output_level = output_level_default
-      if (model%state_variables(i)%output==output_none) output_level = output_level_debug
-      call fm%register(model%state_variables(i)%name, model%state_variables(i)%units, &
-         model%state_variables(i)%long_name, minimum=model%state_variables(i)%minimum, maximum=model%state_variables(i)%maximum, &
-         fill_value=model%state_variables(i)%missing_value, dimensions=(/id_dim_z/), data3d=fabm_pel(_3D_W_,i), category='fabm'//model%state_variables(i)%target%owner%get_path(), output_level=output_level)
-   end do
-   do i=1,size(model%bottom_state_variables)
-      output_level = output_level_default
-      if (model%bottom_state_variables(i)%output==output_none) output_level = output_level_debug
-      call fm%register(model%bottom_state_variables(i)%name, model%bottom_state_variables(i)%units, &
-         model%bottom_state_variables(i)%long_name, minimum=model%bottom_state_variables(i)%minimum, &
-         maximum=model%bottom_state_variables(i)%maximum, fill_value=model%bottom_state_variables(i)%missing_value, &
-         data2d=fabm_ben(_2D_W_,i), category='fabm'//model%bottom_state_variables(i)%target%owner%get_path(), output_level=output_level)
-   end do
-   do i=1,size(model%diagnostic_variables)
-      output_level = output_level_default
-      if (model%diagnostic_variables(i)%output==output_none) output_level = output_level_debug
-      call fm%register(model%diagnostic_variables(i)%name, model%diagnostic_variables(i)%units, &
-         model%diagnostic_variables(i)%long_name, minimum=model%diagnostic_variables(i)%minimum, maximum=model%diagnostic_variables(i)%maximum, &
-         fill_value=model%diagnostic_variables(i)%missing_value, dimensions=(/id_dim_z/), data3d=fabm_diag(_3D_W_,i), category='fabm'//model%diagnostic_variables(i)%target%owner%get_path(), output_level=output_level, used=in_output)
-      if (in_output) model%diagnostic_variables(i)%save = .true.
-   end do
-   do i=1,size(model%horizontal_diagnostic_variables)
-      output_level = output_level_default
-      if (model%horizontal_diagnostic_variables(i)%output==output_none) output_level = output_level_debug
-      call fm%register(model%horizontal_diagnostic_variables(i)%name, model%horizontal_diagnostic_variables(i)%units, &
-         model%horizontal_diagnostic_variables(i)%long_name, minimum=model%horizontal_diagnostic_variables(i)%minimum, maximum=model%horizontal_diagnostic_variables(i)%maximum, &
-         fill_value=model%horizontal_diagnostic_variables(i)%missing_value, data2d=fabm_diag_hz(_2D_W_,i), category='fabm'//model%horizontal_diagnostic_variables(i)%target%owner%get_path(), output_level=output_level, used=in_output)
-      if (in_output) model%horizontal_diagnostic_variables(i)%save = .true.
-   end do
+   LEVEL1 'finalize_register_all_variables()'
 
    return
-   end subroutine register_fabm_variables
+   end subroutine finalize_register_all_variables
 !EOC
-#endif
 
 !-----------------------------------------------------------------------
 
diff --git a/src/output/Makefile b/src/output/Makefile
index 2a88bd303414e0e43a9908387a640ea77e971fd0..ded7cc5927e4b705e08f39e5f71c49bbe2501c7e 100644
--- a/src/output/Makefile
+++ b/src/output/Makefile
@@ -6,7 +6,7 @@ include ../Rules.make
 
 LIB	= $(LIBDIR)/liboutput${buildtype}.a
 
-MODSRC	= output.F90 parser.F90 variable_info.F90 \
+MODSRC	= output.F90 output_processing.F90 \
           ncdf_out.F90 ascii_out.F90 \
           diagnostic_variables.F90 nesting.F90
 
@@ -23,11 +23,9 @@ endif
 MOD += \
 ${LIB}(ncdf_out.o)		\
 ${LIB}(ascii_out.o)		\
+${LIB}(output_processing.o)		\
 ${LIB}(output.o)
 
-#$(LIB)(variable_info.o)		\
-#$(LIB)(parser.o)		\
-
 OBJ     =
 ifneq ($(GETM_NO_3D),true)
 OBJ += \
@@ -42,24 +40,12 @@ modules: $(MOD)
 
 objects: $(OBJ)
 
-tests: test_parser test_varinfo
-
-test_parser: modules objects test_parser.o
-	$(FC) -o $@ $@.o $(LDFLAGS) $(LIB)
-	$(RM) $@.o
-	./$@
-
-test_varinfo: modules objects test_varinfo.o
-	$(FC) -o $@ $@.o $(LDFLAGS) $(LIB)
-	$(RM) $@.o
-	./$@
-
 doc: $(SRC)
 	$(PROTEX) $(SRC) > $(DOCDIR)/output.tex
 	touch doc
 
 clean:
-	$(RM) $(LIB) $(MODDIR)/{output,ncdf_out,ascii_out}.{m,mod}
+	$(RM) $(LIB) $(MODDIR)/{output,output_processing,ncdf_out,ascii_out}.{m,mod}
 
 realclean: clean
 	$(RM) *.o doc
diff --git a/src/output/output_processing.F90 b/src/output/output_processing.F90
new file mode 100644
index 0000000000000000000000000000000000000000..76fdb744080c317aeff70089f5d84f8bf6d53cad
--- /dev/null
+++ b/src/output/output_processing.F90
@@ -0,0 +1,197 @@
+#include "cppdefs.h"
+!-----------------------------------------------------------------------
+!BOP
+!
+! !MODULE: output_processing
+!
+! !INTERFACE:
+   module output_processing
+!
+! !DESCRIPTION:
+!  This modules serves as a container for processing output variables. 
+
+! !USES:
+   use domain, only: imin,imax,jmin,jmax,kmax
+   use field_manager
+   IMPLICIT NONE
+!
+   private
+!
+! !PUBLIC DATA FUNCTIONS:
+   public init_output_processing, register_processed_variables, do_output_processing
+!
+! !PUBLIC DATA MEMBERS:
+   REALTYPE, dimension(:,:),   allocatable, target :: u2d, v2d
+   REALTYPE, dimension(:,:),   allocatable, target :: u2d_destag, v2d_destag
+   REALTYPE, dimension(:,:,:), allocatable, target :: u3d, v3d
+   REALTYPE, dimension(:,:,:), allocatable, target :: u3d_destag, v3d_destag
+!
+! !PRIVATE DATA MEMBERS:
+   logical, target:: u2d_use, v2d_use
+   logical, target:: u2d_destag_use, v2d_destag_use
+   logical, target:: u3d_use, v3d_use
+   logical, target:: u3d_destag_use, v3d_destag_use
+   integer, parameter :: rk = kind(_ONE_)
+!
+! !REVISION HISTORY:
+!  Original author(s): Karsten Bolding & Hans Burchard
+!
+!EOP
+!-----------------------------------------------------------------------
+
+   contains
+
+!-----------------------------------------------------------------------
+!BOP
+!
+! !IROUTINE: init_output_processing - read required variables
+!
+! !INTERFACE:
+   subroutine init_output_processing
+!
+! !USES:
+   IMPLICIT NONE
+!
+! !DESCRIPTION:
+!
+! !INPUT PARAMETERS:
+!    character(len=*), intent(in)        :: filename
+!
+! !REVISION HISTORY:
+!  Original author(s): Karsten Bolding
+!
+! !LOCAL VARIABLES:
+   integer                             :: rc
+!EOP
+!-------------------------------------------------------------------------
+!BOC
+   allocate(u2d(E2DFIELD),stat=rc)
+   if (rc /= 0) stop 'init_output_processing: Error allocating memory (u2d)'
+   u2d = 0._rk
+   allocate(v2d(E2DFIELD),stat=rc)
+   if (rc /= 0) stop 'init_output_processing: Error allocating memory (v2d)'
+   v2d = 0._rk
+   allocate(u2d_destag(E2DFIELD),stat=rc)
+   if (rc /= 0) stop 'init_output_processing: Error allocating memory (u2d_destag)'
+   u2d_destag = 0._rk
+   allocate(v2d_destag(E2DFIELD),stat=rc)
+   if (rc /= 0) stop 'init_output_processing: Error allocating memory (v2d_destag)'
+   v2d_destag = 0._rk
+
+#if 0
+   allocate(u3d(I3DFIELD),stat=rc)
+   if (rc /= 0) stop 'init_output_processing: Error allocating memory (u3d)'
+   allocate(v3d(I3DFIELD),stat=rc)
+   if (rc /= 0) stop 'init_output_processing: Error allocating memory (v3d)'
+   allocate(u3d_destag(I3DFIELD),stat=rc)
+   if (rc /= 0) stop 'init_output_processing: Error allocating memory (u3d_destag)'
+   allocate(v3d_destag(I3DFIELD),stat=rc)
+   if (rc /= 0) stop 'init_output_processing: Error allocating memory (v3d_destag)'
+#endif
+   return
+   end subroutine init_output_processing
+!EOC
+
+!-----------------------------------------------------------------------
+!BOP
+!
+! !ROUTINE: register_processed_variables()
+!
+! !INTERFACE:
+   subroutine register_processed_variables(fm)
+!
+! !DESCRIPTION:
+!
+! !USES:
+   IMPLICIT NONE
+!
+! !INPUT PARAMETERS:
+   type (type_field_manager) :: fm
+!
+! !REVISION HISTORY:
+!  Original author(s): Karsten Bolding & Jorn Bruggeman
+!
+! !LOCAL VARIABLES:
+!EOP
+!-----------------------------------------------------------------------
+!BOC
+   LEVEL2 'register_processed_variables()'
+
+   call fm%register('u2d', 'm/s', 'velocity in local x-direction', standard_name='', data2d=u2d(_2D_W_), fill_value=-9999._rk, category='velocities', used_now=u2d_use)
+   call fm%register('v2d', 'm/s', 'velocity in local y-direction', standard_name='', data2d=v2d(_2D_W_), fill_value=-9999._rk, category='velocities', used_now=v2d_use)
+   call fm%register('u2d-destag', 'm/s', 'velocity in local x-direction(destag)', standard_name='', data2d=u2d_destag(_2D_W_), fill_value=-9999._rk, category='velocities',output_level=output_level_debug, used_now=u2d_destag_use)
+   call fm%register('v2d-destag', 'm/s', 'velocity in local y-direction(destag)', standard_name='', data2d=v2d_destag(_2D_W_), fill_value=-9999._rk, category='velocities',output_level=output_level_debug, used_now=v2d_destag_use)
+
+   return
+   end subroutine register_processed_variables
+!EOC
+
+!-----------------------------------------------------------------------
+!BOP
+! !IROUTINE: do_output_processing - read required variables
+!
+! !INTERFACE:
+  subroutine do_output_processing
+!
+! !USES:
+   use domain, only: az, au, av
+   use variables_2d, only: z,D
+   use variables_2d, only: U,V,DU,DV
+   use variables_3d, only: kmin,hn,uu,hun,vv,hvn
+   IMPLICIT NONE
+!
+! !DESCRIPTION:
+!
+! !INPUT PARAMETERS:
+!
+! !REVISION HISTORY:
+!  Original author(s): Karsten Bolding
+!
+! !LOCAL VARIABLES:
+   REALTYPE, parameter                 :: vel_missing        =-9999.
+!EOP
+!-------------------------------------------------------------------------
+!BOC
+
+!  2D - velocities
+   if (u2d_use .and. v2d_use) then
+      call to_2d_vel(imin,jmin,imax,jmax,au,U,DU,vel_missing,       &
+                     imin,jmin,imax,jmax,u2d)
+      call to_2d_vel(imin,jmin,imax,jmax,av,V,DV,vel_missing,       &
+                     imin,jmin,imax,jmax,v2d)
+   end if
+   if (u2d_destag_use .and. v2d_destag_use) then
+      call to_2d_u(imin,jmin,imax,jmax,az,U,DU,vel_missing,      &
+                   imin,jmin,imax,jmax,u2d_destag)
+      call to_2d_v(imin,jmin,imax,jmax,az,V,DV,vel_missing,      &
+                   imin,jmin,imax,jmax,v2d_destag)
+   end if
+
+#if 0
+!  3D - velocities
+#ifndef NO_3D
+   if (allocated(u3d) .and. allocated(v3d)) then
+      call to_3d_uu(imin,jmin,imax,jmax,kmin,kmax,az, &
+                    hun,uu,vel_missing,u3d)
+      call to_3d_vv (imin,jmin,imax,jmax,kmin,kmax,az, &
+                     hvn,vv,vel_missing,v3d)
+   end if
+
+   if (allocated(u3d_destag) .and. allocated(v3d_destag)) then
+      call to_3d_vel(imin,jmin,imax,jmax,kmin,kmax,au, &
+                     hun,uu,vel_missing,u3d_destag)
+      call to_3d_vel(imin,jmin,imax,jmax,kmin,kmax,av, &
+                     hvn,vv,vel_missing,v3d_destag)
+   end if
+#endif
+#endif
+
+   return
+   end subroutine do_output_processing
+!EOC
+
+   end module output_processing
+
+!-----------------------------------------------------------------------
+! Copyright (C) 2019 - Karsten Bolding & Jorn Bruggeman (BB)           !
+!-----------------------------------------------------------------------
diff --git a/src/output/parser.F90 b/src/output/parser.F90
deleted file mode 100644
index 01efef1a9f5b1b7cc8070ab19e1e2a68023dad55..0000000000000000000000000000000000000000
--- a/src/output/parser.F90
+++ /dev/null
@@ -1,551 +0,0 @@
-#include "cppdefs.h"
-!-----------------------------------------------------------------------
-!BOP
-!
-! !MODULE:
-!
-! !INTERFACE:
-   module parser
-!
-! !DESCRIPTION:
-!  This module define a type $output_set$ defining variables - i.e.
-!  hyperslab and variable list information - for allowing a very flexible
-!  outut scheme in GETM. The actual full output specifications is in a
-!  list of $output_set$'s. Information about output-definitions is read
-!  from an ASCII-file with a relative simple format.
-!  The file contains a number of lines with a given format. Lines starting
-!  with \# or \! and empty lines are discarded. All non-discarded lines must
-!  have a forma like:
-!  \begin{verbatim}
-!   all      |n=50:end:25|elev,temp,salt
-!   test2d   |i=50:60, j=40:70:2, k=10:end,n=10:end:10|metrics,2d
-!   test3d   |i=50:60, j=40:70:2, k=10:end,n=10:end:10|3d,U,V
-!   all      |n=50:end:25|airp
-!  \end{verbatim}
-!  Each line has 3 fields separated by "|". The first field is the
-!  id for an output set. The second field is the hyperslab - in indices
-!  i,j,k,n and the last field is the variable list.\newline
-!  Note that an id can be repeated. In this case the first hyperslab
-!  definition is used and the variable list is the union. For the hyperslap
-!  field it is not necessary to specify for all indices. If a given index
-!  is not defined the minimum and maximum values are used with a stride of
-!  1. Note that the keywords $begin$ and $end$ can be used. The 3 field
-!  contains a comma-separated list of variables to include iin the
-!  output-set. The variable names must match the names of variables given
-!  in $variable_info.F90$
-!  \newline
-!  All implementation details not decided yet.
-!  \newline
-!  Work in progress and some issues still remaining.
-!  A very first version of this routine was made by Alexander Barth.
-!
-! !USES:
-   IMPLICIT NONE
-!
-! !PUBLIC DATA MEMBERS:
-
-! !DEFINED PARAMETERS
-   character, parameter :: field_sep = '|'
-   character, parameter :: var_sep = ','
-   character, parameter :: range_sep = ','
-   character, parameter :: index_sep = ':'
-   integer, parameter   :: max_length = 256
-   integer, parameter   :: max_name   = 16
-   character(len=max_name) :: metrics(8) = [ "dxc", "dyc", "dxu", "dyu", &
-                                             "dxv", "dyv", "dxx", "dyx"  ]
-   character(len=max_name) :: a2d(3) = [ "elev", "U", "V" ]
-#ifdef NO_BAROCLINIC
-   character(len=max_name) :: a3d(3) = [ "elev", "uu", "vv" ]
-#else
-   character(len=max_name) :: a3d(5) = [ "elev", "uu", "vv", "S", "T" ]
-#endif
-!
-! !REVISION HISTORY:
-!  Original author(s): Karsten Bolding
-!
-!  subset of the model domain for output
-   type output_set
-      character(len=max_length) :: id
-!      integer :: nfid
-
-      character(len=max_name), allocatable :: var_list(:)
-      integer, pointer :: varinfo_index(:)
-      integer, pointer :: varids(:)
-
-      integer :: il,ih,is
-      integer :: jl,jh,js
-      integer :: kl,kh,ks
-      integer :: nl,nh,ns
-   end type output_set
-
- ! list of all subsets
-   type(output_set), pointer :: output_list(:)
-
-!EOP
-!-----------------------------------------------------------------------
-
-contains
-
-!-----------------------------------------------------------------------
-!BOP
-! !ROUTINE: parse_output_list() -
-!
-! !INTERFACE:
-   subroutine parse_output_list(fn,imin,imax,jmin,jmax,kmin,kmax,nmin,nmax,ol)
-   IMPLICIT NONE
-!
-! !DESCRIPTION:
-!  This subroutine opens and parses the file with output-set definitions.
-!  On exit of the subroutine a list of output definitions are in the
-!  array of $output_set$-type -  $ol$.
-!
-! !INPUT PARAMETERS:
-   character(len=*), intent(in)   :: fn
-   integer, intent(in)            :: imin,imax
-   integer, intent(in)            :: jmin,jmax
-   integer, intent(in)            :: kmin,kmax
-   integer, intent(in)            :: nmin,nmax
-!
-! !OUTPUT PARAMETERS:
-   type(output_set), pointer :: ol(:)
-!
-! !LOCAL VARIABLES:
-   character(len=max_length) :: line
-   integer                   :: i,j,k=1,l,m,n
-   integer                   :: n_output_set=-1
-   integer                   :: nb_vars
-   integer                   :: Nlines
-   integer                   :: iostat
-   integer                   :: iunit = 10
-   character(len=max_length), allocatable :: id_list(:)
-   character(len=max_length), allocatable :: lines(:)
-   character(len=max_name),   allocatable :: var_names(:)
-!EOP
-!-------------------------------------------------------------------------
-!BOC
-   open(unit=iunit,file=fn,status='old')
-   read(iunit,'(A)',iostat=iostat) line
-
-! loop over all lines in file and get the number of 'valid' lines
-   Nlines = 0
-   do while (iostat == 0)
-     if (line(1:1) == '#' .or. line(1:1) == '!' .or. &
-         len_trim(line) .eq. 0) then
-       ! comments (ignore)
-      else
-         Nlines=Nlines+1
-      end if
-      read(iunit,'(A)',iostat=iostat) line
-   end do
-   allocate(id_list(Nlines))
-   allocate(lines(Nlines))
-!write(*,*) 'Nlines= ',Nlines
-
-!  Get full list of id's and read in lines with information
-   rewind (iunit)
-   read(iunit,'(A)',iostat=iostat) line
-   n=0
-   do while (iostat == 0)
-      if (line(1:1) == '#' .or. line(1:1) == '!' .or. &
-          len_trim(line) .eq. 0) then
-        ! comments (ignore)
-      else
-         i = 1
-         j = next_index(line,field_sep,i)-1
-         n=n+1
-         id_list(n) = trim(line(i:j))
-         lines(n) = line
-      end if
-      read(iunit,'(A)',iostat=iostat) line
-   end do
-   close(iunit)
-
-!  Get unique list of id's
-   n_output_set = n
-   i = n_output_set
-   do n=1,n_output_set
-      if (len_trim(id_list(n)) .gt. 0) then
-         do k=n+1,n_output_set
-            if (len_trim(id_list(k)) .gt. 0) then
-               if (trim(id_list(n)) == trim(id_list(k))) then
-                  id_list(k) = ''
-                  i=i-1
-               end if
-            end if
-         end do
-      end if
-   end do
-!write(*,*) 'Nuniq=  ',i
-   n_output_set=i
-
-!  allocate output list and copy unique id's
-   allocate(ol(n_output_set))
-   do n=1,n_output_set
-      ol(n)%id = trim(id_list(n))
-!write(*,*) trim(ol(n)%id)
-   end do
-
-!  get indices range
-   do n=1,n_output_set
-!     default: output everything
-      ol(n)%il = imin; ol(n)%ih = imax; ol(n)%is = 1
-      ol(n)%jl = jmin; ol(n)%jh = jmax; ol(n)%js = 1
-      ol(n)%kl = kmin; ol(n)%kh = kmax; ol(n)%ks = 1
-      ol(n)%nl = nmin; ol(n)%nh = nmax; ol(n)%ns = 1
-
-!     i and j are start and end index of every token on a line
-      k = 1
-      i = 1
-      j = next_index(lines(k),field_sep,i)-1
-      do while (trim(ol(n)%id) .ne. trim(lines(k)(i:j)))
-         k=k+1
-      end do
-!write(*,*) 'bb ',n,k,trim(lines(k)(i:j))
-      j = next_index(lines(k),field_sep,i)-1
-!      ol(k)%id = trim(line(i:j))
-
-!     get indices range
-      i=j+2
-      j = next_index(lines(k),field_sep,i)-1
-
-      m = next_index(lines(k),range_sep,i,j)-1
-!write (*,*) 'aa ',i,j,m
-!write (*,*) 'aa ',trim(lines(k))
-!write (*,*) 'aa ',lines(k)(i:j)
-!write (*,*) 'aa ',lines(k)(i:m)
-
-      do while (i < j)
-         l = next_index(lines(k),'=',i,m)-1
-
-!write(*,*) 'cc ',n,k
-!write(*,*) 'cc ',i,j,l,m
-!write(*,*) 'cc ',lines(k)(i:l)
-!write(*,*) 'cc ',lines(k)(i:m)
-         if (lines(k)(l:l) == 'i') then
-            i=l+2
-            call parse_range(lines(k)(i:m),ol(n)%il,ol(n)%ih,ol(n)%is,imin,imax)
-         else if (lines(k)(l:l) == 'j') then
-            i=l+2
-            call parse_range(lines(k)(i:m),ol(n)%jl,ol(n)%jh,ol(n)%js,jmin,jmax)
-         else if (lines(k)(l:l) == 'k') then
-            i=l+2
-            call parse_range(lines(k)(i:m),ol(n)%kl,ol(n)%kh,ol(n)%ks,kmin,kmax)
-         else if (lines(k)(l:l) == 'n') then
-            i=l+2
-!write(*,*) 'dd ',lines(k)(i:m)
-            call parse_range(lines(k)(i:m),ol(n)%nl,ol(n)%nh,ol(n)%ns,nmin,nmax)
-         end if
-
-         i = m+2
-         m = next_index(lines(k),range_sep,i,j)-1
-      end do
-   end do
-
-!  concatenate variable lists from dublicated id's
-   do n=1,n_output_set
-
-!     i and j are start and end index of every token on a line
-      do k=1,Nlines
-         i = 1
-         j = next_index(lines(k),field_sep,i)-1
-         if (trim(ol(n)%id) == trim(lines(k)(i:j))) then
-!           skip to variable list
-            j = next_index(lines(k),field_sep,i)
-            i=j+1
-            j = next_index(lines(k),field_sep,i)
-            i=j+1
-!write(*,*) n,k,i,j
-!write(*,*) lines(k)(i:)
-            if (k .eq. n) then
-               line = lines(k)(i:)
-            else
-               line = trim(line)//','//lines(k)(i:)
-            end if
-         end if
-      end do
-!write(*,*) 'vars= ',trim(line)
-
-!     get rid of duplicate variable names
-      nb_vars = count_chars(line,var_sep) + 1
-      allocate(var_names(nb_vars))
-      i = 1
-      j = next_index(line,var_sep,i)-1
-      do l=1,nb_vars
-         var_names(l)=trim(line(i:j))
-         i=j+2
-         j = next_index(line,var_sep,i)-1
-      end do
-#if 0
-do l=1,nb_vars
-  write(*,*) 'aa ',l,trim(var_names(l))
-end do
-#endif
-      m=nb_vars
-      do l=1,nb_vars
-         do k=l+1,nb_vars
-            if (len_trim(var_names(l)) .gt. 0) then
-               if(trim(var_names(l)) == trim(var_names(k))) then
-                  var_names(k) = ''
-                  m=m-1
-               end if
-            end if
-         end do
-      end do
-      nb_vars=m
-#if 0
-i = 0
-do l=1,size(var_names)
-   if (len_trim(var_names(l)) .gt. 0) then
-      i = i+1
-      write(*,*) 'bb ',l,trim(var_names(i))
-   end if
-end do
-#endif
-!     expand predefined definitions of variable lists
-      m=nb_vars
-      do l=1,nb_vars
-         if (trim(var_names(l)) .eq. 'metrics') m=m+size(metrics)-1
-         if (trim(var_names(l)) .eq. '2d')      m=m+size(a2d)-1
-         if (trim(var_names(l)) .eq. '3d')      m=m+size(a3d)-1
-         if (trim(var_names(l)) .eq. 'turb')    m=m+5-1
-      end do
-      nb_vars=m
-!write(*,*) 'm= ',m
-!write(*,*) 'n= ',n
-
-!     allocate(ol(k)%varinfo_index(nb_vars))
-      allocate(ol(n)%var_list(nb_vars))
-      i = 0
-      do l=1,size(var_names)
-         if (len_trim(var_names(l)) .gt. 0) then
-            if (trim(var_names(l)) == "metrics") then
-               do j=1,size(metrics)
-                  i = i+1
-                  ol(n)%var_list(i) = metrics(j)
-               end do
-            else if (trim(var_names(l)) == "2d") then
-               do j=1,size(a2d)
-                  i = i+1
-                  ol(n)%var_list(i) = a2d(j)
-               end do
-            else if (trim(var_names(l)) == "3d") then
-               do j=1,size(a3d)
-                  i = i+1
-                  ol(n)%var_list(i) = a3d(j)
-               end do
-            else
-               i = i+1
-               ol(n)%var_list(i) = trim(var_names(l))
-            end if
-          end if
-      end do
-      deallocate(var_names)
-   end do
-   end subroutine parse_output_list
-!EOC
-
-!-----------------------------------------------------------------------
-!BOP
-! !IROUTINE: print_output_list -
-!
-! !DESCRIPTION:
-!
-! !INTERFACE:
-   subroutine print_output_list(ol)
-   IMPLICIT NONE
-!
-! !INPUT PARAMETERS:
-   type(output_set), pointer      :: ol(:)
-!
-! !LOCAL VARIABLES:
-  integer          :: k,i
-!EOP
-!-------------------------------------------------------------------------
-!BOC
-  do k=1,size(ol)
-    write(6,'(A,A)') 'id: ',trim(ol(k)%id)
-    do i=1,size(ol(k)%var_list)
-!      write(6,'(A,A)') '  variable: ', trim(varinfo_list(ol(k)%varinfo_index(i))%name)
-      write(6,'(A,A)') '  variable: ', trim((ol(k)%var_list(i)))
-    end do
-    write(6,'(A,3I10)') '  i range: ',ol(k)%il,ol(k)%ih,ol(k)%is
-    write(6,'(A,3I10)') '  j range: ',ol(k)%jl,ol(k)%jh,ol(k)%js
-    write(6,'(A,3I10)') '  k range: ',ol(k)%kl,ol(k)%kh,ol(k)%ks
-    write(6,'(A,3I10)') '  n range: ',ol(k)%nl,ol(k)%nh,ol(k)%ns
-  end do
- end subroutine print_output_list
-!EOC
-
-!-----------------------------------------------------------------------
-!BOP
-! !IROUTINE: deallocate_output_list -
-!
-! !DESCRIPTION:
-!
-! !INTERFACE:
-   subroutine deallocate_output_list(ol)
-   IMPLICIT NONE
-!
-! !INPUT/OUTPUT PARAMETERS:
-  type(output_set), pointer       :: ol(:)
-!
-! !LOCAL VARIABLES:
-   integer         :: k
-!EOP
-!-------------------------------------------------------------------------
-!BOC
-  do k=1,size(ol)
-    deallocate(ol(k)%varinfo_index)
-    deallocate(ol(k)%varids)
-  end do
-  deallocate(ol)
-  end subroutine deallocate_output_list
-!EOC
-
-!-----------------------------------------------------------------------
-!BOP
-! !IROUTINE: next_index -
-!
-! !INTERFACE:
-   function next_index(s,sep,start,finish)
-   IMPLICIT NONE
-!
-! !DESCRIPTION:
-!
-! !INPUT PARAMETERS:
-   character(len=*), intent(in)   :: s
-   character, intent(in)          :: sep
-   integer, intent(in)            :: start
-   integer, intent(in), optional  :: finish
-!
-! !LOCAL VARIABLES:
-  integer          :: next_index
-   integer         :: i,f
-!EOP
-!-------------------------------------------------------------------------
-!BOC
-   if (present(finish)) then
-      f = finish
-   else
-      f = len_trim(s)
-   end if
-   do i=start,f
-      if (s(i:i) == sep) then
-         next_index = i
-         return
-      end if
-   end do
-   next_index = f+1
-   end function next_index
-!EOC
-
-!-----------------------------------------------------------------------
-!BOP
-! !IROUTINE: count_chars() -
-!
-! !INTERFACE:
-   function count_chars(s,sep)
-   IMPLICIT NONE
-!
-! !DESCRIPTION:
-!
-! !INPUT PARAMETERS:
-   character(len=*), intent(in)   :: s
-   character, intent(in)          :: sep
-!
-! !LOCAL VARIABLES:
-   integer         :: i
-   integer         :: count_chars
-!EOP
-!-------------------------------------------------------------------------
-!BOC
-   count_chars = 0
-   do i=1,len_trim(s)
-     if (s(i:i) == sep) then
-       count_chars = count_chars+1
-     end if
-   end do
-   end function count_chars
-!EOC
-
-!-----------------------------------------------------------------------
-!BOP
-! !IROUTINE: parse_index() -
-!
-! !DESCRIPTION:
-!
-! !INTERFACE:
-   function parse_index(s,min,max)
-   IMPLICIT NONE
-!
-! !INPUT PARAMETERS:
-   character(len=*), intent(in)   :: s
-   integer, intent(in)            :: min,max
-!
-! !LOCAL VARIABLES:
-   integer parse_index
-!EOP
-!-------------------------------------------------------------------------
-!BOC
-   if (s == 'end') then
-     parse_index = max
-   elseif (s == 'begin') then
-     parse_index = min
-   else
-     read(s,*) parse_index
-   end if
-   end function parse_index
-!EOC
-
-!-----------------------------------------------------------------------
-!BOP
-! !IROUTINE: parse_range() -
-!
-! !INTERFACE:
-   subroutine parse_range(s,first,last,step,min,max)
-   IMPLICIT NONE
-!
-! !DESCRIPTION:
-!
-! !INPUT PARAMETERS:
-   character(len=*), intent(in)   :: s
-   integer, intent(in)            :: min,max
-!
-! !OUTPUT PARAMETERS:
-   integer, intent(out)           :: first,last,step
-!
-! !LOCAL VARIABLES:
-   integer         :: i,j
-!EOP
-!-------------------------------------------------------------------------
-!BOC
-   i = 1
-   j = next_index(s,index_sep,i)-1
-   first = parse_index(s(i:j),min,max)
-   i=j+2
-   j = next_index(s,index_sep,i)-1
-   if (len_trim(s(i:j)) == 0) then
-      last = first
-      step = 1
-   else
-      last = parse_index(s(i:j),min,max)
-
-      i=j+2
-      j = next_index(s,index_sep,i)-1
-
-      if (len_trim(s(i:j)) == 0) then
-         step = 1
-      else
-         step = parse_index(s(i:j),min,max)
-      end if
-   end if
-   end subroutine parse_range
-!EOC
-
-!-----------------------------------------------------------------------
-
-   end module parser
-
-!-----------------------------------------------------------------------
-! Copyright (C) 2009 - Karsten Bolding (BB)                            !
-!-----------------------------------------------------------------------
diff --git a/src/output/test_parser.F90 b/src/output/test_parser.F90
deleted file mode 100644
index 326b6b8c99b319430321fa0ca274968114ac5929..0000000000000000000000000000000000000000
--- a/src/output/test_parser.F90
+++ /dev/null
@@ -1,25 +0,0 @@
-   program test_parser
-!
-! !DESCRIPTION:
-!  program to test the parser routine in parser.F90
-!  To execute:
-!  make test_parser
-!
-! !USES:
-   use parser
-   implicit none
-! !LOCAL VARIABLES
-   integer :: imin=1,imax=100
-   integer :: jmin=100,jmax=200
-   integer :: kmin=1,kmax=15
-   integer :: nmin=1,nmax=10000
-!EOP
-!-----------------------------------------------------------------------
-!BOC
-   call parse_output_list('out.dat',imin,imax,jmin,jmax,kmin,kmax,nmin,nmax, &
-                         output_list)
-   call print_output_list(output_list)
-
-   end program test_parser
-!EOC
-
diff --git a/src/output/test_varinfo.F90 b/src/output/test_varinfo.F90
deleted file mode 100644
index 86c3f1c056da95e0a2b4d335baa2d4d069d18006..0000000000000000000000000000000000000000
--- a/src/output/test_varinfo.F90
+++ /dev/null
@@ -1,19 +0,0 @@
-   program test_varinfo
-!
-! !DESCRIPTION:
-!  Test program for testing routines in variable_info.F90
-!  To execute:
-!  make test_varinfo
-!
-! !USES:
-   use variable_info, only: init_var_info, print_var_info
-   IMPLICIT NONE
-!
-!EOP
-!-----------------------------------------------------------------------
-!BOC
-   call init_var_info()
-   call print_var_info()
-
-   end program test_varinfo
-!EOC
diff --git a/src/output/variable_info.F90 b/src/output/variable_info.F90
deleted file mode 100644
index fd07baf0dd6695cb7cff92ed47862bc45395b100..0000000000000000000000000000000000000000
--- a/src/output/variable_info.F90
+++ /dev/null
@@ -1,569 +0,0 @@
-!#ifdef NETCDF4
-#include "cppdefs.h"
-!-----------------------------------------------------------------------
-!BOP
-! !MODULE:  variable_info - metadata for all variables that can be saved.
-!
-! !INTERFACE:
-   module variable_info
-!
-! !DESCRIPTION:
-!  This module will be used by the GETM output system to obtain information
-!  about variables. Multiple output definitions i.e. hyperslabs and
-!  variable-lists shall all use the information contained here. The key
-!  to the meta-data is the $name$ i.e. search the array until agreement
-!  between $varinfo$->$name$ and variable name searched for.
-!  This module provides meta data information about all variables GETM
-!  that be saved in output files. The specific data for a given variable
-!  is contained in a defined type $varinfo$. $varinfo$ contains a number
-!  of fields including a unique $id$, a $name$, information about the
-!  dimensions of the variable and additional fields that can be used by
-!  NetCDF to give COARDS conforming meta data information.
-!  The $varinfo$ also contains pointer types that can be used to point to
-!  the Fortran variable actually holding the data. For a re-useable type
-!  definition it is necessary to define pointers to 1D, 2D and 3D
-!  variables.
-!  The $varinfo\_list$ can be expanded to also include compound variables
-!  i.e. variables that are functions of other variables. In this case the
-!  unique $id$ shall be used to actually select which variables and what
-!  operations should be done - like:
-!  \begin{verbatim}
-!  select case (varinfo_list(1)%id)
-!     case(?)
-!     case(_SALT_FLUX_)
-!        data = uu*S
-!     case(?)
-!  \end{verbatim}
-!  \newline
-!  The final implementation details not decided yet.\newline
-!  Work in progress.
-!
-! !USES:
-   use exceptions
-#if 0
-!   use domain
-!   use variables_2d
-!   use variables_3d
-!   use meteo
-#endif
-   IMPLICIT NONE
-!
-   private
-!
-! !PUBLIC DATA MEMBERS
-   public init_var_info, print_var_info
-
-   integer, parameter           :: max_length=255
-   type, public                 :: varinfo
-      integer                   :: id=-1
-      character(len=max_length) :: name=''
-      integer                   :: ndims=-1
-      character(len=max_length) :: long_name=''
-      character(len=max_length) :: units=''
-      REALTYPE                  :: mv
-      REALTYPE                  :: vr(2)
-      REALTYPE, pointer         :: data_0d        => null()
-      REALTYPE, pointer         :: data_1d(:)     => null()
-      REALTYPE, pointer         :: data_2d(:,:)   => null()
-      REALTYPE, pointer         :: data_3d(:,:,:) => null()
-!      REALTYPE, pointer         :: save_data_2d(:,:)   => null()
-   end type varinfo
-
-   type(varinfo)                :: varinfo_list(40)
-!
-! !REVISION HISTORY:
-!  Original author(s): Karsten Bolding & Hannes Rennau
-!
-!EOP
-!-----------------------------------------------------------------------
-
-contains
-
-!-----------------------------------------------------------------------
-!BOP
-!
-! !IROUTINE: Optain info on the variable given by var_name
-!
-! !INTERFACE:
-   subroutine init_var_info()
-!
-! !DESCRIPTION:
-!  Here a list meta-data for all possible GETM outputr variables is
-!  constructed. The information is put in an array.
-!
-   IMPLICIT NONE
-#if 1
-#define STATIC
-#ifdef STATIC
-   REALTYPE, target :: z(0:11,-10:10)
-#else
-#endif
-#undef STATIC
-#endif
-
-!
-! !LOCAL VARIABLES
-#define _D_      1
-#define _DU_     2
-#define _DV_     3
-#define _AIRP_   4
-#define _U10_    5
-#define _V10_    6
-#define _T2_     7
-#define _HUM_    8
-#define _TCC_    9
-#define _TAUSX_ 10
-#define _TAUSY_ 11
-#define _SWR_   12
-#define _SHF_   13
-#define _ELEV_  14
-#define _U_     15
-#define _V_     16
-#define _RES_U_ 17
-#define _RES_V_ 18
-#define _UU_    19
-#define _VV_    20
-#define _HCC_   21
-#define _HO_    22
-#define _HN_    23
-#define _HU_    24
-#define _HV_    25
-#define _WW_    26
-#define _T_     27
-#define _S_     28
-#define _TKE_   29
-#define _NUM_   30
-#define _NUH_   31
-#define _EPS_   32
-#define _SS_    33
-#define _NN_    34
-#define _SPM_   35
-   integer         :: n=0
-!EOP
-!-----------------------------------------------------------------------
-!BOC
-
-!   call random_number(z)
-!   n=0
-
-   n=n+1
-!  Domain related variables variables
-   varinfo_list(n)%id=          _D_
-   varinfo_list(n)%name=        'D'
-   varinfo_list(n)%ndims=       2
-   varinfo_list(n)%long_name=   'total water depth'
-   varinfo_list(n)%units=       'm'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ -15.0, 5000.0 ]
-!   varinfo_list(n)%data_2d=>    D
-
-   n=n+1
-   varinfo_list(n)%id=          _DU_
-   varinfo_list(n)%name=        'DU'
-   varinfo_list(n)%ndims=       2
-   varinfo_list(n)%long_name=   'total water depth - U-points'
-   varinfo_list(n)%units=       'm'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ -15.0, 5000.0 ]
-!   varinfo_list(n)%data_2d=>    DU
-
-   n=n+1
-   varinfo_list(n)%id=          _DV_
-   varinfo_list(n)%name=        'DV'
-   varinfo_list(n)%ndims=       2
-   varinfo_list(n)%long_name=   'total water depth - V-points'
-   varinfo_list(n)%units=       'm'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ -15.0, 5000.0 ]
-!   varinfo_list(n)%data_2d=>    DV
-
-!  meteo variables
-   n=n+1
-   varinfo_list(n)%id=          _AIRP_
-   varinfo_list(n)%name=        'airp'
-   varinfo_list(n)%ndims=       2
-   varinfo_list(n)%long_name=   'air pressure'
-   varinfo_list(n)%units='Pascal'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ 90.e3, 110.e3 ]
-!   varinfo_list(n)%data_2d=>    airp
-
-   n=n+1
-   varinfo_list(n)%id=          _U10_
-   varinfo_list(n)%name=        'u10'
-   varinfo_list(n)%ndims=       2
-   varinfo_list(n)%long_name=   'eastward windspeed (10m)'
-   varinfo_list(n)%units='m/s'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ -50.0, 50.0 ]
-!   varinfo_list(n)%data_2d=>    u10
-
-   n=n+1
-   varinfo_list(n)%id=          _V10_
-   varinfo_list(n)%name=        'v10'
-   varinfo_list(n)%ndims=       2
-   varinfo_list(n)%long_name=   'northward windspeed (10m)'
-   varinfo_list(n)%units='m/s'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ -50.0, 50.0 ]
-!   varinfo_list(n)%data_2d=>    v10
-
-   n=n+1
-   varinfo_list(n)%id=          _T2_
-   varinfo_list(n)%name=        't2'
-   varinfo_list(n)%ndims=       2
-   varinfo_list(n)%long_name=   'temperature (2m)'
-   varinfo_list(n)%units='Kelvin'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ 0.0, 325.0 ]
-!   varinfo_list(n)%data_2d=>    t2
-
-   n=n+1
-   varinfo_list(n)%id=          _HUM_
-   varinfo_list(n)%name=        'hum'
-   varinfo_list(n)%ndims=       2
-   varinfo_list(n)%long_name=   'humidity'
-   varinfo_list(n)%units='kg/kg'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ 0.0, 100.0 ]
-!   varinfo_list(n)%data_2d=>    hum
-
-   n=n+1
-   varinfo_list(n)%id=          _TCC_
-   varinfo_list(n)%name=        'tcc'
-   varinfo_list(n)%ndims=       2
-   varinfo_list(n)%long_name=   'total cloud cover'
-   varinfo_list(n)%units=''
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ 0.0, 1.0 ]
-!   varinfo_list(n)%data_2d=>    tcc
-
-   n=n+1
-   varinfo_list(n)%id=          _TAUSX_
-   varinfo_list(n)%name=        'tausx'
-   varinfo_list(n)%ndims=       2
-   varinfo_list(n)%long_name=   'surface stress - x'
-   varinfo_list(n)%units='N/m2'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ -1.0, 1.0 ]
-!   varinfo_list(n)%data_2d=>    tausx
-
-   n=n+1
-   varinfo_list(n)%id=          _TAUSY_
-   varinfo_list(n)%name=        'tausy'
-   varinfo_list(n)%ndims=       2
-   varinfo_list(n)%long_name=   'surface stress - y'
-   varinfo_list(n)%units='N/m2'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ -1.0, 1.0 ]
-!   varinfo_list(n)%data_2d=>    tausy
-
-   n=n+1
-   varinfo_list(n)%id=          _SWR_
-   varinfo_list(n)%name=        'swr'
-   varinfo_list(n)%ndims=       2
-   varinfo_list(n)%long_name=   'short wave radiation'
-   varinfo_list(n)%units='W/m2'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ 0.0, 1500.0 ]
-!   varinfo_list(n)%data_2d=>    swr
-
-   n=n+1
-   varinfo_list(n)%id=          _SHF_
-   varinfo_list(n)%name=        'shf'
-   varinfo_list(n)%ndims=       2
-   varinfo_list(n)%long_name=   'surface heat fluxes'
-   varinfo_list(n)%units='W/m2'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ -1000.0, 1000.0 ]
-!   varinfo_list(n)%data_2d=>    shf
-
-!  2D variables
-   n=n+1
-   varinfo_list(n)%id=          _ELEV_
-   varinfo_list(n)%name=        'elev'
-   varinfo_list(n)%ndims=       2
-   varinfo_list(n)%long_name=   'sea surface elevation'
-   varinfo_list(n)%units=       'm'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ -15.0, 15.0 ]
-!   varinfo_list(n)%data_2d=>    elev
-
-   n=n+1
-   varinfo_list(n)%id=          _U_
-   varinfo_list(n)%name=        'u'
-   varinfo_list(n)%ndims=       2
-   varinfo_list(n)%long_name=   'int. zonal vel.'
-   varinfo_list(n)%units=       'm/s'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ -3, 3 ]
-!   varinfo_list(n)%data_2d=>    u
-
-   n=n+1
-   varinfo_list(n)%id=          _V_
-   varinfo_list(n)%name=        'v'
-   varinfo_list(n)%ndims=       2
-   varinfo_list(n)%long_name=   'int. meridional vel.'
-   varinfo_list(n)%units=       'm/s'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ -3, 3 ]
-!   varinfo_list(n)%data_2d=>    v
-
-   n=n+1
-   varinfo_list(n)%id=          _RES_U_
-   varinfo_list(n)%name=        'res_u'
-   varinfo_list(n)%ndims=       2
-   varinfo_list(n)%long_name=   'residual of U'
-   varinfo_list(n)%units=       'm/s'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ -3, 3 ]
-!   varinfo_list(n)%data_2d=>    res_u
-
-   n=n+1
-   varinfo_list(n)%id=          _RES_V_
-   varinfo_list(n)%name=        'res_v'
-   varinfo_list(n)%ndims=       2
-   varinfo_list(n)%long_name=   'residual of V'
-   varinfo_list(n)%units=       'm/s'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ -3, 3 ]
-!   varinfo_list(n)%data_2d=>    res_v
-
-
-!  3D variables
-   n=n+1
-   varinfo_list(n)%id=          _UU_
-   varinfo_list(n)%name=        'uu'
-   varinfo_list(n)%ndims=       3
-   varinfo_list(n)%long_name=   'zonal velocity'
-   varinfo_list(n)%units=       'm/s'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ -5.0, 5.0 ]
-!   varinfo_list(n)%data_3d=>    uu
-
-   n=n+1
-   varinfo_list(n)%id=          _VV_
-   varinfo_list(n)%name=        'vv'
-   varinfo_list(n)%ndims=       3
-   varinfo_list(n)%long_name=   'meridional velocity'
-   varinfo_list(n)%units=       'm/s'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ -5.0, 5.0 ]
-!   varinfo_list(n)%data_3d=>    vv
-
-   n=n+1
-   varinfo_list(n)%id=          _HCC_
-   varinfo_list(n)%name=        'hcc'
-   varinfo_list(n)%ndims=       3
-   varinfo_list(n)%long_name=   'hydrostatic consistency criterion'
-   varinfo_list(n)%units=       ''
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ 0.0, 1.0 ]
-!   varinfo_list(n)%data_3d=>    hcc
-
-   n=n+1
-   varinfo_list(n)%id=          _HO_
-   varinfo_list(n)%name=        'ho'
-   varinfo_list(n)%ndims=       3
-   varinfo_list(n)%long_name=   'old box height - h-column'
-   varinfo_list(n)%units=       'm'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ 0.0, 1000.0 ]
-!   varinfo_list(n)%data_3d=>    ho
-
-   n=n+1
-   varinfo_list(n)%id=          _HN_
-   varinfo_list(n)%name=        'hn'
-   varinfo_list(n)%ndims=       3
-   varinfo_list(n)%long_name=   'new box height - h-column'
-   varinfo_list(n)%units=       'm'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ 0.0, 1000.0 ]
-!   varinfo_list(n)%data_3d=>    hn
-
-   n=n+1
-   varinfo_list(n)%id=          _HU_
-   varinfo_list(n)%name=        'ho'
-   varinfo_list(n)%ndims=       3
-   varinfo_list(n)%long_name=   'box height - u-column'
-   varinfo_list(n)%units=       'm'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ 0.0, 1000.0 ]
-!   varinfo_list(n)%data_3d=>    hu
-
-   n=n+1
-   varinfo_list(n)%id=          _HV_
-   varinfo_list(n)%name=        'hv'
-   varinfo_list(n)%ndims=       3
-   varinfo_list(n)%long_name=   'new box height - v-column'
-   varinfo_list(n)%units=       'm'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ 0.0, 1000.0 ]
-!   varinfo_list(n)%data_3d=>    hv
-
-   n=n+1
-   varinfo_list(n)%id=          _WW_
-   varinfo_list(n)%name=        'ww'
-   varinfo_list(n)%ndims=       3
-   varinfo_list(n)%long_name=   'vertical velocity'
-   varinfo_list(n)%units=       'm/s'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ -3.0, 3.0 ]
-!   varinfo_list(n)%data_3d=>    ww
-
-   n=n+1
-   varinfo_list(n)%id=          _T_
-   varinfo_list(n)%name=        'T'
-   varinfo_list(n)%ndims=       3
-   varinfo_list(n)%long_name=   'temperature'
-   varinfo_list(n)%units=       'degC'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ -2.0, 40.0 ]
-!   varinfo_list(n)%data_3d=>    T
-
-   n=n+1
-   varinfo_list(n)%id=          _S_
-   varinfo_list(n)%name=        'S'
-   varinfo_list(n)%ndims=       3
-   varinfo_list(n)%long_name=   'salinity'
-   varinfo_list(n)%units=       'psu'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ 0.0, 40.0 ]
-!   varinfo_list(n)%data_3d=>    S
-
-   n=n+1
-   varinfo_list(n)%id=          _TKE_
-   varinfo_list(n)%name=        'tke'
-   varinfo_list(n)%ndims=       3
-   varinfo_list(n)%long_name=   'turbulent kinetic energy'
-   varinfo_list(n)%units=       'm2/s2'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ 0.0, 0.2 ]
-!   varinfo_list(n)%data_3d=>    tke
-
-   n=n+1
-   varinfo_list(n)%id=          _NUM_
-   varinfo_list(n)%name=        'num'
-   varinfo_list(n)%ndims=       3
-   varinfo_list(n)%long_name=   'viscosity'
-   varinfo_list(n)%units=       'm2/s'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ 0.0, 0.2 ]
-!   varinfo_list(n)%data_3d=>    num
-
-   n=n+1
-   varinfo_list(n)%id=          _NUH_
-   varinfo_list(n)%name=        'nuh'
-   varinfo_list(n)%ndims=       3
-   varinfo_list(n)%long_name=   'diffusivity'
-   varinfo_list(n)%units=       'm2/s'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ 0.0, 0.2 ]
-!   varinfo_list(n)%data_3d=>    nuh
-
-   n=n+1
-   varinfo_list(n)%id=          _EPS_
-   varinfo_list(n)%name=        'eps'
-   varinfo_list(n)%ndims=       3
-   varinfo_list(n)%long_name=   'dissipation'
-   varinfo_list(n)%units=       'm2/s3'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ 0.0, 0.2 ]
-!   varinfo_list(n)%data_3d=>    eps
-
-   n=n+1
-   varinfo_list(n)%id=          _SS_
-   varinfo_list(n)%name=        'SS'
-   varinfo_list(n)%ndims=       3
-   varinfo_list(n)%long_name=   'shear-frequency squared'
-   varinfo_list(n)%units=       's-2'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ -3.0, 3.0 ]
-!   varinfo_list(n)%data_3d=>    SS
-
-   n=n+1
-   varinfo_list(n)%id=          _NN_
-   varinfo_list(n)%name=        'NN'
-   varinfo_list(n)%ndims=       3
-   varinfo_list(n)%long_name=   'BVF squared'
-   varinfo_list(n)%units=       's-2'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ -3.0, 3.0 ]
-!   varinfo_list(n)%data_3d=>    NN
-
-   n=n+1
-   varinfo_list(n)%id=          _SPM_
-   varinfo_list(n)%name=        'spm'
-   varinfo_list(n)%ndims=       3
-   varinfo_list(n)%long_name=   'suspended matter concentration'
-   varinfo_list(n)%units=       'kg/m-3'
-   varinfo_list(n)%mv=          -9999.0
-   varinfo_list(n)%vr=          [ -3.0, 3.0 ]
-!   varinfo_list(n)%data_3d=>    spm
-
-   return
-   end subroutine init_var_info
-!EOC
-
-!-----------------------------------------------------------------------
-!BOP
-!
-! !IROUTINE: Optainn info on the variable given by var_name
-!
-! !INTERFACE:
-   subroutine print_var_info()
-!
-! !DESCRIPTION:
-!  Just prints out the meta-data for all variables in $varinfo_list$.
-!  For testing and debugging purposes mainly.
-!
-   IMPLICIT NONE
-!
-! !LOCAL VARIABLES:
-   integer          :: n
-!EOP
-!-----------------------------------------------------------------------
-!BOC
-
-   LEVEL1 'print_var_info() '
-   do n=1,size(varinfo_list)
-      if (varinfo_list(n)%id .gt. 0) then
-      LEVEL3 '----------------------------------------------------------'
-      LEVEL3 'variable id:    ',varinfo_list(n)%id
-      LEVEL3 'variable name:  ',trim(varinfo_list(n)%name)
-      LEVEL3 'number of dims  ',varinfo_list(n)%ndims
-      LEVEL3 'long_name:      ',trim(varinfo_list(n)%long_name)
-      LEVEL3 'units:          ',trim(varinfo_list(n)%units)
-!   LEVEL3 'fill_value=     ',fv
-      LEVEL3 'missing_values= ',varinfo_list(n)%mv
-      LEVEL3 'lower limit=    ',varinfo_list(n)%vr(1)
-      LEVEL3 'upper limit=    ',varinfo_list(n)%vr(2)
-      LEVEL3 'data=         ',varinfo_list(n)%data_2d
-!      LEVEL3 'save_data=    ',varinfo_list(n)%save_data_2d
-      end if
-   end do
-
-   do n=1,size(varinfo_list)
-      if (varinfo_list(n)%id .gt. 0) then
-         select case (varinfo_list(1)%id)
-            case(_D_)
-               LEVEL3 'var ',n,' is ',trim(varinfo_list(n)%long_name)
-            case(_AIRP_)
-               LEVEL3 'var ',n,' is ',trim(varinfo_list(n)%long_name)
-            case(_ELEV_)
-            case(_UU_)
-            case default
-         end select
-      end if
-   end do
-
-   return
-   end subroutine print_var_info
-!EOC
-
-   end module variable_info
-!#endif
-!-----------------------------------------------------------------------
-! Copyright (C) 2009 - Karsten Bolding & Haness Rennau (BB)            !
-!-----------------------------------------------------------------------
-