/wrfv2_fire/phys/module_fr_sfire_core.F

Large files files are truncated, but you can click here to view the full file

!
!*** Jan Mandel August-October 2007 email: jmandel@ucar.edu or Jan.Mandel@gmail.com
!
! With contributions by Minjeong Kim.
#define DEBUG_OUT
#define DEBUG_PRINT
!#define DEBUG_OUT_FUEL_LEFT    

module module_fr_sfire_core

use module_fr_sfire_phys, only: fire_params , fire_ros
use module_fr_sfire_util

implicit none

! The mathematical core of the fire spread model. No physical constants here.
! 
! subroutine sfire_core: only this routine should be called from the outside.
! subroutine fuel_left:  compute remaining fuel from time of ignition.
! subroutine prop_ls: propagation of curve in normal direction.


contains

!
!****************************************
!
    
subroutine init_no_fire(&
    ifds,ifde,jfds,jfde, &
    ifms,ifme,jfms,jfme, &
    ifts,ifte,jfts,jfte, &
    fdx,fdy,time_start,dt,  & ! scalars in
    fuel_frac,fire_area,lfn,tign)    ! arrays out            
implicit none
             
!*** purpose: initialize model to no fire

!*** arguments
integer, intent(in):: ifds,ifde,jfds,jfde   ! fire domain bounds
integer, intent(in):: ifts,ifte,jfts,jfte   ! fire tile bounds
integer, intent(in):: ifms,ifme,jfms,jfme   ! array bounds
real, intent(in) :: fdx,fdy,time_start,dt   ! mesh spacing, time
real, intent(out), dimension (ifms:ifme,jfms:jfme) :: & 
                   fuel_frac,fire_area,lfn,tign       ! model state

!*** calls
intrinsic epsilon
                                                
!*** local
integer:: i,j
real:: lfn_init,time_init,time_now
character(len=128):: msg

time_now = time_start+dt
time_init = time_start+2*dt  ! could be time_start+dt if not for rounding errors
lfn_init = 2*max((ifde-ifds+1)*fdx,(jfde-jfds+1)*fdy)      ! more than domain diameter

if(fire_perimeter_time > 0.) then

  do j=jfts,jfte
     do i=ifts,ifte
        fuel_frac(i,j)=1.          ! fuel at start is 1 by definition
        fire_area(i,j)=0.          ! nothing burning
        lfn(i,j) = tign(i,j)-time_now       ! use specified ignition time as level set function 
        if(.not.lfn(i,j) > 0.)then
           write(msg,*)'i,j=',i,j,' tign=',tign(i,j),' <= time_now =',time_now
           call message(msg,level=0)
           call crash('init_no_fire: ignition time must be after the end of the time step')
        endif
     enddo
  enddo
  call message('init_no_fire: using given ignition time as history',level=1)
 
else

  do j=jfts,jfte
     do i=ifts,ifte
        fuel_frac(i,j)=1.          ! fuel at start is 1 by definition
        fire_area(i,j)=0.          ! nothing burning
        tign(i,j) = time_init      ! ignition in future
        lfn(i,j) = lfn_init        ! no fire 
     enddo
  enddo
  call message('init_no_fire: state set to no fire',level=1)
endif

call check_lfn_tign('init_no_fire',time_now,ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,lfn, tign)

end subroutine init_no_fire

!
!******************
!
 

subroutine ignite_fire( ifds,ifde,jfds,jfde,                    & ! fire domain dims - the whole domain
                        ifms,ifme,jfms,jfme,                      &
                        ifts,ifte,jfts,jfte,                      &
                        ignition_line,                            &
                        start_ts,end_ts,                    &
                        coord_xf,coord_yf,                &     
                        unit_xf,unit_yf,                  &
                        lfn,tign,ignited)
implicit none

!*** purpose: ignite a circular/line fire 

!*** description
! ignite fire in the region within radius r from the line (sx,sy) to (ex,ey).
! the coordinates of nodes are given as the arrays coord_xf and coord_yf
! r is given in m
! one unit of coord_xf is unit_xf m 
! one unit of coord_yf is unit_yf m 
! so a node (i,j) will be ignited iff for some (x,y) on the line
! || ( (coord_xf(i,j) - x)*unit_xf , (coord_yf(i,j) - y)*unit_yf ) || <= r 


!*** arguments
integer, intent(in):: ifds,ifde,jfds,jfde   ! fire domain bounds
integer, intent(in):: ifts,ifte,jfts,jfte   ! fire tile bounds
integer, intent(in):: ifms,ifme,jfms,jfme   ! array bounds
type(line_type), intent(in):: ignition_line    ! description of the ignition line
real, intent(in):: start_ts,end_ts          ! the time step start and end 
real, dimension(ifms:ifme, jfms:jfme), intent(in):: & 
    coord_xf,coord_yf                       !  node coordinates  
real, intent(in):: unit_xf,unit_yf          !  coordinate units in m
real, intent(inout), dimension (ifms:ifme,jfms:jfme) :: & 
                   lfn, tign                ! level function, ignition time (state)
integer, intent(out):: ignited              ! number of nodes newly ignited
                        
!*** local
integer:: i,j
real::lfn_new,tign_new,time_ign,ax,ay,rels,rele,d
! real:: lfn_new_chk,tign_new_chk,tmperr
real:: sx,sy                    ! start of ignition line, from lower left corner
real:: ex,ey                    ! end of ignition line, or zero
real:: st,et                    ! start and end of time of the ignition line
character(len=128):: msg
real::cx2,cy2,dmax,axmin,axmax,aymin,aymax,dmin
real:: start_x,start_y          ! start of ignition line, from lower left corner
real:: end_x,end_y              ! end of ignition line, or zero
real:: radius                   ! all within the radius of the line will ignite
real:: start_time,end_time      ! the ignition time for the start and the end of the line
real:: ros,tos                  ! ignition rate and time of spread
integer:: msglevel=4,smsg=2     ! print at this level
real:: lfn_min

!*** executable

call check_lfn_tign('ignite_fire start',end_ts,ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,lfn, tign)

! copy ignition line fields to local variables
start_x    = ignition_line%start_x ! x coordinate of the ignition line start point (m, or long/lat)
start_y    = ignition_line%start_y ! y coordinate of the ignition line start point
end_x      = ignition_line%end_x   ! x coordinate of the ignition line end point
end_y      = ignition_line%end_y   ! y coordinate of the ignition line end point
start_time = ignition_line%start_time ! ignition time for the start point from simulation start (s)
end_time   = ignition_line%end_time! ignition time for the end poin from simulation start (s)
radius     = ignition_line%radius  ! all within this radius ignites immediately
ros        = ignition_line%ros     ! rate of spread


tos        = radius/ros            ! time of spread to the given radius
st         = start_time            ! the start time of ignition considered in this time step
et         = min(end_ts,end_time)  ! the end time of the ignition segment in this time step

! this should be called whenever (start_ts, end_ts) \subset (start_time, end_time + tos) 
if(start_ts>et+tos .or. end_ts<st)return   ! too late or too early, nothing to do

if(start_time < end_time)then  ! we really want to test start_time .ne. end_time, but avoiding test of floats on equality
        ! segment of nonzero length
        !rels =  (st - start_time) / (end_time - start_time)  ! relative position of st in the segment (start,end)
	!sx = start_x + rels * (end_x - start_x)
	!sy = start_y + rels * (end_y - start_y)
        rels = 0.
        sx = start_x
        sy = start_y
        rele =  (et - start_time) / (end_time - start_time)    ! relative position of et in the segment (start,end)
	ex = start_x + rele * (end_x - start_x)
	ey = start_y + rele * (end_y - start_y)
