# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Risk assessment and error analysis
#
#
# ![Untitled%20presentation%20%281%29.png](attachment:Untitled%20presentation%20%281%29.png)
#
#
# ## First import libraries

import numpy as np
# %matplotlib inline
import numpy as np
from coare3 import coare3
from air_specific_humidity import air_humidity_method_qsat26air
from surface_specific_humidity import sea_humidity_method_qsee
from temperature_structure import temperature_structure_method_psit_30 

# # define variables

# +
#variable	value	units	notes
Cp = 1004.67         #J/K/kg specific heat at constant pressure
Rgas = 287.1

wind = 7.00000          #m/s    average wind speed over ocean
Tsea = 289.25         #K      average SST over ocean
Tair = Tsea-1.3        #K      average surface air temperature over ocean
Tsea_degC = Tsea - 273.15
Tair_degC = Tair - 273.15
rel_hum = 0.80000       #N/A    near surface relative humidity
pres = 1013.00000       #mb     average surface pressure used in Bentamy

#calculate Qsea and Qair from Tsea, pres, Tair, rel_hum from coare
Lv = (2.501-0.00237*Tsea)*1e6 #J/kg   latent heat of vaporization  SB had 2500000 CG changed to coare3 def
sat_vap_pres =  6.1121*np.exp(17.502*Tair_degC/(Tair_degC+240.97))*(1.0007+3.46e-6*pres) #from coare3
Qsea = sea_humidity_method_qsee(Tsea_degC,pres)/1000
Qair = air_humidity_method_qsat26air(Tair_degC,pres,rel_hum*100)/1000
air_density = pres*100./(Rgas*(Tair)*(1.+0.61*Qair))

#cH_over_cE = 0.00120    #N/A    aerodynamic transfer coef
#Qsea=18.0/1000
#Qair=15.0/1000
#Qsea=18.0 #/1000 = .0018
#Qair=15.0       = .0015
#rh=80.0
print('Qair',Qair*1000)
print('Qsea',Qsea*1000)

# -

# #  Observational error for either FluxSat or Current variables
# - errors from FluxSat via Shannon's simulations
# - errors from current variables from publications cited in paper
#

# +
#pick one and comment out the other with #
#ob_err = 'fluxsat'
ob_err = 'current'

#fluxsat
if ob_err == 'fluxsat':
    dTsea = 0.45
    dwind = 0.6
    dTair = 0.7
    dQair = .95
    dQair_computed = Qair * 0.07 #percent from shannon
    dpres = 5
if ob_err == 'current':
    dTsea = 0.5 #ref found
    dwind = 0.8 #ref found
    dTair = 1.42
    dQair = 1.25
    dQair_computed = 1.25/1000.
    dpres = 5

dpres_computed = (pres+dpres)*100 / (287*Tair) - air_density
print('Qair error:', dQair_computed*1000)
# -

#compute Qsea errors from coare3 using uncertainty in SST & pres
Qsea0 = sea_humidity_method_qsee(Tsea_degC,pres)
Qsea1 = sea_humidity_method_qsee(Tsea_degC+dTsea,pres)
Qsea2 = sea_humidity_method_qsee(Tsea_degC-dTsea,pres)
Qsea_dTsea = (np.maximum(np.abs(Qsea0-Qsea2),np.abs(Qsea0-Qsea1)))
Qsea0 = sea_humidity_method_qsee(Tsea_degC,pres)
Qsea1 = sea_humidity_method_qsee(Tsea_degC,pres+dpres)
Qsea2 = sea_humidity_method_qsee(Tsea_degC,pres-dpres)
Qsea_dpres = (np.maximum(np.abs(Qsea0-Qsea2),np.abs(Qsea0-Qsea1)))
dQsea_computed = np.sqrt(Qsea_dTsea**2+Qsea_dpres**2)/1000
print('Qsea and Qsea uncertainty:',Qsea*1000,dQsea_computed*1000)

# # Error propagation for each variable (<NAME>)

# calculate Ch and Ce using coare
inputs = {'u':wind,'t':Tair_degC,'Q':Qair*1000,'Qs':Qsea*1000,'P':pres,'ts':Tsea_degC,'lat':35.}
res = coare3(inputs)
Ch = res['Ch']
Ce = res['Ce']
print(Ch,Ce)


# # Compute Sensible Heat Flux errors

