wrf-fire /wrfv2_fire/phys/module_fr_sfire_core.F

Language Fortran 77 Lines 2141
MD5 Hash 308bf6021acb7d30f760f7a60c13a95a Estimated Cost $31,202 (why?)
Repository git://github.com/jbeezley/wrf-fire.git View Raw File View Project SPDX
   1
   2
   3
   4
   5
   6
   7
   8
   9
  10
  11
  12
  13
  14
  15
  16
  17
  18
  19
  20
  21
  22
  23
  24
  25
  26
  27
  28
  29
  30
  31
  32
  33
  34
  35
  36
  37
  38
  39
  40
  41
  42
  43
  44
  45
  46
  47
  48
  49
  50
  51
  52
  53
  54
  55
  56
  57
  58
  59
  60
  61
  62
  63
  64
  65
  66
  67
  68
  69
  70
  71
  72
  73
  74
  75
  76
  77
  78
  79
  80
  81
  82
  83
  84
  85
  86
  87
  88
  89
  90
  91
  92
  93
  94
  95
  96
  97
  98
  99
 100
 101
 102
 103
 104
 105
 106
 107
 108
 109
 110
 111
 112
 113
 114
 115
 116
 117
 118
 119
 120
 121
 122
 123
 124
 125
 126
 127
 128
 129
 130
 131
 132
 133
 134
 135
 136
 137
 138
 139
 140
 141
 142
 143
 144
 145
 146
 147
 148
 149
 150
 151
 152
 153
 154
 155
 156
 157
 158
 159
 160
 161
 162
 163
 164
 165
 166
 167
 168
 169
 170
 171
 172
 173
 174
 175
 176
 177
 178
 179
 180
 181
 182
 183
 184
 185
 186
 187
 188
 189
 190
 191
 192
 193
 194
 195
 196
 197
 198
 199
 200
 201
 202
 203
 204
 205
 206
 207
 208
 209
 210
 211
 212
 213
 214
 215
 216
 217
 218
 219
 220
 221
 222
 223
 224
 225
 226
 227
 228
 229
 230
 231
 232
 233
 234
 235
 236
 237
 238
 239
 240
 241
 242
 243
 244
 245
 246
 247
 248
 249
 250
 251
 252
 253
 254
 255
 256
 257
 258
 259
 260
 261
 262
 263
 264
 265
 266
 267
 268
 269
 270
 271
 272
 273
 274
 275
 276
 277
 278
 279
 280
 281
 282
 283
 284
 285
 286
 287
 288
 289
 290
 291
 292
 293
 294
 295
 296
 297
 298
 299
 300
 301
 302
 303
 304
 305
 306
 307
 308
 309
 310
 311
 312
 313
 314
 315
 316
 317
 318
 319
 320
 321
 322
 323
 324
 325
 326
 327
 328
 329
 330
 331
 332
 333
 334
 335
 336
 337
 338
 339
 340
 341
 342
 343
 344
 345
 346
 347
 348
 349
 350
 351
 352
 353
 354
 355
 356
 357
 358
 359
 360
 361
 362
 363
 364
 365
 366
 367
 368
 369
 370
 371
 372
 373
 374
 375
 376
 377
 378
 379
 380
 381
 382
 383
 384
 385
 386
 387
 388
 389
 390
 391
 392
 393
 394
 395
 396
 397
 398
 399
 400
 401
 402
 403
 404
 405
 406
 407
 408
 409
 410
 411
 412
 413
 414
 415
 416
 417
 418
 419
 420
 421
 422
 423
 424
 425
 426
 427
 428
 429
 430
 431
 432
 433
 434
 435
 436
 437
 438
 439
 440
 441
 442
 443
 444
 445
 446
 447
 448
 449
 450
 451
 452
 453
 454
 455
 456
 457
 458
 459
 460
 461
 462
 463
 464
 465
 466
 467
 468
 469
 470
 471
 472
 473
 474
 475
 476
 477
 478
 479
 480
 481
 482
 483
 484
 485
 486
 487
 488
 489
 490
 491
 492
 493
 494
 495
 496
 497
 498
 499
 500
 501
 502
 503
 504
 505
 506
 507
 508
 509
 510
 511
 512
 513
 514
 515
 516
 517
 518
 519
 520
 521
 522
 523
 524
 525
 526
 527
 528
 529
 530
 531
 532
 533
 534
 535
 536
 537
 538
 539
 540
 541
 542
 543
 544
 545
 546
 547
 548
 549
 550
 551
 552
 553
 554
 555
 556
 557
 558
 559
 560
 561
 562
 563
 564
 565
 566
 567
 568
 569
 570
 571
 572
 573
 574
 575
 576
 577
 578
 579
 580
 581
 582
 583
 584
 585
 586
 587
 588
 589
 590
 591
 592
 593
 594
 595
 596
 597
 598
 599
 600
 601
 602
 603
 604
 605
 606
 607
 608
 609
 610
 611
 612
 613
 614
 615
 616
 617
 618
 619
 620
 621
 622
 623
 624
 625
 626
 627
 628
 629
 630
 631
 632
 633
 634
 635
 636
 637
 638
 639
 640
 641
 642
 643
 644
 645
 646
 647
 648
 649
 650
 651
 652
 653
 654
 655
 656
 657
 658
 659
 660
 661
 662
 663
 664
 665
 666
 667
 668
 669
 670
 671
 672
 673
 674
 675
 676
 677
 678
 679
 680
 681
 682
 683
 684
 685
 686
 687
 688
 689
 690
 691
 692
 693
 694
 695
 696
 697
 698
 699
 700
 701
 702
 703
 704
 705
 706
 707
 708
 709
 710
 711
 712
 713
 714
 715
 716
 717
 718
 719
 720
 721
 722
 723
 724
 725
 726
 727
 728
 729
 730
 731
 732
 733
 734
 735
 736
 737
 738
 739
 740
 741
 742
 743
 744
 745
 746
 747
 748
 749
 750
 751
 752
 753
 754
 755
 756
 757
 758
 759
 760
 761
 762
 763
 764
 765
 766
 767
 768
 769
 770
 771
 772
 773
 774
 775
 776
 777
 778
 779
 780
 781
 782
 783
 784
 785
 786
 787
 788
 789
 790
 791
 792
 793
 794
 795
 796
 797
 798
 799
 800
 801
 802
 803
 804
 805
 806
 807
 808
 809
 810
 811
 812
 813
 814
 815
 816
 817
 818
 819
 820
 821
 822
 823
 824
 825
 826
 827
 828
 829
 830
 831
 832
 833
 834
 835
 836
 837
 838
 839
 840
 841
 842
 843
 844
 845
 846
 847
 848
 849
 850
 851
 852
 853
 854
 855
 856
 857
 858
 859
 860
 861
 862
 863
 864
 865
 866
 867
 868
 869
 870
 871
 872
 873
 874
 875
 876
 877
 878
 879
 880
 881
 882
 883
 884
 885
 886
 887
 888
 889
 890
 891
 892
 893
 894
 895
 896
 897
 898
 899
 900
 901
 902
 903
 904
 905
 906
 907
 908
 909
 910
 911
 912
 913
 914
 915
 916
 917
 918
 919
 920
 921
 922
 923
 924
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727
1728
1729
1730
1731
1732
1733
1734
1735
1736
1737
1738
1739
1740
1741
1742
1743
1744
1745
1746
1747
1748
1749
1750
1751
1752
1753
1754
1755
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
1775
1776
1777
1778
1779
1780
1781
1782
1783
1784
1785
1786
1787
1788
1789
1790
1791
1792
1793
1794
1795
1796
1797
1798
1799
1800
1801
1802
1803
1804
1805
1806
1807
1808
1809
1810
1811
1812
1813
1814
1815
1816
1817
1818
1819
1820
1821
1822
1823
1824
1825
1826
1827
1828
1829
1830
1831
1832
1833
1834
1835
1836
1837
1838
1839
1840
1841
1842
1843
1844
1845
1846
1847
1848
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884
1885
1886
1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907
1908
1909
1910
1911
1912
1913
1914
1915
1916
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926
1927
1928
1929
1930
1931
1932
1933
1934
1935
1936
1937
1938
1939
1940
1941
1942
1943
1944
1945
1946
1947
1948
1949
1950
1951
1952
1953
1954
1955
1956
1957
1958
1959
1960
1961
1962
1963
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973
1974
1975
1976
1977
1978
1979
1980
1981
1982
1983
1984
1985
1986
1987
1988
1989
1990
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
2001
2002
2003
2004
2005
2006
2007
2008
2009
2010
2011
2012
2013
2014
2015
2016
2017
2018
2019
2020
2021
2022
2023
2024
2025
2026
2027
2028
2029
2030
2031
2032
2033
2034
2035
2036
2037
2038
2039
2040
2041
2042
2043
2044
2045
2046
2047
2048
2049
2050
2051
2052
2053
2054
2055
2056
2057
2058
2059
2060
2061
2062
2063
2064
2065
2066
2067
2068
2069
2070
2071
2072
2073
2074
2075
2076
2077
2078
2079
2080
2081
2082
2083
2084
2085
2086
2087
2088
2089
2090
2091
2092
2093
2094
2095
2096
2097
2098
2099
2100
2101
2102
2103
2104
2105
2106
2107
2108
2109
2110
2111
2112
2113
2114
2115
2116
2117
2118
2119
2120
2121
2122
2123
2124
2125
2126
2127
2128
2129
2130
2131
2132
2133
2134
2135
2136
2137
2138
2139
2140
!
!*** 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)
            s1 = (xt2-xt1)*((1./2.)*ah*(xt2+xt1)+ch)            
            ! s2=int((ch)*(-exp(ae*x+ce)),x,xt1,xt2)
            ceae=ce/ae;
            s2 = -ch*exp(ae*(xt1+ceae))*(xt2-xt1)*intexp(ae*(xt2-xt1))                
            ! s3=int((ah*x)*(-exp(ae*x+ce)),x,xt1,xt2)
            ! s3=int((ah*x)*(-exp(ae*(x+ceae))),x,xt1,xt2)
            ! expand in Taylor series around ae=0
            ! collect(expand(taylor(int(x*(-exp(ae*(x+ceae))),x,xt1,xt2)*ae^2,ae,4)/ae^2),ae)
            ! =(1/8*xt1^4+1/3*xt1^3*ceae+1/4*xt1^2*ceae^2-1/8*xt2^4-1/3*xt2^3*ceae-1/4*xt2^2*ceae^2)*ae^2
            !     + (-1/3*xt2^3-1/2*xt2^2*ceae+1/3*xt1^3+1/2*xt1^2*ceae)*ae 
            !     + 1/2*xt1^2-1/2*xt2^2
            !
            ! coefficient at ae^2 in the expansion, after some algebra            
            a2=(xt1-xt2)*((1./4.)*(xt1+xt2)*ceae**2+(1./3.)* &
               (xt1**2+xt1*xt2+xt2**2)*ceae+(1./8.)* &
               (xt1**3+xt1*(xt2**2)+xt1**2*xt2+xt2**3))               
            d=(ae**4)*a2
            
            if (abs(d)>eps) then
            ! since ae*xt1+ce<=0 ae*xt2+ce<=0 all fine for large ae
            ! for ae, ce -> 0 rounding error approx eps/ae^2
                s3=( exp(ae*(xt1+ceae))*(ae*xt1-1)-&
                     exp(ae*(xt2+ceae))*(ae*xt2-1) )/(ae**2)
                
            !we do not worry about rounding as xt1 -> xt2, then s3 -> 0
            else
                ! coefficient at ae^1 in the expansion
                a1=(xt1-xt2)*((1./2.)*ceae*(xt1+xt2)+(1./3.)*&
                   (xt1**2+xt1*xt2+xt2**2))
                ! coefficient at ae^0 in the expansion for ae->0
                a0=(1./2.)*(xt1-xt2)*(xt1+xt2)
                s3=a0+a1*ae+a2*ae**2; ! approximate the integral
                            end if

            s3=ah*s3                                                
            out=out+s1+s2+s3
            if(out<0. .or. out>1.) then                                  
             write(msg,'(a,g14.4,a)')'WARNING::fuel_fraction ',out,' should be between 0 and 1'
            end if!print
            if(out.ne.out.or..not.out.le.huge(out).or..not.out.ge.-huge(out)) call crash('fuel_fraction out is not a valid number')
                
        end if
    end do ! k     

            out=1-out !fuel_left
