ArithModel types

[jpfairbanks]

From Paul Hein at UA we got the following julia code for PETPT and PETASCE. I think these are good usecases for using ModelingToolkit.Equation to represent all the math. The open question is what to do about the conditionals.

function petpt(msalb, srad, tmax, tmin, xhlai, eo)
	td = (0.6*tmax)+(0.4*tmin)
	if (xhlai <= 0.0)
		albedo = msalb
	else
		albedo = 0.23-((0.23-msalb)*exp(-((0.75*xhlai))))
	end
	slang = srad*23.923
	eeq = (slang*(0.000204-(0.000183*albedo)))*(td+29.0)
	if (tmax > 35.0)
		eo = eeq*(((tmax-35.0)*0.05)+1.1)
	else
		eo = eeq*1.1
	end
	if (tmax < 5.0)
		eo = (eeq*0.01)*exp((0.18*(tmax+20.0)))
	else
		eo = eo
	end
	eo = max(eo, 0.0001)
end

function petasce(canht, doy, msalb, meevp, srad, tdew, tmax, tmin, windht, windrun, xhlai, xlat, xelev, eo)
	tavg = (tmax+tmin)/2.0
	patm = 101.3*(((293.0-(0.0065*xelev))/293.0)^5.26)
	psycon = 0.000665*patm
	udelta = (2503.0*exp(((17.27*tavg)/(tavg+237.3))))/((tavg+237.3)^2.0)
	emax = 0.6108*exp(((17.27*tmax)/(tmax+237.3)))
	emin = 0.6108*exp(((17.27*tmin)/(tmin+237.3)))
	es = (emax+emin)/2.0
	ea = 0.6108*exp(((17.27*tdew)/(tdew+237.3)))
	rhmin = max(20.0, min(80.0, ((ea/emax)*100.0)))
	if (xhlai <= 0.0)
		albedo = msalb
	else
		albedo = 0.23
	end
	rns = (1.0-albedo)*srad
	pie = 3.14159265359
	dr = 1.0+(0.033*cos((((2.0*pie)/365.0)*doy)))
	ldelta = 0.409*sin(((((2.0*pie)/365.0)*doy)-1.39))
	ws = acos(-(((1.0*tan(((xlat*pie)/180.0)))*tan(ldelta))))
	ra1 = (ws*sin(((xlat*pie)/180.0)))*sin(ldelta)
	ra2 = (cos(((xlat*pie)/180.0))*cos(ldelta))*sin(ws)
	ra = (((24.0/pie)*4.92)*dr)*(ra1+ra2)
	rso = (0.75+(2e-05*xelev))*ra
	ratio = srad/rso
	if (ratio < 0.3)
		ratio = 0.3
	else
		ratio = ratio
	end
	if (ratio > 1.0)
		ratio = 1.0
	else
		ratio = ratio
	end
	fcd = (1.35*ratio)-0.35
	tk4 = (((tmax+273.16)^4.0)+((tmin+273.16)^4.0))/2.0
	rnl = ((4.901e-09*fcd)*(0.34-(0.14*sqrt(ea))))*tk4
	rn = rns-rnl
	g = 0.0
	windsp = (((windrun*1000.0)/24.0)/60.0)/60.0
	wind2m = windsp*(4.87/log(((67.8*windht)-5.42)))
	cn = 0.0
	cd = 0.0
	if (meevp == "A")
		cd = 0.38
	else
		cd = cd
	end
	if (meevp == "A")
		cn = 1600.0
	else
		cn = cn
	end
	if (meevp == "G")
		cd = 0.34
	else
		cd = cd
	end
	if (meevp == "G")
		cn = 900.0
	else
		cn = cn
	end
	refet = ((0.408*udelta)*(rn-g))+(((psycon*(cn/(tavg+273.0)))*wind2m)*(es-ea))
	refet = refet/(udelta+(psycon*(1.0+(cd*wind2m))))
	refet = max(0.0001, refet)
	skc = 0.8
	kcbmin = 0.3
	kcbmax = 1.2
	if (xhlai <= 0.0)
		kcb = 0.0
	else
		kcb = max(0.0, (kcbmin+((kcbmax-kcbmin)*(1.0-exp(-(((1.0*skc)*xhlai)))))))
	end
	wnd = max(1.0, min(wind2m, 6.0))
	cht = max(0.001, canht)
	kcmax = 0.0
	if (meevp == "A")
		kcmax = max(1.0, (kcb+0.05))
	else
		kcmax = kcmax
	end
	if (meevp == "G")
		kcmax = max((1.2+(((0.04*(wnd-2.0))-(0.004*(rhmin-45.0)))*((cht/3.0)^0.3))), (kcb+0.05))
	else
		kcmax = kcmax
	end
	if (kcb <= kcbmin)
		fc = 0.0
	else
		fc = ((kcb-kcbmin)/(kcmax-kcbmin))^(1.0+(0.5*canht))
	end
	fw = 1.0
	few = min((1.0-fc), fw)
	ke = max(0.0, min((1.0*(kcmax-kcb)), (few*kcmax)))
	eo = (kcb+ke)*refet
	eo = max(eo, 0.0001)
end