periodicEBcorr Subroutine

  subroutine periodicEBcorr
  ! added by cew216
! ! use a volume sink to counter a  heat/moisture flux from the SEB
! ! sink acts above buildings
  !
  ! use initfac, only :  max_height_index
  use modfields, only : thlp, qtp
  !use modglobal, only: ltempeq, lperiodicEBcorr, ib, ie, jb, je, kb, ke, imax, jtot
  use modglobal, only : ltempeq, lmoist, lperiodicEBcorr, ib, ie, jb, je, ke,&
                          itot, jtot, totheatflux,sinkbase, totqflux, &
                          zh, fraction,xlen,ylen,dzf
  use modmpi, only : comm3d, mpierr, MY_REAL, MPI_SUM
  !
  integer :: i, j, k
  real :: tot_Tflux, tot_qflux, latent_heat_out,R_theta,R_q, H_proj, E_proj, R_theta_scaled,R_q_scaled, abl_height,phi_theta_t,phi_q_t,sink_depth  !, !tot_qflux !, sink_points
  !
  !write(*,*) 'lperiodicEBcorr ', lperiodicEBcorr
  !write(*,*) 'fraction', fraction
  if (lperiodicEBcorr .eqv. .false.) return

  !
  !call MPI_ALLREDUCE(bctfluxsum,   tot_Tflux,1,MY_REAL,MPI_SUM,comm3d,mpierr)
  !call MPI_ALLREDUCE(bcqfluxsum,   tot_qflux,1,MY_REAL,MPI_SUM,comm3d,mpierr)
  call MPI_ALLREDUCE(totheatflux,tot_Tflux,1,MY_REAL,MPI_SUM,comm3d,mpierr)
  call MPI_ALLREDUCE(totqflux,tot_qflux,1,MY_REAL,MPI_SUM,comm3d,mpierr)
  ! Grylls 2021;  R=(phitop-phibot)/l= -phibot/hABL
  ! Since tot_Tflux = phibot/LxLydeltaV
  ! according to  M if use new tot_Tflux then we have phibot = tot_heatflux/lxly
  ! then define phitop = (1-frac)phibot
  ! R = frac*(phitop-phibot)/l
  ! then needs scaling because of the sink based! code this and comment it.
  !name sinkbase with a k so we now its a vertical index
  !!!! The point is to define phitop phibot and R above the loop so the stuff in the loop looks like the eqns.
  ! Do the same for humidity
  ! This follows the work in Grylls 2021
  H_proj = tot_Tflux/(xlen*ylen) ! [Kms^-1]This is total heat flux in divided by the domain cross section.
  E_proj = tot_qflux/(xlen*ylen)
  abl_height = zh(ke+1)/fraction ! We reverse engineer the ABL height from domain height and fraction
  R_theta = H_proj/abl_height ![Ks^-1] This is the forcing F\theta from Grylls 2021
  R_q = E_proj/abl_height
  ! R is only applied above the canopy, between zh(sinkbase+1) and zh(ke+1), so it is scaled up
  ! by (domain depth)/(sink-layer depth) to keep its vertical integral equal to fraction*H_proj.
  ! This is a ratio of heights, not of level counts: the two agree (ke/(ke-sinkbase)) only on an
  ! equidistant grid, and the level count over- or under-removes heat on a stretched grid.
  sink_depth = zh(ke+1) - zh(sinkbase+1)
  R_theta_scaled = R_theta * zh(ke+1)/sink_depth ! [Ks^-1]
  R_q_scaled = R_q * zh(ke+1)/sink_depth
  !phi_theta_t = 0 ! For debugging the flux profile !(1-fraction)*H_proj! The heat flux out the top of the domain.
  phi_theta_t = (1-fraction)*tot_Tflux/(xlen*ylen*dzf(ke)) ![Ks^-1]    
  phi_q_t = (1-fraction)*tot_qflux/(xlen*ylen*dzf(ke))
  !phi_theta_t = (1-fraction)*H_proj
  !phi_q_t = (1-fraction)*E_proj

   if (ltempeq) then
     ! i innermost (column-major; see #330)
     do k = sinkbase +1 , ke!max_height_index +1 , ke ! Only apply the correction over the volume above the buidlings
       do j = jb,je
         do i = ib,ie
           !thlp(i,j,k) = thlp(i,j,k) + fraction*tot_Tflux*(zh(k+1)-zh(k))/(imax*jtot*(zh(ke+1) - zh(max_height_index+1)))
           !thlp(i,j,k) = thlp(i,j,k) - fraction*tot_Tflux/(itot*jtot*(ke-sinkbase)) ! Most recent working version pre M changes cew216 20240112
           thlp(i,j,k) = thlp(i,j,k) + R_theta_scaled
         end do
       end do
     end do
   !end if
  !sensible_heat_out = (1-fraction)*tot_Tflux/(itot*jtot)
    do j = jb,je
      do i = ib,ie
        thlp(i,j,ke) = thlp(i,j,ke) + phi_theta_t
      end do
    end do
  end if
  !
  !
  !
  if (lmoist) then
    ! i innermost (column-major; see #330)
    do k = sinkbase +1,ke ! Only apply the correction over the volume above the buidlings
      do j = jb,je
        do i = ib,ie
          !qtp(i,j,k) = qtp(i,j,k) + fraction*tot_qflux*(zh(k+1)-zh(k))/(imax*jtot*(zh(ke+1) - zh(max_height_index+1)))
          !qtp(i,j,k) = qtp(i,j,k) - fraction*tot_qflux/(itot*jtot*(ke-sinkbase))
          qtp(i,j,k) = qtp(i,j,k) + R_q_scaled
        end do
      end do
    end do
    latent_heat_out = (1-fraction)*tot_qflux/(itot*jtot*(zh(ke+1)-zh(ke)))
    do j = jb,je
      do i = ib,ie
        qtp(i,j,ke) = qtp(i,j,ke) + phi_q_t
      end do
    end do
  end if
  !
  !write(*,*) 'fraction', fraction
  end subroutine periodicEBcorr