end if ! if case0, elseif case4 ,else case123

900 continue 
fuel_frac_left = out
fire_frac_area= area
if(isnotfinite(fuel_frac_left)) call crash('fuel_frac_left is not a valid number')
if(isnotfinite(fire_frac_area)) call crash('fire_frac_area is not a valid number')
end subroutine fuel_left_cell_2






!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real function intexp(ab)
implicit none
DREAL::ab
!calls
intrinsic epsilon

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

!eps = 2.2204*(10.0**(-8))!from matlab
if ( eps < abs(ab)**3/6. ) then
    intexp=(exp(ab)-1)/ab
  else
    intexp=1+ab/2.
end if
end function intexp
!
!****************************************
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!         
!**************************************** 
!       


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


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine sortxt(xytlist,nrow,ncolumn,xt,nxt,nvec)
implicit none
integer::nrow,ncolumn,nxt,nvec
DREAL,dimension(nrow,ncolumn)::xytlist
DREAL,dimension(nxt)::xt

integer::i,j
DREAL::temp

do i=1,nvec
  xt(i)=xytlist(i,1)
end do

do i=1,nvec-1
  do j=i+1,nvec
    if ( xt(i) > xt(j) ) then
         temp = xt(i)
         xt(i)=xt(j)
         xt(j)=temp
    end if
  end do