else
        ! just a point
        rels = 0.
        rele = 1.
	sx = start_x
	sy = start_y
	ex = end_x
	ey = end_y
endif


cx2=unit_xf*unit_xf
cy2=unit_yf*unit_yf

axmin=coord_xf(ifts,jfts)
aymin=coord_yf(ifts,jfts)
axmax=coord_xf(ifte,jfte)
aymax=coord_yf(ifte,jfte)
if(fire_print_msg.ge.smsg)then
!$OMP CRITICAL(SFIRE_CORE_CRIT)
write(msg,'(a,2f11.6,a,2f11.6)')'IGN from ',sx,sy,' to ',ex,ey
call message(msg,level=smsg)
write(msg,'(a,2f10.2,a,2f10.2,a)')'IGN timestep [',start_ts,end_ts,'] in [',start_time,end_time,']'
call message(msg,level=smsg)
write(msg,'(a,2g13.6,a,2g13.6)')'IGN tile coord from  ',axmin,aymin,' to ',axmax,aymax
call message(msg,level=smsg)
!$OMP END CRITICAL(SFIRE_CORE_CRIT)
endif
ignited=0
dmax=0
dmin=huge(dmax)
11      format('IGN ',6(a,g17.7,1x)) 
12      format('IGN ',4(a,2g13.7,1x))
! tmperr=0.
do j=jfts,jfte   
    do i=ifts,ifte
        call check_lfn_tign_ij(i,j,'ignite_fire start',end_ts,lfn(i,j),tign(i,j))
        ax=coord_xf(i,j)
        ay=coord_yf(i,j)

        ! get d= distance from the nearest point on the ignition segment
        ! and time_ign = the ignition time there
        call nearest(d,time_ign,ax,ay,sx,sy,st,ex,ey,et,cx2,cy2)


        ! collect statistics
        dmax=max(d,dmax)
        dmin=min(d,dmin)

        ! tentative new level set function and time of ignition
        ! lfn_new_chk=d - min( radius, ros*(end_ts - time_ign) )  ! lft at end_ts
        ! tign_new_chk = time_ign + d/ros
        if(radius < ros*(end_ts - time_ign))then
            lfn_new = d - radius
            tign_new= time_ign + d/ros
        else
            lfn_new=d-ros*(end_ts - time_ign)
            tign_new=lfn_new/ros + end_ts  ! consistent even with rounding errors
        endif
        ! tmperr=max(tmperr,abs(lfn_new - lfn_new_chk)+abs(tign_new-tign_new_chk))
        !lfn_new=lfn_new_chk
        !tign_new=tign_new_chk
        
        if(fire_print_msg.ge.msglevel)then
!$OMP CRITICAL(SFIRE_CORE_CRIT)
            write(msg,*)'IGN1 i,j=',i,j,' lfn(i,j)=',lfn(i,j),' tign(i,j)=',tign(i,j)
            call message(msg,level=0)
            write(msg,*)'IGN2 i,j=',i,j,' lfn_new= ',lfn_new, ' time_ign= ',time_ign,' d=',d
            call message(msg,level=msglevel)
!$OMP END CRITICAL(SFIRE_CORE_CRIT)
        endif
        if(.not.lfn_new>0.) then
            ignited=ignited+1   ! count
        endif
        if(.not. lfn(i,j)<0. .and. lfn_new < 0.) then
            tign(i,j)=tign_new  ! newly ignited now
            if(fire_print_msg.ge.msglevel)then
!$OMP CRITICAL(SFIRE_CORE_CRIT)
                write(msg,'(a,2i6,a,2g13.6,a,f10.2,a,2f10.2,a)')'IGN ignited cell ',i,j,' at',ax,ay, &
                    ' time',tign(i,j),' in [',start_ts,end_ts,']'
                call message(msg,level=0)
                write(msg,'(a,g10.3,a,f10.2,a,2f10.2,a)')'IGN distance',d,' from ignition line at',time_ign,' in [',st,et,']'
                call message(msg,level=msglevel)
!$OMP END CRITICAL(SFIRE_CORE_CRIT)
            endif
            if(tign(i,j) < start_ts .or. tign(i,j) > end_ts )then
!$OMP CRITICAL(SFIRE_CORE_CRIT)
                write(msg,'(a,2i6,a,f11.6,a,2f11.6,a)')'WARNING ',i,j, &
                ' fixing ignition time ',tign(i,j),' outside of the time step [',start_ts,end_ts,']'
                call message (msg)
!$OMP END CRITICAL(SFIRE_CORE_CRIT)
                tign(i,j) = min(max(tign(i,j),start_ts),end_ts)
            endif
        endif
        lfn(i,j)=min(lfn(i,j),lfn_new)  ! update the level set function
        if(fire_print_msg.ge.msglevel)then
!$OMP CRITICAL(SFIRE_CORE_CRIT)
            write(msg,*)'IGN3 i,j=',i,j,' lfn(i,j)=',lfn(i,j),' tign(i,j)=',tign(i,j)
            call message(msg,level=msglevel)
!$OMP END CRITICAL(SFIRE_CORE_CRIT)
        endif
        call check_lfn_tign_ij(i,j,'ignite_fire end',end_ts,lfn(i,j),tign(i,j))
    enddo
enddo
! write(msg,*)'ignite_fire: max error ',tmperr
! call message(msg)
call check_lfn_tign("ignite_fire end:",end_ts,ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,lfn, tign)
if(fire_print_msg .ge. smsg)then
lfn_min=huge(lfn_min)
do j=jfts,jfte   
    do i=ifts,ifte
        lfn_min=min(lfn_min,lfn(i,j))
    enddo
enddo
!$OMP CRITICAL(SFIRE_CORE_CRIT)
write(msg,'(a,2g13.2,a,g10.2,a,g10.2)')'IGN units ',unit_xf,unit_yf,' m max dist ',dmax,' min',dmin
call message(msg,level=smsg)
write(msg,'(a,f6.1,a,f8.1,a,i10,a,g10.2)')'IGN radius ',radius,' time of spread',tos, &
   ' ignited nodes',ignited,' lfn min',lfn_min
call message(msg,level=smsg)
!$OMP END CRITICAL(SFIRE_CORE_CRIT)
endif
end subroutine ignite_fire

!
!***
!

subroutine check_lfn_tign(s,time_now,ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,lfn, tign)
!*** purpose: check consistency of lfn and ignition
implicit none
!*** arguments
character(len=*),intent(in)::s              ! for print
real, intent(in)::time_now                    ! end of timestep = time_now
integer, intent(in):: ifts,ifte,jfts,jfte   ! fire tile bounds
integer, intent(in):: ifms,ifme,jfms,jfme   ! array bounds
real, intent(in), dimension (ifms:ifme,jfms:jfme) :: & 
                   lfn, tign                ! level function, ignition time (state)
!*** local
integer:: i,j
!***    executable

do j=jfts,jfte   
    do i=ifts,ifte
       call check_lfn_tign_ij(i,j,s,time_now,lfn(i,j),tign(i,j))
    enddo
enddo

end subroutine check_lfn_tign 

