<jupyter_start><jupyter_text>Generate the Wadden Sea forcing data file<jupyter_code>%matplotlib inline 

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import xarray as xr
import src.utilities as ut<jupyter_output><empty_output><jupyter_text>Read a World Ocean Database (WOD) nc file (rectangle near to the eastern part of the Wadden Sea) which is produced by Ocean Data Viewer (ODV) and smooth through years for depths less then 11 m, var1 is depth:<jupyter_code>ds = xr.open_dataset('data/wadden_sea.nc')
df = ds.to_dataframe()
df = df[(df.date_time > '2000-01-01 00:00:00.0') & (df.var1 < 11)]
df = df[['longitude', 'latitude', 'date_time'
        ,'var1', 'var2', 'var3', 'var4', 'var5', 'var6', 'var7'
        ,'var8', 'var9', 'var10', 'var11', 'var12', 'var13', 'var14']]
df.columns = ['longitude', 'latitude', 'date_time'
             ,'Depth', 'Temperature', 'Salinity', 'Oxygen', 'Phosphate', 'Silicate', 'Nitrate'
             ,'Nitrite', 'pH', 'Chlorophyll a', 'Plankton biomass', 'Alkalinity', 'NO2+NO3', 'pCO2']
df['Oxygen'] = df['Oxygen']*44.661 # 1 ml/l = 10^3/22.391 = 44.661 μmol/l
df['Day of year'] = df.date_time.dt.dayofyear.values<jupyter_output><empty_output><jupyter_text>Make an average year from a range of 11 years:<jupyter_code>dfyear = df[(df.date_time > '2000-01-01 00:00:00.0') & (df.date_time < '2011-01-01 00:00:00.0')]
dfyear = dfyear.groupby(['Day of year']).mean()
dfyear = dfyear.reindex(index=range(1,366))<jupyter_output><empty_output><jupyter_text>Interpolate:<jupyter_code>dfyear = dfyear.interpolate(method='linear',limit_direction='forward')
dfyear.loc[0:16,:] = dfyear[-1:].values
dfyear.loc[0:16,'Chlorophyll a'] = 0
dfyear = dfyear.interpolate(method='linear',limit_direction='both')
dfyear = dfyear.rolling(window=16).mean()
dfyear = dfyear.interpolate(method='linear',limit_direction='both')<jupyter_output><empty_output><jupyter_text>Make three successive years from one;other parameters also will be produced for 3 succesive similar years to construct a nc file with input data for a fortran model that resolves water and sediments domains:<jupyter_code>df3years = dfyear.append([dfyear, dfyear])
df3years.index = range(1, len(df3years)+1)
df3years.index.name = "day"<jupyter_output><empty_output><jupyter_text>Chl a data is from Loebl et al. (2007) since the data from WOD is very controversial.We do not need high precision here so we can take some points from the graph from the article and then interpolate.<jupyter_code>chl_array = np.empty(365)
chl_array[:] = np.nan
chl_array[0] = 2; chl_array[15] = 2; chl_array[80] = 3; chl_array[110] = 15;
chl_array[130] = 25; chl_array[170] = 2; chl_array[180] = 2; chl_array[220] = 2;
chl_array[250] = 2; chl_array[300] = 2; chl_array[364] = 2
dfchl = pd.DataFrame({'Chl_a': chl_array})
dfchl = dfchl.interpolate(method='linear',limit_direction='both')
dfchl = dfchl.rolling(window=16).mean()
dfchl = dfchl.interpolate(method='linear',limit_direction='both')
dfchl3years = dfchl.append([dfchl, dfchl])<jupyter_output><empty_output><jupyter_text>Read and process a text file downloaded from http://tsdata.hzg.de/ (Pile Hoernum1 station) with photosynthetically active radiation (PAR) data.We do it here only to check if our theoretical function to generate PAR is legit.<jupyter_code>dfpar = pd.read_csv('data/Pile_Hoernum1_PAR_20080101_20180608.dat', header=None, delim_whitespace=True
                   , skiprows=22)
