      subroutine fcn(n,fi,f,xl,nelmax)
c
c
      use variaveis
c
c     the system of equations
c      
      implicit real *8 (a-h,o-z)
      dimension fi(nelmax),f(nelmax),
     1	xl(nelmax,8,3)
      common /prop/ rhoair,gban,cvair,rkair,alfa,pr,cpair,rkwall
      common /geom/ thick
      common /temp/ tpoute,tpoutw,tpouts,tpoutn,tpoutt,
     1     tground
      common /hwl/ a,b,c
      common /rad/ tpair,alf,epsi,sig,riavg,rhoal,esp,cal
      common /rad1/ iwe,iww,iwn,iws
      common /sol2/ rkiso(5) 
      common /mesh/ nel,nlx,nly,nlz,neleqp
      common /velair/ fanvelo,uhp
c
c	loop to create the temperature equation from
c	a balance of energy for each volume element,
c	either for air or
c	for equipment. Also computes the internal wall 
c	temperatures when necessary
c
c     ieast, seast, etc - counters and summers to compute 
c     average temperatures of the inner surfaces
c
      ieast=0
      seast=0.d0
      iwest=0
      swest=0.d0
      isouth=0
      ssouth=0.d0
      inorth=0
      snorth=0.d0
      itop=0
      sutop=0.d0
c
c     BEGINNING OF THE LOOP OVER ALL ELEMENTS
c
      do i=1,n-5
c
c     mass of air or eqp in the volume element and element geometry
c
	  wid=xl(i,4,2)-xl(i,1,2)
	  hgt=xl(i,5,3)-xl(i,1,3)
	  rlg=xl(i,2,1)-xl(i,1,1)
	  vol=wid*hgt*rlg  ! air or eqp vol
          if (itypel(i).eq.0) then
             rmel=rhoair*vol
             cvel=cvair
             else
                rmel=rhoeqp(i)*vol
                cvel=ceqp(i)
          endif
          areaew=wid*hgt
          areans=hgt*rlg
          areatb=wid*rlg
c
c	temperature derivative f(i) = dT/dtime
c	where fi(i) = Tair or Teqp
c
c	computing each heat transfer interaction based on a given
c	geometry (vector xl and gabinet size)
c
c	East face of element
c
	if (xl(i,1,1).eq.0.d0) then
c	
c       the element could have air or eqp in it and
c	the east face is in contact with the east wall of the gabinet
c       then, natural convection on a vertical wall, if air and
c       conduction if eqp.
c       H/L >> 1, then conduction between the 2 walls is dominant
c       (tall enclosure limit)
c
          ieast=ieast+1
          if(itypel(i).eq.0) then
             c2=rkiso(1)/thick
c             c1=hint
             c1=hnc(rkair,hgt,rlg,wid,gban,qgeqp(i),
     1            fi(i),fi(n-4))
             ce=c2/c1
             tpin=(ce*fi(n-4)+fi(i))/(1.d0+ce)
             qe=-rkiso(1)*areaew*(tpin-fi(n-4))/thick
             seast=seast+tpin
          else
             rlgm=rlg/2
             ue=1.d0/(rlgm/rkeqp(i)+esp/rkwall+thick/rkiso(1))
             qe=-ue*areaew*(fi(i)-fi(n-4))
             tpin=fi(n-4)-qe*thick/rkiso(1)/areaew
             seast=seast+tpin
          endif
          else
c
c      the element could have air or eqp inside, but in both cases
c       conduction only between adjacent elements, through practically
c       no relative motion between air layers or eqp/eqp, or
c     convection when air/eqp
c
             rlgm=rlg/2
             rlgam=(xl(i-1,2,1)-xl(i-1,1,1))/2
             if(itypel(i).eq.0.and.itypel(i-1).eq.0) then
                ue=rkair/(rlgm+rlgam)
                qe=-ue*areaew*(fi(i)-fi(i-1))
                else
                   if(itypel(i).eq.1.and.itypel(i-1).eq.1) then
                      rki=rkeqp(i)
                      rkia=rkeqp(i-1)
                      ue=1.d0/(rlgm/rki+rlgam/rkia)
                      qe=-ue*areaew*(fi(i)-fi(i-1))
                   else
                   if(itypel(i).ne.0.and.itypel(i-1).eq.0) then
                      rg=rlgm/rkeqp(i)
                      if(fanvelo.gt.0.d0) then
