From f92439c9dff6763dd5e6afded59d29d9528f0c8d Mon Sep 17 00:00:00 2001
From: JakobDeutloff <50237419+JakobDeutloff@users.noreply.github.com>
Date: Tue, 19 Apr 2022 18:46:01 +0100
Subject: [PATCH] Timeframe added for TPhiS plots, os directory changes
 improved, exp_heat added

---
 .../Code/build_config_files/build_config_1.py |  14 +-
 .../Code/build_config_files/build_config_2.py |  14 +-
 .../Code/build_config_files/build_config_3.py |  14 +-
 .../Code/build_config_files/build_config_4.py |  14 +-
 .../build_config_files/build_config_MOSAiC.py | 130 ---
 .../build_config_initial_state.py             |  14 +-
 .../build_config_mosaic1.py                   |  18 +-
 Python_files/Code/plotscripts/plot_TPhiS.py   |  93 ++-
 Python_files/Plots/.gitkeep                   |   0
 Run_specifics/config.json                     |  41 +-
 Run_specifics/description.txt                 |   2 +-
 mo_data.f90                                   |   1 +
 mo_grav_drain.f90                             |   3 +-
 mo_grotz.f90                                  |  12 +-
 mo_output.f90                                 | 780 +++++++++---------
 15 files changed, 511 insertions(+), 639 deletions(-)
 delete mode 100644 Python_files/Code/build_config_files/build_config_MOSAiC.py
 create mode 100644 Python_files/Plots/.gitkeep

diff --git a/Python_files/Code/build_config_files/build_config_1.py b/Python_files/Code/build_config_files/build_config_1.py
index 85d241f..0d4b94e 100644
--- a/Python_files/Code/build_config_files/build_config_1.py
+++ b/Python_files/Code/build_config_files/build_config_1.py
@@ -2,15 +2,10 @@ import json
 from func_bound_values import *
 import os
 
-# set wd to build_config_files folder
-wd = os.getcwd()
-# check if already in correct wd
-if wd[-18:] == 'build_config_files':
-    pass
-# if not, change to build_config_files
-else:
-    path = wd + '/Python_files/Code/build_config_files'
-    os.chdir(path)
+# set wd to the directory of the script
+abspath = os.path.abspath(__file__)
+dname = os.path.dirname(abspath)
+os.chdir(dname)
 
 # set constants
 rho_l = 1028
@@ -34,6 +29,7 @@ file.close()
 # **********************************************************************************************************************
 # Initial layer thickness, timestep, output time and simulation time
 # **********************************************************************************************************************
+config['start_time'] = '2022-01-01 00:00:00'  # start time in YYYY-mm-dd HH:MM:SS
 config['dt'] = 1  # time increment [s]]
 config['time'] = 0.0  # initial value of time [s]
 config['time_out'] = 3600  # time between outputs [s]
diff --git a/Python_files/Code/build_config_files/build_config_2.py b/Python_files/Code/build_config_files/build_config_2.py
index 14bda5a..20a52b8 100644
--- a/Python_files/Code/build_config_files/build_config_2.py
+++ b/Python_files/Code/build_config_files/build_config_2.py
@@ -2,15 +2,10 @@ import json
 from func_bound_values import *
 import os
 
-# set wd to build_config_files folder
-wd = os.getcwd()
-# check if already in correct wd
-if wd[-18:] == 'build_config_files':
-    pass
-# if not, change to build_config_files
-else:
-    path = wd + '/Python_files/Code/build_config_files'
-    os.chdir(path)
+# set wd to the directory of the script
+abspath = os.path.abspath(__file__)
+dname = os.path.dirname(abspath)
+os.chdir(dname)
 
 # set constants
 rho_l = 1028
@@ -34,6 +29,7 @@ file.close()
 # **********************************************************************************************************************
 # Initial layer thickness, timestep, output time and simulation time
 # **********************************************************************************************************************
+config['start_time'] = '2022-01-01 00:00:00'  # start time in YYYY-mm-dd HH:MM:SS
 config['dt'] = 30  # time increment [s]]
 config['time'] = 0.0  # initial value of time [s]
 config['time_out'] = 3600 * 6 # time between outputs [s]
diff --git a/Python_files/Code/build_config_files/build_config_3.py b/Python_files/Code/build_config_files/build_config_3.py
index 2502bd0..c1d2943 100644
--- a/Python_files/Code/build_config_files/build_config_3.py
+++ b/Python_files/Code/build_config_files/build_config_3.py
@@ -2,15 +2,10 @@ import json
 from func_bound_values import *
 import os
 
-# set wd to build_config_files folder
-wd = os.getcwd()
-# check if already in correct wd
-if wd[-18:] == 'build_config_files':
-    pass
-# if not, change to build_config_files
-else:
-    path = wd + '/Python_files/Code/build_config_files'
-    os.chdir(path)
+# set wd to the directory of the script
+abspath = os.path.abspath(__file__)
+dname = os.path.dirname(abspath)
+os.chdir(dname)
 
 # set constants
 rho_l = 1028
@@ -34,6 +29,7 @@ file.close()
 # **********************************************************************************************************************
 # Initial layer thickness, timestep, output time and simulation time
 # **********************************************************************************************************************
+config['start_time'] = '2022-01-01 00:00:00'  # start time in YYYY-mm-dd HH:MM:SS
 config['dt'] = 60  # time increment [s]]
 config['time'] = 0.0  # initial value of time [s]
 config['time_out'] = 86400 * 3.5 # time between outputs [s]
diff --git a/Python_files/Code/build_config_files/build_config_4.py b/Python_files/Code/build_config_files/build_config_4.py
index 390eb2d..88142c6 100644
--- a/Python_files/Code/build_config_files/build_config_4.py
+++ b/Python_files/Code/build_config_files/build_config_4.py
@@ -2,15 +2,10 @@ import json
 from func_bound_values import *
 import os
 
-# set wd to build_config_files folder
-wd = os.getcwd()
-# check if already in correct wd
-if wd[-18:] == 'build_config_files':
-    pass
-# if not, change to build_config_files
-else:
-    path = wd + '/Python_files/Code/build_config_files'
-    os.chdir(path)
+# set wd to the directory of the script
+abspath = os.path.abspath(__file__)
+dname = os.path.dirname(abspath)
+os.chdir(dname)
 
 
 # set constants
@@ -35,6 +30,7 @@ file.close()
 # **********************************************************************************************************************
 # Initial layer thickness, timestep, output time and simulation time
 # **********************************************************************************************************************
+config['start_time'] = '2022-01-01 00:00:00'  # start time in YYYY-mm-dd HH:MM:SS
 config['dt'] = 10  # time increment [s]]
 config['time'] = 0.0  # initial value of time [s]
 config['time_out'] = 86400  # time between outputs [s]
diff --git a/Python_files/Code/build_config_files/build_config_MOSAiC.py b/Python_files/Code/build_config_files/build_config_MOSAiC.py
deleted file mode 100644
index 0aab77e..0000000
--- a/Python_files/Code/build_config_files/build_config_MOSAiC.py
+++ /dev/null
@@ -1,130 +0,0 @@
-import json
-from func_bound_values import *
-import os
-
-# set wd to build_config_files folder
-wd = os.getcwd()
-# check if already in correct wd
-if wd[-18:] == 'build_config_files':
-    pass
-# if not, change to build_config_files
-else:
-    path = wd + '/Python_files/Code/build_config_files'
-    os.chdir(path)
-
-# set constants
-rho_l = 1028
-c_l = 3400
-
-# **********************************************************************************************************************
-# Initialize array containig all parameters and set path to its directory
-# **********************************************************************************************************************
-config = {}
-path_config = r'../../../Run_specifics/'
-path_input = '../../../input/'
-
-# **********************************************************************************************************************
-# Set description of run
-# **********************************************************************************************************************
-description = "Testcase_MOSAiC"
-file = open(path_config + r"description.txt", "w")
-file.write(description)
-file.close()
-
-# **********************************************************************************************************************
-# Initial layer thickness, timestep, output time and simulation time
-# **********************************************************************************************************************
-config['dt'] = 20.  # time increment [s]]
-config['time'] = 0.0  # initial value of time [s]
-config['time_out'] = 86400.  # time between outputs [s]
-config['time_total'] = 55468800.  # total length of simulation [s]
-
-# **********************************************************************************************************************
-# Time settings needed when input is given
-# **********************************************************************************************************************
-config['timestep_data'] = 60  # timestep of input data [s]
-config['length_input'] = config['time_total'] / 60. + 1  # Length of your input files. Must match with timestep_data
-
-# **********************************************************************************************************************
-# Layer settings and allocation
-# **********************************************************************************************************************
-config['thick_0'] = 0.02
-config['Nlayer'] = 80
-config['N_active'] = 1
-config['N_top'] = 20
-config['N_bottom'] = 20
-config['N_middle'] = config['Nlayer'] - config['N_top'] - config['N_bottom']
-
-# **********************************************************************************************************************
-# Flags
-# **********************************************************************************************************************
-# ________________________top heat flux____________
-config['boundflux_flag'] = 2
-config['albedo_flag'] = 2
-# ________________________brine_dynamics____________
-config['grav_heat_flag'] = 1
-config['flush_heat_flag'] = 1
-config['flood_flag'] = 2
-config['flush_flag'] = 5
-config['grav_flag'] = 2
-config['harmonic_flag'] = 2
-# ________________________Salinity____________
-config['prescribe_flag'] = 1
-config['salt_flag'] = 1
-# ________________________bottom setting______________________
-config['turb_flag'] = 2  # was on two in 203
-config['bottom_flag'] = 1
-config['tank_flag'] = 1
-# ________________________snow______________________
-config['precip_flag'] = 1
-config['freeboard_snow_flag'] = 0  # < Niels, 2017
-config['snow_flush_flag'] = 1  # < Niels, 2017
-config['styropor_flag'] = 0
-config['lab_snow_flag'] = 0
-# ________________________debugging_____________________
-config['debug_flag'] = 1  # set to 2 for output of all ice layers each timestep
-# ________________________bgc_______________________
-config['bgc_flag'] = 1
-# ________________________initial state_______________________
-config['initial_state_flag'] = 1  # 2 if initial state is given
-
-# **********************************************************************************************************************
-# Tank and turbulent fluxes settings
-# **********************************************************************************************************************
-config['tank_depth'] = 0
-config['alpha_flux_stable'] = 0
-config['alpha_flux_instable'] = 0
-
-# **********************************************************************************************************************
-# BGC Settings
-# **********************************************************************************************************************
-config['N_bgc'] = 2
-config['bgc_bottom_1'] = 400
-config['bgc_bottom_2'] = 500
-
-# **********************************************************************************************************************
-# Construct Input files
-# **********************************************************************************************************************
-
-# Set constant inputs
-const_inputs = {'T_bottom': -1.8, 'S_bu_bottom': 34,
-                'fl_q_bottom': 1., 'precip_s':0.}
-for input in list(const_inputs.keys()):
-    data = np.ones(int(config['length_input'])) * const_inputs[input]
-    np.savetxt(path_input + input + '.txt', data)
-
-# **********************************************************************************************************************
-# setting the initial values of the top and only layer - not needed if initial ice properties are given - set to zero
-# **********************************************************************************************************************
-config['thick_1'] = config['thick_0']
-config['m_1'] = config['thick_0'] * rho_l
-config['S_abs_1'] = config['m_1'] * const_inputs['S_bu_bottom']
-config['H_abs_1'] = config['m_1'] * const_inputs['T_bottom'] * c_l
-
-# **********************************************************************************************************************
-# Write init to .json file
-# **********************************************************************************************************************
-json_object = json.dumps(config, indent=4)
-
-with open(path_config + "config.json", "w") as outfile:
-    outfile.write(json_object)
diff --git a/Python_files/Code/build_config_files/build_config_initial_state.py b/Python_files/Code/build_config_files/build_config_initial_state.py
index 9962425..967f429 100644
--- a/Python_files/Code/build_config_files/build_config_initial_state.py
+++ b/Python_files/Code/build_config_files/build_config_initial_state.py
@@ -2,15 +2,10 @@ import json
 from func_bound_values import *
 import os
 
