 MODULE WATER_DENSITY

 use LAKE_DATATYPES, only : ireals, iintegers
 use PHYS_CONSTANTS, only : &
 & row0


 implicit none

!The coefficients of UNESCO formula
!for water density dependence on temperature

 real(kind=ireals), parameter:: a0 = 800.969d-7
 real(kind=ireals), parameter:: a1 = 588.194d-7
 real(kind=ireals), parameter:: a2 = 811.465d-8
 real(kind=ireals), parameter:: a3 = 476.600d-10

!The coefficients of formula for
!water density dependence on temperature and salinity
!(McCutcheon et al.,Water Quality and Maidment. 
! Handbook on hydrology, 1993)
 
 real(kind=ireals), parameter:: A11 = 288.9414d0
 real(kind=ireals), parameter:: A12 = 508929.2d0
 real(kind=ireals), parameter:: A13 = 68.12963d0
 real(kind=ireals), parameter:: A14 = 3.9863d0
 real(kind=ireals), parameter:: alpha1 = 8.2449d-1
 real(kind=ireals), parameter:: alpha2 = 4.0899d-3
 real(kind=ireals), parameter:: alpha3 = 7.6438d-5
 real(kind=ireals), parameter:: alpha4 = 8.2467d-7
 real(kind=ireals), parameter:: alpha5 = 5.3675d-9
 real(kind=ireals), parameter:: beta1 = 5.724d-3
 real(kind=ireals), parameter:: beta2 = 1.0227d-4
 real(kind=ireals), parameter:: beta3 = 1.6546d-6
 real(kind=ireals), parameter:: gamma1 = 4.8314d-4

!The coefficients for "Kivu" EOS (Schmid et al., 2003)
 real(kind=ireals), parameter :: at1 = 999.843, at2 = 65.4891d-3, at3 = - 8.56272d-3, at4 = 0.059385d-3
 real(kind=ireals), parameter :: beta = 0.75d-3 ! kg/g
 real(kind=ireals), parameter :: beta_co2 = 0.284d-3 ! coefficient for co2, kg/g
 real(kind=ireals), parameter :: beta_ch4 = -1.25d-3 ! coefficient for ch4, kg/g
 real(kind=ireals), parameter :: co2_ref = 0.d0, ch4_ref = 0.d0 ! reference concentrations, kg/g

 real(kind=ireals), save ::  DDENS_DTEMP0, DDENS_DSAL0, DENS_REF

 integer(kind=iintegers), parameter :: eos_Hostetler = 1
 integer(kind=iintegers), parameter :: eos_TEOS = 2
 integer(kind=iintegers), parameter :: eos_Kivu = 3
 integer(kind=iintegers), parameter :: eos_UNESCO = 4
 integer(kind=iintegers), parameter :: eos_McCatcheon = 5

 !Temperature and salinity reference values
 real(kind=ireals), parameter :: temp_ref = 15. ! deg. C
 real(kind=ireals), parameter :: sal_ref = 0. !kg/kg


 contains

 SUBROUTINE SET_DENS_CONSTS(eos)

 implicit none 

 integer(kind=iintegers), intent(in) :: eos
 

 select case (eos)
 case(eos_Hostetler)
   call DDENS_HOSTETLER(temp_ref,DDENS_DTEMP0,DDENS_DSAL0)
   DENS_REF = DENS_HOSTETLER(temp_ref)
 case(eos_TEOS)
   print*, 'Not operational eos: STOP '
   STOP
 case(eos_Kivu)
   call DDENS_KIVU(temp_ref,sal_ref,DDENS_DTEMP0,DDENS_DSAL0)
   DENS_REF = WATER_DENS_TS_KIVU(temp_ref,sal_ref)
 case(eos_UNESCO)
   call DDENS_UNESCO(temp_ref,DDENS_DTEMP0,DDENS_DSAL0)
   DENS_REF = WATER_DENS_T_UNESCO(temp_ref)
 case(eos_McCatcheon)
   DDENS_DTEMP0 = DDENS_DTEMP(temp_ref,sal_ref) 
   DDENS_DSAL0  = DDENS_DSAL (temp_ref,sal_ref)
   DENS_REF = WATER_DENS_TS(temp_ref,sal_ref)
 end select

 END SUBROUTINE SET_DENS_CONSTS


 FUNCTION SET_DDENS_DTEMP(eos,lindens,Temp,Sal)

 implicit none 

 ! Input/output variables
 integer(kind=iintegers), intent(in) :: eos !> Equation of state identifier
 integer(kind=iintegers), intent(in) :: lindens !> Switch for linear density function
 real(kind=ireals), intent(in) :: Temp, Sal

 real(kind=ireals) :: SET_DDENS_DTEMP

 ! Local variables
 real(kind=ireals) :: x

 if (lindens == 0) then
   select case (eos)
   case(eos_Hostetler)
     call DDENS_HOSTETLER(Temp,SET_DDENS_DTEMP,x)
   case(eos_TEOS)
     print*, 'Not operational eos: STOP '
     STOP
   case(eos_Kivu)
     call DDENS_KIVU(Temp,Sal,SET_DDENS_DTEMP,x)
   case(eos_UNESCO)
     call DDENS_UNESCO(Temp,SET_DDENS_DTEMP,x)
   case(eos_McCatcheon)
     SET_DDENS_DTEMP = DDENS_DTEMP(Temp,Sal) 
   end select
 else
   SET_DDENS_DTEMP = DDENS_DTEMP0
 endif

 END FUNCTION SET_DDENS_DTEMP


 FUNCTION SET_DDENS_DSAL(eos,lindens,Temp,Sal)

 implicit none 

 ! Input/output variables
 integer(kind=iintegers), intent(in) :: eos !> Equation of state identifier
 integer(kind=iintegers), intent(in) :: lindens !> Switch for linear density function
 real(kind=ireals), intent(in) :: Temp, Sal

 real(kind=ireals) :: SET_DDENS_DSAL

 ! Local variables
 real(kind=ireals) :: x

 if (lindens == 0) then
   select case (eos)
   case(eos_Hostetler)
     call DDENS_HOSTETLER(Temp,x,SET_DDENS_DSAL)
   case(eos_TEOS)
     print*, 'Not operational eos: STOP '
     STOP
   case(eos_Kivu)
     call DDENS_KIVU(Temp,Sal,x,SET_DDENS_DSAL)
   case(eos_UNESCO)
     call DDENS_UNESCO(Temp,x,SET_DDENS_DSAL)
   case(eos_McCatcheon)
     SET_DDENS_DSAL = DDENS_DSAL(Temp,Sal) 
   end select
 else
   SET_DDENS_DSAL = DDENS_DSAL0
 endif

 END FUNCTION SET_DDENS_DSAL


 !>Subroutine DENSITY_W calculates density profile
 SUBROUTINE DENSITY_W(M,eos,lindens,Tw,Sal,preswat,row)

 integer(kind=iintegers), intent(in) :: eos, M, lindens

 real(kind=ireals), intent(in) :: Tw(1:M+1), Sal(1:M+1), preswat(1:M+1)

 real(kind=ireals), intent(out) :: row(1:M+1)

 integer(kind=iintegers) :: i

 !Calculation of water density profile at the current timestep
 if (lindens == 0) then
   select case (eos)
   case(eos_Hostetler)
     do i = 1, M+1
       row(i) = DENS_HOSTETLER(Tw(i))
     enddo
   case(eos_TEOS)
     do i = 1, M+1
       row(i) = GSW_RHO(Sal(i),Tw(i),preswat(i))
     enddo
   case(eos_Kivu)
     do i = 1, M+1
       row(i) = WATER_DENS_TS_KIVU(Tw(i),Sal(i)) !0.d0
     enddo
   case(eos_UNESCO)
     do i = 1, M+1
       row(i) = WATER_DENS_T_UNESCO(Tw(i)) !0.d0
     enddo
   case(eos_McCatcheon)
     do i = 1, M+1
       row(i) = WATER_DENS_TS(Tw(i),Sal(i)) ! McCutcheon et al., 1993
     enddo
   end select
 else
   do i = 1, M+1
     row(i) = DENS_REF + DDENS_DTEMP0*(Tw(i) - temp_ref) + &
     &                   DDENS_DSAL0 *(Sal(i) - sal_ref) 
   enddo
 endif

 END SUBROUTINE DENSITY_W


 real(kind=ireals) FUNCTION WATER_DENS_T_UNESCO(Temp)

!The function WATER_DENS_T_UNESCO
!returns the density of water, kg/m**3, 
!as a function of temperature 
!according to UNESCO formula
 
 implicit none

 real(kind=ireals), intent(in):: Temp
 
 WATER_DENS_T_UNESCO = &
&  row0*(1+a0+a1*Temp-a2*Temp**2+a3*Temp**3)
 END FUNCTION WATER_DENS_T_UNESCO


 SUBROUTINE DDENS_UNESCO(Temp,ddens_dtemp,ddens_dsal)
 implicit none

 real(kind=ireals), intent(in):: Temp
 real(kind=ireals), intent(out):: ddens_dtemp, ddens_dsal
 
 ddens_dtemp = &
 & row0*(a1 - 2.*a2*Temp + 3.*a3*Temp**2)
 ddens_dsal = 0.

 END SUBROUTINE DDENS_UNESCO


 FUNCTION WATER_DENS_TS_KIVU(Temp,Sal)

! The function WATER_DENS_TS_KIVU
! resturns the density of water, kg/m**3,
! as a function of temperature and salinity
! adjusted for Lake Kivu (Rwanda & DRC) according to
! (Schmid et al. 2003)

! Input variables:
 
! Temp --- water temperature, deg C
! Sal  --- water salinity,    kg/kg
  
  implicit none

  real(kind=ireals), intent(in):: Temp
  real(kind=ireals), intent(in):: Sal

  real(kind=ireals) :: WATER_DENS_TS_KIVU

  real(kind=ireals) :: Sal_g_kg,A,B,C

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

  WATER_DENS_TS_KIVU = (at1 + at2*Temp + at3*Temp**2 + at4*Temp**3)* &
  & (1. + beta*Sal_g_kg + beta_co2*co2_ref + beta_ch4*ch4_ref)

  END FUNCTION WATER_DENS_TS_KIVU


  SUBROUTINE DDENS_KIVU(Temp,Sal,ddens_dtemp,ddens_dsal)
  implicit none

  real(kind=ireals), intent(in) :: Temp, Sal
  real(kind=ireals), intent(out) :: ddens_dtemp, ddens_dsal

  real(kind=ireals) :: Sal_g_kg

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

  ddens_dtemp = (at2 + 2.*at3*Temp + 3.*at4*Temp**2)* &
  & (1. + beta*Sal_g_kg + beta_co2*co2_ref + beta_ch4*ch4_ref)
  ddens_dsal = (at1 + at2*Temp + at3*Temp**2 + at4*Temp**3)* &
  & beta*1.d+3

  END SUBROUTINE DDENS_KIVU


  real(kind=ireals) FUNCTION 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
  
  implicit none

  real(kind=ireals), intent(in):: Temp
  real(kind=ireals), intent(in):: Sal

  real(kind=ireals) Sal_g_kg,A,B,C