end do

end subroutine sortxt
!
!****************************************
!
subroutine matvec(A,m,n,V,nv,out,nout,nrow,ncolumn)
implicit none
integer::m,n,nv,nout,nrow,ncolumn
DREAL,dimension(m,n)::A   ! allocated m by n 
DREAL,dimension(nv)::V    ! allocated nv
DREAL,dimension(nout)::out! allocated nout 

integer::i,j

do i=1,nrow
  out(i)=0.0
  do j=1,ncolumn
    out(i)=out(i)+A(i,j)*V(j)
  end do
end do
end subroutine
!
!****************************************
!
subroutine matmatp(A,mA,nA,B,mB,nB,C,mC,nC,nrow,ncolumn,nP)
implicit none
integer::mA,nA,mB,nB,mC,nC,nrow,ncolumn,nP
DREAL,dimension(mA,nA)::A   ! allocated m by n 
DREAL,dimension(mB,nB)::B   ! allocated m by n 
DREAL,dimension(mC,nC)::C   ! allocated m by n 
integer::i,j,k
do i=1,nrow  
  do j=1,ncolumn
    C(i,j)=0.0
  do k=1,nP
    C(i,j)=C(i,j)+A(i,k)*B(j,k) ! B'
  end do
end do
end do
end subroutine

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