-# set wd to build_config_files folder
-wd = os.getcwd()
-# check if already in correct wd
-if wd[-18:] == 'build_config_files':
-    pass
-# if not, change to build_config_files
-else:
-    path = wd + '/Python_files/Code/build_config_files'
-    os.chdir(path)
+# set wd to the directory of the script
+abspath = os.path.abspath(__file__)
+dname = os.path.dirname(abspath)
+os.chdir(dname)
 
 # set constants
 rho_l = 1028
@@ -34,6 +29,7 @@ file.close()
 # **********************************************************************************************************************
 # Initial layer thickness, timestep, output time and simulation time
 # **********************************************************************************************************************
+config['start_time'] = '2022-01-01 00:00:00'  # start time in YYYY-mm-dd HH:MM:SS
 config['dt'] = 1  # time increment [s]]
 config['time'] = 0.0  # initial value of time [s]
 config['time_out'] = 3600  # time between outputs [s]
diff --git a/Python_files/Code/build_config_files/build_config_mosaic1.py b/Python_files/Code/build_config_files/build_config_mosaic1.py
index 8a4372d..68169dc 100644
--- a/Python_files/Code/build_config_files/build_config_mosaic1.py
+++ b/Python_files/Code/build_config_files/build_config_mosaic1.py
@@ -2,15 +2,10 @@ import json
 from func_bound_values import *
 import os
 
-# set wd to build_config_files folder
-wd = os.getcwd()
-# check if already in correct wd
-if wd[-18:] == 'build_config_files':
-    pass
-# if not, change to build_config_files
-else:
-    path = wd + '/Python_files/Code/build_config_files'
-    os.chdir(path)
+# set wd to the directory of the script
+abspath = os.path.abspath(__file__)
+dname = os.path.dirname(abspath)
+os.chdir(dname)
 
 
 # set constants
@@ -35,6 +30,7 @@ file.close()
 # **********************************************************************************************************************
 # Initial layer thickness, timestep, output time and simulation time
 # **********************************************************************************************************************
+config['start_time'] = '2019-01-01 00:00:00'  # start time in YYYY-mm-dd HH:MM:SS
 config['dt'] = 20  # time increment [s]]
 config['time'] = 0.0  # initial value of time [s]
 config['time_out'] = 86400  # time between outputs [s]
@@ -108,7 +104,7 @@ config['bgc_bottom_2'] = 0
 # **********************************************************************************************************************
 # Set constant inputs
 const_inputs = {'T_bottom': -1, 'S_bu_bottom': 34,
-                'fl_q_bottom': 10, 'precip_s': 0}
+                'fl_q_bottom': 1, 'precip_s': 0}
 for input in list(const_inputs.keys()):
     data = np.ones(config['length_input']) * const_inputs[input]
     np.savetxt(path_input + input + '.txt', data)
@@ -120,7 +116,7 @@ for input in list(const_inputs.keys()):
 config['thick_1'] = config['thick_0']
 config['m_1'] = config['thick_0'] * rho_l
 config['S_abs_1'] = config['m_1'] * const_inputs['S_bu_bottom']
-config['H_abs_1'] = 0
+config['H_abs_1'] = config['m_1'] * const_inputs['T_bottom'] * c_l
 
 # **********************************************************************************************************************
 # Write init to .json file
diff --git a/Python_files/Code/plotscripts/plot_TPhiS.py b/Python_files/Code/plotscripts/plot_TPhiS.py
index 98b3901..3b495cc 100644
--- a/Python_files/Code/plotscripts/plot_TPhiS.py
+++ b/Python_files/Code/plotscripts/plot_TPhiS.py
@@ -11,33 +11,29 @@
 #Loading modules and setting fonts
 import numpy
 import matplotlib.pyplot as plt
+from matplotlib.dates import DateFormatter
 from matplotlib import rc
 import matplotlib
-#rc('text', usetex=True)
-#rc('font', size='10')
-#rc('font', family='serif')
-#Warnings:
-#Contours are interpolated from the middle of each layer. This is most visible in thick layers, and in the snow layer where the contour lines only extend to the the middle. 
-
-#Settings 
-dx           = 0.5 
-timeunit     = '[days]'
-outputfile   = 'pic_TPhiS'
-outputformat = 'png' #e.g. png, jpg, pdf
-free_flag    = 1     #1: freeboard is included, 0:freeboard is not included 
-
-#Contour levels 
-levelsT      = ([-10,-5,-3,-1])
-levelspsi    = ([0.1, 0.2])
-levelsS      = ([3., 8.])
+import os
+import json
+import pandas as pd
 
+# Warnings:
+# Contours are interpolated from the middle of each layer. This is most visible in thick layers, and in the snow layer where the contour lines only extend to the the middle.
 
+# Settings
+outputfile   = 'pic_TPhiS'
+outputformat = 'png' #e.g. png, jpg, pdf
+free_flag    = 1     #1: freeboard is included, 0:freeboard is not included
 
+# set wd to the directory of the script
+abspath = os.path.abspath(__file__)
+dname = os.path.dirname(abspath)
+os.chdir(dname)
 
-#
 
 
-#Loading data
+# Loading data
 S          = numpy.loadtxt("../../../output/dat_S_bu.dat")
 T          = numpy.loadtxt("../../../output/dat_T.dat")
 psi_l      = numpy.loadtxt("../../../output/dat_psi_l.dat")
@@ -45,6 +41,23 @@ thick      = numpy.loadtxt("../../../output/dat_thick.dat")
 snow       = numpy.loadtxt("../../../output/dat_snow.dat")
 freeboard  = numpy.loadtxt("../../../output/dat_freeboard.dat")
 
+# Load config file
+with open('../../../Run_specifics/config.json') as json_file:
+    config = json.load(json_file)
+
+# build time axis
+offset = pd.DateOffset(seconds=config['time_out'])
+time = pd.date_range(config['start_time'], freq=offset, periods=config['time_total']/config['time_out'] + 1).to_series()
+dx           = config['time_out']/(60*60*24)  # get dx in days
+timeunit     = '[days]'
+
+#Contour levels 
+levelsT      = ([-10,-5,-3,-1])
+levelspsi    = ([0.1, 0.2])
+levelsS      = ([3., 8.])
+
+
+
 #Setting freeboard to zero if free_flag = 0
 if   free_flag == 0:
   freeboard[:] = 0.
@@ -68,9 +81,6 @@ T     = numpy.hstack((T_snow,T))
 psi_l = numpy.hstack((psi_l_snow,psi_l))
 S     = numpy.hstack((S_snow,S))
 
-ylen = len(thick[0,:])
-xlen = len(thick[:,0])
-
 #Restructuring the data so it can be ploted by pcolor
 depth = thick*1.
 depth_contour = thick*1.
@@ -83,17 +93,21 @@ i=0
 j=0
 ireal = 0.
 while (i<xlen):
