Skip to content
Snippets Groups Projects
turb_mod.f90 63.4 KiB
Newer Older
  • Learn to ignore specific revisions
  • MODULE TURB
    
    use NUMERICS, only : PROGONKA, IND_STAB_FACT_DB
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    use NUMERIC_PARAMS, only : vector_length, small_value
    
    use LAKE_DATATYPES, only : ireals, iintegers
    
    use DRIVING_PARAMS, only : missing_value
    
    SUBROUTINE K_EPSILON(ix, iy, nx, ny, year, month, day, &
    & hour, kor, a_veg, c_veg, h_veg, dt, &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    & b0, tau_air, tau_wav, tau_i, tau_gr, roughness, fetch)
          
    ! KT_eq calculates eddy diffusivity (ED) in water coloumn following parameterizations:
    
    ! 1) empirical profile of ED;
    ! 2) Semi-empirical formulations for ED;
    ! 2) k-epsilon parameterization for ED, including Kolmogorov relation.
    
    use COMPARAMS, only: &
    & num_ompthr
    
    
    use DRIVING_PARAMS , only : &
    & nstep_keps, &
    & turb_out, &
    & Turbpar, &
    & path, &
    & stabfunc, &
    & kepsbc, &
    & kwe, &
    & M, &
    & omp, &
    
    & eos, lindens, &
    & zserout
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    use ATMOS, only : &
    & uwind, vwind
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    use ARRAYS_TURB , only : &
    & S_integr_positive, &
    & S_integr_negative, &
    & Gen_integr, &
    & Seps_integr_positive, &
    & Seps_integr_negative, &
    & Geneps_integr, &
    & epseps_integr, &
    & eps_integr, &
    & E_integr, &
    & TKE_balance, &
    & eps_balance, &
    & Seps_integr_positive, &
    & Seps_integr_negative, &
    & Geneps_integr, &
    & epseps_integr, &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    & TKE_turb_trans_integr, &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    & veg_sw, &
    & E_it1, E_it2, &
    & E1, E2, &
    & eps_it1, eps_it2, &
    & eps1, eps2, &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    & E12, eps12, &
    & row, row2, &
    & k5_mid, k5, &
    & L, &
    & TF, &
    & k3_mid, &
    & WU_, WV_, &
    & GAMT, GAMU, GAMV, &
    
    & Gen, Gen_seiches, &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    & S, &
    & F, &
    & KT, &
    & TKE_turb_trans, &
    & Re, &
    & C1aup , &
    
    & Ri , Ri_bulk, &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    & k_turb_T_flux, &
    & H_entrainment, &
    & Buoyancy0, signwaveheight, &
    
    & Eseiches, &
    & TKE_budget_terms
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    use ARRAYS_GRID, only : &
    & ddz, ddz2, ddz05, ddz052, ddz054, &
    
    & dzeta_int, dzeta_05int, z_half
    
    use ARRAYS_BATHYM, only : &
    & dhw, dhw0, &
    
    & area_int, area_half, l1, h1, &
    
    & bathymwater, vol, botar, ls
    
    
    use ARRAYS_WATERSTATE, only : &
    & Tw1, Tw2, &
    & Sal1, Sal2, lamw, lamsal
    
    use ARRAYS, only : &
    & um, u2, u1, &
    & vm, v2, v1, &
    & nstep, &
    & time, &
    & uv
    
    use PHYS_CONSTANTS, only: &
    & row0, roa0, &
    & lamw0, &
    & cw, &
    & g, &
    & cw_m_row0, &
    
    & roughness0, niu_wat, &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    use TURB_CONST
    
    
    use PHYS_FUNC, only : &
    & H13
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    use WATER_DENSITY, only: &
    & water_dens_t_unesco, &
    & water_dens_ts, &
    & DDENS_DTEMP0, DDENS_DSAL0, &
    & SET_DDENS_DTEMP, SET_DDENS_DSAL
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    use INOUT_PARAMETERS, only : &
    & lake_subr_unit_min, &
    & lake_subr_unit_max
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    use NUMERICS, only : STEP
    
    use SEICHES_PARAM, only: SEICHE_ENERGY, C_Deff, cnorm
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    implicit none
    
    ! Input variables
    
    ! Reals
    
    real(kind=ireals), intent(in) :: dt
    real(kind=ireals), intent(in) :: kor
    real(kind=ireals), intent(in) :: h_veg, a_veg, c_veg
    real(kind=ireals), intent(in) :: hour
    real(kind=ireals), intent(in) :: b0
    real(kind=ireals), intent(in) :: tau_air, tau_wav, tau_i, tau_gr
    real(kind=ireals), intent(in) :: roughness
    real(kind=ireals), intent(in) :: fetch
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    ! Integers
    
    integer(kind=iintegers), intent(in) :: ix, iy
    integer(kind=iintegers), intent(in) :: nx, ny
    integer(kind=iintegers), intent(in) :: year, month, day
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
            
    ! Local variables     
    ! Reals
    
    real(kind=ireals) :: C_eps1(M), C_eps2(M), C_eps3(M)
    real(kind=ireals) :: CE(M), CEt(M)
    real(kind=ireals) :: lam_E(M+1), lam_eps(M+1)
    
    real(kind=ireals) :: a(vector_length)
    real(kind=ireals) :: b(vector_length)
    real(kind=ireals) :: c(vector_length)
    real(kind=ireals) :: d(vector_length)
    
    real(kind=ireals) :: AG(5)
    
    real(kind=ireals) :: dt05
    
    
    real(kind=ireals) :: dhw_keps, dhw0_keps
    real(kind=ireals) :: pi
    real(kind=ireals) :: ufr
    real(kind=ireals) :: dist_surf
    real(kind=ireals) :: ACC, ACC2, ACCk
    real(kind=ireals) :: dE_it, deps_it
    real(kind=ireals) :: xx(1:2), yy(1:2), zz
    real(kind=ireals) :: lm
    real(kind=ireals) :: ext_lamw
    real(kind=ireals) :: month_sec, day_sec, hour_sec
    real(kind=ireals) :: al_it
    real(kind=ireals) :: maxTKEinput, maxTKEinput_i, maxTKEinput_gr
    
    real(kind=ireals) :: FS, FTKES, TKES, DISS_TKES
    real(kind=ireals) :: FB, FTKEB, TKEB, DISS_TKEB
    
    real(kind=ireals) :: al, alft
    real(kind=ireals) :: GAMUN
    real(kind=ireals) :: urels, vrels
    real(kind=ireals) :: ksurf1(1:2), ksurf2(1:2)
    real(kind=ireals) :: Esurf1(1:2), Esurf2(1:2)
    real(kind=ireals) :: ceps3
    real(kind=ireals) :: vdamp
    real(kind=ireals) :: dE_it_max, deps_it_max
    
    real(kind=ireals) :: E_min, E_min_lamw0, E_min_Racr, eps_min, eps_min0
    
    real(kind=ireals) :: Buoyancy1, dBdz0, z0, bvf2
    real(kind=ireals) :: Eseiches_diss
    
    
    real(kind=ireals) :: tf_integr
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
         
    
    real(kind=ireals) :: GRADU, GRADV, GRADT, GAMVN
    real(kind=ireals) :: TFR
    real(kind=ireals) :: COGRTN, COGRUN, COGRVN
    real(kind=ireals) :: DTDZH
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    real(kind=ireals), allocatable :: work(:,:)
    
    real(kind=ireals), allocatable :: rhotemp(:), rhosal(:)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    ! Integers
    
    integer(kind=iintegers), parameter :: maxiter = 35
    
    integer(kind=iintegers) :: iter, iter_t
    integer(kind=iintegers) :: nl
    integer(kind=iintegers) :: i, j, k
    integer(kind=iintegers) :: keps_coef
    ! integer(kind=iintegers) :: stabfunc
    integer(kind=iintegers) :: time_deriv
    integer(kind=iintegers) :: i_entrain, i_entrain_old
    integer(kind=iintegers) :: nunit_ = lake_subr_unit_min
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    integer(kind=iintegers) :: grav_wave_Gill = 0
    
    integer(kind=iintegers) :: seiches_Goudsmit = 0
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    ! Logicals
    logical, allocatable :: init_keps(:,:)
    logical :: indstab
    logical :: ind_bound
    logical :: iterat
    logical :: firstcall
    logical :: cycle_keps
    logical :: next_tstep_keps
    logical :: smooth
    logical :: perdam
    logical :: semi_implicit_scheme
    logical :: contrgrad
    
    ! Characters
    character :: tp_name*10
    character :: numt*1
    
    
    data firstcall /.true./
    
    ! Externals
    
    real(kind=ireals), external :: DZETA
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    SAVE
    
    
    firstcall_if : if (firstcall) then
    
      if (turb_out%par == 1) then
    
        call CHECK_UNIT(lake_subr_unit_min,lake_subr_unit_max,nunit_)
        open (nunit_,file=path(1:len_trim(path))//'results/debug/err_keps_iter.dat', &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
        & status = 'unknown')
    
        write(unit=nunit_,fmt=*) 'Errors of iterative process in k-epsilon solver'
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      endif
      
      allocate (init_keps(1:nx, 1:ny) )
      init_keps(:,:) = .false.
           
      month_sec   = 30.*24.*60.*60.
      day_sec     = 24*60.*60.
      hour_sec    = 60.*60.
    
    
      pi          = 4.*atan(1.e0_ireals)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
      AL    = g/row0
      ALFT  = 1.d0 
      
    
      semi_implicit_scheme = .true. ! No iterations in k-epsilon, semi-implicit scheme
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      al_it   = 0.2 !0.8
      ACC     = 1.d-20 !1.d-20
      ACC2    = 1.d-20 !1.d-20
      ACCk    = 1.d-20 !1.d-20
      knum    = 0.
    
      z0 = 1.d-2
      
    
      eps_min0    = 1.d-12 ! The range for dissipation rate values 
                           ! 10**(-11) - 10**(-6) m**2/s**3 (Wuest and Lorke, 2003)
      E_min_lamw0 = sqrt(eps_min0*1.E-5*lamw0*(cw_m_row0*CEt0)**(-1) ) ! This expression ensures that
                                                                       ! eddy diffusivity in the decaying
                                                                       ! turbulence (E -> Emin, eps -> eps_min)
                                                                       ! is much less than the molecular diffusivity
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      dE_it_max = 1.d-11 !1.d-9
      deps_it_max = 1.d-15 !1.d-12
    
      smooth = .false.
      perdam = .false.
      vdamp  = 1.d-3 !1.d-3
     
      contrgrad = .false.
    
      ! AG(1)   = 0._ireals
      ! AG(2)   = 0._ireals
      ! AG(3)   = 0._ireals
      ! AG(4)=1. !0.
      ! AG(5)=1. !0.
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      keps_coef = 1
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      ! keps_coef = 1 - standard empirical coefficients in epsilon-equation
      ! keps_coef = 2 - the coefficient by (Aupoix et al., 1989) is used
        
      ! stabfunc = 1  - stability functions CE and CEt are constants
      ! stabfunc = 2  - stability functions CE and CEt are dependent 
      ! on shear and buoyancy (Canuto et al., 2001)
      ! stabfunc = 3  - stability functions CE and CEt are dependent
      ! only on buoyancy (Galperin et al., 1988)
             
      ! FRICTION: VEGETATION
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
      if (h_veg>0) then
        print*, 'Vegetation friction parameterization &
        & is not operational currently: STOP'
        STOP
      endif
    
      veg_sw = 0.
      do i = 1, M+1
        if (h1-dzeta_int(i)*h1 < h_veg) veg_sw(i) = 1.
      enddo
    
    endif firstcall_if
    
    
    if (l1 > 0._ireals) then
      !The minimal values for E and epsilon are determined from two conditions:
      !1.corresponding turbulent diffusion coefficient is much less than molecular coefficient
      !2.free convection occurs in quiscent state under Ra>Racr
      eps_min = 1.E-5*(lamw0/cw_m_row0)**2*Racr*niu_wat/(maxval(ddz)*h1)**4
      E_min_Racr = eps_min*(maxval(ddz)*h1)**2*(CEt0*Racr*niu_wat*lamw0/cw_m_row0)**(-0.5)
      !E_min = max(E_min_lamw0,E_min_Racr)
      E_min = E_min_Racr
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    dt05 = 0.5d0*dt
    
    allocate (work(1:M+1,1:2))
    
    allocate (rhotemp(1:M), rhosal(1:M))
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    ! Implementation of the Crank-Nicolson numerical scheme
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    ! for k-epsilon parameterization using simple iterations to
    ! solve the set of nonlinear finite-difference equations
    
    !  k_turb_T_flux used here from the previous time step
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      ! Mixed layer defined from the gradient of heat turbulent flux
      !i_entrain = M
      !do i = M-1, 1, -1
      !  if (abs(k_turb_T_flux(i+1)-k_turb_T_flux(i))/abs(k_turb_T_flux(i+1) + 1.d-20) &
      !  & > 0.2 ) then ! 0.2 - quite arbitrary value 
      !    i_entrain = i
      !    exit
      !  endif
      !enddo
    
    
      ! Seiche energy evolution, TKE production due to seiche dissipation
      if (seiches_Goudsmit == 1) then
        call SEICHE_ENERGY &
    
        & (M,bathymwater,ls,tau_air-tau_wav,sqrt(uwind*uwind+vwind*vwind), &
        & vol,dt,Eseiches,Eseiches_diss)
    
        !A part of seiche energy dissipation feeding TKE (Goudsmit et al. 2002, eg. (24))
        do i = 1, M
          bvf2 = max(AL*(row(i+1) - row(i))/(h1*ddz(i)),0._ireals)
          Gen_seiches(i) = (1. - 10.*sqrt(C_Deff))*bvf2/(row0*cnorm*botar)*&
          & (area_int(i) - area_int(i+1))/area_half(i)*Eseiches_diss
        enddo
      else
        Gen_seiches(1:M) = 0.
      endif
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      ! Mixed layer defined from TKE gradient
      i_entrain = minloc(E1(1:M),1)
    
      !do i = 1, M
      !  if (E1(i) < 1.E-6) then
      !    i_entrain = i
      !    exit
      !  endif
      !enddo
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      !do i = i_entrain, 1, -1
      !  if (E1(i)/E1(i_entrain) > 5.E+0_ireals) then
      !    i_entrain = i
      !    exit
      !  endif
      !enddo
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      H_entrainment = dzeta_05int(i_entrain) * h1
    
    
      ! Bulk Richardson number for mixed layer
      Ri_bulk = g/row0*(row(i_entrain)-row(1))*(dzeta_int(i_entrain)*h1)/ &
    
      & ( (u1(i_entrain) - u1(1))**2 + (v1(i_entrain) - v1(1))**2 + small_value)
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      ! Significant wave height
      if (kepsbc%par == 2 .or. kepsbc%par == 3) then
        signwaveheight = H13(tau_air,fetch)
      endif
    
      dhw_keps  = dhw  / real(nstep_keps%par)
      dhw0_keps = dhw0 / real(nstep_keps%par)
    
      ! This time interpolation ensures correspondence between mean energy
      ! viscous dissipation in Crank-Nickolson scheme for momentum equation and TKE shear production
      ! when semi-implicit scheme is used (semi_implicit_scheme = .true.)
      do i = 1, M+1
        Um(i) = 0.5d0 * (u2(i) + u1(i))
        Vm(i) = 0.5d0 * (v2(i) + v1(i))
      enddo
    
    
      ! The TKE injection from waves breaking is limited
      ! by the wind kinetic energy input to waves development
    
      maxTKEinput = tau_wav/row0*sqrt(uwind*uwind + vwind*vwind)
    
      ! maxTKEinput_i = tau_i/row0*sqrt(Um(1)*Um(1)+Vm(1)*Vm(1))
      ! maxTKEinput_gr = tau_gr/row0*sqrt(Um(M+1)*Um(M+1)+Vm(M+1)*Vm(M+1))
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      
      next_tstep_keps = .false.
          
      keps_mode: do while (.not.next_tstep_keps)
    
    
      if (.not.init_keps(ix,iy)) then
    
        ! Initialization mode of k-epsilon parameterization: time derivatives are neglected  
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
        time_deriv = 0
      else
    
        ! Evolution mode of k-epsilon parameterization: full equations are solved
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
        time_deriv = 1
      endif
      
      do i = 1, M+1
        E_it1(i)   = E1(i)
        eps_it1(i) = eps1(i)
    
        k2_mid(i)  = 0._ireals
        k2(i)      = 0._ireals
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      enddo
    
      iter_t = 0
    
      cycle_keps = .true.
      
      if (semi_implicit_scheme) then
        iterat = .false.
      else
        iterat = .true.
      endif  
    
    
      ! Setting density derivatives on temperature and salinity
    
        rhotemp(i) = SET_DDENS_DTEMP(eos%par,lindens%par,0.5*(Tw1(i)+Tw1(i+1)), &
                                                         0.5*(Sal1(i)+Sal1(i+1)))
        rhosal(i)  = SET_DDENS_DSAL (eos%par,lindens%par,0.5*(Tw1(i)+Tw1(i+1)), &
                                                         0.5*(Sal1(i)+Sal1(i+1)))
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      iterat_keps: do while (cycle_keps)
    
    
        ! The cycle implements iterations to
        ! solve the set of nonlinear finite-difference equations
        ! of k-epsilon parameterization
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
        do i = 1, M+1
          E12(i)   = 0.5d0*(E1(i)   + E_it1(i))
          eps12(i) = 0.5d0*(eps1(i) + eps_it1(i))
        enddo
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     
    
        ! E12   =  E_it1
        ! eps12 =  eps_it1
        
        ! Calculating the stability functions
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
        do i = 1, M
          if (stabfunc%par == 1) then
    
            lam_eps(i) = lam_eps0
            lam_E  (i) = lam_E0
            CE     (i) = CE0
            CEt    (i) = CEt0
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          else
    
            work(i,1) = E12(i)*E12(i)/((eps12(i) + ACC2)*(eps12(i) + ACC2))* & 
    
            &        (rhotemp(i) * & ! Assuming Pr = Sc
    
            &        ( Sal1(i+1) - Sal1(i)) ) / &
            &        (h1*ddz(i)) *AL ! ~ bouyancy production of TKE
            work(i,2) = E12(i)*E12(i)/((eps12(i) + ACC2)*(eps12(i) + ACC2))* &
            & ( (U1(i+1) - U1(i))*(U1(i+1) - U1(i)) + &
            &   (V1(i+1) - V1(i))*(V1(i+1) - V1(i)) ) / &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
            &   (h1*h1*ddz(i)*ddz(i)) ! ~ shear production of TKE
            if (stabfunc%par == 2) then
              CE(i)  = CE_CANUTO (work(i,2), work(i,1))
              CEt(i) = CEt_CANUTO(work(i,2), work(i,1))
            elseif (stabfunc%par == 3) then
    
              CE(i)  = sqrt(2.d0)*CL_K_KL_MODEL*SMOMENT_GALPERIN(work(i,1))
              CEt(i) = sqrt(2.d0)*CL_K_KL_MODEL*SHEAT_GALPERIN(work(i,1))
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
            endif
            lam_eps(i) = CE(i)/sigmaeps
            lam_E  (i) = CE(i)/sigmaE
          endif
        enddo
    
    
        ! Calculating eddy diffusivities for kinetic energy and its dissipation
        ! in the middle between half levels
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
        do i = 2, M
    
          ! E_mid=0.5d0*(E12(i-1)+E12(i))
          ! eps_mid=0.5d0*(eps12(i-1)+eps12(i))
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          work(i,1) = INTERPOLATE_TO_INTEGER_POINT(E12(i-1),E12(i),ddz(i-1),ddz(i)) ! interpolated TKE
          work(i,2) = INTERPOLATE_TO_INTEGER_POINT(eps12(i-1),eps12(i),ddz(i-1),ddz(i)) ! interpolated dissipation
    
          xx(1) = INTERPOLATE_TO_INTEGER_POINT(lam_eps(i-1),lam_eps(i),ddz(i-1),ddz(i))
          xx(2) = INTERPOLATE_TO_INTEGER_POINT(lam_E(i-1),lam_E(i),ddz(i-1),ddz(i))
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          if (iter_t == 0) then
    
            k2_mid(i) = work(i,1)*work(i,1)/(work(i,2) + ACCk)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          else
    
            k2_mid(i) = k2_mid(i)*al_it + work(i,1)*work(i,1)/(work(i,2) + ACCk)*(1-al_it)
          endif     
          k5_mid(i) = niu_wat + max(k2_mid(i)*xx(1), min_visc/sigmaeps)
          k3_mid(i) = niu_wat + max(k2_mid(i)*xx(2), min_visc/sigmaE)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
        enddo
        do i = 1, M
          if (iter_t == 0) then
    
            k2(i)  = max(E12(i)*E12(i)/(eps12(i) + ACCk), min_visc/CE(i))
            k2t(i) = max(E12(i)*E12(i)/(eps12(i) + ACCk), min_diff/CEt(i))
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          else
    
            k2(i)  = k2(i) *al_it + &
            & max(E12(i)*E12(i)/(eps12(i) + ACCk), min_visc/CE(i))*(1-al_it)
            k2t(i) = k2t(i)*al_it + &
            & max(E12(i)*E12(i)/(eps12(i) + ACCk), min_diff/CEt(i))*(1-al_it)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          endif
    
          k5(i) = niu_wat + k2(i)*lam_eps(i)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          l(i) = E12(i)**(1.5d0)/(eps12(i) + ACC)
    
          TF(i) = sqrt(E12(i))/(l(i) + ACC)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
        enddo
        do i = 2, M-1
          WU_(i) = (E_it1(i+1)-E_it1(i-1))*(U2(i)-U2(i-1))/ &
          & (h1*h1*ddz052(i)*ddz(i-1) )
          WV_(i) = (E_it1(i+1)-E_it1(i-1))*(V2(i)-V2(i-1))/ &
          & (h1*h1*ddz052(i)*ddz(i-1) )
        enddo
     
        if (contrgrad) then
          do i = 2, M
            GRADT = (Tw1(i+1)-Tw1(i-1))/(h1*ddz2(i-1))
            GRADU = (U2(i)-U2(i-1))/(h1*ddz(i-1) )
            GRADV = (V2(i)-V2(i-1))/(h1*ddz(i-1) )
            GAMUN = - CONUV*(WU_(i+1)-WU_(i-1))/(h1*ddz2(i-1))
            GAMVN = - CONUV*(WV_(i+1)-WV_(i-1))/(h1*ddz2(i-1))
            COGRTN = DTDZH * 0.1
            TFR=E_it1(i)/(l(i)*l(i)+ACC) + acc
            COGRUN = GAMUN/TFR
            COGRVN = GAMVN/TFR 
            if(gradT < 0.) then
              GAMT(i) = - AG(3)*cogrtn
              GAMU(i) = - AG(1)*cogrun
              GAMV(i) = - AG(2)*cogrvn
            else
              GAMT(i) = - AG(3)*(AL*GRADT**2+CON0*TFR*COGRTN) / &
              &            (AL*GRADT+CON0*TFR)
              GAMU(i) = - AG(1)*(AL*(CON2*GRADT+(CON2-1.)*GAMT(i))*GRADU + &
              &            CON1*TFR*COGRUN)/(AL*GRADT+CON1*TFR)
              GAMV(i) = - AG(2)*(AL*(CON2*GRADT+(CON2-1.)*GAMT(i))*GRADV + &
              &            CON1*TFR*COGRVN)/(AL*GRADT+CON1*TFR)
            endif
          enddo
          GAMT(1)=0.
          GAMU(1)=0.
          GAMV(1)=0.
        else
          do i = 1, M+1
            GAMT(i) = 0.
            GAMU(i) = 0.
            GAMV(i) = 0.
          enddo
        endif
    
        do i = 1, M
          Gen(i) = ((Um(i+1) - Um(i))*(Um(i+1) - Um(i) + h1*ddz(i)*GAMU(i)) + &
          &        (Vm(i+1) - Vm(i))*(Vm(i+1) - Vm(i) + h1*ddz(i)*GAMV(i))) / &
          &        (h1**2*ddz(i)**2)*CE(i)
    
          ! S(i) = - ( 0.5d0*(row(i+1) + row2(i+1) - row(i) - row2(i)) / &
          ! &       (h1*ddz(i)) + GAMT(i))*ALFT*AL*CEt(i)
    
          S(i) = - ( 0.5d0*( rhotemp(i) * & ! Assuming Pr = Sc
    
          &        (Tw1(i+1) + Tw2(i+1) - Tw1(i) - Tw2(i)) + &
    
          &        (Sal1(i+1) + Sal2(i+1) - Sal1(i) - Sal2(i)) ) / &
          &        (h1*ddz(i)) + GAMT(i) )*ALFT*AL*CEt(i)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          F(i) = Gen(i) + S(i)
    
          ! Adding shear production of turbulence by gravitational waves in stable stratification
          Gen(i) = Gen(i) - grav_wave_Gill*grav_wave_Gill_const*STEP(-S(i))*S(i)*CE(i)/CEt(i) 
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
        enddo
        
        if (iter_t == 0) then
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          dBdz0 = (Buoyancy1 - Buoyancy0)/(0.5*h1*ddz(1))
        endif
          
    
        ! Solution of the equation for TKE (turbulent kinetic energy)
    
        if_tkebc : if (kepsbc%par <= 4) then
          !Neumann boundary conditions for TKE
          !Boundary condition at the top
    
          if (l1 == 0._ireals) then
    
            !Open water
            if (kepsbc%par < 4) then
              FTKES = TKE_FLUX_SHEAR(tau_air,kwe%par,maxTKEinput,kepsbc%par)
            elseif (kepsbc%par == 4) then ! The new boundary condition, 
                                          ! implemented currently only for open water case
              FTKES = TKE_FLUX_BUOY(Buoyancy0,roughness0,H_entrainment)
            endif
    
          elseif (l1 > 0._ireals .and. l1 <= L0) then
    
            !Fractional ice
            FTKES =        b0*TKE_FLUX_SHEAR(tau_i,0._ireals,maxTKEinput,1) + &
            &    (1.d0-b0)*TKE_FLUX_SHEAR(tau_air,kwe%par,maxTKEinput,kepsbc%par)
          elseif (l1 > L0) then
            !Complete ice cover
            FTKES = TKE_FLUX_SHEAR(tau_i,0._ireals,maxTKEinput,1)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          endif
    
          xx(1) = 0.5*(area_int(2)/area_half(1)*k3_mid(2)*dt / &
          &  (h1**2*ddz(1)*ddz05(1)) + &
          &  time_deriv*dzeta_05int(1)*dhw_keps/(h1*ddz05(1) ) - &
          &  time_deriv*dhw0_keps/(h1*ddz05(1) ) )
          yy(1) = dt/(h1*ddz(1))
          b(1) = - xx(1)
          c(1) = - (xx(1) + time_deriv*1.d0) - (abs(S(1)) - S(1))/(TF(1) + ACC)*dt05 - TF(1)*dt
          d(1) = - time_deriv*E1(1) - xx(1)*(E1(2) - E1(1)) - &
          & area_int(1)/area_half(1)*yy(1)*FTKES - &
          & k2t(1)*(S(1) + abs(S(1)))*dt05 - dt*(k2(1)*Gen(1) + Gen_seiches(1))!+eps_it1(i)*dt
          ! Boundary condition at the bottom
          FTKEB = - TKE_FLUX_SHEAR(tau_gr,0._ireals,maxTKEinput,1)
          xx(2) = 0.5*(area_int(M)/area_half(M)*k3_mid(M)*dt / &
          & (h1**2*ddz(M)*ddz05(M-1) ) - &
          & time_deriv*dzeta_05int(M)*dhw_keps/(h1*ddz05(M-1) ) + &
          & time_deriv*dhw0_keps/(h1*ddz05(M-1) ) )
          yy(2) = dt/(h1*ddz(M))
          a(M) = - xx(2)
          c(M) = - (xx(2) + time_deriv*1.d0) - (abs(S(M)) - S(M))/(TF(M) + ACC)*dt05 - &
          & TF(M)*dt
          d(M) = - time_deriv*E1(M) - xx(2)*(E1(M-1) - E1(M)) + &
          & area_int(M+1)/area_half(M)*yy(2)*FTKEB - &
          & k2t(M)*(S(M) + abs(S(M)))*dt05 - dt*(k2(M)*Gen(M) + Gen_seiches(M)) !+eps_it1(i)*dt
        elseif (kepsbc%par == 5) then
          !Dirichlet boundary condition for TKE
          !Boundary condition at the top
    
          if (l1 == 0._ireals) then
    
          elseif (l1 > 0._ireals .and. l1 <= L0) then
    
            !Fractional ice
            TKES = b0*TKE_WALL_SHEAR(tau_i) + (1.-b0)*TKE_WALL_SHEAR(tau_air)
          elseif (l1 > L0) then
            !Complete ice cover
            TKES = TKE_WALL_SHEAR(tau_i)
          endif
          b(1) = 0.
          c(1) = 1.
          d(1) = TKES
          ! Boundary condition at the bottom
          TKEB = TKE_WALL_SHEAR(tau_gr)
          a(M) = 0.
          c(M) = 1.
          d(M) = TKEB
        endif if_tkebc
    
    
        ! The coefficients of linear algebraic equations in the internal points of the mesh
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
        do i = 2, M-1
          a(i) = time_deriv*dzeta_05int(i)*dhw_keps/(ddz054(i)*h1) - &
          & time_deriv*dhw0_keps/(ddz054(i)*h1) - &
          & area_int(i)/area_half(i)*k3_mid(i)*dt/(h1**2*ddz(i)*ddz2(i-1) )
    
          b(i) = -time_deriv*dzeta_05int(i)*dhw_keps/(ddz054(i)*h1) + &
          & time_deriv*dhw0_keps/(ddz054(i)*h1) - &
          & area_int(i+1)/area_half(i)*k3_mid(i+1)*dt/(h1**2*ddz(i)*ddz2(i) )
    
    
          c(i) = a(i) + b(i) - time_deriv*1.d0-(abs(S(i))-S(i))/(TF(i)+ACC)*dt05 - TF(i)*dt
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
          d(i) = -time_deriv*E1(i)-k2t(i)*(S(i)+abs(S(i)))*dt05 - dt*(k2(i)*Gen(i) + Gen_seiches(i))
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          d(i) = d(i)-dt*(h1**2*ddz(i)*ddz2(i))**(-1)* &
          & (area_int(i+1)/area_half(i)*k3_mid(i+1)*(E1(i+1)-E1(i)))   
          d(i) = d(i)+dt*(h1**2*ddz(i)*ddz2(i-1))**(-1)* &
          & (area_int(i)/area_half(i)*k3_mid(i)*(E1(i)-E1(i-1)))
          d(i) = d(i)-time_deriv*dzeta_05int(i)*dhw_keps* &
          & (ddz054(i)*h1)**(-1)*(E1(i+1)-E1(i-1)) 
          d(i) = d(i)+time_deriv*dhw0_keps* &
          & (ddz054(i)*h1)**(-1)*(E1(i+1)-E1(i-1)) 
        enddo
        ind_bound = .false.
        call IND_STAB_FACT_DB(a,b,c,1,M,indstab,ind_bound)
    
        if (indstab .eqv. .false.) then
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          do i = 2, M-1
            a(i) = - k3_mid(i)*dt/(h1**2*ddz(i)*ddz2(i-1) )
            b(i) = - k3_mid(i+1)*dt/(h1**2*ddz(i)*ddz2(i) )
          enddo
        endif
         
        call PROGONKA (a,b,c,d,E_it2,1,M)
        E_it2(M+1) = 0.
                    
        call CHECK_MIN_VALUE(E_it2, M, E_min) 
    
    
        dE_it = maxval(abs(E_it1-E_it2))
    
        ! C1 is given following Satyanarayana et al., 1999; Aupoix et al., 1989
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
        if (keps_coef == 1) then
          do i = 1, M
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
            yy(1) = 0.5*(1. + sign(1._ireals,S(i)))
            zz    = 0.5*(1. - sign(1._ireals,S(i)))
            ceps3 = ceps3_stable*zz + ceps3_unstable*yy(1)
    
            ! ceps3 = 1.14d0
            ! ceps3 < -0.4 should be used for stable stratification (Burchard, 2002)
            ! ceps3 = 1. ! following (Kochergin and Sklyar, 1992)
            ! ceps3 = 1.14d0 ! Baum and Capony (1992)
            ! Other values for ceps3:
            ! 0.<ceps3<0.29 for stable and ceps3 = 1.44 for unstable conditions (Rodi, 1987);
            ! ceps3 >= 1 (Kochergin and Sklyar, 1992)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
            C_eps1(i) = ceps1
            C_eps2(i) = ceps2
            C_eps3(i) = ceps3
          enddo
        elseif (keps_coef == 2) then
          do i = 1, M 
    
            Re(i) = ACC + (2*E_it1(i)/3.)**2/(eps_it1(i)*niu_wat + ACC)
    
            C1aup(i) = Co/(1.d0 + 0.69d0*(2.d0 - Co)/sqrt(Re(i)))
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
            C_eps1(i) = C1aup(i)
            C_eps2(i) = C1aup(i)
            C_eps3(i) = C1aup(i)
          enddo  
        endif
    
    
        ! The solution of the equation for dissipation rate
        
        ! dist_surf is the distance from the surface, 
        ! where the b.c. for epsilon are formulated
        
        ! dist_surf = 0._ireals !0.5d0*ddz(1)*h1
        if_dissbc: if (kepsbc%par <= 4) then
          !Neumann boundary conditions for TKE dissipation
          !Boundary condition at the top
    
          if (l1 == 0._ireals) then
    
            if (kepsbc%par < 4) then
              ! FS = CE0**(0.75d0)*k5(1)*E_it1(1)**(1.5d0)/kar* &
              ! & (roughness0 + dist_surf)**(-2)
              ! FS = DISSIPATION_FLUX_SHEAR(tau_air,roughness0)
    
              ! Extrapolation of the turbulent diffusivity to the boundary level, assuming
              ! logarithmic profile
              ksurf1(1) = max(niu_wat, CE(1)*k2(1)-kar*sqrt(tau_air/row0)*0.5*ddz(1)*h1)
              
              ! ksurf1(1) = kar*sqrt(tau_air/row0)*z0 !roughness0
              ! ksurf1(1) = CE(1)*k2(1)
              Esurf1(1) = E_it1(1) + 0.5*FTKES*sigmaE*h1*ddz(1)/ksurf1(1)
              FS = DISSIPATION_FLUX_SHEAR(tau_air,ksurf1(1),Esurf1(1),z0, &
              & kwe%par,fetch,maxTKEinput,kepsbc%par,signwaveheight)
            elseif (kepsbc%par == 4) then ! The new boundary condition, 
                                          ! implemented currently only for open water case
              Esurf1(1) = E_it1(1)
              FS = DISSIPATION_FLUX_BUOY(Buoyancy0,roughness0,H_entrainment,dBdz0,Esurf1(1))
            endif
    
          elseif (l1 > 0._ireals .and. l1 <= L0) then
    
            ! FS = CE0**(0.75d0)*k5(1)*E_it1(1)**(1.5d0)/kar* &
            ! & (0.0._ireals + dist_surf)**(-2) !0.01 m - roughness of ice 
            ! FS =     b0*DISSIPATION_FLUX_SHEAR(tau_i,1.d-3) + & ! 1.d-3 roughness of ice
            ! & (1.d0-b0)*DISSIPATION_FLUX_SHEAR(tau_air,roughness0)
    
            ! Extrapolation of the turbulent diffusivity to the boundary level, assuming
            ! logarithmic profile
            ksurf1(1) = max(niu_wat, CE(1)*k2(1) - kar*sqrt(tau_air/row0)*0.5*ddz(1)*h1)
            ksurf2(1) = max(niu_wat, CE(1)*k2(1) - kar*sqrt(tau_i/row0)*0.5*ddz(1)*h1)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
            
    
            ! ksurf1(1) = kar*sqrt(tau_air/row0)*z0 !*roughness0
            ! ksurf2(1) = kar*sqrt(tau_i/row0)*z0 !1.d-3 is the roughness of ice bottom
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
            ! ksurf1(1) = CE(1)*k2(1)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
            Esurf1(1) = E_it1(1) + 0.5*FTKES*sigmaE*h1*ddz(1)/ksurf1(1)
    
            Esurf2(1) = E_it1(1) + 0.5*FTKES*sigmaE*h1*ddz(1)/ksurf2(1)
            FS =     b0*DISSIPATION_FLUX_SHEAR(tau_i,ksurf2(1),Esurf2(1),z0, &
            & 0._ireals,fetch,maxTKEinput,1,signwaveheight) + &
            & (1.d0-b0)*DISSIPATION_FLUX_SHEAR(tau_air,ksurf1(1),Esurf1(1),z0, &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
            & kwe%par,fetch,maxTKEinput,kepsbc%par,signwaveheight)
    
          elseif (l1 > L0) then
            ! FS = DISSIPATION_FLUX_SHEAR(tau_i,1.d-3)
    
            ! Extrapolation of the turbulent diffusivity to the boundary level, assuming
            ! logarithmic profile
            ksurf2(1) = max(niu_wat, CE(1)*k2(1) - kar*sqrt(tau_i/row0)*0.5*ddz(1)*h1)
            
            ! ksurf2(1) = kar*sqrt(tau_i/row0)*z0 !1.d-3 is the roughness of ice bottom
            ! ksurf2(1) = CE(1)*k2(1)
            Esurf2(1) = E_it1(1) + 0.5*FTKES*sigmaE*h1*ddz(1)/ksurf2(1)
            FS = DISSIPATION_FLUX_SHEAR(tau_i,ksurf2(1),Esurf2(1),z0, &
            & 0._ireals,fetch,maxTKEinput,1,signwaveheight)
          endif 
          xx(1) = 0.5*(area_int(2)/area_half(1)*k5_mid(2)*dt / &
          & (h1**2*ddz(1)*ddz05(1)) + &
          & time_deriv*dzeta_05int(1)*dhw_keps/(h1*ddz05(1) ) - &
          & time_deriv*dhw0_keps/(h1*ddz05(1) ) )
          yy(1) = dt/(h1*ddz(1))
          b(1) = - xx(1)
          c(1) = - (xx(1) + time_deriv*1.d0) - C_eps2(1)*TF(1)*dt + &
          & C_eps3(1)*(- abs(S(1)) + S(1))/(TF(1) + ACC)*dt05
          d(1) = - time_deriv*eps1(1) - xx(1)*(eps1(2) - eps1(1)) - &
          & area_int(1)/area_half(1)*yy(1)*FS - C_eps3(1) * &
          & (abs(S(1)) + S(1))*TF(1)*k2t(1)*dt05 - &
          & C_eps1(1)*TF(1)*dt*(k2(1)*Gen(1) + Gen_seiches(1)) !+eps_it1(i)*dt
     
          ! Boundary condition at the bottom
          ! dist_surf = 0._ireals ! ddz(M)/2*h1
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          
    
          ! FB = - CE0**(0.75d0)*k5(M)* &
          ! & E_it1(M)**(1.5d0)/kar*(0.0._ireals + dist_surf)**(-2) 
          ! FB = - DISSIPATION_FLUX_SHEAR(tau_gr,1.d-3) ! 1.d-3 m - rougnness of bottom
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          
    
          ! Extrapolation of the turbulent diffusivity to the boundary level, assuming
          ! logarithmic profile
          ksurf2(2) = max(niu_wat, CE(M)*k2(M) - kar*sqrt(tau_gr/row0)*0.5*ddz(M)*h1)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     
    
          ! ksurf2(2) = kar*sqrt(tau_gr/row0)*z0 !1.d-3 is the roughness of bottom
          ! ksurf2(2) = CE(M)*k2(M)
          Esurf2(2) = E_it1(M) - 0.5*FTKEB*sigmaE*h1*ddz(M)/ksurf2(2)
          FB = - DISSIPATION_FLUX_SHEAR(tau_gr,ksurf2(2),Esurf2(2),z0, &
          & 0._ireals,fetch,maxTKEinput,1,signwaveheight)
          xx(2) = 0.5*(area_int(M)/area_half(M)*k5_mid(M)*dt / &
          & (h1**2*ddz(M)*ddz05(M-1) ) - &
          & time_deriv*dzeta_05int(M)*dhw_keps/(h1*ddz05(M-1) ) + &
          & time_deriv*dhw0_keps/(h1*ddz05(M-1) ))
          yy(2) = dt/(h1*ddz(M))
          a(M) = - xx(2)
          c(M) = - (xx(2) + time_deriv*1.d0) - C_eps2(M)*TF(M)*dt + &
          & C_eps3(M)*( - abs(S(M)) + S(M))/(TF(M) + ACC)*dt05
          d(M) = - time_deriv*eps1(M) - xx(2)*(eps1(M-1) - eps1(M)) + &
          & area_int(M+1)/area_half(M)*yy(2)*FB - C_eps3(M) * &
          & (abs(S(M)) + S(M))*TF(M)*k2t(M)*dt05 - &
          & C_eps1(M)*TF(M)*dt*(k2(M)*Gen(M) + Gen_seiches(M))
        elseif (kepsbc%par == 5) then
          !Dirichlet boundary conditions for TKE dissipation
          !Boundary condition at the top
          dist_surf = 0.5*ddz(1)*h1
    
          if (l1 == 0._ireals) then
    
            !Open water
            DISS_TKES = DISSIPATION_WALL_SHEAR(tau_air,Buoyancy0,dist_surf)
    
          elseif (l1 > 0._ireals .and. l1 <= L0) then
    
            !Fractional ice
            DISS_TKES = b0    *DISSIPATION_WALL_SHEAR(tau_i  ,Buoyancy0,dist_surf) + &
            &          (1.-b0)*DISSIPATION_WALL_SHEAR(tau_air,Buoyancy0,dist_surf)
          elseif (l1 > L0) then
            !Complete ice cover
            DISS_TKES = DISSIPATION_WALL_SHEAR(tau_air,Buoyancy0,dist_surf)
          endif
          b(1) = 0.
          c(1) = 1.
          d(1) = DISS_TKES
          !Boundary condition at the bottom
          dist_surf = 0.5*ddz(M)*h1
          DISS_TKEB = DISSIPATION_WALL_SHEAR(tau_gr,0._ireals,dist_surf)
          a(M) = 0.
          c(M) = 1.
          d(M) = DISS_TKEB
        endif if_dissbc
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
        do i = 2, M-1
          a(i) = time_deriv*dzeta_05int(i)*dhw_keps/(ddz054(i)*h1) - &
          & time_deriv*dhw0_keps/(ddz054(i)*h1) - &
          & area_int(i)/area_half(i)*k5_mid(i)*dt/(h1**2*ddz(i)*ddz2(i-1) )
    
          b(i) = -time_deriv*dzeta_05int(i)*dhw_keps/(ddz054(i)*h1) + &
          & time_deriv*dhw0_keps/(ddz054(i)*h1) - &
          & area_int(i+1)/area_half(i)*k5_mid(i+1)*dt/(h1**2*ddz(i)*ddz2(i) )
    
          c(i) = a(i) + b(i) - time_deriv*1.d0 - C_eps2(i)*TF(i)*dt + &
    
          & C_eps3(i)*(-abs(S(i))+S(i))/(TF(i)+ACC)*dt05
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
          d(i) = - time_deriv*eps1(i) - C_eps3(i)*(abs(S(i)) + S(i))*TF(i)*k2t(i)*dt05 - &
    
          & C_eps1(i)*TF(i)*dt*(k2(i)*Gen(i) + Gen_seiches(i))
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          d(i) = d(i)-dt*(h1**2*ddz(i)*ddz2(i))**(-1)* &
          & (area_int(i+1)/area_half(i)*k5_mid(i+1)*(eps1(i+1)-eps1(i)))
          d(i) = d(i)+dt*(h1**2*ddz(i)*ddz2(i-1))**(-1)* &
          & (area_int(i)/area_half(i)*k5_mid(i)*(eps1(i)-eps1(i-1)))
          d(i) = d(i)-time_deriv*dzeta_05int(i)*dhw_keps* &
          & (ddz054(i)*h1)**(-1)*(eps1(i+1)-eps1(i-1))
          d(i) = d(i)+time_deriv*dhw0_keps* &
          & (ddz054(i)*h1)**(-1)*(eps1(i+1)-eps1(i-1))
        enddo
        ind_bound = .false.
        call IND_STAB_FACT_DB(a,b,c,1,M,indstab,ind_bound)
    
        if (indstab .eqv. .false.) then
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          do i = 2, M-1
            a(i) = - k5_mid(i)*dt/(h1**2*ddz(i)*ddz2(i-1))
            b(i) = - k5_mid(i+1)*dt/(h1**2*ddz(i)*ddz2(i))
          enddo
        endif
          
        call PROGONKA (a,b,c,d,eps_it2,1,M)
        eps_it2(M+1) = 0.
     
        call CHECK_MIN_VALUE(eps_it2, M, eps_min)
    
    
        deps_it = maxval(abs(eps_it1-eps_it2))
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
        if (iterat) then
          if ((dE_it > dE_it_max .or. deps_it > deps_it_max) .and. iter_t < maxiter) then
            do i = 1, M+1
              E_it1(i)   = E_it1(i)   * al_it + E_it2(i)   * (1-al_it)
              eps_it1(i) = eps_it1(i) * al_it + eps_it2(i) * (1-al_it)
            enddo
            iter = iter + 1
            iter_t = iter_t + 1
            if (iter == 8) iter = 0
          else
            if (dE_it < 1.E-6 .and. deps_it < 1.E-7) then
              do i = 1, M+1
                E2(i) = E_it2(i)
                eps2(i) = eps_it2(i)
              enddo
              iter = 0
              iter_t = 0
              nl = 2
              cycle_keps = .false.
            else
              do i = 1, M+1
                eps_it1(i) = eps1(i)
                E_it1(i) = E1(i)
              enddo
              iterat = .false.
              iter_t = 0
              nl = 1
            endif
          endif
        else
          do i = 1, M+1
            E2(i) = E_it2(i)
            eps2(i) = eps_it2(i)
          enddo
    
          ! E2   = E_it1   * al_it + E_it2   * (1-al_it)
          ! eps2 = eps_it1 * al_it + eps_it2 * (1-al_it)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          iter = 0
          iter_t = 0
          cycle_keps = .false.
        endif
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
                  
      enddo iterat_keps
    
      
      if (.not.init_keps(ix,iy)) then
    
        ! K-epsilon solver is implemented in initialization mode, i.e.
        ! to obtain initial profiles of E1 and eps1. 
        ! Time derivatives in k-epslilon equations are neglected in this case.
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
        do i = 1, M+1
          E1(i) = E2(i)
          eps1(i) = eps2(i)
        enddo
        init_keps(ix,iy) = .true.
      else
    
        ! The case when k-epsilon solver is implemented in evolution mode, i.e.
        ! E2 and epsilon2 are obtained at the next time step.
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
        next_tstep_keps = .true.
      endif
    
      enddo keps_mode
    
        
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      if (smooth) then
        call VSMOP_LAKE(E2,E2,lbound(E2,1),ubound(E2,1),1,M,vdamp,perdam)
        call CHECK_MIN_VALUE(E2, M, E_min)
        call VSMOP_LAKE(eps2,eps2,lbound(eps2,1),ubound(eps2,1),1,M,vdamp,perdam)
        call CHECK_MIN_VALUE(eps2, M, eps_min)
      endif
    
      do i = 1, M+1
        E12(i) = 0.5d0*(E1(i) + E2(i))
        eps12(i) = 0.5d0*(eps1(i) + eps2(i))
      enddo
    
    
      if (stabfunc%par /= 1) then
        do i = 1, M
            work(i,1) = E2(i)*E2(i)/((eps2(i) + ACC2)*(eps2(i) + ACC2))* & 
    
            &        (rhotemp(i) * & ! Assuming Pr = Sc
    
            &        ( Sal2(i+1) - Sal2(i)) ) / &
            &        (h1*ddz(i)) *AL ! ~ bouyancy production of TKE
            work(i,2) = E2(i)*E2(i)/((eps2(i) + ACC2)*(eps2(i) + ACC2))* &
            & ( (U2(i+1) - U2(i))*(U2(i+1) - U2(i)) + &
            &   (V2(i+1) - V2(i))*(V2(i+1) - V2(i)) ) / &
            &   (h1*h1*ddz(i)*ddz(i)) ! ~ shear production of TKE
            if (stabfunc%par == 2) then
              CEt(i) = CEt_CANUTO(work(i,2), work(i,1))
            elseif (stabfunc%par == 3) then
    
              CEt(i) = sqrt(2.d0)*CL_K_KL_MODEL*SHEAT_GALPERIN(work(i,1))
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      do i = 1, M
    
        KT(i) = max(CEt(i)*E2(i)**2/(eps2(i) + ACCk),min_diff)*cw_m_row0 !E2, eps2
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      enddo