module mod_ppm
   !! The **one-dimensional implementation** of the PPM
   !! (piecewise parabolic method)

   implicit none
   private

#ifdef REAL_KIND
   integer, parameter :: fp = REAL_KIND
#else
   integer, parameter :: fp = 8
#endif

   real(fp), parameter :: eps = epsilon(1.0_fp) !! 极小值

   public :: fp, cal_cfl_time_step, adv1d_by_ppm

contains

   subroutine cal_cfl_time_step(dx, u, dt, dt_, nt)
      !! 计算满足 CFL(Courant–Friedrichs–Lewy) 条件的时间积分步长
      real(fp), intent(in)  :: dx(:) !! 网格分辨率: m
      real(fp), intent(in)  :: u(:)  !! 风速: m/s
      real(fp), intent(in)  :: dt    !! 全局迭代步长: s
      real(fp), intent(out) :: dt_   !! 局部迭代步长: s
      integer, intent(out)  :: nt    !! 迭代多少次

      ! local variable
      real(fp), parameter :: factor = 0.5_fp !! 满足CFL条件的程度 < 0.5

      ! 计算 CFL 约束的最小时间步长
      dt_ = factor * minval(dx/abs(u))

      ! 确保 dt_ 能整除 dt，并计算迭代次数
      if (dt_ < dt) then
         nt = int(dt/dt_) + 1
         dt_ = dt/nt
      else
         dt_ = dt
         nt  = 1
      end if
   end subroutine cal_cfl_time_step

   subroutine adv1d_by_ppm(dt, u, c, advc, dx, ds, volume)
      !! 基于PPM方法实现的一维平流函数

      ! input args
      real(fp), intent(in) :: dt !! 时间间隔: s
      real(fp), intent(in) :: u(:) !! 网格右边界风速: m/s
      ! output args
      real(fp), intent(inout) :: c(:)    !! 网格浓度: ug/m3
      real(fp), intent(out)   :: advc(:) !! 浓度变化: ug/m3
      ! optional args: dx 和 ds必须有一个
      real(fp), optional, intent(in) :: dx    !! 网格分辨率: m
      real(fp), optional, intent(in) :: ds(:) !! 网格分辨率: m
      real(fp), optional, intent(in) :: volume(:) !! 体积校正因子，dx*dy*dz;

      ! local variables
      ! mass grid
      real(fp) :: c6(size(u)) !! 决定了抛物线的凸凹特征
      real(fp) :: deltaC(size(u)) !! rightC - leftC

      ! stag grid: 第i个mass网格的右边为stag网格的i
      real(fp) :: egdeC(0:size(u)) !! 格子的边界处的浓度: ug/m3
      real(fp) :: leftC(size(u))   !! 网格的左边界处的浓度 leftC = (0:n-1): ug/m3
      real(fp) :: rightC(size(u))  !! 网格的右边界处的浓度 rightC = (1:n): ug/m3
      real(fp) :: fluidC(size(u))  !! 边界处的浓度随着风场输出(一小段距离)的平均浓度: ug/m3
      real(fp) :: flux(0:size(u))  !! 边界处的通量: ug/m^2/s

      real(fp) :: CFL !! Courant-Friedrichs-Lewy条件（u*dt/dx）
      real(fp) :: fL  !! 网格的左边界处通量: ug/m^2/s
      real(fp) :: fR  !! 网格的右边界处通量: ug/m^2/s
      integer :: i, n

      n = size(u)

      if (present(dx)) then
         call cal_conc_at_uniform_mesh_edge(n, c, egdeC)
      else if (present(ds)) then
         call cal_ppm_non_uniform_mesh_edge_conc(n, ds, c, egdeC)
      else
         write(*,*) 'You must give a dx'
         stop 1
      end if
   
      leftC = egdeC(0:n-1)
      rightC = egdeC(1:n)

      ! step 3: 计算抛物线参数，并对边界处浓度进行约束
      do i=1, n
         ! 更新边界处浓度
         if ((rightC(i) - c(i))*(c(i) - leftC(i)) > 0.) then
            deltaC(i) = rightC(i) - leftC(i)
            c6(i) = 6.*(c(i) - (rightC(i) + leftC(i))/2.)
            ! Overshoot cases
            if (deltaC(i)*c6(i) > deltaC(i)*deltaC(i)) then
               leftC(i) = 3.*c(i) - 2.*rightC(i)
            else if (-deltaC(i)*deltaC(i) > deltaC(i)*c6(i)) then
               rightC(i) = 3.*c(i) - 2.*leftC(i)
            endif
         else ! 局地极值
            leftC(i) = c(i)
            rightC(i) = c(i)
         end if
         ! 确定每个格子抛物线参数：rightC, leftC, c6;
         deltaC(i) = rightC(i) - leftC(i)
         c6(i) = 6.*(c(i) - (rightC(i) + leftC(i))/2.)
      end do

      ! step 4: Compute fluxes
      ! fluxes from boundary cells
      flux(0) = 0 ! surface boundary: impermeable boundary condition
      if (u(1) >= 0.) flux(1) = u(1) * rightC(1) ! First order
      if (u(1) <  0.) flux(1) = u(1) * leftC(1+1) ! First order
      if (u(n-1) >= 0.) flux(n-1) = u(n-1) * rightC(n-1) ! First order
      if (u(n-1) <  0.) flux(n-1) = u(n-1) * leftC(n) ! First order
      if (u(n) > 0.) flux(n) = u(n) * rightC(n) ! 顶边界, 只出去不进

      ! fluxes from the parabolic distribution
      CFL = 0.5 ! 初始化，没有用
      do i = 2, n-2
         if (present(dx)) then
            CFL = abs((u(i)*dt)/dx) ! Guarantee y = u(i)*dt > 0
         else if (present(ds)) then
            CFL = abs((u(i)*dt)/ds(i)) ! Guarantee y = u(i)*dt > 0
            ! 这里的ds, 是不是应该是 (ds(i) + ds(i))/2 
         end if
         if (u(i) > 0.) then ! 公式 2.12
            fluidC(i) = rightC(i) - CFL/2 * (deltaC(i) - (1. - CFL*2./3.)*c6(i))
         else
            fluidC(i) = leftC(i+1) + CFL/2 * (deltaC(i+1) + (1. - CFL*2./3.)*c6(i+1))
         end if
         flux(i) = u(i) * fluidC(i) ! 公式2.14
      end do

      ! step 5: 时间积分, 前向欧拉
      do i = 2, n-1
         fR = flux(i)
         fL = flux(i-1)
         if (present(volume)) then ! 进行校正体积校正，保证质量守恒
            if (fL > 0 .and. i>1) fL = fL*volume(i-1)/volume(i)
            if (fR < 0 .and. i<n) fR = fR*volume(i+1)/volume(i)
         end if
         if (present(dx)) then
            advc(i) = dt/dx * (fL - fR)
         else if (present(ds)) then
            advc(i) = dt/ds(i) * (fL - fR)
         end if
         c(i) = c(i) + advc(i)
      end do

      ! 极小值，数值误差
      where(c < eps) c = 0.

   end subroutine adv1d_by_ppm

   subroutine cal_conc_at_uniform_mesh_edge(n, c, egdeC)
      !! 计算均匀网格边界处的浓度(二次抛物线)
      ! input args
      integer, intent(in)  :: n     !! 网格数
      real(fp), intent(in) :: c(n)  !! 网格浓度: ug/m3

      ! output args
      real(fp), intent(out) :: egdeC(0:n) !! conc at mesh egde: ug/m3

      ! local variables
      integer :: I
      real(fp) :: slopC(n)  !! average slop

     ! step 1: 计算 2~n-1 个网格的平均坡度 (slopC)
      do i = 2, n-1
         ! Compute average slope in the i'th cell
         slopC(i) = 0.5 * (c(i+1) - c(i-1))
         ! Constraint: 确保网格边界处浓度在左右格子平均浓度之间
         if ((c(i+1) - c(i))*(c(i) - c(i-1)) > 0.) then
            slopC(i) = sign(1._fp, slopC(i))* &
                       min(abs(slopC(i)), &
                       2.*abs(c(i+1) - c(i)), &
                       2.*abs(c(i) - c(i-1)))
         else ! 局地极值情况
            slopC(i) = 0
         end if
      end do

      ! step 2: Compute concentrations at cell boundaries
      ! Zero order polynomial
      egdeC(0) = c(1)
      egdeC(n) = c(n)
      ! First order polynomial
      egdeC(1) = (c(2) + c(1))/2.
      egdeC(n-1) = (c(n) + c(n-1))/2.
      ! 求边界处的值
      do i = 2, n-2 ! 中间网格
         egdeC(i) = c(i) + 0.5*(c(i+1) - c(i)) + (slopC(i) - slopC(i+1))/6.
      end do

   end subroutine cal_conc_at_uniform_mesh_edge

   subroutine cal_ppm_non_uniform_mesh_edge_conc(n, dx, c, egdeC)
      !! 计算变网格边界处浓度((二次抛物线))
      ! input args
      integer, intent(in)  :: n     !! 网格数
      real(fp), intent(in) :: dx(n) !! 分辨率: m
      real(fp), intent(in) :: c(n)  !! 网格浓度: ug/m3

      ! output args
      real(fp), intent(out) :: egdeC(0:n)  !! conc at mesh egde: ug/m3

      ! local variables
      integer :: I
      real(fp) :: x0, x1, x2, x3 ! temp lattice vars.
      real(fp) :: alpha (2:n-1)  ! temp lattice var.
      real(fp) :: chi   (2:n-1)  ! lattice var. for dc
      real(fp) :: psi   (2:n-1)  ! lattice var. for dc
      real(fp) :: mu    (2:n-1)  ! lattice var. for cm
      real(fp) :: nu    (2:n-1)  ! lattice var. for cm
      real(fp) :: lambda(2:n-1)  ! lattice var. for cm
      real(fp) :: slopC(n)

      ! step 1: args
      do i = 2, n-1
         alpha(i) = dx(i) + dx(i+1)

         x0 = dx(i-1) + dx(i)
         x2 = dx(i) / (x0 + dx(i+1))
         chi(i) = x2 * (dx(i-1) + x0) / alpha(i)
         psi(i) = x2 * (alpha(i) + dx(i+1)) / x0
      end do
      do i = 2, n-2
         x1 = dx(i)/alpha(i)
         x2 = 2.0_fp * dx(i+1) / alpha(i)
         x3 = 1.0_fp / (dx(i-1) + alpha(i) + dx(i+2))
         mu(i) = x3 * dx(i) * (dx(i-1) + dx(i)) / (dx(i) + alpha(i))
         nu(i) = x3 * dx(i+1) * (dx(i+1) + dx(i+2)) / (dx(i+1) + alpha(i))
         lambda(i) = x1 + mu(i) * x2 - 2.0*nu(i) * x1
      end do

      ! step 1: 平均坡度 (slopC)
      do i = 2, n-1
         ! Compute average slope in the i'th cell
         slopC(i) = chi(i)*(c(i+1) - c(i)) + psi(i)*(c(i) - c(i-1)) ! equation (1.7)
         ! Constraint: 确保梯度不要太大
         if ((c(i+1) - c(i))*(c(i) - c(i-1)) > 0.) then ! equation (1.8)
            slopC(i) = sign(1._fp, slopC(i))* &
                       min( abs(slopC(i)), &
                            2.*abs(c(i+1) - c(i)), &
                            2.*abs(c(i) - c(i-1)) )
         else ! 局地极值情况
            slopC(i) = 0
         end if
      end do

      ! step 2: Compute concentrations at cell boundaries
      egdeC(0) = c(1) ! Zero order polynomial
      egdeC(n) = c(n)
      egdeC(1) = (dx(1)*c(2) + dx(2)*c(1)) / (dx(1) + dx(2)) ! First order polynomial
      egdeC(n-1) = (dx(n-1)*c(n) + dx(n)*c(n-1)) / (dx(n-1) + dx(n))

      ! second order polynomial inside the domain
      do i = 2, n-1 ! equation (1.6)
         egdeC(i) = c(i) + lambda(i)*(c(i+1) - c(i)) - mu(i)*slopC(i+1) + nu(i)*slopC(i)
      end do

   end subroutine cal_ppm_non_uniform_mesh_edge_conc
   
end module mod_ppm
