'*******************************************************************************************
'*
'* Description: VBA module containing functions to calculate the components of
'*              the Canadian Fire Weather Index system, as described in
'*
'*      Van Wagner, C.E. 1987. Development and structure of the Canadian Forest Fire
'*      Weather Index System. Can. For. Serv., Ottawa, Ont. For. Tech. Rep. 35. 37 p.
'*
'*      Equation numbers from VW87 are listed throughout, to the right of the equations in
'*      in the code.
'*   
'*      An more recent technical description can be found in:
'*      http://www.cawcr.gov.au/publications/technicalreports/CTR_010.pdf
'*
'*      This module is essentially a direct C to VBA translation of Kerry Anderson's
'*      fwi84.c code. The latitude adjustments were developed by Marty Alexander, and used
'*      over Indonesia in the paper:
'*
'*      Field, R.D., Y. Wang, O. Roswintiarti, and Guswanto. A drought-based predictor of
'*      recent haze events in western Indonesia. Atmospheric Environment, 38, 1869-1878, 2004.
'* 
'*      A technical description of the latitude adjustments can be found in Appendix 3 of:
'*      http://cfs.nrcan.gc.ca/pubwarehouse/pdfs/29152.pdf
'*
'*      Robert Field, robert.field@utoronto.ca
'*******************************************************************************************