c
c     compute forced convection heat transfer coefficient
c
                      hi=hfcwall(hgt,rlg,wid)
                      else
c
c     compute natural convection heat transfer coefficient
c
                   hi=hnc(rkair,hgt,rlg,wid,gban,qgeqp(i),
     1                  fi(i),fi(i-1))
                      endif
                   else
                      rg=rlgam/rkeqp(i-1)
                      if(fanvelo.gt.0.d0) then
c
c     compute forced convection heat transfer coefficient
c
                      hi=hfcwall(hgt,rlg,wid)
                      else
c
c     compute natural convection heat transfer coefficient
c
                      hgtia=xl(i-1,5,3)-xl(i-1,1,3)
                      rlgia=xl(i-1,2,1)-xl(i-1,1,1)
                      widia=xl(i-1,4,2)-xl(i-1,1,2)
                   hi=hnc(rkair,hgtia,rlgia,widia,gban,qgeqp(i-1),
     1                  fi(i),fi(i-1))
                      endif  ! end of fanvelo if
                   endif
                   ue=1.d0/(1.d0/hi+rg)
                   qe=ue*areaew*(fi(i-1)-fi(i))
                   endif    
             endif
    	endif
c
c       West face of element
c
        if (xl(i,2,1).eq.a) then
c
c       the element could have air or eqp in it and
c       the west face is in contact with the west wall of the gabinet
c       then, natural convection on a vertical wall, if air and
c       conduction if eqp.
c       H/L >> 1, then conduction between the 2 walls is dominant
c       (tall enclosure limit)
c
          iwest=iwest+1
          if(itypel(i).eq.0) then
             c2=rkiso(2)/thick
c             c1=hint
             c1=hnc(rkair,hgt,rlg,wid,gban,qgeqp(i),
     1            fi(i),fi(n-3))
             cw=c2/c1
             tpin=(cw*fi(n-3)+fi(i))/(1.d0+cw)
             qw=-rkiso(2)*areaew*(tpin-fi(n-3))/thick
             swest=swest+tpin
          else
             rlgm=rlg/2
             uw=1.d0/(rlgm/rkeqp(i)+esp/rkwall+thick/rkiso(2))
             qw=-uw*areaew*(fi(i)-fi(n-3))
             tpin=fi(n-3)-qw*thick/rkiso(2)/areaew
             swest=swest+tpin
          endif
          else
c
c      the element could have air or eqp inside, but in both cases
c       conduction only between adjacent elements, through practically
c       no relative motion between air layers or eqp/eqp, or
c     convection when air/eqp
c
             rlgm=rlg/2
             rlgpm=(xl(i+1,2,1)-xl(i+1,1,1))/2
             if(itypel(i).eq.0.and.itypel(i+1).eq.0) then
                uw=rkair/(rlgm+rlgpm)
                qw=-uw*areaew*(fi(i)-fi(i+1))
                else
                   if(itypel(i).eq.1.and.itypel(i+1).eq.1) then
                      rki=rkeqp(i)
                      rkip=rkeqp(i+1)
                      uw=1.d0/(rlgm/rki+rlgpm/rkip)
                      qw=-uw*areaew*(fi(i)-fi(i+1))
                   else
                   if(itypel(i).ne.0.and.itypel(i+1).eq.0) then
                      rg=rlgm/rkeqp(i)
                      if(fanvelo.gt.0.d0) then
c
c     compute forced convection heat transfer coefficient
c
                      hi=hfcwall(hgt,rlg,wid)
                      else
