scale_atmos_phy_tb_mynn.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
module scale_atmos_phy_tb_mynn
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_index
  use scale_tracer

#ifdef DEBUG
  use scale_debug, only: &
     check
  use scale_const, only: &
     undef  => const_undef, &
     iundef => const_undef2
#endif
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_phy_tb_mynn_setup
  public :: atmos_phy_tb_mynn

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  real(RP), parameter :: oneoverthree = 1.0_rp / 3.0_rp
  real(RP), parameter :: tke_min = 1.0e-10_rp
  real(RP), parameter :: lt_min = 1.0e-6_rp
  real(RP), parameter :: us_min = 1.0e-6_rp

  real(RP)              :: a1
  real(RP)              :: a2
  real(RP), parameter   :: b1 = 24.0_rp
  real(RP), parameter   :: b2 = 15.0_rp
  real(RP)              :: c1
  real(RP), parameter   :: c2 = 0.75_rp
  real(RP), parameter   :: c3 = 0.352_rp
  real(RP), parameter   :: c5 = 0.2_rp
  real(RP), parameter   :: g1 = 0.235_rp
  real(RP)              :: g2
  real(RP)              :: f1
  real(RP)              :: f2
  real(RP)              :: rf1
  real(RP)              :: rf2
  real(RP)              :: rfc
  real(RP)              :: af12
  real(RP), parameter   :: prn = 0.74_rp

  real(RP)              :: sqrt_2pi
  real(RP)              :: rsqrt_2pi
  real(RP)              :: rsqrt_2

  integer :: ke_pbl
  logical :: atmos_phy_tb_mynn_tke_init = .false. 
  real(RP) :: atmos_phy_tb_mynn_n2_max = 1.e3_rp
  real(RP) :: atmos_phy_tb_mynn_nu_min = 1.e-6_rp
  real(RP) :: atmos_phy_tb_mynn_nu_max = 10000._rp
  real(RP) :: atmos_phy_tb_mynn_kh_min = 1.e-6_rp
  real(RP) :: atmos_phy_tb_mynn_kh_max = 10000._rp
  real(RP) :: atmos_phy_tb_mynn_lt_max = 700._rp ! ~ 0.23 * 3 km

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_tb_mynn_setup( &
       TYPE_TB,       &
       CDZ, CDX, CDY, &
       CZ             )
    use scale_process, only: &
       prc_mpistop
    use scale_const, only: &
       pi => const_pi
    implicit none

    character(len=*), intent(in) :: TYPE_TB

    real(RP), intent(in) :: CDZ(ka)
    real(RP), intent(in) :: CDX(ia)
    real(RP), intent(in) :: CDY(ja)
    real(RP), intent(in) :: CZ (ka,ia,ja)

    real(RP) :: ATMOS_PHY_TB_MYNN_PBL_MAX = 1e10_rp 

    namelist / param_atmos_phy_tb_mynn / &
         atmos_phy_tb_mynn_tke_init, &
         atmos_phy_tb_mynn_pbl_max, &
         atmos_phy_tb_mynn_n2_max, &
         atmos_phy_tb_mynn_nu_min, &
         atmos_phy_tb_mynn_nu_max, &
         atmos_phy_tb_mynn_kh_min, &
         atmos_phy_tb_mynn_kh_max, &
         atmos_phy_tb_mynn_lt_max

    integer :: k, i, j
    integer :: ierr
    !---------------------------------------------------------------------------

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '++++++ Module[TURBULENCE] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
    if( io_l ) write(io_fid_log,*) '+++ Mellor-Yamada Nakanishi-Niino Model'

    if ( type_tb /= 'MYNN' ) then
       write(*,*) 'xxx ATMOS_PHY_TB_TYPE is not MYNN. Check!'
       call prc_mpistop
    endif


    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_phy_tb_mynn,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_TB_MYNN. Check!'
       call prc_mpistop
    endif
    if( io_lnml ) write(io_fid_log,nml=param_atmos_phy_tb_mynn)

    do k = ks, ke-1
       do j = js, je
       do i = is, ie
          if ( atmos_phy_tb_mynn_pbl_max >= cz(k,i,j) ) then
             ke_pbl = k
          end if
       end do
       end do
    end do

    a1 = b1 * (1.0_rp - 3.0_rp * g1) / 6.0_rp
    a2 = 1.0_rp / (3.0_rp * g1 * b1**(1.0_rp/3.0_rp) * prn )
    c1 = g1 - 1.0_rp / ( 3.0_rp * a1 * b1**(1.0_rp/3.0_rp) )
    g2 = ( 2.0_rp * a1 * (3.0_rp - 2.0_rp * c2) + b2 * (1.0_rp - c3) ) / b1
    f1 = b1 * (g1 - c1) + 2.0_rp * a1 * (3.0_rp - 2.0_rp * c2) + 3.0_rp * a2 * (1.0_rp - c2) * (1.0_rp - c5)
    f2 = b1 * (g1 + g2) - 3.0_rp * a1 * (1.0_rp - c2)

    rf1 = b1 * (g1 - c1) / f1
    rf2 = b1 * g1 / f2
    rfc = g1 / (g1 + g2)

    af12 = a1 * f1 / ( a2 * f2 )

    sqrt_2pi = sqrt( 2.0_rp * pi )
    rsqrt_2pi = 1.0_rp / sqrt_2pi
    rsqrt_2 = 1.0_rp / sqrt( 2.0_rp )

    return
  end subroutine atmos_phy_tb_mynn_setup

  !-----------------------------------------------------------------------------
  subroutine atmos_phy_tb_mynn( &
       qflx_sgs_momz, qflx_sgs_momx, qflx_sgs_momy, & ! (out)
       qflx_sgs_rhot, qflx_sgs_rhoq,                & ! (out)
       tke,                                         & ! (inout)
       tke_t, nu, ri, pr, n2,                       & ! (out) diagnostic variables
       momz, momx, momy, rhot, dens, qtrc,          & ! (in)
       sflx_mw, sflx_mu, sflx_mv, sflx_sh, sflx_qv, & ! (in)
       gsqrt, j13g, j23g, j33g, mapf, dt            ) ! (in)
    use scale_precision
    use scale_grid_index
    use scale_tracer
    use scale_const, only: &
       grav   => const_grav, &
       r      => const_rdry, &
       rvap   => const_rvap, &
       cp     => const_cpdry, &
       epstvap => const_epstvap
    use scale_grid, only: &
       rcdz => grid_rcdz, &
       rfdz => grid_rfdz, &
       rcdx => grid_rcdx, &
       rcdy => grid_rcdy
    use scale_grid_real, only: &
       cz => real_cz
    use scale_gridtrans, only: &
       i_xyz, &
       i_xyw, &
       i_xvw, &
       i_uyw, &
       i_uyz, &
       i_xvz, &
       i_xy, &
       i_uy, &
       i_xv
    use scale_comm, only: &
       comm_vars8, &
       comm_vars, &
       comm_wait
    use scale_atmos_phy_tb_common, only: &
       diffusion_solver => atmos_phy_tb_diffusion_solver
    use scale_atmos_thermodyn, only: &
       atmos_thermodyn_temp_pres, &
       atmos_thermodyn_templhv, &
       atmos_thermodyn_templhs
    use scale_atmos_saturation, only: &
       atmos_saturation_alpha, &
       atmos_saturation_pres2qsat => atmos_saturation_pres2qsat_all