subroutine check_lfn_tign_ij(i,j,s,time_now,lfnij,tignij)
!*** purpose: check consistency of lfn and ignition
implicit none
!*** arguments
integer, intent(in)::i,j                    ! indices
character(len=*),intent(in)::s              ! for print
real, intent(in)::time_now                  ! end of timestep = time_now
real, intent(in):: lfnij, tignij                ! level function, ignition time (state)
!*** local
character(len=128):: msg
!***    executable

if(.not. lfnij<0. .and. tignij<time_now)then
!$OMP CRITICAL(SFIRE_CORE_CRIT)
   write(msg,*)'i,j=',i,j,' lfn=',lfnij,' tign=',tignij,' time_now=',time_now
   call message(msg,level=0)
   write(msg,*)' tign-time_now=',tignij-time_now
   call message(msg,level=0)
   msg=s//': inconsistent state: lfn>=0 and tign<time_now'
   call crash(msg)
!$OMP END CRITICAL(SFIRE_CORE_CRIT)
endif

end subroutine check_lfn_tign_ij


! called from the inside of a loop, inline if possible
!DEC$ ATTRIBUTES FORCEINLINE
SUBROUTINE nearest(d,t,ax,ay,sx,sy,st,ex,ey,et,cx2,cy2)
        implicit none
!***    arguments
        real, intent(out):: d,t
        real, intent(in):: ax,ay,sx,sy,st,ex,ey,et,cx2,cy2
        ! input:
        ! ax, ay       coordinates of point a
        ! sx,sy,ex,ey  coordinates of endpoints of segment [x,y]
        ! st,et        values at the endpoints x,y
        ! cx2,cy2      x and y unit squared for computing distance 

        ! output
        ! d            the distance of a and the nearest point z on the segment [x,y]
        ! t            linear interpolation from the values st,et to the point z
        !
        ! method: compute d as the distance (ax,ay) from the midpoint (mx,my)
        ! minus a correction (because of rounding errors): |a-c|^2 = |a-m|^2 - |m-c|^2
        ! when |m-c| >= |s-e|/2 the nearest point is one of the endpoints
        ! the computation work also for the case when s=e exactly or approximately 
        !
        !
        !           a    
        !          /| \
        !     s---m-c--e
        !
        ! |m-c| = |a-m| cos (a-m,e-s) 
        !       = |a-m| (a-m).(e-s))/(|a-m|*|e-s|)
!***    local
        real:: mx,my,dam2,dames,am_es,cos2,dmc2,mcrel,mid_t,dif_t,des2,cx,cy
        character(len=128):: msg
        integer::msglevel=4
!***    executable

11      format('IGN ',6(a,g17.7,1x))
12      format('IGN ',4(a,2g13.7,1x))

        ! midpoint m = (mx,my)
        mx = (sx + ex)*0.5
        my = (sy + ey)*0.5
        dam2=(ax-mx)*(ax-mx)*cx2+(ay-my)*(ay-my)*cy2      ! |a-m|^2
        des2 = (ex-sx)*(ex-sx)*cx2+(ey-sy)*(ey-sy)*cy2          ! des2 = |e-s|^2
        dames = dam2*des2
        am_es=(ax-mx)*(ex-sx)*cx2+(ay-my)*(ey-sy)*cy2       ! am_es = (a-m).(e-s)
        if(dames>0)then
            cos2 = (am_es*am_es)/dames                  ! cos2 = cos^2 (a-m,e-s)
        else ! point a already is the midpoint
            cos2 = 0.
        endif
        dmc2 = dam2*cos2                                ! dmc2 = |m-c|^2
        if(4.*dmc2 < des2)then                          ! if |m-c|<=|e-s|/2
            ! d = sqrt(max(dam2 - dmc2,0.))               ! d=|a-m|^2 - |m-c|^2, guard rounding
            mcrel = sign(sqrt(4.*dmc2/des2),am_es)      ! relative distance of c from m 
        elseif(am_es>0)then                             ! if cos > 0, closest is e
            mcrel = 1.0 
        else                                            ! closest is s
            mcrel = -1.0
        endif
	cx = (ex + sx)*0.5 + mcrel*(ex - sx)*0.5     ! interpolate to c by going from m
	cy = (ey + sy)*0.5 + mcrel*(ey - sy)*0.5     ! interpolate to c by going from m
        d=sqrt((ax-cx)*(ax-cx)*cx2+(ay-cy)*(ay-cy)*cy2) ! |a-c|^2
	t = (et + st)*0.5 + mcrel*(et - st)*0.5     ! interpolate to c by going from m
        if(fire_print_msg.ge.msglevel)then
!$OMP CRITICAL(SFIRE_CORE_CRIT)
            write(msg,12)'find nearest to [',ax,ay,'] from [',sx,sy,'] [',ex,ey,']' ! DEB
            call message(msg,level=msglevel)
            write(msg,12)'end times',st,et,' scale squared',cx2,cy2 ! DEB
            call message(msg,level=msglevel)
            write(msg,11)'nearest at mcrel=',mcrel,'from the midpoint, t=',t ! DEB
            call message(msg,level=msglevel)
            write(msg,12)'nearest is [',cx,cy,'] d=',d ! DEB
            call message(msg,level=msglevel)
            write(msg,11)'dam2=',dam2,'des2=',des2,'dames=',dames
            call message(msg,level=msglevel)
            write(msg,11)'am_es=',am_es,'cos2=',cos2,'dmc2=',dmc2 ! DEB
            call message(msg)
!$OMP END CRITICAL(SFIRE_CORE_CRIT)
        endif
END SUBROUTINE nearest


!
!**********************
!            

subroutine fuel_left( &
    ifds,ifde,jfds,jfde, &
    ims,ime,jms,jme, &
    its,ite,jts,jte, &
    ifs,ife,jfs,jfe, &
    lfn, tign, fuel_time, time_now, fuel_frac, fire_area)
implicit none

!*** purpose: determine fraction of fuel remaining
!*** NOTE: because variables are cell centered, need halo/sync width 1 before

!*** Jan Mandel August 2007 email: jmandel@ucar.edu or Jan.Mandel@gmail.com

!*** arguments

integer, intent(in) ::ifds,ifde,jfds,jfde,its,ite,jts,jte,ims,ime &
                      ,jms,jme,ifs,ife,jfs,jfe
real, intent(in), dimension(ims:ime,jms:jme)::lfn,tign,fuel_time
real, intent(in):: time_now
real, intent(out), dimension(ifs:ife,jfs:jfe)::fuel_frac
real, intent(out), dimension(ims:ime,jms:jme):: fire_area

! ims,ime,jms,jme   in   memory dimensions
! its,ite,jts,jte   in   tile dimensions (cells where fuel_frac computed)
! ifs,ife,jfs,jfe   in   fuel_frac memory dimensions
! lfn               in   level function, at nodes at midpoints of cells
! tign              in   ignition time, at nodes at nodes at midpoints of cells
! fuel_time         in   time constant of fuel, per cell
! time_now          in   time now
! fuel_frac         out  fraction of fuel remaining, per cell
! fire_area         out  fraction of cell area on fire

!*** local

integer::i,j,ir,jr,icl,jcl,isubcl,jsubcl,i2,j2,ii,jj,its1,jts1,ite1,jte1
real::fmax,frat,helpsum1,helpsum2,fuel_left_ff,fire_area_ff,rx,ry,tignf(2,2)
real,dimension(3,3)::tff,lff
! help variables instead of arrays fuel_left_f and fire_area_f 
real::lffij,lffi1j,lffij1,lffi1j1,tifij,tifi1j,tifij1,tifi1j1,tx,ty,txx,tyy
! variables for calculation instead of lff(i,j) and tif(i,j)is lffij,tifij etc..#define IFCELLS (ite-its+1)*fuel_left_irl
#define JFCELLS (jte-jts+1)*fuel_left_jrl
character(len=128)::msg
integer::m,omp_get_thread_num
     

