SAMSIM
Functions/Subroutines
mo_thermo_functions Module Reference

Contains subroutines and functions related to multi-phase thermodynamics. More...

Functions/Subroutines

subroutine, public gett (H, S_bu, T_in, T, phi, k)
 Determines equilibrium Temperature of a layer for given S_bu and H as well as solid fraction. More...
 
subroutine, public expulsion (phi, thick, m, psi_s, psi_l, psi_g, V_ex)
 Determines Brine flux expelled from out of a layer due to freezing. More...
 
subroutine, public sub_fl_q (psi_s_1, psi_l_1, psi_g_1, thick_1, T_1, psi_s_2, psi_l_2, psi_g_2, thick_2, T_2, fl_Q)
 Determines conductive heat flux between two layers. More...
 
subroutine, public sub_fl_q_0 (psi_s, psi_l, psi_g, thick, T, T_bound, direct_flag, fl_Q)
 Determines conductive Heat flux between layer and boundary temperatures. More...
 
subroutine, public sub_fl_q_styropor (k_styropor, fl_Q)
 Niels, 2017 add: Determines conductive Heat flux below styropor cover. More...
 
real(wp) function, public func_s_br (T, S_bu)
 Computes salinity of brine pockets for given temperature in Celsius of mushy layer. More...
 
real(wp) function, public func_ddt_s_br (T)
 Computes temperature derivative of brine pocket salinity for given temperature in Celsius of mushy layer. More...
 

Detailed Description

Contains subroutines and functions related to multi-phase thermodynamics.

See the subroutine and function descriptions for details.

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/.

Revision History
Started by Philipp Griewank 2010-07-08
Add function for styropor cover by Niels Fuchs, MPIMET (2017-01-03)
Modified salinity functions by Philipp Griewank, Uni K (2018-08-01)

Function/Subroutine Documentation

subroutine, public mo_thermo_functions::expulsion ( real(wp), intent(in)  phi,
real(wp), intent(in)  thick,
real(wp), intent(inout)  m,
real(wp), intent(out)  psi_s,
real(wp), intent(out)  psi_l,
real(wp), intent(out)  psi_g,
real(wp), intent(out)  V_ex 
)

Determines Brine flux expelled from out of a layer due to freezing.

If the volume of ice and brine exceed the Volume of the layer brine is expelled. The volume of the ejected brine is calculated and exported. The volume fractions are also calculated.

Revision History
first version by Philipp Griewank, (2010-07-19)
changes to mass, Enthalpy and Salinity are now computed in subroutine mass_transfer by Philipp Griewank, (2010-08-24)
real(wp) function, public mo_thermo_functions::func_ddt_s_br ( real(wp), intent(in)  T)

Computes temperature derivative of brine pocket salinity for given temperature in Celsius of mushy layer.

Subroutine computes one ddT_S_br for one given T NaCl solutions and seawater produce slight variations. Which solution is used is specified by salt_flag. Based on Notz 2005 p. 36 ddT_S_br = c2+2*c3*T+3*c4*T^2

For T below -20 we simply use a linear extension based on "Composition of sea ice and its tensile strength". The actual salinity function is much more complicated and depends on the salt composition, but the linear fit is far better than using the polynomial fit.

The NaCL precipitates at -22, leading to ice/salt kristall mix below -22. Ideally the whole code would be modified to take the non-continous transition at -22 but given that there is currently little interest I can't be bothered to put in the effor.

Revision History
First version by Philipp Griewank (2010-07-13)
Added linear bit by Philipp Griewank (2018-07-22)
Parameters
[in]tTemperature in Celsius
Returns
derivative of Brine salinity
real(wp) function, public mo_thermo_functions::func_s_br ( real(wp), intent(in)  T,
real(wp), intent(in), optional  S_bu 
)

Computes salinity of brine pockets for given temperature in Celsius of mushy layer.

Subroutine computes one S_br for one given T in Celsius by third-order polynomial. NaCl solutions and seawater produce slight variations. Which solution is used is determined by salt_flag. S_br = c1+c2*T+c3*T^2+c4*T^3 Based on Notz 2005 p. 36

For T below -20 we simply use a linear extension based on "Composition of sea ice and its tensile strength". The actual salinity function is much more complicated and depends on the salt composition, but the linear fit is far better than using the polynomial fit.