dfpar[0] = dfpar[0] + ' ' + dfpar[1]
dfpar = dfpar.drop(1, axis=1)
dfpar[0] =  pd.to_datetime(dfpar[0])
dfpar = dfpar.rename(columns={0: 'date_time', 2: 'par', 3: 'flag'})
dfpar = dfpar[(dfpar['flag'] < 3) & (dfpar['par'] > 0)]
dfpar = dfpar.drop('flag', axis=1)
dfpar['Day of year'] = dfpar.date_time.dt.dayofyear.values
dfparmean = dfpar.groupby(['Day of year']).mean()
dfparmean = dfparmean.reindex(index=range(1,366))
dfparmean = dfparmean.interpolate(method='linear',limit_direction='both')
dfparmean = dfparmean.rolling(window=16).mean()
dfparmean = dfparmean.interpolate(method='linear',limit_direction='both')
par_generated = np.zeros(365)
par_generated[:] = [ut.surface_radiation(x, 54) for x in range(365)]
par3years_generated = np.concatenate((par_generated, par_generated, par_generated), axis=None)
# visualize measured PAR
name = 'par'; plt.plot(dfparmean.index, dfparmean[name]);
plt.gcf().set_size_inches(10,2); plt.title(name); plt.xlabel('Day'); plt.ylabel('$\mu M\,s^{-1}\,m^{-2}$');
# visualize calculated PAR
name = 'par'; plt.plot(dfparmean.index, par_generated);
plt.gcf().set_size_inches(10,2); plt.title(name); plt.xlabel('Day'); plt.ylabel('$\mu M\,s^{-1}\,m^{-2}$');<jupyter_output><empty_output><jupyter_text>Generate yearly organic matter (OM) influx according to Beusekom (1999):<jupyter_code>from scipy import integrate<jupyter_output><empty_output><jupyter_text>The original flux of OM is 110 $gC\,m^{-2} y^{-1}$ for the total water column, the Wadden sea has an average depth 2.5 meters, so concentration is 110/2.5 for $m^2$.<jupyter_code>annual_flux_om = 110/2.5 * 1000 # mg C m-2 y-1<jupyter_output><empty_output><jupyter_text>Function for parameter maximum 'm' identification with given total yearly value according to simplification that OM influx has an approximately sinusoidal yearly profile:<jupyter_code>def sinusoidal(fr, ti, st, var):
    """fr - lower bound of the parameter
       ti - upper bound of the parameter
       st - step for the parameter selection
       var - the value of the assessed sinusoidal integral"""
    for m in np.arange(fr, ti, st):
        flux = lambda day: m/2 + (m/3)*(1+np.sin(2*np.pi*(day-20)/365))
        sumflux = integrate.quad(flux, 0, 365)
        if var-1 < sumflux[0] < var+1:
            print(m)
sinusoidal(100,1000,0.01,annual_flux_om) # finds a maximum value of yearly OM influx
m = 144; day = np.arange(0,1095,1)
# use sinusoidal function to produce yearly profile of OM influx
y_om = m/2 + (m/3)*(1+np.sin(2*np.pi*(day-20)/365))<jupyter_output><empty_output><jupyter_text>NH4 is from Beusekom et al. (2009).Again we do not need precise data here.<jupyter_code>nh4year = np.zeros(365)
nh4year[0:59] = 7; nh4year[59:90] = 4; nh4year[90:240] = 0
nh4year[240:270] = 2; nh4year[270:305] = 5; nh4year[305:365] = 8
nh4year_beusekom = np.concatenate((nh4year, nh4year, nh4year), axis=None)<jupyter_output><empty_output><jupyter_text>Generate a wadden_sea_out.nc file which consist of 3 successive similar years:<jupyter_code>import itertools as it
from netCDF4 import Dataset

rootgrp = Dataset("data/wadden_sea_out.nc", "w", format='NETCDF3_CLASSIC')

levelface = rootgrp.createDimension("levelface", 3)
levelcntr = rootgrp.createDimension("levelcntr", 2)
time = rootgrp.createDimension("time", 1095)

times = rootgrp.createVariable("time","f8",("time",))
times.units = "seconds since 2010-01-01 00:00:00.0"
times.calendar = "gregorian"

par = rootgrp.createVariable("par","f4",("time",))
par.units = "microM photons m−2 s−1"

hice = rootgrp.createVariable("hice","f4",("time",))
snowthick = rootgrp.createVariable("snowthick","f4",("time",))
icesurft = rootgrp.createVariable("icesurft","f4",("time",))

levelsface = rootgrp.createVariable("levelface","f4",("levelface",))
levelsface.units = "m"
levelscntr = rootgrp.createVariable("levelcntr","f4",("levelcntr",))
levelscntr.units = "m"
temperature = rootgrp.createVariable("temperature","f8",("time","levelcntr",))
temperature.units = "C degrees"
salinity = rootgrp.createVariable("salinity","f8",("time","levelcntr",))
salinity.units = "psu"
sigma = rootgrp.createVariable("sigma","f8",("time","levelcntr",))
sigma.units = "kg m-1"
turbulence = rootgrp.createVariable("turbulence","f8",("time","levelface",))
turbulence.units = "m2 s-1"