#ifdef MORE_HIST
    use scale_history, only: &
       hist_in
#endif
    implicit none

    real(RP), intent(out) :: qflx_sgs_momz(ka,ia,ja,3)
    real(RP), intent(out) :: qflx_sgs_momx(ka,ia,ja,3)
    real(RP), intent(out) :: qflx_sgs_momy(ka,ia,ja,3)
    real(RP), intent(out) :: qflx_sgs_rhot(ka,ia,ja,3)
    real(RP), intent(out) :: qflx_sgs_rhoq(ka,ia,ja,3,qa)

    real(RP), intent(inout) :: tke (ka,ia,ja) ! TKE
    real(RP), intent(out) :: tke_t(ka,ia,ja) ! tendency of TKE
    real(RP), intent(out) :: Nu(ka,ia,ja) ! eddy viscosity (center)
    real(RP), intent(out) :: Ri(ka,ia,ja) ! Richardson number
    real(RP), intent(out) :: Pr(ka,ia,ja) ! Plandtle number
    real(RP), intent(out) :: N2(ka,ia,ja) ! squared Brunt-Vaisala frequency

    real(RP), intent(in)  :: MOMZ(ka,ia,ja)
    real(RP), intent(in)  :: MOMX(ka,ia,ja)
    real(RP), intent(in)  :: MOMY(ka,ia,ja)
    real(RP), intent(in)  :: RHOT(ka,ia,ja)
    real(RP), intent(in)  :: DENS(ka,ia,ja)
    real(RP), intent(in)  :: QTRC(ka,ia,ja,qa)

    real(RP), intent(in)  :: SFLX_MW(ia,ja)
    real(RP), intent(in)  :: SFLX_MU(ia,ja)
    real(RP), intent(in)  :: SFLX_MV(ia,ja)
    real(RP), intent(in)  :: SFLX_SH(ia,ja)
    real(RP), intent(in)  :: SFLX_QV(ia,ja)

    real(RP), intent(in)  :: GSQRT   (ka,ia,ja,7)
    real(RP), intent(in)  :: J13G    (ka,ia,ja,7)
    real(RP), intent(in)  :: J23G    (ka,ia,ja,7)
    real(RP), intent(in)  :: J33G
    real(RP), intent(in)  :: MAPF    (ia,ja,2,4)
    real(DP), intent(in)  :: dt

    real(RP) :: U(ka,ia,ja)
    real(RP) :: V(ka,ia,ja)
    real(RP) :: phiN(ka,ia,ja)


    real(RP) :: sm(ka,ia,ja)
    real(RP) :: sh(ka,ia,ja)
    real(RP) :: l(ka,ia,ja)
    real(RP) :: q(ka,ia,ja)
    real(RP) :: dudz2(ka,ia,ja)
    real(RP) :: q2_2(ka,ia,ja)
    real(RP) :: Kh
    real(RP) :: RHOKh (ka,ia,ja)

    real(RP) :: SFLX_PT(ia,ja) ! surface potential temperature flux

    real(RP) :: a(ka,ia,ja)
    real(RP) :: b(ka,ia,ja)
    real(RP) :: c(ka,ia,ja)
    real(RP) :: d(ka,ia,ja)
    real(RP) :: ap
    real(RP) :: tke_N(ka,ia,ja)

    real(RP) :: POTT(ka,ia,ja)
    real(RP) :: POTV(ka,ia,ja)
    real(RP) :: POTL(ka,ia,ja)
    real(RP) :: TEML(ka,ia,ja)

    real(RP) :: Qw(ka,ia,ja)
    real(RP) :: ql
    real(RP) :: qs
    real(RP) :: qdry

    real(RP) :: temp(ka,ia,ja)
    real(RP) :: pres(ka,ia,ja)

    real(RP) :: lh(ka,ia,ja)
    real(RP) :: lhv(ka,ia,ja)
    real(RP) :: lhs(ka,ia,ja)
    real(RP) :: alpha(ka,ia,ja)

    real(RP) :: ac

    real(RP) :: Q1
    real(RP) :: Qsl(ka,ia,ja)
    real(RP) :: dQsl
    real(RP) :: sigma_s
    real(RP) :: RR
    real(RP) :: Rt
    real(RP) :: betat
    real(RP) :: betaq
    real(RP) :: aa
    real(RP) :: bb
    real(RP) :: cc

    real(RP) :: flux_z(ka,ia,ja)
    real(RP) :: flux_x(ka,ia,ja)
    real(RP) :: flux_y(ka,ia,ja)
    real(RP) :: mflx
    real(RP) :: advc

#ifdef MORE_HIST
    real(RP) :: adv(ka,ia,ja)
    real(RP) :: gen(ka,ia,ja)
#endif

    integer :: k, i, j, iq
    integer :: IIS, IIE, JJS, JJE

    if ( io_l ) write(io_fid_log, *) "*** Physics step: Turbulence (MYNN)"

#ifdef DEBUG
    pott(:,:,:) = undef
    potl(:,:,:) = undef
    temp(:,:,:) = undef
#endif

!OCL XFILL
    do j = js  , je
    do i = is-1, ie
    do k = ks-1, ke
       qflx_sgs_momz(k,i,j,xdir) = 0.0_rp
    end do
    end do
    end do
!OCL XFILL
    do j = js  , je
    do i = is  , ie+1
    do k = ks-1, ke+1
       qflx_sgs_momx(k,i,j,xdir) = 0.0_rp
    end do
    end do
    end do

!OCL XFILL
    do j = js  , je
    do i = is-1, ie
    do k = ks-1, ke+1
       qflx_sgs_momy(k,i,j,xdir) = 0.0_rp
    end do
    end do
    end do
!OCL XFILL
    do j = js  , je
    do i = is-1, ie
    do k = ks-1, ke+1
       qflx_sgs_rhot(k,i,j,xdir) = 0.0_rp
    end do
    end do
    end do