# +
dSH_over_dTsea = air_density * Cp * Ch * wind               * dTsea
dSH_over_dTair = air_density * Cp * Ch * wind               * dTair
dSH_over_dwind = air_density * Cp * Ch        * (Tsea-Tair) * dwind
dSH_over_dpres =               Cp * Ch * wind * (Tsea-Tair) * dpres_computed

total_error = np.sqrt(dSH_over_dTsea**2+
                      dSH_over_dTair**2+
                      dSH_over_dwind**2+
                      dSH_over_dpres**2)

signal = air_density * Cp * Ch * (Tsea-Tair) * wind

percent_error = total_error / signal *100

print('dSH_over_dTsea',dSH_over_dTsea)
print('dSH_over_dTair',dSH_over_dTair)
print('dSH_over_dwind',dSH_over_dwind)
print('dSH_over_dpres',dSH_over_dpres)
print('Sensible Heat Flux total_error',total_error)
print('signal',signal)
print('percent error',percent_error)

sensible_heat_flux_total_error = total_error
# -

# # Compute Latent Heat Flux errors

# +
dLH_over_dqs =   air_density * Ce * wind * Lv               * dQsea_computed
dLH_over_dqair = air_density * Ce * wind * Lv               * dQair_computed
dLH_over_dwind = air_density * Ce *        Lv * (Qsea-Qair) * dwind
dLH_over_dpres =               Ce * wind * Lv * (Qsea-Qair) * dpres_computed

total_error = np.sqrt(dLH_over_dqs**2+
                      dLH_over_dqair**2+
                      dLH_over_dwind**2+
                      dLH_over_dpres**2)

signal = air_density * Lv * Ce * (Qsea-Qair) * wind

percent_error = total_error / signal *100

print('dLH_over_dqs',dLH_over_dqs)
print('dLH_over_dqair',dLH_over_dqair)
print('dLH_over_dwind',dLH_over_dwind)
print('dLH_over_dpres',dLH_over_dpres)
print('Latent Heat Flux total_error',total_error)
print('signal',signal)
print('percent error',percent_error)

latent_heat_flux_total_error = total_error
# -

# # Gradient error calculation from <NAME>.

print('LHF error',latent_heat_flux_total_error)
print('#','gradient_error','spatial_scale (km)')
for number_of_points in range(2,20):
    sampling_distance = 10 # spacing (in swath) between observations
    observational_error = latent_heat_flux_total_error #W/m2 Uncertainty in the observation

    sqpt = 2 * number_of_points**2
    one_over = 1. / np.sqrt(sqpt)
    sqnn2 = np.sqrt( (2*number_of_points) / (2*number_of_points - 2) )
    prod = one_over * sqnn2
    uncertainty_gradient = (observational_error / sampling_distance) * prod
    spatial_scale = number_of_points * sampling_distance
    print(number_of_points,'    ', "%1.2f" % uncertainty_gradient,'        ', spatial_scale, )


# # Below here is develoment area







# # Computed Errors

# +
dQair_computed = (Qsea*rel_hum*dQair)/100

num = 17.67*(Tsea_degC + dTsea) 
dem = (Tsea_degC + 243.5 + dTsea) 

dQsea_computed  = (6.112*np.exp(num/dem)) / pres*0.622 - Qsea
dpres_computed = (pres+dpres)*100 / (287*Tair) - air_density

print('dqair:',dQair_computed)
print('dqs:', dQsea_computed)
print('dpres:',dpres_computed)
# -



# # use COARE model
#
#    Args:
# *         inputs (dic): inputs parameter containing u,us,ts,t,Qs,Q,Rs,Rl,rain,zi,P,zu,zt,zq,lat,jcool,jwave,twave,hwave fields
# -            u (float): wind speed (m/s)  at height zu (m)
# -            us (float): surface current speed in the wind direction (m/s)
# -            ts (float): bulk water temperature (C) if jcool=1, interface water T if jcool=0  
# -            t (float): bulk air temperature (C), height zt
# -            Qs (float): bulk water spec hum (g/kg) if jcool=1, ...
# -            Q (float): bulk air spec hum (g/kg), height zq
# -            Rs (float): downward solar flux (W/m^2)
# -            Rl (float): downard IR flux (W/m^2)
# -            rain (float): rain rate (mm/hr)
# -            zi (float): Planet Boundary Layer depth (m)
# -            P (float): Atmos surface pressure (mb)
# -            zu (float): wind speed measurement height (m)
# -            zt (float): air T measurement height (m)
# -            zq (float): air q measurement height (m)
# -            lat (float): latitude (deg, N=+)
# -            jcool (float): implement cool calculation skin switch, 0=no, 1=yes
# -            jwave (float): implement wave dependent roughness model
# -            twave (float): wave period (s)
# -            hwave (float): wave height (m)
#     
#