'*******************************************************************************************
'* Function Name: FFMC
'* Description: Calculates today's Fine Fuel Moisture Code
'* Parameters:
'*    TEMP is the 12:00 LST temperature in degrees celsius
'*    RH is the 12:00 LST relative humidity in %
'*    WIND is the 12:00 LST wind speed in kph
'*    RAIN is the 24-hour acculumalted rainfall in mm, calculated at 12:00 LST
'*    FFMCo is the previous day's FFMC
'*******************************************************************************************
Public Function FFMC(ByVal TEMP As Double, _
                     ByVal RH As Double, _
                     ByVal WIND As Double, _
                     ByVal RAIN As Double, _
                     ByVal FFMCPrev As Double) As Double
    
    Dim Mo, rf, Mr, Ed, Ew, ko, kd, kl, kw, m, F As Double
          
    Mo = 147.2 * (101# - FFMCPrev) / (59.5 + FFMCPrev)               '''/* 1  '*/
    If (RAIN > 0.5) Then
        rf = RAIN - 0.5                                      '''/* 2  '*/
        If (Mo <= 150#) Then
           Mr = Mo + _
            42.5 * rf * (Exp(-100# / (251# - Mo))) * (1 - Exp(-6.93 / rf)) '''/* 3a '*/
        Else
           Mr = Mo + _
                42.5 * rf * (Exp(-100# / (251# - Mo))) * (1 - Exp(-6.93 / rf)) + _
                0.0015 * (Mo - 150#) ^ 2 * (rf) ^ (0.5) '''/* 3b '*/
        End If
        If (Mr > 250#) Then
           Mr = 250#
        End If
        Mo = Mr
    End If
    Ed = 0.942 * (RH) ^ (0.679) + _
           11# * Exp((RH - 100#) / 10#) + 0.18 * (21.1 - TEMP) * (1# - Exp(-0.115 * RH)) '''/* 4  '*/
    If (Mo > Ed) Then
        ko = 0.424 * (1# - (RH / 100#) ^ 1.7) + _
              0.0694 * (WIND) ^ (0.5) * (1# - (RH / 100#) ^ 8#) '''/* 6a '*/
        kd = ko * 0.581 * Exp(0.0365 * TEMP)                '''/* 6b '*/
        m = Ed + (Mo - Ed) * (10#) ^ (-kd)                '''/* 8  '*/
    Else
        Ew = 0.618 * (RH) ^ (0.753) + _
           10# * Exp((RH - 100#) / 10#) + _
           0.18 * (21.1 - TEMP) * (1# - Exp(-0.115 * RH))   '''/* 5  '*/
        If (Mo < Ew) Then
           kl = 0.424 * (1# - ((100# - RH) / 100#) ^ 1.7) + _
            0.0694 * (WIND) ^ 0.5 * (1 - ((100# - RH) / 100#) ^ 8#) '''/* 7a '*/
           kw = kl * 0.581 * Exp(0.0365 * TEMP)             '''/* 7b '*/
           m = Ew - (Ew - Mo) * (10#) ^ (-kw)             '''/* 9  '*/
        Else
           m = Mo
        End If
    End If
    FFMC = 59.5 * (250# - m) / (147.2 + m)                    '''/* 10 '*/

End Function

'*******************************************************************************************
'* Function Name: DMC
'* Description: Calculates today's Duff Moisture Code
'* Parameters:
'*    TEMP is the 12:00 LST temperature in degrees celsius
'*    RH is the 12:00 LST relative humidity in %
'*    RAIN is the 24-hour acculumalted rainfall in mm, calculated at 12:00 LST
'*    DMCo is the previous day's DMC
'*    Lat is the latitude in decimal degrees of the location for which calculations are being made
'*    Month is the month of Year (1..12) for the current day's calculations.
'*******************************************************************************************
Public Function DMC(ByVal TEMP As Double, _
                    ByVal RH As Double, _
                    ByVal RAIN As Double, _
                    ByVal DMCPrev As Double, _
                    ByVal MONTH As Integer, _
                    ByVal LAT As Double) As Double


    Dim re, Mo, Mr, K, B, P, Pr, Dl As Double
 
    If (RAIN > 1.5) Then
        re = 0.92 * RAIN - 1.27                                '''/* 11  '*/
        Mo = 20# + Exp(5.6348 - DMCPrev / 43.43)                  '''/* 12  '*/
        If (DMCPrev <= 33#) Then
           B = 100# / (0.5 + 0.3 * DMCPrev)                       '''/* 13a '*/
        Else
           If (DMCPrev <= 65#) Then
              B = 14# - 1.3 * (Log(DMCPrev))                      '''/* 13b '*/
           Else
              B = 6.2 * Log(DMCPrev) - 17.2                       '''/* 13c '*/
           End If
        End If
        Mr = Mo + 1000# * re / (48.77 + B * re)              '''/* 14  '*/
        Pr = 244.72 - 43.43 * Log(Mr - 20#)                  '''/* 15  '*/
        If (Pr > 0#) Then
           DMCPrev = Pr
        Else
           DMCPrev = 0#
        End If
    End If
    
    If (TEMP > -1.1) Then
        Dl = DayLength(LAT, MONTH)
        K = 1.894 * (TEMP + 1.1) * (100# - RH) * Dl * 0.000001
    Else
        K = 0#
    End If
    DMC = DMCPrev + 100# * K                                       '''/* 17  '*/
    

End Function


'*******************************************************************************************
'* Function Name: DC
'* Description: Calculates today's Drought Code
'* Parameters:
'*    TEMP is the 12:00 LST temperature in degrees celsius
'*    RAIN is the 24-hour acculumalted rainfall in mm, calculated at 12:00 LST
'*    DCo is the previous day's DC
'*    Lat is the latitude in decimal degrees of the location for which calculations are being made
'*    Month is the month of Year (1..12) for the current day's calculations.
'*******************************************************************************************
Public Function DC(ByVal TEMP As Double, _
                    ByVal RAIN As Double, _
                    ByVal DCPrev As Double, _
                    ByVal MONTH As Integer, _
                    ByVal LAT As Double) As Double

    Dim rd, Qo, Qr, V, D, Dr, Lf As Double
    

    If (RAIN > 2.8) Then
        rd = 0.83 * (RAIN) - 1.27                              '/* 18  */
        Qo = 800# * Exp(-DCPrev / 400#)                          '/* 19  */
        Qr = Qo + 3.937 * rd                                 '/* 20  */
        Dr = 400# * Log(800# / Qr)                           '/* 21  */
        If (Dr > 0#) Then
           DCPrev = Dr
        Else
           DCPrev = 0#
        End If
    End If

    Lf = DayLengthFactor(LAT, MONTH - 1)
     
    If (TEMP > -2.8) Then
        V = 0.36 * (TEMP + 2.8) + Lf                            '/* 22  */
    Else
        V = Lf
    End If
     
    If (V < 0#) Then
        V = 0#
    End If
    D = DCPrev + 0.5 * V                                        '/* 23  */

    DC = D

End Function


'*******************************************************************************************
'* Function Name: ISI
'* Description: Calculates today's Initial Spread Index
'* Parameters:
'*    WIND is the 12:00 LST wind speed in kph
'*    FFMC is the current day's FFMC
'*******************************************************************************************

Public Function ISI(ByVal WIND As Double, _
                    ByVal FFMC As Double) As Double


    Dim fWIND, m, fF, R As Double
     
    fWIND = Exp(0.05039 * WIND)                                   '''/* 24  '*/
    m = 147.2 * (101 - FFMC) / (59.5 + FFMC)                      '''/* 1   '*/
    fF = 91.9 * Exp(-0.1386 * m) * (1# + (m) ^ 5.31 / 49300000#) '''/* 25  '*/
    ISI = 0.208 * fWIND * fF                                     '''/* 26  '*/
     
     
End Function

'*******************************************************************************************
'* Function Name: BUI
'* Description: Calculates today's Buildup Index
'* Parameters:
'*    DMC is the current day's Duff Moisture Code
'*    DC is the current day's Drought Code
'*******************************************************************************************
Public Function BUI(ByVal DMC As Double, _
                    ByVal DC As Double) As Double


    Dim U As Double
    
    If (DMC <= 0.4 * DC) Then
        U = 0.8 * DMC * DC / (DMC + 0.4 * DC)                      '''/* 27a '*/
    Else
        U = DMC - (1# - 0.8 * DC / (DMC + 0.4 * DC)) _
           * (0.92 + (0.0114 * DMC) ^ 1.7)                  '''/* 27b '*/
    End If
    
    BUI = U

End Function

'*******************************************************************************************
'* Function Name: FWI
'* Description: Calculates today's Fire Weather Index
'* Parameters:
'*    ISI is current day's ISI
'*    BUI is the current day's BUI
'*******************************************************************************************

Public Function FWI(ByVal ISI As Double, _
                    ByVal BUI As Double) As Double
    
    Dim fD, B, S As Double

    If (BUI <= 80#) Then
        fD = 0.626 * (BUI) ^ 0.809 + 2#                     '''/* 28a '*/
    Else
        fD = 1000# / (25# + 108.64 * Exp(-0.023 * BUI))        '''/* 28b '*/
    End If
    B = 0.1 * ISI * fD                                        '''/* 29  '*/
    If (B > 1#) Then
        S = Exp(2.72 * (0.434 * Log(B)) ^ 0.647)         '''/* 30a '*/
    Else
        S = B                                                '''/* 30b '*/
    End If
    
    FWI = S


End Function


'*******************************************************************************************
'* Function Name: DSR
'* Description: Calculates today's Daily Severity Rating
'* Parameters:
'*    FWI is current day's FWI
'*******************************************************************************************

Public Function DSR(ByVal FWI As Double) As Double
    
    DSR = 0.0272 * (FWI ^ 1.77)             '''/* 41 '*/

End Function


'* The following two functions refer to the MEA daylength adjustment 'note'.

'*******************************************************************************************
'* Function Name: DayLengthFactor
'* Description: Calculates latitude/date dependent day length factor for Drought Code
'* Parameters:
'*      Latitude is latitude in decimal degrees of calculation location
'*      Month is integer (1..12) value of month of year for which calculation is being made
'*
'*******************************************************************************************

Private Function DayLengthFactor(ByVal Latitude As Double, _
                              ByVal MONTH As Double) As Double

    LfN = Array(-1.6, -1.6, -1.6, 0.9, 3.8, 5.8, 6.4, 5#, 2.4, 0.4, -1.6, -1.6)
    LfS = Array(6.4, 5#, 2.4, 0.4, -1.6, -1.6, -1.6, -1.6, -1.6, 0.9, 3.8, 5.8)


'    '/* Use Northern hemisphere numbers */
 '   '/* something goes wrong with >= */
    If (Latitude > 15#) Then
        retVal = LfN(MONTH)
    
'    '/* Use Equatorial numbers */
    ElseIf (Latitude <= 15 And Latitude > -15) Then
        retVal = 1.39

'    '/* Use Southern hemisphere numbers */
    ElseIf (Latitude <= -15#) Then
        retVal = LfS(MONTH)
    End If

    DayLengthFactor = retVal

End Function


'*******************************************************************************************
'* Function Name: DayLength
'* Description: Calculates latitude/date dependent day length for DMC calculation
'* Parameters:
'*      Latitude is latitude in decimal degrees of calculation location
'*      Month is integer (1..12) value of month of year for which calculation is being made
'*
'*******************************************************************************************

Private Function DayLength(ByVal Latitude As Double, _
                 ByVal MONTH As Long) As Double


    '''/* Day Length Arrays for four diff't latitude ranges '*/
     DayLength46N = Array(6.5, 7.5, 9#, 12.8, 13.9, 13.9, 12.4, 10.9, 9.4, 8#, 7#, 6#)
    DayLength20N = Array(7.9, 8.4, 8.9, 9.5, 9.9, 10.2, 10.1, 9.7, 9.1, 8.6, 8.1, 7.8)
    DayLength20S = Array(10.1, 9.6, 9.1, 8.5, 8.1, 7.8, 7.9, 8.3, 8.9, 9.4, 9.9, 10.2)
    DayLength40S = Array(11.5, 10.5, 9.2, 7.9, 6.8, 6.2, 6.5, 7.4, 8.7, 10#, 11.2, 11.8)

    ''/* default to return error code '*/
    retVal = errVal
    
    ''/*
    '    Use four ranges which respectively span:
    '        - 90N to 33 N
    '        - 33 N to 0
    '        - 0 to -30
    '        - -30 to -90
    '*/
    If ((Latitude <= 90) And (Latitude > 33#)) Then
    
        retVal = DayLength46N(MONTH - 1)
    
    ElseIf ((Latitude <= 33#) And (Latitude > 15#)) Then
        retVal = DayLength20N(MONTH - 1)
        
    ElseIf ((Latitude <= 15#) And (Latitude > -15#)) Then
        retVal = 9#
    
    ElseIf ((Latitude <= -15#) And (Latitude > -30#)) Then
    
        retVal = DayLength20S(MONTH - 1)
    ElseIf ((Latitude <= -30#) And (Latitude >= -90#)) Then
        retVal = DayLength40S(MONTH - 1)
    End If

    DayLength = retVal
    

End Function








'* Daylength utility functions


'*******************************************************************************************
'* Function Name: Date2Julian
'* Description: Calculates Julian from calendar date
'* Parameters: Month, Day and Year of date to convert
'*
'*******************************************************************************************

Private Function Date2Julian(ByVal MONTH As Long, _
                             ByVal Day As Long, _
                             ByVal Year As Long) As Long


    DaysInMonth = Array(0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334)


    Julian = DaysInMonth(MONTH - 1) + Day

    If ((Year Mod 4) = 0 And MONTH > 1) Then
        Julian = Julian + 1
    End If

    Date2Julian = Julian
   
End Function


'''/*-------------------------------------------------------------------
'   This subroutine was modified from a C program called sunrise.c
'   written by Bear Giles (bear@fsl.noaa.gov) with algorithms for solar
'   declination, equation of time, and length of day provided by
'   Joseph Bartlo (jabartlo@delphi.com)
'   -------------------------------------------------------------------'*/
Private Function DayLengthExact(ByVal LAT As Double, _
                             ByVal lng As Double, _
                             ByVal jday As Double) As Double


    Dim dtr, x, decl, eot, noon, T As Double

    If (LAT < -80# Or LAT > 80# Or lng < -180# Or lng > 180#) Then
        DayLengthExact = errVal

    Else
        '''/* degrees to radians '*/
        dtr = 3.1415926 / 180#

        '''/* fraction of a year '*/
        x = dtr * 360# / 366# * jday

        '''/* solar declination '*/
        decl = dtr * (0.33029 _
            - 22.9717 * Cos(x) + 3.8346 * Sin(x) _
            - 0.3495 * Cos(2# * x) + 0.0261 * Sin(2# * x) _
            - 0.1392 * Cos(3# * x) + 0.0727 * Sin(3# * x))

        '''/* equation of time ? I don't recognize the significance of this '*/
        eot = -0.0001062 + 0.009554 * Cos(x) - 0.122671 * Sin(x) _
            - 0.05335 * Cos(2# * x) - 0.156121 * Sin(2# * x)
        noon = 12# - (lng / 15#) - eot

        '''/* ... so we'll just use Solar noon for now '*/
        noon = 12#

        '''/* fraction of a day in sunlight '*/
        T = -Tan(dtr * (LAT)) * Tan(decl) _
            + -Sin(dtr * 0.8) / (Cos(dtr * (LAT)) * Cos(decl)) _
        
        If (T < -1# Or T > 1#) Then
            DayLengthExact = errVal
        Else
            
            DayLengthExact = 2# * (12# / 3.1415926) * Acos(T)
        End If
    
    End If



End Function
                             
Private Function Acos(ByVal x As Double) As Double
    Acos = Atn(-x / Sqr(-x * x + 1)) + 2 * Atn(1)
End Function