subroutine initEB
!initialise everything necessary to calculate the energy balance
use modglobal, only:AM, BM,CM,DM,EM,FM,GM, HM, IDM, inAM, bb,w,dumv,Tdash, bldT, nfcts,nwalllayers
use initfac, only:facdi, faccp, faclami, fackappa, netsw, facem, fachf, facef, fachfi, facT, facLWin, facain,facefi,facwsoil,facf,facets,facTdash,facqsat,facf,fachurel
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: open_nc, define_nc,ncinfo,writestat_dims_nc
integer :: i,j,k,l,m,n
real :: dum
allocate(AM(1:nwalllayers+1,1:nwalllayers+1))
allocate(inAM(1:nwalllayers+1,1:nwalllayers+1))
allocate(CM(1:nwalllayers+1,1:nwalllayers+1))
allocate(bb(1:nwalllayers+1))
allocate(BM(1:nwalllayers+1,1:nwalllayers+1))
allocate(DM(1:nwalllayers+1,1:nwalllayers+1))
allocate(EM(1:nwalllayers+1,1:nwalllayers+1))
allocate(FM(1:nwalllayers+1,1:nwalllayers+1))
allocate(GM(1:nwalllayers+1,1:nwalllayers+1))
allocate(HM(1:nwalllayers+1,1:nwalllayers+1))
allocate(IDM(1:nwalllayers+1,1:nwalllayers+1))
allocate(w(1:nwalllayers+1))
allocate(dumv(1:nwalllayers+1))
allocate(Tdash(1:nwalllayers+1))
write(*,*) "nwalllayers",nwalllayers
BM=0.;DM=0.;EM=0.;FM=0.;GM=0.;HM=0.;w=0.;dumv=0.;Tdash=0.;
AM=0.;inAM=0.;CM=0.;IDM=0.;bb=0.
do j=1,nwalllayers+1
IDM(j,j)=1.0
end do
!Fortran is column major, i.e. left dimensions should be iterated first
! e.g. (1,1)->(2,1)->(3,1)->(1,2)->... since they are next to each other on memory
!first index moves "up and down" second "left and right" (as always)
m=1; !position along columns
do j=2,nwalllayers+1
AM(j,m)=0.5
AM(j,m+1)=0.5
m=m+1
end do
AM(1,1)=1.0
if (nwalllayers == 3) then
inAM = matinv4(AM)
!!alternatively
!inAM=matinv3(AM)
!!or
else
inAM=gaussji(AM,IDM,nwalllayers+1)
end if
! write facet temperatures to facT.xxx.nc, and energies to facEB.xxx.nc
if (lwriteEBfiles) then
Tname(6:8) = cexpnr
EBname(7:9) = cexpnr
allocate(ncstatT(nstatT,4))
call ncinfo(tncstatT(1,:),'t', 'Time', 's', 'time')
call ncinfo(ncstatT( 1,:),'T' ,'Temperature', 'K','flt')
call ncinfo(ncstatT( 2,:),'dTdz','Temperature gradient','K/m','flt' )
allocate(ncstatEB(nstatEB,4))
call ncinfo(tncstatEB(1,:),'t', 'Time', 's', 'time')
call ncinfo(ncstatEB( 1,:),'netsw', 'Shortwave radiation', 'W','ft')
call ncinfo(ncstatEB( 2,:),'LWin', 'Longwave radiation', 'W','ft')
call ncinfo(ncstatEB( 3,:),'hf', 'Sensible heat flux', 'W','ft')
call ncinfo(ncstatEB( 4,:),'ef', 'Latent heat flux', 'W','ft')
call ncinfo(ncstatEB( 5,:),'WGR','Moisture?', 'W','ft')
if (myid==0) then
call open_nc(Tname, ncidT, nrecT, nfcts=nfcts, nlyrs=nwalllayers+1)
call open_nc(EBname, ncidEB, nrecEB, nfcts=nfcts)
if (nrecT==0) then
call define_nc( ncidT, 1, tncstatT)
call writestat_dims_nc(ncidT)
end if
if (nrecT==0) then
call define_nc( ncidEB, 1, tncstatEB)
call writestat_dims_nc(ncidEB)
end if
call define_nc( ncidT, nstatT, ncstatT)
call define_nc( ncidEB, nstatEB, ncstatEB)
endif !myid==0
end if
end subroutine initEB