# +
inputs = {'u':wind,'t':Tair_degC,'Q':Qair*1000,'Qs':Qsea*1000,'P':pres,'ts':Tsea_degC,'lat':35.}
#print(inputs)
res = coare3(inputs)
inputs['P']=pres+dpres
res2 = coare3(inputs)
dif2 = abs(res['hsb']-res2['hsb'])
dif2a = abs(res['hlb']-res2['hlb'])
inputs['P']=pres-dpres
res2 = coare3(inputs)
dif3 = abs(res['hsb']-res2['hsb'])
dif3a = abs(res['hlb']-res2['hlb'])
SH_perror = max(dif2,dif3)
LH_perror = max(dif2a,dif3a)
print('error due to pressure SENS:',SH_perror)

inputs = {'u':wind,'t':Tair_degC,'Q':Qair*1000,'Qs':Qsea*1000,'P':pres,'ts':Tsea_degC,'lat':35.}
res = coare3(inputs)
inputs['u']=wind+dwind
res2 = coare3(inputs)
dif2 = abs(res['hsb']-res2['hsb'])
dif2a = abs(res['hlb']-res2['hlb'])
inputs['u']=wind-dwind
res2 = coare3(inputs)
dif3 = abs(res['hsb']-res2['hsb'])
dif3a = abs(res['hlb']-res2['hlb'])
SH_winderror = max(dif2,dif3)
LH_winderror = max(dif2a,dif3a)
print('error due to wind SENS:',SH_winderror)

inputs = {'u':wind,'t':Tair_degC,'Q':Qair*1000,'Qs':Qsea*1000,'P':pres,'ts':Tsea_degC,'lat':35.}
res = coare3(inputs)
inputs['ts']=Tsea_degC+dTsea
res2 = coare3(inputs)
dif2 = abs(res['hsb']-res2['hsb'])
dif2a = abs(res['hlb']-res2['hlb'])
inputs['ts']=Tsea_degC-dTsea
res2 = coare3(inputs)
dif3 = abs(res['hsb']-res2['hsb'])
dif3a = abs(res['hlb']-res2['hlb'])
#print('DIF:', res['hsb']-res3['hsb'],res['hlb']-res3['hlb'])
SH_Tseaerror = max(dif2,dif3)
print('error due to Tsea SENS:',SH_Tseaerror)

inputs = {'u':wind,'t':Tair_degC,'Q':Qair*1000,'Qs':Qsea*1000,'P':pres,'ts':Tsea_degC,'lat':35.}
res = coare3(inputs)
inputs['t']=Tair_degC+dTair
res2 = coare3(inputs)
dif2 = abs(res['hsb']-res2['hsb'])
dif2a = abs(res['hlb']-res2['hlb'])
inputs['t']=Tair_degC-dTair
res2 = coare3(inputs)
dif3 = abs(res['hsb']-res2['hsb'])
dif3a = abs(res['hlb']-res2['hlb'])
#print('DIF:', res['hsb']-res3['hsb'],res['hlb']-res3['hlb'])
SH_Tairerror = max(dif2,dif3)
#LH_Tairerror = max(dif2a,dif3a)
print('error due to Tair SENS:',SH_Tairerror)

total_error = np.sqrt(SH_Tairerror**2+
                      SH_Tseaerror**2+
                      SH_winderror**2+
                      SH_perror**2)
print('Sensible total error',total_error)

# +
inputs = {'u':wind,'t':Tair_degC,'Q':Qair*1000,'Qs':Qsea*1000,'P':pres,'ts':Tsea_degC,'lat':35.}
#print(inputs)
res = coare3(inputs)
inputs['P']=pres+dpres
res2 = coare3(inputs)
dif2 = abs(res['hsb']-res2['hsb'])
dif2a = abs(res['hlb']-res2['hlb'])
inputs['P']=pres-dpres
res2 = coare3(inputs)
dif3 = abs(res['hsb']-res2['hsb'])
dif3a = abs(res['hlb']-res2['hlb'])
SH_perror = max(dif2,dif3)
LH_perror = max(dif2a,dif3a)
print('error due to pressure LAT:',LH_perror)

