#!/usr/bin/perl -I/nwprod/lib/incmod/perl

# $Id: sub_burn_data.pl,v 1.8 2004/06/22 00:54:57 bluesky Exp $

#------------------------------------------------------------------------------
# Filename:  sub_burn_data.pl
#
# Purpose:  Fill-in missing data for each burn.
#
# Subroutines:
#    fill_burn_data
#    fill_location_specific_data
#    fill_time_specific_data
#    fill_location_time_specific_data
#    fill_burn_specific_data
#    fill_fuel_specific_data
#    get_elev_netCDF
#    remove_duplicate_burns
#
#  USDA Forest Service - Pacific Northwest Wildland Fire Labratory
#  Copyright (C) 2003-2004
#------------------------------------------------------------------------------

$DIR_HOME = exists $ENV{BSKY_HOME} ? $ENV{BSKY_HOME} : '/bluesky';

use Math::Trig ;
use Math::Complex ;
# rrd (3/2/2005)
# use NetCDF ;
use Text::CSV_XS;
use IO::File;

$USHsmoke = $ENV{USHsmoke};

require "$USHsmoke/smoke_sub_emissions_model_baseclass.pl";
require "$USHsmoke/smoke_sub_emissions_model_EPM.pl";

require "$USHsmoke/smoke_sub_fuelload_baseclass.pl";
require "$USHsmoke/smoke_sub_fuelload_HARDY.pl";
require "$USHsmoke/smoke_sub_fuelload_NFDRS.pl";

sub fill_burn_data
{

    # need to fill in data for each burn
    # the data can be broken down by type as listed below

    # id key
    # (assume that this is always given)
    #     ID

    # location specific
    #     LAT,LON,TWN,STWNSub,TWND,RNG,SRNG,RNGD,SECT,ELEV,OWN,STATE,CNTY,SLOPE,SNOW

    # time specific 
    #     DATE,TIME
    
    # time and location specific
    #     RAIN,WIND,HARV

    # burn specific 
    # (requires location and time first to determine likely burn type and area)
    #     TYPE,AREA

    # fuel specific
    # (requires location, time, burn type to determine likely loadings)
    #     1HR,10HR,100HR,1kHR,10kHR,10k+HR,SHRUB,GRASS,ROT,DUFF,LITTER,FMM,FM10,FM1k,VEG

    # emissions specific
    # (use emissions model to get these)
    #     DUR,PM25,PM10,PM,CO,CO2,CH4,NMHC

    chdir("$MODELS_EMISSIONS_DIR");
    $cwd = `pwd`;
    print ("luke- in smoke_sub_burn_data.pl the cwd is $cwd \n");

    $EMISSIONS_MODEL = $EMISSIONS_TYPE->new();

    # loop over all the fires

    if ( $BSKY_NUMFIRES > 0 ) {  # check for more than zero fires so that loop doesn't act on fire zero

        # check for duplicate burns and remove
        remove_duplicate_burns ();

 	for (my $j = 0; $j < $BSKY_NUMFIRES; $j++) {

            my $delete_burn = 0;
	    my $hash_ref = $BURN[$j];
 	    debug(100, "Fill defaults for BURN[j]{ID} $$hash_ref{ID}");

	    # first fill the location specific data
            # If burn location is outside of domain, then delete burn
	    $delete_burn = fill_location_specific_data($hash_ref);
	    if ( $delete_burn ) {
                debug ( 100, "Remove $$hash_ref{ID} from list of burns");
                splice @BURN, $j, 1;
                $j--;
                $BSKY_NUMFIRES--;
            }
            next if $delete_burn;

	    # fill time specific data
            # If burn start time is outside of time period, then delete burn
	    $delete_burn = fill_time_specific_data($hash_ref);
	    if ( $delete_burn ) {
                debug ( 100, "Remove $$hash_ref{ID} from list of burns");
                splice @BURN, $j, 1;
                $j--;
                $BSKY_NUMFIRES--;
            }
            next if $delete_burn;

	    # next fill time and location specific data
	    fill_location_time_specific_data($hash_ref);

	    # next fill burn specific data
	    fill_burn_specific_data($hash_ref);

	    # next fill fuel specific data
	    fill_fuel_specific_data($hash_ref);

	    # finally compute emissions data
	    $EMISSIONS_MODEL->run_model($hash_ref, $j);
	}
    }
    debug ( 100, "New bsky_numfires = $BSKY_NUMFIRES");
    chdir("$WORKING_DIR");
    $cwd = `pwd`;
    print ("luke- in smoke_sub_burn_data.pl the cwd is $cwd \n");
};