! Converting salinity units from kg/kg to g/kg
  Sal_g_kg = Sal*1.d+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
  WATER_DENS_TS = WATER_DENS_TS + &
 & A*Sal_g_kg + B*Sal_g_kg**1.5 + C*Sal_g_kg**2.

  END FUNCTION WATER_DENS_TS


  real(kind=ireals) FUNCTION DDENS_DTEMP(Temp,Sal)

! The function DDENS_DTEMP
! returns the derivative of water density
! on temperature, kg/(m**3*C)
  
! Input variables:
 
! Temp --- water temperature, deg C
! Sal  --- water salinity,    kg/kg

  implicit none

  real(kind=ireals), intent(in):: Temp
  real(kind=ireals), intent(in):: Sal      

  real(kind=ireals) Sal_g_kg,A,B,C

  DDENS_DTEMP = &
  & -(row0/A12)*(Temp-A14)/( (Temp+A13)**2 ) * &
  & ( (Temp+A13)*( (Temp-A14)+2.*(Temp+A11) ) - &
  &  (Temp+A11)*(Temp-A14) )

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

  A = - alpha2 + 2.*alpha3*Temp - 3.*alpha4*Temp**2 + 4.*alpha5*Temp**3

  B = + beta2  - 2.*beta3*Temp

  DDENS_DTEMP = DDENS_DTEMP + A*Sal_g_kg + B*Sal_g_kg**1.5

  END FUNCTION DDENS_DTEMP
  

  real(kind=ireals) FUNCTION DDENS_DSAL(Temp,Sal)

! The function DDENS_DSAL
! returns the derivative of water density
! on salinity, , kg/(m**3*kg/kg)
  
! Input variables:
 
! Temp --- water temperature, deg C
! Sal  --- water salinity,    kg/kg

  implicit none

  real(kind=ireals), intent(in):: Temp
  real(kind=ireals), intent(in):: Sal

  real(kind=ireals) Sal_g_kg,A,B,C

  A =   alpha1         - alpha2*Temp    + alpha3*Temp**2 &
  &   - alpha4*Temp**3 + alpha5*Temp**4

  B = - beta1          + beta2*Temp     - beta3*Temp**2

  C = gamma1

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

  DDENS_DSAL = &
 & A + 1.5*B*Sal_g_kg**0.5 + 2.*C*Sal_g_kg

  DDENS_DSAL = DDENS_DSAL*1.d+3

  END FUNCTION DDENS_DSAL
  
  
  FUNCTION DENS_HOSTETLER(temp)
  
! Function DENS_HOSTETLER calculates the density of water
! dependent on temperature only, following Hostetler and Bartlein, 1990

  implicit none

  real(kind=ireals) :: DENS_HOSTETLER
  
  real(kind=ireals), intent(in) :: temp
  
! Parameters      
  real(kind=ireals), parameter :: row0 = 1.d+3
  real(kind=ireals), parameter :: temp0 = 3.85
  real(kind=ireals), parameter :: const = 1.9549d-5
  real(kind=ireals), parameter :: const_power = 1.68
  
  DENS_HOSTETLER = row0*(1.-const*abs(temp-temp0)**const_power)
  
  END FUNCTION DENS_HOSTETLER


  SUBROUTINE DDENS_HOSTETLER(Temp,ddens_dtemp,ddens_dsal)
  implicit none

! Parameters      
  real(kind=ireals), parameter :: row0 = 1.d+3
  real(kind=ireals), parameter :: temp0 = 3.85
  real(kind=ireals), parameter :: const = 1.9549d-5
  real(kind=ireals), parameter :: const_power = 1.68

  real(kind=ireals), intent(in):: Temp
  real(kind=ireals), intent(out):: ddens_dtemp, ddens_dsal
  
  ddens_dtemp = row0*( - const_power*const*abs(temp-temp0)**(const_power-1)* &
  & sign(1._ireals,temp-temp0))
  ddens_dsal = 0.

  END SUBROUTINE DDENS_HOSTETLER


!--------------------------------------------------------------------------

!--------------------------------------------------------------------------
! density and enthalpy, based on the 48-term expression for density from TEOS-2010
!--------------------------------------------------------------------------

!==========================================================================
FUNCTION GSW_RHO(sa_,ct_,p_) 
!==========================================================================

!  Calculates in-situ density from Absolute Salinity and Conservative 
!  Temperature, using the computationally-efficient 48-term expression for
!  density in terms of SA, CT and p (McDougall et al., 2011).
!
! sa     : Absolute Salinity                               [kg/kg]
! ct     : Conservative Temperature                        [deg C]
! p      : sea pressure                                    [Pa]
! 
! gsw_rho  : in-situ density (48 term equation)

implicit none

integer, parameter :: r14 = selected_real_kind(14,30)

real (r14), parameter :: v01 =  9.998420897506056d2, v02 =  2.839940833161907d0
real (r14), parameter :: v03 = -3.147759265588511d-2, v04 =  1.181805545074306d-3
real (r14), parameter :: v05 = -6.698001071123802d0, v06 = -2.986498947203215d-2
real (r14), parameter :: v07 =  2.327859407479162d-4, v08 = -3.988822378968490d-2
real (r14), parameter :: v09 =  5.095422573880500d-4, v10 = -1.426984671633621d-5
real (r14), parameter :: v11 =  1.645039373682922d-7, v12 = -2.233269627352527d-2
real (r14), parameter :: v13 = -3.436090079851880d-4, v14 =  3.726050720345733d-6
real (r14), parameter :: v15 = -1.806789763745328d-4, v16 =  6.876837219536232d-7
real (r14), parameter :: v17 = -3.087032500374211d-7, v18 = -1.988366587925593d-8
real (r14), parameter :: v19 = -1.061519070296458d-11, v20 =  1.550932729220080d-10
real (r14), parameter :: v21 =  1.0d0, v22 =  2.775927747785646d-3, v23 = -2.349607444135925d-5
real (r14), parameter :: v24 =  1.119513357486743d-6, v25 =  6.743689325042773d-10
real (r14), parameter :: v26 = -7.521448093615448d-3, v27 = -2.764306979894411d-5
real (r14), parameter :: v28 =  1.262937315098546d-7, v29 =  9.527875081696435d-10
real (r14), parameter :: v30 = -1.811147201949891d-11, v31 = -3.303308871386421d-5
real (r14), parameter :: v32 =  3.801564588876298d-7, v33 = -7.672876869259043d-9
real (r14), parameter :: v34 = -4.634182341116144d-11, v35 =  2.681097235569143d-12
real (r14), parameter :: v36 =  5.419326551148740d-6, v37 = -2.742185394906099d-5
real (r14), parameter :: v38 = -3.212746477974189d-7, v39 =  3.191413910561627d-9
real (r14), parameter :: v40 = -1.931012931541776d-12, v41 = -1.105097577149576d-7
real (r14), parameter :: v42 =  6.211426728363857d-10, v43 = -1.119011592875110d-10
real (r14), parameter :: v44 = -1.941660213148725d-11, v45 = -1.864826425365600d-14
real (r14), parameter :: v46 =  1.119522344879478d-14, v47 = -1.200507748551599d-15
real (r14), parameter :: v48 =  6.057902487546866d-17 

real (kind=ireals), intent(in) :: sa_, ct_, p_
real (r14) :: sa, ct, p, sqrtsa, v_hat_denominator, &
& v_hat_numerator, gsw_rho

sa = sa_*1.d+3 ! Converting kg/kg to g/kg
p = p_*1.d-4   ! Converting pressure from Pa to dbar
ct = ct_

sqrtsa = sqrt(sa)

v_hat_denominator = v01 + ct*(v02 + ct*(v03 + v04*ct))  &
         + sa*(v05 + ct*(v06 + v07*ct) &
         + sqrtsa*(v08 + ct*(v09 + ct*(v10 + v11*ct)))) &
         + p*(v12 + ct*(v13 + v14*ct) + sa*(v15 + v16*ct) &
         + p*(v17 + ct*(v18 + v19*ct) + v20*sa))

v_hat_numerator = v21 + ct*(v22 + ct*(v23 + ct*(v24 + v25*ct))) &
       + sa*(v26 + ct*(v27 + ct*(v28 + ct*(v29 + v30*ct))) + v36*sa &
       + sqrtsa*(v31 + ct*(v32 + ct*(v33 + ct*(v34 + v35*ct)))))  &
       + p*(v37 + ct*(v38 + ct*(v39 + v40*ct))  &
       + sa*(v41 + v42*ct) &
       + p*(v43 + ct*(v44 + v45*ct + v46*sa) &
       + p*(v47 + v48*ct)))

GSW_RHO = v_hat_denominator/v_hat_numerator

return
END FUNCTION GSW_RHO


  END MODULE WATER_DENSITY


  MODULE PHYS_FUNC

  use LAKE_DATATYPES, only : ireals, iintegers

  contains
  real(kind=ireals) FUNCTION TURB_DENS_FLUX(tempflux,salflux,Temp,Sal)

! Function TURB_DENS_FLUX            _____
! returns the turbulent density flux w'ro' in water
! at given temperature, 
!          salinity,         ____
!          temperature flux  w'T'
!                            ____
!      and salinity    flux  w's'

! Input variables:
! tempflux --- kinematic heat flux, m*C/s
! salflux  --- salinity       flux, m*(kg/kg)/s
! Temp     --- temperature        , C
! Sal      --- salinity           , kg/kg

  use water_density, only : &
 & ddens_dtemp, &
 & ddens_dsal
  
  implicit none
  
  real(kind=ireals), intent(in):: tempflux
  real(kind=ireals), intent(in):: salflux
  real(kind=ireals), intent(in):: Temp
  real(kind=ireals), intent(in):: Sal

  TURB_DENS_FLUX = &
 & ddens_dtemp(Temp,Sal)*tempflux + &
 & ddens_dsal (Temp,Sal)*salflux

  END FUNCTION TURB_DENS_FLUX

 
!  real(kind=ireals) FUNCTION dirdif()
!  use atmos, only:
! & shortwave
! implicit none
!  real(kind=ireals) cloud
!  real(kind=ireals) dirdif0,b,S0,sinh0

!  SAVE

!  b=1./3.
!  S0=1367.
!  cloud = 0.
 
!  dirdif0 = (shortwave-b*S0*sinh0())/(b*(S0*sinh0()-shortwave))
!  dirdif = dirdif0*(1.-cloud)
!  dirdif = max(dirdif,0.d0)
  
!  END FUNCTION dirdif


  FUNCTION SINH0(year,month,day,hour,phi)