ammonium = rootgrp.createVariable("ammonium","f8",("time","levelcntr",))
ammonium.units = "mM N m-3"
nitrate = rootgrp.createVariable("nitrate","f8",("time","levelcntr",))
nitrate.units = "mM N m-3"
sulfate = rootgrp.createVariable("sulfate","f8",("time","levelcntr",))
nitrate.units = "mM S m-3"
silicate = rootgrp.createVariable("silicate","f8",("time","levelcntr",))
silicate.units = "mM Si m-3"
phosphate = rootgrp.createVariable("phosphate","f8",("time","levelcntr",))
phosphate.units = "mM P m-3"
oxygen = rootgrp.createVariable("oxygen","f8",("time","levelcntr",))
oxygen.units = "mM O2 m-3"
om = rootgrp.createVariable("om_flux","f8",("time","levelcntr",))
om.units = "mg C m-2 s-1"
chl_a = rootgrp.createVariable("chl_a","f8",("time","levelcntr",))
chl_a.units = "mg Chl a"

# lvlface = np.array([0.5, 0.25, 0])
# lvlcntr = np.array([0.375, 0.125])
# lvlface = np.array([1.5, 0.75, 0])
# lvlcntr = np.array([1.125, 0.375])
# lvlface = np.array([2.0, 1.0, 0])
# lvlcntr = np.array([1.5, 0.5])
lvlface = np.array([2.5, 1.25, 0])
lvlcntr = np.array([1.875, 0.625])
# lvlface = np.array([3.5, 1.75, 0])
# lvlcntr = np.array([2.625, 0.875])

levelsface[:] = lvlface
levelscntr[:] = lvlcntr

times[:] = df3years.index.values * 24 * 60 * 60
par[:] = par3years_generated
hice[:] = 0
snowthick[:] = 0
icesurft[:] = 0
temperature[:,:] = np.vstack((df3years.Temperature.values, df3years.Temperature.values)).T
salinity[:,:] = np.vstack((df3years.Salinity.values, df3years.Salinity.values)).T<jupyter_output><empty_output><jupyter_text>Calculate eddy diffusivity, see utilities.py:<jupyter_code>sigma[:,0] = [ut.svan(s,t,p) for s,t,p in 
              zip(df3years.Salinity.values, df3years.Temperature.values, it.count(lvlcntr[0]))]
sigma[:,1] = [ut.svan(s,t,p) for s,t,p in 
              zip(df3years.Salinity.values, df3years.Temperature.values, it.count(lvlcntr[1]))]
sigmaface = np.zeros((1095,3))
sigmaface[:,0] = [ut.svan(s,t,p) for s,t,p in 
              zip(df3years.Salinity.values, df3years.Temperature.values, it.count(lvlface[0]))]
sigmaface[:,1] = [ut.svan(s,t,p) for s,t,p in 
              zip(df3years.Salinity.values, df3years.Temperature.values, it.count(lvlface[1]))]
sigmaface[:,2] = [ut.svan(s,t,p) for s,t,p in 
              zip(df3years.Salinity.values, df3years.Temperature.values, it.count(lvlface[2]))]
turbulence[:,0] = ut.gargett(sigmaface[:,0], 1)
turbulence[:,1] = ut.gargett(sigmaface[:,1], 1)
turbulence[:,2] = ut.gargett(sigmaface[:,2], 1)
ammonium[:,:] = np.vstack((nh4year_beusekom, nh4year_beusekom)).T
sulfate[:,:] = 25000
nitrate[:,:] = np.vstack((df3years.Nitrate.values, df3years.Nitrate.values)).T
silicate[:,:] = np.vstack((df3years.Silicate.values, df3years.Silicate.values)).T
phosphate[:,:] = np.vstack((df3years.Phosphate.values, df3years.Phosphate.values)).T
oxygen[:,:] = np.vstack((df3years.Oxygen.values, df3years.Oxygen.values)).T
chl_a[:,:] = np.vstack((dfchl3years['Chl_a'].values, dfchl3years['Chl_a'].values)).T
om[:,:] = np.vstack((y_om/86400, y_om/86400)).T
rootgrp.close()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1765797531780644, 'fraction_numerical': 0.05103461074510532, 'mean_word_length': 6.003248862897986, 'pattern_counts': {'":': 0, '<': 58, '<?xml version=': 0, '>': 56, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '321247', 'license': 'CC-BY-3.0', 'n_tokens_mistral': 4531, 'n_tokens_neox': 4166, 'n_words': 823, 'path': 's_1_generate_netcdf.ipynb', 'repo_name': 'limash/ws_notebook', 'size': '10777'}