!OCL XFILL
    do iq = 1, qa
    do j = js  , je
    do i = is-1, ie
    do k = ks-1, ke+1
       qflx_sgs_rhoq(k,i,j,xdir,iq) = 0.0_rp
    end do
    end do
    end do
    end do

!OCL XFILL
    do j = js-1, je
    do i = is  , ie
    do k = ks-1, ke
       qflx_sgs_momz(k,i,j,ydir) = 0.0_rp
    end do
    end do
    end do
!OCL XFILL
    do j = js-1, je
    do i = is  , ie
    do k = ks-1, ke+1
       qflx_sgs_momx(k,i,j,ydir) = 0.0_rp
    end do
    end do
    end do
!OCL XFILL
    do j = js  , je+1
    do i = is  , ie
    do k = ks-1, ke+1
       qflx_sgs_momy(k,i,j,ydir) = 0.0_rp
    end do
    end do
    end do
!OCL XFILL
    do j = js-1, je
    do i = is  , ie
    do k = ks-1, ke+1
       qflx_sgs_rhot(k,i,j,ydir) = 0.0_rp
    end do
    end do
    end do
!OCL XFILL
    do iq = 1, qa
    do j = js-1, je
    do i = is  , ie
    do k = ks-1, ke+1
       qflx_sgs_rhoq(k,i,j,ydir,iq) = 0.0_rp
    end do
    end do
    end do
    end do
!OCL XFILL
    do j = js, je
    do i = is, ie
    do k = ks-1, ke+1
       qflx_sgs_momz(k,i,j,zdir) = 0.0_rp
    end do
    end do
    end do
!OCL XFILL
    do iq = i_qv+1, qa
    do j = js-1, je
    do i = is  , ie
    do k = ks-1, ke+1
       qflx_sgs_rhoq(k,i,j,zdir,iq) = 0.0_rp
    end do
    end do
    end do
    end do



    do jjs = js, je, jblock
    jje = jjs+jblock-1
    do iis = is, ie, iblock
    iie = iis+iblock-1

!OCL XFILL
       do j = jjs  , jje+1
       do i = iis-1, iie+1
       do k = ks, ke_pbl+1
          u(k,i,j) = 0.5_rp * ( momx(k,i,j) + momx(k,i-1,j) ) / dens(k,i,j)
       end do
       end do
       end do

!OCL XFILL
       do j = jjs-1, jje+1
       do i = iis  , iie+1
       do k = ks, ke_pbl+1
          v(k,i,j) = 0.5_rp * ( momy(k,i,j) + momy(k,i,j-1) ) / dens(k,i,j)
       end do
       end do
       end do

    end do
    end do

    call atmos_thermodyn_temp_pres( temp, pres, & ! (out)
                                    dens, rhot, qtrc ) ! (in)


    call atmos_thermodyn_templhv( lhv, temp )
    call atmos_thermodyn_templhs( lhs, temp )
    call atmos_saturation_alpha( alpha, temp )

!OCL LOOP_NOFUSION,PREFETCH_SEQUENTIAL(SOFT),SWP
    do j = js, je
    do i = is, ie
       do k = ks, ke_pbl+1
          ql = 0.0_rp
          if ( i_qc > 0 ) ql = qtrc(k,i,j,i_qc)
!          do iq = QWS, QWE
!             ql = ql + QTRC(k,i,j,iq)
!          end do
          qs = 0.0_rp
          if ( i_qi > 0 ) qs = qtrc(k,i,j,i_qi)
!          do iq = QIS, QIE
!             qs = qs + QTRC(k,i,j,iq)
!          end do
          qdry = 1.0_rp
          do iq = qqs, qqe
             qdry = qdry - qtrc(k,i,j,iq)
          end do

          qw(k,i,j) = qtrc(k,i,j,i_qv) + ql + qs

          lh(k,i,j) = lhv(k,i,j) * alpha(k,i,j) + lhs(k,i,j) * ( 1.0_rp-alpha(k,i,j) )

          pott(k,i,j) = rhot(k,i,j) / dens(k,i,j)
          ! liquid water potential temperature
          potl(k,i,j) = pott(k,i,j) * (1.0_rp - 1.0_rp * (lhv(k,i,j) * ql + lhs(k,i,j) * qs) / ( temp(k,i,j) * cp ) )
          teml(k,i,j) = potl(k,i,j) * temp(k,i,j) / pott(k,i,j)

          ! virtual potential temperature for derivertive
!          POTV(k,i,j) = ( 1.0_RP + EPSTvap * Qw(k,i,j) ) * POTL(k,i,j)
          potv(k,i,j) = ( qdry + (epstvap+1.0_rp) * qtrc(k,i,j,i_qv) ) * potl(k,i,j)

       end do

       sflx_pt(i,j) = sflx_sh(i,j) / ( cp * dens(ks,i,j) ) &
                    * pott(ks,i,j) / temp(ks,i,j)

       n2(ks,i,j) = min(atmos_phy_tb_mynn_n2_max, &
                        grav * ( potv(ks+1,i,j) - potv(ks,i,j) ) &
                             / ( ( cz(ks+1,i,j)-cz(ks,i,j) ) * potv(ks,i,j) ) )
       do k = ks+1, ke_pbl
          n2(k,i,j) = min(atmos_phy_tb_mynn_n2_max, &
                          grav * ( potv(k+1,i,j) - potv(k-1,i,j) ) &
                               / ( ( cz(k+1,i,j)-cz(k-1,i,j) ) * potv(k,i,j) ) )
       end do

       dudz2(ks,i,j) = ( ( u(ks+1,i,j) - u(ks,i,j) )**2 + ( v(ks+1,i,j) - v(ks,i,j) )**2 ) &
                     / ( cz(ks+1,i,j) - cz(ks,i,j) )**2
       do k = ks+1, ke_pbl
          dudz2(k,i,j) = ( ( u(k+1,i,j) - u(k-1,i,j) )**2 + ( v(k+1,i,j) - v(k-1,i,j) )**2 ) &
                       / ( cz(k+1,i,j) - cz(k-1,i,j) )**2
       end do

       do k = ks, ke_pbl
          ri(k,i,j) = n2(k,i,j) / max(dudz2(k,i,j), 1e-10_rp)
       end do

    end do
    end do

    if ( atmos_phy_tb_mynn_tke_init ) then
!OCL XFILL
       do j = js, je
       do i = is, ie
       do k = ks, ke_pbl
          q(k,i,j) = sqrt( tke_min*2.0_rp )
       end do
       end do
       end do
    else