subroutine prop_ls( id, &                                ! for debug
                ids,ide,jds,jde, &                       ! domain dims
                ims,ime,jms,jme, &                       ! memory dims
                ips,ipe,jps,jpe, &                ! patch - nodes owned by this process 
                its,ite,jts,jte, &                       ! tile dims
                ts,dt,dx,dy,     &                       ! scalars in
                tbound,          &                       ! scalars out
                lfn_in,lfn_out,tign,ros,  &              ! arrays inout          
                fp               &
                   )
implicit none

!*** purpose: advance level function in time

! Jan Mandel August 2007 - February 2008

!*** description
!
! Propagation of closed curve by a level function method. The level function
! lfn is defined by its values at the nodes of a rectangular grid. 
! The area where lfn < 0 is inside the curve. The curve is 
! described implicitly by lfn=0. Points where the curve intersects gridlines
! can be found by linear interpolation from nodes.
!
! The level function is advanced from time ts to time ts + dt. 
!
! The level function should be initialized to (an approximation of) the signed
! distance from the curve. If the initial curve is a circle, the initial level
! function is simply the distance from the center minus the radius.
! 
! The curve moves outside with speed given by function speed_func.
!   
! Method: Godunov/ENO method for the normal motion. The timestep is checked for
! CFL condition. For a straight segment in a constant field and locally linear
! level function, the method reduces to the exact normal motion. The advantage of 
! the level set method is that it treats automatically special cases such as
! the curve approaching itself and merging components of the area inside the curve.
!
! Based on S. Osher and R. Fedkiw, Level set methods and dynamic implicit surfaces,
! Springer, 2003, Sec. 6.4, as implemented in toolboxLS for Matlab by 
! I. Mitchell, A toolbox of Level Set Methods (Version 1.1), TR-2007-11,
! Dept. Computer Science, University of British Columbia, 2007
! http://www.cs.ubc.ca/\~mitchell/Toolbo\LS
! 
  
!*** arguments 

! id                in    unique identification for prints and dumps
! ids,ide,jds,jde   in    domain dimensions
! ims,ime,jms,jme   in    memory dimensions
! its,ite,jts,jte   in    tile dimensions
! ts                in    start time
! dt                in    time step
! dx,dy             in    grid spacing
! tbound            out   bound on stable time step from CFL condition, if tbound>=dt then OK
! lfn_in,lfn_out    inout,out the level set function at nodes
! tign              inout the ignition time at nodes

! The dimensions are cell-based, the nodal value is associated with the south-west corner.
! The whole computation is on domain indices ids:ide+1,jds:jde+1.
!
! The region where new lfn and tign are computed is the tile its:ite,jts:jte 
! except when the tile is at domain upper boundary, an extra band of points is added:
! if ite=ide then region goes up to ite+1, if jte=jde then region goes up to jte+1.

! The time step requires values from 2 rows of nodes beyond the region except when at the 
! domain boundary one-sided derivatives are used. This is implemented by extending the input
! beyond the domain boundary so sufficient memory bounds must be allocated. 
! The update on all tiles can be done in parallel. To avoid the race condition (different regions
! of the same array updated by different threads), the in and out versions of the
! arrays lft and tign are distinct. If the time step dt is larger
! that the returned tbound, the routine should be called again with timestep td<=tbound, and then
! having distinct inputs and outputs comes handy.

!*** calls
!
! tend_ls
!

integer,intent(in)::id,ims,ime,jms,jme,ids,ide,jds,jde,its,ite,jts,jte,ips,ipe,jps,jpe 
real,dimension(ims:ime,jms:jme),intent(inout)::lfn_in,tign
real,dimension(ims:ime,jms:jme),intent(out)::lfn_out,ros
real,intent(in)::dx,dy,ts,dt
real,intent(out)::tbound
type(fire_params),intent(in)::fp

!*** local 
! arrays
#define IMTS its-1
#define IMTE ite+1
#define JMTS jts-1
#define JMTE jte+1
real,dimension(IMTS:IMTE,JMTS:JMTE):: tend, lfn1 ! region-sized with halo
! scalars
real::grad2,rr,tbound2,a,a1 ! a=0 euler, a=0.5 heun

real::gradx,grady,aspeed,err,aerr,time_future,time_now
real::tmp
integer::ihs,ihe,jhs,jhe
integer::ihs2,ihe2,jhs2,jhe2
integer::i,j,its1,ite1,jts1,jte1,k,kk,id1
character(len=128)::msg
integer::nfirenodes,nfireline
real::sum_err,min_err,max_err,sum_aerr,min_aerr,max_aerr   

! constants
integer,parameter :: mstep=1000, printl=1
real, parameter:: zero=0.,one=1.,eps=epsilon(zero),tol=100*eps, &
    safe=2.,rmin=safe*tiny(zero),rmax=huge(zero)/safe

! f90 intrinsic function

intrinsic max,min,sqrt,nint,epsilon,tiny,huge
  
!*** executable

!$OMP CRITICAL(SFIRE_CORE_CRIT)
write(msg,'(a8,i5,a6,i5,3(a1,i5))')'prop_ls:',id,' tile ',its,':',ite,',',jts,':',jte
!$OMP END CRITICAL(SFIRE_CORE_CRIT)
call message(msg)