sub fill_location_specific_data
{
    my $fire = shift;
    # keys to fill
    #     LAT,LON,TWN,STWNSub,TWND,RNG,SRNG,RNGD,SECT,ELEV,OWN,STATE,CNTY,SLOPE,SNOW

    ToDo('code to fill in location specific data needs to be revised');

    my %defaults = ( 'LAT',     '',
		     'LON',     '',
		     'TWN',     '',  
		     'STWNSub', '',  
		     'TWND',    '',
		     'RNG',     '',
		     'SRNG',    '',
		     'RNGD',    '',
		     'SECT',    '',
		     'ELEV',    '',
		     'OWN',     'Unknown',
		     'STATE',   '',
		     'CNTY',    '',
		     'SLOPE',   '10',
		     'SNOW',    '05'
		      );

    if ( ($$fire{LAT} eq "") && ($$fire{LON} eq "") ) {

       if ( ($$fire{TWN}   ne "") &&
            ($$fire{TWND}  ne "") &&
            ($$fire{RNG}   ne "") &&
            ($$fire{RNGD}  ne "") &&
            ($$fire{SECT}  ne "") &&
            ($$fire{STATE} ne "") ) {

          open (TRS2LL_INP, "> trs2ll.inp")
             or die "Could not open trs2ll.inp\n";
          printf ( TRS2LL_INP "%d%s%d%s%d,%s\n",
             $$fire{TWN}, $$fire{TWND},
             $$fire{RNG}, $$fire{RNGD},
             $$fire{SECT}, $$fire{STATE} );

          run_command ( "running trs2ll", "$MODELS_UTILITY_DIR/trs2ll/trs2ll" );

          open (TRS2LL_OUT, "< trs2ll.out")
             or die "Could not open trs2ll.out\n";
          my (@trs, @state, @lat, @lon);
          $_ = <TRS2LL_OUT>;
          ($trs, $state, $lat, $lon) = split;

          $$fire{LAT} = $lat;
          $$fire{LON} = $lon;
          debug(100, "after trs2ll, lat = $$fire{LAT}, lon = $$fire{LON}, id = $$fire{ID}");

       }

       else {
          debug (1, "BIG TIME ERROR - Insufficient location data given (Need: TRS and State, or Lat/Lon)");
          debug (100, "Fire can not be located in the domain. Remove fire ID = $$fire{ID}");
          return 1;
       }

    }

    # smo 8/1/2002 - add new WESTHEMS env var
    #              - use to check that LON is neg

    if ( $WESTHEMS eq 'TRUE' ) {
       if ( $$fire{LON} gt 0 ) {
          debug (100, "Positive burn longitude in western hemisphere for $$fire{ID}" );
          $$fire{LON} = -$$fire{LON};
          debug (100, "Burn longitude is now negative $$fire{LON}" );
       }
    }

    # Get the terrain elevation where the burn is located. Using the elevation
    # given by the burn reporting system results in inconsistencies.

    my $geofile = $MODELS_MET_DIR."/geo_CRS_2d";
    if ( -e $geofile ) {
        run_command ( "ln -s $MODELS_MET_DIR/geo_CRS_2d geo.cdf",
                      "ln -s $MODELS_MET_DIR/geo_CRS_2d geo.cdf" );
        $$fire{ELEV} = get_elev_netCDF ( $MET_MODEL, $$fire{LAT}, $$fire{LON} );
        debug (100, "Elevation from get_elev_netCDF is $$fire{ELEV}");
        run_command ( "rm ./geo.cdf", "rm ./geo.cdf" );

        if ( (0 + $$fire{ELEV}) < 0 ) {
            debug (100, "Fire is outside the domain. Remove fire ID = $$fire{ID}");
            return 1;
        }
    }

    default_fill(\%{$fire}, \%defaults);
    return 0;
    
};