!OCL XFILL
       do j = js, je
       do i = is, ie
       do k = ks, ke_pbl
          q(k,i,j) = sqrt( max(tke(k,i,j), tke_min)*2.0_rp )
       end do
       end do
       end do
    end if

    ! length
    call get_length( &
         l, & ! (out)
         dens, & ! (in)
         q, n2, & ! (in)
         sflx_mu, sflx_mv, sflx_pt, & ! (in)
         pott ) ! (in)

    call get_q2_level2( &
         q2_2, & ! (out)
         dudz2, ri, l ) ! (in)

    if ( atmos_phy_tb_mynn_tke_init ) then
!OCL XFILL
       do j = js, je
       do i = is, ie
       do k = ks, ke_pbl
          tke(k,i,j) = max(q2_2(k,i,j) * 0.5_rp, tke_min)
          q(k,i,j) = sqrt( tke(k,i,j) * 2.0_rp )
       end do
       end do
       end do
!OCL XFILL
       do j = js, je
       do i = is, ie
       do k = ke_pbl+1, ke
          tke(k,i,j) = 0.0_rp
       end do
       end do
       end do
       call comm_vars8(tke, 1)
       call comm_wait (tke, 1)

       atmos_phy_tb_mynn_tke_init = .false.
    end if

    call get_smsh( sm, sh, & ! (out)
                   q, q2_2, & ! (in)
                   l, n2, dudz2 ) ! (in)

    call atmos_saturation_pres2qsat( qsl(ks:ke_pbl,:,:), & ! (out)
                                     teml(ks:ke_pbl,:,:), pres(ks:ke_pbl,:,:), & ! (in)
                                     ke_pbl-ks+1 ) ! (in)

!OCL LOOP_NOFUSION,PREFETCH_SEQUENTIAL(SOFT),SWP
    do j = js, je
    do i = is, ie
       do k = ks+1, ke_pbl

          dqsl = qsl(k,i,j) * lh(k,i,j) / ( rvap * potl(k,i,j)**2 )
          aa = 1.0_rp / ( 1.0_rp + lh(k,i,j)/cp * dqsl )
          bb = temp(k,i,j) / pott(k,i,j) * dqsl
          ac = min( q(k,i,j)/sqrt(q2_2(k,i,j)), 1.0_rp )
          sigma_s = max( sqrt( 0.25_rp * aa**2 * l(k,i,j)**2 * ac * b2 * sh(k,i,j) ) &
                       * abs( qw(k+1,i,j) - qw(k-1,i,j) - bb * ( potl(k+1,i,j)-potl(k-1,i,j) ) ) &
                       / ( cz(k+1,i,j) - cz(k-1,i,j) ), &
                       1e-10_rp )
          q1 = aa * ( qw(k,i,j) - qsl(k,i,j) ) * 0.5_rp / sigma_s
          rr = min( max( 0.5_rp * ( 1.0_rp + erf(q1*rsqrt_2) ), 0.0_rp ), 1.0_rp )
          ql = min( max( 2.0_rp * sigma_s * ( rr * q1 + rsqrt_2pi * exp(-0.5_rp*q1**2) ), &
                    0.0_rp ), &
                    qw(k,i,j) * 0.5_rp )
          cc = ( 1.0_rp + epstvap * qw(k,i,j) - (1.0_rp+epstvap) * ql ) * pott(k,i,j)/temp(k,i,j) * lh(k,i,j) / cp &
               - (1.0_rp+epstvap) * pott(k,i,j)
          rt = min( max( rr - ql / (2.0_rp*sigma_s*sqrt_2pi) * exp(-q1**2 * 0.5_rp), 0.0_rp ), 1.0_rp )
          betat = 1.0_rp + epstvap * qw(k,i,j) - (1.0_rp+epstvap) * ql - rt * aa * bb * cc
          betaq = epstvap * pott(k,i,j) + rt * aa * cc
          n2(k,i,j) = min(atmos_phy_tb_mynn_n2_max, &
                          grav * ( ( potl(k+1,i,j) - potl(k-1,i,j) ) * betat &
                                 + ( qw(k+1,i,j) - qw(k-1,i,j) ) * betaq ) &
                               / ( cz(k+1,i,j) - cz(k-1,i,j) ) / potv(k,i,j) )
       end do
       n2(ks,i,j) = n2(ks+1,i,j)

       do k = ks, ke_pbl
          ri(k,i,j) = n2(k,i,j) / max(dudz2(k,i,j), 1e-10_rp)
       end do
       do k = ke_pbl+1, ke
          ri(k,i,j) = 0.0_rp
          n2(k,i,j) = 0.0_rp
       end do

    end do
    end do

    ! length
    call get_length( &
         l, & ! (out)
         dens, & ! (in)
         q, n2, & ! (in)
         sflx_mu, sflx_mv, sflx_pt, & ! (in)
         pott ) ! (in)

    call get_q2_level2( &
         q2_2, & ! (out)
         dudz2, ri, l ) ! (in)

    call get_smsh( &
         sm, sh, & ! (out)
         q, q2_2, & ! (in)
         l, n2, dudz2 ) ! (in)


    do j = js, je
    do i = is, ie
       do k = 1, ks-1
          nu(k,i,j) = 0.0_rp
       end do
       do k = ks, ke_pbl
          nu(k,i,j) = max( min( l(k,i,j) * q(k,i,j) * sm(k,i,j), &
                           atmos_phy_tb_mynn_nu_max ), &
                           atmos_phy_tb_mynn_nu_min )
       end do
       do k = ke_pbl+1, ka
          nu(k,i,j) = 0.0_rp
       end do

       do k = ks, ke_pbl
          kh = max( min( l(k,i,j) * q(k,i,j) * sh(k,i,j), &
                    atmos_phy_tb_mynn_kh_max ), &
                    atmos_phy_tb_mynn_kh_min )
          rhokh(k,i,j) = dens(k,i,j) * kh
          pr(k,i,j) = nu(k,i,j) / kh
       end do
       do k = ke_pbl+1, ke
          pr(k,i,j) = 1.0_rp
       end do
    end do
    end do

    call comm_vars( nu,   1 )

    ! time integration

    !  for velocities
