.

.

.

.

.

# HealthyNumerics

HealthPoliticsEconomics | Quant Analytics | Numerics

# Computational Fluid Dynamics : Estimation of the atmospheric surface fluxes

```import datetime
from dateutil.relativedelta import relativedelta

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator

float_formatter = lambda x: "%8.2f" % x
np.set_printoptions(formatter={'float_kind':float_formatter})
```

We give here a procedure to compute the elevation and irraditaion of the sun

```def solarRadiation(day):
Latit = 47.
NumberOfDay = day.timetuple()[7]                                           # number of days since 1st January
TMST = day.timetuple()[3]*3600.0 + day.timetuple()[4]* 60.0 + day.timetuple()[5] # True mean solar time in sec
Ht = -np.pi*(43200.0-TMST)/43200.0;                                        # hourley angel
sElevation = np.arctan(sinHs/np.sqrt(1.0-sinHs*sinHs));                    # elevation of sun
sAzimut = np.cos(declin)*np.sin(Ht)/np.cos(sElevation);                    # azimuth of sun
kn = NumberOfDay * 2.0*np.pi/365.0;
I0 = 1353.0 + 45.33*np.cos(kn)   + 0.88* np.cos(2*kn) \
+ 0.005*np.cos(3*kn) + 1.8 * np.sin(kn)   \
+ 0.10 *np.sin(2*kn)+ 0.184* np.sin(3*kn)                       # solar constant
#S0:= I0 * sinHs                                                            # extraterrestrial radiation
return I0*sinHs, Ht
```

## Daily course of solar radition for some selected months

For March, June, September and December of a year we take day 21 and plot the daily corse of the solar extraterestrial raditaion

```firstDate = datetime.datetime(2017, 3,21,0,0,0)
lastDate  = datetime.datetime(2018, 3,20,0,0,0)
daySteps  = 1
monthSteps= 3
secSteps  = 0.5 * 3600

dt = datetime.timedelta(seconds = secSteps)

fig = plt.figure(figsize=(42,11))

date = firstDate
while date <= lastDate:
nextDay = date.timetuple()[2] + 1
print('date : ',date)

d0 = date.replace(hour= 0, minute= 0, second= 0)
d1 = d0.replace(day = nextDay)
Time= np.arange(d0, d1, dt).astype(datetime.datetime)

A  = np.array( [solarRadiation(Time[t]) for t in range(0,len(Time))] )
I    = np.maximum(A[:,0], 0.)   #B[B < 0] = 0
Ht = A[:,1]
ax1.plot(Ht,I,lw=20, alpha=0.2)

#date += datetime.timedelta(days=daySteps)
date += relativedelta(months= monthSteps)

plt.xlabel(r'Hour angel', fontsize=40)
plt.ylabel(r'\$[W/m^2]\$',  fontsize=40)
tickFontsize=40
plt.xticks(fontsize = tickFontsize, rotation=320);
plt.yticks(fontsize = tickFontsize);
plt.tick_params(which='major', length=20, width=10)
plt.legend(['Mar', 'Jun', 'Sep', 'Dec'], loc='upper right',prop={'size': 60})
plt.show()
```
```date :  2017-03-21 00:00:00
date :  2017-06-21 00:00:00
date :  2017-09-21 00:00:00
date :  2017-12-21 00:00:00
```

## Estimation of the atmospheric surface fluxes

`I` represents the extraterrestrial radiation at the outer bound of the earth atmosphere. So the transmission through the atmosphere still has to be included. For the moment we ommit this step. To get a first impression of the energy balance at the earth surface we use the following model, which was designed to calculate the energy balance from routine weather data.

``` A.A.M.Holtslag, A.P.VanUlden:
A simple scheme for daytime estimates of the surface fluxes from routine weather data
J Clim Appl Met 22 (1983) 517-529

A.P.VanUlden, A.A.M.Holtslag:
Estimation of atmospheric boundary layer parameters for diffusion applications
J Cli Appl Met 24 (1985) 1196-1207
```

We will use `I` as a parameter for the global raditaion, and we will simulate roughly the temperature `T` based on `I`. `T` shall ondulate between 15 and 25 degree centigrade following the form of `I` , shifted by 3 hours from noon to 3 p.m.

```T = 272. + 15. + 10. * np.roll(I , int((3.*3600.)/secSteps))/max(I)
```
```#----- Energy balance model---------------------------------
def EB2(Gl,T10):
''' Input:
z10:  measurement hight,          z0:  roughness length,
a  :  albedo,                     Gl:  global radiation,
T10:  temperature at 10 m,        N :  degree of cloudines
LEfrac: part of evaporating area
'''
#Gl = 750;        T10 = 15. + 273.15;
z10 = 10.0;        z0 = 0.5;
a   =  0.7;         N = 0.5;
LEfrac = 0.9;   sigma = 5.67E-8; #Stefan-Boltzmann constant

c1  = 5.31E-13;                c2 = 60.0;     c3 = 0.12;
T3  = np.exp(3.0*np.log(T10)); T4 = T3*T10;   T6 = T3*T3;
Qstar = ((1.0-a)*Gl + c1*T6 - sigma*T4 + c2*N)/(1+c3);

#---- ground heat flux -----------------------------------
Ag   = 1.0;     S  = np.exp(0.055*(T10-279.0));
alfa = 1.0;     Ch = 0.38*  ((1.0-alfa)*S + 1.0) / (S+1.0);
G    = Ag/(4.0*sigma*T3)*Ch*Qstar;

#---- latent heat flux ---------------------------------------
LE   =   (S/(S+1.0)*(Qstar-G) + 20.0) * LEfrac

#---- sensible heat flux  -----------------------------------
Ho   = Qstar - LE - G;
return Qstar, G, LE, Ho
```

We run the model for a mid summer day

```d0 = datetime.datetime(2018,6,21,0,0,0)
d1 = datetime.datetime(2018,6,22,0,0,0)
secSteps  = 0.5 * 3600

dt = datetime.timedelta(seconds = secSteps)
Time  = np.arange(d0, d1, dt).astype(datetime.datetime)

A  = np.array( [solarRadiation(Time[t]) for t in range(0,len(Time))] )
I  = np.maximum(A[:,0], 0.)   # A[A< 0] = 0; sets negative values to zero
T  = 272. + 15. + 10. * np.roll(I,int((3.*3600.)/secSteps))/max(I)
Qstar, G, LE, Ho = EB2(I,T)
```

#### Graphics

```fig = plt.figure(figsize=(42,11))

t = 24 * (Ht-min(Ht))/(max(Ht)-min(Ht))
ax1.plot(t,Qstar, t,2.*G, t,LE, t,Ho,lw=20, alpha=0.4)

x_majorLoc = MultipleLocator(3)
ax1.set_xlim(0, 24)
ax1.xaxis.set_major_locator(x_majorLoc)

plt.legend(['Q*', 'G', 'LE', 'Ho'], loc='upper right',prop={'size': 60})
plt.suptitle('Atmospheric Surface Layer Energy Balance', fontsize=50, fontweight='bold')
plt.title('21. June',fontsize=40)
plt.xlabel(r'Local solar time of the site', fontsize=40)
plt.ylabel(r'\$[W/m^2]\$',  fontsize=40)

tickFontsize=40
plt.xticks(fontsize = tickFontsize, rotation=320);
plt.yticks(fontsize = tickFontsize);
plt.tick_params(which='major', length=20, width=10)
plt.show()
```