Skip to content
Snippets Groups Projects
phys_func.f90 66.7 KiB
Newer Older
  • Learn to ignore specific revisions
  • Victor Stepanenko's avatar
    Victor Stepanenko committed
     MODULE WATER_DENSITY
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     use LAKE_DATATYPES, only : ireals, iintegers
     use PHYS_CONSTANTS, only : &
     & row0
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     implicit none
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    !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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    !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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    !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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
     real(kind=ireals), save ::  DDENS_DTEMP0, DDENS_DSAL0, DENS_REF
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
     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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
     !Temperature and salinity reference values
     real(kind=ireals), parameter :: temp_ref = 15. ! deg. C
     real(kind=ireals), parameter :: sal_ref = 0. !kg/kg
    
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     contains
    
     SUBROUTINE SET_DENS_CONSTS(eos)
    
     implicit none 
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     
    
     select case (eos)
     case(eos_Hostetler)
       call DDENS_HOSTETLER(temp_ref,DDENS_DTEMP0,DDENS_DSAL0)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     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)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     case(eos_UNESCO)
       call DDENS_UNESCO(temp_ref,DDENS_DTEMP0,DDENS_DSAL0)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     case(eos_McCatcheon)
    
       DDENS_DTEMP0 = DDENS_DTEMP(temp_ref,sal_ref) 
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
       DDENS_DSAL0  = DDENS_DSAL (temp_ref,sal_ref)
    
       DENS_REF = WATER_DENS_TS(temp_ref,sal_ref)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     end select
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     END SUBROUTINE SET_DENS_CONSTS
    
    
    
     FUNCTION SET_DDENS_DTEMP(eos,lindens,Temp,Sal)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
     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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
     ! Local variables
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
     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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
     END FUNCTION SET_DDENS_DTEMP
    
    
    
     FUNCTION SET_DDENS_DSAL(eos,lindens,Temp,Sal)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
     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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
     ! Local variables
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
     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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
     END FUNCTION SET_DDENS_DSAL
    
    
    
     !>Subroutine DENSITY_W calculates density profile
     SUBROUTINE DENSITY_W(M,eos,lindens,Tw,Sal,preswat,row)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
     integer(kind=iintegers), intent(in) :: eos, M, lindens
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
     real(kind=ireals), intent(in) :: Tw(1:M+1), Sal(1:M+1), preswat(1:M+1)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
     !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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
       do i = 1, M+1
    
         row(i) = DENS_REF + DDENS_DTEMP0*(Tw(i) - temp_ref) + &
         &                   DDENS_DSAL0 *(Sal(i) - sal_ref) 
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
       enddo
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     END SUBROUTINE DENSITY_W
    
     real(kind=ireals) FUNCTION WATER_DENS_T_UNESCO(Temp)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    !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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     
     WATER_DENS_T_UNESCO = &
    &  row0*(1+a0+a1*Temp-a2*Temp**2+a3*Temp**3)
     END FUNCTION WATER_DENS_T_UNESCO
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     
     ddens_dtemp = &
     & row0*(a1 - 2.*a2*Temp + 3.*a3*Temp**2)
     ddens_dsal = 0.
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     END SUBROUTINE DDENS_UNESCO
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     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:
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    ! Temp --- water temperature, deg C
    ! Sal  --- water salinity,    kg/kg
      
      implicit none
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      real(kind=ireals), intent(in):: Temp
      real(kind=ireals), intent(in):: Sal
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    ! Converting salinity units from kg/kg to g/kg
      Sal_g_kg = Sal*1.d+3
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      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)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      END FUNCTION WATER_DENS_TS_KIVU
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    ! Converting salinity units from kg/kg to g/kg
      Sal_g_kg = Sal*1.d+3
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      END SUBROUTINE DDENS_KIVU
    
      real(kind=ireals) FUNCTION WATER_DENS_TS(Temp,Sal)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    ! 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)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    ! Input variables:
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    ! Temp --- water temperature, deg C
    ! Sal  --- water salinity,    kg/kg
      
      implicit none
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      real(kind=ireals), intent(in):: Temp
      real(kind=ireals), intent(in):: Sal
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    ! Converting salinity units from kg/kg to g/kg
      Sal_g_kg = Sal*1.d+3
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    ! The first term, dependent on temperature      
      WATER_DENS_TS = row0 * &
     & ( 1.-(Temp+A11)*(Temp-A14)**2/(A12*(Temp+A13) ) )
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      A =  alpha1         - alpha2*Temp    + alpha3*Temp**2 &
     &    -alpha4*Temp**3 + alpha5*Temp**4
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      B = -beta1          + beta2*Temp     - beta3*Temp**2
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      C = gamma1
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    ! 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.
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      END FUNCTION WATER_DENS_TS
    
      real(kind=ireals) FUNCTION DDENS_DTEMP(Temp,Sal)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    ! The function DDENS_DTEMP
    ! returns the derivative of water density
    ! on temperature, kg/(m**3*C)
      
    ! Input variables:
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    ! Temp --- water temperature, deg C
    ! Sal  --- water salinity,    kg/kg
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      implicit none
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      real(kind=ireals), intent(in):: Temp
      real(kind=ireals), intent(in):: Sal      
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      real(kind=ireals) Sal_g_kg,A,B,C
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      END FUNCTION DDENS_DTEMP
      
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      real(kind=ireals) FUNCTION DDENS_DSAL(Temp,Sal)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    ! The function DDENS_DSAL
    ! returns the derivative of water density
    ! on salinity, , kg/(m**3*kg/kg)
      
    ! Input variables:
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    ! Temp --- water temperature, deg C
    ! Sal  --- water salinity,    kg/kg
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      implicit none
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      real(kind=ireals), intent(in):: Temp
      real(kind=ireals), intent(in):: Sal
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      A =   alpha1         - alpha2*Temp    + alpha3*Temp**2 &
      &   - alpha4*Temp**3 + alpha5*Temp**4
    
      B = - beta1          + beta2*Temp     - beta3*Temp**2
    
      C = gamma1
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    ! Converting salinity units from kg/kg to g/kg
      Sal_g_kg = Sal*1.d+3
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      DDENS_DSAL = &
     & A + 1.5*B*Sal_g_kg**0.5 + 2.*C*Sal_g_kg
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      DDENS_DSAL = DDENS_DSAL*1.d+3
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      implicit none
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      
    ! 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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      
    
      DENS_HOSTETLER = row0*(1.-const*abs(temp-temp0)**const_power)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      
      END FUNCTION DENS_HOSTETLER
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      real(kind=ireals), intent(in):: Temp
      real(kind=ireals), intent(out):: ddens_dtemp, ddens_dsal
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      
    
      ddens_dtemp = row0*( - const_power*const*abs(temp-temp0)**(const_power-1)* &
    
      & sign(1._ireals,temp-temp0))
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      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_
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      END MODULE WATER_DENSITY
    
    
      MODULE PHYS_FUNC
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      contains
    
      real(kind=ireals) FUNCTION TURB_DENS_FLUX(tempflux,salflux,Temp,Sal)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    ! 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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
      TURB_DENS_FLUX = &
     & ddens_dtemp(Temp,Sal)*tempflux + &
     & ddens_dsal (Temp,Sal)*salflux
    
      END FUNCTION TURB_DENS_FLUX
    
     
    
    !  real(kind=ireals) FUNCTION dirdif()
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    !  use atmos, only:
    ! & shortwave
    ! implicit none
    
    !  real(kind=ireals) cloud
    !  real(kind=ireals) dirdif0,b,S0,sinh0
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    !  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)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      
    !  END FUNCTION dirdif
    
    
      FUNCTION SINH0(year,month,day,hour,phi)
    
    ! SINH0 is sine of solar angle ( = cosine of zenith angle)   
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      implicit none
      
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      
    
      integer(kind=iintegers), intent(in) :: year
      integer(kind=iintegers), intent(in) :: month
      integer(kind=iintegers), intent(in) :: day
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      real(kind=ireals)   , intent(in) :: hour
      real(kind=ireals)   , intent(in) :: phi
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      real(kind=ireals) delta
      real(kind=ireals) theta
      real(kind=ireals) pi
       real(kind=ireals) phi_rad
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     
      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)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      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) 
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
      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)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      
    ! 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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    !phase = 1 - liquid, 0 - ice
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    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
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    real(kind=ireals) FUNCTION ES(phase,t,p)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    !     ES is pressure of saturated water vapour, Pa    
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    implicit none
    
    real(kind=ireals) t,p,aw,bw,ai,bi,a,b
    integer(kind=iintegers) phase
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    !phase = 1 - liquid, 0 - ice
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    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)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    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 
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    real(kind=ireals), parameter :: dtdc = 66.7, pnorm = 101325
    real(kind=ireals), parameter :: Pa2dbar = 1.d-4
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    end select
    
    END FUNCTION MELTPNT
    
    FUNCTION FP_THETA(s,p,sat)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    !   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)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      
      FUNCTION WATER_FREEZE_MELT(temperature, grid_size, &
     & melting_point, switch)
     
    
    ! The function WATER_FREEZE_MELT checks, if the
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    ! 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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      
    
      integer(kind=iintegers), intent(in) :: switch ! +1 for water -> ice transition
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
                                       ! -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 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 ARRAYS_TURB, only : turb_type
      use WATER_DENSITY, only : SET_DDENS_DTEMP, SET_DDENS_DSAL
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
      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
    
      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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    ! 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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      
    
    ! External functions
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      
    
      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) :: zv, zm, x, y, x1, x2
    
      real(kind=ireals) :: u1, u2, v1, v2, Sal1, Sal2
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    ! 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)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      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) 
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
      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/&
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      & ((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