sub fill_time_specific_data
{
    my $fire = shift;
    # keys to fill
    #     DATE,TIME
    my %defaults =   (
		      'DATE', $FORECAST_DATE_STR,
                      'TIME', '1000'
		      );
                      # note:  make DATE fill be from file date
    default_fill(\%{$fire}, \%defaults);

    # smo 8/1/2002 - add new TIMEZONE env var
    #              - use to convert to GMT if tz not in burn data

    debug(100, "Fire Local Date, Time = $$fire{DATE}, $$fire{TIME}");
    $my_timezone = get_timezone ($$fire{STATE});
    $my_datetime_in_seconds = convert_datetime_to_seconds ($$fire{DATE}, $$fire{TIME});
    $my_datetime = make_date_string($my_datetime_in_seconds+$my_timezone*3600,
       "%04d%02d%02d%02d%02d", "gmt");
    $$fire{DATE} = substr ( $my_datetime, 0, 8 );
    $$fire{TIME} = substr ( $my_datetime, 8, 4 );
    debug(100, "Fire GMT Date, Time = $$fire{DATE}, $$fire{TIME}");

    debug (100, "MODEL_DATEHR_STR = $MODEL_DATEHR_STR, my_datetime = $my_datetime, MODEL_END_DATEHR_STR = $MODEL_END_DATEHR_STR");
    if ( ($my_datetime/100 < $MODEL_DATEHR_STR) or ($my_datetime/100 > $MODEL_END_DATEHR_STR) ) {
        debug (100, "DateTime outside of forecast period. Remove fire ID = $$fire{ID}");
        return 1;
    } else {
        return 0;
    }
};  

sub fill_location_time_specific_data
{
    my $fire = shift;
    # keys to fill
    # (requires location specific and time specific data first)
    #     RAIN,WIND,HARV
    ToDo('code to fill in location and time data needs to be revised');
    my %defaults = qw(
		      RAIN    8 
		      WIND    6
		      HARV    199808
		      );
    default_fill(\%{$fire}, \%defaults);
};


sub fill_burn_specific_data
{
    my $fire = shift;
    ToDo('code to fill in burn specific data needs to be revised');
    # keys to fill: TYPE, AREA, DUR

    if ( ($$fire{DUR} eq "") and ($$fire{AREA} > 0) ) {

       if ( ($$fire{TYPE} eq "Landings") or
            ($$fire{TYPE} eq "Dozer Piles") or
            ($$fire{TYPE} eq "Hand Piles") ) {
           $$fire{DUR} = 90;
           debug(100, "fill DUR: fire type = $$fire{TYPE}, dur = $$fire{DUR}");
       }

       elsif ( ($$fire{TYPE} eq "Mechanical Piles") ) {
           $$fire{DUR} = 1.08*($$fire{AREA}) + 86;
           debug(100, "fill DUR: fire type = $$fire{TYPE}, dur = $$fire{DUR}");
       }

       elsif ( ($$fire{TYPE} eq "Unspecified") or
                 ($$fire{TYPE} eq "Wildlife Habitat") or
                 ($$fire{TYPE} eq "Range" ) ) {
           $$fire{DUR} = 35*log($$fire{AREA}) + 60;
           debug(100, "fill DUR: fire type = $$fire{TYPE}, dur = $$fire{DUR}");
       }

       elsif ( ($$fire{TYPE} eq "Understory") or
                 ($$fire{TYPE} eq "Broadcast") ) {
           $$fire{DUR} = 24*log($$fire{AREA}) + 131;
           debug(100, "fill DUR: fire type = $$fire{TYPE}, dur = $$fire{DUR}");
       }
       else {
#          $$fire{DUR} = 35*log($$fire{AREA}) + 60;
           $$fire{DUR} = 0.15*$$fire{AREA};
           debug(100, "fill DUR: fire type = $$fire{TYPE}, dur = $$fire{DUR}");
       }

    }

    if ( $$fire{DUR} < 1 ) {
       $$fire{DUR} = 1.5;
    }

    my %defaults = qw(
		      TYPE Underburn
		      AREA 10
		      DUR  1.5
		      );
    default_fill(\%{$fire}, \%defaults);    
};