c
c     compute natural convection heat transfer coefficient
c
                   hi=hnc(rkair,hgt,rlg,wid,gban,qgeqp(i),
     1                  fi(i),fi(i+1))
                      endif
                   else
                      rg=rlgpm/rkeqp(i+1)
                      if(fanvelo.gt.0.d0) then
c
c     compute forced convection heat transfer coefficient
c
                      hi=hfcwall(hgt,rlg,wid)
                      else
c
c     compute natural convection heat transfer coefficient
c
                      hgtip=xl(i+1,5,3)-xl(i+1,1,3)
                      rlgip=xl(i+1,2,1)-xl(i+1,1,1)
                      widip=xl(i+1,4,2)-xl(i+1,1,2)
                hi=hnc(rkair,hgtip,rlgip,widip,gban,qgeqp(i+1),
     1                  fi(i),fi(i+1))
                      endif
                   endif
                   uw=1.d0/(1.d0/hi+rg)
                   qw=uw*areaew*(fi(i+1)-fi(i))
                   endif
             endif
    	endif
c
c     Top face of element
c
        if (xl(i,5,3).eq.c) then 
c
c     The element could have air or eqp in it and
c     The top face is in contact with the top wall of the gabinet
c
          itop=itop+1
          if(itypel(i).eq.0) then
             hgtm=hgt/2
             hi=hnc(rkair,hgt,rlg,wid,gban,qgeqp(i),
     1            fi(i),fi(n))
             ut=1.d0/(1.d0/hi+esp/rkwall+thick/rkiso(5))
             qt=-ut*areatb*(fi(i)-fi(n))
             tpin=fi(n)-qt*thick/rkiso(5)/areatb
             sutop=sutop+tpin
          else
             hgtm=hgt/2
             ut=1.d0/(hgtm/rkeqp(i)+esp/rkwall+thick/rkiso(5))
             qt=-ut*areatb*(fi(i)-fi(n))
             tpin=fi(n)-qt*thick/rkiso(5)/areatb
             sutop=sutop+tpin
          endif
           else
c
c     The element might have air or eqp inside.
c     The top face is in contact with another volume element.
c     We have to analize 4 possibilities air/air, airabove/eqpbelow,
c     airbelow/eqpabove and eqp/eqp.
c     First, we check if the element uses a heat pipe or not
c     If there is a heat pipe, then we compute the appropriate 
c     vaporization temperature and the heat rejected to the heat pipe
c     by the heat source
c
              if (itypel(i).eq.1.and.qgeqp(i).ne.0.d0
     *             .and.ifhp(i).ne.0) then
		 ihx=ifhp(i)
		 tvap=(uhp*atc(ihx)*fi(ihx)+
     *	rkeqp(i)*areatb*fi(i))/(uhp*atc(ihx)+rkeqp(i)*areatb)
		qgeqp(ihx)=uhp*atc(ihx)*(tvap-fi(ihx))
		qt=-qgeqp(ihx)
              else