! SINH0 is sine of solar angle ( = cosine of zenith angle)   
  use TIME_LAKE_MOD, only : JULIAN_DAY  
  implicit none
  
  real(kind=ireals) :: sinh0
  
  integer(kind=iintegers), intent(in) :: year
  integer(kind=iintegers), intent(in) :: month
  integer(kind=iintegers), intent(in) :: day

  real(kind=ireals)   , intent(in) :: hour
  real(kind=ireals)   , intent(in) :: phi

  real(kind=ireals) delta
  real(kind=ireals) theta
  real(kind=ireals) pi
   real(kind=ireals) phi_rad

  integer(kind=iintegers) nday
 
  pi=4.*datan(1.d0)

  nday  = JULIAN_DAY(year,month,day)

  delta = 23.5d0*pi/180.d0*cos(2*pi*(float(nday)-173.d0)/365.d0)
  theta = pi*(hour-12.d0)/12.d0
  phi_rad = phi*pi/180.d0
  sinh0 = sin(phi_rad)*sin(delta) + &
  & cos(phi_rad)*cos(delta)*cos(theta)
  sinh0=max(sinh0,0.d0) 

  END FUNCTION SINH0


  !> Function LONGWAVE_RADIATION returns incident 
  !! longwave (atmospheric) radiation on the earth surface
  FUNCTION LONGWAVE_RADIATION(tempair,humair,pressure,cloud_sky_fraction) result(LR)
  use UNITS, only : Pa_to_hPa
  use PHYS_CONSTANTS, only : Kelvin0, Rd_d_Rwv, sigma
  implicit none
  real(kind=ireals), intent(in) :: tempair !> Surface air temperature, Celsius
  real(kind=ireals), intent(in) :: humair  !> Surface air specific humidity, kg/kg
  real(kind=ireals), intent(in) :: pressure!> Surface atmospheric pressure, Pa
  real(kind=ireals), intent(in) :: cloud_sky_fraction !> Cloud-sky fraction, 0-1, n/d
  !> Output variables
  real(kind=ireals) :: LR
  !> Local variables
  real(kind=ireals) :: emis, a1, a2, a3, parpres, Rwv_d_Rd

  integer(kind=iintegers), parameter :: CAN = 1 !(Anstrom, 1915; Marthews et al., Theor Appl Climatol, 2012)
  integer(kind=iintegers), parameter :: CBR = 2 !(Brunt, 1932; Marthews et al., Theor Appl Climatol, 2012)
  integer(kind=iintegers), parameter :: CSW = 3 !(Swinbank, 1963; Marthews et al., Theor Appl Climatol, 2012)
  integer(kind=iintegers), parameter :: CIJ = 4 !(Idso and Jackson, 1969; Marthews et al., Theor Appl Climatol, 2012)
  integer(kind=iintegers), parameter :: CBT = 5 !(Brutsaert, 1975, 1982; Marthews et al., Theor Appl Climatol, 2012)
  integer(kind=iintegers), parameter :: CID = 6 !(Idso, 1981; Marthews et al., Theor Appl Climatol, 2012)
  integer(kind=iintegers), parameter :: CKZ = 7 !(Konzelmann et al., 1994; Marthews et al., Theor Appl Climatol, 2012)
  integer(kind=iintegers), parameter :: nCSE = CAN !switch for clear-sky emissivity parameterization

  Rwv_d_Rd = 1./Rd_d_Rwv
  parpres = humair*pressure*Rwv_d_Rd/(1. + humair*(Rwv_d_Rd - 1.))*Pa_to_hPa

  select case(nCSE)
  case(CAN)
    a1 = 0.83
    a2 = 0.18
    a3 = 0.067
    emis = a1 - a2*10**( - a3*parpres)
  case(CBR)
    a1 = 0.51
    a2 = 0.066
    emis = a1 + a2*sqrt(parpres)
  case(CSW)
    a1 = 0.0000092
    emis = a1*(tempair + Kelvin0)**2
  case(CIJ)
    a1 = 0.261
    a2 = 0.000777
    emis = 1. - a1*exp( - a2*tempair**2)
  case(CBT)
    a1 = 1.24
    emis = a1*(parpres/(tempair + Kelvin0))**(1./7.)
  case(CID)
    a1 = 0.7
    a2 = 0.0000595
    a3 = 1500.
    emis = a1 + a2*parpres*exp(a3/(tempair + Kelvin0))
  case(CKZ)
    a1 = 0.23
    a2 = 0.484
    emis = a1 + a2*(parpres/(tempair + Kelvin0))**(1./8.)
  end select

  !Correction of emissivity for cloudiness
  a1 = 0.22
  emis = emis*(1. + a1*cloud_sky_fraction)
  emis = min(emis,1._ireals)
  LR = emis*sigma*(tempair + Kelvin0)**4

  END FUNCTION LONGWAVE_RADIATION

  !> Function SHORTWAVE_RADIATION returns incident 
  !! shortwave (solar) radiation (direct+diffuse) on the earth surface
  FUNCTION SHORTWAVE_RADIATION(time_group,phi,cloud_sky_fraction) result(SR)
  use PHYS_CONSTANTS, only : solar_constant
  use ARRAYS, only : time_type
  implicit none
  !> Input variables
  type(time_type), intent(in) :: time_group !>The group of time/date variables, time assumed local, not GMT
  real(kind=ireals), intent(in) :: phi !>The latitude, deg.
  real(kind=ireals), intent(in) :: cloud_sky_fraction !> Cloud-sky fraction, 0-1, n/d
  !> Output variables
  real(kind=ireals) :: SR
  !> Local variables
  real(kind=ireals) :: sinsolh, trans_atm, SR_dir, SR_dif, a1, a2, SR_atten

  sinsolh = SINH0(time_group%year,time_group%month,time_group%day,time_group%hour,phi)
  trans_atm = 0.6 ! atmospheric transmissivity for direct solar radiation
  SR_atten = trans_atm**(1./sinsolh)
  !Direct solar radiation at the horizontal surface (Burger's law)
  SR_dir = solar_constant*SR_atten*sinsolh
  !Diffuse solar radiation at the horizontal surface (Liu and Jordan, Solar Energy, 1960; 
  !Jemaa et al., Energy Procedia, 2013)
  a1 = 0.271
  a2 = 0.294
  SR_dif = solar_constant*(a1 - a2*SR_atten)*sinsolh
  SR = SR_dir + SR_dif
  !Correction of global radiation by cloudiness (Kasten and Czeplak, Solar Energy, 1980; 
  !Jemaa et al., Energy Procedia, 2013)
  a1 = 0.75
  a2 = 3.4
  SR = SR*(1.-a1*cloud_sky_fraction**a2)

  END FUNCTION SHORTWAVE_RADIATION


  real(kind=ireals) FUNCTION QS(phase,t,p)
  
! QS - specific humidity, kg/kg, for saturated water vapour

  implicit none
  real(kind=ireals) t,p,aw,bw,ai,bi,a,b,es
  integer(kind=iintegers) phase
!phase = 1 - liquid, 0 - ice
aw = 7.6326    
bw = 241.9    
ai = 9.5  
bi = 265.5 
a  = phase*aw+(1-phase)*ai
b  = phase*bw+(1-phase)*bi
es=610.7*10.**(a*t/(b+t))
QS=0.622*es/p
END FUNCTION

!> Derivative of saturated specific humidity on temperature
FUNCTION DQSDT(phase,t,p)
! QS - specific humidity, kg/kg, for saturated water vapour
implicit none
real(kind=ireals) :: t,p,aw,bw,ai,bi,a,b,qs,es
real(kind=ireals) :: DQSDT
integer(kind=iintegers) :: phase
!phase = 1 - liquid, 0 - ice
aw = 7.6326    
bw = 241.9    
ai = 9.5  
bi = 265.5 
a  = phase*aw+(1-phase)*ai
b  = phase*bw+(1-phase)*bi
es = 610.7*10.**(a*t/(b+t))
qs = 0.622*es/p
DQSDT = qs*a*b*log(10.)/(b+t)**2
END FUNCTION DQSDT


real(kind=ireals) FUNCTION ES(phase,t,p)

!     ES is pressure of saturated water vapour, Pa    

implicit none
real(kind=ireals) t,p,aw,bw,ai,bi,a,b
integer(kind=iintegers) phase
!phase = 1 - liquid, 0 - ice
aw = 7.6326    
bw = 241.9    
ai = 9.5  
bi = 265.5 
a  = phase*aw+(1-phase)*ai
b  = phase*bw+(1-phase)*bi
ES = 610.7*10.**(a*t/(b+t))
END FUNCTION


real(kind=ireals) FUNCTION MELTPNT(C,p,npar)
implicit none

real(kind=ireals), intent(in) :: C ! salinity, kg/kg
real(kind=ireals), intent(in) :: p ! pressure, Pa
integer(kind=iintegers), intent(in) :: npar 

real(kind=ireals), parameter :: dtdc = 66.7, pnorm = 101325
real(kind=ireals), parameter :: Pa2dbar = 1.d-4

select case(npar)
case(0)
  MELTPNT = 0.
case(1)
  MELTPNT = 0. - C*dtdc
case(2)
  ! Jackett et al. 2004 formula
  MELTPNT = FP_THETA(C*1.E+3_ireals, (p-pnorm)*Pa2dbar,'air-sat') ! convertion to 0/00
end select

END FUNCTION MELTPNT

FUNCTION FP_THETA(s,p,sat)

!   potential temperature freezing point of seawater, as in
!   Jackett, McDougall, Feistel, Wright and Griffies (2004), submitted JAOT
!
!   s                : salinity                               (psu)
!   p                : gauge pressure                         (dbar)
!                      (absolute pressure - 10.1325 dbar)
!   sat              : string variable
!                       'air-sat'  - air saturated water       
!                       'air-free' - air free water
!
!   fp_theta         : potential freezing temperature         (deg C, ITS-90)
!
!   check value      : fp_theta(35,200,'air-sat')   = -2.076426227617581 deg C
!                      fp_theta(35,200,'air-free') = -2.074408175943127 deg C
!
!   DRJ on 2/6/04


implicit real(kind=ireals)(a-h,o-z)

character*(*) sat


sqrts = sqrt(s)

tf_num =                          2.5180516744541290d-03    +     &
                              s*(-5.8545863698926184d-02    +     &
                          sqrts*( 2.2979985780124325d-03    -     &
                         sqrts *  3.0086338218235500d-04))  +     &
                              p*(-7.0023530029351803d-04    +     &
                              p*( 8.4149607219833806d-09    +     &
                             s *  1.1845857563107403d-11));

tf_den =                          1.0000000000000000d+00    +     &
                              p*(-3.8493266309172074d-05    +     &
                             p *  9.1686537446749641d-10)   +     &
                      s*s*sqrts*  1.3632481944285909d-06 


fp_theta = tf_num/tf_den;


if(sat.eq.'air-sat') then 
    fp_theta = fp_theta          -2.5180516744541290d-03    +     &
                             s *  1.428571428571429d-05
elseif(sat.eq.'air-free') then  
    continue
else  
    print *, '***   Error in fp_theta.f90: invalid third argument   ***'
    print *
    stop
endif


return
END FUNCTION FP_THETA
  
  FUNCTION WATER_FREEZE_MELT(temperature, grid_size, &
 & melting_point, switch)
 
! The function WATER_FREEZE_MELT checks, if the
! heat storage of a grid cell is larger or equal
! the latent heat, necessary for phase transition
  
  use PHYS_CONSTANTS, only : &
 & cw_m_row0, &
 & ci_m_roi, &
 & row0_m_Lwi, &
 & roi_m_Lwi
  use NUMERIC_PARAMS, only : &
 & min_ice_thick, &
 & min_water_thick, &
 & T_phase_threshold
 
  implicit none
  
! Input variables
  real(kind=ireals), intent(in) :: temperature ! Temperature, Celsius
  real(kind=ireals), intent(in) :: grid_size
  real(kind=ireals), intent(in) :: melting_point
  
  integer(kind=iintegers), intent(in) :: switch ! +1 for water -> ice transition
                                   ! -1 for ice -> water transition
  logical :: WATER_FREEZE_MELT
  
  if (switch == +1) then
    if ( ( melting_point - T_phase_threshold - temperature) * &
 &   cw_m_row0*grid_size > min_ice_thick*roi_m_Lwi .or. &
 &   ( melting_point - T_phase_threshold - temperature) > 0.2d0) &
 &   then
      WATER_FREEZE_MELT = .true.
    else
      WATER_FREEZE_MELT = .false.
    endif
  elseif (switch == -1) then
    if ( (temperature - T_phase_threshold - melting_point) * &
 &   ci_m_roi*grid_size > min_water_thick*row0_m_Lwi .or. &
 &   (temperature - T_phase_threshold - melting_point) > 0.2d0) &
 &   then
      WATER_FREEZE_MELT = .true.
    else
      WATER_FREEZE_MELT = .false.
    endif     
  else
    write(*,*) 'Wrong switch in WATER_FREEZE_MELT: STOP'
    STOP
  endif
  
  END FUNCTION WATER_FREEZE_MELT


  !> Subroutine TURB_SCALES calculates turbulent scales, both bulk and local
  SUBROUTINE TURB_SCALES(gsp, ls, bathymwater, wst, RadWater, &
  & k_turb_T_flux, T_massflux, row, eflux0_kinem, &
  & turb_density_flux, Buoyancy0, tau, kor, i_maxN, H_mixed_layer, maxN, w_conv_scale, &
  & T_conv_scale, Wedderburn, LakeNumber, Rossby_rad, ThermThick, ReTherm, RiTherm, trb)
 
  use DRIVING_PARAMS, only : M, eos, lindens
  use PHYS_CONSTANTS, only : g, row0, cw_m_row0, g_d_Kelvin0, niu_wat
  use ARRAYS_BATHYM, only : bathym, layers_type
  use ARRAYS_GRID, only : gridspacing_type
  use WATER_DENSITY, only : DDENS_DTEMP0, DDENS_DSAL0
  use NUMERIC_PARAMS, only : small_value
  use ARRAYS_WATERSTATE, only : waterstate_type
  use RADIATION, only : rad_type
  use ARRAYS_TURB, only : turb_type
  use WATER_DENSITY, only : SET_DDENS_DTEMP, SET_DDENS_DSAL

  implicit none

! Input variables
  type(gridspacing_type), intent(in) :: gsp
  type(layers_type), intent(in) :: ls
  type(bathym), intent(in) :: bathymwater(1:M+1)
  type(waterstate_type), intent(in) :: wst
  type(rad_type), intent(in) :: RadWater

  real(kind=ireals), intent(in) :: k_turb_T_flux(1:M)
  real(kind=ireals), intent(in) :: T_massflux(1:M)
  real(kind=ireals), intent(in) :: row(1:M+1)
  real(kind=ireals), intent(in) :: eflux0_kinem
  real(kind=ireals), intent(in) :: turb_density_flux
  real(kind=ireals), intent(in) :: Buoyancy0
  real(kind=ireals), intent(in) :: tau
  real(kind=ireals), intent(in) :: kor

! Output variables
  integer(kind=iintegers), intent(out) :: i_maxN

  real(kind=ireals), intent(out) :: H_mixed_layer
  real(kind=ireals), intent(out) :: maxN
  real(kind=ireals), intent(out) :: w_conv_scale
  real(kind=ireals), intent(out) :: T_conv_scale
  real(kind=ireals), intent(out) :: Wedderburn !Wedderburn number
  real(kind=ireals), intent(out) :: LakeNumber
  real(kind=ireals), intent(out) :: Rossby_rad !Internal Rossby radius
  real(kind=ireals), intent(out) :: ThermThick !Thermocline thickness
  real(kind=ireals), intent(out) :: ReTherm    !Reynolds number in thermocline
  real(kind=ireals), intent(out) :: RiTherm    !Richardson number in thermocline

  type(turb_type), intent(inout) :: trb
  
! External functions
  real(kind=ireals), external:: DZETA
  
  real(kind=ireals), parameter :: Ttherm0 = 14., Ttherm1 = 8. !temperatures of top and bottom of thermocline, C

! Auxilliary variables
  integer(kind=iintegers) :: maxlocat, i
  real(kind=ireals), allocatable :: bvf(:)
  real(kind=ireals) :: zv, zm, x, y, x1, x2
  real(kind=ireals) :: u1, u2, v1, v2, Sal1, Sal2

! Mixed-layer depth - depth of maximal heat flux
!  maxlocat = 1
!  do i = 2, M
!    if (k_turb_T_flux(i) + T_massflux(i) > &
!    & k_turb_T_flux(maxlocat) + T_massflux(maxlocat)) maxlocat = i
!  enddo
!  H_mixed_layer = dzeta_05int(maxlocat)*h1

! Mixed-layer depth - depth of maximal Bruent-Vaisala frequency
!  allocate (bvf(2:M-1))      
!  do i = 2, M-1
!    bvf(i) = g/row0*(row(i+1) - row(i-1))/(h1*(ddz(i-1) + ddz(i)))
!  enddo
!  maxlocat = maxloc(bvf,1) + lbound(bvf,1) - 1
!  deallocate (bvf)
!  H_mixed_layer = dzeta_int(maxlocat)*h1

  call MIXED_LAYER_CALC(row,gsp%ddz,gsp%ddz2,gsp%dzeta_int,gsp%dzeta_05int, &
  & ls%h1,M,i_maxN,H_mixed_layer,maxN)

  w_conv_scale = max( (Buoyancy0 - &
  & (RadWater%integr(1) - RadWater%integr(i_maxN))*g_d_Kelvin0/cw_m_row0) * &
  & H_mixed_layer, 0._ireals)**(1._ireals/3._ireals) 

  T_conv_scale = -eflux0_kinem/(w_conv_scale + small_value)

  Wedderburn = max(min(g*(row(M+1) - row(1))*H_mixed_layer*H_mixed_layer/&
  & ((tau+small_value)*max(bathymwater(1)%Lx,bathymwater(1)%Ly)),&
  & 1.e+2_ireals),1.e-2_ireals)

  ! Depth of volume center (zv) and of center of mass (zm)
  zm = 0.; x1 = 0.
  zv = 0.; x2 = 0.
  do i = 1, M+1
    y = gsp%ddz05(i-1)*ls%h1*bathymwater(i)%area_int
    x = y*gsp%z_full(i)
    zm = zm + row(i)*x
    x1 = x1 + row(i)*y
    zv = zv + x
    x2 = x2 + y
  enddo
  zm = zm/x1
  zv = zv/x2

  LakeNumber = (zm - zv)*ls%vol*row0*g*2.*H_mixed_layer / &
  &            (zv*(tau+small_value)*bathymwater(1)%area_int* &
  &            max(bathymwater(1)%Lx,bathymwater(1)%Ly))

  LakeNumber = max(min(LakeNumber,1.e+2_ireals),1.e-2_ireals)

  Rossby_rad = sqrt(g*max(row(M+1) - row(1),small_value)/row0*H_mixed_layer)/max(kor,small_value)

  ! Thermocline thickness
  x1 = 0.; x2 = ls%h1
  u1 = 0.; u2 = 0.
  v1 = 0.; v2 = 0.
  Sal1 = 0.; Sal2 = 0.
  do i = 1, M
    if (wst%Tw2(i) >= Ttherm0 .and. wst%Tw2(i+1) < Ttherm0) then
      y = gsp%ddz(i)*ls%h1
      x1 = gsp%z_full(i) + (wst%Tw2(i) - Ttherm0)*y/(wst%Tw2(i) - wst%Tw2(i+1))
      x = ( x1 - gsp%z_full(i) )/y
      u1 = wst%u2(i) + x*( wst%u2(i+1) - wst%u2(i) )
      v1 = wst%v2(i) + x*( wst%v2(i+1) - wst%v2(i) )
      Sal1 = wst%Sal2(i) + x*( wst%Sal2(i+1) - wst%Sal2(i) )
    endif
    if (wst%Tw2(i) >= Ttherm1 .and. wst%Tw2(i+1) < Ttherm1) then
      y = gsp%ddz(i)*ls%h1
      x2 = gsp%z_full(i) + (wst%Tw2(i) - Ttherm1)*y/(wst%Tw2(i) - wst%Tw2(i+1))
      x = ( x2 - gsp%z_full(i) )/y
      u2 = wst%u2(i) + x*( wst%u2(i+1) - wst%u2(i) )
      v2 = wst%v2(i) + x*( wst%v2(i+1) - wst%v2(i) )
      Sal2 = wst%Sal2(i) + x*( wst%Sal2(i+1) - wst%Sal2(i) )
    endif
  enddo
  ThermThick = x2 - x1 

  ReTherm = sqrt((u2 - u1)*(u2 - u1) + (v2 - v1)*(v2 - v1))*ThermThick / niu_wat

  RiTherm = g/row0*(DDENS_DTEMP0*(Ttherm1 - Ttherm0) + DDENS_DSAL0*(Sal2 - Sal1))*ThermThick / &
  & ( (u2 - u1)*(u2 - u1) + (v2 - v1)*(v2 - v1) + small_value)

  ! Rp: ratio of density increments by salinity and temperature,
  ! Rpdens : ratio of density fluxes caused by salinity and temperature fluxes
  do i = 1, M
    x1 = SET_DDENS_DTEMP(eos%par,lindens%par,0.5*(wst%Tw2(i+1) + wst%Tw2(i)), &
    & 0.5*(wst%Sal2(i+1) + wst%Sal2(i)))
    x2 = SET_DDENS_DSAL(eos%par,lindens%par, 0.5*(wst%Tw2(i+1) + wst%Tw2(i)), &
    & 0.5*(wst%Sal2(i+1) + wst%Sal2(i)))
    trb%Rp(i) = - x2/x1*(wst%Sal2(i+1) - wst%Sal2(i))/(wst%Tw2(i+1) - wst%Tw2(i) + small_value)
    trb%Rpdens(i) = - wst%lamsal(i)/wst%lamw(i)*x2/x1* &
    & (wst%Sal2(i+1) - wst%Sal2(i))/(wst%Tw2(i+1) - wst%Tw2(i) + small_value)
  enddo

  !if (Wedderburn < 0) then
  !print*, Wedderburn, (row(M+1) - row(1)), tau
  !read*
  !endif
 
  END SUBROUTINE TURB_SCALES


  SUBROUTINE MIXED_LAYER_CALC(row,ddz,ddz2,dzeta_int,dzeta_05int,h,M, &
  & i_mixed_layer,H_mixed_layer,maxN)

  ! Subroutine calculates the mixed-layer depth

  use PHYS_CONSTANTS, only: &
  & g, &
  & row0

  implicit none

! Input variables
  real(kind=ireals), intent(in) :: row (1:M+1)
  real(kind=ireals), intent(in) :: ddz(1:M)
  real(kind=ireals), intent(in) :: ddz2(1:M-1)
  real(kind=ireals), intent(in) :: dzeta_int(1:M+1)
  real(kind=ireals), intent(in) :: dzeta_05int(1:M)
 
  real(kind=ireals), intent(in) :: h

  integer(kind=iintegers), intent(in) :: M

! Output variables
  integer(kind=iintegers), intent(out) :: i_mixed_layer
  real(kind=ireals), intent(out) :: H_mixed_layer
  real(kind=ireals), intent(out) :: maxN

! Local variables
  integer(kind=iintegers) :: i, ml
  real(kind=ireals) :: a, b, c, x1, x2

  real(kind=ireals), allocatable :: bvf(:)

! Mixed-layer depth - depth of maximal heat flux
!  maxlocat = 1
!  do i = 2, M
!    if (k_turb_T_flux(i) + T_massflux(i) > &
!    & k_turb_T_flux(maxlocat) + T_massflux(maxlocat)) maxlocat = i
!  enddo
!  H_mixed_layer = dzeta_05int(maxlocat)*h1

  allocate (bvf(1:M))      
  do i = 1, M
    bvf(i) = g/row0*(row(i+1) - row(i))/(h*ddz(i))
  enddo

  ! Mixed-layer depth - depth of maximal Bruent-Vaisala frequency
  ml = maxloc(bvf,1) !+ lbound(bvf,1) - 1
  maxN = sqrt(max(bvf(ml),0._ireals))

  if ((row(M) - row(1))/row0 < 1.E-3) ml = M

  !! Mixed-layer depth as a depth of sharp increase of static stability
  !ml = M
  !do i = ml, 1, -1
  !  if (abs(bvf(i)/bvf(ml)) < 7.E-1_ireals) then
  !    ml = i
  !    exit
  !  endif
  !enddo

  ! Maximum of quadratic function
  !if (ml >= 2 .and. ml <= M-1) then
  !  x1 = (bvf(ml+1) - bvf(ml)  )/(h*(dzeta_05int(ml+1) - dzeta_05int(ml  )))
  !  x2 = (bvf(ml)   - bvf(ml-1))/(h*(dzeta_05int(ml  ) - dzeta_05int(ml-1)))
  !  a =   ( x1 - x2 ) / (h*(dzeta_05int(ml+1) - dzeta_05int(ml-1)))
  !  b =   x1 - a*(h*(dzeta_05int(ml+1) + dzeta_05int(ml  )))
  !  H_mixed_layer = - 0.5*b/a
  !else
    H_mixed_layer = dzeta_05int(ml)*h
  !endif
  deallocate (bvf)

  i_mixed_layer = ml

  !if (H_mixed_layer < 0.05) then
  !  H_mixed_layer = dzeta_05int(M)*h
  !  i_mixed_layer = M
  !endif

  ! H_mixed_layer = dzeta_int(ml)*h

  END SUBROUTINE MIXED_LAYER_CALC


  !> The function W_SEDIM calculates the sedimentation speed
  !! of organic particles, following Song et al., 2008, 
  !! DOI: 10.3882/j. issn. 1674-2370.2008.01.005
  FUNCTION W_SEDIM(d,ind)
     
  use PHYS_CONSTANTS, only : g, niu_wat

  implicit none

  ! Input variables
  !> Particle diameter, m
  real(kind=ireals), intent(in) :: d
  !> Switch between \f$Re\rightarrow 0\f$ and \f$Re \geq 10^5\f$ limits
  integer(kind=iintegers), intent(in) :: ind

  ! Output variables
  !> Sedimentation speed, m/s
  real(kind=ireals) :: W_SEDIM

  ! Local variables
  real(kind=ireals), parameter :: A = 30., B = 1.25 
  real(kind=ireals), parameter :: delta = 0.25 ! Density parameter for organic particles, Avnimelech et al. 2001

  real(kind=ireals), save :: A1, B1
  logical, save :: firstcall = .true.

  if (firstcall) then
    A1 = 4.*delta*g/(3.*A*niu_wat) 
    B1 = 4.*delta*g/(3.*B)
  endif

  if     (ind == 1) then
    W_SEDIM = A1*d*d
  elseif (ind == 2) then
    W_SEDIM = sqrt(B1*d)
  endif

  if (firstcall) firstcall = .false.
  END FUNCTION W_SEDIM


  FUNCTION SHORTRAD(shortrad0,z)

  implicit none

  real(kind=ireals), intent(in) :: shortrad0 ! Radiation at the surface
  real(kind=ireals), intent(in) :: z ! depth, m

! Shortwave radiation in water column from Arctic ocean model by N.Yakovlev
! Penetrating into ocean radiation parameters (Paulson&Simpson, 1977).
  real(kind=ireals), parameter :: Ra = 0.68 ! Fraction of light frequency radiation
  real(kind=ireals), parameter :: dzi1 = 1.2 ! Penetration of high frequency, m
  real(kind=ireals), parameter :: dzi2 = 28. ! Penetration of low frequency, m

  real(kind=ireals) :: SHORTRAD

  SHORTRAD = shortrad0*(Ra*exp(-z/dzi1)+(1.-Ra)*exp(-z/dzi2))

  END FUNCTION SHORTRAD


  FUNCTION WATER_ALBEDO(sinh0)
  implicit none

  real(kind=ireals) :: WATER_ALBEDO

! Input variables
  real(kind=ireals), intent(in) :: sinh0

! Local variables
  real(kind=ireals), parameter :: const1 = 0.05d0
  real(kind=ireals), parameter :: const2 = 0.15d0

  WATER_ALBEDO = const1/(sinh0 + const2)

  END FUNCTION WATER_ALBEDO


  FUNCTION F_AI(T,h)
! Ice albedo, from Arctic ocean model by N.Yakovlev
  implicit none 

! Input variables
  real(kind=ireals), intent(in) :: T, h

  real(kind=ireals) :: F_AI
  real(kind=ireals), parameter :: aow = 0.10, & ! AOMIP
  &                     ai = 0.65 ! Boulder Ice CCSM3 (Ice version 4, 2002) 

  if (T .LT. -1.0)then ! -1C - see BoulderIce 
    if (h .LE. 50.)then ! see BoulderIce - this is a good approx.
      F_ai=aow +(ai-aow)*h/50.
    else
      F_ai= ai
    end if
  else
    if(h.LE.50.)then ! see BoulderIce - this is a good approx.
      F_ai=aow +(ai-0.075*(T+1.)-aow)*h/50. 
    else
      F_ai= ai -0.075*(T+1.)
    end if
  end if

  return
  END FUNCTION F_AI
 
  
  FUNCTION F_AS(T)
! Snow albedo, from Arctic ocean model by N.Yakovlev
  if(T .LT. -1.0)then ! -1C - see BoulderIce

    F_as=0.80 ! Old Aomip, close to CCSM3

  else

    F_as=0.80 -0.1*(T+1.0) ! Pinto, 1999, SHEBA + Boulder Ice, CCSM3

  end if

  return
  END FUNCTION F_AS


  FUNCTION NETLWRAD(T,Ta,epsa,cloud,emis)
! Net longwave radiation at the surface, Rosati&Miyakoda, 1988
  use PHYS_CONSTANTS, only: &
  & sigma

  implicit none

  real(kind=ireals), intent(in) :: T ! Surface temperature, K
  real(kind=ireals), intent(in) :: Ta ! Air temperature, K
  real(kind=ireals), intent(in) :: epsa ! Water vapor pressure, Pa
  real(kind=ireals), intent(in) :: cloud ! Cloudiness, fraction
  real(kind=ireals), intent(in) :: emis ! Emissivity, fraction

  real(kind=ireals) :: NETLWRAD

  real(kind=ireals) :: AA, BB

  AA = emis*sigma*(0.39-0.005*sqrt(epsa))*(1.-0.8*cloud)
  BB = 4.*emis*sigma
  NETLWRAD = - (AA*T + BB*(T-Ta))*T**3 ! TC - corrected

  END FUNCTION NETLWRAD


  FUNCTION EXTINCT_SNOW(snow_density)
  implicit none

  real(kind=ireals) :: EXTINCT_SNOW

! Input variables
  real(kind=ireals), intent(in) :: snow_density ! Snow density, kg/m**3

! Local variables
  real(kind=ireals), parameter :: const1 = 0.13d0
  real(kind=ireals), parameter :: const2 = 3.4d0
  real(kind=ireals), parameter :: extinct_snow_max = 1.d+7

  EXTINCT_SNOW = exp(-(const1*snow_density+const2) )
  
! For Sparkling2002-2005 experiment
! EXTINCT_SNOW = exp(-extinct_snow_max)


  END FUNCTION EXTINCT_SNOW


  real(kind=ireals) FUNCTION UNFRWAT(T,i)
  use DRIVING_PARAMS
  use ARRAYS_SOIL, only : WLM0, WLM7
  
!CALCULATION OF LIQUID WATER CONTENT IN FREEZING SOIL
!T - DEG. C; WLM0,WLM7,WLMAX - KG/KG

  implicit none
  real(kind=ireals) T
  integer(kind=iintegers) i

  unfrwat = (WLM0(i)-WLM7(i))*exp(T/3.) + WLM7(i)
  !unfrwat = 0 
  
  RETURN
  END FUNCTION UNFRWAT


  real(kind=ireals) FUNCTION WL_MAX(por,rosdry,wi,mh)
  use PHYS_CONSTANTS, only: &
  & row0, &
  & roi
  use METH_OXYG_CONSTANTS, only : &
  & densmh
  implicit none
  real(kind=ireals), intent(in):: por,rosdry,wi,mh
  real(kind=ireals) :: work

  if (mh > 0.) then
    work = densmh
  else
    work = roi
  endif

  WL_MAX = por*row0/(rosdry*(1 - por)) - wi*row0/work
  END FUNCTION WL_MAX


  real(kind=ireals) FUNCTION WI_MAX(por,rosdry,mh)
  use PHYS_CONSTANTS, only: &
  & roi
  use METH_OXYG_CONSTANTS, only : &
  & densmh
  implicit none
  real(kind=ireals), intent(in):: por,rosdry,mh
  real(kind=ireals) :: work

  if (mh > 0.) then
    work = densmh
  else
    work = roi
  endif

  WI_MAX = por*work/(rosdry*(1 - por))
  END FUNCTION WI_MAX

  
  FUNCTION SOIL_COND_JOHANSEN(wl,wi,rosdry,por)
  
  use PHYS_CONSTANTS, only: &
  & row0, roi, &
  & row0_d_roi, roi_d_row0, &
  & lamw0, lami
  
  implicit none
  
  real(kind=ireals) :: SOIL_COND_JOHANSEN
  
! Input variables
  real(kind=ireals), intent(in) :: wl, wi ! liquid water and ice content 
                                ! (in respect to dry soil mass), kg/kg
  real(kind=ireals), intent(in) :: rosdry ! dry soil (soil particles) density, kg/m**3
  real(kind=ireals), intent(in) :: por ! soil porosity, m**3/m**3
  
! Local variables
  real(kind=ireals) :: water_vol_ratio, ice_vol_ratio ! water and ice volume ratio 
                                            ! (in respect to bulk soil volume), m**3/m**3
  real(kind=ireals) :: waterice_vol_ratio ! total volume ratio of water and ice, m**3/m**3
  real(kind=ireals) :: Kersten ! Kersten number, n/d
  real(kind=ireals) :: CK_const, CK_const1, CK_const2
  real(kind=ireals) :: lambda_sat, lambda_dry ! heat conduction coefficients for saturated
                                    ! and dry soil, W/(m*K)
                                    
  real(kind=ireals) :: water_sat_ratio, ice_sat_ratio                       
                   
  real(kind=ireals), save :: CK_consts(1:4) ! Cote and Konrad (2005) constants, n/d
  real(kind=ireals), save :: CK_consts1(1:3) ! Cote and Konrad (2005) constants, W/(m*K)
  real(kind=ireals), save :: CK_consts2(1:3) ! Cote and Konrad (2005) constants, n/d
  real(kind=ireals), save :: quartz_ratio ! quartz content of the total solids content
  real(kind=ireals), save :: lambda_quartz ! heat conduction coefficient for quartz, W/(m*K)
  real(kind=ireals), save :: lambda_othmin ! heat conduction coefficient for non-quartz
                                 ! minerals, W/(m*K)
  real(kind=ireals), save :: lambda_solids ! heat conduction coefficient for soild part
                                 ! of the soil, W/(m*K)
  
  logical, save :: firstcall = .true.
  
  if (firstcall) then
    CK_consts(1) = 4.6d0 ! for gravel and coarse sand
    CK_consts(2) = 3.55d0 ! for medium and fine sand
    CK_consts(3) = 1.9d0 ! silty and clay soils
    CK_consts(4) = 0.6d0 ! organic fibrous soils
    
    CK_consts1(1) = 1.7d0 ! for crashed rocks
    CK_consts1(2) = 0.75d0 ! for mineral soils
    CK_consts1(3) = 0.3d0 ! organic fibrous soils
    
    CK_consts2(1) = 1.8d0 ! for crashed rocks
    CK_consts2(2) = 1.2d0 ! for mineral soils
    CK_consts2(3) = 0.87d0 ! organic fibrous soils
    
    quartz_ratio = 0.1d0
    lambda_quartz = 7.7d0
    if (quartz_ratio > 0.2d0) then
      lambda_othmin = 2.d0
    else
      lambda_othmin = 3.d0
    endif
    lambda_solids = lambda_quartz**quartz_ratio * &
    &               lambda_othmin**(1.d0-quartz_ratio)
  endif
  
! Convertation from mass ratios to volume ratios      
  water_vol_ratio = wl / &
  & (por*(wl + row0/roi*wi + row0/rosdry))
  ice_vol_ratio = wi / &
  & (por*(wi + roi/row0*wl + roi/rosdry))
  
  waterice_vol_ratio = water_vol_ratio + ice_vol_ratio
  
  CK_const = CK_consts(3) ! silty and clay soils are assumed
  Kersten = CK_const*waterice_vol_ratio / &
  & (1.d0 + (CK_const - 1.d0)*waterice_vol_ratio)
  
  water_sat_ratio = por*water_vol_ratio/waterice_vol_ratio
  ice_sat_ratio = por*ice_vol_ratio/waterice_vol_ratio
  
! This is the original formula from Johansen (1975) extended for the case
! with the ice content
  lambda_sat = lambda_solids**(1.d0 - por) * &
  & lamw0**(water_sat_ratio)*lami**(ice_sat_ratio)
  
  CK_const1 = CK_consts1(2) ! mineral soils are assumed
  CK_const2 = CK_consts2(2) ! mineral soils are assumed
! Cote and Konrad (2005) formula for heat conduction coefficient of dry soil
  lambda_dry = CK_const1*10.d0**(-CK_const2*por)
  
! The normalized soil conductivity concept by Johansen (1975)      
  SOIL_COND_JOHANSEN = (lambda_sat - lambda_dry)*Kersten + lambda_dry
  
  if (firstcall) firstcall = .false.
  END FUNCTION SOIL_COND_JOHANSEN


  FUNCTION PRESMHYDRDISS(T,S)

  implicit none

  real(kind=ireals), intent(in) :: T ! Temperature, K
  real(kind=ireals), intent(in) :: S ! Salinity, ppt

  real(kind=ireals) :: PRESMHYDRDISS ! Pressure of lower limit of methane hydrate stability, Pa 
                           ! (Tishchenko et al. 2005)

  real(kind=ireals), parameter :: c1 = 1.6444866d+3, c2 = 0.1374178, c3 = 5.4979866d+4
  real(kind=ireals), parameter :: c4 = 2.64118188d+2, c5 = 1.1178266d+4, c6 = 7.67420344
  real(kind=ireals), parameter :: c7 = 4.515213d-3, c8 = 2.04872879d+5, c9 = 2.17246046d+3
  real(kind=ireals), parameter :: c10 = 1.70484431d+2, c11 = 0.118594073, c12 = 7.0581304d-5
  real(kind=ireals), parameter :: c13 = 3.09796169d+3, c14 = 33.2031996 

  real(kind=ireals) :: logpres

  logpres =  - c1 - c2*T + &
  & c3/T + c4*log(T) + S* &
  & (c5 + c6*T - c7*T**2 - c8/T - c9*log(T)) + S**2* &
  & (c10 + c11*T - c12*T**2 - c13/T - c14*log(T))
  PRESMHYDRDISS = exp(logpres)*1.d+6 !Converting to Pa

  END FUNCTION PRESMHYDRDISS



  FUNCTION TEMPMHYDRDISS(p)

! The temperature (K) of methane hydrate dissociation in seawater 
! at 35 ppt salinity according to (Dickens and Quinby-Hunt, 1994, GRL)
  use PHYS_CONSTANTS, only : &
  & Kelvin0

  implicit none

  real(kind=ireals), intent(in) :: p ! Pressure, Pa

  real(kind=ireals) :: TEMPMHYDRDISS

  real(kind=ireals), parameter ::  c1 = 3.79E-3, c2 = 2.83E-4

  TEMPMHYDRDISS = 1./(c1 - c2*log10(p*1.E-6)) - Kelvin0

  END FUNCTION TEMPMHYDRDISS


  FUNCTION MELTINGPOINT(S,p,methhydr)

  use DRIVING_PARAMS, only : nmeltpoint

  implicit none

  real(kind=ireals), intent(in) :: S, p 
  real(kind=ireals), intent(in) :: methhydr
  real(kind=ireals) :: MELTINGPOINT

  if (methhydr > 0.) then
    MELTINGPOINT = TEMPMHYDRDISS(p)
  else
    MELTINGPOINT = MELTPNT(S,p,nmeltpoint%par)
  endif

  END FUNCTION MELTINGPOINT


  !> Function calculates ice salinity at the phase change front
  !! according to V.L.Tsurikov (Nazintsev and Panov, 2000)
  FUNCTION SALICEBOT(dhbdt,salwat,por)
  use PHYS_CONSTANTS, only : row0, roi
  implicit none
  !>Input/output variables
  real(kind=ireals), intent(in) :: dhbdt !> the rate of water freeze at the phase change front, m/s
  real(kind=ireals), intent(in) :: salwat !> water salinity below the phase change front, kg/kg
  real(kind=ireals), intent(in) :: por  !> ice porosity, m**3/m**3
  real(kind=ireals) :: SALICEBOT

  !>Local variables
  real(kind=ireals), parameter :: alpha = 7., beta = 7., gamma_ = 10.3
  real(kind=ireals) :: w
  
  w = max(dhbdt,0._ireals)*3.6E+6 !Converting from m/s to mm/h

  SALICEBOT = (row0*por + (1.-por)*roi) * &
  &           salwat*alpha*sqrt(w)/(beta*sqrt(w) + gamma_) 

  END FUNCTION SALICEBOT


  FUNCTION REACPOT_ARRHEN &
  & (delta_E, temp, temp0, reacpot_arrhen0)

! The FUNCTION REACPOT_ARRHEN calculates the reaction potential
! according to Arrhenius equation

  use PHYS_CONSTANTS, only : &
  & R_univ

  implicit none      
  
! Input variables
  real(kind=ireals), intent(in) :: delta_E ! activation energy, J/mol
  real(kind=ireals), intent(in) :: temp    ! temperature, Kelvin
  real(kind=ireals), intent(in) :: temp0   ! reference temperature, Kelvin
  real(kind=ireals), intent(in) :: reacpot_arrhen0 ! reaction potential at the
                                         ! reference temeperature, mol/(m**3*s)
                                         
! Output variables
  real(kind=ireals) :: REACPOT_ARRHEN
  
! Local variables and constants
  
  REACPOT_ARRHEN = &
  & reacpot_arrhen0*exp(delta_E/R_univ*(1./temp0 - 1./temp))
  
  END FUNCTION REACPOT_ARRHEN
  
  
  FUNCTION HENRY_CONST(henry_const0, temp_dep, temp_ref, temp) !, radius)
  
  ! Function HENRY_CONST calculates the Henry constant of a substance for a given temperature
  
  use PHYS_CONSTANTS, only : &
  & surf_tension_wat, R_univ

  implicit none
  
  real(kind=ireals) :: HENRY_CONST
  
  ! Input variables
  real(kind=ireals), intent(in) :: henry_const0 ! Henry constant at the reference temperature, mol/(m**3*Pa)
  real(kind=ireals), intent(in) :: temp_dep ! Temperature dependence (enthalpy solution devided by universal gas constant), K
  real(kind=ireals), intent(in) :: temp_ref ! Reference temperature, K
  real(kind=ireals), intent(in) :: temp ! Temperature, K