sub fill_fuel_specific_data
{
    my $fire = shift;
    # keys to fill
    # (requires location, time, burn specific data to determine likely loadings)
    #     1HR,10HR,100HR,1kHR,10kHR,10k+HR,SHRUB,GRASS,ROT,DUFF,LITTER,FMM,FM10,FM1k,VEG

    # TODO: BUT, what if user specifically did not enter a fuel loading for one of these categories?
    debug(100, "fill_fuel_specific_data:  area of fire = $$fire{AREA}");

    my %fuel;

    if ( ($$fire{AREA}      > 0) &&
        (($$fire{'1HR'}    eq "") ||
         ($$fire{'10HR'}   eq "") ||
         ($$fire{'100HR'}  eq "") ||
         ($$fire{'1kHR'}   eq "") ||
         ($$fire{'10kHR'}  eq "") ||
         ($$fire{'10k+HR'} eq "") ||
         ($$fire{DUFF}     eq "") ||
         ($$fire{SHRUB}    eq "")) ) {

       # Only get fuel loadings if fuel loadings are not provided
       # and AREA of fire > 0.

       if ( $FUELLOAD_NFDRS_ALWAYS eq 'TRUE' ) {

          chdir "$MODELS_FUELLOAD_DIR/NFDRS";
          $cwd = `pwd`;
          print ("luke- in smoke_sub_burn_data.pl the cwd is $cwd \n");
          %fuel = $NFDRS_FUELLOAD->run_model ( $$fire{ID}, $$fire{LAT}, $$fire{LON} );
       }
       elsif ( $HARDY_FUELLOAD->is_valid_state ($$fire{STATE}) ) {

          chdir "$MODELS_FUELLOAD_DIR/HARDY";
          $cwd = `pwd`;
          print ("luke- in smoke_sub_burn_data.pl the cwd is $cwd \n");
          %fuel = $HARDY_FUELLOAD->run_model ( $$fire{ID}, $$fire{LAT}, $$fire{LON} );
       }
       else {

          chdir "$MODELS_FUELLOAD_DIR/NFDRS";
          $cwd = `pwd`;
          print ("luke- in smoke_sub_burn_data.pl the cwd is $cwd \n");
          %fuel = $NFDRS_FUELLOAD->run_model ( $$fire{ID}, $$fire{LAT}, $$fire{LON} );
       }

    }

    if ( ($fuel{'1HR'}    ne "") &&
         ($fuel{'10HR'}   ne "") &&
         ($fuel{'100HR'}  ne "") &&
         ($fuel{'1kHR'}   ne "") &&
         ($fuel{'10kHR'}  ne "") &&
         ($fuel{'10k+HR'} ne "") &&
         ($fuel{DUFF}   ne "") &&
         ($fuel{SHRUB}  ne "") ) {

        debug(100, "Fill default fuel loadings from Hardy/NFDRS, vegtype = $fuel{'VEG'}");

        $def_1hr    = $fuel{'1HR'};
        $def_10hr   = $fuel{'10HR'};
        $def_100hr  = $fuel{'100HR'};
        $def_1khr   = $fuel{'1kHR'};
        $def_10khr  = $fuel{'10kHR'};
        $def_10kphr = $fuel{'10k+HR'};
        $def_shrub  = $fuel{SHRUB};
        $def_duff   = $fuel{DUFF};
        $def_grass  = $fuel{GRASS};
        $def_veg    = $fuel{VEG};
    } else {
        if ( $$fire{AREA} > 0 ) {
           debug(100, "Fill default fuel loadings from Quartz vegtype = $fuel{'VEG'}");
           $def_1hr    = 1;
           $def_10hr   = 2.2;
           $def_100hr  = 1.6;
           $def_1khr   = 5.4;
           $def_10khr  = 24.6;
           $def_10kphr = 0.1;
           $def_shrub  = 1;
           $def_duff   = 8;
           $def_grass  = 0;
           $def_veg    = 'Quartz Complex';
        } else {
           debug(100, "AREA of fire is $$fire{AREA}, do not fill with fuel loadings");
           $def_1hr    = 0;
           $def_10hr   = 0;
           $def_100hr  = 0;
           $def_1khr   = 0;
           $def_10khr  = 0;
           $def_10kphr = 0;
           $def_shrub  = 0;
           $def_duff   = 0;
           $def_grass  = 0;
           $def_veg    = 'No Veg';
       }
    }

    my %defaults =   (
		      '1HR',     $def_1hr,
		      '10HR',    $def_10hr,
		      '100HR',   $def_100hr,
		      '1kHR',    $def_1khr,
		      '10kHR',   $def_10khr,
		      '10k+HR',  $def_10kphr,
		      'SHRUB',   $def_shrub,
		      'GRASS',   $def_grass,
		      'ROT',     '1',
		      'DUFF',    $def_duff,
		      'LITTER',  '14',
		      'FMM',     'A',
		      'FM10',    '9',
		      'FM1k',    '12',
		      'VEG',     $def_veg
		      );
    default_fill(\%{$fire}, \%defaults);    
    chdir("$MODELS_EMISSIONS_DIR");
    $cwd = `pwd`;
    print ("luke- in smoke_sub_burn_data.pl the cwd is $cwd \n");
};