c
c     Next, we identify the element number over the current element
c
              it=0
              do j=1,n-5
                 if (xl(i,5,1).eq.xl(j,1,1).and.
     1			xl(i,5,2).eq.xl(j,1,2).and.
     1			xl(i,5,3).eq.xl(j,1,3)) then
                    it=j
                    goto 322
                 endif
              enddo
 322          continue
              delt=fi(it)-fi(i)
             hgtm=hgt/2
             hgttm=(xl(it,5,3)-xl(it,1,3))/2
             if(itypel(i).ne.1.and.itypel(it).eq.0) then
                if (fanvelo.gt.0.d0) then
                   velo=fanvelo
                else
                   velo=alfa*(gban*abs(delt)*c)**0.5d0
                endif
                rmdot=rhoair*velo*areatb/2
                qt=rmdot*cpair*delt
                else
                   if(itypel(i).eq.1.and.itypel(it).eq.0) then
                      hi=hnc(rkair,hgt,rlg,wid,gban,qgeqp(i),
     1                     fi(i),fi(it))
                         qt=1.d0/(1.d0/hi+hgtm/rkeqp(i))
     1                     *areatb*(fi(it)-fi(i))
                      else
                      if(itypel(i).eq.0.and.itypel(it).eq.1) then
                         rkit=rkeqp(it)
                      hgtit=xl(it,5,3)-xl(it,1,3)
                      rlgit=xl(it,2,1)-xl(it,1,1)
                      widit=xl(it,4,2)-xl(it,1,2)
                   hi=hnc(rkair,hgtit,rlgit,widit,gban,qgeqp(it),
     1                  fi(i),fi(it))
                         ut=1.d0/(1.d0/hi+hgttm/rkit)
                         qt=ut*areatb*delt
                         else
                         rki=rkeqp(i)
                         rkit=rkeqp(it)
                         ut=1.d0/(hgtm/rki+hgttm/rkit)
                         qt=ut*areatb*delt
                      endif    
                   endif
             endif
          endif
    	endif
c
c     Bottom face of element
c
        if (xl(i,1,3).eq.0.d0) then 
c
c     The element could have air or eqpt in it and
c     The bottom face is in contact with the bottom wall of the gabinet
c     Pure conduction if the lower temperature is
c     smaller than the upper temperature, else heat transfer is by
c     Benard convection.
c
           if(itypel(i).eq.0) then
c              tu=2*fi(i)-tground
              delt=tground-fi(i)
              hi=hnc(rkair,hgt,rlg,wid,gban,qgeqp(i),
     1             fi(i),tground)              
              qb=hi*areatb*delt
           else
              hgtm=hgt/2
              qb=-rkeqp(i)*areatb*(fi(i)-tground)/hgtm
           endif
          else
c
c     The element might have air or eqp inside.
c     The bottom face is in contact with another volume element.
c     We have to analize 4 possibilities air/air, airabove/eqpbelow,
c     airbelow/eqpabove and eqp/eqp.
c     First, we identify the element number below the current element
c
              ib=0
              do j=1,n-5
                 if (xl(i,1,1).eq.xl(j,5,1).and.
     1			xl(i,1,2).eq.xl(j,5,2).and.
     1			xl(i,1,3).eq.xl(j,5,3)) then
                    ib=j
                    goto 325
                 endif
              enddo
 325          continue
              delt=fi(ib)-fi(i)
             hgtm=hgt/2
             hgtbm=(xl(ib,5,3)-xl(ib,1,3))/2
             if(itypel(i).ne.1.and.itypel(ib).eq.0) then
                if(fanvelo.gt.0.d0) then
                   velo=fanvelo
                else
                   velo=alfa*(gban*abs(delt)*c)**0.5d0
                endif
                rmdot=rhoair*velo*areatb/2
                qb=rmdot*cpair*delt
                else
                   if(itypel(ib).eq.1.and.itypel(i).eq.0) then
                      hgtib=xl(ib,5,3)-xl(ib,1,3)
                      rlgib=xl(ib,2,1)-xl(ib,1,1)
                      widib=xl(ib,4,2)-xl(ib,1,2)
                   hi=hnc(rkair,hgtib,rlgib,widib,gban,qgeqp(ib),
     1                  fi(i),fi(ib))
                      qb=1.d0/(1.d0/hi+hgtbm/rkeqp(ib))*areatb*
     1                     delt
                      else
                      if(itypel(ib).eq.0.and.itypel(i).eq.1) then
                         rki=rkeqp(i)
                   hi=hnc(rkair,hgt,rlg,wid,gban,qgeqp(i),
     1                  fi(i),fi(ib))
                         ub=1.d0/(1.d0/hi+hgtm/rki)
                         qb=ub*areatb*delt
                         else
                         rki=rkeqp(i)
                         rkib=rkeqp(ib)
                         ub=1.d0/(hgtm/rki+hgtbm/rkib)
                         qb=ub*areatb*delt
                      endif    
                   endif
             endif
    	endif