call check_lfn_tign('prop_ls start',ts,its,ite,jts,jte,ims,ime,jms,jme,lfn_in,tign)

    a=fire_back_weight ! from module_fr_sfire_util
    a1=1. - a
    
    ! tend = F(lfn)

    ihs2=max(its-2,ids)   ! need lfn two beyond the tile but not outside the domain 
    ihe2=min(ite+2,ide)
    jhs2=max(jts-2,jds) 
    jhe2=min(jte+2,jde)

    ihs=max(its-1,ids)   ! compute tend one beyond the tile but not outside the domain 
    ihe=min(ite+1,ide)
    jhs=max(jts-1,jds) 
    jhe=min(jte+1,jde)


#ifdef DEBUG_OUT    
    call write_array_m(ihs,ihe,jhs,jhe,ims,ime,jms,jme,lfn_in,'lfn_in',id)
#endif

    ! check array dimensions
    call check_mesh_2dim(ihs2,ihe2,jhs2,jhe2,ims,ime,jms,jme)
    call print_2d_stats(ihs2,ihe2,jhs2,jhe2,ims,ime,jms,jme, &
                   lfn_in,'prop_ls: lfn in')
    
    ! NOTE: tend_ls will extrapolate to one node strip at domain boundaries
    ! so that it can compute gradient at domain boundaries.
    ! To avoid copying, lfn_in is declared inout.
    ! At tile boundaries that are not domain boundaries values of lfn_in two nodes
    ! outside of the tile are needed.
    id1 = id  ! for debug prints
    if(id1.ne.0)id1=id1+1000
    call  tend_ls( id1, &
    ims,ime,jms,jme, &                       ! memory dims for lfn_in
    IMTS,IMTE,JMTS,JMTE, &                   ! memory dims for tend 
    ids,ide,jds,jde, &                       ! domain dims - where lfn exists
    ips,ipe,jps,jpe, &                       ! patch - nodes owned by this process 
    ihs,ihe,jhs,jhe, &                       ! where tend computed
    ims,ime,jms,jme, &                       ! memory dims for ros 
    its,ite,jts,jte, &                       ! tile dims - where to set ros
    ts,dt,dx,dy,      &                      ! scalars in
    lfn_in, &                                ! arrays in
    tbound, &                                ! scalars out 
    tend, ros, &                              ! arrays out        
    fp         &                             ! params
)

#ifdef DEBUG_OUT    
    call write_array_m(ihs,ihe,jhs,jhe,IMTS,IMTE,JMTS,JMTE,tend,'tend1',id)
#endif

    ! Euler method, the half-step, same region as ted
    do j=jhs,jhe
        do i=ihs,ihe
            lfn1(i,j) = lfn_in(i,j) + dt*tend(i,j)
        enddo
    enddo
    
    call print_2d_stats(ihs,ihe,jhs,jhe,IMTS,IMTE,JMTS,JMTE, &
                   lfn1,'prop_ls: lfn1')
    ! tend = F(lfn1) on the tile (not beyond)

    if(id1.ne.0)id1=id1+1000
    call  tend_ls( id1,&
    IMTS,IMTE,JMTS,JMTE, &                   ! memory dims for lfn
    IMTS,IMTE,JMTS,JMTE, &                   ! memory dims for tend 
    ids,ide,jds,jde,     &                   ! domain dims - where lfn exists
    ips,ipe,jps,jpe, &                       ! patch - nodes owned by this process 
    its,ite,jts,jte, &                       ! tile dims - where is tend computed
    ims,ime,jms,jme, &                       ! memory dims for ros 
    its,ite,jts,jte, &                       ! tile dims - where is ros set
    ts+dt,dt,dx,dy,      &                   ! scalars in
    lfn1, &                                  ! arrays in
    tbound2, &                               ! scalars out 
    tend,ros, &                               ! arrays out        
    fp &
)

#ifdef DEBUG_OUT    
    call write_array_m(its,ite,jts,jte,IMTS,IMTE,JMTS,JMTE,tend,'tend2',id)
#endif

    call print_2d_stats(its,ite,jts,jte,IMTS,IMTE,JMTS,JMTE,tend,'prop_ls: tend2')
        
    tbound=min(tbound,tbound2)

!$OMP CRITICAL(SFIRE_CORE_CRIT)
    write(msg,'(a,f10.2,4(a,f7.2))')'prop_ls: time',ts,' dt=',dt,' bound',min(tbound,999.99), &
        ' dx=',dx,' dy=',dy
!$OMP END CRITICAL(SFIRE_CORE_CRIT)
    call message(msg)
    if(dt>tbound)then
!$OMP CRITICAL(SFIRE_CORE_CRIT)
        write(msg,'(2(a,f10.2))')'prop_ls: WARNING: time step ',dt, &
        ' > bound =',tbound
