1.Different setting for underice and open-water minimal K and epsilon...

1.Different setting for underice and open-water minimal K and epsilon 2.Velocity amplitude changed in constant-density layers identification algorithm

source/model/lake.f90

@@ -2495,6 +2495,7 @@ if (dyn_pgrad%par == 5) then !diagnostic computation of TKE dissipation in a riv
   & g*tan(pi*alphax/180.),z_half,eps_ziriv)
 endif
 
+
 !print*,h1,l1,hs1,ls1,Tw2(1),Sal2(1),Sal2(M+1)
 !if (ix == 10 .and. iy == 10) print*, 'Lake', &
 !& 'T1=',Tw1(1),'T2=',Tw1(2),'h1=',h1,'S=',shortwave, &
@@ -2786,6 +2787,7 @@ if (flag_print) then ! The output in ASCII files
     &  eps1, z_half, M, &
     &  zgrid_out, ngrid_out%par, & !ngrid_out
     &  outpath, 'epsilon', i, ndec, .false.)
+    if (dyn_pgrad%par == 5) then
       ndec = -1
       i = i + 2
       call PROFILE_OUTPUT &
@@ -2795,6 +2797,7 @@ if (flag_print) then ! The output in ASCII files
       &  eps_ziriv, z_half, M, &
       &  zgrid_out, ngrid_out%par, & !ngrid_out
       &  outpath, 'epsziriv', i, ndec, .false.)
+    endif
     ndec = -1
     i = i + 2
     call PROFILE_OUTPUT &

source/model/momentum_mod.f90

@@ -1846,7 +1846,7 @@ real(kind=ireals), intent(out) :: LxLyCDL(1:M+1,1:2) !Horizontal dimensions of c
 
 
 !Local variables
-real(kind=ireals), parameter :: a = 1.e-1, u = 0.1
+real(kind=ireals), parameter :: a = 1.e-1, u = 0.01
 real(kind=ireals) :: omega, dro, dro_, lnumx, lnumy,ldenx,ldeny, maxN, Lxi, Lyi
 real(kind=ireals), allocatable :: Hlars(:), work(:), depths(:)
 !integer(kind=iintegers) :: iwork(1:M+2)

source/model/turb_mod.f90

@@ -321,6 +321,7 @@ firstcall_if : if (firstcall) then
 
 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
@@ -328,6 +329,7 @@ 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
+endif
 
 dt05 = 0.5d0*dt