sub get_elev_netCDF ()
{
    $projectionHash = shift ;
    $ilat = shift;
    $ilon = shift;
    debug (100, "get ELEV:  $ilat $ilon");

    $pi = 3.14159265 ;
    $radius = 20000000/$pi ;

    # Projection data - Lambert Conformal Conic (LCC)
    $sp1       = $projectionHash->{alpha} * $pi / 180 ; # 1st std parallel for LCC

    $sp2       = $projectionHash->{beta}  * $pi / 180 ; # 2nd std parallel for LCC
    # if std parallels are equal then divide by zero error occurs in eqn 15-3
    # check if equal and add 1 degree to beta if they are
    if ( $projectionHash->{alpha} eq $projectionHash->{beta} ) {
        $sp2 = ($projectionHash->{beta}+1)  * $pi / 180;
    }

    $olat      = $projectionHash->{latC}  * $pi / 180 ; # center lat
    $olon      = $projectionHash->{lonC}  * $pi / 180 ; # center lon
    $xllcorner = $projectionHash->{xllM}  ;             # LCC lower left corner x coord (m)
    $yllcorner = $projectionHash->{yllM}  ;             # LCC lower left corner y coord (m)
    $cellsize  = $projectionHash->{dxKM}  * 1000 ;      # grid cell size (m)
    $ncols     = $projectionHash->{nxCRS} ;             # number of columns
    $nrows     = $projectionHash->{nyCRS} ;             # number of rows

    $ncfn = "geo.cdf" ;
    $ncid = NetCDF::open($ncfn, NetCDF::NC_NOWRITE) ;
    $varid = NetCDF::varid($ncid, "ELEV") ;

    $lat = $ilat * $pi / 180;
    $lon = $ilon * $pi / 180;

    # equation (15-3) (Snyder)
    $n = log(cos($sp1)/cos($sp2)) /
          log(tan($pi/4.0 + $sp2/2.0)/tan($pi/4.0 + $sp1/2.0)) ;

    # equation (15-2)
    $F = cos($sp1) * tan($pi/4 + $sp1/2.0) ** $n / $n ;

    # equations (15-1), (14-4), (15-1a)
    $rho = $radius * $F / tan($pi/4 + $lat/2.0) ** $n ;
    $theta = $n * ($lon - $olon) ;
    $rho0 = $radius * $F / tan($pi/4 + $olat/2.0) ** $n ;

    # equations (14-1), (14-2)
    ($x, $y) = ($rho * sin($theta), $rho0 - $rho * cos($theta)) ;

    # Calculate the grid row and column.
    $gridy = int(($y-$yllcorner)/$cellsize) ;
    $gridx = int(($x-$xllcorner)/$cellsize) ;
    debug (100, "get_elev_netCDF:  gridx = $gridx, gridy = $gridy");

    if ( ($gridy < 0) || ($gridy > $nrows) ) {
       debug (1, "Error get_elev_netCDF:  y-dim outside domain");
       return -1;
    }
    elsif ( ($gridx < 0) || ($gridx > $ncols) ) {
       debug (1, "Error get_elev_netCDF:  x-dim outside domain");
       return -1;
    }
    else {
       @coords = (0, 0, $gridy, $gridx) ;
       $data = "" ;
       $junk = NetCDF::varget1($ncid, $varid, \@coords, $data) ;
       $elev = $data ;
       NetCDF::close($ncid) ;
       return $elev;
    }
}