c
c	North face of element
c
	if (xl(i,1,2).eq.0.d0) then
c	
c       the element could have air or eqp in it and
c	the north face is in contact with the north wall of the gabinet
c       then, natural convection on a vertical wall, if air and
c       conduction if eqp.
c       H/L >> 1, then conduction between the 2 walls is dominant
c       (tall enclosure limit)
c
          inorth=inorth+1
          if(itypel(i).eq.0) then
             c2=rkiso(4)/thick
c             c1=hint
             c1=hnc(rkair,hgt,rlg,wid,gban,qgeqp(i),
     1            fi(i),fi(n-1))
             cn=c2/c1
             tpin=(cn*fi(n-1)+fi(i))/(1.d0+cn)
             qn=-rkiso(4)*areans*(tpin-fi(n-1))/thick
             snorth=snorth+tpin
          else
             widm=wid/2
             un=1.d0/(widm/rkeqp(i)+esp/rkwall+thick/rkiso(4))
             qn=-un*areans*(fi(i)-fi(n-1))
             tpin=fi(n-1)-qn*thick/rkiso(4)/areans
             snorth=snorth+tpin
          endif
          else
c
c      the element could have air or eqp inside, but in both cases
c       conduction only between adjacent elements, through practically
c       no relative motion between air layers or eqp/eqp, or
c     convection when air/eqp.
c	Determine which is the adjacent element
c
              itn=0
              do j=1,n-5
                 if (xl(i,1,1).eq.xl(j,4,1).and.
     1			xl(i,1,2).eq.xl(j,4,2).and.
     1			xl(i,1,3).eq.xl(j,4,3)) then
                    itn=j
                    goto 422
                 endif
              enddo
 422          continue
             widm=wid/2
             widnm=(xl(itn,4,2)-xl(itn,1,2))/2
             if(itypel(i).eq.0.and.itypel(itn).eq.0) then
                un=rkair/(widm+widnm)
                qn=-un*areans*(fi(i)-fi(itn))
                else
                   if(itypel(i).eq.1.and.itypel(itn).eq.1) then
                      rki=rkeqp(i)
                      rkin=rkeqp(itn)
                      un=1.d0/(widm/rki+widnm/rkin)
                      qn=-un*areans*(fi(i)-fi(itn))
                   else
c************************************************************
                   if(itypel(i).ne.0.and.itypel(itn).eq.0) then
                      wg=widm/rkeqp(i)
                      if(fanvelo.gt.0.d0) then
c
c     compute forced convection heat transfer coefficient
c
                      hi=hfcwall(hgt,rlg,wid)
                      else
c
c     compute natural convection heat transfer coefficient
c
                   hi=hnc(rkair,hgt,rlg,wid,gban,qgeqp(i),
     1                  fi(i),fi(itn))
                      endif
!c**********************************************************
                   else
                      wg=widnm/rkeqp(itn)
                      if(fanvelo.gt.0.d0) then
!c
!c     compute forced convection heat transfer coefficient
!c
                      hi=hfcwall(hgt,rlg,wid)
                      else
!c
!c     compute natural convection heat transfer coefficient
!c
                      hgtitn=xl(itn,5,3)-xl(itn,1,3)
                      rlgitn=xl(itn,2,1)-xl(itn,1,1)
                      widitn=xl(itn,4,2)-xl(itn,1,2)
             hi=hnc(rkair,hgtitn,rlgitn,widitn,gban,qgeqp(itn),
     1                  fi(i),fi(itn))
                      endif
                   endif
                   un=1.d0/(1.d0/hi+wg)
                   qn=un*areans*(fi(itn)-fi(i))
                   endif
             endif
    	endif
!c
!c	South face of element
!c
	if (xl(i,4,2).eq.b) then