!$OMP END CRITICAL(SFIRE_CORE_CRIT)
        call message(msg)
    endif
    
    ! combine lfn1 and lfn_in + dt*tend -> lfn_out
    
    do j=jts,jte
        do i=its,ite
            lfn_out(i,j) = a1*lfn1(i,j) + a*(lfn_in(i,j) + dt*tend(i,j))
            lfn_out(i,j) = min(lfn_out(i,j),lfn_in(i,j)) ! fire area can only increase
        enddo
    enddo      

    ! compute ignition time by interpolation
    ! the node was not burning at start but it is burning at end
    ! interpolate from the level functions at start and at end
    ! lfn_in   is the level set function value at time ts
    ! lfn_out  is the level set function value at time ts+dt
    ! 0        is the level set function value at time tign(i,j)
    ! thus assuming the level function is approximately linear =>
    ! tign(i,j)= ts + ((ts + td) - ts) * lfn_in / (lfn_in - lfn_out)
    !        =        ts + dt * lfn_in  / (lfn_in - lfn_out)
    !        = (ts + dt) + dt * lfn_out / (lfn_in - lfn_out)
    ! the second form is better because we want make sure tign < ts+dt (used only for lfn_out < 0)

    time_now=ts+dt
    time_future=ts+2*dt
    do j=jts,jte
        do i=its,ite
            call check_lfn_tign_ij(i,j,'prop_ls before',ts,lfn_in(i,j),tign(i,j))
            if (lfn_out(i,j) < 0. ) then  ! now on fire
               if( .not. lfn_in(i,j) < 0)  then ! was not before, set ignition time by interpolation
                  tign(i,j) = time_now + dt * lfn_out(i,j) / (lfn_in(i,j) - lfn_out(i,j))
                  ! tmp = max(tmp,abs(ts + dt * lfn_in(i,j) / (lfn_in(i,j) - lfn_out(i,j)) - tign(i,j)))
               endif  !  was already burning, leave ignition time unchanged
            else ! set the ignition time outside of burning region
               tign(i,j)=time_future
            endif
            call check_lfn_tign_ij(i,j,'prop_ls after',time_now,lfn_out(i,j),tign(i,j))
        enddo
    enddo
    ! write(msg,*)'prop_ls: tign err',tmp
    ! call message(msg)
    
    ! check local speed error and stats 
    ! may not work correctly in parallel
    ! init stats
    nfirenodes=0
    nfireline=0
    sum_err=0.
    min_err=rmax
    max_err=rmin     
    sum_aerr=0.
    min_aerr=rmax
    max_aerr=rmin    
    its1=its+1
    jts1=jts+1
    ite1=ite-1
    jte1=jte-1
    ! loop over right inside of the domain
    ! cannot use values outside of the domain, would have to reflect and that
    ! would change lfn_out
    ! cannot use values outside of tile, not synchronized yet
    ! so in parallel mode the statistics is not accurate
    do j=jts1,jte1
        do i=its1,ite1
            if(lfn_out(i,j)>0.0)then   ! a point out of burning region
                if(lfn_out(i+1,j)<=0.or.lfn_out(i,j+1)<=0.or. & ! neighbor in burning region
                   lfn_out(i-1,j)<=0.or.lfn_out(i,j-1)<=0)then ! point next to fireline
                   gradx=(lfn_out(i+1,j)-lfn_out(i-1,j))/(2.0*dx) ! central differences
                   grady=(lfn_out(i,j+1)-lfn_out(i,j-1))/(2.0*dy)
                   grad2=sqrt(gradx*gradx+grady*grady)
                   aspeed = (lfn_in(i,j)-lfn_out(i,j))/(dt*max(grad2,rmin))                   
                    rr = speed_func(gradx,grady,dx,dy,i,j,fp)
                   err=aspeed-rr
                   sum_err=sum_err+err
                   min_err=min(min_err,err)
                   max_err=max(max_err,err)     
                   aerr=abs(err)
                   sum_aerr=sum_aerr+aerr
                   min_aerr=min(min_aerr,aerr)
                   max_aerr=max(max_aerr,aerr)
                   nfireline=nfireline+1
                endif
            else
                nfirenodes=nfirenodes+1
            endif
        enddo
    enddo
!!$OMP CRITICAL(SFIRE_CORE_CRIT)
!    write(msg,'(2(a,i6,f8.4))')'prop_ls: nodes burning',nfirenodes, &
!      (100.*nfirenodes)/((ite1-its1+1)*(jte1-jts1+1)),'% next to fireline',nfireline
!!$OMP END CRITICAL(SFIRE_CORE_CRIT)
!   call message(msg)
!    if(nfireline>0)then
!        call print_stat_line('speed error',its1,ite1,jts1,jte1,min_err,max_err,sum_err/nfireline)
!        call print_stat_line('abs(speed error)',its1,ite1,jts1,jte1,min_aerr,max_aerr,sum_aerr/nfireline)
!    endif

    ! check if the fire did not get to the domain boundary
    do k=-1,1,2
        ! do kk=1,(boundary_guard*(ide-ids+1))/100  ! in %
        do kk=1,boundary_guard   ! measured in cells
            i=ids+k*kk
            if(i.ge.its.and.i.le.ite)then
                do j=jts,jte
                    if(lfn_out(i,j)<=0.)goto 9
                enddo
            endif
    enddo
        ! do kk=1,(boundary_guard*(jde-jds+1))/100
        do kk=1,boundary_guard    ! measured in cells
            j=jds+k*kk
            if(j.ge.jts.and.j.le.jte)then
                do i=its,ite
                    if(lfn_out(i,j)<=0.)goto 9
                enddo
            endif
        enddo
    enddo
    goto 10
