diff --git a/mo_data.f90 b/mo_data.f90
index fbb19f2f1fc21d0ba9a63ceafdff54b56bac120f..335e8446656c267166789621606fb42fc2523aa9 100644
--- a/mo_data.f90
+++ b/mo_data.f90
@@ -1,219 +1,220 @@
-!>
-!! Sets data and contains all flag descriptions.
-!!
-!! All data needed by mo_grotz are set in this module.
-!! Most arrays are allocated after the needed dimension is specified for each testcase in mo_init.f90.
-!!
-!!
-!! @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
-!! Initialized by Philipp Griewank, IMPRS (2010-07-14) \n
-!! Add several variables by Niels Fuchs, MPIMET (2017-03-01)
-
-!!
-MODULE mo_data
-
- USE mo_parameters, ONLY:wp
-
- IMPLICIT NONE
-
-
- !----Arrays
- REAL(wp), DIMENSION(:), ALLOCATABLE :: H !< Enthalpy [J]
- REAL(wp), DIMENSION(:), ALLOCATABLE :: H_abs !< specific Enthalpy [J/kg]
- REAL(wp), DIMENSION(:), ALLOCATABLE :: Q !< Heat in layer [J]
- REAL(wp), DIMENSION(:), ALLOCATABLE :: fl_Q !< Heat flux between layers [J/s]
- REAL(wp), DIMENSION(:), ALLOCATABLE :: T !< Temperature [C]
- REAL(wp), DIMENSION(:), ALLOCATABLE :: S_bu !< Bulk Salinity [g/kg]
- REAL(wp), DIMENSION(:), ALLOCATABLE :: fl_S !< Salinity flux [(g/s]
- REAL(wp), DIMENSION(:), ALLOCATABLE :: S_abs !< Absolute Salinity [g]
- REAL(wp), DIMENSION(:), ALLOCATABLE :: S_br !< Brine salinity [g/kg]
- REAL(wp), DIMENSION(:), ALLOCATABLE :: thick !< Layer thickness [m]
- REAL(wp), DIMENSION(:), ALLOCATABLE :: m !< Mass [kg]
- REAL(wp), DIMENSION(:), ALLOCATABLE :: fl_m !< Mass fluxes between layers [kg]
- REAL(wp), DIMENSION(:), ALLOCATABLE :: V_s !< Volume [m^3] of solid
- REAL(wp), DIMENSION(:), ALLOCATABLE :: V_l !< Volume [m^3] of liquid
- REAL(wp), DIMENSION(:), ALLOCATABLE :: V_g !< Volume [m^3] of gas
- REAL(wp), DIMENSION(:), ALLOCATABLE :: V_ex !< Volume of brine due expelled due to freezing [m^3] of solid, gas & liquid
- REAL(wp), DIMENSION(:), ALLOCATABLE :: phi !< Solid mass fraction
- REAL(wp), DIMENSION(:), ALLOCATABLE :: psi_s !< Solid volume fraction
- REAL(wp), DIMENSION(:), ALLOCATABLE :: psi_l !< Liquid volume fraction
- REAL(wp), DIMENSION(:), ALLOCATABLE :: psi_g !< Gas volume fraction
- REAL(wp), DIMENSION(:), ALLOCATABLE :: ray !< Rayleigh number of each layer
- REAL(wp), DIMENSION(:), ALLOCATABLE :: perm,flush_v,flush_h,flush_v_old,flush_h_old !< Permeability [?]
-
- REAL(wp) :: dt !< Timestep [s]
- REAL(wp) :: timestep_data !< Timestep of input data [s]
- REAL(wp) :: thick_0 !< Initial layer thickness [m]
- REAL(wp) :: time !< Time [s]
- REAL(wp) :: freeboard !< Height of ice surface above (or below) waterlevel [m]
- REAL(wp) :: T_freeze !< Freezing temperature [C]
- INTEGER :: Nlayer !< Number of layers
- INTEGER :: N_bottom !< Number of bottom layers
- INTEGER :: N_middle !< Number of middle layers
- INTEGER :: N_top !< Number of top layers
- INTEGER :: N_active !< Number of Layers active in the present
- INTEGER :: i !< Index, normally used for time
- INTEGER :: k !< Index, normally used for layer
- INTEGER :: styropor_flag
- REAL(wp) :: time_out !< Time between outputs [s]
- REAL(wp) :: time_total !< Time of simulation [s]
- INTEGER :: i_time !< Number of timesteps
- INTEGER :: i_time_out !< Number of timesteps between each output
- INTEGER :: n_time_out !< Counts number of timesteps between output
- CHARACTER*12000 :: format_T,format_psi,format_thick,format_snow,format_integer,format_T2m_top,format_bgc,&
- &format_melt !< Format strings for output. Niels(2017) add: melt output
- CHARACTER*12000 :: format_perm !< Niels(2017) add: permeability output
-
- !----Boundary conditions
- REAL(wp) :: T_bottom !< Temperature of water beneath the ice [C]
- REAL(wp) :: T_top !< Temperature at the surface [C]
- REAL(wp) :: S_bu_bottom !< Salinity beneath the ice [g/kg]
- REAL(wp) :: T2m !< Two meter Temperature [C]
- REAL(wp) :: fl_q_bottom !< Bottom heat flux [J*s]
-
- !----Snow & Precip
- REAL(wp) :: psi_s_snow !< Solid volume fraction of snow layer
- REAL(wp) :: psi_l_snow !< Liquid volume fraction of snow layer
- REAL(wp) :: psi_g_snow !< Gas volume fraction of snow layer
- REAL(wp) :: phi_s !< Solid mass fraction of snow layer
- REAL(wp) :: S_abs_snow !< Absolute salinity of snow layer [g]
- REAL(wp) :: H_abs_snow !< Absolute enthalpy of snow layer [J]
- REAL(wp) :: m_snow !< Mass of snow layer [kg]
- REAL(wp) :: T_snow !< Temperature of snow layer [C]
- REAL(wp) :: thick_snow,test !< Thickness of snow layer [m]
- REAL(wp) :: liquid_precip !< Liquid precip, [meter of water/s]
- REAL(wp) :: solid_precip !< Solid precip, [meter of water /s]
- REAL(wp) :: fl_q_snow !< flow of heat into the snow layer
-
- !----Vital signs
- REAL(wp) :: energy_stored !< Total amount of energy stored, control is freezing point temperature of S_bu_bottom [J]
- REAL(wp) :: total_resist !< Thermal resistance of the whole column []
- REAL(wp) :: surface_water !< Percentage of water fraction in the top 5cm [%]
- REAL(wp) :: freshwater !< Meters of freshwater stored in column [m]
- REAL(wp) :: thickness !< Meters of ice [m]
- REAL(wp) :: bulk_salin !< Salt/Mass [ppt]
-
- !----Model and numerics
- REAL(wp) :: thick_min !< Parameter for snow, determines when snow is in thermal equilibrium with the ice and when it is totally neglected
- REAL(wp),SAVE :: T_test !< First guess for getT subroutine
-
- !----Radiation
- REAL(wp) :: albedo !< Amount of short wave radiation which is reflected at the top surface
- REAL(wp) :: fl_sw !< Incoming shortwave radiation [W/m**2]
- REAL(wp) :: fl_lw !< Incoming longwave radiation [W/m**2]
- REAL(wp) :: fl_sen !< Sensitive heat flux [W/m**2]
- REAL(wp) :: fl_lat !< Latent heat flux [W/m**2]
- REAL(wp) :: fl_rest !< Bundled longwave,sensitive and latent heat flux [W/m**2]
- REAL(wp), DIMENSION(:), ALLOCATABLE :: fl_rad !< Energy flux of absorbed sw radiation of each layer [J/s]
-
- !----Gravity drainage
- REAL(wp) :: grav_drain !< brine flux of gravity drainage between two outputs [kg/s]
- REAL(wp) :: grav_salt !< salt flux moved by gravity drainage between two outputs [kg*ppt/s]
- REAL(wp) :: grav_temp !< average temperature of gravity drainage brine between two outputs [T]
-
- !----Flushing
- REAL(wp) :: melt_thick !< thickness of fully liquid part of top layer [m]
- REAL(wp) :: melt_thick_snow, melt_thick_snow_old !< Niels(2017) add: thickness of excess fully liquid part from snow_melt_processes [m]
- REAL(wp), DIMENSION(3) :: melt_thick_output !< Niels, 2017 add: output field of surface liquid meltwater sizes
-
- !----Lab fluxes
- REAL(wp) :: alpha_flux_instable !< Proportionality constant which determines energy flux by the temperature difference T_top>T2m [W/C]
- REAL(wp) :: alpha_flux_stable !< Proportionality constant which determines energy flux by the temperature difference T_top<T2m [W/C]
-
- !----Flags
- INTEGER :: grav_flag !< 1: no gravity drainage, 2: Gravity drainage, 3: Simple Drainage
- INTEGER :: prescribe_flag !< 1: nothing happens, 2: prescribed Salinity profile is prescribed at each timestep (does not disable brine dynamics, just overwrites the salinity!)
- INTEGER :: grav_heat_flag !< 1: nothing happens, 2: compensates heatfluxes in grav_flag = 2
- INTEGER :: flush_heat_flag !< 1: nothing happens, 2: compensates heatfluxes in flush_flag = 5
- INTEGER :: turb_flag !< 1: No bottom turbulence, 2: Bottom mixing
- INTEGER :: salt_flag !< 1: Sea salt, 2: NaCL
- INTEGER :: boundflux_flag !< 1: top and bottom cooling plate, 2:top Notz fluxes, bottom cooling plate 3: top flux=a*(T-T_s)
- INTEGER :: flush_flag !< 1: no flushing, 4:meltwater is removed artificially, 5:vert and horiz flushing, 6: simplified
- INTEGER :: flood_flag !< 1: no flooding, 2:normal flooding, 3:simple flooding
- INTEGER :: bottom_flag !< 1: nothing changes, 2: deactivates all bottom layer dynamics, useful for some debugging and idealized tests
- INTEGER :: debug_flag !< 1: no raw layer output, 2: each layer is output at every timestep (warning, file size can be very large)
- INTEGER :: precip_flag !< 0: solid and liquid precipitation, 1:phase determined by T2m
- INTEGER :: harmonic_flag !< 1: minimal permeability is used to calculate Rayleigh number, 2:harmonic mean is used for Rayleigh number
- INTEGER :: tank_flag !< 1: nothing, 2: S_bu_bottom and bgc_bottom are calculated as if the experiment is conducted in a tank
- INTEGER :: lab_snow_flag !< Niels, 2017 add: 0: lab setup without snow covers, 1: lab setup include snow influence on heat fluxes
- INTEGER :: albedo_flag !< 1: simple albedo, 2: normal albedo, see func_albedo for details
- INTEGER :: freeboard_snow_flag !< Niels, 2017 add: 0: respect the mass of snow in the freeboard calculation, 1: don't
- INTEGER :: snow_flush_flag !< Niels, 2017 add: 0: all meltwater from snow forms slush, 1: meltwater partly leads to flushing, ratio defined by "k_snow_flush"
- INTEGER :: initial_state_flag !< Jakob 2022 add: 1: default behaviour - starts from ice thickness = 0, 2: H_abs, S_abs, m, thick of initial ice given
-
-
- !##########################################################################################
- !Variables used to import data
- !##########################################################################################
- INTEGER :: Length_Input !< Sets the input length for atmoflux_flag==2, common value of 13169
- REAL(wp), DIMENSION(:), ALLOCATABLE :: Tinput !< Niels, 2017 add: used to read in top temperature for field experiment tests, dimension needs to be set in the code
- !REAL(wp), DIMENSION(:), ALLOCATABLE :: precipinput !< Niels, 2017 add: used to read in precipation for field experiment tests, dimension needs to be set in the code
- REAL(wp), DIMENSION(:), ALLOCATABLE :: ocean_T_input !< Niels, 2017 add: used to read in ocean temperature for field experiment tests, dimension needs to be set in the code
- REAL(wp), DIMENSION(:), ALLOCATABLE :: ocean_flux_input !< Niels, 2017 add: used to read in oceanic heat flux for field experiment tests, dimension needs to be set in the code
- REAL(wp), DIMENSION(:), ALLOCATABLE :: styropor_input !< Niels, 2017 add: if styropor is used in the lab on top of the ice to simulate snow heat fluxes
- REAL(wp), DIMENSION(:), ALLOCATABLE :: Ttop_input !< Niels, 2017 add: used for testcase 111, comparison with greenland harp data, uppermost harp temperature is seen as Ttop
- REAL(wp), DIMENSION(:), ALLOCATABLE :: fl_sw_input !< Used to read in sw fluxes
- REAL(wp), DIMENSION(:), ALLOCATABLE :: fl_lw_input !< Used to read in lw fluxes
- REAL(wp), DIMENSION(:), ALLOCATABLE :: fl_sen_input !< Used to read in sensible heat fluxes
- REAL(wp), DIMENSION(:), ALLOCATABLE :: fl_lat_input !< Used to read in latent heat fluxes
- REAL(wp), DIMENSION(:), ALLOCATABLE :: T2m_input !< Used to read in 2Tm
- REAL(wp), DIMENSION(:), ALLOCATABLE :: precip_l_input !< Used to read in liquid precipitation
- REAL(wp), DIMENSION(:), ALLOCATABLE :: precip_s_input !< Used to read in solid precipitation
- REAL(wp), DIMENSION(:), ALLOCATABLE :: T_top_input !< Used to read in T_top
- REAL(wp), DIMENSION(:), ALLOCATABLE :: T_bottom_input !< Used to read in T_bottom
- REAL(wp), DIMENSION(:), ALLOCATABLE :: S_bu_bottom_input !< Used to read in S_bu_bottom
- REAL(wp), DIMENSION(:), ALLOCATABLE :: fl_q_bottom_input !< Used to read in fl_q_bottom
- REAL(wp), DIMENSION(:), ALLOCATABLE :: time_input !< Used to read in time from ERA for atmoflux_flag==2
- INTEGER :: time_counter !< Keeps track of input data
- INTEGER :: days !< Days of simulation, specified in init.py
-
-
- !REAL(wp), DIMENSION(:), ALLOCATABLE :: perm !< permeability values
-
-
- !##########################################################################################
- !Chemicals Baby!
- !All arrays needed to support bigogeochemical tracers
- !Chemical matrixes: index 1 defines the chemical, index 2 the layer
- !##########################################################################################
- INTEGER :: bgc_flag !< 1: no bgc, 2:bgc
- INTEGER :: N_bgc !< Number of chemicals
- REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fl_brine_bgc !< Brine fluxes in a matrix, [kg/s], first index is the layer of origin, and the second index is the layer of arrival
- REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bgc_abs !< Absolute amount of chemicals [kmol] for each tracer
- REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bgc_bu !< Bulk amounts of chemicals [kmol/kg]
- REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bgc_br !< Brine concentrations of chems [kmol/kg]
- REAL(wp), DIMENSION(:), ALLOCATABLE :: bgc_bottom !< Bulk concentrations of chems below the ice [kmol/kg]
-
-
- !##########################################################################################
- !Variables needed for tank experiments in which concentrations below the ice change over time
- !##########################################################################################
- REAL(wp), DIMENSION(:), ALLOCATABLE :: bgc_total !< Total of chems, for lab experiments with a fixed total amount
- REAL(wp) :: m_total !< Total initial water mass, for lab experiments with a fixed total amount
- REAL(wp) :: S_total !< Total initial salt mass, for lab experiments with a fixed total amount
- REAL(wp) :: tank_depth !< water depth in meters, used to calculate concentrations below ice for tank experiments
-
- !#######
- !Misc
- !#######
- CHARACTER*3 :: flush_question='No!' !< Niels, 2017 add: used to indicate in stdout wether flushing occurs at this moment or not
- REAL(wp) :: melt_err=0._wp !< Niels, 2017 add: used to check how much meltwater vanishes in flushing routine
- 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
-
+!>
+!! Sets data and contains all flag descriptions.
+!!
+!! All data needed by mo_grotz are set in this module.
+!! Most arrays are allocated after the needed dimension is specified for each testcase in mo_init.f90.
+!!
+!!
+!! @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
+!! Initialized by Philipp Griewank, IMPRS (2010-07-14) \n
+!! Add several variables by Niels Fuchs, MPIMET (2017-03-01)
+!! edited by Niels Fuchs, UHH (2024-05-03), corrected Enthalpy declaration (before specific and total enthalpy were confused)
+
+!!
+MODULE mo_data
+
+ USE mo_parameters, ONLY:wp
+
+ IMPLICIT NONE
+
+
+ !----Arrays
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: H !< specific Enthalpy [J/kg]
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: H_abs !< total Enthalpy [J]
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: Q !< Heat in layer [J]
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: fl_Q !< Heat flux between layers [J/s]
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: T !< Temperature [C]
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: S_bu !< Bulk Salinity [g/kg]
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: fl_S !< Salinity flux [(g/s]
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: S_abs !< Absolute Salinity [g]
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: S_br !< Brine salinity [g/kg]
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: thick !< Layer thickness [m]
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: m !< Mass [kg]
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: fl_m !< Mass fluxes between layers [kg]
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: V_s !< Volume [m^3] of solid
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: V_l !< Volume [m^3] of liquid
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: V_g !< Volume [m^3] of gas
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: V_ex !< Volume of brine due expelled due to freezing [m^3] of solid, gas & liquid
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: phi !< Solid mass fraction
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: psi_s !< Solid volume fraction
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: psi_l !< Liquid volume fraction
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: psi_g !< Gas volume fraction
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: ray !< Rayleigh number of each layer
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: perm,flush_v,flush_h,flush_v_old,flush_h_old !< Permeability [?]
+
+ REAL(wp) :: dt !< Timestep [s]
+ REAL(wp) :: timestep_data !< Timestep of input data [s]
+ REAL(wp) :: thick_0 !< Initial layer thickness [m]
+ REAL(wp) :: time !< Time [s]
+ REAL(wp) :: freeboard !< Height of ice surface above (or below) waterlevel [m]
+ REAL(wp) :: T_freeze !< Freezing temperature [C]
+ INTEGER :: Nlayer !< Number of layers
+ INTEGER :: N_bottom !< Number of bottom layers
+ INTEGER :: N_middle !< Number of middle layers
+ INTEGER :: N_top !< Number of top layers
+ INTEGER :: N_active !< Number of Layers active in the present
+ INTEGER :: i !< Index, normally used for time
+ INTEGER :: k !< Index, normally used for layer
+ INTEGER :: styropor_flag
+ REAL(wp) :: time_out !< Time between outputs [s]
+ REAL(wp) :: time_total !< Time of simulation [s]
+ INTEGER :: i_time !< Number of timesteps
+ INTEGER :: i_time_out !< Number of timesteps between each output
+ INTEGER :: n_time_out !< Counts number of timesteps between output
+ CHARACTER*12000 :: format_T,format_psi,format_thick,format_snow,format_integer,format_T2m_top,format_bgc,&
+ &format_melt !< Format strings for output. Niels(2017) add: melt output
+ CHARACTER*12000 :: format_perm !< Niels(2017) add: permeability output
+
+ !----Boundary conditions
+ REAL(wp) :: T_bottom !< Temperature of water beneath the ice [C]
+ REAL(wp) :: T_top !< Temperature at the surface [C]
+ REAL(wp) :: S_bu_bottom !< Salinity beneath the ice [g/kg]
+ REAL(wp) :: T2m !< Two meter Temperature [C]
+ REAL(wp) :: fl_q_bottom !< Bottom heat flux [J*s]
+
+ !----Snow & Precip
+ REAL(wp) :: psi_s_snow !< Solid volume fraction of snow layer
+ REAL(wp) :: psi_l_snow !< Liquid volume fraction of snow layer
+ REAL(wp) :: psi_g_snow !< Gas volume fraction of snow layer
+ REAL(wp) :: phi_s !< Solid mass fraction of snow layer
+ REAL(wp) :: S_abs_snow !< Absolute salinity of snow layer [g]
+ REAL(wp) :: H_abs_snow !< Absolute enthalpy of snow layer [J]
+ REAL(wp) :: m_snow !< Mass of snow layer [kg]
+ REAL(wp) :: T_snow !< Temperature of snow layer [C]
+ REAL(wp) :: thick_snow,test !< Thickness of snow layer [m]
+ REAL(wp) :: liquid_precip !< Liquid precip, [meter of water/s]
+ REAL(wp) :: solid_precip !< Solid precip, [meter of water /s]
+ REAL(wp) :: fl_q_snow !< flow of heat into the snow layer
+
+ !----Vital signs
+ REAL(wp) :: energy_stored !< Total amount of energy stored, control is freezing point temperature of S_bu_bottom [J]
+ REAL(wp) :: total_resist !< Thermal resistance of the whole column []
+ REAL(wp) :: surface_water !< Percentage of water fraction in the top 5cm [%]
+ REAL(wp) :: freshwater !< Meters of freshwater stored in column [m]
+ REAL(wp) :: thickness !< Meters of ice [m]
+ REAL(wp) :: bulk_salin !< Salt/Mass [ppt]
+
+ !----Model and numerics
+ REAL(wp) :: thick_min !< Parameter for snow, determines when snow is in thermal equilibrium with the ice and when it is totally neglected
+ REAL(wp),SAVE :: T_test !< First guess for getT subroutine
+
+ !----Radiation
+ REAL(wp) :: albedo !< Amount of short wave radiation which is reflected at the top surface
+ REAL(wp) :: fl_sw !< Incoming shortwave radiation [W/m**2]
+ REAL(wp) :: fl_lw !< Incoming longwave radiation [W/m**2]
+ REAL(wp) :: fl_sen !< Sensitive heat flux [W/m**2]
+ REAL(wp) :: fl_lat !< Latent heat flux [W/m**2]
+ REAL(wp) :: fl_rest !< Bundled longwave,sensitive and latent heat flux [W/m**2]
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: fl_rad !< Energy flux of absorbed sw radiation of each layer [J/s]
+
+ !----Gravity drainage
+ REAL(wp) :: grav_drain !< brine flux of gravity drainage between two outputs [kg/s]
+ REAL(wp) :: grav_salt !< salt flux moved by gravity drainage between two outputs [kg*ppt/s]
+ REAL(wp) :: grav_temp !< average temperature of gravity drainage brine between two outputs [T]
+
+ !----Flushing
+ REAL(wp) :: melt_thick !< thickness of fully liquid part of top layer [m]
+ REAL(wp) :: melt_thick_snow, melt_thick_snow_old !< Niels(2017) add: thickness of excess fully liquid part from snow_melt_processes [m]
+ REAL(wp), DIMENSION(3) :: melt_thick_output !< Niels, 2017 add: output field of surface liquid meltwater sizes
+
+ !----Lab fluxes
+ REAL(wp) :: alpha_flux_instable !< Proportionality constant which determines energy flux by the temperature difference T_top>T2m [W/C]
+ REAL(wp) :: alpha_flux_stable !< Proportionality constant which determines energy flux by the temperature difference T_top<T2m [W/C]
+
+ !----Flags
+ INTEGER :: grav_flag !< 1: no gravity drainage, 2: Gravity drainage, 3: Simple Drainage
+ INTEGER :: prescribe_flag !< 1: nothing happens, 2: prescribed Salinity profile is prescribed at each timestep (does not disable brine dynamics, just overwrites the salinity!)
+ INTEGER :: grav_heat_flag !< 1: nothing happens, 2: compensates heatfluxes in grav_flag = 2
+ INTEGER :: flush_heat_flag !< 1: nothing happens, 2: compensates heatfluxes in flush_flag = 5
+ INTEGER :: turb_flag !< 1: No bottom turbulence, 2: Bottom mixing
+ INTEGER :: salt_flag !< 1: Sea salt, 2: NaCL
+ INTEGER :: boundflux_flag !< 1: top and bottom cooling plate, 2:top Notz fluxes, bottom cooling plate 3: top flux=a*(T-T_s)
+ INTEGER :: flush_flag !< 1: no flushing, 4:meltwater is removed artificially, 5:vert and horiz flushing, 6: simplified
+ INTEGER :: flood_flag !< 1: no flooding, 2:normal flooding, 3:simple flooding
+ INTEGER :: bottom_flag !< 1: nothing changes, 2: deactivates all bottom layer dynamics, useful for some debugging and idealized tests
+ INTEGER :: debug_flag !< 1: no raw layer output, 2: each layer is output at every timestep (warning, file size can be very large)
+ INTEGER :: precip_flag !< 0: solid and liquid precipitation, 1:phase determined by T2m
+ INTEGER :: harmonic_flag !< 1: minimal permeability is used to calculate Rayleigh number, 2:harmonic mean is used for Rayleigh number
+ INTEGER :: tank_flag !< 1: nothing, 2: S_bu_bottom and bgc_bottom are calculated as if the experiment is conducted in a tank
+ INTEGER :: lab_snow_flag !< Niels, 2017 add: 0: lab setup without snow covers, 1: lab setup include snow influence on heat fluxes
+ INTEGER :: albedo_flag !< 1: simple albedo, 2: normal albedo, see func_albedo for details
+ INTEGER :: freeboard_snow_flag !< Niels, 2017 add: 0: respect the mass of snow in the freeboard calculation, 1: don't
+ INTEGER :: snow_flush_flag !< Niels, 2017 add: 0: all meltwater from snow forms slush, 1: meltwater partly leads to flushing, ratio defined by "k_snow_flush"
+ INTEGER :: initial_state_flag !< Jakob 2022 add: 1: default behaviour - starts from ice thickness = 0, 2: H_abs, S_abs, m, thick of initial ice given
+
+
+ !##########################################################################################
+ !Variables used to import data
+ !##########################################################################################
+ INTEGER :: Length_Input !< Sets the input length for atmoflux_flag==2, common value of 13169
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: Tinput !< Niels, 2017 add: used to read in top temperature for field experiment tests, dimension needs to be set in the code
+ !REAL(wp), DIMENSION(:), ALLOCATABLE :: precipinput !< Niels, 2017 add: used to read in precipation for field experiment tests, dimension needs to be set in the code
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: ocean_T_input !< Niels, 2017 add: used to read in ocean temperature for field experiment tests, dimension needs to be set in the code
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: ocean_flux_input !< Niels, 2017 add: used to read in oceanic heat flux for field experiment tests, dimension needs to be set in the code
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: styropor_input !< Niels, 2017 add: if styropor is used in the lab on top of the ice to simulate snow heat fluxes
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: Ttop_input !< Niels, 2017 add: used for testcase 111, comparison with greenland harp data, uppermost harp temperature is seen as Ttop
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: fl_sw_input !< Used to read in sw fluxes
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: fl_lw_input !< Used to read in lw fluxes
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: fl_sen_input !< Used to read in sensible heat fluxes
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: fl_lat_input !< Used to read in latent heat fluxes
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: T2m_input !< Used to read in 2Tm
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: precip_l_input !< Used to read in liquid precipitation
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: precip_s_input !< Used to read in solid precipitation
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: T_top_input !< Used to read in T_top
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: T_bottom_input !< Used to read in T_bottom
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: S_bu_bottom_input !< Used to read in S_bu_bottom
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: fl_q_bottom_input !< Used to read in fl_q_bottom
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: time_input !< Used to read in time from ERA for atmoflux_flag==2
+ INTEGER :: time_counter !< Keeps track of input data
+ INTEGER :: days !< Days of simulation, specified in init.py
+
+
+ !REAL(wp), DIMENSION(:), ALLOCATABLE :: perm !< permeability values
+
+
+ !##########################################################################################
+ !Chemicals Baby!
+ !All arrays needed to support bigogeochemical tracers
+ !Chemical matrixes: index 1 defines the chemical, index 2 the layer
+ !##########################################################################################
+ INTEGER :: bgc_flag !< 1: no bgc, 2:bgc
+ INTEGER :: N_bgc !< Number of chemicals
+ REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fl_brine_bgc !< Brine fluxes in a matrix, [kg/s], first index is the layer of origin, and the second index is the layer of arrival
+ REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bgc_abs !< Absolute amount of chemicals [kmol] for each tracer
+ REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bgc_bu !< Bulk amounts of chemicals [kmol/kg]
+ REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bgc_br !< Brine concentrations of chems [kmol/kg]
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: bgc_bottom !< Bulk concentrations of chems below the ice [kmol/kg]
+
+
+ !##########################################################################################
+ !Variables needed for tank experiments in which concentrations below the ice change over time
+ !##########################################################################################
+ REAL(wp), DIMENSION(:), ALLOCATABLE :: bgc_total !< Total of chems, for lab experiments with a fixed total amount
+ REAL(wp) :: m_total !< Total initial water mass, for lab experiments with a fixed total amount
+ REAL(wp) :: S_total !< Total initial salt mass, for lab experiments with a fixed total amount
+ REAL(wp) :: tank_depth !< water depth in meters, used to calculate concentrations below ice for tank experiments
+
+ !#######
+ !Misc
+ !#######
+ CHARACTER*3 :: flush_question='No!' !< Niels, 2017 add: used to indicate in stdout wether flushing occurs at this moment or not
+ REAL(wp) :: melt_err=0._wp !< Niels, 2017 add: used to check how much meltwater vanishes in flushing routine
+ 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_grotz.f90 b/mo_grotz.f90
index 008deed85c95d3d2c4f63fe0c9370f5122bf5e38..b0f69ca22ff5b2b3bd26425d26b22d02012edfdc 100644
--- a/mo_grotz.f90
+++ b/mo_grotz.f90
@@ -79,6 +79,7 @@ CONTAINS
!! Basic thermodynamics and layer_dynamics for fixed boundaries seem stable, backup made. by griewank (2010-08-10) \n
!! Add some more outputs, changed routine names and arguments with respect to newly introduces flags by Niels Fuchs, MPIMET (2017-03-01) \n
!! Added a bit of description with the run down of what happends by Philipp Griewank, Uni K (2018-08-08)
+ !! Edited by Niels Fuchs, UHH (2024-05-03), correct layer enthalpy balance in brine expulsion routine
SUBROUTINE grotz ()
USE mo_parameters
@@ -287,7 +288,7 @@ CONTAINS
!##########################################################################################
!Brine flux due to Expulsion
!##########################################################################################
- CALL expulsion_flux (thick, V_ex, Nlayer, N_active, psi_g, fl_m, m)
+ CALL expulsion_flux (thick, V_ex, Nlayer, N_active, psi_g, fl_m, m, H_abs)
IF (i .NE. 1) THEN !is disabled on first timestep to avoid expulsion due to unbalanced initial conditions
!This ensures that the salinity is not effected if the initial mass of the layers is to high
CALL mass_transfer (Nlayer, N_active, T, H_abs, S_abs, S_bu, T_bottom, S_bu_bottom, fl_m)
diff --git a/mo_layer_dynamics.f90 b/mo_layer_dynamics.f90
index 1ed56b7e67981e0321c1443eedff55d678c33f66..27d1c53cebeeb02ad307f7640b003a9cfcf6da03 100644
--- a/mo_layer_dynamics.f90
+++ b/mo_layer_dynamics.f90
@@ -1,718 +1,718 @@
-!>
-!! Mo_layer_dynamics contains all subroutines for the growth and shrinking of layer thickness.
-!!
-!! The middle layers have flexible thickness in contrast to the lower and upper layers which have static thickness.
-!! The details are provided in the separate subroutines.
-!!
-!! @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
-!! Shrinking and growth at the bottom are started by Philipp Griewank, IMPRS (2010-07-28) \n
-!! add melt_thick_output by Niels Fuchs, MPIMET (2017-03-01)
-!!
-MODULE mo_layer_dynamics
-
- USE mo_parameters
- USE mo_output, ONLY:output_raw_lay
-
-
- IMPLICIT NONE
-
- PRIVATE
-
- PUBLIC :: layer_dynamics,top_melt,top_grow
-
-
-
-CONTAINS
-
-
- !>
- !! Organizes the Semi-Adaptive grid SAMSIM uses.
- !!
- !! Modifies the grid and all core variables due to growth or melt.
- !! Calls the different subroutines according to current conditions.
- !! All subroutines can be called with or without biogeochemical tracers active, which is triggered by providing bgc_abs when calling the subroutine.
- !! See Griewank PhD thesis for a full description of the grid.
- !!
- !! Conditions under which following layer dynamics subroutines are called:
- !! - bottom_melt: lowest layer is ice free, second lowest layer has a solid fraction smaller than phi_s_min/2, and all Nlayer layers are active.
- !! - bottom_melt_simple: lowest layer is ice free, second lowest layer has a solid fraction smaller than phi_s_min/2, and not all Nlayer layers are active.
- !! - bottom_melt_simple: lowest layer is ice free, second lowest layer has a solid fraction smaller than phi_s_min/2, all Nlayer layers are active, and the thickness of the middle layers equals thick_0
- !! - bottom_growth_simple: lowest layer has a solid fraction higher than psi_s_min, and not all Nlayer layers are active
- !! - bottom_growth: lowest layer has a solid fraction higher than psi_s_min, and all Nlayer layers are active
- !! - top_grow: top layer thicker than 3/2 * thick_0
- !! - top_melt: top layer thinner than 1/2 * thick_0
- !!
- !! If debug_flag is set to 2 the layer values will be written into the debug output (thermoXX.dat) before and after layer dynamics with a string to identify which subroutine was called
- !!
- !! @par Revision History
- !! created by Philipp Griewank, IMPRS (2010-07-29) \n
- !! first complete and hopefully stable version by Philipp Griewank, IMPRS (2010-08-10)
-
- SUBROUTINE layer_dynamics (phi,N_active,Nlayer,N_bottom,N_middle,N_top,m,S_abs,H_abs,thick,thick_0,T_bottom,S_bu_bottom,&
- & bottom_flag,debug_flag,melt_thick_output,N_bgc,bgc_abs,bgc_bottom)
-
- INTEGER, INTENT(in) :: Nlayer,N_bottom,N_middle,N_top,bottom_flag,debug_flag
- INTEGER, INTENT(inout) :: N_active
- REAL(wp), DIMENSION(Nlayer), INTENT(in) :: phi
- REAL(wp), INTENT(in) :: T_bottom,S_bu_bottom,thick_0
- REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: m,S_abs,H_abs,thick
- INTEGER, INTENT(in) :: N_bgc
- REAL(wp), INTENT(inout) :: melt_thick_output !< Niels, 2017 add: melt_thick_output !OBS: only 3rd element in standard melt_thick_output vector!
- REAL(wp), DIMENSION(Nlayer,N_bgc), INTENT(inout), OPTIONAL :: bgc_abs
- REAL(wp), DIMENSION(N_bgc), INTENT(in), OPTIONAL :: bgc_bottom
-
-
- IF (debug_flag==2) THEN
- CALL output_raw_lay(Nlayer,N_active,H_abs,m,S_abs,thick,'layer_')
- END IF
-
- !_____________________________________________________________________________________________________________
- !___Bottom melt_______________________________________________________________________________________________
- !_____________________________________________________________________________________________________________
- IF (phi(Nlayer-1)<=psi_s_min/2._wp .AND. phi(N_active)<0.00001_wp .AND. N_active==Nlayer &
- & .AND. thick(N_top+1)/thick_0>1.000001_wp .AND. bottom_flag==1) THEN
- IF ( PRESENT(bgc_abs)) THEN
- CALL bottom_melt(Nlayer,N_middle,N_top,m,S_abs,H_abs,thick,N_bgc,bgc_abs)
- ELSE
- CALL bottom_melt(Nlayer,N_middle,N_top,m,S_abs,H_abs,thick,N_bgc)
- END IF
- IF (debug_flag==2) THEN
- CALL output_raw_lay(Nlayer,N_active,H_abs,m,S_abs,thick,'BoMelt')
- END IF
- ELSE IF ( N_active > 1 .AND. N_active<Nlayer .AND. phi(N_active)<0.00001_wp &
- & .AND. phi(MAX(N_active-1,1))<=psi_s_min/2._wp .AND. bottom_flag==1) THEN
- IF ( PRESENT(bgc_abs)) THEN
- CALL bottom_melt_simple(N_active,Nlayer,m,S_abs,H_abs,thick,N_bgc,bgc_abs)
- ELSE
- CALL bottom_melt_simple(N_active,Nlayer,m,S_abs,H_abs,thick,N_bgc)
- END IF
-
- IF (debug_flag==2) THEN
- CALL output_raw_lay(Nlayer,N_active,H_abs,m,S_abs,thick,'BoMeS1')
- END IF
- ELSE IF ( N_active > 1 .AND. phi(N_active)<0.00001_wp .AND. phi(MAX(N_active-1,1))<=psi_s_min/2._wp &
- .AND. (thick(N_top+1)/thick_0)<1.01_wp .AND. bottom_flag==1 ) THEN
- IF ( PRESENT(bgc_abs)) THEN
- CALL bottom_melt_simple(N_active,Nlayer,m,S_abs,H_abs,thick,N_bgc,bgc_abs)
- ELSE
- CALL bottom_melt_simple(N_active,Nlayer,m,S_abs,H_abs,thick,N_bgc)
- END IF
- IF (debug_flag==2) THEN
- CALL output_raw_lay(Nlayer,N_active,H_abs,m,S_abs,thick,'BoMeS2')
- END IF
-
-
- !_____________________________________________________________________________________________________________
- !___Bottom growth_____________________________________________________________________________________________
- !_____________________________________________________________________________________________________________
-
- ELSE IF (phi(N_active)>psi_s_min .AND. N_active<Nlayer .AND. bottom_flag==1) THEN
- IF ( PRESENT(bgc_abs)) THEN
- CALL bottom_growth_simple(N_active,Nlayer,m,S_abs,H_abs,thick,thick_0,T_bottom,S_bu_bottom,N_bgc,bgc_abs,bgc_bottom)
- ELSE
- CALL bottom_growth_simple(N_active,Nlayer,m,S_abs,H_abs,thick,thick_0,T_bottom,S_bu_bottom,N_bgc)
- END IF
- IF (debug_flag==2) THEN
- CALL output_raw_lay(Nlayer,N_active,H_abs,m,S_abs,thick,'BoGrSi')
- END IF
-
- ELSE IF (phi(Nlayer)>psi_s_min .AND. bottom_flag==1) THEN
- IF ( PRESENT(bgc_abs)) THEN
- CALL bottom_growth (Nlayer,N_bottom,N_middle,N_top,m,S_abs,H_abs,thick,T_bottom,S_bu_bottom,N_bgc,bgc_abs,bgc_bottom)
- ELSE
- CALL bottom_growth (Nlayer,N_bottom,N_middle,N_top,m,S_abs,H_abs,thick,T_bottom,S_bu_bottom,N_bgc)
- END IF
- IF (debug_flag==2) THEN
- CALL output_raw_lay(Nlayer,N_active,H_abs,m,S_abs,thick,'BoGrow')
- END IF
-
- !_____________________________________________________________________________________________________________
- !___Top growth________________________________________________________________________________________________
- !_____________________________________________________________________________________________________________
- ELSE IF (thick(1)>1.5_wp*thick_0) THEN
- melt_thick_output = melt_thick_output - thick(1) !< Niels, 2017 add: subtract top growth from melt thick output
- IF ( PRESENT(bgc_abs)) THEN
- CALL top_grow(Nlayer,N_active,N_middle,N_top,thick_0,m,S_abs,H_abs,thick,N_bgc,bgc_abs)
- ELSE
- CALL top_grow(Nlayer,N_active,N_middle,N_top,thick_0,m,S_abs,H_abs,thick,N_bgc)
- END IF
- IF (debug_flag==2) THEN
- CALL output_raw_lay(Nlayer,N_active,H_abs,m,S_abs,thick,'ToGro2')
- END IF
- melt_thick_output = melt_thick_output + thick(1)
-
- !_____________________________________________________________________________________________________________
- !___Top melt__________________________________________________________________________________________________
- !_____________________________________________________________________________________________________________
- ELSE IF (thick(1)<0.5_wp*thick_0) THEN
- melt_thick_output = melt_thick_output - thick(1)
- IF ( PRESENT(bgc_abs)) THEN
- CALL top_melt(Nlayer,N_active,N_middle,N_top,thick_0,m,S_abs,H_abs,thick,N_bgc,bgc_abs)
- ELSE
- CALL top_melt(Nlayer,N_active,N_middle,N_top,thick_0,m,S_abs,H_abs,thick,N_bgc)
- END IF
- IF (debug_flag==2) THEN
- CALL output_raw_lay(Nlayer,N_active,H_abs,m,S_abs,thick,'ToMel2')
- END IF
- melt_thick_output = melt_thick_output + thick(1)
-
-
- END IF
-
- END SUBROUTINE layer_dynamics
-
-
- !!>
- !! Controls melting of the surface layer.
- !!
- !! Subroutine is called when thick(1)<0.5*thick_0(1).
- !! If N_active<Nlayer then N_active = N_active-1.
- !! If N_active=Nlayer and middle layer thickness = thick_0 than N_active = N_active-1.
- !! Otherwise if N_active=Nlayer the middle layers are shrunk.
- !!
- !! @par Revision History
- !! Pasted by Philipp Griewank, IMPRS (2011-05-10)
- !!
-
-
- SUBROUTINE top_melt (Nlayer,N_active,N_middle,N_top,thick_0,m,S_abs,H_abs,thick,N_bgc,bgc_abs)
-
- INTEGER, INTENT(in) :: Nlayer,N_middle,N_top
- INTEGER, INTENT(inout) :: N_active
- REAL(wp), INTENT(in) :: thick_0
- REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: m,S_abs,H_abs,thick
- INTEGER, INTENT(in) :: N_bgc
- REAL(wp), DIMENSION(Nlayer,N_bgc), INTENT(inout), OPTIONAL :: bgc_abs
- INTEGER :: k
- REAL(wp) :: loss_m,loss_S_abs,loss_H_abs
- REAL(wp) :: shift !> distance the border between layer moves
- REAL(wp), DIMENSION(Nlayer) :: rho !> density of layer
- REAL(wp), DIMENSION(Nlayer) :: H,S_bu
- REAL(wp), DIMENSION(N_bgc) :: loss_bgc
- REAL(wp), DIMENSION(Nlayer,N_bgc) :: bgc_temp, bgc_bulk
-
-
- IF (PRESENT(bgc_abs)) THEN
- bgc_temp = bgc_abs
- ELSE
- bgc_temp = 0._wp
- END IF
-
- !Calculating rho, H and S_bu
- DO k = 1,N_active
- rho(k) = m(k)/thick(k)
- S_bu(k) = S_abs(k)/m(k)
- H(k) = H_abs(k)/m(k)
- bgc_bulk(k,:) = bgc_temp(k,:)/m(k)
- END DO
-
- !Adjusting the top layer from thick(1)<0.5*thick_0(1) to thick(1)=thick(1)+thick_0(1)
- loss_m = thick_0 *rho(1)
- loss_S_abs = loss_m *S_bu(1)
- loss_H_abs = loss_m *H(1)
- loss_bgc(:) = loss_m *bgc_bulk(1,:)
-
- m(1) = m(1) +m(2)
- S_abs(1) = S_abs(1) +S_abs(2)
- H_abs(1) = H_abs(1) +H_abs(2)
- bgc_temp(1,:) = bgc_temp(1,:) +bgc_temp(2,:)
- thick(1) = thick(1) +thick(2)
-
- !Adjusting top layers
- DO k = 2,MIN(N_top-1,N_active-1)
- m(k) = rho (k+1) *thick_0
- S_abs(k) = S_bu(k+1) *rho(k+1)*thick_0
- H_abs(k) = H (k+1) *rho(k+1)*thick_0
- bgc_temp(k,:) = bgc_bulk (k+1,:) *rho(k+1)*thick_0
- END DO
-
-
- !Removing bottom layer if N_active<N_top
- IF (N_active<=N_top) THEN
- m(N_active) = 0.0_wp
- S_abs(N_active) = 0.0_wp
- H_abs(N_active) = 0.0_wp
- thick(N_active) = 0.0_wp
- bgc_temp(N_active,:) = 0.0_wp
-
- N_active = N_active-1
-
- ELSE IF (N_active>N_top .AND. N_active<=Nlayer .AND. thick(N_top+1)/thick_0<1.00001_wp) THEN
- !Removing bottom layer and moving all other layers if N_active<Nlayer .and. N_active>N_top.
- !Is also activated if thick(N_top+1)=thick_0
- DO k = N_top,N_active-1
- m(k) = rho (k+1) *thick_0
- S_abs(k) = S_bu(k+1) *rho(k+1)*thick_0
- H_abs(k) = H (k+1) *rho(k+1)*thick_0
- bgc_temp(k,:) = bgc_bulk (k+1,:) *rho(k+1)*thick_0
- END DO
- m(N_active) = 0.0_wp
- S_abs(N_active) = 0.0_wp
- H_abs(N_active) = 0.0_wp
- bgc_temp(N_active,:) = 0.0_wp
- thick(N_active) = 0.0_wp
-
-
- N_active = N_active-1
- END IF
-
- IF (N_active==Nlayer .AND. thick(N_top+1)-thick_0>=0.000001_wp ) THEN
- !Middle layers are adjusted__________________________
- loss_m = thick_0*rho(N_top+1)
- loss_S_abs = loss_m*S_bu(N_top+1)
- loss_H_abs = loss_m*H(N_top+1)
- loss_bgc = loss_m*bgc_bulk(N_top+1,:)
-
- !layer N_top needs new values________________________
- m(N_top) = loss_m
- S_abs(N_top) = loss_S_abs
- H_abs(N_top) = loss_H_abs
- bgc_temp(N_top,:) = loss_bgc(:)
-
-
- DO k=N_top+1,N_middle+N_top
- m(k) = m(k) -loss_m
- H_abs(k) = H_abs(k) -loss_H_abs
- S_abs(k) = S_abs(k) -loss_S_abs
- bgc_temp(k,:) = bgc_temp(k,:) -loss_bgc(:)
-
- shift = thick_0*REAL(N_middle-k+N_top)/REAL(N_middle) !distance the border moves
-
- loss_m = shift *rho(k+1)
- loss_S_abs = loss_m *S_bu(k+1)
- loss_H_abs = loss_m *H(k+1)
- loss_bgc = loss_m *bgc_bulk(k+1,:)
-
- m(k) = m(k) +loss_m
- H_abs(k) = H_abs(k) +loss_H_abs
- S_abs(k) = S_abs(k) +loss_S_abs
- bgc_temp(k,:) = bgc_temp(k,:) +loss_bgc(:)
- END DO
-
-
-
-
-
- !Middle layer thickness is adjusted________________________
- DO k = N_top+1,N_top+N_middle
- thick(k) = thick(k)-thick_0/REAL(N_middle)
- END DO
- END IF
-
-
- !Checks if current thick vector is possible for current N_active and thick_0
- !Uses 0.501 instead of 0.500 to give a bit of wriggle room for numerical noise
- IF(thick_0*(N_active+0.501_wp)<=SUM(thick) .AND. N_active<Nlayer) THEN
- PRINT*,'wtf layer issue',thick_0,SUM(thick),N_active,Nlayer
- STOP 7889
- END IF
-
- IF (PRESENT(bgc_abs)) THEN
- bgc_abs = bgc_temp
- END IF
- END SUBROUTINE top_melt
-
-
- !>
- !! Controls melting at the bottom.
- !!
- !! Subroutine is called if the second lowest layer is fully liquid and the middle layers are thicker then the bottom layers.
- !! When this occurs the thickness of the middle layers is reduced by 1/N_middle*thick_0.
- !! The lower layers are adjusted upwards, layer(k)=layer(k-1).
- !!
- !!
- !! @par Revision History
- !! Started by Philipp Griewank, IMPRS (2010-07-28) \n
- !! Linear profiles replaced with simple upstream by Philipp Griewank, IMPRS (2010-08-25)
- !!
- SUBROUTINE bottom_melt (Nlayer,N_middle,N_top,m,S_abs,H_abs,thick,N_bgc,bgc_abs)
-
- INTEGER, INTENT(in) :: Nlayer,N_middle,N_top
- REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: m,S_abs,H_abs,thick
- INTEGER :: k
- REAL(wp) :: loss_m,loss_S_abs,loss_H_abs,shift
- REAL(wp), DIMENSION(Nlayer) :: rho !>density of layer
- REAL(wp), DIMENSION(Nlayer) :: H,S_bu
- INTEGER, INTENT(in) :: N_bgc
- REAL(wp), DIMENSION(Nlayer,N_bgc), INTENT(inout), OPTIONAL :: bgc_abs
- REAL(wp), DIMENSION(N_bgc) :: loss_bgc
- REAL(wp), DIMENSION(Nlayer,N_bgc) :: bgc_temp, bgc_bulk
-
-
- IF (PRESENT(bgc_abs)) THEN
- bgc_temp = bgc_abs
- ELSE
- bgc_temp = 0._wp
- END IF
-
-
- !Calculating rho, H and S_bu
- DO k = N_top+1,Nlayer
- rho(k) = m(k)/thick(k)
- S_bu(k) = S_abs(k)/m(k)
- H(k) = H_abs(k)/m(k)
- bgc_bulk(k,:) = bgc_temp(k,:)/m(k)
-
- END DO
-
- loss_m = 0._wp
- loss_S_abs = 0._wp
- loss_H_abs = 0._wp
- loss_bgc = 0._wp
-
- !_________________________Middle layers are adjusted__________________________
- DO k = N_top+1,N_top+N_middle
-
- !each layer receives all that the layer above loses
- m(k) = m(k) +loss_m
- H_abs(k) = H_abs(k) +loss_H_abs
- S_abs(k) = S_abs(k) +loss_S_abs
- bgc_temp(k,:) = bgc_temp(k,:) +loss_bgc(:)
-
-
-
- shift = thick(Nlayer)*(k-N_top)/ REAL(N_middle) !distance the border moves
-
- loss_m = shift*rho(k)
- loss_H_abs = loss_m*H(k)
- loss_S_abs = loss_m*S_bu(k)
- loss_bgc(:) = loss_m*bgc_bulk(k,:)
-
-
- m(k) = m(k) -loss_m
- H_abs(k) = H_abs(k) -loss_H_abs
- S_abs(k) = S_abs(k) -loss_S_abs
- bgc_temp(k,:) = bgc_temp(k,:) -loss_bgc(:)
-
-
- END DO
-
-
- !________________________Middle layer thickness is adjusted
-
- DO k = N_top+1,N_top+N_middle
- thick(k) = thick(k)-thick(Nlayer)/REAL(N_middle)
- END DO
- !_________________________Bottom layers are adjusted__________________________
-
- DO k = N_top+N_middle+1,Nlayer
- H_abs(k) = rho(k-1)*thick(k)*H(k-1)
- S_abs(k) = rho(k-1)*thick(k)*S_bu(k-1)
- m(k) = rho(k-1)*thick(k)
- bgc_temp(k,:) = rho(k-1)*thick(k)*bgc_bulk(k-1,:)
- END DO
-
- IF (PRESENT(bgc_abs)) THEN
- bgc_abs = bgc_temp
- END IF
-
- END SUBROUTINE bottom_melt
-
-
-
- !>
- !! Controls bottom growth.
- !!
- !! Subroutine is called if the lowest layer is not fully liquid and the middle layers are thicker then thick_0.
- !! When this occurs the thickness of the middle layers is increased by 1/N_middle*thick(bottom).
- !! The lower layers are adjusted downwards, layer(k)=layer(k+1).
- !!
- !!
- !! @par Revision History
- !! Started by Philipp Griewank, IMPRS (2010-07-29) \n
- !! Linear profiles removed, replaced with simple upwind by Philipp Griewank, IMPRS (2010-08-25)
- !!
- SUBROUTINE bottom_growth (Nlayer,N_bottom,N_middle,N_top,m,S_abs,H_abs,thick,T_bottom,S_bu_bottom,N_bgc,bgc_abs,bgc_bottom)
-
- INTEGER, INTENT(in) :: Nlayer,N_bottom,N_middle,N_top
- REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: m,S_abs,H_abs,thick
- REAL(wp), INTENT(in) :: T_bottom,S_bu_bottom
- INTEGER :: k
- REAL(wp) :: gain_m,gain_S_abs,gain_H_abs,shift
- REAL(wp), DIMENSION(Nlayer) :: rho !>density of layer
- REAL(wp), DIMENSION(Nlayer) :: H,S_bu
- INTEGER, INTENT(in) :: N_bgc
- REAL(wp), DIMENSION(Nlayer,N_bgc), INTENT(inout), OPTIONAL :: bgc_abs
- REAL(wp), DIMENSION(N_bgc), INTENT(in), OPTIONAL :: bgc_bottom
- REAL(wp), DIMENSION(N_bgc) :: gain_bgc
- REAL(wp), DIMENSION(Nlayer,N_bgc) :: bgc_temp, bgc_bulk
-
-
- IF (PRESENT(bgc_abs)) THEN
- bgc_temp = bgc_abs
- ELSE
- bgc_temp = 0._wp
- END IF
-
-
- !Calculating rho, H and S_bu
- DO k = N_top+1,N_top+N_middle+1
- rho(k) = m(k)/thick(k)
- S_bu(k) = S_abs(k)/m(k)
- H(k) = H_abs(k)/m(k)
- bgc_bulk(k,:) = bgc_temp(k,:)/m(k)
- END DO
-
- gain_m = 0._wp
- gain_S_abs = 0._wp
- gain_H_abs = 0._wp
- gain_bgc = 0._wp
-
- !_________________________Middle layers are adjusted__________________________
- DO k = N_top+1,N_top+N_middle
-
- !each layer loses all that the upper layer gained
- m(k) = m(k) -gain_m
- H_abs(k) = H_abs(k) -gain_H_abs
- S_abs(k) = S_abs(k) -gain_S_abs
- bgc_temp(k,:) = bgc_temp(k,:) -gain_bgc(:)
-
- shift = thick(Nlayer)*(k-N_top)/ REAL(N_middle) !distance the border moves
-
- gain_m = shift*rho(k+1)
- gain_H_abs = gain_m*H(k+1)
- gain_S_abs = gain_m*S_bu(k+1)
- gain_bgc(:)= gain_m*bgc_bulk(k+1,:)
-
- m(k) = m(k) +gain_m
- H_abs(k) = H_abs(k) +gain_H_abs
- S_abs(k) = S_abs(k) +gain_S_abs
- bgc_temp(k,:) = bgc_temp(k,:) +gain_bgc(:)
-
- END DO
-
- !________________________Middle layer thickness is adjusted
- DO k=N_top+1,N_top+N_middle
- thick(k) = thick(k)+thick(Nlayer)/REAL(N_middle)
- END DO
-
- !_________________________Bottom layers are adjusted__________________________
- DO k = Nlayer-N_bottom+1,Nlayer-1
- H_abs(k) = H_abs(k+1)
- S_abs(k) = S_abs(k+1)
- m(k) = m(k+1)
- bgc_temp(k,:) = bgc_temp(k+1,:)
- END DO
-
- !_________________________Lowest layer is set using T_bottom and S_bu_bottom
- m(Nlayer) = thick(Nlayer)*rho_l
- H_abs(Nlayer) = m(Nlayer)*T_bottom*c_l
- S_abs(Nlayer) = m(Nlayer)*S_bu_bottom
-
- IF (PRESENT(bgc_abs)) THEN
- bgc_temp(Nlayer,:) = m(Nlayer)*bgc_bottom(:)
- bgc_abs = bgc_temp
- END IF
-
- END SUBROUTINE bottom_growth
-
-
-
-
- !>
- !! Controls growth if the middle layers are not all active.
- !!
- !! Subroutine is called if the lowest layer is not fully liquid and the middle layers are not all activated.
- !! A new layer is activated.
- !! The new lowest layer is assigned standard values of T_bottom and S_bu_bottom.
- !!
- !!
- !! @par Revision History
- !! Started by Philipp Griewank, IMPRS (2010-07-29) \n
- !! Expanded to deal with not full top layers by Philipp Griewank, IMPRS (2011-01-13) \n
- !! Expansion removed and bgc added by Philipp Griewank, IMPRS (2014-02-04) \n
- SUBROUTINE bottom_growth_simple (N_active,Nlayer,m,S_abs,H_abs,thick,thick_0,T_bottom,S_bu_bottom,N_bgc,bgc_abs,bgc_bottom)
-
- INTEGER, INTENT(in) :: Nlayer
- INTEGER, INTENT(inout) :: N_active
- REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: m,S_abs,H_abs,thick
- REAL(wp), INTENT(in) :: thick_0
- REAL(wp), INTENT(in) :: T_bottom,S_bu_bottom
- INTEGER, INTENT(in) :: N_bgc
- REAL(wp), DIMENSION(Nlayer,N_bgc), INTENT(inout), OPTIONAL :: bgc_abs
- REAL(wp), DIMENSION(N_bgc), INTENT(in), OPTIONAL :: bgc_bottom
-
-
- N_active = N_active+1
- thick(N_active) = thick_0
-
- m(N_active) = thick(N_active)*rho_l
- H_abs(N_active) = m(N_active) *T_bottom*c_l
- S_abs(N_active) = m(N_active) *S_bu_bottom
-
- IF (PRESENT(bgc_abs)) THEN
- bgc_abs(N_active,:) = bgc_bottom(:)*m(N_active)
- END IF
-
- END SUBROUTINE bottom_growth_simple
-
- !>
- !! Controls bottom melting if the middle layers are not all active.
- !!
- !! Subroutine is called if the second lowest active layer is fully liquid and the middle layers are not all activated, or if the middle layers are thick_0 thick.
- !! The lowest layer is simply deactivated.
- !!
- !!
- !! @par Revision History
- !! Started by Philipp Griewank, IMPRS (2010-07-30>)
- !!
- SUBROUTINE bottom_melt_simple(N_active,Nlayer,m,S_abs,H_abs,thick,N_bgc,bgc_abs)
-
- INTEGER, INTENT(in) :: Nlayer
- INTEGER, INTENT(inout) :: N_active
- REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: m,S_abs,H_abs,thick
- INTEGER, INTENT(in) :: N_bgc
- REAL(wp), DIMENSION(Nlayer,N_bgc), INTENT(inout), OPTIONAL :: bgc_abs
-
- thick(N_active) = 0._wp
- m(N_active) = 0._wp
- S_abs(N_active) = 0._wp
- H_abs(N_active) = 0._wp
- IF(PRESENT(bgc_abs)) THEN
- bgc_abs(N_active,:) = 0._wp
- END IF
- N_active = N_active-1
-
- END SUBROUTINE bottom_melt_simple
-
-
-
- !>
- !! Top grow subroutine.
- !!
- !! Should be called when the top layer is thicker then 1.5 *thick_0.
- !! If N_active=Nlayer middle layers are expanded by thick_0/N_middle and top layers are moved one down.
- !! IF N_active<Nlayer then N_active=N_active+1 and all layers are shifted downwards.
- !!
- !!
- !! @par Revision History
- !! Started by Philipp Griewank, IMPRS (2011-05-10>)
- !!
-
-
- SUBROUTINE top_grow (Nlayer,N_active,N_middle,N_top,thick_0,m,S_abs,H_abs,thick,N_bgc,bgc_abs)
-
- INTEGER, INTENT(in) :: Nlayer,N_middle,N_top
- INTEGER, INTENT(inout) :: N_active
- REAL(wp), INTENT(in) :: thick_0
- REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: m,S_abs,H_abs,thick
- INTEGER :: k
- REAL(wp) :: loss_m,loss_S_abs,loss_H_abs,shift
- REAL(wp), DIMENSION(Nlayer) :: rho !>density of layer
- REAL(wp), DIMENSION(Nlayer) :: H,S_bu
- INTEGER, INTENT(in) :: N_bgc
- REAL(wp), DIMENSION(Nlayer,N_bgc), INTENT(inout), OPTIONAL :: bgc_abs
- REAL(wp), DIMENSION(N_bgc) :: loss_bgc
- REAL(wp), DIMENSION(Nlayer,N_bgc) :: bgc_temp, bgc_bulk
-
- IF (PRESENT(bgc_abs)) THEN
- bgc_temp = bgc_abs
- ELSE
- bgc_temp = 0._wp
- END IF
- !Calculating rho, H and S_bu
- DO k = 1,N_active
- rho(k) = m(k)/thick(k)
- S_bu(k) = S_abs(k)/m(k)
- H(k) = H_abs(k)/m(k)
- bgc_bulk(k,:) = bgc_temp(k,:)/m(k)
- END DO
-
- !Adjusting the top layer from thick(1)>1.5*thick_0(1) to thick(1)=thick(1)-thick_0(1)
- loss_m = thick_0*rho(1)
- loss_S_abs = loss_m*S_bu(1)
- loss_H_abs = loss_m*H(1)
- loss_bgc(:) = loss_m*bgc_bulk(1,:)
-
- m(1) = m(1) -loss_m
- S_abs(1) = S_abs(1) -loss_S_abs
- H_abs(1) = H_abs(1) -loss_H_abs
- bgc_temp(1,:) = bgc_temp(1,:) -loss_bgc(:)
- thick(1) = thick(1) -thick_0
-
- !Adjusting top layers
- DO k = 2,MIN(N_top,N_active)
- m(k) = rho (k-1) *thick_0
- S_abs(k) = S_bu(k-1) *rho(k-1)*thick_0
- H_abs(k) = H (k-1) *rho(k-1)*thick_0
- bgc_temp(k,:) = bgc_bulk(k-1,:) *rho(k-1)*thick_0
- END DO
-
- !Adding new bottom layer if N_active<=N_top+1
- IF (N_active<=N_top) THEN
- N_active = N_active+1
- m(N_active) = rho (N_active-1) *thick_0
- S_abs(N_active) = S_bu(N_active-1) *thick_0*rho(N_active-1)
- H_abs(N_active) = H (N_active-1) *thick_0*rho(N_active-1)
- bgc_temp(N_active,:) = bgc_bulk(N_active-1,:) *thick_0*rho(N_active-1)
- thick(N_active) = thick_0
-
- !Adding new bottom layer and moving all other layers if N_active<Nlayer .and. N_active>N_top
- ELSE IF (N_active>N_top .AND. N_active<Nlayer) THEN
- DO k = N_top+1,N_active
- m(k) = rho (k-1) *thick_0
- S_abs(k) = S_bu(k-1) *rho(k-1)*thick_0
- H_abs(k) = H (k-1) *rho(k-1)*thick_0
- bgc_temp(k,:) = bgc_bulk(k-1,:) *rho(k-1)*thick_0
- END DO
- N_active = N_active+1
- m(N_active) = rho (N_active-1) *thick_0
- S_abs(N_active) = S_bu(N_active-1) *thick_0 *rho(N_active-1)
- H_abs(N_active) = H (N_active-1) *thick_0 *rho(N_active-1)
- bgc_temp(N_active,:) = bgc_bulk(N_active-1,:) *thick_0 *rho(N_active-1)
- thick(N_active) = thick_0
-
- ELSE IF (N_active==Nlayer) THEN
- !Middle layers are adjusted__________________________
- loss_m = thick_0*rho(N_top)
- loss_S_abs = loss_m*S_bu(N_top)
- loss_H_abs = loss_m*H(N_top)
- loss_bgc(:) = loss_m*bgc_bulk(N_top,:)
-
- DO k = N_top+1,N_middle+N_top
- m(k) = m(k) +loss_m
- H_abs(k) = H_abs(k) +loss_H_abs
- S_abs(k) = S_abs(k) +loss_S_abs
- bgc_temp(k,:) = bgc_temp(k,:) +loss_bgc(:)
-
- shift = thick_0*REAL(N_middle-k+N_top)/REAL(N_middle) !distance the border moves
-
- loss_m = shift *rho(k)
- loss_S_abs = loss_m *S_bu(k)
- loss_H_abs = loss_m *H(k)
- loss_bgc(:) = loss_m *bgc_bulk(k,:)
-
- m(k) = m(k) -loss_m
- H_abs(k) = H_abs(k) -loss_H_abs
- S_abs(k) = S_abs(k) -loss_S_abs
- bgc_temp(k,:) = bgc_temp(k,:) -loss_bgc(:)
- END DO
-
- !Middle layer thickness is adjusted
- DO k = N_top+1,N_top+N_middle
- thick(k) = thick(k)+thick_0/REAL(N_middle)
- END DO
- END IF
-
-
- IF (PRESENT(bgc_abs)) THEN
- bgc_abs = bgc_temp
- END IF
-
- END SUBROUTINE top_grow
-
-
-END MODULE mo_layer_dynamics
+!>
+!! Mo_layer_dynamics contains all subroutines for the growth and shrinking of layer thickness.
+!!
+!! The middle layers have flexible thickness in contrast to the lower and upper layers which have static thickness.
+!! The details are provided in the separate subroutines.
+!!
+!! @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
+!! Shrinking and growth at the bottom are started by Philipp Griewank, IMPRS (2010-07-28) \n
+!! add melt_thick_output by Niels Fuchs, MPIMET (2017-03-01)
+!!
+MODULE mo_layer_dynamics
+
+ USE mo_parameters
+ USE mo_output, ONLY:output_raw_lay
+
+
+ IMPLICIT NONE
+
+ PRIVATE
+
+ PUBLIC :: layer_dynamics,top_melt,top_grow
+
+
+
+CONTAINS
+
+
+ !>
+ !! Organizes the Semi-Adaptive grid SAMSIM uses.
+ !!
+ !! Modifies the grid and all core variables due to growth or melt.
+ !! Calls the different subroutines according to current conditions.
+ !! All subroutines can be called with or without biogeochemical tracers active, which is triggered by providing bgc_abs when calling the subroutine.
+ !! See Griewank PhD thesis for a full description of the grid.
+ !!
+ !! Conditions under which following layer dynamics subroutines are called:
+ !! - bottom_melt: lowest layer is ice free, second lowest layer has a solid fraction smaller than phi_s_min/2, and all Nlayer layers are active.
+ !! - bottom_melt_simple: lowest layer is ice free, second lowest layer has a solid fraction smaller than phi_s_min/2, and not all Nlayer layers are active.
+ !! - bottom_melt_simple: lowest layer is ice free, second lowest layer has a solid fraction smaller than phi_s_min/2, all Nlayer layers are active, and the thickness of the middle layers equals thick_0
+ !! - bottom_growth_simple: lowest layer has a solid fraction higher than psi_s_min, and not all Nlayer layers are active
+ !! - bottom_growth: lowest layer has a solid fraction higher than psi_s_min, and all Nlayer layers are active
+ !! - top_grow: top layer thicker than 3/2 * thick_0
+ !! - top_melt: top layer thinner than 1/2 * thick_0
+ !!
+ !! If debug_flag is set to 2 the layer values will be written into the debug output (thermoXX.dat) before and after layer dynamics with a string to identify which subroutine was called
+ !!
+ !! @par Revision History
+ !! created by Philipp Griewank, IMPRS (2010-07-29) \n
+ !! first complete and hopefully stable version by Philipp Griewank, IMPRS (2010-08-10)
+
+ SUBROUTINE layer_dynamics (phi,N_active,Nlayer,N_bottom,N_middle,N_top,m,S_abs,H_abs,thick,thick_0,T_bottom,S_bu_bottom,&
+ & bottom_flag,debug_flag,melt_thick_output,N_bgc,bgc_abs,bgc_bottom)
+
+ INTEGER, INTENT(in) :: Nlayer,N_bottom,N_middle,N_top,bottom_flag,debug_flag
+ INTEGER, INTENT(inout) :: N_active
+ REAL(wp), DIMENSION(Nlayer), INTENT(in) :: phi
+ REAL(wp), INTENT(in) :: T_bottom,S_bu_bottom,thick_0
+ REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: m,S_abs,H_abs,thick
+ INTEGER, INTENT(in) :: N_bgc
+ REAL(wp), INTENT(inout) :: melt_thick_output !< Niels, 2017 add: melt_thick_output !OBS: only 3rd element in standard melt_thick_output vector!
+ REAL(wp), DIMENSION(Nlayer,N_bgc), INTENT(inout), OPTIONAL :: bgc_abs
+ REAL(wp), DIMENSION(N_bgc), INTENT(in), OPTIONAL :: bgc_bottom
+
+
+ IF (debug_flag==2) THEN
+ CALL output_raw_lay(Nlayer,N_active,H_abs,m,S_abs,thick,'layer_')
+ END IF
+
+ !_____________________________________________________________________________________________________________
+ !___Bottom melt_______________________________________________________________________________________________
+ !_____________________________________________________________________________________________________________
+ IF (phi(Nlayer-1)<=psi_s_min/2._wp .AND. phi(N_active)<0.00001_wp .AND. N_active==Nlayer &
+ & .AND. thick(N_top+1)/thick_0>1.000001_wp .AND. bottom_flag==1) THEN
+ IF ( PRESENT(bgc_abs)) THEN
+ CALL bottom_melt(Nlayer,N_middle,N_top,m,S_abs,H_abs,thick,N_bgc,bgc_abs)
+ ELSE
+ CALL bottom_melt(Nlayer,N_middle,N_top,m,S_abs,H_abs,thick,N_bgc)
+ END IF
+ IF (debug_flag==2) THEN
+ CALL output_raw_lay(Nlayer,N_active,H_abs,m,S_abs,thick,'BoMelt')
+ END IF
+ ELSE IF ( N_active > 1 .AND. N_active<Nlayer .AND. phi(N_active)<0.00001_wp &
+ & .AND. phi(MAX(N_active-1,1))<=psi_s_min/2._wp .AND. bottom_flag==1) THEN
+ IF ( PRESENT(bgc_abs)) THEN
+ CALL bottom_melt_simple(N_active,Nlayer,m,S_abs,H_abs,thick,N_bgc,bgc_abs)
+ ELSE
+ CALL bottom_melt_simple(N_active,Nlayer,m,S_abs,H_abs,thick,N_bgc)
+ END IF
+
+ IF (debug_flag==2) THEN
+ CALL output_raw_lay(Nlayer,N_active,H_abs,m,S_abs,thick,'BoMeS1')
+ END IF
+ ELSE IF ( N_active > 1 .AND. phi(N_active)<0.00001_wp .AND. phi(MAX(N_active-1,1))<=psi_s_min/2._wp &
+ .AND. (thick(N_top+1)/thick_0)<1.01_wp .AND. bottom_flag==1 ) THEN
+ IF ( PRESENT(bgc_abs)) THEN
+ CALL bottom_melt_simple(N_active,Nlayer,m,S_abs,H_abs,thick,N_bgc,bgc_abs)
+ ELSE
+ CALL bottom_melt_simple(N_active,Nlayer,m,S_abs,H_abs,thick,N_bgc)
+ END IF
+ IF (debug_flag==2) THEN
+ CALL output_raw_lay(Nlayer,N_active,H_abs,m,S_abs,thick,'BoMeS2')
+ END IF
+
+
+ !_____________________________________________________________________________________________________________
+ !___Bottom growth_____________________________________________________________________________________________
+ !_____________________________________________________________________________________________________________
+
+ ELSE IF (phi(N_active)>psi_s_min .AND. N_active<Nlayer .AND. bottom_flag==1) THEN
+ IF ( PRESENT(bgc_abs)) THEN
+ CALL bottom_growth_simple(N_active,Nlayer,m,S_abs,H_abs,thick,thick_0,T_bottom,S_bu_bottom,N_bgc,bgc_abs,bgc_bottom)
+ ELSE
+ CALL bottom_growth_simple(N_active,Nlayer,m,S_abs,H_abs,thick,thick_0,T_bottom,S_bu_bottom,N_bgc)
+ END IF
+ IF (debug_flag==2) THEN
+ CALL output_raw_lay(Nlayer,N_active,H_abs,m,S_abs,thick,'BoGrSi')
+ END IF
+
+ ELSE IF (phi(Nlayer)>psi_s_min .AND. bottom_flag==1) THEN
+ IF ( PRESENT(bgc_abs)) THEN
+ CALL bottom_growth (Nlayer,N_bottom,N_middle,N_top,m,S_abs,H_abs,thick,T_bottom,S_bu_bottom,N_bgc,bgc_abs,bgc_bottom)
+ ELSE
+ CALL bottom_growth (Nlayer,N_bottom,N_middle,N_top,m,S_abs,H_abs,thick,T_bottom,S_bu_bottom,N_bgc)
+ END IF
+ IF (debug_flag==2) THEN
+ CALL output_raw_lay(Nlayer,N_active,H_abs,m,S_abs,thick,'BoGrow')
+ END IF
+
+ !_____________________________________________________________________________________________________________
+ !___Top growth________________________________________________________________________________________________
+ !_____________________________________________________________________________________________________________
+ ELSE IF (thick(1)>1.5_wp*thick_0) THEN
+ melt_thick_output = melt_thick_output - thick(1) !< Niels, 2017 add: subtract top growth from melt thick output
+ IF ( PRESENT(bgc_abs)) THEN
+ CALL top_grow(Nlayer,N_active,N_middle,N_top,thick_0,m,S_abs,H_abs,thick,N_bgc,bgc_abs)
+ ELSE
+ CALL top_grow(Nlayer,N_active,N_middle,N_top,thick_0,m,S_abs,H_abs,thick,N_bgc)
+ END IF
+ IF (debug_flag==2) THEN
+ CALL output_raw_lay(Nlayer,N_active,H_abs,m,S_abs,thick,'ToGro2')
+ END IF
+ melt_thick_output = melt_thick_output + thick(1)
+
+ !_____________________________________________________________________________________________________________
+ !___Top melt__________________________________________________________________________________________________
+ !_____________________________________________________________________________________________________________
+ ELSE IF (thick(1)<0.5_wp*thick_0) THEN
+ melt_thick_output = melt_thick_output - thick(1)
+ IF ( PRESENT(bgc_abs)) THEN
+ CALL top_melt(Nlayer,N_active,N_middle,N_top,thick_0,m,S_abs,H_abs,thick,N_bgc,bgc_abs)
+ ELSE
+ CALL top_melt(Nlayer,N_active,N_middle,N_top,thick_0,m,S_abs,H_abs,thick,N_bgc)
+ END IF
+ IF (debug_flag==2) THEN
+ CALL output_raw_lay(Nlayer,N_active,H_abs,m,S_abs,thick,'ToMel2')
+ END IF
+ melt_thick_output = melt_thick_output + thick(1)
+
+
+ END IF
+
+ END SUBROUTINE layer_dynamics
+
+
+ !!>
+ !! Controls melting of the surface layer.
+ !!
+ !! Subroutine is called when thick(1)<0.5*thick_0(1).
+ !! If N_active<Nlayer then N_active = N_active-1.
+ !! If N_active=Nlayer and middle layer thickness = thick_0 than N_active = N_active-1.
+ !! Otherwise if N_active=Nlayer the middle layers are shrunk.
+ !!
+ !! @par Revision History
+ !! Pasted by Philipp Griewank, IMPRS (2011-05-10)
+ !! Edited by Niels Fuchs, UHH (2024-05-03), changed thick_0 to thick(k) in middle layers to work better with ice core initialization, since middle layer thickness can vary
+
+
+ SUBROUTINE top_melt (Nlayer,N_active,N_middle,N_top,thick_0,m,S_abs,H_abs,thick,N_bgc,bgc_abs)
+
+ INTEGER, INTENT(in) :: Nlayer,N_middle,N_top
+ INTEGER, INTENT(inout) :: N_active
+ REAL(wp), INTENT(in) :: thick_0
+ REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: m,S_abs,H_abs,thick
+ INTEGER, INTENT(in) :: N_bgc
+ REAL(wp), DIMENSION(Nlayer,N_bgc), INTENT(inout), OPTIONAL :: bgc_abs
+ INTEGER :: k
+ REAL(wp) :: loss_m,loss_S_abs,loss_H_abs
+ REAL(wp) :: shift !> distance the border between layer moves
+ REAL(wp), DIMENSION(Nlayer) :: rho !> density of layer
+ REAL(wp), DIMENSION(Nlayer) :: H,S_bu
+ REAL(wp), DIMENSION(N_bgc) :: loss_bgc
+ REAL(wp), DIMENSION(Nlayer,N_bgc) :: bgc_temp, bgc_bulk
+
+
+ IF (PRESENT(bgc_abs)) THEN
+ bgc_temp = bgc_abs
+ ELSE
+ bgc_temp = 0._wp
+ END IF
+
+ !Calculating rho, H and S_bu
+ DO k = 1,N_active
+ rho(k) = m(k)/thick(k)
+ S_bu(k) = S_abs(k)/m(k)
+ H(k) = H_abs(k)/m(k)
+ bgc_bulk(k,:) = bgc_temp(k,:)/m(k)
+ END DO
+
+ !Adjusting the top layer from thick(1)<0.5*thick_0(1) to thick(1)=thick(1)+thick_0(1)
+ loss_m = thick_0 *rho(1)
+ loss_S_abs = loss_m *S_bu(1)
+ loss_H_abs = loss_m *H(1)
+ loss_bgc(:) = loss_m *bgc_bulk(1,:)
+
+ m(1) = m(1) +m(2)
+ S_abs(1) = S_abs(1) +S_abs(2)
+ H_abs(1) = H_abs(1) +H_abs(2)
+ bgc_temp(1,:) = bgc_temp(1,:) +bgc_temp(2,:)
+ thick(1) = thick(1) +thick(2)
+
+ !Adjusting top layers
+ DO k = 2,MIN(N_top-1,N_active-1)
+ m(k) = rho (k+1) *thick_0
+ S_abs(k) = S_bu(k+1) *rho(k+1)*thick_0
+ H_abs(k) = H (k+1) *rho(k+1)*thick_0
+ bgc_temp(k,:) = bgc_bulk (k+1,:) *rho(k+1)*thick_0
+ END DO
+
+
+ !Removing bottom layer if N_active<N_top
+ IF (N_active<=N_top) THEN
+ m(N_active) = 0.0_wp
+ S_abs(N_active) = 0.0_wp
+ H_abs(N_active) = 0.0_wp
+ thick(N_active) = 0.0_wp
+ bgc_temp(N_active,:) = 0.0_wp
+
+ N_active = N_active-1
+
+ ELSE IF (N_active>N_top .AND. N_active<=Nlayer .AND. thick(N_top+1)/thick_0<1.00001_wp) THEN
+ !Removing bottom layer and moving all other layers if N_active<Nlayer .and. N_active>N_top.
+ !Is also activated if thick(N_top+1)=thick_0
+ DO k = N_top,N_active-1
+ m(k) = rho (k+1) *thick(k+1) ! thick_0 changed by NF, 2024
+ S_abs(k) = S_bu(k+1) *rho(k+1)*thick(k+1) ! thick_0 changed by NF, 2024
+ H_abs(k) = H (k+1) *rho(k+1)*thick(k+1) ! thick_0 changed by NF, 2024
+ bgc_temp(k,:) = bgc_bulk (k+1,:) *rho(k+1)*thick(k+1) ! thick_0 changed by NF, 2024
+ END DO
+ m(N_active) = 0.0_wp
+ S_abs(N_active) = 0.0_wp
+ H_abs(N_active) = 0.0_wp
+ bgc_temp(N_active,:) = 0.0_wp
+ thick(N_active) = 0.0_wp
+
+
+ N_active = N_active-1
+ END IF
+
+ IF (N_active==Nlayer .AND. thick(N_top+1)-thick_0>=0.000001_wp ) THEN
+ !Middle layers are adjusted__________________________
+ loss_m = thick_0*rho(N_top+1)
+ loss_S_abs = loss_m*S_bu(N_top+1)
+ loss_H_abs = loss_m*H(N_top+1)
+ loss_bgc = loss_m*bgc_bulk(N_top+1,:)
+
+ !layer N_top needs new values________________________
+ m(N_top) = loss_m
+ S_abs(N_top) = loss_S_abs
+ H_abs(N_top) = loss_H_abs
+ bgc_temp(N_top,:) = loss_bgc(:)
+
+
+ DO k=N_top+1,N_middle+N_top
+ m(k) = m(k) -loss_m
+ H_abs(k) = H_abs(k) -loss_H_abs
+ S_abs(k) = S_abs(k) -loss_S_abs
+ bgc_temp(k,:) = bgc_temp(k,:) -loss_bgc(:)
+
+ shift = thick_0*REAL(N_middle-k+N_top)/REAL(N_middle) !distance the border moves
+
+ loss_m = shift *rho(k+1)
+ loss_S_abs = loss_m *S_bu(k+1)
+ loss_H_abs = loss_m *H(k+1)
+ loss_bgc = loss_m *bgc_bulk(k+1,:)
+
+ m(k) = m(k) +loss_m
+ H_abs(k) = H_abs(k) +loss_H_abs
+ S_abs(k) = S_abs(k) +loss_S_abs
+ bgc_temp(k,:) = bgc_temp(k,:) +loss_bgc(:)
+ END DO
+
+
+
+
+
+ !Middle layer thickness is adjusted________________________
+ DO k = N_top+1,N_top+N_middle
+ thick(k) = thick(k)-thick_0/REAL(N_middle)
+ END DO
+ END IF
+
+
+ !Checks if current thick vector is possible for current N_active and thick_0
+ !Uses 0.501 instead of 0.500 to give a bit of wriggle room for numerical noise
+ IF(thick_0*(N_active+0.501_wp)<=SUM(thick) .AND. N_active<Nlayer) THEN
+ PRINT*,'wtf layer issue',thick_0,SUM(thick),N_active,Nlayer
+ STOP 7889
+ END IF
+
+ IF (PRESENT(bgc_abs)) THEN
+ bgc_abs = bgc_temp
+ END IF
+ END SUBROUTINE top_melt
+
+
+ !>
+ !! Controls melting at the bottom.
+ !!
+ !! Subroutine is called if the second lowest layer is fully liquid and the middle layers are thicker then the bottom layers.
+ !! When this occurs the thickness of the middle layers is reduced by 1/N_middle*thick_0.
+ !! The lower layers are adjusted upwards, layer(k)=layer(k-1).
+ !!
+ !!
+ !! @par Revision History
+ !! Started by Philipp Griewank, IMPRS (2010-07-28) \n
+ !! Linear profiles replaced with simple upstream by Philipp Griewank, IMPRS (2010-08-25)
+ !!
+ SUBROUTINE bottom_melt (Nlayer,N_middle,N_top,m,S_abs,H_abs,thick,N_bgc,bgc_abs)
+
+ INTEGER, INTENT(in) :: Nlayer,N_middle,N_top
+ REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: m,S_abs,H_abs,thick
+ INTEGER :: k
+ REAL(wp) :: loss_m,loss_S_abs,loss_H_abs,shift
+ REAL(wp), DIMENSION(Nlayer) :: rho !>density of layer
+ REAL(wp), DIMENSION(Nlayer) :: H,S_bu
+ INTEGER, INTENT(in) :: N_bgc
+ REAL(wp), DIMENSION(Nlayer,N_bgc), INTENT(inout), OPTIONAL :: bgc_abs
+ REAL(wp), DIMENSION(N_bgc) :: loss_bgc
+ REAL(wp), DIMENSION(Nlayer,N_bgc) :: bgc_temp, bgc_bulk
+
+
+ IF (PRESENT(bgc_abs)) THEN
+ bgc_temp = bgc_abs
+ ELSE
+ bgc_temp = 0._wp
+ END IF
+
+
+ !Calculating rho, H and S_bu
+ DO k = N_top+1,Nlayer
+ rho(k) = m(k)/thick(k)
+ S_bu(k) = S_abs(k)/m(k)
+ H(k) = H_abs(k)/m(k)
+ bgc_bulk(k,:) = bgc_temp(k,:)/m(k)
+
+ END DO
+
+ loss_m = 0._wp
+ loss_S_abs = 0._wp
+ loss_H_abs = 0._wp
+ loss_bgc = 0._wp
+
+ !_________________________Middle layers are adjusted__________________________
+ DO k = N_top+1,N_top+N_middle
+
+ !each layer receives all that the layer above loses
+ m(k) = m(k) +loss_m
+ H_abs(k) = H_abs(k) +loss_H_abs
+ S_abs(k) = S_abs(k) +loss_S_abs
+ bgc_temp(k,:) = bgc_temp(k,:) +loss_bgc(:)
+
+
+
+ shift = thick(Nlayer)*(k-N_top)/ REAL(N_middle) !distance the border moves
+
+ loss_m = shift*rho(k)
+ loss_H_abs = loss_m*H(k)
+ loss_S_abs = loss_m*S_bu(k)
+ loss_bgc(:) = loss_m*bgc_bulk(k,:)
+
+
+ m(k) = m(k) -loss_m
+ H_abs(k) = H_abs(k) -loss_H_abs
+ S_abs(k) = S_abs(k) -loss_S_abs
+ bgc_temp(k,:) = bgc_temp(k,:) -loss_bgc(:)
+
+
+ END DO
+
+
+ !________________________Middle layer thickness is adjusted
+
+ DO k = N_top+1,N_top+N_middle
+ thick(k) = thick(k)-thick(Nlayer)/REAL(N_middle)
+ END DO
+ !_________________________Bottom layers are adjusted__________________________
+
+ DO k = N_top+N_middle+1,Nlayer
+ H_abs(k) = rho(k-1)*thick(k)*H(k-1)
+ S_abs(k) = rho(k-1)*thick(k)*S_bu(k-1)
+ m(k) = rho(k-1)*thick(k)
+ bgc_temp(k,:) = rho(k-1)*thick(k)*bgc_bulk(k-1,:)
+ END DO
+
+ IF (PRESENT(bgc_abs)) THEN
+ bgc_abs = bgc_temp
+ END IF
+
+ END SUBROUTINE bottom_melt
+
+
+
+ !>
+ !! Controls bottom growth.
+ !!
+ !! Subroutine is called if the lowest layer is not fully liquid and the middle layers are thicker then thick_0.
+ !! When this occurs the thickness of the middle layers is increased by 1/N_middle*thick(bottom).
+ !! The lower layers are adjusted downwards, layer(k)=layer(k+1).
+ !!
+ !!
+ !! @par Revision History
+ !! Started by Philipp Griewank, IMPRS (2010-07-29) \n
+ !! Linear profiles removed, replaced with simple upwind by Philipp Griewank, IMPRS (2010-08-25)
+ !!
+ SUBROUTINE bottom_growth (Nlayer,N_bottom,N_middle,N_top,m,S_abs,H_abs,thick,T_bottom,S_bu_bottom,N_bgc,bgc_abs,bgc_bottom)
+
+ INTEGER, INTENT(in) :: Nlayer,N_bottom,N_middle,N_top
+ REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: m,S_abs,H_abs,thick
+ REAL(wp), INTENT(in) :: T_bottom,S_bu_bottom
+ INTEGER :: k
+ REAL(wp) :: gain_m,gain_S_abs,gain_H_abs,shift
+ REAL(wp), DIMENSION(Nlayer) :: rho !>density of layer
+ REAL(wp), DIMENSION(Nlayer) :: H,S_bu
+ INTEGER, INTENT(in) :: N_bgc
+ REAL(wp), DIMENSION(Nlayer,N_bgc), INTENT(inout), OPTIONAL :: bgc_abs
+ REAL(wp), DIMENSION(N_bgc), INTENT(in), OPTIONAL :: bgc_bottom
+ REAL(wp), DIMENSION(N_bgc) :: gain_bgc
+ REAL(wp), DIMENSION(Nlayer,N_bgc) :: bgc_temp, bgc_bulk
+
+
+ IF (PRESENT(bgc_abs)) THEN
+ bgc_temp = bgc_abs
+ ELSE
+ bgc_temp = 0._wp
+ END IF
+
+
+ !Calculating rho, H and S_bu
+ DO k = N_top+1,N_top+N_middle+1
+ rho(k) = m(k)/thick(k)
+ S_bu(k) = S_abs(k)/m(k)
+ H(k) = H_abs(k)/m(k)
+ bgc_bulk(k,:) = bgc_temp(k,:)/m(k)
+ END DO
+
+ gain_m = 0._wp
+ gain_S_abs = 0._wp
+ gain_H_abs = 0._wp
+ gain_bgc = 0._wp
+
+ !_________________________Middle layers are adjusted__________________________
+ DO k = N_top+1,N_top+N_middle
+
+ !each layer loses all that the upper layer gained
+ m(k) = m(k) -gain_m
+ H_abs(k) = H_abs(k) -gain_H_abs
+ S_abs(k) = S_abs(k) -gain_S_abs
+ bgc_temp(k,:) = bgc_temp(k,:) -gain_bgc(:)
+
+ shift = thick(Nlayer)*(k-N_top)/ REAL(N_middle) !distance the border moves
+
+ gain_m = shift*rho(k+1)
+ gain_H_abs = gain_m*H(k+1)
+ gain_S_abs = gain_m*S_bu(k+1)
+ gain_bgc(:)= gain_m*bgc_bulk(k+1,:)
+
+ m(k) = m(k) +gain_m
+ H_abs(k) = H_abs(k) +gain_H_abs
+ S_abs(k) = S_abs(k) +gain_S_abs
+ bgc_temp(k,:) = bgc_temp(k,:) +gain_bgc(:)
+
+ END DO
+
+ !________________________Middle layer thickness is adjusted
+ DO k=N_top+1,N_top+N_middle
+ thick(k) = thick(k)+thick(Nlayer)/REAL(N_middle)
+ END DO
+
+ !_________________________Bottom layers are adjusted__________________________
+ DO k = Nlayer-N_bottom+1,Nlayer-1
+ H_abs(k) = H_abs(k+1)
+ S_abs(k) = S_abs(k+1)
+ m(k) = m(k+1)
+ bgc_temp(k,:) = bgc_temp(k+1,:)
+ END DO
+
+ !_________________________Lowest layer is set using T_bottom and S_bu_bottom
+ m(Nlayer) = thick(Nlayer)*rho_l
+ H_abs(Nlayer) = m(Nlayer)*T_bottom*c_l
+ S_abs(Nlayer) = m(Nlayer)*S_bu_bottom
+
+ IF (PRESENT(bgc_abs)) THEN
+ bgc_temp(Nlayer,:) = m(Nlayer)*bgc_bottom(:)
+ bgc_abs = bgc_temp
+ END IF
+
+ END SUBROUTINE bottom_growth
+
+
+
+
+ !>
+ !! Controls growth if the middle layers are not all active.
+ !!
+ !! Subroutine is called if the lowest layer is not fully liquid and the middle layers are not all activated.
+ !! A new layer is activated.
+ !! The new lowest layer is assigned standard values of T_bottom and S_bu_bottom.
+ !!
+ !!
+ !! @par Revision History
+ !! Started by Philipp Griewank, IMPRS (2010-07-29) \n
+ !! Expanded to deal with not full top layers by Philipp Griewank, IMPRS (2011-01-13) \n
+ !! Expansion removed and bgc added by Philipp Griewank, IMPRS (2014-02-04) \n
+ SUBROUTINE bottom_growth_simple (N_active,Nlayer,m,S_abs,H_abs,thick,thick_0,T_bottom,S_bu_bottom,N_bgc,bgc_abs,bgc_bottom)
+
+ INTEGER, INTENT(in) :: Nlayer
+ INTEGER, INTENT(inout) :: N_active
+ REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: m,S_abs,H_abs,thick
+ REAL(wp), INTENT(in) :: thick_0
+ REAL(wp), INTENT(in) :: T_bottom,S_bu_bottom
+ INTEGER, INTENT(in) :: N_bgc
+ REAL(wp), DIMENSION(Nlayer,N_bgc), INTENT(inout), OPTIONAL :: bgc_abs
+ REAL(wp), DIMENSION(N_bgc), INTENT(in), OPTIONAL :: bgc_bottom
+
+
+ N_active = N_active+1
+ thick(N_active) = thick_0
+
+ m(N_active) = thick(N_active)*rho_l
+ H_abs(N_active) = m(N_active) *T_bottom*c_l
+ S_abs(N_active) = m(N_active) *S_bu_bottom
+
+ IF (PRESENT(bgc_abs)) THEN
+ bgc_abs(N_active,:) = bgc_bottom(:)*m(N_active)
+ END IF
+
+ END SUBROUTINE bottom_growth_simple
+
+ !>
+ !! Controls bottom melting if the middle layers are not all active.
+ !!
+ !! Subroutine is called if the second lowest active layer is fully liquid and the middle layers are not all activated, or if the middle layers are thick_0 thick.
+ !! The lowest layer is simply deactivated.
+ !!
+ !!
+ !! @par Revision History
+ !! Started by Philipp Griewank, IMPRS (2010-07-30>)
+ !!
+ SUBROUTINE bottom_melt_simple(N_active,Nlayer,m,S_abs,H_abs,thick,N_bgc,bgc_abs)
+
+ INTEGER, INTENT(in) :: Nlayer
+ INTEGER, INTENT(inout) :: N_active
+ REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: m,S_abs,H_abs,thick
+ INTEGER, INTENT(in) :: N_bgc
+ REAL(wp), DIMENSION(Nlayer,N_bgc), INTENT(inout), OPTIONAL :: bgc_abs
+
+ thick(N_active) = 0._wp
+ m(N_active) = 0._wp
+ S_abs(N_active) = 0._wp
+ H_abs(N_active) = 0._wp
+ IF(PRESENT(bgc_abs)) THEN
+ bgc_abs(N_active,:) = 0._wp
+ END IF
+ N_active = N_active-1
+
+ END SUBROUTINE bottom_melt_simple
+
+
+
+ !>
+ !! Top grow subroutine.
+ !!
+ !! Should be called when the top layer is thicker then 1.5 *thick_0.
+ !! If N_active=Nlayer middle layers are expanded by thick_0/N_middle and top layers are moved one down.
+ !! IF N_active<Nlayer then N_active=N_active+1 and all layers are shifted downwards.
+ !!
+ !!
+ !! @par Revision History
+ !! Started by Philipp Griewank, IMPRS (2011-05-10>)
+ !! edited by Niels Fuchs, UHH (2024-05-03), changed according to changes in function top_melt
+
+
+ SUBROUTINE top_grow (Nlayer,N_active,N_middle,N_top,thick_0,m,S_abs,H_abs,thick,N_bgc,bgc_abs)
+
+ INTEGER, INTENT(in) :: Nlayer,N_middle,N_top
+ INTEGER, INTENT(inout) :: N_active
+ REAL(wp), INTENT(in) :: thick_0
+ REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: m,S_abs,H_abs,thick
+ INTEGER :: k
+ REAL(wp) :: loss_m,loss_S_abs,loss_H_abs,shift
+ REAL(wp), DIMENSION(Nlayer) :: rho !>density of layer
+ REAL(wp), DIMENSION(Nlayer) :: H,S_bu
+ INTEGER, INTENT(in) :: N_bgc
+ REAL(wp), DIMENSION(Nlayer,N_bgc), INTENT(inout), OPTIONAL :: bgc_abs
+ REAL(wp), DIMENSION(N_bgc) :: loss_bgc
+ REAL(wp), DIMENSION(Nlayer,N_bgc) :: bgc_temp, bgc_bulk
+
+ IF (PRESENT(bgc_abs)) THEN
+ bgc_temp = bgc_abs
+ ELSE
+ bgc_temp = 0._wp
+ END IF
+ !Calculating rho, H and S_bu
+ DO k = 1,N_active
+ rho(k) = m(k)/thick(k)
+ S_bu(k) = S_abs(k)/m(k)
+ H(k) = H_abs(k)/m(k)
+ bgc_bulk(k,:) = bgc_temp(k,:)/m(k)
+ END DO
+
+ !Adjusting the top layer from thick(1)>1.5*thick_0(1) to thick(1)=thick(1)-thick_0(1)
+ loss_m = thick_0*rho(1)
+ loss_S_abs = loss_m*S_bu(1)
+ loss_H_abs = loss_m*H(1)
+ loss_bgc(:) = loss_m*bgc_bulk(1,:)
+
+ m(1) = m(1) -loss_m
+ S_abs(1) = S_abs(1) -loss_S_abs
+ H_abs(1) = H_abs(1) -loss_H_abs
+ bgc_temp(1,:) = bgc_temp(1,:) -loss_bgc(:)
+ thick(1) = thick(1) -thick_0
+
+ !Adjusting top layers
+ DO k = 2,MIN(N_top,N_active)
+ m(k) = rho (k-1) *thick_0
+ S_abs(k) = S_bu(k-1) *rho(k-1)*thick_0
+ H_abs(k) = H (k-1) *rho(k-1)*thick_0
+ bgc_temp(k,:) = bgc_bulk(k-1,:) *rho(k-1)*thick_0
+ END DO
+
+ !Adding new bottom layer if N_active<=N_top+1
+ IF (N_active<=N_top) THEN
+ N_active = N_active+1
+ m(N_active) = rho (N_active-1) *thick_0
+ S_abs(N_active) = S_bu(N_active-1) *thick_0*rho(N_active-1)
+ H_abs(N_active) = H (N_active-1) *thick_0*rho(N_active-1)
+ bgc_temp(N_active,:) = bgc_bulk(N_active-1,:) *thick_0*rho(N_active-1)
+ thick(N_active) = thick_0
+
+ !Adding new bottom layer and moving all other layers if N_active<Nlayer .and. N_active>N_top
+ ELSE IF (N_active>N_top .AND. N_active<Nlayer) THEN
+ DO k = N_top+1,N_active
+ m(k) = rho (k-1) *thick(k-1) !thick_0
+ S_abs(k) = S_bu(k-1) *rho(k-1)*thick(k-1) !thick_0
+ H_abs(k) = H (k-1) *rho(k-1)*thick(k-1) !thick_0
+ bgc_temp(k,:) = bgc_bulk(k-1,:) *rho(k-1)*thick(k-1) !thick_0
+ END DO
+ N_active = N_active+1
+ m(N_active) = rho (N_active-1) *thick_0
+ S_abs(N_active) = S_bu(N_active-1) *thick_0 *rho(N_active-1)
+ H_abs(N_active) = H (N_active-1) *thick_0 *rho(N_active-1)
+ bgc_temp(N_active,:) = bgc_bulk(N_active-1,:) *thick_0 *rho(N_active-1)
+ thick(N_active) = thick_0
+
+ ELSE IF (N_active==Nlayer) THEN
+ !Middle layers are adjusted__________________________
+ loss_m = thick_0*rho(N_top)
+ loss_S_abs = loss_m*S_bu(N_top)
+ loss_H_abs = loss_m*H(N_top)
+ loss_bgc(:) = loss_m*bgc_bulk(N_top,:)
+
+ DO k = N_top+1,N_middle+N_top
+ m(k) = m(k) +loss_m
+ H_abs(k) = H_abs(k) +loss_H_abs
+ S_abs(k) = S_abs(k) +loss_S_abs
+ bgc_temp(k,:) = bgc_temp(k,:) +loss_bgc(:)
+
+ shift = thick_0*REAL(N_middle-k+N_top)/REAL(N_middle) !distance the border moves
+
+ loss_m = shift *rho(k)
+ loss_S_abs = loss_m *S_bu(k)
+ loss_H_abs = loss_m *H(k)
+ loss_bgc(:) = loss_m *bgc_bulk(k,:)
+
+ m(k) = m(k) -loss_m
+ H_abs(k) = H_abs(k) -loss_H_abs
+ S_abs(k) = S_abs(k) -loss_S_abs
+ bgc_temp(k,:) = bgc_temp(k,:) -loss_bgc(:)
+ END DO
+
+ !Middle layer thickness is adjusted
+ DO k = N_top+1,N_top+N_middle
+ thick(k) = thick(k)+thick_0/REAL(N_middle)
+ END DO
+ END IF
+
+
+ IF (PRESENT(bgc_abs)) THEN
+ bgc_abs = bgc_temp
+ END IF
+
+ END SUBROUTINE top_grow
+
+
+END MODULE mo_layer_dynamics
diff --git a/mo_mass.f90 b/mo_mass.f90
index fb0ea511b389aa4cb07e618f0d6e2708b1fdcd39..7cc3f158e66ef461014c247fa751912d9b23d8b9 100644
--- a/mo_mass.f90
+++ b/mo_mass.f90
@@ -1,215 +1,218 @@
-!>
-!! Regulates mass transfers and their results.
-!!
-!! Ultimately all processes which involve a mass flux should be stored here.
-!!
-!!
-!! @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
-!! Begin implementing Expulsion by Philipp Griewank, IMPRS (2010-08-24)
-!!
-MODULE mo_mass
-
- USE mo_parameters, ONLY:wp,rho_l,c_l
- USE mo_thermo_functions, ONLY:func_S_br
- USE mo_functions, ONLY:func_sat_O2
- IMPLICIT NONE
-
- PUBLIC :: mass_transfer,expulsion_flux,bgc_advection
-
-
-CONTAINS
-
- !>
- !! Calculates the effects of mass transfers on H_abs and S_abs.
- !!
- !! The effects of brine displaced by expulsion, flushing or drainage expansion lead to changes in mass, salt ans enthalpy.
- !! This subroutine calculates the effects on S_abs and H_abs.
- !! A very simple upwind strategy is employed, Brine from below has T and S_br of the lower layer, and brine from above T and S_br of the upper layer.
- !! To avoid negative salinity, the maximum amount of advective salt is the total salt content of the layer.
- !! The amount of mass transfered is calculated in other subroutines.
- !!
- !! This subroutine was started as a quick and dirty way to simulate the bottom freezing experiment described in Notz 2005 p. 85
- !! IMPORTANT: Before this subroutine expelled brine was removed from the system and its effects were determined in subroutine expulsion.
- !! S_bu must be up to date!
- !!
- !! @par Revision History
- !! Brought to life by Philipp Griewank, IMPRS (2010-08-24) \n
- !! Modified to work with all processes by Philipp Griewank, IMPRS (2010-11-27)
- SUBROUTINE mass_transfer (Nlayer,N_active,T,H_abs,S_abs,S_bu,T_bottom,S_bu_bottom,fl_m)
-
- IMPLICIT NONE
-
- INTEGER, INTENT(in) :: Nlayer, N_active
- REAL(wp), INTENT(in) :: T_bottom,S_bu_bottom
- REAL(wp), DIMENSION(Nlayer), INTENT(in) :: T,S_bu
- REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: H_abs,S_abs
- REAL(wp), DIMENSION(Nlayer+1), INTENT(in) :: fl_m
- REAL(wp), DIMENSION(Nlayer+1) :: TT,SS_bu,SS_abs !Same as T and S_bu but expanded to include bottom values
- INTEGER :: k
-
-
- TT(1:N_active) = T(1:N_active)
- SS_bu(1:N_active) = S_bu(1:N_active)
- SS_abs(1:N_active) = S_abs(1:N_active)
-
- TT(N_active+1) = T_bottom
- SS_bu(N_active+1) = S_bu_bottom
- SS_abs(N_active+1) = S_bu_bottom*2000._wp
-
-
-
- DO k = 1,N_active
-
- IF (fl_m(k+1)>0.) THEN
- H_abs(k) = H_abs(k) +fl_m(k+1)*TT(k+1)*c_l
- S_abs(k) = S_abs(k) +MIN(fl_m(k+1)*func_S_br(TT(k+1),SS_bu(k+1)), SS_abs(k+1))
- ELSE IF (fl_m(k+1)<0.) THEN
- H_abs(k) = H_abs(k) +fl_m(k+1)*TT(k)*c_l
- S_abs(k) = S_abs(k) +MAX(fl_m(k+1)*func_S_br(TT(k),SS_bu(k)), -S_abs(k))
- END IF
-
- IF (fl_m(k)>0.) THEN
- H_abs(k) = H_abs(k) -fl_m(k)*TT(k)*c_l
- S_abs(k) = S_abs(k) -MIN(fl_m(k)*func_S_br(TT(k),SS_bu(k)), S_abs(k))
- END IF
- END DO
-
- DO k = 2, N_active
-
- IF (fl_m(k)<0) THEN
- H_abs(k) = H_abs(k) -fl_m(k)*TT(k-1)*c_l
- S_abs(k) = S_abs(k) -MAX(fl_m(k)*func_S_br( TT(k-1),SS_bu(k-1)), -S_abs(k-1))
- END IF
- END DO
- END SUBROUTINE mass_transfer
-
-
-
- !>
- !! Generates the fluxes caused by expulsion.
- !!
- !! Brine displaced by expansion of a freezing mushy layer lead to a mass, enthalpy and salt flux.
- !! This subroutine calculates the amount of brine which moves between the layers caused by V_ex and how the mass in the layers changes.
- !! Vary basic assumptions are made. Brine always moves downward (negative), no horizontal movement are allowed and gas pockets can be filled.
- !! The upper boundary layer is not permeable but the bottom one is.
- !! This subroutine was started as a quick and dirty way to simulate the bottom freezing experiment described in Notz 2005 p. 85
- !!
- !! @par Revision History
- !! Brought to life by Philipp Griewank, IMPRS (2010-08-24) \n
- !! Simplified by Philipp Griewank, IMPRS (2010-11-27)
- SUBROUTINE expulsion_flux (thick,V_ex,Nlayer,N_active,psi_g,fl_m,m)
-
- INTEGER, INTENT(in) :: Nlayer, N_active
- REAL(wp), DIMENSION(Nlayer), INTENT(in) :: V_ex,thick
- REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: psi_g,m
- REAL(wp), DIMENSION(Nlayer+1), INTENT(out) :: fl_m
-
- INTEGER::k
-
- fl_m(1:Nlayer+1) = 0._wp
- fl_m(2) = -V_ex(1)*rho_l
- DO k = 2,N_active
- IF (psi_g(k)<0.001) THEN
- fl_m(k+1) = -V_ex(k)*rho_l+fl_m(k)
- ELSE
- fl_m(k+1) = -MAX((V_ex(k)-psi_g(k)*thick(k))*rho_l ,0.0_wp)
- psi_g(k) = MAX((psi_g(k)*thick(k)-V_ex(k))/thick(k) ,0.0_wp)
- END IF
- END DO
-
- DO k = 1,N_active
- m(k) = m(k) +fl_m(k+1)-fl_m(k)
- END DO
-
- END SUBROUTINE expulsion_flux
-
- !>
- !! Calculates how the brine fluxes stored in fl_brine_bgc advect bgc tracers
- !!
- !! A very simple upwind strategy is employed.
- !! To avoid negative tracer densities, the maximum amount of advection is restricted to the current tracer content in a layer divided by three.
- !! Three is chosen as a limit as currently each layer can have a maximum of three flows leaving the layer (to the layer above, the layer below, and the lowest layer).
- !! The advection scheme is likely overly diffusive, but given the limitations we are working with (e.g. changing brine volumes) nothing more sophisticated can be applied easily.
- !!
- !! For gases it might make sense to limit the brine density to saturation value in advecting brine, to take bubble formation into account. This needs to be specified in bgc_advection, and is a first attempt (both scientifically and code wise) which should be used with caution!
- !!
- !! @par Revision History
- !! Brought to life by Philipp Griewank, IMPRS (2014-02-10)
- SUBROUTINE bgc_advection (Nlayer,N_active,N_bgc,fl_brine_bgc,bgc_abs,psi_l,thick,bgc_bottom)
-
- IMPLICIT NONE
-
- INTEGER, INTENT(in) :: Nlayer,N_active,N_bgc
- REAL(wp), DIMENSION(Nlayer), INTENT(in) :: psi_l,thick
- REAL(wp), DIMENSION(N_bgc), INTENT(in) :: bgc_bottom
- REAL(wp), DIMENSION(Nlayer+1,Nlayer+1), INTENT(in) :: fl_brine_bgc
- REAL(wp), DIMENSION(Nlayer,N_bgc), INTENT(inout) :: bgc_abs
- REAL(wp), DIMENSION(Nlayer,N_bgc) :: bgc_temp
- REAL(wp), DIMENSION(Nlayer,N_bgc) :: bgc_br
- REAL(wp), DIMENSION(N_bgc) :: flux
- INTEGER ::k,i,j
-
- bgc_temp = bgc_abs
-
- !Calculating brine concentration of advecting brine, saturation limitations must be included manually!
- DO k = 1,N_active
- bgc_br(k,:) = bgc_abs(k,:)/(MAX(psi_l(k)*thick(k)*rho_l,0.000000000000001_wp))
-
- !How to limit the first tracer to the oxygen saturation limit
- !sat_O2 = func_sat_O2(T(k),S_abs(k)/m(k))
- !bgc_br(k,1) = MIN(bgc_br(k,1),1.25*sat_O2)
- END DO
-
- !Internal flows
- DO i = 1,N_active
- DO j = 1,N_active
- DO k = 1,N_bgc
- IF (fl_brine_bgc(i,j)*bgc_br(i,k) > bgc_abs(i,k)/3._wp) THEN
- END IF
- flux(k) = MIN(fl_brine_bgc(i,j)*bgc_br(i,k),bgc_abs(i,k)/3._wp)
- END DO
- bgc_temp(i,:) = bgc_temp(i,:) - flux(:)
- bgc_temp(j,:) = bgc_temp(j,:) + flux(:)
-
- END DO
- END DO
-
- !Flows which leave the domain
- DO i = 1,N_active
- j = N_active+1
- DO k = 1,N_bgc
- flux(k) = MIN(fl_brine_bgc(i,j)*bgc_br(i,k),bgc_abs(i,k)/3._wp)
- END DO
- bgc_temp(i,:) = bgc_temp(i,:) - flux(:)
- END DO
-
- !Flows which enter the domain
- i = N_active+1
- DO j = 1,N_active
- DO k = 1,N_bgc
- flux(k) = fl_brine_bgc(i,j)*bgc_bottom(k)
- END DO
- bgc_temp(j,:) = bgc_temp(j,:) + flux(:)
- END DO
-
- bgc_abs = bgc_temp
- END SUBROUTINE bgc_advection
-
-
-
-END MODULE mo_mass
+!>
+!! Regulates mass transfers and their results.
+!!
+!! Ultimately all processes which involve a mass flux should be stored here.
+!!
+!!
+!! @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
+!! Begin implementing Expulsion by Philipp Griewank, IMPRS (2010-08-24)
+!!
+MODULE mo_mass
+
+ USE mo_parameters, ONLY:wp,rho_l,c_l
+ USE mo_thermo_functions, ONLY:func_S_br
+ USE mo_functions, ONLY:func_sat_O2
+ IMPLICIT NONE
+
+ PUBLIC :: mass_transfer,expulsion_flux,bgc_advection
+
+
+CONTAINS
+
+ !>
+ !! Calculates the effects of mass transfers on H_abs and S_abs.
+ !!
+ !! The effects of brine displaced by expulsion, flushing or drainage expansion lead to changes in mass, salt ans enthalpy.
+ !! This subroutine calculates the effects on S_abs and H_abs.
+ !! A very simple upwind strategy is employed, Brine from below has T and S_br of the lower layer, and brine from above T and S_br of the upper layer.
+ !! To avoid negative salinity, the maximum amount of advective salt is the total salt content of the layer.
+ !! The amount of mass transfered is calculated in other subroutines.
+ !!
+ !! This subroutine was started as a quick and dirty way to simulate the bottom freezing experiment described in Notz 2005 p. 85
+ !! IMPORTANT: Before this subroutine expelled brine was removed from the system and its effects were determined in subroutine expulsion.
+ !! S_bu must be up to date!
+ !!
+ !! @par Revision History
+ !! Brought to life by Philipp Griewank, IMPRS (2010-08-24) \n
+ !! Modified to work with all processes by Philipp Griewank, IMPRS (2010-11-27)
+ SUBROUTINE mass_transfer (Nlayer,N_active,T,H_abs,S_abs,S_bu,T_bottom,S_bu_bottom,fl_m)
+
+ IMPLICIT NONE
+
+ INTEGER, INTENT(in) :: Nlayer, N_active
+ REAL(wp), INTENT(in) :: T_bottom,S_bu_bottom
+ REAL(wp), DIMENSION(Nlayer), INTENT(in) :: T,S_bu
+ REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: H_abs,S_abs
+ REAL(wp), DIMENSION(Nlayer+1), INTENT(in) :: fl_m
+ REAL(wp), DIMENSION(Nlayer+1) :: TT,SS_bu,SS_abs !Same as T and S_bu but expanded to include bottom values
+ INTEGER :: k
+
+
+ TT(1:N_active) = T(1:N_active)
+ SS_bu(1:N_active) = S_bu(1:N_active)
+ SS_abs(1:N_active) = S_abs(1:N_active)
+
+ TT(N_active+1) = T_bottom
+ SS_bu(N_active+1) = S_bu_bottom
+ SS_abs(N_active+1) = S_bu_bottom*2000._wp
+
+
+
+ DO k = 1,N_active
+
+ IF (fl_m(k+1)>0.) THEN
+ H_abs(k) = H_abs(k) +fl_m(k+1)*TT(k+1)*c_l
+ S_abs(k) = S_abs(k) +MIN(fl_m(k+1)*func_S_br(TT(k+1),SS_bu(k+1)), SS_abs(k+1))
+ ELSE IF (fl_m(k+1)<0.) THEN
+ H_abs(k) = H_abs(k) +fl_m(k+1)*TT(k)*c_l
+ S_abs(k) = S_abs(k) +MAX(fl_m(k+1)*func_S_br(TT(k),SS_bu(k)), -S_abs(k))
+ END IF
+
+ IF (fl_m(k)>0.) THEN
+ H_abs(k) = H_abs(k) -fl_m(k)*TT(k)*c_l
+ S_abs(k) = S_abs(k) -MIN(fl_m(k)*func_S_br(TT(k),SS_bu(k)), S_abs(k))
+ END IF
+ END DO
+
+ DO k = 2, N_active
+
+ IF (fl_m(k)<0) THEN
+ H_abs(k) = H_abs(k) -fl_m(k)*TT(k-1)*c_l
+ S_abs(k) = S_abs(k) -MAX(fl_m(k)*func_S_br( TT(k-1),SS_bu(k-1)), -S_abs(k-1))
+ END IF
+ END DO
+ END SUBROUTINE mass_transfer
+
+
+
+ !>
+ !! Generates the fluxes caused by expulsion.
+ !!
+ !! Brine displaced by expansion of a freezing mushy layer lead to a mass, enthalpy and salt flux.
+ !! This subroutine calculates the amount of brine which moves between the layers caused by V_ex and how the mass in the layers changes.
+ !! Vary basic assumptions are made. Brine always moves downward (negative), no horizontal movement are allowed and gas pockets can be filled.
+ !! The upper boundary layer is not permeable but the bottom one is.
+ !! This subroutine was started as a quick and dirty way to simulate the bottom freezing experiment described in Notz 2005 p. 85
+ !!
+ !! @par Revision History
+ !! Brought to life by Philipp Griewank, IMPRS (2010-08-24) \n
+ !! Simplified by Philipp Griewank, IMPRS (2010-11-27)
+ !! edited by Niels Fuchs, UHH (2024-05-03), added brine mass removal from H_abs
+ SUBROUTINE expulsion_flux (thick,V_ex,Nlayer,N_active,psi_g,fl_m,m, H_abs)
+
+ INTEGER, INTENT(in) :: Nlayer, N_active
+ REAL(wp), DIMENSION(Nlayer), INTENT(in) :: V_ex,thick
+ REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: psi_g,m, H_abs
+ REAL(wp), DIMENSION(Nlayer+1), INTENT(out) :: fl_m
+
+ INTEGER::k
+
+ H = H_abs(:) / m(:)
+ fl_m(1:Nlayer+1) = 0._wp
+ fl_m(2) = -V_ex(1)*rho_l
+ DO k = 2,N_active
+ IF (psi_g(k)<0.001) THEN
+ fl_m(k+1) = -V_ex(k)*rho_l+fl_m(k)
+ ELSE
+ fl_m(k+1) = -MAX((V_ex(k)-psi_g(k)*thick(k))*rho_l ,0.0_wp)
+ psi_g(k) = MAX((psi_g(k)*thick(k)-V_ex(k))/thick(k) ,0.0_wp)
+ END IF
+ END DO
+
+ DO k = 1,N_active
+ m(k) = m(k) +fl_m(k+1)-fl_m(k)
+ H_abs(k) = H_abs(k) + (fl_m(k+1)-fl_m(k))*H(k) ! exp. brine mass must be removed from H_abs, temperature loss is calculated in mass transfer subs.
+ END DO
+
+ END SUBROUTINE expulsion_flux
+
+ !>
+ !! Calculates how the brine fluxes stored in fl_brine_bgc advect bgc tracers
+ !!
+ !! A very simple upwind strategy is employed.
+ !! To avoid negative tracer densities, the maximum amount of advection is restricted to the current tracer content in a layer divided by three.
+ !! Three is chosen as a limit as currently each layer can have a maximum of three flows leaving the layer (to the layer above, the layer below, and the lowest layer).
+ !! The advection scheme is likely overly diffusive, but given the limitations we are working with (e.g. changing brine volumes) nothing more sophisticated can be applied easily.
+ !!
+ !! For gases it might make sense to limit the brine density to saturation value in advecting brine, to take bubble formation into account. This needs to be specified in bgc_advection, and is a first attempt (both scientifically and code wise) which should be used with caution!
+ !!
+ !! @par Revision History
+ !! Brought to life by Philipp Griewank, IMPRS (2014-02-10)
+ SUBROUTINE bgc_advection (Nlayer,N_active,N_bgc,fl_brine_bgc,bgc_abs,psi_l,thick,bgc_bottom)
+
+ IMPLICIT NONE
+
+ INTEGER, INTENT(in) :: Nlayer,N_active,N_bgc
+ REAL(wp), DIMENSION(Nlayer), INTENT(in) :: psi_l,thick
+ REAL(wp), DIMENSION(N_bgc), INTENT(in) :: bgc_bottom
+ REAL(wp), DIMENSION(Nlayer+1,Nlayer+1), INTENT(in) :: fl_brine_bgc
+ REAL(wp), DIMENSION(Nlayer,N_bgc), INTENT(inout) :: bgc_abs
+ REAL(wp), DIMENSION(Nlayer,N_bgc) :: bgc_temp
+ REAL(wp), DIMENSION(Nlayer,N_bgc) :: bgc_br
+ REAL(wp), DIMENSION(N_bgc) :: flux
+ INTEGER ::k,i,j
+
+ bgc_temp = bgc_abs
+
+ !Calculating brine concentration of advecting brine, saturation limitations must be included manually!
+ DO k = 1,N_active
+ bgc_br(k,:) = bgc_abs(k,:)/(MAX(psi_l(k)*thick(k)*rho_l,0.000000000000001_wp))
+
+ !How to limit the first tracer to the oxygen saturation limit
+ !sat_O2 = func_sat_O2(T(k),S_abs(k)/m(k))
+ !bgc_br(k,1) = MIN(bgc_br(k,1),1.25*sat_O2)
+ END DO
+
+ !Internal flows
+ DO i = 1,N_active
+ DO j = 1,N_active
+ DO k = 1,N_bgc
+ IF (fl_brine_bgc(i,j)*bgc_br(i,k) > bgc_abs(i,k)/3._wp) THEN
+ END IF
+ flux(k) = MIN(fl_brine_bgc(i,j)*bgc_br(i,k),bgc_abs(i,k)/3._wp)
+ END DO
+ bgc_temp(i,:) = bgc_temp(i,:) - flux(:)
+ bgc_temp(j,:) = bgc_temp(j,:) + flux(:)
+
+ END DO
+ END DO
+
+ !Flows which leave the domain
+ DO i = 1,N_active
+ j = N_active+1
+ DO k = 1,N_bgc
+ flux(k) = MIN(fl_brine_bgc(i,j)*bgc_br(i,k),bgc_abs(i,k)/3._wp)
+ END DO
+ bgc_temp(i,:) = bgc_temp(i,:) - flux(:)
+ END DO
+
+ !Flows which enter the domain
+ i = N_active+1
+ DO j = 1,N_active
+ DO k = 1,N_bgc
+ flux(k) = fl_brine_bgc(i,j)*bgc_bottom(k)
+ END DO
+ bgc_temp(j,:) = bgc_temp(j,:) + flux(:)
+ END DO
+
+ bgc_abs = bgc_temp
+ END SUBROUTINE bgc_advection
+
+
+
+END MODULE mo_mass