!OCL INDEPENDENT
    do j = js, je
    do i = is, ie

       ap = - dt * 0.5_rp * ( dens(ks  ,i,j)*nu(ks  ,i,j) &
                            + dens(ks+1,i,j)*nu(ks+1,i,j) ) &
                          * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
       a(ks,i,j) = ap * rcdz(ks) / ( dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
       c(ks,i,j) = 0.0_rp
       b(ks,i,j) = - a(ks,i,j) + 1.0_rp
       do k = ks+1, ke_pbl-1
          c(k,i,j) = ap * rcdz(k) / ( dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
          ap = - dt * 0.5_rp * ( dens(k  ,i,j)*nu(k  ,i,j) &
                               + dens(k+1,i,j)*nu(k+1,i,j) ) &
                             * rfdz(k) / gsqrt(k,i,j,i_xyw)
          a(k,i,j) = ap * rcdz(k) / ( dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
          b(k,i,j) = - a(k,i,j) - c(k,i,j) + 1.0_rp
       end do
       a(ke_pbl,i,j) = 0.0_rp
       c(ke_pbl,i,j) = ap * rcdz(ke_pbl) / ( dens(ke_pbl,i,j) * gsqrt(ke_pbl,i,j,i_xyz) )
       b(ke_pbl,i,j) = - c(ke_pbl,i,j) + 1.0_rp

       d(ks,i,j) = u(ks,i,j) + dt * sflx_mu(i,j) * rcdz(ks) / ( dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
       do k = ks+1, ke_pbl
          d(k,i,j) = u(k,i,j)
       end do
       call diffusion_solver( &
            phin(:,i,j),                     & ! (out)
            a(:,i,j), b(:,i,j), c(:,i,j), d(:,i,j), & ! (in)
            ke_pbl                           ) ! (in)
    end do
    end do

    call comm_vars( phin, 2 )
    call comm_wait( nu,   1 )
    call comm_wait( phin, 2 )

!OCL INDEPENDENT
    do jjs = js, je, jblock
    jje = jjs+jblock-1
    do iis = is, ie, iblock
    iie = iis+iblock-1

       do j = jjs, jje
       do i = iis, iie
       do k = ks, ke_pbl-1
          qflx_sgs_momx(k,i,j,zdir) = - 0.03125_rp & ! 1/4/4/2
               * ( dens(k,i,j) + dens(k+1,i,j) + dens(k,i+1,j) + dens(k+1,i+1,j) ) &
               * ( nu(k,i,j) + nu(k+1,i,j) + nu(k,i+1,j) + nu(k+1,i+1,j) ) &
               * ( (phin(k+1,i,j)+phin(k+1,i+1,j)) - (phin(k,i,j)+phin(k,i+1,j)) ) &
               * j33g * rfdz(k) / gsqrt(k,i,j,i_uyw)
       end do
       end do
       end do
       do j = jjs, jje
       do i = iis, iie
          qflx_sgs_momx(ks-1,i,j,zdir) = 0.0_rp
       end do
       end do
       do j = jjs, jje
       do i = iis, iie
       do k = ke_pbl, ke
          qflx_sgs_momx(k,i,j,zdir) = 0.0_rp
       end do
       end do
       end do


       ! integration V
!OCL INDEPENDENT
       do j = jjs, jje
       do i = iis, iie
          d(ks,i,j) = v(ks,i,j) + dt * sflx_mv(i,j) * rcdz(ks) / ( dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
          do k = ks+1, ke_pbl
             d(k,i,j) = v(k,i,j)
          end do
          call diffusion_solver( &
               phin(:,i,j),                     & ! (out)
               a(:,i,j), b(:,i,j), c(:,i,j), d(:,i,j), & ! (in)
               ke_pbl                           ) ! (in)
       end do
       end do

    end do
    end do

    call comm_vars( phin, 1 )
    call comm_wait( phin, 1 )

    do jjs = js, je, jblock
    jje = jjs+jblock-1
    do iis = is, ie, iblock
    iie = iis+iblock-1

       do j = jjs, jje
       do i = iis, iie
          qflx_sgs_momy(ks-1,i,j,zdir) = 0.0_rp
          do k = ks, ke_pbl-1
             qflx_sgs_momy(k,i,j,zdir) = - 0.03125_rp & ! 1/4/4/2
               * ( dens(k,i,j) + dens(k+1,i,j) + dens(k,i,j+1) + dens(k+1,i,j+1) ) &
               * ( nu(k,i,j) + nu(k+1,i,j) + nu(k,i,j+1) + nu(k+1,i,j+1) ) &
               * ( (phin(k+1,i,j)+phin(k+1,i,j+1)) - (phin(k,i,j)+phin(k,i,j+1)) ) &
               * j33g * rfdz(k) / gsqrt(k,i,j,i_xvw)
          end do
          do k = ke_pbl, ke
             qflx_sgs_momy(k,i,j,zdir) = 0.0_rp
          end do
       end do
       end do


       ! for scalars
       ! integration POTT
!OCL INDEPENDENT
       do j = jjs, jje
       do i = iis, iie
          ap = - dt * 0.5_rp * ( rhokh(ks,i,j) + rhokh(ks+1,i,j) ) &
             * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
          a(ks,i,j) = ap * rcdz(ks) / (dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
          c(ks,i,j) = 0.0_rp
          b(ks,i,j) = - a(ks,i,j) + 1.0_rp
          do k = ks+1, ke_pbl-1
             c(k,i,j) = ap * rcdz(k) / (dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
             ap = - dt * 0.5_rp * ( rhokh(k,i,j) + rhokh(k+1,i,j) ) &
                * rfdz(k) / gsqrt(k,i,j,i_xyw)
             a(k,i,j) = ap * rcdz(k) / ( dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
             b(k,i,j) = - a(k,i,j) - c(k,i,j) + 1.0_rp
          end do
          a(ke_pbl,i,j) = 0.0_rp
          c(ke_pbl,i,j) = ap * rcdz(ke_pbl) / (dens(ke_pbl,i,j) * gsqrt(ke_pbl,i,j,i_xyz) )
          b(ke_pbl,i,j) = - c(ke_pbl,i,j) + 1.0_rp
          d(ks,i,j) = potl(ks,i,j) + dt * sflx_pt(i,j) * rcdz(ks) / gsqrt(ks,i,j,i_xyz)
          do k = ks+1, ke_pbl
             d(k,i,j) = potl(k,i,j)
          end do
          call diffusion_solver( &
               phin(:,i,j),                     & ! (out)
               a(:,i,j), b(:,i,j), c(:,i,j), d(:,i,j), & ! (in)
               ke_pbl                           ) ! (in)
       end do
       end do

       do j = jjs, jje
       do i = iis, iie
          qflx_sgs_rhot(ks-1,i,j,zdir) = 0.0_rp
          do k = ks, ke_pbl-1
             qflx_sgs_rhot(k,i,j,zdir) = - 0.5_rp & ! 1/2
               * ( rhokh(k,i,j) + rhokh(k+1,i,j) ) &
               * j33g * ( phin(k+1,i,j) - phin(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,i_xyw)
          end do
          do k = ke_pbl, ke
             qflx_sgs_rhot(k,i,j,zdir) = 0.0_rp
          end do
       end do
       end do


       ! integration QV
       iq = i_qv
!OCL INDEPENDENT
       do j = jjs, jje
       do i = iis, iie
          d(ks,i,j) = qw(ks,i,j) + dt * sflx_qv(i,j) * rcdz(ks) / ( dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
          do k = ks+1, ke_pbl
             d(k,i,j) = qw(k,i,j)
          end do
          call diffusion_solver( &
               phin(:,i,j),                     & ! (out)
               a(:,i,j), b(:,i,j), c(:,i,j), d(:,i,j), & ! (in)
               ke_pbl                           ) ! (in)
       end do
       end do

       do j = jjs, jje
       do i = iis, iie
          qflx_sgs_rhoq(ks-1,i,j,zdir,iq) = 0.0_rp
          do k = ks, ke_pbl-1
             qflx_sgs_rhoq(k,i,j,zdir,iq) = - 0.5_rp & ! 1/2
                  * ( rhokh(k,i,j) + rhokh(k+1,i,j) ) &
                  * j33g * ( phin(k+1,i,j) - phin(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,i_xyw)
          end do
          do k = ke_pbl, ke
             qflx_sgs_rhoq(k,i,j,zdir,iq) = 0.0_rp
          end do
       end do
       end do


       ! time integration tke

       do j = jjs, jje
       do i = iis, iie
       do k = ks , ke_pbl-1
          mflx = j33g * momz(k,i,j) / ( mapf(i,j,1,i_xy)*mapf(i,j,2,i_xy) ) &
               + j13g(k,i,j,i_xyw) * 0.25_rp * ( momx(k,i-1,j)+momx(k,i,j)+momx(k+1,i-1,j)+momx(k+1,i,j) ) / mapf(i,j,2,i_xy) &
               + j23g(k,i,j,i_xyw) * 0.25_rp * ( momy(k,i,j-1)+momy(k,i,j)+momy(k+1,i,j-1)+momy(k+1,i,j) ) / mapf(i,j,1,i_xy)
          flux_z(k,i,j) = 0.5_rp * (    mflx  * ( tke(k+1,i,j)+tke(k,i,j) ) &
                                   -abs(mflx) * ( tke(k+1,i,j)-tke(k,i,j) ) )
       end do
       end do
       end do
       do j = jjs  , jje
       do i = iis-1, iie
       do k = ks   , ke_pbl
          mflx = momx(k,i,j) * gsqrt(k,i,j,i_uyz) / mapf(i,j,2,i_uy)
          flux_x(k,i,j) = 0.5_rp * (    mflx  * ( tke(k,i+1,j)+tke(k,i,j) ) &
                                   -abs(mflx) * ( tke(k,i+1,j)-tke(k,i,j) ) )
       end do
       end do
       end do
       do j = jjs-1, jje
       do i = iis  , iie
       do k = ks   , ke_pbl
          mflx = momy(k,i,j) * gsqrt(k,i,j,i_xvz) / mapf(i,j,1,i_xv)
          flux_y(k,i,j) = 0.5_rp * (    mflx  * ( tke(k,i,j+1)+tke(k,i,j) ) &
                                   -abs(mflx) * ( tke(k,i,j+1)-tke(k,i,j) ) )
       end do
       end do
       end do

       do j = jjs, jje
       do i = iis, iie
          advc = ( ( flux_z(ks,i,j)                    ) * rcdz(ks) &
                 + ( flux_x(ks,i,j) - flux_x(ks,i-1,j) ) * rcdx(i) &
                 + ( flux_y(ks,i,j) - flux_y(ks,i,j-1) ) * rcdy(j) ) &
                 * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / ( gsqrt(ks,i,j,i_xyz) * dens(ks,i,j) )
          d(ks,i,j) = tke(ks,i,j) + dt * ( nu(ks,i,j) * (dudz2(ks,i,j) - n2(ks,i,j)/pr(ks,i,j)) &
                                     - advc )
#ifdef MORE_HIST
          adv(ks,i,j) = advc
          gen(ks,i,j) = nu(ks,i,j) * (dudz2(ks,i,j) - n2(ks,i,j)/pr(ks,i,j))
#endif
          do k = ks+1, ke_pbl-1
             advc = ( ( flux_z(k,i,j) - flux_z(k-1,i,j) ) * rcdz(k) &
                    + ( flux_x(k,i,j) - flux_x(k,i-1,j) ) * rcdx(i) &
                    + ( flux_y(k,i,j) - flux_y(k,i,j-1) ) * rcdy(j) ) &
                    * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / ( gsqrt(k,i,j,i_xyz) * dens(k,i,j) )
             d(k,i,j) = tke(k,i,j) &
                  + dt * ( nu(k,i,j) * (dudz2(k,i,j) - n2(k,i,j)/pr(k,i,j)) &
                         - advc )
#ifdef MORE_HIST
             adv(k,i,j) = advc
             gen(k,i,j) = nu(k,i,j) * (dudz2(k,i,j) - n2(k,i,j)/pr(k,i,j))
#endif
          end do
          advc = ( (                    - flux_z(ke_pbl-1,i,j) ) * rcdz(ke_pbl) &
                 + ( flux_x(ke_pbl,i,j) - flux_x(ke_pbl,i-1,j) ) * rcdx(i) &
                 + ( flux_y(ke_pbl,i,j) - flux_y(ke_pbl,i,j-1) ) * rcdy(j) ) &
                 * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / ( gsqrt(ke_pbl,i,j,i_xyz) * dens(ke_pbl,i,j) )
          d(ke_pbl,i,j) = tke(ke_pbl,i,j) + dt * ( nu(ke_pbl,i,j) * (dudz2(ke_pbl,i,j) - n2(ke_pbl,i,j)/pr(ke_pbl,i,j)) &
                                             - advc )
#ifdef MORE_HIST
          adv(ke_pbl,i,j) = advc
          gen(ke_pbl,i,j) = nu(ke_pbl,i,j) * (dudz2(ke_pbl,i,j) - n2(ke_pbl,i,j)/pr(ke_pbl,i,j))
#endif
       end do
       end do

!OCL INDEPENDENT
       do j = jjs, jje
       do i = iis, iie
          ap = - dt * 1.5_rp * ( dens(ks  ,i,j)*nu(ks  ,i,j) &
                               + dens(ks+1,i,j)*nu(ks+1,i,j) ) &
             * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
          a(ks,i,j) = ap * rcdz(ks) / ( dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
          c(ks,i,j) = 0.0_rp
          b(ks,i,j) = - a(ks,i,j) + 1.0_rp + 2.0_rp * dt * q(ks,i,j) / ( b1 * l(ks,i,j) )
          do k = ks+1, ke_pbl-1
             c(k,i,j) = ap * rcdz(k) / ( dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
             ap = - dt * 1.5_rp * ( dens(k  ,i,j)*nu(k  ,i,j) &
                                  + dens(k+1,i,j)*nu(k+1,i,j))  &
                * rfdz(k) / gsqrt(k,i,j,i_xyw)
             a(k,i,j) = ap * rcdz(k) / ( dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
             b(k,i,j) = - a(k,i,j) - c(k,i,j) + 1.0_rp + 2.0_rp * dt * q(k,i,j) / ( b1 * l(k,i,j) )
          end do

          a(ke_pbl,i,j) = 0.0_rp
          c(ke_pbl,i,j) = ap * rcdz(ke_pbl) / ( dens(ke_pbl,i,j) * gsqrt(ke_pbl,i,j,i_xyz) )
          b(ke_pbl,i,j) = - c(ke_pbl,i,j) + 1.0_rp + 2.0_rp * dt * q(ke_pbl,i,j) / ( b1 * l(ke_pbl,i,j) )
          call diffusion_solver( &
               tke_n(:,i,j),     & ! (out)
               a(:,i,j), b(:,i,j), c(:,i,j), d(:,i,j), & ! (in)
               ke_pbl                           ) ! (in)
#ifdef DEBUG
          do k = ks+1, ke_pbl-1
             tke(1,i,j) = d(k,i,j) &
               + ( ( (tke_n(k+1,i,j)-tke_n(k,i,j)) * rfdz(k) / gsqrt(k,i,j,i_xyw) &
                   * (nu(k+1,i,j)*dens(k+1,i,j)+nu(k,i,j)*dens(k,i,j))*1.5_rp &
                   - (tke_n(k,i,j)-tke_n(k-1,i,j)) * rfdz(k-1) / gsqrt(k-1,i,j,i_xyw) &
                   * (nu(k,i,j)*dens(k,i,j)+nu(k-1,i,j)*dens(k-1,i,j))*1.5_rp ) &
                 * rcdz(k) / gsqrt(k,i,j,i_xyz) / dens(k,i,j) &
                 - 2.0_rp*tke_n(k,i,j) * q(k,i,j) / (b1 * l(k,i,j)) ) &
                 * dt
             if ( tke_n(k,i,j) > 0.1_rp .AND. abs(tke(1,i,j) - tke_n(k,i,j))/tke_n(k,i,j) > 1e-10_rp ) then
                advc = ( tke(k,i,j) - d(k,i,j) ) / dt + nu(k,i,j) * (dudz2(k,i,j) - n2(k,i,j)/pr(k,i,j))
                write(*,*)k,i,j,tke(1,i,j),tke_n(k,i,j), tke(k,i,j),nu(k,i,j),dudz2(k,i,j),n2(k,i,j),pr(k,i,j),advc
                open(90, file="mynn.dat")
                write(90,*)ke_pbl-ks+1
                write(90,*)dt, b1
                write(90,*)tke(ks:ke_pbl,i,j)
                write(90,*)q(ks:ke_pbl,i,j)
                write(90,*)l(ks:ke_pbl,i,j)
                write(90,*)rcdz(ks:ke_pbl)
                write(90,*)rfdz(ks:ke_pbl)
                write(90,*)nu(ks:ke_pbl,i,j)
                write(90,*)dens(ks:ke_pbl,i,j)
                write(90,*)dudz2(ks:ke_pbl,i,j)
                write(90,*)n2(ks:ke_pbl,i,j)
                write(90,*)pr(ks:ke_pbl,i,j)
                close(90)
                call abort
             end if
          end do
#endif
          do k = ks, ke_pbl
             tke_t(k,i,j) = ( max(tke_n(k,i,j), tke_min) - tke(k,i,j) ) / dt
          end do
          do k = ke_pbl+1, ke
             tke_t(k,i,j) = 0.0_rp
          end do

       end do
       end do
    end do
    end do

#ifdef MORE_HIST
    adv(ke_pbl+1:ke,:,:) = 0.0_rp
    gen(ke_pbl+1:ke,:,:) = 0.0_rp
    dudz2(ke_pbl+1:ke,:,:) = 0.0_rp
    l(ke_pbl+1:ke,:,:) = 0.0_rp
    potv(ke_pbl+1:ke,:,:) = 0.0_rp
    potl(ke_pbl+1:ke,:,:) = 0.0_rp
    call hist_in(adv, 'TKE_advc', 'advection of TKE', 'm2/s3', nohalo=.true.)
    call hist_in(gen, 'TKE_gen', 'generation of TKE', 'm2/s3', nohalo=.true.)
    call hist_in(dudz2, 'dUdZ2', 'dudz2', 'm2/s2', nohalo=.true.)
    call hist_in(l, 'L_mix', 'minxing length', 'm', nohalo=.true.)
    call hist_in(potv, 'POTV', 'virtual potential temperature', 'K', nohalo=.true.)
    call hist_in(potl, 'POTL', 'liquid potential temperature', 'K', nohalo=.true.)
#endif

    return
  end subroutine atmos_phy_tb_mynn

  subroutine get_length( &
       l, &
       DENS, &
       q, n2, &
       SFLX_MU, SFLX_MV, SFLX_PT, &
       PT0 )
    use scale_grid_real, only: &
       fz => real_fz
    use scale_const, only: &
       grav   => const_grav, &
       karman => const_karman, &
       cp     => const_cpdry, &
       eps    => const_eps
    implicit none

    real(RP), intent(out) :: l(ka,ia,ja)
    real(RP), intent(in) :: DENS(ka,ia,ja)
    real(RP), intent(in) :: q(ka,ia,ja)
    real(RP), intent(in) :: n2(ka,ia,ja)
    real(RP), intent(in) :: SFLX_MU(ia,ja)
    real(RP), intent(in) :: SFLX_MV(ia,ja)
    real(RP), intent(in) :: SFLX_PT(ia,ja)
    real(RP), intent(in) :: PT0(ka,ia,ja)

    real(RP) :: ls
    real(RP) :: lt
    real(RP) :: lb
    real(RP) :: rlm
    real(RP) :: rlt

    real(RP) :: qc
    real(RP) :: int_q
    real(RP) :: int_qz
    real(RP) :: rn2sr
    real(RP) :: us
    real(RP) :: zeta

    real(RP) :: z
    real(RP) :: qdz

    real(RP) :: sw
    integer :: k, i, j

!OCL LOOP_NOFUSION,PREFETCH_SEQUENTIAL(SOFT),SWP
    do j = js, je
    do i = is, ie
       int_qz = 0.0_rp
       int_q = 0.0_rp
       do k = ks, ke_pbl
          qdz = q(k,i,j) * ( fz(k,i,j) - fz(k-1,i,j) )
          int_qz = int_qz + ((fz(k,i,j)+fz(k-1,i,j))*0.5_rp-fz(ks-1,i,j)) * qdz
          int_q  = int_q + qdz
       end do
       ! LT
       lt = min( max(0.23_rp * int_qz / (int_q + eps), &
                    lt_min), &
                    atmos_phy_tb_mynn_lt_max )
       rlt = 1.0_rp / lt

       us = ( sflx_mu(i,j)**2 + sflx_mv(i,j)**2 )**0.25_rp / dens(ks,i,j) ! friction velocity
       us = max(us, us_min)
       rlm = - karman * grav * sflx_pt(i,j) / (pt0(ks,i,j) * us**3 )

       qc = (grav/pt0(ks,i,j)*max(sflx_pt(i,j),0.0_rp)*lt)**oneoverthree

       do k = ks, ke_pbl
          z = ( fz(k,i,j)+fz(k-1,i,j) )*0.5_rp - fz(ks-1,i,j)
          zeta = z * rlm

          ! LS
          sw = sign(0.5_rp, zeta) + 0.5_rp ! 1 for zeta >= 0, 0 for zeta < 0
          ls = karman * z &
             * ( sw / (1.0_rp + 2.7_rp*min(zeta,1.0_rp)*sw ) &
               + ( (1.0_rp - 100.0_rp*zeta)*(1.0_rp-sw) )**0.2_rp )

          ! LB
          sw  = sign(0.5_rp, n2(k,i,j)) + 0.5_rp ! 1 for dptdz >0, 0 for dptdz < 0
          rn2sr = 1.0_rp / ( sqrt(n2(k,i,j)*sw) + 1.0_rp-sw)
          lb = (1.0_rp + 5.0_rp * sqrt(qc*rn2sr/lt)) * q(k,i,j) * rn2sr * sw & ! qc=0 when SFLX_PT < 0
             +  999.e10_rp * (1.0_rp-sw)

          ! L
          l(k,i,j) = 1.0_rp / ( 1.0_rp/ls + rlt + 1.0_rp/lb )
       end do
    end do
    end do

    return
  end subroutine get_length

  subroutine get_q2_level2( &
       q2_2,        &
       dudz2, Ri, l )
    implicit none

    real(RP), intent(out) :: q2_2(ka,ia,ja)
    real(RP), intent(in)  :: dudz2(ka,ia,ja)
    real(RP), intent(in)  :: Ri(ka,ia,ja)
    real(RP), intent(in)  :: l(ka,ia,ja)

    real(RP) :: rf
    real(RP) :: sm_2
    real(RP) :: sh_2

    real(RP) :: q2
    integer :: k, i, j

!OCL LOOP_NOFUSION,PREFETCH_SEQUENTIAL(SOFT),SWP
    do j = js, je
    do i = is, ie
    do k = ks, ke_pbl
       rf = min(0.5_rp / af12 * ( ri(k,i,j) &
                              + af12*rf1 &
                              - sqrt(ri(k,i,j)**2 + 2.0_rp*af12*(rf1-2.0_rp*rf2)*ri(k,i,j) + (af12*rf1)**2) ), &
                rfc)
       sh_2 = 3.0_rp * a2 * (g1+g2) * (rfc-rf) / (1.0_rp-rf)
       sm_2 = sh_2 * af12 * (rf1-rf) / (rf2-rf)
       q2 = b1 * l(k,i,j)**2 * sm_2 * (1.0_rp-rf) * dudz2(k,i,j)
       q2_2(k,i,j) = max( q2, 1.e-10_rp )
    end do
    end do
    end do

    return
  end subroutine get_q2_level2

  subroutine get_smsh( &
         sm, sh, & ! (out)
         q, q2_2, & ! (in)
         l, n2, dudz2 ) ! (in)
    implicit none
    real(RP), intent(out) :: sm(ka,ia,ja)
    real(RP), intent(out) :: sh(ka,ia,ja)
    real(RP), intent(in)  :: q(ka,ia,ja)
    real(RP), intent(in)  :: q2_2(ka,ia,ja)
    real(RP), intent(in)  :: l(ka,ia,ja)
    real(RP), intent(in)  :: n2(ka,ia,ja)
    real(RP), intent(in)  :: dudz2(ka,ia,ja)

    real(RP) :: l2q2
    real(RP) :: ac
    real(RP) :: ac2
    real(RP) :: p1
    real(RP) :: p2
    real(RP) :: p3
    real(RP) :: p4
    real(RP) :: p5
    real(RP) :: rd25
    real(RP) :: gh

    integer :: k, i, j

!OCL LOOP_NOFUSION,PREFETCH_SEQUENTIAL(SOFT),SWP
    do j = js, je
    do i = is, ie
    do k = ks, ke_pbl

       ! level 2.5
       ac = min(q(k,i,j)/sqrt(q2_2(k,i,j)), 1.0_rp)
       ac2 = ac**2
       l2q2 = ( l(k,i,j) / q(k,i,j) )**2
       gh = - n2(k,i,j) * l2q2

       p1 = 1.0_rp - 3.0_rp * ac2 * a2 * b2 * (1.0_rp-c3) * gh
       p2 = 1.0_rp - 9.0_rp * ac2 * a1 * a2 * (1.0_rp-c2) * gh
       p3 = p1 + 9.0_rp * ac2 * a2**2 * (1.0_rp-c2) * (1.0_rp-c5) * gh
       p4 = p1 - 12.0_rp * ac2 * a1 * a2 * (1.0_rp-c2) * gh
       p5 = 6.0_rp * ac2 * a1**2 * dudz2(k,i,j) * l2q2

       rd25 = 1.0_rp / max(p2 * p4 + p5 * p3, 1.e-20_rp)
       sm(k,i,j) = max( ac * a1 * (p3 - 3.0_rp * c1 * p4) * rd25, 0.0_rp )
       sh(k,i,j) = max( ac * a2 * (p2 + 3.0_rp * c1 * p5) * rd25, 0.0_rp )

    end do
    end do
    end do

    return
  end subroutine get_smsh

  function erf(x)
    real(RP), parameter :: a = 0.1400122886866665_rp     ! -8(pi-3)/(3pi(pi-4))
    real(RP), parameter :: fourpi = 1.2732395447351628_rp ! 4/pi

    real(RP), intent(in) :: x
    real(RP) :: erf

    real(RP) :: x2

    x2 = x**2
    erf = sign( sqrt( 1.0_rp - exp(-x2 * (fourpi+a*x2)/(1.0_rp+a*x2) ) ), x )

    return
  end function erf

end module scale_atmos_phy_tb_mynn