-  while (j<ylen):
-    depth[i,j]=-sum(thick[i,0:j])+freeboard[i]+thick_snow[i]
-    #Contour depth is slightly different
-    depth_contour[i,j]=-sum(thick[i,0:j])-thick[i,j]/2.+freeboard[i]+thick_snow[i]
-    Xgrid[i,j]=ireal
-    j=j+1
-  i=i+1
-  j=0
-  ireal=ireal+dx
+    while (j<ylen):
+        depth[i,j]=-sum(thick[i,0:j])+freeboard[i]+thick_snow[i]
+        #Contour depth is slightly different
+        depth_contour[i,j]=-sum(thick[i,0:j])-thick[i,j]/2.+freeboard[i]+thick_snow[i]
+        Xgrid[i,j]=ireal
+        j=j+1
+    i=i+1
+    j=0
+    ireal=ireal+dx
 
+depth = numpy.column_stack((depth,depth[:,-1]-(depth[:,-1]-depth[:,-2])))
+Xgrid = numpy.column_stack((Xgrid,Xgrid[:,-1]))
 
+depth = numpy.vstack((depth, depth[-1,:]))
+Xgrid = numpy.vstack((Xgrid, Xgrid[-1,:]))
 
 
 #Custom colormaps
@@ -202,7 +216,7 @@ c1 = plt.colorbar(pad=0.01)
 c1.set_label(r'T')
 plt.axis([Xgrid.min(), Xgrid.max(), ymin, ymax])
 #ax1.fill_between(Xaxis[:],freeboard[:], snow[:,0]+freeboard[:,],facecolor='white',edgecolor='white')
-CS1 = plt.contour(Xgrid, depth_contour, T, levelsT,colors='k')
+CS1 = plt.contour(Xgrid[:-1,:-1], depth_contour, T, levelsT,colors='k')
 plt.clabel(CS1, fontsize=9, inline=1,fmt='%1.0f')
 plt.plot(Xaxis[:],freeboard[:],'k--')
 plt.ylabel(r'depth [m]')
@@ -223,7 +237,7 @@ plt.pcolor(Xgrid,depth,psi_l,cmap=psi_cmap,vmin=zmin,vmax=zmax)
 plt.axis([Xgrid.min(), Xgrid.max(), ymin, ymax])
 c2 = plt.colorbar(pad=0.01)
 c2.set_label(r'$\phi_l$')
-CS2 = plt.contour(Xgrid, depth_contour, psi_l, levelspsi, colors='k')
+CS2 = plt.contour(Xgrid[:-1,:-1], depth_contour, psi_l, levelspsi, colors='k')
 plt.clabel(CS2, fontsize=9, inline=1,fmt='%1.1f')
 #ax2.fill_between(Xaxis[:],freeboard[:], snow[:,0]+freeboard[:,],facecolor='white',edgecolor='white')
 plt.plot(Xaxis[:],freeboard[:],'k--')
@@ -244,16 +258,23 @@ plt.pcolor(Xgrid,depth,S,cmap=S_cmap,vmin=zmin,vmax=zmax)
 plt.axis([Xgrid.min(), Xgrid.max(), ymin, ymax])
 c3 = plt.colorbar(pad=0.01)
 c3.set_label(r'$S_{bu}$')
-CS3 = plt.contour(Xgrid, depth_contour, S, levelsS,colors='k')
+CS3 = plt.contour(Xgrid[:-1,:-1], depth_contour, S, levelsS,colors='k')
 plt.clabel(CS3, fontsize=9, inline=1,fmt='%2.0f')
 #ax3.fill_between(Xaxis[:],freeboard[:], snow[:,0]+freeboard[:,],facecolor='white',edgecolor='white')
 plt.plot(Xaxis[:],freeboard[:],'k--')
 plt.ylabel(r'depth [m]')
-plt.xlabel(r'time '+timeunit)
 ax3.set_facecolor((0.5,0.5,0.5))
 
+# Set x ticks with times
+ticks = ax3.get_xticks() * (60*60*24)
+index = (ticks / config['time_out']).astype(int)
+times_plot = [str(time[idx].strftime(format='%Y-%m-%d %Hh')) for idx in index]
+ax3.set_xticklabels(times_plot)
+plt.xticks(rotation=20)
+
+plt.tight_layout()
 
 #Saving and exporting
-plt.savefig('../../Plots/'+outputfile+'.'+outputformat,dpi=1000)
+plt.savefig('../../Plots/'+outputfile+'.'+outputformat, dpi=1000)
 plt.close()
 
diff --git a/Python_files/Plots/.gitkeep b/Python_files/Plots/.gitkeep
new file mode 100644
index 0000000..e69de29
diff --git a/Run_specifics/config.json b/Run_specifics/config.json
index 8bdcc22..a2c92d2 100644
--- a/Run_specifics/config.json
+++ b/Run_specifics/config.json
@@ -1,36 +1,37 @@
 {
-    "dt": 20.0,
+    "start_time": "2020-04-19 12:00:00",
+    "dt": 1,
     "time": 0.0,
-    "time_out": 86400.0,
-    "time_total": 55468800.0,
-    "timestep_data": 60,
-    "length_input": 924481.0,
-    "thick_0": 0.02,
-    "Nlayer": 80,
+    "time_out": 3600,
+    "time_total": 259200,
+    "length_input": 259200.0,
+    "timestep_data": 1,
+    "thick_0": 0.002,
+    "Nlayer": 90,
     "N_active": 1,
-    "N_top": 20,
-    "N_bottom": 20,
-    "N_middle": 40,
-    "boundflux_flag": 2,
+    "N_top": 5,
+    "N_bottom": 5,
+    "N_middle": 80,
+    "boundflux_flag": 1,
     "albedo_flag": 2,
     "grav_heat_flag": 1,
     "flush_heat_flag": 1,
     "flood_flag": 2,
-    "flush_flag": 5,
+    "flush_flag": 1,
     "grav_flag": 2,
     "harmonic_flag": 2,
     "prescribe_flag": 1,
-    "salt_flag": 1,
-    "turb_flag": 2,
+    "salt_flag": 2,
+    "turb_flag": 1,
     "bottom_flag": 1,
     "tank_flag": 1,
-    "precip_flag": 1,
+    "precip_flag": 0,
     "freeboard_snow_flag": 0,
     "snow_flush_flag": 1,
     "styropor_flag": 0,
     "lab_snow_flag": 0,
     "debug_flag": 1,
-    "bgc_flag": 1,
+    "bgc_flag": 2,
     "initial_state_flag": 1,
     "tank_depth": 0,
     "alpha_flux_stable": 0,
@@ -38,8 +39,8 @@
     "N_bgc": 2,
     "bgc_bottom_1": 400,
     "bgc_bottom_2": 500,
-    "thick_1": 0.02,
-    "m_1": 20.56,
-    "S_abs_1": 699.04,
-    "H_abs_1": -125827.19999999998
+    "thick_1": 0.002,
+    "m_1": 2.056,
+    "S_abs_1": 69.904,
+    "H_abs_1": -6990.400000000001
 }
\ No newline at end of file
diff --git a/Run_specifics/description.txt b/Run_specifics/description.txt
index fd739a7..7badf45 100644
--- a/Run_specifics/description.txt
+++ b/Run_specifics/description.txt
@@ -1 +1 @@
-Testcase_MOSAiC
\ No newline at end of file
+Testcase_1
\ No newline at end of file
diff --git a/mo_data.f90 b/mo_data.f90
index 280fe10..fbb19f2 100644
--- a/mo_data.f90
+++ b/mo_data.f90
@@ -214,5 +214,6 @@ MODULE mo_data
   INTEGER                               :: length_input_lab     !< Niels, 2017 add: used to allocate lab testcase input arrays in mo_init, set value in testcases
   REAL(wp)                              :: grav_heat_flux_down = 0._wp  !< output variable for heat loss at the ice bottom caused by gravity drainage (negative = gain)
   REAL(wp)                              :: grav_heat_flux_up = 0._wp !< output variable for heat gain at the ice bottom caused by gravity drainage (sea water replacing brine, positive = loss)
+  REAL(wp)                              :: exp_heat = 0._wp
 END MODULE mo_data
 
diff --git a/mo_grav_drain.f90 b/mo_grav_drain.f90
index b68f0a9..26da096 100644
--- a/mo_grav_drain.f90
+++ b/mo_grav_drain.f90
@@ -171,7 +171,7 @@ CONTAINS
        END IF
     END DO
 
-    grav_salt  = grav_salt - SUM(S_abs(:))
+    
 
     !Influence of summed up upward transport resulting from the downward drainage
     fl_m(1)=0.0_wp
@@ -188,6 +188,7 @@ CONTAINS
 
     CALL mass_transfer (Nlayer,N_active,T,H_abs,S_abs,S_bu,T_bottom,S_bu_bottom,fl_m)
 
+    grav_salt  = grav_salt - SUM(S_abs(:))
     grav_drain = grav_drain+fl_m(N_active+1)
 
     !Heatflux correction
diff --git a/mo_grotz.f90 b/mo_grotz.f90
index 1f26ba1..008deed 100644
--- a/mo_grotz.f90
+++ b/mo_grotz.f90
@@ -269,7 +269,7 @@ CONTAINS
                 melt_thick_snow = 0.0_wp  !< Niels, 2017 add: melt_thick_snow
             END IF
 
-
+            exp_heat=exp_heat+SUM(H_abs(:))
             !##########################################################################################
             !Inner layer thermodynamics and Expulsion
             !##########################################################################################