The NaCL precipitates at -22, leading to ice/salt kristall mix below -22. Ideally the whole code would be modified to take the non-continous transition at -22 but given that there is currently little interest I can't be bothered to put in the effor.

Revision History
First version by Philipp Griewank (2010-07-12)
Changed to go through 0 by Philipp Griewank (2014-05-07)
Added linear bit by Philipp Griewank (2018-07-22)
Parameters
[in]tTemperature in Celsius
Returns
Brine salinity
subroutine, public mo_thermo_functions::gett ( real(wp), intent(in)  H,
real(wp), intent(in)  S_bu,
real(wp), intent(in)  T_in,
real(wp), intent(out)  T,
real(wp), intent(out)  phi,
integer, intent(in), optional  k 
)

Determines equilibrium Temperature of a layer for given S_bu and H as well as solid fraction.

The temperature of a fully liquid layer is used to see if the resulting brine salinity is lower than the bulk salinity. After checking if the layer is a fluid or a mushy layer the temperature is calculated by solving f(T) = 0 using the Newton method. f(T) = -latent_heat-H+latent_heat*S_bu/S_br(T) + c_s*T+c_s_beta*T^2/2 f'(T) = c_s+c_s_beta*T-latent_heat*S_bu*S_br'(T)/S_br^2 Described in Notz2005, subsubsection 5.6.1. See func_S_br(T) and func_ddT_S_br(T). First guess T_0 must be given, low first guess lead to overshooting which would lead to very high Temperatures. To avoid this, an if loop sets T to freezing T when T>0. Freezing T is also calculated at the beginning using the Newton-Method. If S_bu<0.001 then it is treated as pure ice.

Revision History
first version by Philipp Griewank (2010-07-13)
Freezing temperature is calculated and introduced if T goes above 0 by Philipp Griewank (2010-07-13)
Added if loops to deal with saltless ice by Philipp Griewank (2010-11-27)
Parameters
[in]hEnthalpy [J/kg]
[in]s_buBulk Salinity [g/kg]
[in]t_ininput Temperature for T_0 [C]
[out]tTemperature [C]
[out]phisolid fraction
subroutine, public mo_thermo_functions::sub_fl_q ( real(wp), intent(in)  psi_s_1,
real(wp), intent(in)  psi_l_1,
real(wp), intent(in)  psi_g_1,
real(wp), intent(in)  thick_1,
real(wp), intent(in)  T_1,
real(wp), intent(in)  psi_s_2,
real(wp), intent(in)  psi_l_2,
real(wp), intent(in)  psi_g_2,
real(wp), intent(in)  thick_2,
real(wp), intent(in)  T_2,
real(wp), intent(out)  fl_Q 
)

Determines conductive heat flux between two layers.

Details can be found in Notz 2005, especially equation 5.7. The gas volume is assumed to have no thermal properties at all. First the thermal resistance R is calculated using the approximated thermal conductivity of the mushy layer (see Notz 2005 eq. 3.41.). Then the heat flux Q is simply (T_1-T_2)/R "_1" denotes the upper layer and "_2" the lower layer. A positive heat flux is from lower to upper layer.

Revision History
First version by Philipp Griewank (2010-07-21)
subroutine, public mo_thermo_functions::sub_fl_q_0 ( real(wp), intent(in)  psi_s,
real(wp), intent(in)  psi_l,
real(wp), intent(in)  psi_g,
real(wp), intent(in)  thick,
real(wp), intent(in)  T,
real(wp), intent(in)  T_bound,
integer  direct_flag,
real(wp), intent(out)  fl_Q 
)

Determines conductive Heat flux between layer and boundary temperatures.

Details can be found in Notz 2005, especially equation 5.10 and 5.11. The gas volume is assumed to have no thermal properties. direct_flag denotes if the boundary layer is above or below the layer. 1 : = layer above boundary -1: = layer below boundary

Revision History
first version by Philipp Griewank (2010-07-21)
Parameters
[in]t_boundT_bound temperature of boundary layer
subroutine, public mo_thermo_functions::sub_fl_q_styropor ( real(wp), intent(in)  k_styropor,
real(wp), intent(inout)  fl_Q 
)

Niels, 2017 add: Determines conductive Heat flux below styropor cover.

Standard approach.

Revision History
first version by Niels Fuchs, MPIMET (2017-01-03)