!  real(kind=ireals), intent(in) :: radius ! Curvature radius, m, positive for drops, negative for bubbles and zero for flat surface
  
  HENRY_CONST = henry_const0*exp(temp_dep*(1./temp-1./temp_ref))
! Effect of bubble surface curvatue on saturation pressure
!  if (radius /= 0.d0) HENRY_CONST = HENRY_CONST* &
!  & exp(-2.*surf_tension_wat*mol_vol_water/(radius*R_univ*temp))
  
  END FUNCTION HENRY_CONST


  !> Function HC_CORR_CARBEQUIL returns a correction multiplier for 
  !! CO_2 Henry constant to get solubulity, taking into account
  !! carbonate equilibrium.
  !! \f[
  !! HC\_CORR\_CARBEQUIL = (1+k_1 10^{pH}+k_1*k_2 10^{2pH})
  !! \f]
  !! Concomitantly, it is a ratio of DIC to CO_2 molar concenatrations.
  !! DIC molar concentration is a molar concentration of C atoms in DIC.
  FUNCTION HC_CORR_CARBEQUIL(temp,pH)

  implicit none

  !Input variables
  !>Water temperature, deg. Kelvin
  real(kind=ireals), intent(in) :: temp 
  !>pH
  real(kind=ireals), intent(in) :: pH 

  real(kind=ireals) :: HC_CORR_CARBEQUIL

  !Local variables
  real(kind=ireals), parameter :: k1_0 = 4.3E-7, k2_0 = 4.7E-11
  real(kind=ireals), parameter :: Eact_d_RT_1 = 921.4, Eact_d_RT_2 = 1787.4
  real(kind=ireals), parameter :: temp_ref = 298.

  real(kind=ireals) :: k1,k2

  k1 = k1_0*exp( - Eact_d_RT_1 * (1./temp - 1./temp_ref) )
  k2 = k2_0*exp( - Eact_d_RT_2 * (1./temp - 1./temp_ref) )
  HC_CORR_CARBEQUIL = (1. + k1*10._ireals**(pH) + k1*k2*10._ireals**(2*pH))

  END FUNCTION HC_CORR_CARBEQUIL


  FUNCTION DIFF_WATER_OXYGEN(temp_C)
  
  ! Function DIFF_WATER_OXYGEN calculates the molecular diffusivity 
  ! of oxygen dissolved in liquid water, m**2/s
  ! (Broecker and Peng, 1974)
  
  implicit none
  
  real(kind=ireals) :: DIFF_WATER_OXYGEN
  
  ! Input variables
  real(kind=ireals), intent(in) :: temp_C ! temperature, degrees Celsius

  ! Local variables
  real(kind=ireals), parameter :: const1 = 1.11d-9
  real(kind=ireals), parameter :: const2 = 4.96d-11
  
  DIFF_WATER_OXYGEN = const1 + const2*temp_C
   
  END FUNCTION DIFF_WATER_OXYGEN

  
  FUNCTION DIFF_WATER_METHANE(temp_C)
  
  ! Function DIFF_WATER_METHANE calculates the molecular diffusivity 
  ! of methane dissolved in liquid water, m**2/s
  ! (Broecker and Peng, 1974)
  
  implicit none
  
  real(kind=ireals) :: DIFF_WATER_METHANE
  
  ! Input variables
  real(kind=ireals), intent(in) :: temp_C ! temperature, degrees Celsius

  ! Local variables
  real(kind=ireals), parameter :: const1 = 9.798d-10
  real(kind=ireals), parameter :: const2 = 2.986d-11
  real(kind=ireals), parameter :: const3 = 4.381d-13
  
  DIFF_WATER_METHANE = const1 + const2*temp_C + const3*temp_C*temp_C
   
  END FUNCTION DIFF_WATER_METHANE
  
  
  FUNCTION DIFF_AIR_METHANE(temp_C)
  
  ! Function DIFF_WATER_METHANE calculates the molecular diffusivity 
  ! of methane dissolved in liquid water, m**2/s
  ! (Lerman, 1979)
  
  implicit none
  
  real(kind=ireals) :: DIFF_AIR_METHANE
  
  ! Input variables
  real(kind=ireals), intent(in) :: temp_C ! temperature, degrees Celsius

  ! Local variables
  real(kind=ireals), parameter :: const1 = 1.875d-5
  real(kind=ireals), parameter :: const2 = 1.3d-7
  
  DIFF_AIR_METHANE = const1 + const2*temp_C
   
  END FUNCTION DIFF_AIR_METHANE

  
  FUNCTION DIFF_WATER_CARBDI(temp_C)
  
  ! Function DIFF_WATER_CARBDI calculates the molecular diffusivity 
  ! of crabon dioxide dissolved in liquid water, m**2/s
  ! (Broecker and Peng, 1974)
  
  implicit none
  
  real(kind=ireals) :: DIFF_WATER_CARBDI
  
  ! Input variables
  real(kind=ireals), intent(in) :: temp_C ! temperature, degrees Celsius

  ! Local variables
  real(kind=ireals), parameter :: const1 = 9.39d-10
  real(kind=ireals), parameter :: const2 = 2.671d-11
  real(kind=ireals), parameter :: const3 = 4.095d-13
  
  DIFF_WATER_CARBDI = const1 + const2*temp_C + const3*temp_C*temp_C
   
  END FUNCTION DIFF_WATER_CARBDI

  
  FUNCTION DIFF_AIR_CARBDI(temp_C)
  
  ! Function DIFF_WATER_CARBDI calculates the molecular diffusivity 
  ! of carbon dioxide dissolved in liquid water, m**2/s
  ! (Lerman, 1979)
  
  implicit none
  
  real(kind=ireals) :: DIFF_AIR_CARBDI
  
  ! Input variables
  real(kind=ireals), intent(in) :: temp_C ! temperature, degrees Celsius

  ! Local variables
  real(kind=ireals), parameter :: const1 = 1.35d-5
  real(kind=ireals), parameter :: const2 = 0.9d-7
        
  DIFF_AIR_CARBDI = const1 + const2*temp_C
   
  END FUNCTION DIFF_AIR_CARBDI
  
  
  FUNCTION SCHMIDT_NUMBER_METHANE(tempC)
  
  implicit none
  
  real(kind=ireals) :: SCHMIDT_NUMBER_METHANE
  
  ! Input variables      
  real(kind=ireals), intent(in) :: tempC
  
  ! Local variables
  real(kind=ireals), parameter :: const1 = 1.898d+3
  real(kind=ireals), parameter :: const2 = -1.101d+2
  real(kind=ireals), parameter :: const3 = 2.834
  real(kind=ireals), parameter :: const4 = -2.791d-2
  
  SCHMIDT_NUMBER_METHANE = &
  & const1 + const2*tempC + const3*tempC**2 + const4*tempC**3
  
  END FUNCTION SCHMIDT_NUMBER_METHANE
  
  
  FUNCTION SCHMIDT_NUMBER_OXYGEN(tempC)
  
  implicit none
  
  real(kind=ireals) :: SCHMIDT_NUMBER_OXYGEN
  
  ! Input variables      
  real(kind=ireals), intent(in) :: tempC
  
  ! Local variables
  real(kind=ireals), parameter :: const1 = 1.8006d+3
  real(kind=ireals), parameter :: const2 = -1.201d+2
  real(kind=ireals), parameter :: const3 = 3.7818
  real(kind=ireals), parameter :: const4 = -4.7608d-2
  
  SCHMIDT_NUMBER_OXYGEN = &
  & const1 + const2*tempC + const3*tempC**2 + const4*tempC**3
  
  END FUNCTION SCHMIDT_NUMBER_OXYGEN
  
  
  FUNCTION SCHMIDT_NUMBER_CARBDI(tempC)
  
  implicit none
  
  real(kind=ireals) :: SCHMIDT_NUMBER_CARBDI
  
  ! Input variables      
  real(kind=ireals), intent(in) :: tempC
  
  ! Local variables
  real(kind=ireals), parameter :: const1 = 1.911d+3
  real(kind=ireals), parameter :: const2 = -1.137d+2
  real(kind=ireals), parameter :: const3 = 2.967
  real(kind=ireals), parameter :: const4 = -2.943d-2
  
  SCHMIDT_NUMBER_CARBDI = &
  & const1 + const2*tempC + const3*tempC**2 + const4*tempC**3
  
  END FUNCTION SCHMIDT_NUMBER_CARBDI
  
  
  FUNCTION GAS_WATATM_FLUX &
  & (tempC,wind10,surf_conc,partial_pressure,henry_const,gasindic,eps)

  use PHYS_CONSTANTS, only: &
  & niu_wat
        
  ! Subroutine GAS_WATATM_FLUX calculates the upward gas flux at the water-air interface, mol/(m**2*s)
  ! following formulations by (McGillis et al., 2000; Cole and Caraco, 1998; Riera et al., 1999) 
  ! and others (described in Wania, 2007)
        
  implicit none
  
  real(kind=ireals) :: GAS_WATATM_FLUX
  
  ! Input variables
  real(kind=ireals), intent(in) :: tempC ! water surface temperature, Celsius
  real(kind=ireals), intent(in) :: wind10 ! wind speed at 10 m above the surface, m/s
  real(kind=ireals), intent(in) :: surf_conc ! gas concentration at the water surface, mol/m**3
  real(kind=ireals), intent(in) :: partial_pressure ! partial pressure of a gas in the atmosphere, Pa
  real(kind=ireals), intent(in) :: henry_const ! Henry constant of a gas, mol/(m**3*Pa)
  real(kind=ireals), intent(in) :: eps  ! kinetic energy dissipation
  
  integer(kind=iintegers), intent(in) :: gasindic ! gas indicator: methane - 8, oxygen - 9, carbon dioxide - 10
        
  ! Local constants
  real(kind=ireals), parameter :: constvel1 = 5.75d-6 !(Cole and Coraco et al.,1998) ! 3.78d-8 ! SI units
  real(kind=ireals), parameter :: constvel2 = 5.97d-7 !(Cole and Coraco et al.,1998) ! 1.36d-9 ! SI units
  real(kind=ireals), parameter :: c1 = 0.5 ! (JOUNI J. HEISKANEN et al., 2014)
  ! real(kind=ireals), parameter :: c1 = 1.2 (MacIntyre et al., 2010)
  real(kind=ireals), parameter :: Schmidt_scale = 600.
  real(kind=ireals), parameter :: Schmidt_number_min = 1.
  real(kind=ireals), parameter :: wind10_power = 1.7
  
  integer(kind=iintegers), parameter :: methane_indic = 8
  integer(kind=iintegers), parameter :: oxygen_indic = 9
  integer(kind=iintegers), parameter :: carbdi_indic = 10
  integer(kind=iintegers), parameter :: COLE_CARACO = 1 ! model depending on wind-speed (Cole and Coraco et al.,1998)
  integer(kind=iintegers), parameter :: SURFACE_RENEWAL = 2 ! surface renewal model (Soloviev and Schluessel, 1994;&
  ! & MacIntyre et al.,1995, 2001; Zappa et al., 2007; Turney and Banerjee,2008)
  
  ! Local variables
  real(kind=ireals) :: Schmidt_number
  real(kind=ireals) :: piston_velocity, piston_velocity600
  real(kind=ireals) :: surf_conc_atmequil

  integer(kind=iintegers) :: piston_velocity_model ! choice of model

  piston_velocity_model = SURFACE_RENEWAL

  if (piston_velocity_model == COLE_CARACO) then
    piston_velocity600 = constvel1 + constvel2*wind10**wind10_power
  elseif (piston_velocity_model == SURFACE_RENEWAL) then
    piston_velocity600 = c1*(eps*niu_wat)**0.25/sqrt(Schmidt_scale)
  endif
  
  ! piston_velocity600 = constvel1 + constvel2*wind10**wind10_power      
  if     (gasindic == methane_indic) then
    Schmidt_number = max(SCHMIDT_NUMBER_METHANE(tempC),Schmidt_number_min)
  elseif (gasindic == oxygen_indic) then
    Schmidt_number = max(SCHMIDT_NUMBER_OXYGEN(tempC),Schmidt_number_min)
  elseif (gasindic == carbdi_indic) then
    Schmidt_number = max(SCHMIDT_NUMBER_CARBDI(tempC),Schmidt_number_min)
  endif
  piston_velocity = piston_velocity600*sqrt(Schmidt_scale/Schmidt_number)
  surf_conc_atmequil = partial_pressure*henry_const
  GAS_WATATM_FLUX = piston_velocity*(surf_conc-surf_conc_atmequil)
  
  END FUNCTION GAS_WATATM_FLUX
  
  
  FUNCTION CHARNOCK_Z0(velfrict)

  use PHYS_CONSTANTS, only : &
  & g, niu_atm, const2 => charnock_const
  
  implicit none

  real(kind=ireals) :: CHARNOCK_Z0

  ! Input variables
  real(kind=ireals), intent(in) :: velfrict

  ! Local variables
  real(kind=ireals), parameter :: const1 = 0.111
  ! real(kind=ireals), parameter :: const2 = 0.0144
  ! real(kind=ireals), parameter :: const2 = 0.35
  real(kind=ireals), parameter :: roughness_min = 1.d-5
  real(kind=ireals), parameter :: roughness_max = 1.1d-1

  CHARNOCK_Z0 = &
  & min(max(const1*niu_atm/velfrict + const2*velfrict**2/g, &
  & roughness_min),roughness_max)

  END FUNCTION CHARNOCK_Z0


  FUNCTION DOMWAVESPEED(ufr,fetch)