@@ -297,7 +297,7 @@ CONTAINS
                     END DO
                 END IF
             END IF
-
+            exp_heat=exp_heat-SUM(H_abs(:))
 
             !##########################################################################################
             !Raw Output, if debug flag is set to 2 then all layer data is output
@@ -336,13 +336,14 @@ CONTAINS
 
                 grav_heat_flux_down = grav_heat_flux_down / time_out
                 grav_heat_flux_up = grav_heat_flux_up / time_out
+                exp_heat = exp_heat / time
 
                 !Calling standard output
                 CALL output(Nlayer, T, psi_s, psi_l, thick, S_bu, ray, format_T, format_psi, &
                         format_thick, format_snow, freeboard, thick_snow, T_snow, psi_l_snow, psi_s_snow, &
                         energy_stored, freshwater, total_resist, thickness, bulk_salin, &
                         grav_drain, grav_salt, grav_temp, T2m, T_top, perm, format_perm, flush_v, flush_h, psi_g, &
-                        & melt_thick_output, grav_heat_flux_down, grav_heat_flux_up, format_melt)
+                        & melt_thick_output, grav_heat_flux_down, grav_heat_flux_up, exp_heat, format_melt)
 
                 IF (bgc_flag==2) THEN
                     CALL output_bgc(Nlayer, N_active, bgc_bottom, N_bgc, bgc_abs, psi_l, thick, m, format_bgc)
@@ -375,8 +376,8 @@ CONTAINS
                         'T_top: ', T_top, &
                         'fl_sw: ', fl_sw, &
                         'fl_rest; ', fl_rest, &
-						', Habs_top:', H_abs(2), &
-						', Habs_bot:', H_abs(N_active)
+                        ', Habs_top:', H_abs(2), &
+                        ', Habs_bot:', H_abs(N_active)
 
 
 
@@ -385,6 +386,7 @@ CONTAINS
                 grav_drain = 0.0_wp
                 grav_salt = 0.0_wp
                 grav_temp = 0.0_wp
+                exp_heat =0.0_wp
                 grav_heat_flux_down = 0._wp
                 grav_heat_flux_up = 0._wp
                 melt_thick_output(:) = 0._wp
