def dens_hostetler(temp):
	#Water density
	#temp - water temperature, Celsius
	row0 = 1.e+3
	temp0 = 3.85
	const = 1.9549e-5
	const_power = 1.68
	return row0*(1.-const*abs(temp-temp0)**const_power)

def water_dens_ts(Temp,Sal):
	# The function WATER_DENS_TS
	# resturns the density of water, kg/m**3,
	# as a function of temperature and salinity
	# according to
	# (McCutcheon et al.,Water Quality and Maidment. 
	#  Handbook on Hydrology, 1993)
	
	# Input variables:
	# Temp --- water temperature, deg C
	# Sal  --- water salinity,    kg/kg
	A11 = 288.9414e0
	A12 = 508929.2e0
	A13 = 68.12963e0
	A14 = 3.9863e0
	alpha1 = 8.2449e-1
	alpha2 = 4.0899e-3
	alpha3 = 7.6438e-5
	alpha4 = 8.2467e-7
	alpha5 = 5.3675e-9
	beta1 = 5.724e-3
	beta2 = 1.0227e-4
	beta3 = 1.6546e-6
	gamma1 = 4.8314e-4

	row0 = 1.e+3

	# Converting salinity units from kg/kg to g/kg
	Sal_g_kg = Sal*1.e+3

	# The first term, dependent on temperature      
	water_dens_ts = row0 * \
 	( 1.-(Temp+A11)*(Temp-A14)**2/(A12*(Temp+A13) ) )
	
	A =  alpha1         - alpha2*Temp    + alpha3*Temp**2 - \
	alpha4*Temp**3 + alpha5*Temp**4
	
	B = - beta1          + beta2*Temp     - beta3*Temp**2
	
	C = gamma1

	# Adding the second term, dependent on temperature and salinity
	return water_dens_ts + A*Sal_g_kg + B*Sal_g_kg**1.5 + C*Sal_g_kg**2.


def bvf(temp1,temp2,sal1,sal2,dz):
	#Buoyancy frequency squared
	g=9.814
	row0=1e+3
	#return (max(g/row0*(dens_hostetler(temp2)-dens_hostetler(temp1))/dz,0.))**0.5
	return (max(g/row0*(water_dens_ts(temp2,sal2) - \
	water_dens_ts(temp1,sal1))/dz,0.))**0.5
	#return (max(g/row0*(temp2-temp1)*(-0.15)/dz,0.))**0.5

def mldepth(temp,z):
	#Mixed-layer depth
	n = len(temp)
	bvfarr=[0]*(n-1)
	for i in range(n-1):
		bvfarr[i] = bvf(temp[i],temp[i+1],z[i+1]-z[i])
	j = bvfarr.index(max(bvfarr))
	if j == 0:
		return float('nan')
	else:
		return z[j] + 0.5*(z[j+1] - z[j])

def maxbvf(temp,z):
	#Maximal bouyancy frequency
	n = len(temp)
	bvfarr=[0]*(n-1)
	for i in range(n-1):
		bvfarr[i] = bvf(temp[i],temp[i+1],z[i+1]-z[i])
	return max(bvfarr)

def thermthick(temp,z,thermt1,thermt2):
	#Thermocline thickness
	n = len(temp)
	#Upper boundary
	z1 = 0
	for i in range(n-1):
		if temp[i] > thermt1 and temp[i+1] < thermt1:
			dT = temp[i+1] - temp[i]
			dz = z[i+1] - z[i]
			z1 = z[i] + dz*(thermt1 - temp[i])/dT
	#Lower boundary
	z2 = z[n-1]
	for i in range(n-1):
		if temp[i] > thermt2 and temp[i+1] < thermt2:
			dT = temp[i+1] - temp[i]
			dz = z[i+1] - z[i]
			z2 = z[i] + dz*(thermt2 - temp[i])/dT
	if (z1 == 0 or z2 == 0):
		return float('nan')
	else:
		return z2 - z1

def radabs(S0,z):
        import numpy as np
        row0 = 1.e+3
        cp = 3.99E+3
        alpha0 = 0.5
        albedo = 0.06
        beta_abs = 0.35
        H = 2.#mixed-layer depth
        H_b = 7.#lake depth
        C = 1./(H_b - H)# coefficient for linear increase of ext. coef. in the pycnocline
        alphaint = alpha0*(z + 0.5*C*(z-H)**2)
        alphaz = alpha0*(1. + C*(z-H))
        return alphaz*S0*(1.-albedo)*(1.-beta_abs)/(row0*cp)*np.exp(-alphaint)


def KatoPhilDepth(time):
	#Mixed-layer depth evolution in Kato-Phillips experiment 
	#according to Price (1979), see Burchard, 2002, p.119
	#time, seconds
	g = 9.814; row0 = 1.E+3
	ddens_dtemp0 = -0.15 #-8.78E-2
	dtdz = -15./5.
	drodz = ddens_dtemp0*dtdz
	N0 = (g/row0 * drodz)**0.5
	u_s = 10**(-2.5)
	depth = 1.05*u_s*N0**(-0.5)*time**0.5
	return depth




