EB Subroutine

  subroutine EB
    !calculates the energy balance for every facet
    use modglobal, only: nfcts, boltz, tEB, AM, BM,CM,DM,EM,FM,GM,HM, inAM, bb,w, dumv,Tdash, timee, dtEB, tnextEB, rk3step, rhoa, cp, lEB, ntrun, lwriteEBfiles,nfaclyrs
    use initfac, only: faclam, faccp, netsw, facem, fachf, facef, fachfi, facT, facLWin, faca,facefi,facf,facets,facTdash,facqsat,facwsoil,facf,fachurel,facd,fackappa
    use modmpi, only: myid, comm3d, mpierr, MY_REAL, nprocs, cmyid
    use modstat_nc, only : writestat_nc, writestat_1D_nc, writestat_2D_nc
    real  :: ca = 0., cb = 0., cc = 0., cd = 0., ce = 0., cf = 0.
    real  :: ab = 0.
    integer :: l, n, m,i,j
    character(19) name

    if (.not. (lEB)) return
    !calculate latent heat flux from vegetation and soil
    call intqH
    !calculate energy balance, update facet temperature and soil moisture
    if ((rk3step .eq. 3) .and. (timee .ge. tnextEB)) then

      if (myid .eq. 0) then
        tEB = timee - tEB !time since last calculation of energy balance
        !write (*, *) "doing EB, time since last EB:", tEB

        !calculate time mean, facet area mean latent heat flux and update green roof
        !ILS13 02.05.18 ABOUT updateGR: convert latent heatflux E properly should be done before temperature calculatation. BUT the rest of updateGR should be done after!
        !update green roof
        call updateGR

        !get longwave fluxes for all facets
        call calclw

        !get time mean, facet area mean sensible heat flux
        do n = 1, nfcts
          fachfi(n) = fachfi(n)/tEB/faca(n)*rhoa*cp !mean heat flux since last EB calculation (time average)
          !since fachf is the sum over all cells making up a facet we need to divide by the number of cells, assuming a given density to convert to W/m2
        end do

        !solve the system:
        !see Suter 2018
        !A * T'= bb + B * T,   where T' = dT/dz
        !C * d/dtT + D d/dtT'= e * T'
        !
        !-> T(n+1)=(F-G*dt)^-1*(F*T+w*dt)
        !where F=(C + D*A^-1*B), G=(E*A^-1*B), w=(E*A^-1*bb)


        do n = 1, nfcts
          if (facets(n) < -100) cycle

          !calculate wallflux and update surface temperature
          !! define time dependent fluxes
          ab = boltz*facem(n)*(facT(n, 1)**3)/faclam(n, 1) ! ab*T is the Stefan-Boltzman law
          bb(1) = -(netsw(n) + facLWin(n) + fachfi(n) + facefi(n))/faclam(n, 1) !net surface flux

          !!define the matrices to solve wall heat flux
          !! CREATE MATRICES BASED ON WALL PROPERTIES
          i=1;m=0; !position along columns, placeholder for layerindex since only 3 layers implemented (initfac.f90)
          do j=1,nfaclyrs
            m=j  !!CARE!!! ONLY 3 LAYERS ARE CURRENTLY BEING READ FROM INPUT FILES. PROPERTIES OF LAYER 3 ARE USED FOR SUBSEQUENT LAYERS!!!
            ca=1./facd(n,m)
            BM(j+1,i)=-ca
            BM(j+1,i+1)=ca
            EM(j,i)=-faclam(n,m)
            EM(j,i+1)=faclam(n,m+1)
            cb=faccp(n,m)*facd(n,m)/2.
            CM(j,i)=cb
            CM(j,i+1)=cb
            ca=faccp(n,m)*facd(n,m)**2/12.
            DM(j,i)=ca
            DM(j,i+1)=-ca
            i=i+1
          end do
          CM(nfaclyrs+1,nfaclyrs+1)=1.
          BM(1,1)=ab


          w = matmul(EM, matmul(inAM,bb))*tEB !easier than loop and sum
          HM = matmul(inAM,BM)
          FM = CM + matmul(DM,HM)
          GM = matmul(EM,HM)
          HM = FM-GM*tEB
          if (nfaclyrs == 3) then
            GM = matinv4(HM)
          else
            GM = gaussji(HM,IDM,nfaclyrs+1)
          end if
          !instead of inverting matrix HM and multiplying by GM (=HM^-1) it would be waster to do a  left matrix division HM\x is faster than (HM^-1)*x
          dumv = matmul(GM, (matmul(FM,facT(n,:))+w))

          facT(n, :) = dumv
          !calculate Temperature gradient dT/dz=>Tdash so we can output it
          !ground heat flux = lambda dT/dz
          w = matmul(BM, dumv)
          facTdash(n, :) = matmul(inAM, (bb + w))

          !end if
        end do

        if (lwriteEBfiles) then
          if (myid == 0) then
            allocate(varsT(nfcts,nfaclyrs+1,nstatT))
            varsT(:,:,1) = facT(1:nfcts,1:nfaclyrs+1)
            varsT(:,:,2) = facTdash(1:nfcts,1:nfaclyrs+1)
            call writestat_nc(ncidT,1,tncstatT,(/timee/),nrecT,.true.)
            call writestat_2D_nc(ncidT,nstatT,ncstatT,varsT,nrecT,nfcts,nfaclyrs+1)
            deallocate(varsT)

            allocate(varsEB(nfcts,nstatEB))
            varsEB(:,1) = netsw(1:nfcts)
            varsEB(:,2) = facLWin(1:nfcts)
            varsEB(:,3) = boltz*facem(1:nfcts)*facT(1:nfcts,1)**4
            varsEB(:,4) = fachfi(1:nfcts)
            varsEB(:,5) = facefi(1:nfcts)
            varsEB(:,6) = facwsoil(1:nfcts)
            ! add longwave out
            call writestat_nc(ncidEB,1,tncstatEB,(/timee/),nrecEB,.true.)
            call writestat_1D_nc(ncidEB,nstatEB,ncstatEB,varsEB,nrecEB,nfcts)
            deallocate(varsEB)

          end if !myid
        end if

        tEB = timee !set time of last calculation of energy balance to current time
        tnextEB = NINT((timee + dtEB))*1.0  !rounded to nearest integer  (e.g. if current time is 10.013s and dtEb=10s, then the next energy balance will be calculated at t>=20s)
        !write (*, *) "time, time next EB", timee, tnextEB

        do n = 1, nfcts
          fachfi(n) = 0.
          facefi(n) = 0.
        end do
      end if !myid==0

      !write (*, *) "bcasting facT"
      call MPI_BCAST(facT(0:nfcts, 1:nfaclyrs+1), (nfaclyrs+1)*(nfcts + 1), MY_REAL, 0, comm3d, mpierr)
      call MPI_BCAST(tnextEB, 1, MY_REAL, 0, comm3d, mpierr)
      call MPI_BCAST(facqsat(0:nfcts), nfcts + 1, MY_REAL, 0, comm3d, mpierr)
      call MPI_BCAST(facf(0:nfcts, 1:5), (nfcts + 1)*5, MY_REAL, 0, comm3d, mpierr)
      call MPI_BCAST(fachurel(0:nfcts), nfcts + 1, MY_REAL, 0, comm3d, mpierr)
      !call MPI_BCAST(facwsoil(0:nfcts), nfcts + 1, MY_REAL, 0, comm3d, mpierr)
    end if !time>tnextEB

  end subroutine EB