diff --git a/mo_output.f90 b/mo_output.f90
index 92cf3d8..6393e62 100644
--- a/mo_output.f90
+++ b/mo_output.f90
@@ -1,390 +1,390 @@
-!>
-!! All things output.
-!!
-!! Used to clean up root.f90 and make it easier to implement changes to the output. 
-!!
-!!
-!! @author Philipp Griewank
-!!
-!!  COPYRIGHT
-!!
-!! This file is part of SAMSIM.
-!!
-!!  SAMSIM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
-!!
-!!  SAMSIM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
-!!
-!!  You should have received a copy of the GNU General Public License along with SAMSIM. If not, see <http://www.gnu.org/licenses/>.
-!
-!!
-!! @par Revision History
-!! Brought from the womb by Philipp Griewank, IMPRS (2010-10-11) \n
-!! add some output values by Niels Fuchs, MPIMET (2017-03-01)
-!!
-MODULE mo_output
-
-  USE mo_parameters
-  IMPLICIT NONE
-
-  PRIVATE
-  PUBLIC :: output_begin,output_begin_bgc,output_raw,output_raw_snow,output_raw_lay,output,output_bgc,output_settings!,output_end
-
-CONTAINS
-  !>
-  !! Settings output.
-  !!
-  !! Writes important values to latter identify run.
-  !!
-  !! @par Revision History
-  !! created by Philipp Griewank, IMPRS (2011-02-12)
-  !!
-  SUBROUTINE output_settings(description,N_top,N_bottom,Nlayer,fl_q_bottom,T_bottom,S_bu_bottom,thick_0,time_out,    &
-    time_total,dt,boundflux_flag,albedo_flag,grav_flag,flush_flag,flood_flag,grav_heat_flag,flush_heat_flag,    &
-    harmonic_flag,prescribe_flag,salt_flag,turb_flag,bottom_flag,tank_flag,precip_flag,bgc_flag,N_bgc,k_snow_flush)
-    INTEGER,         INTENT(in) :: N_top,N_bottom,Nlayer, boundflux_flag,albedo_flag,grav_flag,flush_flag, &
-                                   flood_flag,grav_heat_flag , flush_heat_flag,harmonic_flag,prescribe_flag,salt_flag,turb_flag,  &
-                                   bottom_flag,tank_flag,precip_flag,bgc_flag,N_bgc
-    REAL(wp),        INTENT(in) :: fl_q_bottom,T_bottom,S_bu_bottom,thick_0,time_out,time_total,dt,k_snow_flush
-    CHARACTER*12000, INTENT(in) :: description
-    
-    
-    OPEN(1234,file='./output/dat_settings.dat' ,STATUS='replace')
-    WRITE(1234,*)'################  Description  ###############'
-    WRITE(1234,*) TRIM(description)
-
-       
-    !*************************************************************************************************************************
-    !Layer settings and allocation
-    !*************************************************************************************************************************
-    WRITE(1234,*) '##############  Basic settings  ##############'
-    WRITE(1234,'(A16,F15.3)') 'dt              =',dt
-    WRITE(1234,'(A16,F15.3)') 'thick_0         =',thick_0
-    WRITE(1234,'(A16,F15.3)') 'time_out        =',time_out
-    WRITE(1234,'(A16,F15.3)') 'time_total      =',time_total
-    
-    
-    WRITE(1234,'(A16,F15.3)') 'fl_q_bottom     =',fl_q_bottom
-    WRITE(1234,'(A16,F15.3)') 'T_bottom        =',T_bottom
-    WRITE(1234,'(A16,F15.3)') 'S_bu_bottom     =',S_bu_bottom
-    
-    
-    WRITE(1234,'(A16,I9.0)')  'N_top           =',N_top
-    WRITE(1234,'(A16,I9.0)')  'N_middle        =',Nlayer-N_top-N_bottom
-    WRITE(1234,'(A16,I9.0)')  'N_bottom        =',N_bottom
-    WRITE(1234,'(A16,I9.0)')  'Nlayer          =',Nlayer
-    
-    !*************************************************************************************************************************
-    !Flags
-    !*************************************************************************************************************************
-    WRITE(1234,*) '##################  Flags  ###################'
-    !________________________top heat flux____________
-    WRITE(1234,'(A16,I9.0)')  'boundflux_flag  =',boundflux_flag
-    WRITE(1234,'(A16,I9.0)')  'albedo_flag     =',albedo_flag
-    !________________________brine_dynamics____________
-    WRITE(1234,'(A16,I9.0)')  'grav_flag       =',grav_flag
-    WRITE(1234,'(A16,I9.0)')  'flush_flag      =',flush_flag
-    WRITE(1234,'(A16,I9.0)')  'flood_flag      =',flood_flag
-    WRITE(1234,'(A16,I9.0)')  'grav_heat_flag  =',grav_heat_flag
-    WRITE(1234,'(A16,I9.0)')  'flush_heat_flag =',flush_heat_flag
-    WRITE(1234,'(A16,I9.0)')  'harmonic_flag   =',harmonic_flag
-    WRITE(1234,'(A16,F15.3)')  'k_snow_flush    =',k_snow_flush !< Niels, 2017
-    !________________________Salinity____________
-    WRITE(1234,'(A16,I9.0)')  'prescribe_flag  =',prescribe_flag
-    WRITE(1234,'(A16,I9.0)')  'salt_flag       =',salt_flag
-    !________________________bottom setting______________________
-    WRITE(1234,'(A16,I9.0)')  'turb_flag       =',turb_flag
-    WRITE(1234,'(A16,I9.0)')  'bottom_flag     =',bottom_flag
-    WRITE(1234,'(A16,I9.0)')  'tank_flag       =',tank_flag
-    WRITE(1234,'(A16,I9.0)')  'precip_flag     =',precip_flag
-    WRITE(1234,'(A16,I9.0)')  'bgc_flag        =',bgc_flag
-    WRITE(1234,'(A16,I9.0)')  'N_bgc           =',N_bgc
-
-
-    CLOSE(1234)
-  END SUBROUTINE output_settings
-
-  !>
-  !! Standard output.
-  !!
-  !! For time=n*time_out data is exported.
-  !!
-  !! @par Revision History
-  !! created by Philipp Griewank, IMPRS (2010-10-11)
-  !!
-  SUBROUTINE output(Nlayer,T,psi_s,psi_l,thick,S_bu,ray,format_T,format_psi, &
-       format_thick,format_snow,freeboard,thick_snow,T_snow,psi_l_snow,psi_s_snow,           &
-       energy_stored,freshwater,total_resist,thickness,bulk_salin,                &
-       grav_drain,grav_salt,grav_temp,T2m,T_top,perm,format_perm,flush_v,flush_h,psi_g,melt_thick_output, &
-       grav_heat_flux_down, grav_heat_flux_up,format_melt)
-    INTEGER,                       INTENT(in) :: Nlayer
-    REAL(wp), DIMENSION(Nlayer),   INTENT(in) :: T,psi_s,psi_l,thick,S_bu,perm,flush_v,flush_h,psi_g
-    REAL(wp), DIMENSION(Nlayer-1), INTENT(in) :: ray
-    REAL(wp),                      INTENT(in) :: freeboard,thick_snow,T_snow,psi_l_snow,psi_s_snow,energy_stored,&
-                                                 freshwater,thickness,bulk_salin, &
-                                                 total_resist,grav_drain,grav_salt,grav_temp,T2m,T_top, &
-                                                 grav_heat_flux_down, grav_heat_flux_up
-    REAL(wp), DIMENSION(3),        INTENT(in) :: melt_thick_output  !< Niels, 2017: 1: accumulated melt_thick, 2: accumulated melt_thick_snow, 3: accumulated top ice thickness variations (recheck 3: in mo_layer_dynamics)
-    CHARACTER*12000,               INTENT(in) :: format_T,format_psi,format_thick,format_snow,format_perm,format_melt
-
-    WRITE(30,format_T)     T
-    WRITE(31,format_psi)   psi_s
-    WRITE(32,format_thick) thick
-    WRITE(33,format_T)     S_bu
-    WRITE(34,format_T)     ray
-    WRITE(35,format_psi)   psi_l
-    WRITE(40,'(F9.3)')     freeboard
-    WRITE(41,format_snow)  thick_snow,T_snow,psi_l_snow,psi_s_snow
-    WRITE(42,'(F15.1,2x,F10.5,2x,F10.5,2x,F10.5,2x,F10.5)')energy_stored,freshwater,total_resist,thickness,bulk_salin
-    WRITE(43,'(F9.6,2X,F9.5,2X,F7.3)')grav_drain,grav_salt,grav_temp
-    WRITE(45,*)            T2m,T_top
-    WRITE(46,format_perm) perm  !< Niels, 2017 add: output permeability
-    WRITE(47,format_perm) flush_v   !< Niels, 2017 add: output vertical flushing
-    WRITE(48,format_perm) flush_h   !< Niels, 2017 add: output horizontal flushing
-    WRITE(49,format_psi) psi_g  !< Niels, 2017 add: output gas fraction !OBS: not simulated physically in SAMSIM
-    WRITE(50,format_melt) melt_thick_output !< Niels, 2017
-    WRITE(51,format_thick) grav_heat_flux_down !< Niels, 2017
-    WRITE(52,format_thick) grav_heat_flux_up !< Niels, 2017
-  
-  END SUBROUTINE output
-
-  !>
-  !! Standard bgc output.
-  !!
-  !! For time=n*time_out data is exported.
-  !!
-  !! @par Revision History
-  !! created by Philipp Griewank, IMPRS (2014-02-06)
-  !!
-  SUBROUTINE output_bgc(Nlayer,N_active,bgc_bottom,N_bgc,bgc_abs,psi_l,thick,m,format_bgc)
-    INTEGER,                             INTENT(in) :: Nlayer, N_bgc, N_active
-    REAL(wp), DIMENSION(N_bgc),          INTENT(in) :: bgc_bottom
-    REAL(wp), DIMENSION(Nlayer),         INTENT(in) :: psi_l,m,thick
-    REAL(wp), DIMENSION(Nlayer,N_bgc),   INTENT(in) :: bgc_abs
-    REAL(wp), DIMENSION(Nlayer)                     :: bgc_bu,bgc_br
-    CHARACTER*12000,                     INTENT(in) :: format_bgc
-    INTEGER :: k,kk
-    
-    DO k=1,N_bgc
-      DO kk=1,Nlayer
-        IF (kk<=N_active) THEN
-           IF (abs(m(kk)) > tolerance) THEN
-              bgc_bu(kk) = bgc_abs(kk,k)/m(kk)
-              IF (abs(psi_l(kk)) > tolerance .AND. abs(thick(kk)) > tolerance) THEN
-                 bgc_br(kk) = bgc_abs(kk,k)/psi_l(kk)/thick(kk)/rho_l
-              ELSE
-                 bgc_br(kk) = 0000._wp
-              END IF
-           ELSE
-              bgc_bu(kk) = 0._wp
-              bgc_br(kk) = 0._wp
-           END IF
-        ELSE
-              bgc_bu(kk) = bgc_bottom(k)
-              bgc_br(kk) = bgc_bottom(k)
-        END IF
-      END DO
-      WRITE(2*k+400,format_bgc)bgc_bu
-      WRITE(2*k+401,format_bgc)bgc_br
-    END DO
-  
-  END SUBROUTINE output_bgc
-
-  !>
-  !! Output for debugging purposes.
-  !!
-  !! Data for each layer is written out each time step to aid in finding errors or understanding model behavior.
-  !!
-  !! @par Revision History
-  !! created by Philipp Griewank, IMPRS (2010-10-11)
-  !!
-  SUBROUTINE output_raw(Nlayer,N_active,time,T,thick,S_bu,psi_s,psi_l,psi_g)
-    INTEGER,                     INTENT(in) :: Nlayer,N_active
-    INTEGER                                 :: k
-    REAL(wp),                    INTENT(in) :: time
-    REAL(wp), DIMENSION(Nlayer), INTENT(in) :: T,thick,S_bu,psi_s,psi_l,psi_g
-
-    DO k = 1,Nlayer
-       IF ( k<=N_active ) THEN
-          WRITE(k+200, '(F8.4,2x,f10.3,2x,f7.5,2x,f9.5,2x,f4.2,2x,f4.2,2x,f4.2)') &
-               time/(86400._wp),T(k),thick(k),S_bu(k),psi_s(k),psi_l(k),psi_g(k)
-       ELSE
-          WRITE(k+200,'(F6.2,2x,f7.3,2x,f5.3,2x,f8.5,2x,f4.2,2x,f4.2,2x,f4.2)')time/(86400._wp),0.0,0.0,0.0,0.0,0.0,0.0
-       END IF
-    END DO
-
-  END SUBROUTINE output_raw
-
-  !>
-  !! Output for debugging purposes.
-  !!
-  !! Data of snow layer is written out at each time step to aid in finding errors or understanding model behavior.
-  !!
-  !! @par Revision History
-  !! created by Philipp Griewank, IMPRS (2010-10-11)
-  !!
-  SUBROUTINE output_raw_snow(time,T_snow,thick_snow,S_abs_snow,m_snow,psi_s_snow,psi_l_snow,psi_g_snow)
-    REAL(wp), INTENT(in) :: time
-    REAL(wp), INTENT(in) :: T_snow,thick_snow,S_abs_snow,m_snow,psi_s_snow,psi_l_snow,psi_g_snow
-
-    IF (thick_snow>0.0_wp) THEN
-       WRITE(66, '(F8.4,2x,f10.3,2x,f5.3,2x,f4.1,2x,f4.2,2x,f4.2,2x,f4.2)') &
-            time/(30._wp*86400._wp),T_snow,thick_snow,S_abs_snow/MAX(m_snow,0.001_wp),psi_s_snow,psi_l_snow,psi_g_snow
-    ELSE 
-       WRITE(66, '(F8.4,2x,f10.3,2x,f5.3,2x,f4.1,2x,f4.2,2x,f4.2,2x,f4.2)') &
-            time/(30._wp*86400._wp),0.0,0.0,0.0,0.0,0.0,0.0
-    END IF
-
-  END SUBROUTINE output_raw_snow
-
-  !>
-  !! Output for debugging layer dynamics..
-  !!
-  !! Is used when debug_flag = 2 to track when which layer dynamics occur (see mo_layer_dynamics).
-  !!
-  !!
-
-  SUBROUTINE output_raw_lay(Nlayer,N_active,H_abs,m,S_abs,thick,string)
-    INTEGER,                     INTENT(in) :: Nlayer,N_active
-    INTEGER                                 :: k
-    REAL(wp), DIMENSION(Nlayer), INTENT(in) :: H_abs,S_abs,thick,m 
-    REAL(wp)                                :: mm 
-    CHARACTER*6,                 INTENT(in) :: string
-
-    DO k = 1,Nlayer
-       IF ( k<=N_active ) THEN
-          IF (abs(m(k)) < tolerance) THEN
-             mm = 99999999._wp
-          ELSE
-             mm = m(k)
-          END IF
-          WRITE(k+200, '(A6,2x,f11.0,2x,f6.3,2x,f8.5,2x,f9.3,2x,I2)') &
-               string,h_abs(k),thick(k),s_abs(k)/mm,mm/MAX(thick(k),0.0000000000000000001_wp),N_active
-       ELSE 
-          WRITE(k+200, '(A6,2x,f11.0,2x,f6.3,2x,f8.1,2x,f9.3,2x,I2)')string,0.0,0.0,0.0,0.0,N_active
-       END IF
-    END DO
-  END SUBROUTINE output_raw_lay
-
-  !>
-  !! Output files are opened and format strings are created.
-  !!
-  !! Format strings are defined according to the number of layers used which define the output format.
-  !! Files are opened.
-  !!
-  !! @par Revision History
-  !! created by Philipp Griewank, IMPRS (2010-10-11)
-  !! moved by Philipp Griewank, IMPRS (2011-03-09)
-  !!
-  SUBROUTINE output_begin(Nlayer,debug_flag,format_T,format_psi,format_thick,format_snow,format_T2m_top,format_perm,&
-                          &format_melt)
-    INTEGER,         INTENT(in)  :: Nlayer,debug_flag
-    CHARACTER*12000, INTENT(out) :: format_T,format_psi,format_thick,format_snow,format_T2m_top,format_perm,&
-                                    &format_melt
-    INTEGER                      :: k
-    CHARACTER*12000              :: output_string
-
-
-    IF (debug_flag==2) THEN
-       DO k=1,Nlayer
-          IF (k<10) THEN
-             WRITE(output_string,'(A18,I1,A4)')'./output/thermo00',k,'.txt'
-             OPEN(k+200,file=TRIM(ADJUSTL(output_string)),STATUS='replace',Recl=12288)
-          ELSE IF (k<100) THEN
-             WRITE(output_string,'(A17,I2,A4)')'./output/thermo0',k,'.txt'
-             OPEN(k+200,file=TRIM(ADJUSTL(output_string)),STATUS='replace',Recl=12288)
-          ELSE
-             WRITE(output_string,'(A16,I3,A4)')'./output/thermo',k,'.txt'
-             OPEN(k+200,file=TRIM(ADJUSTL(output_string)),STATUS='replace',Recl=12288)
-          END IF
-       END DO
-    END IF
-
-    format_T     = '(F9.3,2x'
-    format_perm  = '(ES14.7,2x'
-    format_psi   = '(F9.3,2x'
-    format_thick = '(F9.5,2x'
-
-    format_snow    = '(F9.3,2x,F9.3,2x,F9.3,2x,F9.3)'
-    format_T2m_top = '(F9.3,2x,F9.3,2x,F9.3)'
-    format_melt    = '(ES14.7,2x,ES14.7,2x,ES14.7)'
-
-    DO k=1,Nlayer
-       format_T     = TRIM(format_T)//',f9.3,2x'
-       format_psi   = TRIM(format_psi)//',f9.3,2x'
-       format_thick = TRIM(format_thick)//',f9.5,2x'
-       format_perm  = TRIM(format_perm)//',ES14.7,2x'
-    END DO
-
-    format_T     = TRIM(format_T)//')'
-    format_psi   = TRIM(format_psi)//')'
-    format_thick = TRIM(format_thick)//')'
-    format_perm  = TRIM(format_perm)//')'
-
-    
-    OPEN(30,file='./output/dat_T.dat'           ,STATUS='replace',Recl=12288)
-    OPEN(31,file='./output/dat_psi_s.dat'       ,STATUS='replace',Recl=12288)
-    OPEN(32,file='./output/dat_thick.dat'       ,STATUS='replace',Recl=12288)
-    OPEN(33,file='./output/dat_S_bu.dat'        ,STATUS='replace',Recl=12288)
-    OPEN(34,file='./output/dat_ray.dat'         ,STATUS='replace',Recl=12288)
-    OPEN(35,file='./output/dat_psi_l.dat'       ,STATUS='replace',Recl=12288)
-    OPEN(40,file='./output/dat_freeboard.dat'   ,STATUS='replace',Recl=12288)
-    OPEN(41,file='./output/dat_snow.dat'        ,STATUS='replace',Recl=12288)
-    OPEN(42,file='./output/dat_vital_signs.dat' ,STATUS='replace',Recl=12288)
-    OPEN(43,file='./output/dat_grav_drain.dat'  ,STATUS='replace',Recl=12288)
-    OPEN(45,file='./output/dat_T2m_T_top.dat'   ,STATUS='replace',Recl=12288)
-    OPEN(46,file='./output/dat_perm.dat'        ,STATUS='replace',Recl=12288)
-    OPEN(47,file='./output/dat_flush_v.dat'        ,STATUS='replace',Recl=12288)
-    OPEN(48,file='./output/dat_flush_h.dat'        ,STATUS='replace',Recl=12288)
-    OPEN(49,file='./output/dat_psi_g.dat'        ,STATUS='replace',Recl=12288)
-    OPEN(50,file='./output/dat_melt.dat'        ,STATUS='replace',Recl=12288)
-    OPEN(51,file='./output/dat_grav_heat_flux_down.dat'        ,STATUS='replace',Recl=12288)
-    OPEN(52,file='./output/dat_grav_heat_flux_up.dat'        ,STATUS='replace',Recl=12288)
-
-    OPEN(66,file='./output/snow.txt'            ,STATUS='replace')
-
-  END SUBROUTINE output_begin
-  
-  
-  !>
-  !! Output files for bgc are opened and format strings are created.
-  !!
-  !! Same thing as out_begin but for bgc
-  !! Each tracer is outputted in bulk and in brine concentration in a separate file.
-  !! Added ADJUSTL to the output strings because they got wierd 
-  !!
-  !! @par Revision History
-  !! created by Dr. Philipp Griewank, MPI (2014-02-07)
-  !! fix by Dr. Philipp Griewank, UniK (2018-05-18)
-  SUBROUTINE output_begin_bgc(Nlayer,N_bgc,format_bgc)
-    INTEGER,         INTENT(in)  :: Nlayer,N_bgc
-    CHARACTER*12000, INTENT(out) :: format_bgc
-    INTEGER                      :: k
-    CHARACTER*12000              :: output_string
-
-     
-
-    format_bgc = '(F16.8,2x'
-
-    DO k=1,Nlayer
-      format_bgc = TRIM(format_bgc)//',f16.8,2x'
-    END DO
-    
-    format_bgc = TRIM(format_bgc)//')'
-
-    DO k=1,N_bgc
-       IF (k<10) THEN
-          WRITE(output_string,'(A18,I1,A7)')'./output/dat_bgc0',k,'.bu.dat'
-          OPEN(2*k+400,file=TRIM(ADJUSTL(output_string)),STATUS='replace',Recl=12288)
-          WRITE(output_string,'(A18,I1,A7)')'./output/dat_bgc0',k,'.br.dat'
-          OPEN(2*k+401,file=TRIM(ADJUSTL(output_string)),STATUS='replace',Recl=12288)
-       ELSE 
-          WRITE(output_string,'(A17,I2,A7)')'./output/dat_bgc',k,'.bu.dat'
-          OPEN(2*k+400,file=TRIM(ADJUSTL(output_string)),STATUS='replace',Recl=12288)
-          WRITE(output_string,'(A17,I2,A7)')'./output/dat_bgc',k,'.br.dat'
-          OPEN(2*k+401,file=TRIM(ADJUSTL(output_string)),STATUS='replace',Recl=12288)
-       END IF
-    END DO
-
-  END SUBROUTINE output_begin_bgc
-
-END MODULE mo_output
+!>
+!! All things output.
+!!
+!! Used to clean up root.f90 and make it easier to implement changes to the output. 
+!!
+!!
+!! @author Philipp Griewank
+!!
+!!  COPYRIGHT
+!!
+!! This file is part of SAMSIM.
+!!
+!!  SAMSIM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
+!!
+!!  SAMSIM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
+!!
+!!  You should have received a copy of the GNU General Public License along with SAMSIM. If not, see <http://www.gnu.org/licenses/>.
+!
+!!
+!! @par Revision History
+!! Brought from the womb by Philipp Griewank, IMPRS (2010-10-11) \n
+!! add some output values by Niels Fuchs, MPIMET (2017-03-01)
+!!
+MODULE mo_output
+
+  USE mo_parameters
+  IMPLICIT NONE
+
+  PRIVATE
+  PUBLIC :: output_begin,output_begin_bgc,output_raw,output_raw_snow,output_raw_lay,output,output_bgc,output_settings!,output_end
+
+CONTAINS
+  !>
+  !! Settings output.
+  !!
+  !! Writes important values to latter identify run.
+  !!
+  !! @par Revision History
+  !! created by Philipp Griewank, IMPRS (2011-02-12)
+  !!
+  SUBROUTINE output_settings(description,N_top,N_bottom,Nlayer,fl_q_bottom,T_bottom,S_bu_bottom,thick_0,time_out,    &
+    time_total,dt,boundflux_flag,albedo_flag,grav_flag,flush_flag,flood_flag,grav_heat_flag,flush_heat_flag,    &
+    harmonic_flag,prescribe_flag,salt_flag,turb_flag,bottom_flag,tank_flag,precip_flag,bgc_flag,N_bgc,k_snow_flush)
+    INTEGER,         INTENT(in) :: N_top,N_bottom,Nlayer, boundflux_flag,albedo_flag,grav_flag,flush_flag, &
+                                   flood_flag,grav_heat_flag , flush_heat_flag,harmonic_flag,prescribe_flag,salt_flag,turb_flag,  &
+                                   bottom_flag,tank_flag,precip_flag,bgc_flag,N_bgc
+    REAL(wp),        INTENT(in) :: fl_q_bottom,T_bottom,S_bu_bottom,thick_0,time_out,time_total,dt,k_snow_flush
+    CHARACTER*12000, INTENT(in) :: description
+    
+    
+    OPEN(1234,file='./output/dat_settings.dat' ,STATUS='replace')
+    WRITE(1234,*)'################  Description  ###############'
+    WRITE(1234,*) TRIM(description)
+
+       
+    !*************************************************************************************************************************
+    !Layer settings and allocation
+    !*************************************************************************************************************************
+    WRITE(1234,*) '##############  Basic settings  ##############'
+    WRITE(1234,'(A16,F15.3)') 'dt              =',dt
+    WRITE(1234,'(A16,F15.3)') 'thick_0         =',thick_0
+    WRITE(1234,'(A16,F15.3)') 'time_out        =',time_out
+    WRITE(1234,'(A16,F15.3)') 'time_total      =',time_total
+    
+    
+    WRITE(1234,'(A16,F15.3)') 'fl_q_bottom     =',fl_q_bottom
+    WRITE(1234,'(A16,F15.3)') 'T_bottom        =',T_bottom
+    WRITE(1234,'(A16,F15.3)') 'S_bu_bottom     =',S_bu_bottom
+    
+    
+    WRITE(1234,'(A16,I9.0)')  'N_top           =',N_top
+    WRITE(1234,'(A16,I9.0)')  'N_middle        =',Nlayer-N_top-N_bottom
+    WRITE(1234,'(A16,I9.0)')  'N_bottom        =',N_bottom
+    WRITE(1234,'(A16,I9.0)')  'Nlayer          =',Nlayer
+    
+    !*************************************************************************************************************************
+    !Flags
+    !*************************************************************************************************************************
+    WRITE(1234,*) '##################  Flags  ###################'
+    !________________________top heat flux____________
+    WRITE(1234,'(A16,I9.0)')  'boundflux_flag  =',boundflux_flag
+    WRITE(1234,'(A16,I9.0)')  'albedo_flag     =',albedo_flag
+    !________________________brine_dynamics____________
+    WRITE(1234,'(A16,I9.0)')  'grav_flag       =',grav_flag
+    WRITE(1234,'(A16,I9.0)')  'flush_flag      =',flush_flag
+    WRITE(1234,'(A16,I9.0)')  'flood_flag      =',flood_flag
+    WRITE(1234,'(A16,I9.0)')  'grav_heat_flag  =',grav_heat_flag
+    WRITE(1234,'(A16,I9.0)')  'flush_heat_flag =',flush_heat_flag
+    WRITE(1234,'(A16,I9.0)')  'harmonic_flag   =',harmonic_flag
+    WRITE(1234,'(A16,F15.3)')  'k_snow_flush    =',k_snow_flush !< Niels, 2017
+    !________________________Salinity____________
+    WRITE(1234,'(A16,I9.0)')  'prescribe_flag  =',prescribe_flag
+    WRITE(1234,'(A16,I9.0)')  'salt_flag       =',salt_flag
+    !________________________bottom setting______________________
+    WRITE(1234,'(A16,I9.0)')  'turb_flag       =',turb_flag
+    WRITE(1234,'(A16,I9.0)')  'bottom_flag     =',bottom_flag
+    WRITE(1234,'(A16,I9.0)')  'tank_flag       =',tank_flag
+    WRITE(1234,'(A16,I9.0)')  'precip_flag     =',precip_flag
+    WRITE(1234,'(A16,I9.0)')  'bgc_flag        =',bgc_flag
+    WRITE(1234,'(A16,I9.0)')  'N_bgc           =',N_bgc
+
+
+    CLOSE(1234)
+  END SUBROUTINE output_settings
+
+  !>
+  !! Standard output.
+  !!
+  !! For time=n*time_out data is exported.
+  !!
+  !! @par Revision History
+  !! created by Philipp Griewank, IMPRS (2010-10-11)
+  !!
+  SUBROUTINE output(Nlayer,T,psi_s,psi_l,thick,S_bu,ray,format_T,format_psi, &
+       format_thick,format_snow,freeboard,thick_snow,T_snow,psi_l_snow,psi_s_snow,           &
+       energy_stored,freshwater,total_resist,thickness,bulk_salin,                &
+       grav_drain,grav_salt,grav_temp,T2m,T_top,perm,format_perm,flush_v,flush_h,psi_g,melt_thick_output, &
+       grav_heat_flux_down, grav_heat_flux_up, exp_heat, format_melt)
+    INTEGER,                       INTENT(in) :: Nlayer
+    REAL(wp), DIMENSION(Nlayer),   INTENT(in) :: T,psi_s,psi_l,thick,S_bu,perm,flush_v,flush_h,psi_g
+    REAL(wp), DIMENSION(Nlayer-1), INTENT(in) :: ray
+    REAL(wp),                      INTENT(in) :: freeboard,thick_snow,T_snow,psi_l_snow,psi_s_snow,energy_stored,&
+                                                 freshwater,thickness,bulk_salin, &
+                                                 total_resist,grav_drain,grav_salt,grav_temp,T2m,T_top, &
+                                                 grav_heat_flux_down, grav_heat_flux_up, exp_heat
+    REAL(wp), DIMENSION(3),        INTENT(in) :: melt_thick_output  !< Niels, 2017: 1: accumulated melt_thick, 2: accumulated melt_thick_snow, 3: accumulated top ice thickness variations (recheck 3: in mo_layer_dynamics)
+    CHARACTER*12000,               INTENT(in) :: format_T,format_psi,format_thick,format_snow,format_perm,format_melt
+
+    WRITE(30,format_T)     T
+    WRITE(31,format_psi)   psi_s
+    WRITE(32,format_thick) thick
+    WRITE(33,format_T)     S_bu
+    WRITE(34,format_T)     ray
+    WRITE(35,format_psi)   psi_l
+    WRITE(40,'(F9.3)')     freeboard
+    WRITE(41,format_snow)  thick_snow,T_snow,psi_l_snow,psi_s_snow
+    WRITE(42,'(F15.1,2x,F10.5,2x,F10.5,2x,F10.5,2x,F10.5)')energy_stored,freshwater,total_resist,thickness,bulk_salin
+    WRITE(43,'(F9.6,2X,F9.5,2X,F7.3,2X,F13.6)')grav_drain,grav_salt,grav_temp, exp_heat
+    WRITE(45,*)            T2m,T_top
+    WRITE(46,format_perm) perm  !< Niels, 2017 add: output permeability
+    WRITE(47,format_perm) flush_v   !< Niels, 2017 add: output vertical flushing
+    WRITE(48,format_perm) flush_h   !< Niels, 2017 add: output horizontal flushing
+    WRITE(49,format_psi) psi_g  !< Niels, 2017 add: output gas fraction !OBS: not simulated physically in SAMSIM
+    WRITE(50,format_melt) melt_thick_output !< Niels, 2017
+    WRITE(51,format_thick) grav_heat_flux_down !< Niels, 2017
+    WRITE(52,format_thick) grav_heat_flux_up !< Niels, 2017
+  
+  END SUBROUTINE output
+
+  !>
+  !! Standard bgc output.
+  !!
+  !! For time=n*time_out data is exported.
+  !!
+  !! @par Revision History
+  !! created by Philipp Griewank, IMPRS (2014-02-06)
+  !!
+  SUBROUTINE output_bgc(Nlayer,N_active,bgc_bottom,N_bgc,bgc_abs,psi_l,thick,m,format_bgc)
+    INTEGER,                             INTENT(in) :: Nlayer, N_bgc, N_active
+    REAL(wp), DIMENSION(N_bgc),          INTENT(in) :: bgc_bottom
+    REAL(wp), DIMENSION(Nlayer),         INTENT(in) :: psi_l,m,thick
+    REAL(wp), DIMENSION(Nlayer,N_bgc),   INTENT(in) :: bgc_abs
+    REAL(wp), DIMENSION(Nlayer)                     :: bgc_bu,bgc_br
+    CHARACTER*12000,                     INTENT(in) :: format_bgc
+    INTEGER :: k,kk
+    
+    DO k=1,N_bgc
+      DO kk=1,Nlayer
+        IF (kk<=N_active) THEN
+           IF (abs(m(kk)) > tolerance) THEN
+              bgc_bu(kk) = bgc_abs(kk,k)/m(kk)
+              IF (abs(psi_l(kk)) > tolerance .AND. abs(thick(kk)) > tolerance) THEN
+                 bgc_br(kk) = bgc_abs(kk,k)/psi_l(kk)/thick(kk)/rho_l
+              ELSE
+                 bgc_br(kk) = 0000._wp
+              END IF
+           ELSE
+              bgc_bu(kk) = 0._wp
+              bgc_br(kk) = 0._wp
+           END IF
+        ELSE
+              bgc_bu(kk) = bgc_bottom(k)
+              bgc_br(kk) = bgc_bottom(k)
+        END IF
+      END DO
+      WRITE(2*k+400,format_bgc)bgc_bu
+      WRITE(2*k+401,format_bgc)bgc_br
+    END DO
+  
+  END SUBROUTINE output_bgc
+
+  !>
+  !! Output for debugging purposes.
+  !!
+  !! Data for each layer is written out each time step to aid in finding errors or understanding model behavior.
+  !!
+  !! @par Revision History
+  !! created by Philipp Griewank, IMPRS (2010-10-11)
+  !!
+  SUBROUTINE output_raw(Nlayer,N_active,time,T,thick,S_bu,psi_s,psi_l,psi_g)
+    INTEGER,                     INTENT(in) :: Nlayer,N_active
+    INTEGER                                 :: k
+    REAL(wp),                    INTENT(in) :: time
+    REAL(wp), DIMENSION(Nlayer), INTENT(in) :: T,thick,S_bu,psi_s,psi_l,psi_g
+
+    DO k = 1,Nlayer
+       IF ( k<=N_active ) THEN
+          WRITE(k+200, '(F8.4,2x,f10.3,2x,f7.5,2x,f9.5,2x,f4.2,2x,f4.2,2x,f4.2)') &
+               time/(86400._wp),T(k),thick(k),S_bu(k),psi_s(k),psi_l(k),psi_g(k)
+       ELSE
+          WRITE(k+200,'(F6.2,2x,f7.3,2x,f5.3,2x,f8.5,2x,f4.2,2x,f4.2,2x,f4.2)')time/(86400._wp),0.0,0.0,0.0,0.0,0.0,0.0
+       END IF
+    END DO
+
+  END SUBROUTINE output_raw
+
+  !>
+  !! Output for debugging purposes.
+  !!
+  !! Data of snow layer is written out at each time step to aid in finding errors or understanding model behavior.
+  !!
+  !! @par Revision History
+  !! created by Philipp Griewank, IMPRS (2010-10-11)
+  !!
+  SUBROUTINE output_raw_snow(time,T_snow,thick_snow,S_abs_snow,m_snow,psi_s_snow,psi_l_snow,psi_g_snow)
+    REAL(wp), INTENT(in) :: time
+    REAL(wp), INTENT(in) :: T_snow,thick_snow,S_abs_snow,m_snow,psi_s_snow,psi_l_snow,psi_g_snow
+
+    IF (thick_snow>0.0_wp) THEN
+       WRITE(66, '(F8.4,2x,f10.3,2x,f5.3,2x,f4.1,2x,f4.2,2x,f4.2,2x,f4.2)') &
+            time/(30._wp*86400._wp),T_snow,thick_snow,S_abs_snow/MAX(m_snow,0.001_wp),psi_s_snow,psi_l_snow,psi_g_snow
+    ELSE 
+       WRITE(66, '(F8.4,2x,f10.3,2x,f5.3,2x,f4.1,2x,f4.2,2x,f4.2,2x,f4.2)') &
+            time/(30._wp*86400._wp),0.0,0.0,0.0,0.0,0.0,0.0
+    END IF
+
+  END SUBROUTINE output_raw_snow
+
+  !>
+  !! Output for debugging layer dynamics..
+  !!
+  !! Is used when debug_flag = 2 to track when which layer dynamics occur (see mo_layer_dynamics).
+  !!
+  !!
+
+  SUBROUTINE output_raw_lay(Nlayer,N_active,H_abs,m,S_abs,thick,string)
+    INTEGER,                     INTENT(in) :: Nlayer,N_active
+    INTEGER                                 :: k
+    REAL(wp), DIMENSION(Nlayer), INTENT(in) :: H_abs,S_abs,thick,m 
+    REAL(wp)                                :: mm 
+    CHARACTER*6,                 INTENT(in) :: string
+
+    DO k = 1,Nlayer
+       IF ( k<=N_active ) THEN
+          IF (abs(m(k)) < tolerance) THEN
+             mm = 99999999._wp
+          ELSE
+             mm = m(k)
+          END IF
+          WRITE(k+200, '(A6,2x,f11.0,2x,f6.3,2x,f8.5,2x,f9.3,2x,I2)') &
+               string,h_abs(k),thick(k),s_abs(k)/mm,mm/MAX(thick(k),0.0000000000000000001_wp),N_active
+       ELSE 
+          WRITE(k+200, '(A6,2x,f11.0,2x,f6.3,2x,f8.1,2x,f9.3,2x,I2)')string,0.0,0.0,0.0,0.0,N_active
+       END IF
+    END DO
+  END SUBROUTINE output_raw_lay
+
+  !>
+  !! Output files are opened and format strings are created.
+  !!
+  !! Format strings are defined according to the number of layers used which define the output format.
+  !! Files are opened.
+  !!
+  !! @par Revision History
+  !! created by Philipp Griewank, IMPRS (2010-10-11)
+  !! moved by Philipp Griewank, IMPRS (2011-03-09)
+  !!
+  SUBROUTINE output_begin(Nlayer,debug_flag,format_T,format_psi,format_thick,format_snow,format_T2m_top,format_perm,&
+                          &format_melt)
+    INTEGER,         INTENT(in)  :: Nlayer,debug_flag
+    CHARACTER*12000, INTENT(out) :: format_T,format_psi,format_thick,format_snow,format_T2m_top,format_perm,&
+                                    &format_melt
+    INTEGER                      :: k
+    CHARACTER*12000              :: output_string
+
+
+    IF (debug_flag==2) THEN
+       DO k=1,Nlayer
+          IF (k<10) THEN
+             WRITE(output_string,'(A18,I1,A4)')'./output/thermo00',k,'.txt'
+             OPEN(k+200,file=TRIM(ADJUSTL(output_string)),STATUS='replace',Recl=12288)
+          ELSE IF (k<100) THEN
+             WRITE(output_string,'(A17,I2,A4)')'./output/thermo0',k,'.txt'
+             OPEN(k+200,file=TRIM(ADJUSTL(output_string)),STATUS='replace',Recl=12288)
+          ELSE
+             WRITE(output_string,'(A16,I3,A4)')'./output/thermo',k,'.txt'
+             OPEN(k+200,file=TRIM(ADJUSTL(output_string)),STATUS='replace',Recl=12288)
+          END IF
+       END DO
+    END IF
+
+    format_T     = '(F9.3,2x'
+    format_perm  = '(ES14.7,2x'
+    format_psi   = '(F9.3,2x'
+    format_thick = '(F9.5,2x'
+
+    format_snow    = '(F9.3,2x,F9.3,2x,F9.3,2x,F9.3)'
+    format_T2m_top = '(F9.3,2x,F9.3,2x,F9.3)'
+    format_melt    = '(ES14.7,2x,ES14.7,2x,ES14.7)'
+
+    DO k=1,Nlayer
+       format_T     = TRIM(format_T)//',f9.3,2x'
+       format_psi   = TRIM(format_psi)//',f9.3,2x'
+       format_thick = TRIM(format_thick)//',f9.5,2x'
+       format_perm  = TRIM(format_perm)//',ES14.7,2x'
+    END DO
+
+    format_T     = TRIM(format_T)//')'
+    format_psi   = TRIM(format_psi)//')'
+    format_thick = TRIM(format_thick)//')'
+    format_perm  = TRIM(format_perm)//')'
+
+    
+    OPEN(30,file='./output/dat_T.dat'           ,STATUS='replace',Recl=12288)
+    OPEN(31,file='./output/dat_psi_s.dat'       ,STATUS='replace',Recl=12288)
+    OPEN(32,file='./output/dat_thick.dat'       ,STATUS='replace',Recl=12288)
+    OPEN(33,file='./output/dat_S_bu.dat'        ,STATUS='replace',Recl=12288)
+    OPEN(34,file='./output/dat_ray.dat'         ,STATUS='replace',Recl=12288)
+    OPEN(35,file='./output/dat_psi_l.dat'       ,STATUS='replace',Recl=12288)
+    OPEN(40,file='./output/dat_freeboard.dat'   ,STATUS='replace',Recl=12288)
+    OPEN(41,file='./output/dat_snow.dat'        ,STATUS='replace',Recl=12288)
+    OPEN(42,file='./output/dat_vital_signs.dat' ,STATUS='replace',Recl=12288)
+    OPEN(43,file='./output/dat_grav_drain.dat'  ,STATUS='replace',Recl=12288)
+    OPEN(45,file='./output/dat_T2m_T_top.dat'   ,STATUS='replace',Recl=12288)
+    OPEN(46,file='./output/dat_perm.dat'        ,STATUS='replace',Recl=12288)
+    OPEN(47,file='./output/dat_flush_v.dat'        ,STATUS='replace',Recl=12288)
+    OPEN(48,file='./output/dat_flush_h.dat'        ,STATUS='replace',Recl=12288)
+    OPEN(49,file='./output/dat_psi_g.dat'        ,STATUS='replace',Recl=12288)
+    OPEN(50,file='./output/dat_melt.dat'        ,STATUS='replace',Recl=12288)
+    OPEN(51,file='./output/dat_grav_heat_flux_down.dat'        ,STATUS='replace',Recl=12288)
+    OPEN(52,file='./output/dat_grav_heat_flux_up.dat'        ,STATUS='replace',Recl=12288)
+
+    OPEN(66,file='./output/snow.txt'            ,STATUS='replace')
+
+  END SUBROUTINE output_begin
+  
+  
+  !>
+  !! Output files for bgc are opened and format strings are created.
+  !!
+  !! Same thing as out_begin but for bgc
+  !! Each tracer is outputted in bulk and in brine concentration in a separate file.
+  !! Added ADJUSTL to the output strings because they got wierd 
+  !!
+  !! @par Revision History
+  !! created by Dr. Philipp Griewank, MPI (2014-02-07)
+  !! fix by Dr. Philipp Griewank, UniK (2018-05-18)
+  SUBROUTINE output_begin_bgc(Nlayer,N_bgc,format_bgc)
+    INTEGER,         INTENT(in)  :: Nlayer,N_bgc
+    CHARACTER*12000, INTENT(out) :: format_bgc
+    INTEGER                      :: k
+    CHARACTER*12000              :: output_string
+
+     
+
+    format_bgc = '(F16.8,2x'
+
+    DO k=1,Nlayer
+      format_bgc = TRIM(format_bgc)//',f16.8,2x'
+    END DO
+    
+    format_bgc = TRIM(format_bgc)//')'
+
+    DO k=1,N_bgc
+       IF (k<10) THEN
+          WRITE(output_string,'(A18,I1,A7)')'./output/dat_bgc0',k,'.bu.dat'
+          OPEN(2*k+400,file=TRIM(ADJUSTL(output_string)),STATUS='replace',Recl=12288)
+          WRITE(output_string,'(A18,I1,A7)')'./output/dat_bgc0',k,'.br.dat'
+          OPEN(2*k+401,file=TRIM(ADJUSTL(output_string)),STATUS='replace',Recl=12288)
+       ELSE 
+          WRITE(output_string,'(A17,I2,A7)')'./output/dat_bgc',k,'.bu.dat'
+          OPEN(2*k+400,file=TRIM(ADJUSTL(output_string)),STATUS='replace',Recl=12288)
+          WRITE(output_string,'(A17,I2,A7)')'./output/dat_bgc',k,'.br.dat'
+          OPEN(2*k+401,file=TRIM(ADJUSTL(output_string)),STATUS='replace',Recl=12288)
+       END IF
+    END DO
+
+  END SUBROUTINE output_begin_bgc
+
+END MODULE mo_output
-- 
GitLab