call check_mesh_2dim(its-1,ite+1,jts-1,jte+1,ims,ime,jms,jme)
call check_mesh_2dim(its,ite,jts,jte,ifs,ife,jfs,jfe)
call check_lfn_tign('fuel_left start',time_now,its,ite,jts,jte,ims,ime,jms,jme,lfn, tign)

! refinement
ir=fuel_left_irl
jr=fuel_left_jrl

if ((ir.ne.2).or.(jr.ne.2)) then 
   call crash('fuel_left: ir.ne.2 or jr.ne.2 ')
endif

rx=1./ir 
ry=1./jr

! interpolate level set function to finer grid
#ifdef DEBUG_OUT_FUEL_LEFT    
    call write_array_m(1,IFCELLS+1,1,JFCELLS+1,1,IFCELLS+1,1,JFCELLS+1,lff,'lff',0)
    call write_array_m(1,IFCELLS+1,1,JFCELLS+1,1,IFCELLS+1,1,JFCELLS+1,tif,'tif',0)
#endif

!
! example for ir=2:
!
!                     |      coarse cell      |
!      its-1                     its                                   ite             ite+1
! -------X------------|-----.-----X-----.-----|--........----|----------X----------|---------X
!           fine node 1           2           3                   2*(ite-its+1) 
!                fine cell  1           2                             cell 2*(ite-its+1)



!  Loop over cells in Tile 
!  Changes made by Volodymyr Kondratenko 09/24/2009
its1=max(its,ifds+1)
ite1=min(ite,ifde-1)
jts1=max(jts,jfds+1)
jte1=min(jte,jfde-1)

! jm: initialize fuel_frac - we do not compute it near the domain boundary becau!se we cannot interpolate!
do j=jts,jte
   do i=its,ite
      fuel_frac(i,j)=1.
      fire_area(i,j)=0.
   enddo
enddo

do icl=its1,ite1
  do jcl=jts1,jte1
    helpsum1=0
    helpsum2=0

    call tign_lfn_interpolation(time_now,icl,jcl,ims,ime,jms,jme, &
                                  tign,lfn,tff,lff)

    do isubcl=1,ir
      do jsubcl=1,jr 
        if(fuel_left_method.eq.1)then
          call fuel_left_cell_1( fuel_left_ff, fire_area_ff, &
     lff(isubcl,jsubcl),lff(isubcl,jsubcl+1),lff(isubcl+1,jsubcl),lff(isubcl+1,jsubcl+1), &
     tff(isubcl,jsubcl),tff(isubcl,jsubcl+1),tff(isubcl+1,jsubcl),tff(isubcl+1,jsubcl+1), &
     time_now, fuel_time(icl,jcl))
        elseif(fuel_left_method.eq.2)then
          call fuel_left_cell_2( fuel_left_ff, fire_area_ff, &
     lff(isubcl,jsubcl),lff(isubcl,jsubcl+1),lff(isubcl+1,jsubcl),lff(isubcl+1,jsubcl+1), &
     tff(isubcl,jsubcl),tff(isubcl,jsubcl+1),tff(isubcl+1,jsubcl),tff(isubcl+1,jsubcl+1), &
     time_now, fuel_time(icl,jcl)) 
! dont forget to change fire_area_ff here
        else
          call crash('fuel_left: unknown fuel_left_method')
        endif

        ! consistency check
        if(fire_area_ff.lt.-1e-6 .or.  &
          (fire_area_ff.eq.0. .and. fuel_left_ff.lt.1.-1e-6))then
!$OMP CRITICAL(SFIRE_CORE_CRIT)
           write(msg,'(a,2i6,2(a,f11.8))')'fuel_left: at node',i,j, &
              ' of refined mesh fuel burnt',1-fuel_left_ff,' fire area',fire_area_ff
!$OMP END CRITICAL(SFIRE_CORE_CRIT)
           call crash(msg)
        endif

        helpsum1=helpsum1+fuel_left_ff
        helpsum2=helpsum2+fire_area_ff
        ! write(*,*)icl,jcl,isubcl,jsubcl,fuel_left_ff,fire_area_ff
      enddo
    enddo
    ! write(*,*)icl,jcl,helpsum1,helpsum2
    fuel_frac(icl,jcl)=helpsum1 / (ir*jr) ! multiply by weight for averaging over coarse cell
    fire_area(icl,jcl)=helpsum2 / (ir*jr) ! multiply by weight for averaging over coarse cell
  enddo 
enddo



#ifdef DEBUG_OUT_FUEL_LEFT
    call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,fire_area,'fire_area',0)
    call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,fuel_frac,'fuel_frac',0)
    call write_array_m(1,IFCELLS,1,JFCELLS,1,IFCELLS,1,JFCELLS,fuel_left_f,'fuel_left_f',0)
    call write_array_m(1,IFCELLS,1,JFCELLS,1,IFCELLS,1,JFCELLS,fire_area_f,'fire_area_f',0)
#endif

! consistency check after sum
!fmax=0
!do j=jts,jte
!    do i=its,ite        
!       if(fire_area(i,j).eq.0.)then
!           if(fuel_frac(i,j).lt.1.-1e-6)then
!!$OMP CRITICAL(SFIRE_CORE_CRIT)
!               write(msg,'(a,2i6,2(a,f11.8))')'fuel_left: at node',i,j, &
!                   ' fuel burnt',1-fuel_frac(i,j),' but fire area',fire_area(i,j)
!!$OMP END CRITICAL(SFIRE_CORE_CRIT)
!               call crash(msg)
!           endif
!       else
!           frat=(1-fuel_frac(i,j))/fire_area(i,j)
!           fmax=max(fmax,frat)
!       endif
!    enddo
!enddo
!$OMP CRITICAL(SFIRE_CORE_CRIT)
write(msg,'(a,4i6,a,e15.7)')'fuel_left: tile',its,ite,jts,jte,' max fuel burnt/area',fmax 
!$OMP END CRITICAL(SFIRE_CORE_CRIT)
call message(msg)

return


end subroutine fuel_left

!
!************************
!


!
!*************************
!Subroutine that is calculating tign and lfn of four endpoints of the subcell
! that is located at isubcl,jsubcl of the cell -(icl,jcl)
!
subroutine tign_lfn_interpolation(time_now,icl,jcl,ims,ime,jms,jme, &
                                  tign,lfn,tff,lff)
real, intent(in):: time_now    ! not ignited nodes will have tign set to >= time_now
integer, intent(in) :: icl,jcl
integer, intent(in) :: ims,ime,jms,jme  ! memory dimensions of all arrays
real, intent(in), dimension(ims:ime,jms:jme)::lfn,tign
real, intent(out),dimension(3,3)::tff,lff

