#!/usr/bin/python3

# compute cloud fraction using HSPF equation
#
# inputs:
#	GLW	= DOWNWARD LONG WAVE FLUX AT GROUND SURFACE (W/m2)
#	T2	= TEMP at 2 M (C)
#
# output:
#	C	= Cloud fraction (range: 0-1)

# gmauger Jun 2017
import numpy as np

# constants:
SIGMA  = 5.67e-8        # W/m2/K4
KATRAD = 9.0;		# unitless

def cloud_fraction(GLW,T2):
    T2 = T2 + 273.15
    ratio = GLW / (SIGMA * KATRAD * 1E-6 * T2**6)
    
    Csq = (ratio - 1.0) / 0.0017
    
    if ( Csq < 0 ):
        C = 0
    elif ( Csq >= 100 ):
        C = 1
    else:
        C = Csq**(0.5) / 10

    return C

def cloud_fraction_vector(GLW,T2):
    T2 = T2 + 273.15    
    ratio = GLW / (SIGMA * KATRAD * 1E-6 * T2**6)

    #print(np.argwhere(np.isnan(T2)))    

    Csq = (ratio - 1.0) / 0.0017
    Csq = np.where(Csq < 0, 0, Csq)
    
    C = np.where(Csq >= 100, 1, Csq**(0.5)/10)
    
    return C