!c	
!c       the element must have air in it and
!c	the south face is in contact with the south wall of the gabinet
!c       then, natural convection on a vertical wall, if air and
!c       conduction if eqp.
!c       H/L >> 1, then conduction between the 2 walls is dominant
!c       (tall enclosure limit)
!c
          isouth=isouth+1
          if(itypel(i).eq.0) then
             c2=rkiso(3)/thick
!c             c1=hint
             c1=hnc(rkair,hgt,rlg,wid,gban,qgeqp(i),
     1            fi(i),fi(n-2))
             cs=c2/c1
             tpin=(cs*fi(n-2)+fi(i))/(1.d0+cs)
             qs=-rkiso(3)*areans*(tpin-fi(n-2))/thick
             ssouth=ssouth+tpin
          else
             widm=wid/2
             us=1.d0/(widm/rkeqp(i)+esp/rkwall+thick/rkiso(3))
             qs=-us*areans*(fi(i)-fi(n-2))
             tpin=fi(n-2)-qs*thick/rkiso(3)/areans
             ssouth=ssouth+tpin
          endif
          else
!c
!c      the element could have air or eqp inside, but in both cases
!c       conduction only between adjacent elements, through practically
!c       no relative motion between air layers or eqp/eqp, or
!c     convection when air/eqp.
!c	Determine which is the adjacent element
!c
              its=0
              do j=1,n-5
                 if (xl(i,4,1).eq.xl(j,1,1).and.
     1			xl(i,4,2).eq.xl(j,1,2).and.
     1			xl(i,4,3).eq.xl(j,1,3)) then
                    its=j
                    goto 423
                 endif
              enddo
 423          continue
             widm=wid/2
             widsm=(xl(its,4,2)-xl(its,1,2))/2
             if(itypel(i).eq.0.and.itypel(its).eq.0) then
                us=rkair/(widm+widsm)
                qs=-us*areans*(fi(i)-fi(its))
                else
                   if(itypel(i).eq.1.and.itypel(its).eq.1) then
                      rki=rkeqp(i)
                      rkis=rkeqp(its)
                      us=1.d0/(widm/rki+widsm/rkis)
                      qs=-us*areans*(fi(i)-fi(its))
                   else
!c************************************************************
                   if(itypel(i).ne.0.and.itypel(its).eq.0) then
                      wg=widm/rkeqp(i)
                      if(fanvelo.gt.0.d0) then
!c
!c     compute forced convection heat transfer coefficient
!c
                      hi=hfcwall(hgt,rlg,wid)
                      else
!c
!c     compute natural convection heat transfer coefficient
!c
                   hi=hnc(rkair,hgt,rlg,wid,gban,qgeqp(i),
     1                  fi(i),fi(its))
                      endif
!c**********************************************************
                   else
                      wg=widsm/rkeqp(its)
                      if(fanvelo.gt.0.d0) then
!c
!c     compute forced convection heat transfer coefficient
!c
                      hi=hfcwall(hgt,rlg,wid)
                      else
!c
!c     compute natural convection heat transfer coefficient
!c
                      hgtits=xl(its,5,3)-xl(its,1,3)
                      rlgits=xl(its,2,1)-xl(its,1,1)
                      widits=xl(its,4,2)-xl(its,1,2)
            hi=hnc(rkair,hgtits,rlgits,widits,gban,qgeqp(its),
     1                  fi(i),fi(its))
                      endif
                   endif
                   us=1.d0/(1.d0/hi+wg)
                   qs=us*areans*(fi(its)-fi(i))
                   endif
             endif
    	endif
!c*****************************
!c
      	f(i)=(1.d0/rmel/cvel)*(qe+qn+qw+qs+qt+qb+qgeqp(i))
!c
      enddo
!c
!c     now compute time derivatives of external wall
!c     temperatures from radiation heat transfer
!c
      alri=alf*riavg  
!c
!c     this is only for the top and 2 side walls. The other 2 walls do
!c     not receive sun radiosity at the same time
!c
      epsig=epsi*sig
      tpa=tpair-6.d0
      rec=rhoal*esp*cal
