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, tnextEB, rk3step, rhoa, cp, lEB, ntrun, lwriteEBfiles,nwalllayers
use initfac, only: faclami, netsw, facem, fachf, facef, fachfi, facT, facLWin, facain,facefi,facf,facets,facTdash,facqsat,facwsoil,facf,fachurel,facd,facdi,fackappa
use modmpi, only: myid, comm3d, mpierr, MPI_INTEGER, MPI_DOUBLE_PRECISION, MY_REAL, nprocs, cmyid, MPI_REAL8, MPI_REAL4, MPI_SUM
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/facain(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, 2) < -100) then !it's a bounding wall, no reason to do energy balance
cycle
else
!calculate wallflux and update surface temperature
!! define time dependent fluxes
ab = faclami(n, 1)*boltz*facem(n)*(facT(n, 1)**3) ! ab*T is the Stefan-Boltzman law
bb(1) = -faclami(n, 1)*(netsw(n) + facLWin(n) + fachfi(n) + facefi(n)) !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,nwalllayers
m=min(j,3) !!CARE!!! ONLY 3 LAYERS ARE CURRENTLY BEING READ FROM INPUT FILES. PROPERTIES OF LAYER 3 ARE USED FOR SUBSEQUENT LAYERS!!!
ca=facdi(n,m)
BM(j+1,i)=-ca
BM(j+1,i+1)=ca
EM(j,i)=-fackappa(n,m)
EM(j,i+1)=fackappa(n,m+1)
cb=facd(n,m)/2
CM(j,i)=cb
CM(j,i+1)=cb
ca=cb**2*0.33333333
DM(j,i)=ca
DM(j,i+1)=-ca
i=i+1
end do
CM(nwalllayers+1,nwalllayers+1)=1.0
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 (nwalllayers == 3) then
GM = matinv4(HM)
else
GM = gaussji(HM,IDM,nwalllayers+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,nwalllayers+1,nstatT))
varsT(:,:,1) = facT(1:nfcts,1:nwalllayers+1)
varsT(:,:,2) = facTdash(1:nfcts,1:nwalllayers+1)
call writestat_nc(ncidT,1,tncstatT,(/timee/),nrecT,.true.)
call writestat_2D_nc(ncidT,nstatT,ncstatT,varsT,nrecT,nfcts,nwalllayers+1)
deallocate(varsT)
allocate(varsEB(nfcts,nstatEB))
varsEB(:,1) = netsw(1:nfcts)
varsEB(:,2) = facLWin(1:nfcts)
varsEB(:,3) = fachfi(1:nfcts)
varsEB(:,4) = facefi(1:nfcts)
varsEB(:,5) = facwsoil(1:nfcts)
call writestat_nc(ncidEB,1,tncstatEB,(/timee/),nrecEB,.true.)
call writestat_1D_nc(ncidEB,nstatEB,ncstatEB,varsEB,nrecEB,nfcts)
deallocate(varsEB)
endif !myid
! if (lwriteEBfiles) then
! name = 'tEB____________.txt'
! open (unit=11, file=name, position='append')
! write (11, '(1(F10.4,:,","))') tEB
! close (11)
! name = 'dummy__________.csv'
! write (name(6:15), '(F10.4)') timee
! write (name(1:5), '(A5)') 'netsw'
! open (unit=11, file=name, position='append')
! write (11, '(249(F10.4,:,","))') netsw(1:nfcts)
! close (11)
! write (name(1:5), '(A5)') 'LWin_'
! open (unit=11, file=name, position='append')
! write (11, '(249(F10.4,:,","))') facLWin(1:nfcts)
! close (11)
! write (name(1:5), '(A5)') 'hf___'
! open (unit=11, file=name, position='append')
! write (11, '(249(F10.4,:,","))') fachfi(1:nfcts)
! close (11)
! write (name(1:5), '(A5)') 'ef___'
! open (unit=11, file=name, position='append')
! write (11, '(249(F10.4,:,","))') facefi(1:nfcts)
! close (11)
! write (name(1:5), '(A5)') 'WGR__'
! open (unit=11, file=name, position='append')
! write (11, '(249(F10.4,:,","))') facwsoil(1:nfcts)
! close (11)
! write (name(1:5), '(A5)') 'facT1'
! open (unit=11, file=name, position='append')
! write (11, '(249(F10.4,:,","))') facT(1:nfcts, 1)
! close (11)
! write (name(1:5), '(A5)') 'facT2'
! open (unit=11, file=name, position='append')
! write (11, '(249(F10.4,:,","))') facT(1:nfcts, 2)
! close (11)
! write (name(1:5), '(A5)') 'facT3'
! open (unit=11, file=name, position='append')
! write (11, '(249(F10.4,:,","))') facT(1:nfcts, 3)
! close (11)
! write (name(1:5), '(A5)') 'faTd1'
! open (unit=11, file=name, position='append')
! write (11, '(249(F10.4,:,","))') facTdash(1:nfcts, 1)
! write (name(1:5), '(A5)') 'faTd2'
! open (unit=11, file=name, position='append')
! write (11, '(249(F10.4,:,","))') facTdash(1:nfcts, 2)
! write (name(1:5), '(A5)') 'faTd3'
! open (unit=11, file=name, position='append')
! write (11, '(249(F10.4,:,","))') facTdash(1:nfcts, 3)
! close (11)
! write (name(1:5), '(A5)') 'faTd4'
! open (unit=11, file=name, position='append')
! write (11, '(249(F10.4,:,","))') facTdash(1:nfcts, 4)
! close (11)
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:nwalllayers+1), (nwalllayers+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