initEB Subroutine

  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,nfaclyrs
    use initfac, only:facd, faccp, faclam, fackappa, netsw, facem, fachf, facef, fachfi, facT, facLWin,facefi,facwsoil,facf,facets,facTdash,facqsat,facf,fachurel
    use modmpi, only:myid, comm3d, mpierr, MY_REAL, nprocs, cmyid
    use modstat_nc,only: open_nc, define_nc,ncinfo,writestat_dims_nc
    integer :: i,j,k,l,m,n
    real :: dum

    if (.not. lEB) return

    allocate(AM(1:nfaclyrs+1,1:nfaclyrs+1))
    allocate(inAM(1:nfaclyrs+1,1:nfaclyrs+1))
    allocate(CM(1:nfaclyrs+1,1:nfaclyrs+1))
    allocate(bb(1:nfaclyrs+1))
    allocate(BM(1:nfaclyrs+1,1:nfaclyrs+1))
    allocate(DM(1:nfaclyrs+1,1:nfaclyrs+1))
    allocate(EM(1:nfaclyrs+1,1:nfaclyrs+1))
    allocate(FM(1:nfaclyrs+1,1:nfaclyrs+1))
    allocate(GM(1:nfaclyrs+1,1:nfaclyrs+1))
    allocate(HM(1:nfaclyrs+1,1:nfaclyrs+1))
    allocate(IDM(1:nfaclyrs+1,1:nfaclyrs+1))
    allocate(w(1:nfaclyrs+1))
    allocate(dumv(1:nfaclyrs+1))
    allocate(Tdash(1:nfaclyrs+1))

    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,nfaclyrs+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,nfaclyrs+1
      AM(j,m)=0.5
      AM(j,m+1)=0.5
      m=m+1
    end do
    AM(1,1)=1.0
    if (nfaclyrs == 3) then
      inAM = matinv4(AM)
      !!alternatively
      !inAM=matinv3(AM)
      !!or
    else
      inAM=gaussji(AM,IDM,nfaclyrs+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', 'Net shortwave', 'W/m^2','ft')
      call ncinfo(ncstatEB( 2,:),'LWin', 'Incoming longwave', 'W/m^2','ft')
      call ncinfo(ncstatEB( 3,:),'LWout', 'Outgoing longwave', 'W/m^2','ft')
      call ncinfo(ncstatEB( 4,:),'hf', 'Sensible heat', 'W/m^2','ft')
      call ncinfo(ncstatEB( 5,:),'ef', 'Latent heat', 'W/m^2','ft')
      call ncinfo(ncstatEB( 6,:),'WGR','Water content', '?','ft')


      if (myid==0) then
        call open_nc(Tname, ncidT, nrecT, nfcts=nfcts, nlyrs=nfaclyrs+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 (nrecEB==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