!     |                 |                 |       
!  -(3,1)-------------(3,2)-------------(3,3)
!     |                 |                 |
!     |   (2,1)         |      (2,2)      |       
!     |                 |                 |
!     |                 |                 |
!     |                 |                 |
!     |                 |                 |
!   (2,1)--------node-(icl,jcl)---------(2,3)-----------(icl,jcl+1)-------------|
!     |          sub-node (2,2)           |
!     |                 |                 |
!     |   (1,1)         |      (1,2)      |       each fire mesh cell is decomposed in 4
!     |                 |                 |       tff,lff is computed at the nodes of 
!     |                 |                 |       the subcells, numbered (1,1)...(3,3)
!   (1,1)-------------(1,2)-------------(1,3)--     
!     |                 |                 |       
!

!**********************
	
!       Direct calculation tif and lff, avoiding big auxiliary arrays, just for case ir=jr=2
! Checking whether icl or jcl is on the boundary

      call tign_lfn_four_pnts_interp(tign(icl-1,jcl-1),tign(icl-1,jcl),tign(icl,jcl-1),  & 
                                     tign(icl,jcl),lfn(icl-1,jcl-1),lfn(icl-1,jcl),      &
                                     lfn(icl,jcl-1),lfn(icl,jcl),lff(1,1),tff(1,1),time_now)

      call tign_lfn_line_interp(tign(icl-1,jcl),tign(icl,jcl),lfn(icl-1,jcl),lfn(icl,jcl), &
                             lff(1,2),tff(1,2),time_now)


      call tign_lfn_four_pnts_interp(tign(icl-1,jcl),tign(icl-1,jcl+1),tign(icl,jcl),  & 
                                     tign(icl,jcl+1),lfn(icl-1,jcl),lfn(icl-1,jcl+1),      &
                                     lfn(icl,jcl),lfn(icl,jcl+1),lff(1,3),tff(1,3),time_now)

      call tign_lfn_line_interp(tign(icl,jcl-1),tign(icl,jcl),lfn(icl,jcl-1),lfn(icl,jcl), &
                             lff(2,1),tff(2,1),time_now)

      lff(2,2)=lfn(icl,jcl)
      tff(2,2)=tign(icl,jcl)

      call tign_lfn_line_interp(tign(icl,jcl+1),tign(icl,jcl),lfn(icl,jcl+1),lfn(icl,jcl), &
                             lff(2,3),tff(2,3),time_now)

      call tign_lfn_four_pnts_interp(tign(icl,jcl-1),tign(icl,jcl),tign(icl+1,jcl-1),  & 
                                     tign(icl+1,jcl),lfn(icl,jcl-1),lfn(icl,jcl),      &
                                     lfn(icl+1,jcl-1),lfn(icl+1,jcl),lff(3,1),tff(3,1),time_now)

      call tign_lfn_line_interp(tign(icl+1,jcl),tign(icl,jcl),lfn(icl+1,jcl),lfn(icl,jcl), &
                             lff(3,2),tff(3,2),time_now)

      call tign_lfn_four_pnts_interp(tign(icl,jcl),tign(icl,jcl+1),tign(icl+1,jcl),  & 
                                     tign(icl+1,jcl+1),lfn(icl,jcl),lfn(icl,jcl+1),      &
                                     lfn(icl+1,jcl),lfn(icl+1,jcl+1),lff(3,3),tff(3,3),time_now)

end subroutine tign_lfn_interpolation

!
!************************
!

subroutine tign_lfn_line_interp(tign1,tign2,lfn1,lfn2,lfn_subcl,tign_subcl,time_now)
!***purpose: computes time of ignition of the point(*_subcl) that lies on the line
! between two points, whose lfn and tign is known

!*** arguments
!
real,intent(in) :: tign1,tign2   ! ignition times at the two endpoints
real,intent(in) :: lfn1,lfn2     ! level set function at the two endpoints
real,intent(in) :: time_now      ! tign>=time_now => no fire
real,intent(out) :: lfn_subcl,tign_subcl ! interpolated to midpoint

! Case 1: both points are on fire -> tign is interpolated linearly
! Case 2: both points are not on fire -> tign_subcl=time_now, not burning
! Case 3: One point is not fire, another is not -> tign_subcl set from
!         the equation  tign - time_now = c*lfn, where the proportionality
!         constant is found from the values of tign and lfn at the burning endpoint.
!         In particular, lfn(x)=0 implies tign(x)=time_now, so the representations of
!         the fireline by lfn and tign are consistent.
!         This is a special case of Case 3 from subroutine tign_lfn_four_pnts_interp.

!*** local
real :: c

!*** executable

! check consistency of inputs
call check_lfn_tign_ij(0,1,'tign_lfn_line_interp',time_now,lfn1,tign1) 
call check_lfn_tign_ij(0,2,'tign_lfn_line_interp',time_now,lfn2,tign2) 

lfn_subcl=0.5*(lfn1+lfn2)  ! interpolate lfn linearly

if (.not. lfn_subcl < 0.) then ! never test floats on equality
   tign_subcl=time_now         ! midpoint not on fire => none on fire
elseif ((lfn1 < 0.).AND.(lfn2 < 0.)) then
   tign_subcl=0.5*(tign1+tign2) ! both on fire, interpolate linearly
else  ! one is on fire one is not
   if (lfn1 < 0) then  ! 1 is on fire, 2 is not
      c = (tign1-time_now)/lfn1
   elseif (lfn2 < 0) then
      c = (tign2-time_now)/lfn2
   else
      call crash('tign_lfn_line_interp: one of lfn1 or lfn2 should be < 0')
   endif
   if( c < 0.)call crash('tign_lfn_line_interp: bad ignition times, c<0')
   tign_subcl=c*lfn_subcl+time_now;
endif
end subroutine tign_lfn_line_interp
!
!************************
!

subroutine tign_lfn_four_pnts_interp(tign1,tign2,tign3,tign4, &
  lfn1,lfn2,lfn3,lfn4,lfn_subcl,tign_subcl,time_now)

!***purpose: computes time of ignition of the point(*_subcl) that lies on the middle
! of square with given 4 points on its ends, whose lfn and tign is known
! since lfn is interpolated linearly, lfn_subcl is known 

!***arguments
real,intent(in) :: tign1,tign2,tign3,tign4 ! time of ignition of all corner points
real,intent(in) :: lfn1,lfn2,lfn3,lfn4 ! lfn of central and all corner points
real,intent(in) :: time_now
real,intent(out) :: lfn_subcl,tign_subcl

