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

# $Id: sub_fuelload_HARDY.pl,v 1.2 2004/02/17 22:53:27 smo Exp $

#-------------------------------------------------------------------------------
# Filename:  sub_fuelload_HARDY.pl
#
# Purpose:  Package implementing interface with the Hardy et al. fuel loadings.
#
# Original version by Mitchell Johnson.  Revisions by S. ONeill.
#
#  USDA Forest Service - Pacific Northwest Wildland Fire Labratory
#  Copyright (C) 2003-2004
#-------------------------------------------------------------------------------

package HARDY;

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

require "$USHsmoke/smoke_sub_fuelload_baseclass.pl";

@ISA = (Fuelload_Class);

$cwd = `pwd`;
# $cwd = system(pwd);
print ("luke in smoke_sub_fuelload_HARDY.pl the cwd is $cwd \n");

use Math::Trig ;
use Math::Complex ;
#use NetCDF ; #eliminated in DIC 2004 (NetCDF files changed to CSV files)
use Text::CSV_XS;
use IO::File;

sub new
{
    ::debug(1, "using HARDY fuelload system");
    my $self = Fuelload_Class->new;
    bless($self);
    return $self;
};


sub run_model ()
{
    $hash = shift ;
    $burnID = shift ;
    $ilat = shift;
    $ilon = shift;
    ::debug (100, "Hardy:  $burnID $ilat $ilon");

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

    # Projection data - Lambert Conformal Conic (LCC)
    $sp1 = 33 * $pi / 180 ;     # 1st std parallel for LCC
    $sp2 = 45 * $pi / 180 ;     # 2nd std parallel for LCC
    $olat = 40 * $pi / 180 ;    # center lat
    $olon = -120 * $pi / 180 ;  # center lon
    $xllcorner = -400000;       # LCC lower left corner x coord (m)
    $yllcorner = -960000;       # LCC lower left corner y coord (m)
    $cellsize = 1000;           # grid cell size (m)
    $ncols = 2030;              # number of columns
    $nrows = 2140;              # number of rows

    $vflinenumber[1]  = [-1, 36, 37, 35] ;  # Spruce/Fir
    $vflinenumber[2]  = [-1, 18, 19, 17] ;  # Mixed Conifer
    $vflinenumber[3]  = [-1, 12, 13, 11] ;  # Douglas-Fir Community
    $vflinenumber[4]  = [-1, 30, 31, 29] ;  # Ponderosa Pine
    $vflinenumber[5]  = [-1, 15, 16, 14] ;  # Lodgepole Pine
    $vflinenumber[6]  = [-1, 27, 28, 26] ;  # Pinyon/Juniper Woodland
    $vflinenumber[7]  = [-1,  5,  6,  4] ;  # Chaparral
    $vflinenumber[8]  = [-1,  3,  3,  3] ;  # Aspen/Hardwoods
    $vflinenumber[9]  = [-1,  7,  7,  7] ;  # Cottonwoods/Willow/Riparian
    $vflinenumber[10] = [-1, 21, 22, 20] ;  # Oak Brush
    $vflinenumber[11] = [-1, 33, 34, 32] ;  # Sage
    $vflinenumber[12] = [-1,  9, 10,  8] ;  # Desert Shrub
    $vflinenumber[13] = [-1,  2,  1,  1] ;  # Annual Grass
    $vflinenumber[14] = [-1, 24, 25, 23] ;  # Perennial Grass
    $vflinenumber[15] = [-1, 38, 38, 38] ;  # Alpine Water
    $vflinenumber[16] = [-1, 39, 39, 39] ;  # Water
    $vflinenumber[17] = [-1, 40, 40, 40] ;  # Barren
    $vflinenumber[18] = [-1, 41, 41, 41] ;  # Agriculture


#--------	START: PART 1 OF CHANGES TO LOAD CSV MATRIX INSTEAD OF NetCDF (DIC 2004)---------------------------------------


#     $ncfn = "hardyfuel.cdf" ;
#     $ncid = NetCDF::open($ncfn, NetCDF::NC_NOWRITE) ;
#     $veg_varid = NetCDF::varid($ncid, "VEG") ;
#     $wfl_varid = NetCDF::varid($ncid, "WFL") ;

    my $ifile_VEG = "$FIX_DIR/smoke_VEG_hardyfuel.txt";
       die "FATAL ERROR: $! unable to open file VEG_hardyfuel.txt\n"
       unless defined $ifile_VEG;

    open(INFO_VEG,$ifile_VEG);
    @matrixVEG = <INFO_VEG>;


    my $ifile_WFL = "$FIX_DIR/smoke_WFL_hardyfuel.txt";
       die "FATAL ERROR: $! unable to open file WFL_hardyfuel.txt\n"
       unless defined $ifile_WFL;

    open(INFO_WFL,$ifile_WFL);
    @matrixWFL = <INFO_WFL>;


#--------	END: PART 1 OF CHANGES TO LOAD CSV MATRIX INSTEAD OF NetCDF (DIC 2004)------------------------

    $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)) ;

    # Remembering that the X coordinate ends up being the column number.
    $vegy = int(($y-$yllcorner)/$cellsize) ;
    $vegx = int(($x-$xllcorner)/$cellsize) ;

    if ( ($vegy < 0) || ($vegy > $nrows) ) {
       ::debug (1, "Error hardy_fuelload.pl:  y-dim outside domain, specifying barren fuelload");
       $linenumber = 40 ;
    }
    elsif ( ($vegx < 0) || ($vegx > $ncols) ) {
       ::debug (1, "Error hardy_fuelload.pl:  x-dim outside domain, specifying barren fuelload");
       $linenumber = 40 ;
    }
    else {

# --------	START: PART 2 OF CHANGES TO LOAD CSV MATRIX INSTEAD OF NetCDF (DIC 2004)---------------------------------------

#        @coords = (0, 0, $vegy, $vegx) ;
#        $data = "" ;
#        $junk = NetCDF::varget1($ncid, $veg_varid, \@coords, $data) ;
#        $veg_type = $data ;
#        NetCDF::varget1($ncid, $wfl_varid, \@coords, $data) ;
#        $fuel_load = $data ;



         $row_veg=$matrixVEG[$vegy];
         split(/ /,$row_veg);
         $veg_type=$_[$vegx];

         $row_wfl=$matrixWFL[$vegy];
         split(/ /,$row_wfl);
         $fuel_load=$_[$vegx];

# --------	END: PART 2 OF CHANGES TO LOAD CSV MATRIX INSTEAD OF NetCDF (DIC 2004)---------------
       if ( $fuel_load eq 0 ) { $fuel_load = 2; }

       $linenumber = $vflinenumber[$veg_type][$fuel_load] ;

    }

    ::debug (100, "linenumber = $linenumber, veg_type = $veg_type, fuel_load = $fuel_load") ;

    if ($linenumber == -1 || $linenumber == 0) {

       ::debug (1, "Error Running hardy_fuelload.pl! - NO DATA") ;

    } else {

       # open Fuelload_lookup.csv and get the fuel loadings
       my $ifile_handle = new IO::File "< $FIX_DIR/smoke_fuelload_lookup_hardy.csv";
          die "FATAL ERROR: $! unable to open file Fuelload_lookup.csv\n"
          unless defined $ifile_handle;

       my $csv = Text::CSV_XS->new();
       my $header = $csv->getline($ifile_handle);
       my $line;
       my $lineNum;
       until($ifile_handle->eof()) {

          $lineNum++;
          $line = $csv->getline($ifile_handle);

          if ( $lineNum eq $linenumber ) {
             my %fuel;
             @fuel{@$header} = @$line;
             ::debug (100, "Hardy fuelloads = @$line");
             return (%fuel);
          }
       }
       ::debug (1, "Error obtaining Hardy Fuel loadings!");
    }

    #NetCDF::close($ncid) ; #Unnecessary with changes in NetCDF to CSV files
};


sub is_valid_state ()
{

#   method to check that the given state is covered by Hardy et al fuel loadings

    $hash = shift;
    $state = shift;

    if ( ($state eq 'WA') || ($state eq 'wa') ||
         ($state eq 'OR') || ($state eq 'or') ||
         ($state eq 'ID') || ($state eq 'id') ||
         ($state eq 'MT') || ($state eq 'mt') ||
         ($state eq 'CA') || ($state eq 'ca') ||
         ($state eq 'NV') || ($state eq 'nv') ||
         ($state eq 'AZ') || ($state eq 'az') ||
         ($state eq 'NM') || ($state eq 'nm') ||
         ($state eq 'CO') || ($state eq 'co') ||
         ($state eq 'WY') || ($state eq 'wy') ||
         ($state eq 'UT') || ($state eq 'ut') ) {

        return 1;
    }
    else {

        return 0;
    }
};