! Calculating fetch-dependent dominant wave speed,
! following Wuest and Lorke (2003)

  use PHYS_CONSTANTS, only : &
  & g

  implicit none
  real(kind=ireals) :: DOMWAVESPEED

  ! Input variables
  real(kind=ireals), intent(in) :: ufr ! Friction velocity, m/s
  real(kind=ireals), intent(in) :: fetch ! Wind fetch, m

  ! Locals
  real(kind=ireals), parameter :: C1 = 0.051, C2 = 0.96
  real(kind=ireals), save :: C3
  logical, save :: firstcall = .true.

  if (firstcall) then
    C3 = (C1/C2)**(2./3.)
  endif

  DOMWAVESPEED = C3*(ufr*g*fetch)**(1./3.)

  if (firstcall) firstcall = .false.
  END FUNCTION DOMWAVESPEED


  FUNCTION H13(momentum_flux,fetch)

! Significant wave height, Hasselman et al. (1973)

  use PHYS_CONSTANTS, only : &
  & g, roa0

  implicit none
  
  real(kind=ireals) :: H13
 
! Input variables
  real(kind=ireals), intent(in) :: momentum_flux, fetch

! Locals
  real(kind=ireals), parameter :: const = 0.051

  H13 = const*sqrt(momentum_flux*fetch/(g*roa0))

  END FUNCTION H13



  FUNCTION LOGFLUX(vel, ds, z, z0, z0s, mult, ind)

  ! Computes a flux of substance in neutral (logarythmic) surface layer
  ! LOGFLUX is positive in the direction in respect to which the difference
  ! ds is taken

  use PHYS_CONSTANTS, only : &
  & kappa

  implicit none

  ! Input variables
  real(kind=ireals), intent(in) :: vel, ds, z, z0, z0s, mult
  integer(kind=iintegers), intent(in) :: ind
  
  ! Output varaible
  real(kind=ireals) :: LOGFLUX

  ! Local variables
  real(kind=ireals), save :: kappa2
  logical, save :: firstcall = .true.

  if (firstcall) kappa2 = kappa*kappa

  if (ind == 1) then ! return flux
    LOGFLUX = - mult*kappa2*vel*ds/(log(z/z0)*log(z/z0s))
  elseif (ind == 2) then ! return exchange coefficient
    LOGFLUX =   mult*kappa2*vel/(log(z/z0)*log(z/z0s))
  endif

  if (firstcall) firstcall = .false.
  END FUNCTION LOGFLUX


  FUNCTION TEMPWATR(kturb,rad,Tsoil,exchcoef,Tw,Tflux)

  ! Water temperature near lateral bottom from radial temperature
  ! profile scaling

  implicit none

  ! Input variables
  real(kind=ireals), intent(in) :: kturb,rad,Tsoil,exchcoef,Tw,Tflux

  ! Output variable
  real(kind=ireals) :: TEMPWATR

  ! Local variables
  real(kind=ireals) :: C = 1.e+1

  TEMPWATR = Tw + C*rad/kturb*(Tsoil*exchcoef - Tflux)/ & ! 
  &          (1 + C*rad/kturb*       exchcoef)

  !print*, Tsoil*exchcoef, Tflux

  END FUNCTION TEMPWATR


  !> Subroutine calculates background eddy diffusivity after (Hondzo and Stefan, 1993)
  SUBROUTINE DIFFMIN_HS(area,wst,gsp,ls,M)
  use PHYS_CONSTANTS, only : row0, g, cw_m_row0
  use ARRAYS_GRID, only : gridspacing_type
  use ARRAYS_WATERSTATE, only : waterstate_type
  use ARRAYS_BATHYM, only : layers_type
  use DRIVING_PARAMS, only : backdiff0
  implicit none
  !Input/output variables
  real(kind=ireals)     , intent(in) :: area
  type(waterstate_type) , intent(in) :: wst
  type(gridspacing_type), intent(in) :: gsp
  type(layers_type)     , intent(in) :: ls
  integer(kind=iintegers), intent(in) :: M
  !Local variables
  integer(kind=iintegers) :: i
  real(kind=ireals) :: x, y

  real(kind=ireals), parameter :: C1 = 8.17E-4, C2 = 0.56, C3 = -0.43, xmin = 7.E-5
  real(kind=ireals), save :: z !A multiplier in expression for background diffusivity
  logical, save :: firstcall = .true.

  if (firstcall) then
    if (backdiff0%par >= 0.) then
      z = backdiff0%par 
    else
      z = C1 !default value
    endif
  endif

  do i = 1, M
    x = g/row0*(wst%row(i+1)-wst%row(i))/(ls%h1*gsp%ddz(i))
    if (x <= 0.) then
      wst%lamw_back(i) = 0.
    else
      x = max(x, xmin)
      y = area*1.E-6 ! converting from m**2 to km**2
      wst%lamw_back(i) = z*(y)**C2 * x**C3
      wst%lamw_back(i) = wst%lamw_back(i) * 1.E-4 !converting from cm**2/s to m**2/s
      wst%lamw_back(i) = wst%lamw_back(i) * cw_m_row0
      wst%lamw_back(i) = wst%lamw_back(i) * 5.E-1 !Calibration multiplyer
    endif
  enddo

  if (firstcall) firstcall = .false.
  END SUBROUTINE DIFFMIN_HS


  !> Subroutine calculates background eddy diffusivity from KPP model 
  !! (Zhang et al., HESS, 2019 and references therein)
  SUBROUTINE DIFFMIN_KPP(wst,gsp,ls,M) 
  use PHYS_CONSTANTS, only : row0, g, cw_m_row0
  use ARRAYS_GRID, only : gridspacing_type
  use ARRAYS_WATERSTATE, only : waterstate_type
  use ARRAYS_BATHYM, only : layers_type
  use DRIVING_PARAMS, only : backdiff0
  implicit none
  !Input/output variables
  type(waterstate_type) , intent(inout) :: wst
  type(gridspacing_type), intent(in) :: gsp
  type(layers_type)     , intent(in) :: ls
  integer(kind=iintegers), intent(in) :: M
  !Local variables
  logical, save :: firstcall = .true.
  real(kind=ireals), parameter :: k0_default = 1.E-5, Ri0 = 7.E-1, p = 3.
  real(kind=ireals), save :: k0
  real(kind=ireals) :: x, y, Ri, dz
  integer(kind=iintegers) :: i !Loop index

  if (firstcall) then
    if (backdiff0%par >= 0.) then
      k0 = backdiff0%par 
    else
      k0 = k0_default !default value
    endif
  endif

  do i = 1, M
    dz = ls%h1*gsp%ddz(i)
    !Brunt-Vaisala frequency squared
    x = g/row0*(wst%row(i+1) - wst%row(i)) / dz
    !Shear frequency squared
    y = ( (wst%u1(i+1) - wst%u1(i))**2 + (wst%v1(i+1) - wst%v1(i))**2 ) / dz**2
    !Gradient Richardson number
    Ri = x / y
    if     (Ri <= 0.) then
      wst%lamw_back(i) = k0
    elseif (Ri < Ri0) then
      wst%lamw_back(i) = k0*(1. - (Ri/Ri0)**2)**p
    else
      wst%lamw_back(i) = 0.
    endif
  enddo
  wst%lamw_back(:) = wst%lamw_back(:) * cw_m_row0

  if (firstcall) firstcall = .false.
  END SUBROUTINE DIFFMIN_KPP

  !> Function HORIZCONC returns a value of dissolved concentration in the mixed layer
  !! at distance r from center of circular lake; according to 
  !! (DelSontro et al., Ecosystems, 2018) extended by inclusion of linear sink term
  FUNCTION HORIZCONC(Caver,r,rL,kdiff,kexch,ksink,hML) result(Cr)
  use NUMERICS, only : BESSI !Bessel function
  implicit none
  !Input/output variables
  real(kind=ireals), intent(in) :: Caver !> horizontally averaged concentration in the ML
  real(kind=ireals), intent(in) :: r     !> distance form lake center, m
  real(kind=ireals), intent(in) :: rL    !> lake radius, m
  real(kind=ireals), intent(in) :: kdiff !> effective diffusion coefficient in horizontal, m**2/s
  real(kind=ireals), intent(in) :: kexch !> gas exchange koefficient with the atmosphere (piston velocity), m/s
  real(kind=ireals), intent(in) :: ksink !> the rate of concentration sink in the ML, s**(-1)
  real(kind=ireals), intent(in) :: hML   !> mixed-layer (ML) depth
  real(kind=ireals)             :: Cr    !> concentration at the distance r from lake center
  !Local variables
  real(kind=ireals) :: x
  x = sqrt((kexch/hML+ksink)/kdiff)
  Cr = 0.5*Caver*rL*BESSI(0_4,x*r)/BESSI(1_4,x*rL)
  END FUNCTION HORIZCONC

  END MODULE PHYS_FUNC