! Case 1: all 4 points are on fire -> tign is interpolated linearly
!   tign_subcl=0.25*(tign1+tign2+tign3+tign4)
! Case 2: all points are not on fire -> subcl - not burning
! Case 3: some points are on fire, others are not
! Here we assume that when lfn(x)=0 -> tign(x)=time_now (values interpolated to point x)
! For this case, tign of central point is approximated by
! tign~=c*lfn+time_now, which for our case is
! tign_subcl~=c*lfn_subcl+time_now
! where c is computed by least squares from the values at the points that are on fire
!
! sum(c*lfn(Ai)+time_now-tign(Ai))^2 ---> min
! for all lfn(Ai)<0
!
! solution for 'c' would be 
!
!         sum(lfn(Ai)*lfn(Ai)
!    c=   ------------------------------, both sums are over Ai, s.t lfn(Ai)<0
!         sum(tign(Ai)-time_now)*lfn(Ai)
!


real :: a,b,c,err
err=0.0001

call check_lfn_tign_ij(0,1,'tign_lfn_four_pnts_interp',time_now,lfn1,tign1) 
call check_lfn_tign_ij(0,2,'tign_lfn_four_pnts_interp',time_now,lfn2,tign2) 
call check_lfn_tign_ij(0,3,'tign_lfn_four_pnts_interp',time_now,lfn3,tign3) 
call check_lfn_tign_ij(0,4,'tign_lfn_four_pnts_interp',time_now,lfn4,tign4) 
 
lfn_subcl=0.25*(lfn1+lfn2+lfn3+lfn4)

if(.not. lfn_subcl < 0.) then ! midpoint not on fire, most frequent
      ! Case 2  
      tign_subcl=time_now
elseif((lfn1 < 0.).AND.(lfn2 < 0.).AND.(lfn3 < 0.).AND.(lfn4 < 0.)) then
      ! Case 1
      tign_subcl=0.25*(tign1+tign2+tign3+tign4)  ! all on fire, interpolate
else ! some on fire
      ! Case 3
      ! tign_subcl~=c*lfn_subcl+time_now
      a=0; 
      b=0;
      if (lfn1 < 0.) then
         a=a+lfn1*lfn1
         b=b+(tign1-time_now)*lfn1
      endif
      if (lfn2 < 0.) then
           a=a+lfn2*lfn2 
           b=b+(tign2-time_now)*lfn2
      endif
      if (lfn3 < 0.) then
           a=a+lfn3*lfn3
           b=b+(tign3-time_now)*lfn3
      endif
      if (lfn4 < 0.) then
           a=a+lfn4*lfn4
           b=b+(tign4-time_now)*lfn4
      endif 
      if (.not. a>0.) call crash('tign_lfn_four_pnts_interp: none of lfn < 0')
      if( b < 0.)call crash('tign_lfn_four_pnts_interp: bad ignition times, b<0') ! can have 0 because of rounding
      c=b/a;
      tign_subcl=c*lfn_subcl+time_now;
endif

end subroutine tign_lfn_four_pnts_interp


!
!************************
!


subroutine fuel_left_cell_1( fuel_frac_left, fire_frac_area, &
    lfn00,lfn01,lfn10,lfn11, &
    tign00,tign01,tign10,tign11,&
    time_now, fuel_time_cell)
!*** purpose: compute the fuel fraction left in one cell
implicit none
!*** arguments
real, intent(out):: fuel_frac_left, fire_frac_area ! 
real, intent(in)::lfn00,lfn01,lfn10,lfn11    ! level set function at 4 corners of the cell
real, intent(in)::tign00,tign01,tign10,tign11! ignition time at the  4 corners of the cell
real, intent(in)::time_now                   ! the time now
real, intent(in)::fuel_time_cell            ! time to burns off to 1/e
!*** Description
! The area burning is given by the condition L <= 0, where the function P is
! interpolated from lfn(i,j)
!
! The time since ignition is the function T, interpolated in from tign(i,j)-time_now.
! The values of tign(i,j) where lfn(i,j)>=0 are ignored, tign(i,j)=0 is taken 
! when lfn(i,j)=0.
!
! The function computes an approxmation  of the integral
!
!
!                                  /\
!                                  |              
! fuel_frac_left  =      1   -     | 1 -  exp(-T(x,y)/fuel_time_cell)) dxdy
!                                  |            
!                                 \/
!                                0<x<1
!                                0<y<1
!                             L(x,y)<=0
!
! When the cell is not burning at all (all lfn>=0), then fuel_frac(i,j)=1.
! Because of symmetries, the result should not depend on the mesh spacing dx dy
! so dx=1 and dy=1 assumed.
!
! Example:
!
!        lfn<0         lfn>0
!      (0,1) -----O--(1,1)            O = points on the fireline, T=tnow
!            |      \ |               A = the burning area for computing
!            |       \|                        fuel_frac(i,j)
!            |   A    O 
!            |        |
!            |        |
!       (0,0)---------(1,0)
!       lfn<0          lfn<0
!
! Approximations allowed: 
! The fireline can be approximated by straight line(s).
! When all cell is burning, approximation by 1 point Gaussian quadrature is OK.
! 
! Requirements:
! 1. The output should be a continuous function of the arrays lfn and
!  tign whenever lfn(i,j)=0 implies tign(i,j)=tnow.  
! 2. The output should be invariant to the symmetries of the input in each cell.
! 3. Arbitrary combinations of the signs of lfn(i,j) should work.
! 4. The result should be at least 1st order accurate in the sense that it is
!    exact if the time from ignition is a linear function.
!
! If time from ignition is approximated by polynomial in the burnt
! region of the cell, this is integral of polynomial times exponential
! over a polygon, which can be computed exactly.
!
! Requirement 4 is particularly important when there is a significant decrease
! of the fuel fraction behind the fireline on the mesh scale, because the
! rate of fuel decrease right behind the fireline is much larger 
! (exponential...). This will happen when
!
! change of time from ignition within one mesh cell / fuel_time_cell is not << 1
!
! This is the same as
!
!               mesh cell size
!  X =    -------------------------      is not << 1
!       fireline speed * fuel_time_cell
!         
!
! When X is large then the fuel burnt in one timestep in the cell is
! approximately proportional to length of  fireline in that cell.
!
! When X is small then the fuel burnt in one timestep in the cell is
! approximately proportional to the area of the burning region.
!

!*** calls
intrinsic tiny

!*** local
real::ps,aps,area,ta,out
real::t00,t01,t10,t11
real,parameter::safe=tiny(aps)
character(len=128)::msg

! the following algorithm is a very crude approximation

! minus time since ignition, 0 if no ignition yet
! it is possible to have 0 in fire region when ignitin time falls in 
! inside the time step because lfn is updated at the beginning of the time step

t00=tign00-time_now
if(lfn00>0. .or. t00>0.)t00=0.
t01=tign01-time_now
if(lfn01>0. .or. t01>0.)t01=0.
t10=tign10-time_now
if(lfn10>0. .or. t10>0.)t10=0.
t11=tign11-time_now
if(lfn11>0. .or. t11>0.)t11=0.

! approximate burning area, between 0 and 1   
ps = lfn00+lfn01+lfn10+lfn11   
aps = abs(lfn00)+abs(lfn01)+abs(lfn10)+abs(lfn11)
aps=max(aps,safe)
area =(-ps/aps+1.)/2.
area = max(area,0.) ! make sure area is between 0 and 1
area = min(area,1.)
    
! average negative time since ignition
ta=0.25*(t00+t01+t10+t11)

! exp decay in the burning area
out=1.
!if(area>0.)out=1. - area*(1. - exp(ta/fuel_time_cell))
if(area>0)out=area*exp(ta/fuel_time_cell) + (1. - area)

if(out>1.)then
!$OMP CRITICAL(SFIRE_CORE_CRIT)
    write(msg,*)'out=',out,'>1 area=',area,' ta=',ta
    call message(msg)
    write(msg,*)'tign=',tign00,tign01,tign10,tign11,' time_now=',time_now
!$OMP END CRITICAL(SFIRE_CORE_CRIT)
    call message(msg)
    !call message('WARNING: fuel_left_cell_1: fuel fraction > 1')
    call crash('fuel_left_cell_1: fuel fraction > 1')
endif

!out = max(out,0.) ! make sure out is between 0 and 1
!out = min(out,1.)

fuel_frac_left = out
fire_frac_area = area

end subroutine fuel_left_cell_1

!
!****************************************
!
! function calculation fuel_frac made by Volodymyr Kondratenko on the base of
! the Matlab code made by Jan Mandel and the fortran port by Minjeong Kim