!c
!c     compute heat transfer coefficient outside the gabinet
!c     in the absence of winds (natural convection) - Lienhardt correlation
!c     worst possible situation (critical condition), for each wall
!c     nevertheless, the average hair = 10 - 15 W/(m2.K) 
!c     (Heat Transfer, Bejan, 1993)
!c
!c     east wall
!c
      hair=hnc(rkair,c,a,b,gban,0.d0,
     1     fi(n-4),tpair)
      teavg=seast/ieast
      if(iwe.eq.0) then
         ale=0.d0
         else
            ale=alri
      endif
      f(n-4)=(ale+epsig*(tpa**4-fi(n-4)**4)+hair*(tpair-fi(n-4))
     1     -rkiso(1)*(fi(n-4)-teavg)/thick)/rec
!c
!c     west wall
!c
      hair=hnc(rkair,c,a,b,gban,0.d0,
     1     fi(n-3),tpair)
      twavg=swest/iwest
      if(iww.eq.0) then
         alw=0.d0
         else
            alw=alri
      endif
      f(n-3)=(alw+epsig*(tpa**4-fi(n-3)**4)+hair*(tpair-fi(n-3))
     1     -rkiso(2)*(fi(n-3)-twavg)/thick)/rec
!c
!c     south wall
!c
      hair=hnc(rkair,c,a,b,gban,0.d0,
     1     fi(n-2),tpair)
      tsavg=ssouth/isouth
      if(iws.eq.0) then
         als=0.d0
         else
            als=alri
      endif
      f(n-2)=(als+epsig*(tpa**4-fi(n-2)**4)+hair*(tpair-fi(n-2))
     1     -rkiso(3)*(fi(n-2)-tsavg)/thick)/rec
!c
!c     north wall
!c
      hair=hnc(rkair,c,a,b,gban,0.d0,
     1     fi(n-1),tpair)
      tnavg=snorth/inorth
      if(iwn.eq.0) then
         aln=0.d0
         else
            aln=alri
      endif
      f(n-1)=(aln+epsig*(tpa**4-fi(n-1)**4)+hair*(tpair-fi(n-1))
     1     -rkiso(4)*(fi(n-1)-tnavg)/thick)/rec
!c
!c     top wall
!c
!c
!c     heat transfer coefficient computed here from Benard convection
!c     Note that the height between the two driving temperatures (tpair
!c     and fi(n)) cancels out, since Ra**(1/3)
!c
      hair=hnc(rkair,c,a,b,gban,0.d0,
     1     fi(n),tpair)
      ttavg=sutop/itop
      f(n)=(alri+epsig*(tpa**4-fi(n)**4)+hair*(tpair-fi(n))
     1     -rkiso(5)*(fi(n)-ttavg)/thick)/rec
!c
      return
      end
!c------------------------------------------------------------------
      function hnc(rkair,hgt,rlg,wid,gban,qgeqp,tel,tside)
      implicit real *8 (a-h,o-z)
      common /hwl/ a,b,c
      common /heat/ ct,cq
      if(qgeqp.eq.0.d0) then
         hnc=ct*rkair*(3.47d0+(gban*(hgt/3)**3*abs(tel-tside))
     1        **(1.d0/4.d0))/(hgt/3)
      else
         face=2*(rlg*hgt+wid*hgt+wid*rlg)
         hnc=cq*rkair*(gban*qgeqp/face/rkair)
     1        **0.25d0
      endif
      return
      end
!c------------------------------------------------------------------
      function hfcwall(hgt,rlg,wid)
      implicit real *8 (a-h,o-z)
      common /velair/ fanvelo,uhp
      common /prop/ rhoair,gban,cvair,rkair,alfa,pr,cpair,rkwall
      rnu=pr*alfa
      re=fanvelo*hgt/rnu
      hfcwall=rkair*0.037d0*pr**(1.d0/3.d0)*(re**0.8d0-23550.d0)/hgt
      return
      end
!c------------------------------------------------------------------