9   continue
!$OMP CRITICAL(SFIRE_CORE_CRIT)
    write(msg,'(a,i2,a,2i8)')'prop_ls: fire',boundary_guard, &
        ' cells from domain boundary at node ',i,j
!$OMP END CRITICAL(SFIRE_CORE_CRIT)
    call message(msg)     
    call crash('prop_ls: increase the fire region')
10  continue

    call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme, &
                   lfn_out,'prop_ls: lfn out')

    call check_lfn_tign('prop_ls end',ts+dt,its,ite,jts,jte,ims,ime,jms,jme,lfn_out, tign)
end subroutine prop_ls

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

subroutine tend_ls( id, &
    lims,lime,ljms,ljme, &                   ! memory dims for lfn
    tims,time,tjms,tjme, &                   ! memory dims for tend 
    ids,ide,jds,jde, &                       ! domain - nodes where lfn defined
    ips,ipe,jps,jpe, &                       ! patch - nodes owned by this process 
    ints,inte,jnts,jnte, &                   ! region - nodes where tend computed
    ims,ime,jms,jme, &                       ! memory dims for ros 
    its,ite,jts,jte, &                       ! tile dims - where is ros set
    t,dt,dx,dy,      &                       ! scalars in
    lfn, &                                   ! arrays in
    tbound, &                                ! scalars out 
    tend, ros,  &                              ! arrays out
    fp &
)

implicit none
! purpose
! compute the right hand side of the level set equation

!*** arguments
integer,intent(in)::id,lims,lime,ljms,ljme,tims,time,tjms,tjme
integer,intent(in)::ims,ime,jms,jme,its,ite,jts,jte
integer, intent(in)::ids,ide,jds,jde,ints,inte,jnts,jnte,ips,ipe,jps,jpe 
real,intent(in)::t                                     ! time
real,intent(in)::dt,dx,dy                                 ! mesh step
real,dimension(lims:lime,ljms:ljme),intent(inout)::lfn ! level set function
real,intent(out)::tbound                               ! max allowed time step
real,dimension(tims:time,tjms:tjme),intent(out)::tend  ! tendency (rhs of the level set pde)
real,dimension(ims:ime,jms:jme),intent(out)::ros  ! rate of spread 
type(fire_params),intent(in)::fp

!*** local 
real:: te,diffLx,diffLy,diffRx,diffRy, & 
   diffCx,diffCy,diff2x,diff2y,grad,rr, &
   ros_back,ros_wind,ros_slope,advx,advy,scale,nvx,nvy,speed
integer::i,j,itso,iteo,jtso,jteo
character(len=128)msg

! constants
real, parameter:: eps=epsilon(0.0)
!intrinsic epsilon
real, parameter:: zero=0.,one=1.,tol=100*eps, &
    safe=2.,rmin=safe*tiny(zero),rmax=huge(zero)/safe


! f90 intrinsic function

intrinsic max,min,sqrt,nint,tiny,huge

!*** executable
    
    ! check array dimensions
    call check_mesh_2dim(ints-1,inte+1,jnts-1,jnte+1,lims,lime,ljms,ljme)
    call check_mesh_2dim(ints,inte,jnts,jnte,tims,time,tjms,tjme)
    
    call continue_at_boundary(1,1,fire_lfn_ext_up, &   !extend by extrapolation but never down 
    lims,lime,ljms,ljme, &                ! memory dims
    ids,ide,jds,jde, &                    ! domain - nodes where lfn defined
    ips,ipe,jps,jpe, &                    ! patch - nodes owned by this process 
    ints,inte,jnts,jnte, &                ! tile - nodes update by this thread
    itso,iteo,jtso,jteo, &                ! where set now
    lfn)                                  ! array

    call print_2d_stats(itso,iteo,jtso,jteo,lims,lime,ljms,ljme, &
                   lfn,'tend_ls: lfn cont')

#ifdef DEBUG_OUT
    call write_array_m(ints-1,inte+1,jnts-1,jnte+1,lims,lime,ljms,ljme,lfn,'tend_lfn_in',id)