subroutine fuel_left_cell_2( fuel_frac_left, fire_frac_area, &
                             lfn00,lfn01,lfn10,lfn11, &
                             tign00,tign01,tign10,tign11,&
                             time_now, fuel_time_cell)
!*** purpose: compute the fuel fraction left in one cell
implicit none
!*** arguments
real, intent(out):: fuel_frac_left, fire_frac_area ! 
real, intent(in)::lfn00,lfn01,lfn10,lfn11    ! level set function at 4 corners of the cell
real, intent(in)::tign00,tign01,tign10,tign11! ignition time at the  4 corners of the cell
real, intent(in)::time_now                   ! the time now
real, intent(in)::fuel_time_cell             ! burns time to 1/e
!*** Description
! The burning area is given by the condition L <= 0, where the function L is
! interpolated from lfn(i,j)
!
! The time since ignition is the function T, interpolated in from tign(i,j)-time_now.
! The values of tign(i,j) where lfn(i,j)>=0 are ignored, tign(i,j)=0 is taken 
! when lfn(i,j)=0.
!
! The function computes an approxmation  of the integral
!
!
!                                  /\
!                                  |              
! fuel_frac_left  =      1   -     | 1 -  exp(-T(x,y)/fuel_time_cell)) dxdy
!                                  |            
!                                 \/
!                                0<x<1
!                                0<y<1
!                             L(x,y)<=0
!
! When the cell is not burning at all (all lfn>=0), then fuel_frac(i,j)=1.
! Because of symmetries, the result should not depend on the mesh spacing dx dy
! so dx=1 and dy=1 assumed.
!
! Example:
!
!        lfn<0         lfn>0
!      (0,1) -----O--(1,1)            O = points on the fireline, T=tnow
!            |      \ |               A = the burning area for computing
!            |       \|                        fuel_frac(i,j)
!            |   A    O 
!            |        |
!            |        |
!       (0,0)---------(1,0)
!       lfn<0          lfn<0
!
! Approximations allowed: 
! The fireline can be approximated by straight line(s).
! When all cell is burning, approximation by 1 point Gaussian quadrature is OK.
! 
! Requirements:
! 1. The output should be a continuous function of the arrays lfn and
!  tign whenever lfn(i,j)=0 implies tign(i,j)=tnow.  
! 2. The output should be invariant to the symmetries of the input in each cell.
! 3. Arbitrary combinations of the signs of lfn(i,j) should work.
! 4. The result should be at least 1st order accurate in the sense that it is
!    exact if the time from ignition is a linear function.
!
! If time from ignition is approximated by polynomial in the burnt
! region of the cell, this is integral of polynomial times exponential
! over a polygon, which can be computed exactly.
!
! Requirement 4 is particularly important when there is a significant decrease
! of the fuel fraction behind the fireline on the mesh scale, because the
! rate of fuel decrease right behind the fireline is much larger 
! (exponential...). This will happen when
! 
! change of time from ignition within one mesh cell * fuel speed is not << 1
!
! This is the same as
!
!         mesh cell size*fuel_speed 
!         -------------------------      is not << 1
!             fireline speed
!         
!
! When X is large then the fuel burnt in one timestep in the cell is
! approximately proportional to length of  fireline in that cell.
!
! When X is small then the fuel burnt in one timestep in the cell is
! approximately proportional to the area of the burning region.

!*** calls
intrinsic tiny
#define DREAL real(kind=8)
!*** local
DREAL::ps,aps,area,ta,out
DREAL::t00,t01,t10,t11
DREAL,parameter::safe=tiny(aps)
DREAL::dx,dy ! mesh sizes
integer::i,j,k

DREAL,dimension(3)::u

DREAL::tweight,tdist
integer::kk,ll,ss
DREAL::rnorm
DREAL,dimension(8,2)::xylist,xytlist
DREAL,dimension(8)::tlist,llist,xt
DREAL,dimension(5)::xx,yy
DREAL,dimension(5)::lfn,tign

integer:: npoint
DREAL::tt,x0,y0,xts,xte,yts,yte,xt1,xt2
DREAL::lfn0,lfn1,dist,nr,s,errQ,ae,ce,ceae,a0,a1,a2,d,cet
DREAL::s1,s2,s3
DREAL::upper,lower,ah,ch,aa,cc,aupp,cupp,alow,clow
DREAL,dimension(2,2)::mQ
DREAL,dimension(2)::ut
character(len=128)::msg

!calls
intrinsic epsilon

DREAL, parameter:: zero=0.,one=1.,eps=epsilon(zero)

!!!! For finite differences by VK
DREAL::tign_middle,dt_dx,dt_dy,lfn_middle,a,b,c
DREAL:: alg_err

!*** executable

call crash('fuel_left_method=2 not working, please use fuel_left_method=1')

! check consistency
call check_lfn_tign_ij(0,0,'fuel_left_cell_2',time_now,lfn00,tign00) 
call check_lfn_tign_ij(0,1,'fuel_left_cell_2',time_now,lfn01,tign01) 
call check_lfn_tign_ij(1,0,'fuel_left_cell_2',time_now,lfn10,tign10) 
call check_lfn_tign_ij(1,1,'fuel_left_cell_2',time_now,lfn11,tign11) 

alg_err=0
dx=1
dy=1
t00=time_now-tign00
if(lfn00>=0. .or. t00<0.)t00=0.
t01=time_now-tign01
if(lfn01>=0. .or. t01<0.)t01=0.
t10=time_now-tign10
if(lfn10>=0. .or. t10<0.)t10=0.
t11=time_now-tign11
if(lfn11>=0. .or. t11<0.)t11=0.

! approximate burning area, between 0 and 1  
! was taken from fuel_left_cell_1 made by Jan, will need to be changed 
ps = lfn00+lfn01+lfn10+lfn11
aps = abs(lfn00)+abs(lfn01)+abs(lfn10)+abs(lfn11)
aps=max(aps,safe)
area =(-ps/aps+1.)/2.
area = max(area,zero) ! make sure area is between 0 and 1
area = min(area,one)

!*** case0 Do nothing
if ( lfn00>=0 .and. lfn10>=0 .and. lfn01>=0 .and. lfn11>=0 ) then
    out = 1.0 !  fuel_left, no burning
    area= 0.
!*** case4 all four coners are burning
else if (lfn00<=0 .and. lfn10<=0 .and. lfn01<=0 .and. lfn11<=0) then
! All burning
! T=u(1)*x+u(2)*y+u(3)
! t(0,0)=tign(1,1)
! t(0,fd(2))=tign(1,2)
! t(fd(1),0)=tign(2,1)
! t(fd(1),fd(2))=tign(2,2)
! t(g/2,h/2)=sum(tign(i,i))/4
! dt/dx=(1/2h)*(t10-t00+t11-t01)
! dt/dy=(1/2h)*(t01-t00+t11-t10)
! approximate T(x,y)=u(1)*x+u(2)*y+u(3) Using finite differences
! t(x,y)=t(h/2,h/2)+(x-h/2)*dt/dx+(y-h/2)*dt/dy
! u(1)=dt/dx
! u(2)=dt/dy
! u(3)=t(h/2,h/2)-h/2(dt/dx+dt/dy)

 tign_middle=(t00+t01+t10+t11)/4

 ! since mesh_size is 1 we replace fd(1) and fd(2) by 1
 dt_dx=(t10-t00+t11-t01)/2
 dt_dy=(t01-t00+t11-t10)/2

 u(1)=dt_dx
 u(2)=dt_dy
 u(3)=tign_middle-(dt_dx+dt_dy)/2

