!WRF:MEDIATION_LAYER:IND_HEAT
! Module and subroutines added by dmbrown1
! This is the industrial heat addition section

module module_ind_heat

use module_model_constants, only: cp,xlv

contains

SUBROUTINE ind_heat(grid , config_flags                             &
          ,ids,ide, kds,kde, jds,jde                                &
            ,ims,ime, kms,kme, jms,jme                              &
            ,ips,ipe, kps,kpe, jps,jpe                              &
            ,rho, dz8w, rthindten)

     USE module_domain, only: domain , get_ijk_from_subgrid
     USE module_configure, only: grid_config_rec_type

     implicit none

     TYPE(domain) , TARGET  :: grid ! The main grid that holds everything.  

     TYPE (grid_config_rec_type) , INTENT(IN) :: config_flags ! namelist.input

     integer, intent(in):: &
             ids,ide, kds,kde, jds,jde                              &
            ,ims,ime, kms,kme, jms,jme                              &
            ,ips,ipe, kps,kpe, jps,jpe ! array indices

     ! input variables: density and thickness.  
     real,dimension(ims:ime, kms:kme, jms:jme), intent(in)::rho,dz8w
     real,dimension(ims:ime, kms:kme, jms:jme), intent(out) :: rthindten

     integer :: its,ite,jts,jte,kts,kte   ! atmospheric tiles
     integer :: ij, i, j, k, kk           ! General array indices
     integer :: lu, lu_idx   ! The current land use type and the land use type to add the heat to.
     real :: heat            ! The amount of industrial heat to add to the system. 

     real cp_i, rho_i, mut ! inverse of specific heat and density and total mass.


    heat = config_flags%ind_heat_amt    ! The heating comes from namelist.input.
    lu_idx = config_flags%lu_type       ! As does the lu type to add heat to.  

    do ij=1, grid%num_tiles

      ! Here we arrange the tiles
      its = max(grid%i_start(ij), ids)
      ite = min(grid%i_end(ij), ide-1)
      jts = max(grid%j_start(ij), jds)
      jte = min(grid%j_end(ij), jde-1)
      kts = kds
      kte = kde
      
      cp_i = 1./cp      ! Specific heat inverse

      ! We initialise the rthindten variable to zero everywhere
      do j=jts,jte
        do k=kts,kte
          do i=its,ite
            rthindten(i,k,j) = 0.
          enddo
        enddo
      enddo

      ! We add all heat to the first atmospheric level.  
      kk = 1
      
      do j = jts,jte
          do i = its,ite
              lu = INT(grid%LU_INDEX(i,j)+0.1) ! Land use type  
              rho_i = 1./rho(i,kk,j) ! The inverse of the density
              mut = grid%mut(i,j)    ! The total atmospheric column mass
              if (lu .eq. lu_idx) then
                 rthindten(i,kk,j) = mut * cp_i * rho_i * heat / dz8w(i,kk,j)
                 ! We check to make sure that the calculation is not undefined.
                 ! If it is, we set the value to zero.  
                 if (rthindten(i,kk,j) .ne. rthindten(i,kk,j)) then
                     rthindten(i,kk,j) = 0
                 endif
              endif
          enddo
      enddo

    enddo

end subroutine ind_heat
end module module_ind_heat

! End subroutine by dmbrown1.