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 85d241f5572f1794ac85a0d0b0a0d1f6ccac248a..0d4b94ef71663a3d6c31be573de0f353596eb363 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 14bda5a0b81279900a2cabccd2634df2789228f7..20a52b82ac3d030153a3337e0dc6e5d3e6db40b7 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 2502bd00734b2b3e3ef266f10d34fba75a95979f..c1d294398dd64efe383cc5acd19f1e00ffc49ecb 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 390eb2dadc8f2ab6b5964e264b63861424fb7a0d..88142c6c04d992254f2a5ea0de08fee48fca9d60 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 0aab77eae56481b3293d8e1da07ccf9e73e98a30..0000000000000000000000000000000000000000
--- 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 996242548d6c28ad6ab449d60b4d0cacc243c6c2..967f429309828b829c10330cf8e419b060366c22 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 8a4372d34a766b8be8f5639cda6539b85f3c9c8c..68169dc7febe1cc9dc857fde6459f91085e2577c 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 98b3901909442f52bcccdc19dd4ded4d33d4073e..3b495cc852c5c0bc8a5c858362c45057cb5e625e 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 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/Run_specifics/config.json b/Run_specifics/config.json
index 8bdcc22482567f19f683ab34edc6e762a2eb9363..a2c92d2d1a118f3c7c449cfa45f29e1b726a2759 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 fd739a7179d8b19c7bbf7aae2ccfc2cf65f3e747..7badf450dcb4c4c8c2cccb624627e84b467dbd8f 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 280fe10eeaa95c4633aa850837e472735cda1244..fbb19f2f1fc21d0ba9a63ceafdff54b56bac120f 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 b68f0a92d16b90040dd13ed550f4b486bd60f640..26da0962b51df3ffa0dcaa9ea996260d88045455 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 1f26ba1bfcf524a2cf2328436d64e2c25cbe6ea7..008deed85c95d3d2c4f63fe0c9370f5122bf5e38 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 92cf3d8c2a97866690b2141d28fb54b56141293a..6393e62d2a68c44d24a42a2b3470ee1d4bdcb2ed 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