! integrate
u(1)=-u(1)/fuel_time_cell
u(2)=-u(2)/fuel_time_cell
u(3)=-u(3)/fuel_time_cell
    s1=u(1)
    s2=u(2)            
    out=exp(u(3))*intexp(s1)*intexp(s2)
    area=1
    if ( out<0 .or. out>1.0 ) then
        call message('WARNING: fuel_left_cell: case all burning: out should be between 0 and 1')
    end if
!*** case 1,2,3- other cases
!*** part of cell is burning 
else
    ! set xx, yy for the coner points
    ! move these values out of i and j loop to speed up
    ! comments for xx, yy - make center [0,0], cyclic, counterclockwise
    ! comments for lfn,tign - cyclic, counterclockwise
    xx(1) = -0.5
    xx(2) =  0.5
    xx(3) =  0.5
    xx(4) = -0.5
    xx(5) = -0.5
    yy(1) = -0.5
    yy(2) = -0.5
    yy(3) =  0.5
    yy(4) =  0.5
    yy(5) = -0.5     
    lfn(1)=lfn00
    lfn(2)=lfn10
    lfn(3)=lfn11
    lfn(4)=lfn01
    lfn(5)=lfn00
    tign(1)=t00
    tign(2)=t10
    tign(3)=t11
    tign(4)=t01
    tign(5)=t00
    npoint = 0 ! number of points in polygon

    do k=1,4
        lfn0=lfn(k  )
        lfn1=lfn(k+1)
        if ( lfn0 <= 0.0 ) then
            npoint = npoint + 1
            xylist(npoint,1)=xx(k) ! add corner to list
            xylist(npoint,2)=yy(k)
            tlist(npoint)=-tign(k) ! time since ignition
            llist(npoint)=lfn0
        end if
        if ( lfn0*lfn1 < 0 ) then
            npoint = npoint + 1
! coordinates of intersection of the fire line with segment k k+1
! lfn(t)=lfn0+t*(lfn1-lfn0)=0            
            tt=lfn0/(lfn0-lfn1)
            x0=xx(k)+( xx(k+1)-xx(k) )*tt
            y0=yy(k)+( yy(k+1)-yy(k) )*tt
            xylist(npoint,1)=x0
            xylist(npoint,2)=y0
            tlist(npoint)=0 ! now at ignition
            llist(npoint)=0 ! on fireline
        end if
    end do

    ! make the list circular and trim to size
    tlist(npoint+1)=tlist(1)
    llist(npoint+1)=llist(1)   
    xylist(npoint+1,1)=xylist(1,1)
    xylist(npoint+1,2)=xylist(1,2)

! approximate L(x,y)=u(1)*x+u(2)*y+u(3)
 lfn_middle=(lfn00+lfn01+lfn10+lfn11)/4
 dt_dx=(lfn10-lfn00+lfn11-lfn01)/2
 dt_dy=(lfn01-lfn00+lfn11-lfn10)/2
 u(1)=dt_dx
 u(2)=dt_dy
 u(3)=lfn_middle-(dt_dx+dt_dy)/2
! finding the coefficient c, reminder we work over one subcell only
! T(x,y)=c*L(x,y)+time_now
    a=0
    b=0

    if (lfn00 <= 0) then
                        a=a+lfn00*lfn00
                        if (t00 < 0) then
                        call crash('fuel_burnt_fd: tign(i1) should be less then time_now')
                        else
                        b=b+t00*lfn00
                        end if
    end if


    if (lfn01 <= 0) then
                        a=a+lfn01*lfn01
                        if (t01< 0) then
                        call crash('fuel_burnt_fd: tign(i1) should be less then time_now')
                        else
                        b=b+t01*lfn01
                        end if
    end if


    if (lfn10<=0) then
                        a=a+lfn10*lfn10
                        if (t10<0) then
                        call crash('fuel_burnt_fd: tign(i1) should be less then time_now')
                        else
                        b=b+t10*lfn10
                        end if
  end if

    if (lfn11<=0) then
                        a=a+lfn11*lfn11
                        if (t11<0) then
                        call crash('fuel_burnt_fd: tign(i1) should be less then time_now')
                        else
                        b=b+t11*lfn11
                        end if
    end if


                     if (a==0) then
                        call crash('fuel_burnt_fd: if c is on fire then one of cells should be on fire')
                     end if
                        c=b/a
                        u(1)=u(1)*c
                        u(2)=u(2)*c
                        u(3)=u(3)*c

    ! rotate to gradient on x only
    nr = sqrt(u(1)**2+u(2)**2)
    if(.not.nr.gt.eps)then
        out=1.
        goto 900
    endif
    c = u(1)/nr
    s = u(2)/nr
! rotation such that Q*u(1:2)-[something;0]
    mQ(1,1)=c
    mQ(1,2)=s
    mQ(2,1)=-s
    mQ(2,2)=c            
    ! mat vec multiplication
    call matvec(mQ,2,2,u,3,ut,2,2,2)            
    errQ = ut(2) ! should be zero            
    ae = -ut(1)/fuel_time_cell
    ce = -u(3)/fuel_time_cell ! -T(xt,yt)/fuel_time=ae*xt+ce      
    cet=ce!keep ce for later
    call matmatp(xylist,8,2,mQ,2,2,xytlist,8,2,npoint+1,2,2)            
    call sortxt( xytlist, 8,2, xt,8,npoint ) !sort ascending in x           
    out=0.0
    aupp=0.0
    cupp=0.0
    alow=0.0
    clow=0.0
    do k=1,npoint-1
        xt1=xt(k)
        xt2=xt(k+1)
        upper=0
        lower=0
        ah=0
        ch=0
        if ( xt2-xt1 > eps*100 ) then
! slice of nonzero width                
! find slice height as h=ah*x+ch
            do ss=1,npoint ! pass counterclockwise
                xts=xytlist(ss,1) ! start point of the line
                yts=xytlist(ss,2)
                xte=xytlist(ss+1,1) ! end point of the line
                yte=xytlist(ss+1,2)
                  
                if ( (xts>xt1 .and. xte>xt1) .or. &
                     (xts<xt2 .and. xte<xt2) ) then
                    aa = 0 ! do nothing
                    cc = 0
                else ! line y=a*x+c through (xts,yts), (xte,yte)
                    aa = (yts-yte)/(xts-xte)
                    cc = (xts*yte-xte*yts)/(xts-xte)                    
                    if (xte<xts) then ! upper boundary
                        aupp = aa
                        cupp = cc
                        ah=ah+aa
                        ch=ch+cc
                        upper=upper+1
                    else              ! lower boundary
                        alow = aa
                        clow = cc
                        lower=lower+1
                    end if
                end if!(xts>xt1 .and. xte>xt1)              
            end do ! ss
            ce=cet !use stored ce
!shift small amounts exp(-**) to avoid negative fuel burnt
            if (ae*xt1+ce > 0 ) then
              ce=ce-(ae*xt1+ce)!
            end if
            if (ae*xt2+ce > 0) then
            ce=ce-(ae*xt2+ce)
            end if

            ah = aupp-alow
            ch = cupp-clow  
            ! integrate (ah*x+ch)*(1-exp(ae*x+ce) from xt1 to xt2
            ! numerically sound for ae->0, ae -> infty
            ! this can be important for different model scales
            ! esp. if someone runs the model in single precision!!
            ! s1=int((ah*x+ch),x,xt1,xt2)
     …