inputs = {'u':wind,'t':Tair_degC,'Q':Qair*1000,'Qs':Qsea*1000,'P':pres,'ts':Tsea_degC,'lat':35.}
res = coare3(inputs)
inputs['u']=wind+dwind
res2 = coare3(inputs)
dif2 = abs(res['hsb']-res2['hsb'])
dif2a = abs(res['hlb']-res2['hlb'])
inputs['u']=wind-dwind
res2 = coare3(inputs)
dif3 = abs(res['hsb']-res2['hsb'])
dif3a = abs(res['hlb']-res2['hlb'])
SH_winderror = max(dif2,dif3)
LH_winderror = max(dif2a,dif3a)
print('error due to wind LAT:',LH_winderror)

inputs = {'u':wind,'t':Tair_degC,'Q':Qair*1000,'Qs':Qsea*1000,'P':pres,'ts':Tsea_degC,'lat':35.}
#print(inputs)
print('qair',Qair*1000,dQair_computed*1000)
res = coare3(inputs)
inputs['Q']=Qair*1000+dQair_computed*1000
res2 = coare3(inputs)
dif2 = abs(res['hsb']-res2['hsb'])
dif2a = abs(res['hlb']-res2['hlb'])
inputs['Q']=Qair*1000-dQair_computed*1000
res2 = coare3(inputs)
dif3 = abs(res['hsb']-res2['hsb'])
dif3a = abs(res['hlb']-res2['hlb'])
LH_Qairerror = max(dif2a,dif3a)
print('error due to Qair LAT:',LH_Qairerror)

inputs = {'u':wind,'t':Tair_degC,'Q':Qair*1000,'Qs':Qsea*1000,'P':pres,'ts':Tsea_degC,'lat':35.}
res = coare3(inputs)
inputs['Qs']=Qsea*1000+dQsea_computed*1000
res2 = coare3(inputs)
dif2 = abs(res['hsb']-res2['hsb'])
dif2a = abs(res['hlb']-res2['hlb'])
inputs['Qs']=Qsea*1000-dQsea_computed*1000
res2 = coare3(inputs)
dif3 = abs(res['hsb']-res2['hsb'])
dif3a = abs(res['hlb']-res2['hlb'])
LH_Qseaerror = max(dif2a,dif3a)
print('error due to Qsea LAT:',LH_Qseaerror)

#calculate totals
total_error = np.sqrt(LH_Qairerror**2+
                      LH_Qseaerror**2+
                      LH_winderror**2+
                      LH_perror**2)
print('Latent total error',total_error)

# -

Qair0 = air_humidity_method_qsat26air(Tair,pres,rel_hum)
Qair1 = air_humidity_method_qsat26air(Tair+dTair,pres,rel_hum)
Qair2 = air_humidity_method_qsat26air(Tair-dTair,pres,rel_hum)
Qair_dTair = max(np.abs(Qsea0-Qsea2),np.abs(Qsea0-Qsea1))
Qair0 = air_humidity_method_qsat26air(Tair,pres,rel_hum)
Qair1 = air_humidity_method_qsat26air(Tair,pres+dpres,rel_hum)
Qair2 = air_humidity_method_qsat26air(Tair,pres-dpres,rel_hum)
Qair_dpres = max(np.abs(Qsea0-Qsea2),np.abs(Qsea0-Qsea1))
Qair0 = air_humidity_method_qsat26air(Tair,pres,rel_hum)
Qair1 = air_humidity_method_qsat26air(Tair,pres,rel_hum+drh)
Qair2 = air_humidity_method_qsat26air(Tair,pres,rel_hum-drh)
Qair_drh = max(np.abs(Qsea0-Qsea2),np.abs(Qsea0-Qsea1))
dQair_computed = np.sqrt(Qair_dTair**2+Qair_dpres**2+Qair_drh**2)/1000
print('Qair and Qair uncertainty:',Qair*1000,dQair_computed*1000)


StarcoderdataJupyterScriptsDedupFiltered