#endif
    
    tbound=0    
    do j=jnts,jnte
        do i=ints,inte
            ! one sided differences
            diffRx = (lfn(i+1,j)-lfn(i,j))/dx
            diffLx = (lfn(i,j)-lfn(i-1,j))/dx
            diffRy = (lfn(i,j+1)-lfn(i,j))/dy
            diffLy = (lfn(i,j)-lfn(i,j-1))/dy
            diffCx = diffLx+diffRx   ! 	TWICE CENTRAL DIFFERENCE
            diffCy = diffLy+diffRy
    
            !upwinding - select right or left derivative
            select case(fire_upwinding)
            case(0)  ! none
                grad=sqrt(diffCx**2 + diffCy**2)
            case(1) ! standard
                diff2x=select_upwind(diffLx,diffRx)
                diff2y=select_upwind(diffLy,diffRy)
                grad=sqrt(diff2x*diff2x + diff2y*diff2y)
            case(2) ! godunov per osher/fedkiw
                diff2x=select_godunov(diffLx,diffRx)
                diff2y=select_godunov(diffLy,diffRy)
                grad=sqrt(diff2x*diff2x + diff2y*diff2y)
            case(3) ! eno
                diff2x=select_eno(diffLx,diffRx)
                diff2y=select_eno(diffLy,diffRy)
                grad=sqrt(diff2x*diff2x + diff2y*diff2y)
            case(4) ! Sethian - twice stronger pushdown of bumps
                grad=sqrt(max(diffLx,0.)**2+min(diffRx,0.)**2   &
                        + max(diffLy,0.)**2+min(diffRy,0.)**2)
            case default
                grad=0.
            end select
  
            ! normal direction, from central differences
            scale=sqrt(diffCx*diffCx+diffCy*diffCy+eps) 
            nvx=diffCx/scale
            nvy=diffCy/scale
                      
            ! wind speed in direction of spread
            ! speed =  vx(i,j)*nvx + vy(i,j)*nvy
        
    
            ! get rate of spread from wind speed and slope

            call fire_ros(ros_back,ros_wind,ros_slope, &
            nvx,nvy,i,j,fp)

            rr=ros_back + ros_wind + ros_slope
            if(fire_grows_only.gt.0)rr=max(rr,0.)

            ! set ros for output
            if(i.ge.its.and.i.le.ite.and.j.ge.jts.and.j.le.jte)ros(i,j)=rr

            ! get rate of spread
            te = -rr*grad   ! normal term 

            ! cfl condition
            if (grad > 0.) then
                 tbound = max(tbound,rr*(abs(diff2x)/dx+abs(diff2y)/dy)/grad)
            endif

            ! add numerical viscosity
            te=te + fire_viscosity*abs(rr)*((diffRx-diffLx)+(diffRy-diffLy))

            tend(i,j)=te
        enddo
    enddo        

    call print_2d_stats(ints,inte,jnts,jnte,tims,time,tjms,tjme, &
                   tend,'tend_ls: tend out')

    ! the final CFL bound
    tbound = 1/(tbound+tol)

end subroutine tend_ls

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

real function select_upwind(diffLx,diffRx)
implicit none
real, intent(in):: diffLx, diffRx
real diff2x

! upwind differences, L or R if bith same sign, otherwise zero    

diff2x=0
if (diffLx>0.and.diffRx>0.)diff2x=diffLx
if (diffLx<0.and.diffRx<0.)diff2x=diffRx

select_upwind=diff2x
end function select_upwind


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


real function select_godunov(diffLx,diffRx)
implicit none
real, intent(in):: diffLx, diffRx
real diff2x,diffCx

! Godunov scheme: upwind differences, L or R or none    
! always test on > or < never = , much faster because of IEEE
! central diff >= 0 => take left diff if >0, ortherwise 0
! central diff <= 0 => take right diff if <0, ortherwise 0

diff2x=0
diffCx=diffRx+diffLx
if (diffLx>0.and..not.diffCx<0)diff2x=diffLx
if (diffRx<0.and.     diffCx<0)diff2x=diffRx

select_godunov=diff2x
end function select_godunov

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

real function select_eno(diffLx,diffRx)
implicit none
real, intent(in):: diffLx, diffRx
real diff2x

! 1st order ENO scheme

if    (.not.diffLx>0 .and. .not.diffRx>0)then
    diff2x=diffRx
elseif(.not.diffLx<0 .and. .not.diffRx<0)then
    diff2x=diffLx
elseif(.not.diffLx<0 .and. .not.diffRx>0)then
    if(.not. abs(diffRx) < abs(diffLx))then
        diff2x=diffRx
    else
        diff2x=diffLx
    endif
else
    diff2x=0.
endif

select_eno=diff2x
end function select_eno
      
!
!**************************
!

real function speed_func(diffCx,diffCy,dx,dy,i,j,fp)
!*** purpose
!    the level set method speed function
implicit none
!*** arguments
real, intent(in)::diffCx,diffCy  ! x and y coordinates of the direction of propagation
real, intent(in)::dx,dy  ! x and y coordinates of the direction of propagation
integer, intent(in)::i,j         ! indices of the node to compute the speed at
type(fire_params),intent(in)::fp
!*** local
real::scale,nvx,nvy,r
real::ros_back , ros_wind , ros_slope
real, parameter:: eps=epsilon(0.0)
!*** executable
            ! normal direction, from central differences
            scale=sqrt(diffCx*diffCx+diffCy*diffCy+eps) 
            nvx=diffCx/scale
            nvy=diffCy/scale
                      
            ! get rate of spread from wind speed and slope

            call fire_ros(ros_back,ros_wind,ros_slope, &
            nvx,nvy,i,j,fp)

            r=ros_back + ros_wind + ros_slope
            if(fire_grows_only.gt.0)r=max(r,0.)
            speed_func=r

end function speed_func

end module module_fr_sfire_core
Back to Top