  !*****************************************************************************
  !                                                                            *
  ! 2021 L. Bakels: This module contains all turbulence related subroutines    *
  ! 2023 PS: include psih, psim into this module
  !                                                                            *
  !*****************************************************************************

module turbulence_mod
  use par_mod
  use com_mod
  use particle_mod
  use pbl_profile_mod

  implicit none

  real :: ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw
  real :: sigw,dsigwdz,dsigw2dz

!$OMP THREADPRIVATE(ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw, &
!$OMP sigw,dsigwdz,dsigw2dz)

contains

subroutine turbulence_pbl(ipart,nrand,dt,zts,rhoa,rhograd,thread)

  use cbl_mod

  implicit none

  integer, intent(in) ::   &
    ipart,                 & ! particle index
    thread                   ! number of the omp thread
  integer, intent(inout) ::&
    nrand                    ! random number used for turbulence
  real,intent(in) ::       &
    dt,                    & ! real(ldt)
    rhoa,                  & ! air density, used in CBL
    rhograd                  ! vertical gradient of the air density, used in CBL
  real,intent(inout) ::    &
    zts                      ! local 'real' copy of the particle position
  real ::                     &
    delz,                     & ! change in vertical position due to turbulence
    ru,rv,rw,wp,icbt_r,       & ! used for computing turbulence
    dtf,rhoaux,dtftlw,ath,bth,& ! CBL related
    ptot_lhh,Q_lhh,phi_lhh,   & ! CBL related
    old_wp_buf                  ! CBL related
  integer ::                  &
    flagrein,                 & ! flag used in CBL scheme
    i                           ! loop variable
  integer(kind=2) :: icbt

  ! tlw,dsigwdz and dsigw2dz are defined in hanna
    if (turbswitch) then
      call hanna(zts)
    else
      call hanna1(zts)
    endif

  !*****************************************
  ! Determine the new diffusivity velocities
  !*****************************************

  ! Horizontal components
  !**********************
    if (nrand+1.gt.maxrand) nrand=1
    if (dt/tlu.lt..5) then
      part(ipart)%turbvel%u=(1.-dt/tlu)*part(ipart)%turbvel%u+rannumb(nrand)*sigu*sqrt(2.*dt/tlu)
    else
      ru=exp(-dt/tlu)
      part(ipart)%turbvel%u=ru*part(ipart)%turbvel%u+rannumb(nrand)*sigu*sqrt(1.-ru**2)
    endif
    if (dt/tlv.lt..5) then
      part(ipart)%turbvel%v=(1.-dt/tlv)*part(ipart)%turbvel%v+rannumb(nrand+1)*sigv*sqrt(2.*dt/tlv)
    else
      rv=exp(-dt/tlv)
      part(ipart)%turbvel%v=rv*part(ipart)%turbvel%v+rannumb(nrand+1)*sigv*sqrt(1.-rv**2)
    endif
    nrand=nrand+2

    if (nrand+ifine.gt.maxrand) nrand=1
    rhoaux=rhograd/rhoa
    dtf=dt*fine

    dtftlw=dtf/tlw

  ! Loop over ifine short time steps for vertical component
  !********************************************************
    wp=part(ipart)%turbvel%w
    icbt=part(ipart)%icbt
    do i=1,ifine
      icbt_r=real(icbt)
  ! Determine the drift velocity and density correction velocity
  !*************************************************************

      if (turbswitch) then
        if (dtftlw.lt..5) then
  !*************************************************************
  !************** CBL options added by mc see routine cblf90 ***
          ! LB needs to be checked if this works with openmp
          if (cblflag.eq.1) then  !modified by mc
            if (-h/ol.gt.5) then  !modified by mc
              flagrein=0
              nrand=nrand+1
              old_wp_buf=wp
              call cbl(wp,zts,wst,h,rhoa,rhograd,&
                sigw,dsigwdz,tlw,ptot_lhh,Q_lhh,phi_lhh,ath,bth,ol,flagrein) !inside the routine for inverse time
              wp=(wp+ath*dtf+&
                bth*rannumb(nrand)*sqrt(dtf))*icbt_r
              delz=wp*dtf
              if ((flagrein.eq.1).or.(wp.ne.wp).or.((wp-1.).eq.wp)) then
                call reinit_particle(zts,wst,h,sigw,old_wp_buf,nrand,ol)
                wp=old_wp_buf
                delz=wp*dtf
                nan_count(thread+1)=nan_count(thread+1)+1
              end if
            else
              nrand=nrand+1
              old_wp_buf=wp
              ath=-wp/tlw+sigw*dsigwdz+&
                wp*wp/sigw*dsigwdz+sigw*sigw/rhoa*rhograd  !1-note for inverse time should be -wp/tlw*ldirect+... calculated for wp=-wp
                                                                                  !2-but since ldirect =-1 for inverse time and this must be calculated for (-wp) and
                                                                                  !3-the gaussian pdf is symmetric (i.e. pdf(w)=pdf(-w) ldirect can be discarded
              bth=sigw*rannumb(nrand)*sqrt(2.*dtftlw)
              wp=(wp+ath*dtf+bth)*icbt_r
              delz=wp*dtf
              if ((wp.ne.wp).or.((wp-1.).eq.wp)) then ! Catch infinity or NaN
                nrand=nrand+1
                wp=sigw*rannumb(nrand)
                delz=wp*dtf
                nan_count(thread+1)=nan_count(thread+1)+1
              end if
            end if
  !******************** END CBL option *******************************
  !*******************************************************************
          else
               wp=((1.-dtftlw)*wp+rannumb(nrand+i)*sqrt(2.*dtftlw) &
               +dtf*(dsigwdz+rhoaux*sigw))*icbt_r
               delz=wp*sigw*dtf
          end if
        else
          rw=exp(-dtftlw)
          wp=(rw*wp+rannumb(nrand+i)*sqrt(1.-rw**2) &
               +tlw*(1.-rw)*(dsigwdz+rhoaux*sigw))*icbt_r
          delz=wp*sigw*dtf
        endif

      else
        rw=exp(-dtftlw)
        wp=(rw*wp+rannumb(nrand+i)*sqrt(1.-rw**2)*sigw &
             +tlw*(1.-rw)*(dsigw2dz+rhoaux*sigw**2))*icbt_r
        delz=wp*dtf
      endif

  !****************************************************
  ! Compute turbulent vertical displacement of particle
  !****************************************************

      if (abs(delz).gt.h) delz=mod(delz,h)

  ! Determine if particle transfers to a "forbidden state" below the ground
  ! or above the mixing height
  !************************************************************************

      if (delz.lt.-zts) then         ! reflection at ground
        icbt=-1
        call set_z(ipart,-zts-delz)
      else if (delz.gt.(h-zts)) then ! reflection at h
        icbt=-1
        call set_z(ipart,-zts-delz+2.*h)
      else                         ! no reflection
        icbt=1
        call set_z(ipart,zts+delz)
      endif

      if (i.ne.ifine) then
        zeta=zts/h
        call hanna_short(zts)
      endif
      zts=real(part(ipart)%z)
    end do
    part(ipart)%turbvel%w=wp
    part(ipart)%icbt=icbt
    if (cblflag.ne.1) nrand=nrand+i
end subroutine turbulence_pbl

subroutine turbulence_above_pbl(dt,nrand,ux,vy,wp,tropop,zts)

  implicit none

  integer, intent(inout) ::       &
    nrand                           ! random number used for turbulence
  real, intent(inout) ::          &
    ux,vy,wp                        ! random turbulent velocities above PBL
  real, intent(in) ::             &
    tropop,                       & ! height of troposphere
    zts,                          & ! height of particle
    dt                              ! real(ldt)
  real ::                         &
    uxscale,wpscale,              & ! factor used in calculating turbulent perturbations above PBL
    weight                          ! transition above the tropopause

  if (zts.lt.tropop) then  ! in the troposphere
    uxscale=sqrt(2.*d_trop/dt)
    if (nrand+1.gt.maxrand) nrand=1
    ux=rannumb(nrand)*uxscale
    vy=rannumb(nrand+1)*uxscale
    nrand=nrand+2
    wp=0.
  else if (zts.lt.tropop+1000.) then     ! just above the tropopause: make transition
    weight=(zts-tropop)/1000.
    uxscale=sqrt(2.*d_trop/dt*(1.-weight))
    if (nrand+2.gt.maxrand) nrand=1
    ux=rannumb(nrand)*uxscale
    vy=rannumb(nrand+1)*uxscale
    wpscale=sqrt(2.*d_strat/dt*weight)
    wp=rannumb(nrand+2)*wpscale+d_strat/1000.
    nrand=nrand+3
  else                 ! in the stratosphere
    if (nrand.gt.maxrand) nrand=1
    ux=0.
    vy=0.
    wpscale=sqrt(2.*d_strat/dt)
    wp=rannumb(nrand)*wpscale
    nrand=nrand+1
  endif
end subroutine turbulence_above_pbl

subroutine turbulence_mesoscale(nrand,dxsave,dysave,ipart,usig,vsig,wsig,wsigeta,eps_eta)

  implicit none

  integer, intent(inout) ::       &
    nrand                           ! random number used for turbulence
  integer, intent(in) ::          &
    ipart                              ! particle index
  real, intent(in) ::             &
    eps_eta,usig,vsig,wsig,wsigeta
  real, intent(inout) ::          &
    dxsave,dysave                   ! accumulated displacement in long and lat
  real ::                         &
    r,rs                            ! mesoscale related

  r=exp(-2.*real(abs(lsynctime))/real(lwindinterv))
  rs=sqrt(1.-r**2)
  if (nrand+2.gt.maxrand) nrand=1
  part(ipart)%mesovel%u=r*part(ipart)%mesovel%u+rs*rannumb(nrand)*usig*fturbmeso
  part(ipart)%mesovel%v=r*part(ipart)%mesovel%v+rs*rannumb(nrand+1)*vsig*fturbmeso
  dxsave=dxsave+part(ipart)%mesovel%u*real(lsynctime)
  dysave=dysave+part(ipart)%mesovel%v*real(lsynctime)

#ifdef ETA
  part(ipart)%mesovel%w=r*part(ipart)%mesovel%w+rs*rannumb(nrand+2)*wsigeta*fturbmeso
  call update_zeta(ipart,part(ipart)%mesovel%w*real(lsynctime))
  if (part(ipart)%zeta.ge.1.) call set_zeta(ipart,1.-(part(ipart)%zeta-1.))
  if (part(ipart)%zeta.eq.1.) call update_zeta(ipart,-eps_eta)

#else
  part(ipart)%mesovel%w=r*part(ipart)%mesovel%w+rs*rannumb(nrand+2)*wsig*fturbmeso
  call update_z(ipart,part(ipart)%mesovel%w*real(lsynctime))
  if (part(ipart)%z.lt.0.) call set_z(ipart,-1.*part(ipart)%z)    ! if particle below ground -> refletion
#endif
end subroutine turbulence_mesoscale

subroutine hanna(z)
  !                 i
  !*****************************************************************************
  !                                                                            *
  !   Computation of \sigma_i and \tau_L based on the scheme of Hanna (1982)   *
  !   Source: 'Atmospheric Turbulence and Air Polution', chapter 4 and 7,      *
  !   J.A. Businger, edited by F.T.M. Nieuwstadt and H. van Dop                *            *
  !                                                                            *
  !   Author: A. Stohl                                                         *
  !                                                                            *
  !   4 December 1997                                                          *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  ! dsigwdz [1/s]     vertical gradient of sigw                                *
  ! ol [m]            Obukhov length                                           *
  ! sigu, sigv, sigw  standard deviations of turbulent velocity fluctuations   *
  ! tlu [s]           Lagrangian time scale for the along wind component.      *
  ! tlv [s]           Lagrangian time scale for the cross wind component.      *
  ! tlw [s]           Lagrangian time scale for the vertical wind component.   *
  ! ust, ustar [m/s]  friction velocity                                        *
  ! wst, wstar [m/s]  convective velocity scale                                *
  !                                                                            *
  !*****************************************************************************

  implicit none

  real :: corr,z

  !**********************
  ! 1. Neutral conditions
  !**********************

  ! The addition of 1.e-2 in sigu,sigv,sigw comes from ???

  if (h/abs(ol).lt.1.) then
    ust=max(1.e-4,ust)
    corr=z/ust

    ! Eq. 7.25 Hanna 1982: sigu/ust=2.0*exp(-3*f*z/ust),
    ! where f, the Coriolis parameter, is set to 1e-4
    sigu=1.e-2+2.0*ust*exp(-3.e-4*corr)

    ! Eq. 7.26 Hanna 1982: sigv/ust=sigw/ust=1.3*exp(-2*f*z/ust),
    ! where f, the Coriolis parameter, is set to 1e-4
    sigw=1.3*ust*exp(-2.e-4*corr)

    ! ???
    dsigwdz=-2.e-4*sigw
    sigw=sigw+1.e-2
    sigv=sigw

    ! Eq.7.27 Hanna 1982: TL=0.5*z/sigw/(1+15*f*z/ust) assumed to be valid
    ! for all three components
    tlu=0.5*z/sigw/(1.+1.5e-3*corr)
    tlv=tlu
    tlw=tlu

  !***********************
  ! 2. Unstable conditions
  !***********************

  else if (ol.lt.0.) then

  ! Determine sigmas
  !*****************

    ! Eq. 4.15 Caughey 1982
    sigu=1.e-2+ust*(12.-0.5*h/ol)**0.33333
    sigv=sigu

    ! Ryall & Maryon 1998
    sigw=sqrt(1.2*wst**2*(1.-.9*zeta)*zeta**0.66666+ &
         (1.8-1.4*zeta)*ust**2)+1.e-2
    ! ???
    dsigwdz=0.5/sigw/h*(-1.4*ust**2+wst**2* &
         (0.8*max(zeta,1.e-3)**(-.33333)-1.8*zeta**0.66666))

  ! Determine average Lagrangian time scale
  !****************************************

    ! Eq. 7.17 Hanna  1982
    tlu=0.15*h/sigu
    tlv=tlu
    if (z.lt.abs(ol)) then
      tlw=0.1*z/(sigw*(0.55-0.38*abs(z/ol)))
    else if (zeta.lt.0.1) then
      tlw=0.59*z/sigw
    else
      tlw=0.15*h/sigw*(1.-exp(-5*zeta))
    endif

  !*********************
  ! 3. Stable conditions
  !*********************

  else
    sigu=1.e-2+2.*ust*(1.-zeta)   ! Eq. 7.20 Hanna 1982
    sigv=1.e-2+1.3*ust*(1.-zeta)  ! Eq. 7.19 Hanna 1982
    sigw=sigv
    dsigwdz=-1.3*ust/h            ! ???
    tlu=0.15*h/sigu*(sqrt(zeta))  ! Eq. 7.22 Hanna 1982
    tlv=0.467*tlu                 ! Eq. 7.23 Hanna 1982
    tlw=0.1*h/sigw*zeta**0.8      ! Eq. 7.24 Hanna 1982
  endif

  tlu=max(10.,tlu)
  tlv=max(10.,tlv)
  tlw=max(30.,tlw)

  if (dsigwdz.eq.0.) dsigwdz=1.e-10
end subroutine hanna

subroutine hanna1(z)
  !                  i
  !*****************************************************************************
  !                                                                            *
  !   Computation of \sigma_i and \tau_L based on the scheme of Hanna (1982)   *
  !   Source: 'Atmospheric Turbulence and Air Polution', chapter 4 and 7,      *
  !   J.A. Businger, edited by F.T.M. Nieuwstadt and H. van Dop                *
  !                                                                            *
  !   Author: A. Stohl                                                         *
  !                                                                            *
  !   4 December 1997                                                          *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  ! dsigwdz [1/s]     vertical gradient of sigw                                *
  ! ol [m]            Obukhov length                                           *
  ! sigu, sigv, sigw  standard deviations of turbulent velocity fluctuations   *
  ! tlu [s]           Lagrangian time scale for the along wind component.      *
  ! tlv [s]           Lagrangian time scale for the cross wind component.      *
  ! tlw [s]           Lagrangian time scale for the vertical wind component.   *
  ! ust, ustar [m/s]  friction velocity                                        *
  ! wst, wstar [m/s]  convective velocity scale                                *
  !                                                                            *
  !*****************************************************************************

  implicit none

  real :: z,s1,s2

  !**********************
  ! 1. Neutral conditions
  !**********************

  if (h/abs(ol).lt.1.) then

    ust=max(1.e-4,ust)

    ! Eq. 7.25 Hanna 1982: sigu/ust=2.0*exp(-3*f*z/ust),
    ! where f, the Coriolis parameter, is set to 1e-4
    sigu=2.0*ust*exp(-3.e-4*z/ust)
    sigu=max(sigu,1.e-5)

    ! Eq. 7.26 Hanna 1982: sigv/ust=sigw/ust=1.3*exp(-2*f*z/ust),
    ! where f, the Coriolis parameter, is set to 1e-4
    sigv=1.3*ust*exp(-2.e-4*z/ust)
    sigv=max(sigv,1.e-5)
    sigw=sigv

    ! ???
    dsigw2dz=-6.76e-4*ust*exp(-4.e-4*z/ust)

    ! Eq.7.27 Hanna 1982: TL=0.5*z/sigw/(1+15*f*z/ust) assumed to be valid
    ! for all three components
    tlu=0.5*z/sigw/(1.+1.5e-3*z/ust)
    tlv=tlu
    tlw=tlu

  !***********************
  ! 2. Unstable conditions
  !***********************

  else if (ol.lt.0.) then

  ! Determine sigmas
  !*****************

    ! Eq. 4.15 Caughey 1982
    sigu=ust*(12.-0.5*h/ol)**0.33333
    sigu=max(sigu,1.e-6)
    sigv=sigu

    ! Eq. 7.15 Hanna 1982
    if (zeta.lt.0.03) then
      sigw=0.96*wst*(3*zeta-ol/h)**0.33333
      dsigw2dz=1.8432*wst*wst/h*(3*zeta-ol/h)**(-0.33333)
    else if (zeta.lt.0.4) then
      s1=0.96*(3*zeta-ol/h)**0.33333
      s2=0.763*zeta**0.175
      if (s1.lt.s2) then
        sigw=wst*s1
        dsigw2dz=1.8432*wst*wst/h*(3*zeta-ol/h)**(-0.33333)
      else
        sigw=wst*s2
        dsigw2dz=0.203759*wst*wst/h*zeta**(-0.65)
      endif
    else if (zeta.lt.0.96) then
      sigw=0.722*wst*(1-zeta)**0.207
      dsigw2dz=-.215812*wst*wst/h*(1-zeta)**(-0.586)
    else if (zeta.lt.1.00) then
      sigw=0.37*wst
      dsigw2dz=0.
    endif
    sigw=max(sigw,1.e-6)

  ! Determine average Lagrangian time scale
  !****************************************

    ! Eq. 7.17 Hanna  1982
    tlu=0.15*h/sigu
    tlv=tlu
    if (z.lt.abs(ol)) then
      tlw=0.1*z/(sigw*(0.55-0.38*abs(z/ol)))
    else if (zeta.lt.0.1) then
      tlw=0.59*z/sigw
    else
      tlw=0.15*h/sigw*(1.-exp(-5*zeta))
    endif

  !*********************
  ! 3. Stable conditions
  !*********************

  else
    sigu=2.*ust*(1.-zeta)
    sigv=1.3*ust*(1.-zeta)
    sigu=max(sigu,1.e-6)
    sigv=max(sigv,1.e-6)
    sigw=sigv
    dsigw2dz=3.38*ust*ust*(zeta-1.)/h
    tlu=0.15*h/sigu*(sqrt(zeta))
    tlv=0.467*tlu
    tlw=0.1*h/sigw*zeta**0.8
  endif

  tlu=max(10.,tlu)
  tlv=max(10.,tlv)
  tlw=max(30.,tlw)
end subroutine hanna1

subroutine hanna_short(z)
  !                       i
  !*****************************************************************************
  !                                                                            *
  !   Computation of \sigma_i and \tau_L based on the scheme of Hanna (1982)   *
  !                                                                            *
  !   Author: A. Stohl                                                         *
  !                                                                            *
  !   4 December 1997                                                          *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  ! dsigwdz [1/s]     vertical gradient of sigw                                *
  ! ol [m]            Obukhov length                                           *
  ! sigu, sigv, sigw  standard deviations of turbulent velocity fluctuations   *
  ! tlu [s]           Lagrangian time scale for the along wind component.      *
  ! tlv [s]           Lagrangian time scale for the cross wind component.      *
  ! tlw [s]           Lagrangian time scale for the vertical wind component.   *
  ! ust, ustar [m/s]  friction velocity                                        *
  ! wst, wstar [m/s]  convective velocity scale                                *
  !                                                                            *
  !*****************************************************************************

  implicit none

  real :: z

  !**********************
  ! 1. Neutral conditions
  !**********************

  if (h/abs(ol).lt.1.) then
    ust=max(1.e-4,ust)
    sigw=1.3*exp(-2.e-4*z/ust)
    dsigwdz=-2.e-4*sigw
    sigw=sigw*ust+1.e-2
    tlw=0.5*z/sigw/(1.+1.5e-3*z/ust)

  !***********************
  ! 2. Unstable conditions
  !***********************

  else if (ol.lt.0.) then

  ! Determine sigmas
  !*****************

    sigw=sqrt(1.2*wst**2*(1.-.9*zeta)*zeta**0.66666+ &
         (1.8-1.4*zeta)*ust**2)+1.e-2
    dsigwdz=0.5/sigw/h*(-1.4*ust**2+wst**2* &
         (0.8*max(zeta,1.e-3)**(-.33333)-1.8*zeta**0.66666))

  ! Determine average Lagrangian time scale
  !****************************************

    if (z.lt.abs(ol)) then
      tlw=0.1*z/(sigw*(0.55-0.38*abs(z/ol)))
    else if (zeta.lt.0.1) then
      tlw=0.59*z/sigw
    else
      tlw=0.15*h/sigw*(1.-exp(-5*zeta))
    endif

  !*********************
  ! 3. Stable conditions
  !*********************

  else
    sigw=1.e-2+1.3*ust*(1.-zeta)
    dsigwdz=-1.3*ust/h
    tlw=0.1*h/sigw*zeta**0.8
  endif

  tlu=max(10.,tlu)
  tlv=max(10.,tlv)
  tlw=max(30.,tlw)
  if (dsigwdz.eq.0.) dsigwdz=1.e-10
end subroutine hanna_short

subroutine windalign(u,v,ffap,ffcp,ux,vy)
  !                     i i  i    i   o  o
  !*****************************************************************************
  !                                                                            *
  !  Transformation from along- and cross-wind components to u and v           *
  !  components.                                                               *
  !                                                                            *
  !     Author: A. Stohl                                                       *
  !                                                                            *
  !     3 June 1996                                                            *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  ! ffap  turbulent wind in along wind direction                               *
  ! ffcp  turbulent wind in cross wind direction                               *
  ! u     main wind component in x direction                                   *
  ! ux    turbulent wind in x direction                                        *
  ! v     main wind component in y direction                                   *
  ! vy    turbulent wind in y direction                                        *
  !                                                                            *
  !*****************************************************************************

  implicit none

  real :: u,v,ffap,ffcp,ux,vy,ffinv,ux1,ux2,vy1,vy2,sinphi,cosphi
  real,parameter :: eps=1.e-30

  ! Transform along wind components
  !********************************

  ffinv=1./max(sqrt(u*u+v*v),eps)
  sinphi=v*ffinv
  vy1=sinphi*ffap
  cosphi=u*ffinv
  ux1=cosphi*ffap

  ! Transform cross wind components
  !********************************

  ux2=-sinphi*ffcp
  vy2=cosphi*ffcp

  ! Add contributions from along and cross wind components
  !*******************************************************

  ux=ux1+ux2
  vy=vy1+vy2
end subroutine windalign

end module turbulence_mod