!==========================================================================================!
!==========================================================================================!
! Subroutine nakanishi                                                                     !
! Developed by Marcos Longo (Lab. MASTER - Univ. Sao Paulo/EPS Harvard University)         !
! Sao Paulo, April 29, 2005                                                                !
!                                                                                          !
! This routine is intended to find the TKE tendency by using a 2.5-level model             !
! based on Mellor-Yamada scheme. It's actually an update of the existent M-Y code          !
! in RAMS, outputing also the Sig-W (vertical velocity standard-deviation) and             !
! the Lagrangian time scale, Obukhov Length and PBL depth, using the following references: !
!                                                                                          !
! JANJIC, Z. I. Nonsingular implementation of the Mellor-Yamada level 2.5 scheme in the    !
!   NCEP meso model, office note # 437, National Centers for Environmental Prediction,     !
!   2001, 61 pp.                                                                           !
!                                                                                          !
! NAKANISHI, M. Improvement of the Mellor-Yamada turbulence closure model based on         !
!   large-eddy simulation data. Boundary-Layer Meteor., v. 99, p. 349-378, 2001.           !
!                                                                                          !
! NAKANISHI, M.; NIINO, H. An improved Mellor-Yamada level-3 model with condens-           !
!   ation physics: its design and verification. Boundary-Layer Meteor., v.112,             !
!   p. 1-31, 2004.                                                                         !
!                                                                                          !
! NAKANISHI, M.; NIINO, H. An improved Mellor-Yamada level-3 model with condensation       !
!   physics: its numerical stability and application to a regional prediction of advection !
!   fog. Boundary-Layer Meteor., vol. 119, p. 397-407, 2006.                               !
!                                                                                          !
! HANNA, S. R. Application in air pollution modeling. In: NIEUSWSTADT, F. M. T.;           !
!   VAN DOP, H. Atmospheric turbulence and air pollution modelling. Dordrecht: D.          !
!   Reidel Publishing Company, 1982, chap. 7, p. 275-310.                                  !
!                                                                                          !
! VOGEZELANG, D. H. P.; HOLTSLAG, A. M. Evaluation and model impacts of alternati-         !
!   ve boundary-layer height formulations. Boundary-Layer Meteor., v. 81, p. 245-          !
!   269, 1996.                                                                             !
!==========================================================================================!
!==========================================================================================!
subroutine nakanishi(m1,m2,m3,m4,ia,iz,ja,jz,jd,tkep,tket,vt3dd,vt3de,vt3dh,vt3di          &
                    ,vt3dj,scr1,rtgt,theta,rv,rtp,dn0,up,vp,vegz0,patchz0,tstar,ustar      &
                    ,patch_area,sflux_u,sflux_v,sflux_t,flpu_R,flpv_R,flpw_R,kpbl,pblhgt,lmo     &
                    ,tl,sigw)
   !---------------------------------------------------------------------------------------!
   !   Reference for variables, following Nakanishi (2001) and Nakanishi and Niino         !
   ! (2004,2006), which are based on Helfand and Labraga (1988) and Mellor and Yamada      !
   ! (1974,1982).                                                                          !
   !                                                                                       !
   ! dzloc  -> dz for integrals                                                            !
   ! scr1   -> Km x rho0                                                                     !
   ! sumtkz -> int_0^H z sqrt(e) dz                                                        !
   ! sumtk  -> int_0^H sqrt(e) dz                                                          !
   ! tkep   -> 1/2 q^2 = e (TKE)                                                            !
   ! vctr1  -> z                                                                           !
   ! vctr5  -> zeta
   ! vctr9  -> L                                                                           !
   ! vctr19 -> Gm                                                                          !
   ! vctr20 -> Gh                                                                          !
   ! vctr21 -> THETAV0                                                                     !
   ! vctr22 -> THETA V', used for PBL in the unstable case                                 !
   ! vctr23 -> Sm                                                                          !
   ! vctr24 -> Sh                                                                          !
   ! vctr29 -> zonal wind in thermodynamic grid                                            !
   ! vctr30 -> sqrt(2e) = q                                                                !
   ! vctr31 -> meridional wind in thermodynamic grid                                       !
   ! vctr32 -> q dz                                                                        !
   ! vctr33 -> q z dz                                                                      !
   ! vctr35 -> q^3                                                                          !
   ! vctr36 -> Ps+Pb=(q^3/L) x (Sm x Gm+Sh x Gh)                                                  !
   ! vctr37 -> epsilon=(q^3/B1 x L)                                                           !
   ! vctr38 -> 2e = q^2                                                                     !
   ! vt3dd  -> dU/dz at the u points                                                       !
   ! vt3de  -> dV/dz at the v points                                                       !
   ! vt3dh  -> Kh*rho0                                                                     !
   ! vt3di  -> (dU/dz)^2+(dV/dz)^2 (M^2) at the beginning, then Kq x rho0                       !
   ! vt3dj  -> g/theta_v d(THETA_v)/dz (N^2)                                                   !
   !---------------------------------------------------------------------------------------!

   !---------------------------------------------------------------------------------------!
   ! Loading modules                                                                       !
   !---------------------------------------------------------------------------------------!
   use mem_grid,     only: nstbot,zm,zt
   use mem_scratch,  only:  vctr1, vctr5, vctr9,vctr19,vctr20,vctr21,vctr22,vctr23,vctr24  &
                           ,vctr29,vctr30,vctr31,vctr32,vctr33,vctr38,vctr35,vctr36,vctr37

   use rconstants,     only: abslmomin,abswltlmin,cp,grav,ltscalemax,sigwmin,vonk,lturbmin &
                            ,tkmin,onethird

   use turb_constants, only:  nna1,nna2,nnb1,nnb2,nnc1,nnc2,nnc3,nnc4,nnc5                 &
                             ,nnq1,nnq2,nnq3,nnq4,nnq5,nnq6                                &
                             ,nngama1,nngama2,nnf1,nnf2,nnrf1,nnrf2,nnri1,nnri2,nnri3      &
                             ,nnrfc,nnce1a,nnce1b,nnce2,nnce3,nnce4,nncr1                  &
                             ,nno1,nno2,nno3,nno4, nno5, nno6, nno7, nno8, nnreq,nnrsl     &
                             ,nnmacheps

   use mem_stilt     , only: virtt
   implicit none

   !---------------------------------------------------------------------------------------!
   !   Variable declaration section.                                                       !
   !---------------------------------------------------------------------------------------!
   !----- Arguments (Input/Output/Both) ---------------------------------------------------!
   integer  , intent(in)                          :: m1,m2,m3,m4,ia,iz,ja,jz,jd
   real  , intent(in)    , dimension(m2,m3)    :: flpu_R,flpv_R,flpw_R
   real     , intent(in)    , dimension(m2,m3)    :: rtgt, sflux_u, sflux_v, sflux_t
   real     , intent(in)    , dimension(m1,m2,m3) :: tkep, vt3dd, vt3de, theta
   real     , intent(in)    , dimension(m1,m2,m3) :: vt3dj, rv, rtp, dn0, up, vp
   real     , intent(in)    , dimension(m2,m3,m4) :: vegz0, patchz0, ustar
   real     , intent(in)    , dimension(m2,m3,m4) :: tstar, patch_area 
   real     , intent(inout) , dimension(m2,m3)    :: kpbl
   real     , intent(inout) , dimension(m2,m3)    :: pblhgt, lmo
   real     , intent(inout) , dimension(m1,m2,m3) :: vt3dh, scr1, tl, sigw
   real     , intent(inout) , dimension(m1,m2,m3) :: vt3di,tket
   !----- Internal variables --------------------------------------------------------------!
   integer                                        :: i,j,k,p,k2u,k2um1,k2v,k2vm1,k2w 
   real                                           :: weightsurf,du0dz,dv0dz,sumtkz,sumtk       
   real                                           :: tket2,aux,dzloc,rf,ri,sh2,sm2,gm2
   real                                           :: qlevel2,oneminusalpha,oneminusalpha2
   real                                           :: e1,e2,e3,e4,r1,r2,wltl0,lt,ls,lb
   real                                           :: suminv,z0w,ustarw,tstarw,wstarw,qc
   real                                           :: dzloc1up,dzloc1dn,dzloc2,dq3dz
   real                                           :: dlsmdz,d2q3dz2,janjc,janjd,janjg
   real                                           :: janjh,janji,janjp1,janjt1
   logical                                        :: stable,neutral
   !----- Constants -----------------------------------------------------------------------!
   real, parameter                                :: ustarmin=0.1
   !----- Function ------------------------------------------------------------------------!
   real, external                                 :: cbrt
   !---------------------------------------------------------------------------------------!
   integer  , dimension(m2,m3)    :: flpu,flpv,flpw
   !---------------------------------------------------------------------------------------!
   !   Big loops on horizontal.                                                            !
   !---------------------------------------------------------------------------------------!

   !lfr: Solving a  problem with integer inside vtables
   flpu=int(flpu_R);flpv=int(flpv_R);flpw=int(flpw_R)

   tl = 0.
   
   do j=ja,jz
      do i=ia,iz
         !k2u   = nint(flpu(i,j))
         !k2um1 = nint(flpu(i-1,j))
         !k2v   = nint(flpv(i,j))
         !k2vm1 = nint(flpv(i,j-jd))
         !k2w   = nint(flpw(i,j))

         k2u   = (flpu(i,j))
         k2um1 = (flpu(i-1,j))
         k2v   = (flpv(i,j))
         k2vm1 = (flpv(i,j-jd))
         k2w   = (flpw(i,j))


         !---------------------------------------------------------------------------------!
         !  k2w is the first useful level                                                  !
         !---------------------------------------------------------------------------------!
         do k=k2w,m1-1
            vctr30(k) = max(sqrt(2.0 * tkep(k,i,j)),sqrt(2*tkmin))    ! q 
            vctr38(k) = max(2.0 * tkep(k,i,j),2*tkmin)                ! q^2
            vctr35(k) = vctr38(k)*vctr30(k)                           ! q^3

            !------------------------------------------------------------------------------!
            !  Restrictions in M^2, just to avoid singularities (following Janjic 2001)     !
            !------------------------------------------------------------------------------!
            if (vt3di(k,i,j) < nnmacheps)  vt3di(k,i,j)=(1.+nnmacheps)*nnmacheps*nnreq

            !------------------------------------------------------------------------------!
            !    Determining some integrated variables which will be necessary for length  !
            ! scale derivation.                                                            !
            !------------------------------------------------------------------------------!
            vctr1(k)=(zt(k)-zm(k2w-1))*rtgt(i,j) 
            dzloc=(zm(k)-zm(k-1))*rtgt(i,j)
            vctr33(k)=vctr30(k)*dzloc
            vctr32(k)=vctr33(k)*vctr1(k)
            !----- Deriving some variables that are needed for PBL depth estimation -------!
            vctr21(k)=virtt(theta(k,i,j),rv(k,i,j),rtp(k,i,j))
            vctr29(k)=0.5*(up(k,i,j)+up(k,i-1,j)) 
            vctr31(k)=0.5*(vp(k,i,j)+vp(k,i,j-jd))
         end do !----- k=k2,m1-1 !
    
         !----- The sum is between k2w and m1-1: k1 and m1 are just boundaries ------------!
         sumtkz=sum(vctr32(k2w:(m1-1)))
         sumtk =sum(vctr33(k2w:(m1-1)))
         !---------------------------------------------------------------------------------!

         !----- Finding the average dU/dz and dV/dz at the thermodynamic level k2 ---------!
         du0dz=0.25 * ( vt3dd(k2um1,i-1,j)    + vt3dd(k2u,i,j)                             &
                      + vt3dd(k2um1-1,i-1,j)  + vt3dd(k2u-1,i,j) )
         dv0dz=0.25 * ( vt3de(k2vm1,i,j-jd)   + vt3de(k2v,i,j)                             &
                      + vt3de(k2vm1-1,i,j-jd) + vt3de(k2v-1,i,j) )
         !---------------------------------------------------------------------------------!


         !---------------------------------------------------------------------------------!
         !    Here I'll add a weighted average to z0, ustar and tstar, to avoid patch      !
         ! dependance.                                                                     !
         !---------------------------------------------------------------------------------!
         z0w=0.                                                      
         ustarw=0.
         tstarw=0.
         do p=1,m4
            z0w    =    z0w + max(vegz0(i,j,p),patchz0(i,j,p)) * patch_area(i,j,p)
            ustarw = ustarw + ustar(i,j,p)                     * patch_area(i,j,p)
            tstarw = tstarw + tstar(i,j,p)                     * patch_area(i,j,p)
         end do
         ustarw = max(ustarw,ustarmin)
         !---------------------------------------------------------------------------------!
         ! Here I'll find the <w'Theta'>g and truncate for small values                    !
         !---------------------------------------------------------------------------------!
         if (abs(sflux_t(i,j)/dn0(k2w,i,j)) < abswltlmin) then
            wltl0=sign(abswltlmin,sflux_t(i,j)/dn0(k2w,i,j))
         else
            wltl0=sflux_t(i,j)/dn0(k2w,i,j)
         end if

         !---------------------------------------------------------------------------------!
         ! Finding the Obukhov length (LMO)                                                !
         ! From Nakanishi (2001), Eq. 4:                                                   !
         !              THETA0 u*^3                                                         !
         ! LMO = - --------------------                                                    !
         !           k g <w theta>(g)                                                      !
         !                                                                                 !
         !   Since LMO can be close to zero, but either positive or negative, it must be   !
         ! truncated to a small value, but keeping the same sign                           !
         !---------------------------------------------------------------------------------!
         lmo(i,j)= - vctr21(k2w) * ustarw * ustarw * ustarw / (vonk * grav * wltl0)
         if (abs(lmo(i,j)) < abslmomin) lmo(i,j)=sign(abslmomin,lmo(i,j))

         do k=k2w,m1-1
            !------------------------------------------------------------------------------!
            !    Calculating the length scale L, based on Nakanishi (2001), and GH (nega-  !
            ! tive of dimensionless square of Brunt Vaisala frequency) and GM (dimension-  !
            ! less square of mean shear)                                                   !
            !------------------------------------------------------------------------------!
            vctr5(k) = vctr1(k)/lmo(i,j) ! zeta
            !------------------------------------------------------------------------------!
            ! Finding Ls, following equation A2 of Nakanishi and Niino (2004):             !
            !------------------------------------------------------------------------------!
            if ( vctr5(k) >= 1 ) then
               ls = vonk * vctr1(k) / nnq1                       
            elseif (vctr5(k) >= 0 ) then                       
               ls = vonk * vctr1(k) / (1 + nnq2 * vctr5(k))       
            else !if (vctr5(k) < 0) then                       
               ls = vonk * vctr1(k) * (1 - nnq3 * vctr5(k))**nnq4   
            end if                                            
            !------------------------------------------------------------------------------!
            ! Finding Lt, following equation A3 of Nakanishi and Niino (2004):             !
            !------------------------------------------------------------------------------!
            lt = nnq5 * sumtkz / sumtk
            !------------------------------------------------------------------------------!
            ! Finding Lb, following equation A4 of Nakanishi and Niino (2004):             !
            !------------------------------------------------------------------------------!
            if (vt3dj(k,i,j) > 0 .and. vctr5(k) >= 0) then
               lb = vctr30(k) / max(sqrt(vt3dj(k,i,j)),1.e-10)
            elseif (vt3dj(k,i,j) > 0 .and. vctr5(k) < 0) then
               qc= cbrt((grav/vctr21(k2w))*wltl0*lt)
               lb = (1. + nnq6 *sqrt(qc/max((lt*sqrt(vt3dj(k,i,j))),1.e-20)))*             &
                          vctr30(k)/max(sqrt(vt3dj(k,i,j)),1.e-10)
            else
               lb = 1.e20 !----- Making it very large but without risking overflow. -------!
            end if
            suminv=(1./ls)+(1./lt)+(1./lb) !-----Equation A1, Nakanishi and Niino (2004) --!
            vctr9(k)=1./suminv
            !------------------------------------------------------------------------------!
            !    Restriction on L: from Nakanishi and Niino (2006), based on Janjic(2001). !
            ! This is to guarantee that Cw, defined as <w^2>/q^2 < Rsl (about 0.14)          !
            !------------------------------------------------------------------------------!
            if (vt3dj(k,i,j) > 0) then
               janjc=(nno1*vt3dj(k,i,j)+nno2*vt3di(k,i,j))*vt3dj(k,i,j)
               janjd=(nno3*vt3di(k,i,j)+nno4*vt3dj(k,i,j))
               janjg=(nno5*vt3dj(k,i,j)+nno6*vt3di(k,i,j))*vt3dj(k,i,j)-3.*nnrsl*janjc
               janjh=nno7*vt3di(k,i,j)+nno8*vt3dj(k,i,j)-3.*nnrsl*janjd
               janji=1-3.*nnrsl
               janjt1=(-janjh+sqrt(janjh*janjh-4.*janjg*janji))/(2.*janji)
               if (janjt1 > nnmacheps*nnmacheps) then
                  vctr9(k)=min(vctr9(k),vctr30(k)/sqrt(janjt1))
               end if
            end if
         end do !----- k=k2,m1-1 !

         !---------------------------------------------------------------------------------!
         !   Obtaining an estimation of the PBL height. First, I need to check whether the !
         ! PBL is stable, neutral, or convective. This is obtained by verifying whether    !
         ! z/LMO is > 0, =0 or <0                                                          !
         !---------------------------------------------------------------------------------!
         stable = vctr5(k2w)  >=  1.e-4
         neutral = abs(vctr5(k2w)) < 1.e-4 
         if (stable .or. neutral) then 
            !------------------------------------------------------------------------------!
            !   If the boundary layer is STABLE or NEUTRAL, then we follow Vogelezang and  !
            ! Holtslag (1996) equation (3); (SBL is defined when (w'theta')g < 0). The PBL !
            ! depth is defined as the first level where average Ri > 0.25.                 !
            !------------------------------------------------------------------------------!
            kpbl(i,j) = k2w
            pblhgt(i,j)=vctr1(k2w)
            sboundlay: do k=k2w+1,m1-1
               ri = grav * (vctr21(k)-vctr21(k2w)) * (vctr1(k)-vctr1(k2w))                 &
                  / ( vctr21(k2w) * ( (vctr29(k)-vctr29(k2w)) * (vctr29(k)-vctr29(k2w))    &
                                    + (vctr31(k)-vctr31(k2w)) * (vctr31(k)-vctr31(k2w))    &
                                    + 100.*ustarw*ustarw ) )
               if (ri >= 0.25) then
                  kpbl(i,j) = k
                  pblhgt(i,j) = 0.5 * (vctr1(k-1)+vctr1(k))
                  exit sboundlay
               end if
            end do sboundlay
         else
            !------------------------------------------------------------------------------!
            !    -> Or, if the PBL is convective, then the PBL is defined as the level     !
            !       where the first minimum of w'theta'                                    !
            !------------------------------------------------------------------------------!
            kpbl(i,j) = k2w
            pblhgt(i,j)=vctr1(k2w)
            aux= wltl0 * grav / (cp * vctr21(k2w))
            convmixlay: do k=k2w+1,m1-1
               kpbl(i,j) = k
               pblhgt(i,j)=0.5*(vctr1(k)+vctr1(k-1))
               wstarw=cbrt(aux * pblhgt(i,j) )
               vctr22(k)=vctr21(k2w)+8.5*wltl0/(wstarw*cp)
               ri = grav * (vctr21(k)-vctr22(k)) * (vctr1(k)-vctr1(k2w))                   &
                  / ( vctr22(k) * ((vctr29(k)-vctr29(k2w)) * (vctr29(k)-vctr29(k2w))       &
                    + (vctr31(k)-vctr31(k2w)) * (vctr31(k)-vctr31(k2w))                    &
                    + 100.*ustarw*ustarw) )
               if (ri >= 0.25) exit convmixlay
            end do convmixlay
         end if
         !---------------------------------------------------------------------------------!

         do k=k2w,m1-1

            !------------------------------------------------------------------------------!
            !    Calculating both gradient and flux Richardson numbers (ri and rf, respec- !
            ! tively):                                                                     !
            !------------------------------------------------------------------------------!
            ri=vt3dj(k,i,j)/vt3di(k,i,j)
            rf=nnri1*(ri+nnri2-sqrt(ri*ri-nnri3*ri+nnri2*nnri2))
            !------------------------------------------------------------------------------!
            ! Finding the SH and SM of the Level 2 model (Nakanish, 2001):                 !
            !------------------------------------------------------------------------------!
            if (rf < nnrfc) then
               sh2=3.*nna2*(nngama1+nngama2)*(nnrfc-rf)/(1.-rf)
               sm2=nna1*nnf1*(nnrf1-rf)*sh2/(nna2*nnf2*(nnrf2-rf))
               qlevel2=max(vctr9(k)*sqrt(nnb1*vt3di(k,i,j)*sm2*(1.-rf)),sqrt(2.*tkmin))
            else
               qlevel2=sqrt(2.*tkmin)
            end if
            !------------------------------------------------------------------------------!
            ! Finding Gm and Gh                                                            !
            !------------------------------------------------------------------------------!
            aux=vctr9(k)*vctr9(k)/vctr38(k)
            vctr19(k)=aux*vt3di(k,i,j)       !----- Gm ------------------------------------!
            vctr20(k)=-aux*vt3dj(k,i,j)      !----- Gh ------------------------------------!

            !------------------------------------------------------------------------------!
            !     Case of growing turbulence, it must use the modified Level 2,5 model      !
            ! (Nakanishi and Niino, 2004), otherwise it just uses the original Level 2,5.   !
            !------------------------------------------------------------------------------!
            if (vctr30(k) < qlevel2) then
               oneminusalpha=vctr30(k)/max(qlevel2,1.e-10)
               oneminusalpha2=oneminusalpha*oneminusalpha
            else
               oneminusalpha=1.
               oneminusalpha2=1.
            end if
            e1=1.+oneminusalpha2*(nnce1a*vctr19(k)-nnce1b*vctr20(k))
            e2=-oneminusalpha2*nnce2*vctr20(k)
            e3=oneminusalpha2*nnce3*vctr19(k)
            e4=1.-oneminusalpha2*nnce4*vctr20(k)
            r1=oneminusalpha*nncr1
            r2=oneminusalpha*nna2
            vctr23(k)=(r2*e2-r1*e4)/(e2*e3-e1*e4) !----- Sm -------------------------------!
            vctr24(k)=(r1*e3-r2*e1)/(e2*e3-e1*e4) !----- Sh -------------------------------!

            !------------------------------------------------------------------------------!
            !    Deriving the convective velocity standard-deviation.                      !
            !    The max was inserted because there is no physical limitation which effec- !
            ! tively prevents the variance to be negative (although in preliminary tests   !
            ! the negative values were significantly smaller than the positive ones.       !
            !------------------------------------------------------------------------------!
            sigw(k,i,j) = sqrt(max(vctr38(k)*(onethird-2*nna1*vctr23(k)*vctr19(k)          &
                        + 4.*nna1*(1. - nnc2)*vctr24(k)*vctr20(k)),sigwmin*sigwmin))

         end do
         !---------------------------------------------------------------------------------!

         !----- Finding the surface weight ------------------------------------------------!
         weightsurf=vctr1(k2w)/vctr1(k2w+1)

         !----- Saving the TKE tendency at the bottom for later ---------------------------!
         tket2=tket(k2w,i,j)
         
         do k=k2w,m1-1
            !------------------------------------------------------------------------------!
            ! The vertical Lagrangian timescale is found following Hanna (1982)            !
            !------------------------------------------------------------------------------!
            if (stable) then
            !------------------------------------------------------------------------------!
            ! Stable PBL, so we shall use equation 7.24                                    !
            !------------------------------------------------------------------------------!
               tl(k,i,j)= 0.10 * (pblhgt(i,j)**0.20) * (vctr1(k)**0.80) / sigw(k,i,j)
            elseif (neutral) then
            !------------------------------------------------------------------------------!
            !    This is the formula proposed by Hanna (1982), but without the Coriolis    !
            ! term: later papers, such as Wilson (200, JAM), also disconsider this.        !
            !------------------------------------------------------------------------------!
               tl(k,i,j)=0.5*vctr1(k)/sigw(k,i,j)
            !------------------------------------------------------------------------------!
            !   For unstable PBLs we should follow equation 7.17, considering all possible !
            ! conditions                                                                   !
            !------------------------------------------------------------------------------!
            elseif (vctr1(k) >= 0.1* pblhgt(i,j)) then ! if z/h > 0.1
               ! If PBL is too low, this number can cause underflow in exp. ---------------!
               aux       = max(40.,5.0*vctr1(k)/pblhgt(i,j))
               tl(k,i,j) = 0.15*pblhgt(i,j)*(1-exp(-aux))/sigw(k,i,j)
            elseif ((z0w - vctr1(k)) <= lmo(i,j)) then ! -(z-z0)/Lmo > 1 and z/h < 0.1
               tl(k,i,j)= 0.59*vctr1(k)/sigw(k,i,j)
            else                                       ! -(z-z0)/Lmo < 1 and z/h < 0.1
               tl(k,i,j)= 0.10*vctr1(k)/(sigw(k,i,j)*(0.55+0.38*(vctr1(k)-z0w)/lmo(i,j)))
            end if
            tl(k,i,j)=min(tl(k,i,j),ltscalemax) ! Avoiding exaggerated values...
            !------------------------------------------------------------------------------!

            !------------------------------------------------------------------------------!
            ! Finding the TKE tendency.                                                    !
            !------------------------------------------------------------------------------!
            !----- Finding Ps+Pb-epsilon=(q^3/L) x (SmGm+ShGh-1/B1) --------------------------!
            vctr36(k)=vctr35(k)*(vctr23(k)*vctr19(k)+vctr24(k)*vctr20(k))/vctr9(k)
            vctr37(k)=vctr35(k)/(nnb1*vctr9(k))
            !------------------------------------------------------------------------------!
            !    Finding the tendency term, based on Nakanishi and Niino equation 5. Note  !
            ! that all the 2 factors vanish because here it's the TKE tendency, and        !
            ! TKE=q^2/2.                                                                    !
            !------------------------------------------------------------------------------!
            tket(k,i,j)=tket(k,i,j)+vctr36(k)-vctr37(k)

            !----- Km=L x q x Sm (scr1) -------------------------------------------------------!
            scr1(k,i,j)=vctr9(k)*vctr30(k)*vctr23(k)*dn0(k,i,j)
            !----- Kh=L x q x Sh (vt3dh) ------------------------------------------------------!
            vt3dh(k,i,j)=vctr9(k)*vctr30(k)*vctr24(k)*dn0(k,i,j)
            !-----  Kq=L x q x Sq (vt3di). N&N (2004), proposed Sq=2*Sm instead of Sq=0.2 -----!
            vt3di(k,i,j)=2.*scr1(k,i,j) ! Since Km=L x q x Sm... 
         end do

         !---------------------------------------------------------------------------------!
         !    Different closure for the surface, using surface fluxes. If weightsurf =0    !
         ! it'll ignore.                                                                   !
         !---------------------------------------------------------------------------------!
         if (nstbot == 1 .and. weightsurf > 0) then
            tket(k2w,i,j) = tket2 -vctr37(k2w) + (1.-weightsurf)*vctr36(k2w+1)             &
                          + weightsurf*(-sflux_u(i,j)*du0dz-sflux_v(i,j)*dv0dz             &
                          + grav * sflux_t(i,j) / theta(k2w,i,j)) / dn0(k2w,i,j)
         end if
      end do
   end do
   return
end subroutine nakanishi
!==========================================================================================!
!==========================================================================================!

!------------------------------------------------------------------------------------------!
!------------------------------------------------------------------------------------------!
! Function cbrt                                                                            !
! Developed by Marcos Longo (Lab. MASTER/IAG/USP)                                          !
!                                                                                          !
!   The aim of this function is simply to find the cubic root of all numbers,              !
! including the negative ones.                                                             !
!------------------------------------------------------------------------------------------!
real function cbrt(x)
implicit none
real, intent(in) :: x
if (x > 0) then
  cbrt=x**(1./3.)
else
  cbrt=-((-x)**(1./3.))
end if 
end function cbrt
!------------------------------------------------------------------------------------------!
subroutine nakanishi_light(m1,m2,m3,m4,ia,iz,ja,jz,jd,tkep,tket,vt3dd,vt3de,vt3dh,vt3di    &
                    ,vt3dj,scr1,rtgt,theta,rv,rtp,dn0,up,vp,vegz0,patchz0,tstar,ustar      &
                    ,patch_area,sflux_u,sflux_v,sflux_t,flpu_R,flpv_R,flpw_R, &
		    !kpbl,pblhgt,!srf
		    lmo     )
                    !,tl,sigw)!srf
   !---------------------------------------------------------------------------------------!
   !   Reference for variables, following Nakanishi (2001) and Nakanishi and Niino         !
   ! (2004,2006), which are based on Helfand and Labraga (1988) and Mellor and Yamada      !
   ! (1974,1982).                                                                          !
   !                                                                                       !
   ! dzloc  -> dz for integrals                                                            !
   ! scr1   -> Km x rho0                                                                     !
   ! sumtkz -> int_0^H z sqrt(e) dz                                                        !
   ! sumtk  -> int_0^H sqrt(e) dz                                                          !
   ! tkep   -> 1/2 q^2 = e (TKE)                                                            !
   ! vctr1  -> z                                                                           !
   ! vctr5  -> zeta
   ! vctr9  -> L                                                                           !
   ! vctr19 -> Gm                                                                          !
   ! vctr20 -> Gh                                                                          !
   ! vctr21 -> THETAV0                                                                     !
   ! vctr22 -> THETA V', used for PBL in the unstable case                                 !
   ! vctr23 -> Sm                                                                          !
   ! vctr24 -> Sh                                                                          !
   ! vctr29 -> zonal wind in thermodynamic grid                                            !
   ! vctr30 -> sqrt(2e) = q                                                                !
   ! vctr31 -> meridional wind in thermodynamic grid                                       !
   ! vctr32 -> q dz                                                                        !
   ! vctr33 -> q z dz                                                                      !
   ! vctr35 -> q^3                                                                          !
   ! vctr36 -> Ps+Pb=(q^3/L) x (Sm x Gm+Sh x Gh)                                                  !
   ! vctr37 -> epsilon=(q^3/B1 x L)                                                           !
   ! vctr38 -> 2e = q^2                                                                     !
   ! vt3dd  -> dU/dz at the u points                                                       !
   ! vt3de  -> dV/dz at the v points                                                       !
   ! vt3dh  -> Kh*rho0                                                                     !
   ! vt3di  -> (dU/dz)^2+(dV/dz)^2 (M^2) at the beginning, then Kq x rho0                       !
   ! vt3dj  -> g/theta_v d(THETA_v)/dz (N^2)                                                   !
   !---------------------------------------------------------------------------------------!

   !---------------------------------------------------------------------------------------!
   ! Loading modules                                                                       !
   !---------------------------------------------------------------------------------------!
   use mem_grid,     only: nstbot,zm,zt
   use mem_scratch,  only:  vctr1, vctr5, vctr9,vctr19,vctr20,vctr21,vctr22,vctr23,vctr24  &
                           ,vctr29,vctr30,vctr31,vctr32,vctr33,vctr38,vctr35,vctr36,vctr37

   use rconstants,     only: abslmomin,abswltlmin,cp,grav,ltscalemax,sigwmin,vonk,lturbmin &
                            ,tkmin,onethird

   use turb_constants, only:  nna1,nna2,nnb1,nnb2,nnc1,nnc2,nnc3,nnc4,nnc5                 &
                             ,nnq1,nnq2,nnq3,nnq4,nnq5,nnq6                                &
                             ,nngama1,nngama2,nnf1,nnf2,nnrf1,nnrf2,nnri1,nnri2,nnri3      &
                             ,nnrfc,nnce1a,nnce1b,nnce2,nnce3,nnce4,nncr1                  &
                             ,nno1,nno2,nno3,nno4, nno5, nno6, nno7, nno8, nnreq,nnrsl     &
                             ,nnmacheps

   use mem_stilt     , only: virtt
   implicit none

   !---------------------------------------------------------------------------------------!
   !   Variable declaration section.                                                       !
   !---------------------------------------------------------------------------------------!
   !----- Arguments (Input/Output/Both) ---------------------------------------------------!
   integer  , intent(in)                          :: m1,m2,m3,m4,ia,iz,ja,jz,jd
   real  , intent(in)    , dimension(m2,m3)    :: flpu_R,flpv_R,flpw_R
   real     , intent(in)    , dimension(m2,m3)    :: rtgt, sflux_u, sflux_v, sflux_t
   real     , intent(in)    , dimension(m1,m2,m3) :: tkep, vt3dd, vt3de, theta
   real     , intent(in)    , dimension(m1,m2,m3) :: vt3dj, rv, rtp, dn0, up, vp
   real     , intent(in)    , dimension(m2,m3,m4) :: vegz0, patchz0, ustar
   real     , intent(in)    , dimension(m2,m3,m4) :: tstar, patch_area 
   !integer  , intent(inout) , dimension(m2,m3)    :: kpbl
   !real     , intent(inout) , dimension(m2,m3)    :: pblhgt
   real     , intent(inout) , dimension(m2,m3)    :: lmo
   real     , intent(inout) , dimension(m1,m2,m3) :: vt3dh, scr1 !, tl, sigw
   real     , intent(inout) , dimension(m1,m2,m3) :: vt3di,tket
   !----- Internal variables --------------------------------------------------------------!
   integer                                        :: i,j,k,p,k2u,k2um1,k2v,k2vm1,k2w 
   real                                           :: weightsurf,du0dz,dv0dz,sumtkz,sumtk       
   real                                           :: tket2,aux,dzloc,rf,ri,sh2,sm2,gm2
   real                                           :: qlevel2,oneminusalpha,oneminusalpha2
   real                                           :: e1,e2,e3,e4,r1,r2,wltl0,lt,ls,lb
   real                                           :: suminv,z0w,ustarw,tstarw,wstarw,qc
   real                                           :: dzloc1up,dzloc1dn,dzloc2,dq3dz
   real                                           :: dlsmdz,d2q3dz2,janjc,janjd,janjg
   real                                           :: janjh,janji,janjp1,janjt1
   logical                                        :: stable,neutral
   !----- Constants -----------------------------------------------------------------------!
   real, parameter                                :: ustarmin=0.1
   !----- Function ------------------------------------------------------------------------!
   real, external                                 :: cbrt
   !---------------------------------------------------------------------------------------!

   !---------------------------------------------------------------------------------------!
   !   Big loops on horizontal.                                                            !
   !---------------------------------------------------------------------------------------!
   integer, dimension(m2,m3)    :: flpu,flpv,flpw
   
   !lfr: Solving a  problem with integer inside vtables
   flpu=int(flpu_R);flpv=int(flpv_R);flpw=int(flpw_R)
   
   !tl = 0.
   
   do j=ja,jz
      do i=ia,iz
         !k2u   = nint(flpu(i,j))
         !k2um1 = nint(flpu(i-1,j))
         !k2v   = nint(flpv(i,j))
         !k2vm1 = nint(flpv(i,j-jd))
         !k2w   = nint(flpw(i,j))

         k2u   = (flpu(i,j))
         k2um1 = (flpu(i-1,j))
         k2v   = (flpv(i,j))
         k2vm1 = (flpv(i,j-jd))
         k2w   = (flpw(i,j))


         !---------------------------------------------------------------------------------!
         !  k2w is the first useful level                                                  !
         !---------------------------------------------------------------------------------!
         do k=k2w,m1-1
            vctr30(k) = max(sqrt(2.0 * tkep(k,i,j)),sqrt(2*tkmin))    ! q 
            vctr38(k) = max(2.0 * tkep(k,i,j),2*tkmin)                ! q^2
            vctr35(k) = vctr38(k)*vctr30(k)                           ! q^3

            !------------------------------------------------------------------------------!
            !  Restrictions in M^2, just to avoid singularities (following Janjic 2001)     !
            !------------------------------------------------------------------------------!
            if (vt3di(k,i,j) < nnmacheps)  vt3di(k,i,j)=(1.+nnmacheps)*nnmacheps*nnreq

            !------------------------------------------------------------------------------!
            !    Determining some integrated variables which will be necessary for length  !
            ! scale derivation.                                                            !
            !------------------------------------------------------------------------------!
            vctr1(k)=(zt(k)-zm(k2w-1))*rtgt(i,j) 
            dzloc=(zm(k)-zm(k-1))*rtgt(i,j)
            vctr33(k)=vctr30(k)*dzloc
            vctr32(k)=vctr33(k)*vctr1(k)
            !----- Deriving some variables that are needed for PBL depth estimation -------!
            vctr21(k)=virtt(theta(k,i,j),rv(k,i,j),rtp(k,i,j))
            !vctr29(k)=0.5*(up(k,i,j)+up(k,i-1,j)) 
            !vctr31(k)=0.5*(vp(k,i,j)+vp(k,i,j-jd))
         end do !----- k=k2,m1-1 !
    
         !----- The sum is between k2w and m1-1: k1 and m1 are just boundaries ------------!
         sumtkz=sum(vctr32(k2w:(m1-1)))
         sumtk =sum(vctr33(k2w:(m1-1)))
         !---------------------------------------------------------------------------------!

         !----- Finding the average dU/dz and dV/dz at the thermodynamic level k2 ---------!
         du0dz=0.25 * ( vt3dd(k2um1,i-1,j)    + vt3dd(k2u,i,j)                             &
                      + vt3dd(k2um1-1,i-1,j)  + vt3dd(k2u-1,i,j) )
         dv0dz=0.25 * ( vt3de(k2vm1,i,j-jd)   + vt3de(k2v,i,j)                             &
                      + vt3de(k2vm1-1,i,j-jd) + vt3de(k2v-1,i,j) )
         !---------------------------------------------------------------------------------!


         !---------------------------------------------------------------------------------!
         !    Here I'll add a weighted average to z0, ustar and tstar, to avoid patch      !
         ! dependance.                                                                     !
         !---------------------------------------------------------------------------------!
         z0w=0.                                                      
         ustarw=0.
         tstarw=0.
         do p=1,m4
            z0w    =    z0w + max(vegz0(i,j,p),patchz0(i,j,p)) * patch_area(i,j,p)
            ustarw = ustarw + ustar(i,j,p)                     * patch_area(i,j,p)
            tstarw = tstarw + tstar(i,j,p)                     * patch_area(i,j,p)
         end do
         ustarw = max(ustarw,ustarmin)
         !---------------------------------------------------------------------------------!
         ! Here I'll find the <w'Theta'>g and truncate for small values                    !
         !---------------------------------------------------------------------------------!
         if (abs(sflux_t(i,j)/dn0(k2w,i,j)) < abswltlmin) then
            wltl0=sign(abswltlmin,sflux_t(i,j)/dn0(k2w,i,j))
         else
            wltl0=sflux_t(i,j)/dn0(k2w,i,j)
         end if

         !---------------------------------------------------------------------------------!
         ! Finding the Obukhov length (LMO)                                                !
         ! From Nakanishi (2001), Eq. 4:                                                   !
         !              THETA0 u*^3                                                         !
         ! LMO = - --------------------                                                    !
         !           k g <w theta>(g)                                                      !
         !                                                                                 !
         !   Since LMO can be close to zero, but either positive or negative, it must be   !
         ! truncated to a small value, but keeping the same sign                           !
         !---------------------------------------------------------------------------------!
         lmo(i,j)= - vctr21(k2w) * ustarw * ustarw * ustarw / (vonk * grav * wltl0)
         if (abs(lmo(i,j)) < abslmomin) lmo(i,j)=sign(abslmomin,lmo(i,j))

         do k=k2w,m1-1
            !------------------------------------------------------------------------------!
            !    Calculating the length scale L, based on Nakanishi (2001), and GH (nega-  !
            ! tive of dimensionless square of Brunt Vaisala frequency) and GM (dimension-  !
            ! less square of mean shear)                                                   !
            !------------------------------------------------------------------------------!
            vctr5(k) = vctr1(k)/lmo(i,j) ! zeta
            !------------------------------------------------------------------------------!
            ! Finding Ls, following equation A2 of Nakanishi and Niino (2004):             !
            !------------------------------------------------------------------------------!
            if ( vctr5(k) >= 1 ) then
               ls = vonk * vctr1(k) / nnq1                       
            elseif (vctr5(k) >= 0 ) then                       
               ls = vonk * vctr1(k) / (1 + nnq2 * vctr5(k))       
            else !if (vctr5(k) < 0) then                       
               ls = vonk * vctr1(k) * (1 - nnq3 * vctr5(k))**nnq4   
            end if                                            
            !------------------------------------------------------------------------------!
            ! Finding Lt, following equation A3 of Nakanishi and Niino (2004):             !
            !------------------------------------------------------------------------------!
            lt = nnq5 * sumtkz / sumtk
            !------------------------------------------------------------------------------!
            ! Finding Lb, following equation A4 of Nakanishi and Niino (2004):             !
            !------------------------------------------------------------------------------!
            if (vt3dj(k,i,j) > 0 .and. vctr5(k) >= 0) then
               lb = vctr30(k) / max(sqrt(vt3dj(k,i,j)),1.e-10)
            elseif (vt3dj(k,i,j) > 0 .and. vctr5(k) < 0) then
               qc= cbrt((grav/vctr21(k2w))*wltl0*lt)
               lb = (1. + nnq6 *sqrt(qc/max((lt*sqrt(vt3dj(k,i,j))),1.e-20)))*             &
                          vctr30(k)/max(sqrt(vt3dj(k,i,j)),1.e-10)
            else
               lb = 1.e20 !----- Making it very large but without risking overflow. -------!
            end if
            suminv=(1./ls)+(1./lt)+(1./lb) !-----Equation A1, Nakanishi and Niino (2004) --!
            vctr9(k)=1./suminv
            !------------------------------------------------------------------------------!
            !    Restriction on L: from Nakanishi and Niino (2006), based on Janjic(2001). !
            ! This is to guarantee that Cw, defined as <w^2>/q^2 < Rsl (about 0.14)          !
            !------------------------------------------------------------------------------!
            if (vt3dj(k,i,j) > 0) then
               janjc=(nno1*vt3dj(k,i,j)+nno2*vt3di(k,i,j))*vt3dj(k,i,j)
               janjd=(nno3*vt3di(k,i,j)+nno4*vt3dj(k,i,j))
               janjg=(nno5*vt3dj(k,i,j)+nno6*vt3di(k,i,j))*vt3dj(k,i,j)-3.*nnrsl*janjc
               janjh=nno7*vt3di(k,i,j)+nno8*vt3dj(k,i,j)-3.*nnrsl*janjd
               janji=1-3.*nnrsl
               janjt1=(-janjh+sqrt(janjh*janjh-4.*janjg*janji))/(2.*janji)
               if (janjt1 > nnmacheps*nnmacheps) then
                  vctr9(k)=min(vctr9(k),vctr30(k)/sqrt(janjt1))
               end if
            end if
         end do !----- k=k2,m1-1 !

         !---------------------------------------------------------------------------------!
         !   Obtaining an estimation of the PBL height. First, I need to check whether the !
         ! PBL is stable, neutral, or convective. This is obtained by verifying whether    !
         ! z/LMO is > 0, =0 or <0                                                          !
         !---------------------------------------------------------------------------------!
         !stable = vctr5(k2w)  >=  1.e-4
         !neutral = abs(vctr5(k2w)) < 1.e-4 
         !if (stable .or. neutral) then 
         !   !------------------------------------------------------------------------------!
         !   !   If the boundary layer is STABLE or NEUTRAL, then we follow Vogelezang and  !
         !   ! Holtslag (1996) equation (3); (SBL is defined when (w'theta')g < 0). The PBL !
         !   ! depth is defined as the first level where average Ri > 0.25.                 !
         !   !------------------------------------------------------------------------------!
         !   kpbl(i,j) = k2w
         !   pblhgt(i,j)=vctr1(k2w)
         !   sboundlay: do k=k2w+1,m1-1
         !      ri = grav * (vctr21(k)-vctr21(k2w)) * (vctr1(k)-vctr1(k2w))                 &
         !         / ( vctr21(k2w) * ( (vctr29(k)-vctr29(k2w)) * (vctr29(k)-vctr29(k2w))    &
         !                           + (vctr31(k)-vctr31(k2w)) * (vctr31(k)-vctr31(k2w))    &
         !                           + 100.*ustarw*ustarw ) )
         !      if (ri >= 0.25) then
         !         kpbl(i,j) = k
         !         pblhgt(i,j) = 0.5 * (vctr1(k-1)+vctr1(k))
         !         exit sboundlay
         !      end if
         !   end do sboundlay
         !else
         !   !------------------------------------------------------------------------------!
         !   !    -> Or, if the PBL is convective, then the PBL is defined as the level     !
         !   !       where the first minimum of w'theta'                                    !
         !   !------------------------------------------------------------------------------!
         !   kpbl(i,j) = k2w
         !   pblhgt(i,j)=vctr1(k2w)
         !   aux= wltl0 * grav / (cp * vctr21(k2w))
         !   convmixlay: do k=k2w+1,m1-1
         !      kpbl(i,j) = k
         !      pblhgt(i,j)=0.5*(vctr1(k)+vctr1(k-1))
         !      wstarw=cbrt(aux * pblhgt(i,j) )
         !      vctr22(k)=vctr21(k2w)+8.5*wltl0/(wstarw*cp)
         !      ri = grav * (vctr21(k)-vctr22(k)) * (vctr1(k)-vctr1(k2w))                   &
         !         / ( vctr22(k) * ((vctr29(k)-vctr29(k2w)) * (vctr29(k)-vctr29(k2w))       &
         !           + (vctr31(k)-vctr31(k2w)) * (vctr31(k)-vctr31(k2w))                    &
         !           + 100.*ustarw*ustarw) )
         !      if (ri >= 0.25) exit convmixlay
         !   end do convmixlay
         !end if
         !---------------------------------------------------------------------------------!

         do k=k2w,m1-1

            !------------------------------------------------------------------------------!
            !    Calculating both gradient and flux Richardson numbers (ri and rf, respec- !
            ! tively):                                                                     !
            !------------------------------------------------------------------------------!
            ri=vt3dj(k,i,j)/vt3di(k,i,j)
            rf=nnri1*(ri+nnri2-sqrt(ri*ri-nnri3*ri+nnri2*nnri2))
            !------------------------------------------------------------------------------!
            ! Finding the SH and SM of the Level 2 model (Nakanish, 2001):                 !
            !------------------------------------------------------------------------------!
            if (rf < nnrfc) then
               sh2=3.*nna2*(nngama1+nngama2)*(nnrfc-rf)/(1.-rf)
               sm2=nna1*nnf1*(nnrf1-rf)*sh2/(nna2*nnf2*(nnrf2-rf))
               qlevel2=max(vctr9(k)*sqrt(nnb1*vt3di(k,i,j)*sm2*(1.-rf)),sqrt(2.*tkmin))
            else
               qlevel2=sqrt(2.*tkmin)
            end if
            !------------------------------------------------------------------------------!
            ! Finding Gm and Gh                                                            !
            !------------------------------------------------------------------------------!
            aux=vctr9(k)*vctr9(k)/vctr38(k)
            vctr19(k)=aux*vt3di(k,i,j)       !----- Gm ------------------------------------!
            vctr20(k)=-aux*vt3dj(k,i,j)      !----- Gh ------------------------------------!

            !------------------------------------------------------------------------------!
            !     Case of growing turbulence, it must use the modified Level 2,5 model      !
            ! (Nakanishi and Niino, 2004), otherwise it just uses the original Level 2,5.   !
            !------------------------------------------------------------------------------!
            if (vctr30(k) < qlevel2) then
               oneminusalpha=vctr30(k)/max(qlevel2,1.e-10)
               oneminusalpha2=oneminusalpha*oneminusalpha
            else
               oneminusalpha=1.
               oneminusalpha2=1.
            end if
            e1=1.+oneminusalpha2*(nnce1a*vctr19(k)-nnce1b*vctr20(k))
            e2=-oneminusalpha2*nnce2*vctr20(k)
            e3=oneminusalpha2*nnce3*vctr19(k)
            e4=1.-oneminusalpha2*nnce4*vctr20(k)
            r1=oneminusalpha*nncr1
            r2=oneminusalpha*nna2
            vctr23(k)=(r2*e2-r1*e4)/(e2*e3-e1*e4) !----- Sm -------------------------------!
            vctr24(k)=(r1*e3-r2*e1)/(e2*e3-e1*e4) !----- Sh -------------------------------!

            !------------------------------------------------------------------------------!
            !    Deriving the convective velocity standard-deviation.                      !
            !    The max was inserted because there is no physical limitation which effec- !
            ! tively prevents the variance to be negative (although in preliminary tests   !
            ! the negative values were significantly smaller than the positive ones.       !
            !------------------------------------------------------------------------------!
            !sigw(k,i,j) = sqrt(max(vctr38(k)*(onethird-2*nna1*vctr23(k)*vctr19(k)          &
            !            + 4.*nna1*(1. - nnc2)*vctr24(k)*vctr20(k)),sigwmin*sigwmin))

         end do
         !---------------------------------------------------------------------------------!

         !----- Finding the surface weight ------------------------------------------------!
         weightsurf=vctr1(k2w)/vctr1(k2w+1)

         !----- Saving the TKE tendency at the bottom for later ---------------------------!
         tket2=tket(k2w,i,j)
         
         do k=k2w,m1-1
            !------------------------------------------------------------------------------!
            ! The vertical Lagrangian timescale is found following Hanna (1982)            !
            !------------------------------------------------------------------------------!
            !if (stable) then
            !------------------------------------------------------------------------------!
            ! Stable PBL, so we shall use equation 7.24                                    !
            !------------------------------------------------------------------------------!
            !   tl(k,i,j)= 0.10 * (pblhgt(i,j)**0.20) * (vctr1(k)**0.80) / sigw(k,i,j)
            !elseif (neutral) then
            !------------------------------------------------------------------------------!
            !    This is the formula proposed by Hanna (1982), but without the Coriolis    !
            ! term: later papers, such as Wilson (200, JAM), also disconsider this.        !
            !------------------------------------------------------------------------------!
            !   tl(k,i,j)=0.5*vctr1(k)/sigw(k,i,j)
            !------------------------------------------------------------------------------!
            !   For unstable PBLs we should follow equation 7.17, considering all possible !
            ! conditions                                                                   !
            !------------------------------------------------------------------------------!
            !elseif (vctr1(k) >= 0.1* pblhgt(i,j)) then ! if z/h > 0.1
            !   ! If PBL is too low, this number can cause underflow in exp. ---------------!
            !   aux       = max(40.,5.0*vctr1(k)/pblhgt(i,j))
            !   tl(k,i,j) = 0.15*pblhgt(i,j)*(1-exp(-aux))/sigw(k,i,j)
            !elseif ((z0w - vctr1(k)) <= lmo(i,j)) then ! -(z-z0)/Lmo > 1 and z/h < 0.1
            !   tl(k,i,j)= 0.59*vctr1(k)/sigw(k,i,j)
            !else                                       ! -(z-z0)/Lmo < 1 and z/h < 0.1
            !   tl(k,i,j)= 0.10*vctr1(k)/(sigw(k,i,j)*(0.55+0.38*(vctr1(k)-z0w)/lmo(i,j)))
            !end if
            !tl(k,i,j)=min(tl(k,i,j),ltscalemax) ! Avoiding exaggerated values...
            !------------------------------------------------------------------------------!

            !------------------------------------------------------------------------------!
            ! Finding the TKE tendency.                                                    !
            !------------------------------------------------------------------------------!
            !----- Finding Ps+Pb-epsilon=(q^3/L) x (SmGm+ShGh-1/B1) --------------------------!
            vctr36(k)=vctr35(k)*(vctr23(k)*vctr19(k)+vctr24(k)*vctr20(k))/vctr9(k)
            vctr37(k)=vctr35(k)/(nnb1*vctr9(k))
            !------------------------------------------------------------------------------!
            !    Finding the tendency term, based on Nakanishi and Niino equation 5. Note  !
            ! that all the 2 factors vanish because here it's the TKE tendency, and        !
            ! TKE=q^2/2.                                                                    !
            !------------------------------------------------------------------------------!
            tket(k,i,j)=tket(k,i,j)+vctr36(k)-vctr37(k)

            !----- Km=L x q x Sm (scr1) -------------------------------------------------------!
            scr1(k,i,j)=vctr9(k)*vctr30(k)*vctr23(k)*dn0(k,i,j)
            !----- Kh=L x q x Sh (vt3dh) ------------------------------------------------------!
            vt3dh(k,i,j)=vctr9(k)*vctr30(k)*vctr24(k)*dn0(k,i,j)
            !-----  Kq=L x q x Sq (vt3di). N&N (2004), proposed Sq=2*Sm instead of Sq=0.2 -----!
            vt3di(k,i,j)=2.*scr1(k,i,j) ! Since Km=L x q x Sm... 
         end do

         !---------------------------------------------------------------------------------!
         !    Different closure for the surface, using surface fluxes. If weightsurf =0    !
         ! it'll ignore.                                                                   !
         !---------------------------------------------------------------------------------!
         if (nstbot == 1 .and. weightsurf > 0) then
	    tket(k2w,i,j) = tket2 -vctr37(k2w) + (1.-weightsurf)*vctr36(k2w+1)             &
                          + weightsurf*(-sflux_u(i,j)*du0dz-sflux_v(i,j)*dv0dz             &
                          + grav * sflux_t(i,j) / theta(k2w,i,j)) / dn0(k2w,i,j)
         end if
      end do
   end do
   return
end subroutine nakanishi_light
!==========================================================================================!
!==========================================================================================!