sub remove_duplicate_burns ()
{
    # Remove duplicate burns.
    # If two fires have same ID, DATE, and TIME then delete one of them.
    # If the fires are wildfires, then delete the 209 report fire, not the manually input fire.
    # Want manual input of wildfire data to override the automated data
    # since the manually input data is (probably, hopefully!) based on better information.

    debug (200, "Entering remove_duplicate_burns");

    for (my $j = 0; $j < $BSKY_NUMFIRES; $j++) {

        my $bref1 = $BURN[$j];
        debug (200, "bref1 = $$bref1{ID}, $$bref1{DATE}, $$bref1{TIME}, $$bref1{TYPE}");

        for (my $jj = $j+1; $jj < $BSKY_NUMFIRES; $jj++) {

            my $bref2 = $BURN[$jj];
#           debug (200, "bref2 = $$bref2{ID}, $$bref2{DATE}, $$bref2{TIME}, $$bref2{TYPE}");

            if ( ($$bref1{ID}   eq $$bref2{ID})   &&
                 ($$bref1{DATE} eq $$bref2{DATE}) &&
                 ($$bref1{TIME} eq $$bref2{TIME}) ) {

                # check if fire is a wildfire, don't delete WILDFIREMANUAL, only WILDFIRE209
                if ( (substr($$bref1{TYPE}, 0) =~ /WILDFIRE/) and
                     (substr($$bref2{TYPE}, 0) =~ /WILDFIRE/) ) {

                    if ( $$bref1{TYPE} eq "WILDFIRE209" ) {
                        # delete $bref1
                        debug ( 100, "Remove bref1 $$bref1{ID}, $$bref1{TYPE} from list of burns");
                        splice @BURN, $j, 1;
                        $j--;
                        $BSKY_NUMFIRES--;
                    }
                    elsif ( $$bref2{TYPE} eq "WILDFIRE209" ) {
                        # delete $bref2
                        debug ( 100, "Remove bref2 $$bref2{ID}, $$bref2{TYPE} from list of burns");
                        splice @BURN, $jj, 1;
                        $jj--;
                        $BSKY_NUMFIRES--;
                    }
                    else {
                        # PROBLEMS!!
                        debug ( 100, "ERROR:  Duplicate wildfires identified for $$bref2{ID} but unable to delete one");
                        debug ( 100, "ERROR:  Wildfire types are:  $$bref1{TYPE} $$bref2{TYPE}" );
                    }
                }
                else {
                    # delete bref2
                    debug ( 100, "Remove bref2 $$bref2{ID} from list of burns");
                    splice @BURN, $jj, 1;
                    $jj--;
                    $BSKY_NUMFIRES--;
                }
            }
        }
    }
}

return 1;