phys_func.f90

 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
!  common /cloud/ cloud
! data cloud /0./

!  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)   
  
  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

  integer(kind=iintegers), external:: JULIAN_DAY
 
  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


  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