<jupyter_start><jupyter_text>Idealized diurnal cycleThis notebook runs [GOTM](https://gotm.net/) with idealized diurnal heating, where the solar radiation $F_s$ is$$F_s(t) = F_0 \max\left\{\cos\left[2\pi\left(\frac{t}{86400}-0.5\right)\right], 0\right\},$$where $F_0$ is the maximum solar radiation (W m$^{-2}$) at noon and $t$ is the time of a day (s).<jupyter_code>import sys
import copy
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
sys.path.append("../../../gotmtool")
from gotmtool import *
def diurnal_cycle_factor(time):
    """Get the factor of an idealized diurnal cycle (0, 1) following
       f = max(cos(2*pi(t-0.5)), 0.)

    :time:   (datetime like) time
    :return: (numpy.array) diurnal cycle factor

    """
    tod = (3600.*time.hour+60.*time.minute+time.second)/86400.
    return np.maximum(np.cos(2.*np.pi*(tod-0.5)), 0.0)<jupyter_output><empty_output><jupyter_text>Create a modelCreate a model with environment file `../../.gotm_env.yaml`, which is created by `gotm_env_init.py`.<jupyter_code>m = Model(name='Idealized_diurnal_cycle', environ='../../.gotm_env.yaml')<jupyter_output><empty_output><jupyter_text>Take a look at what are defined in the environment file.<jupyter_code>for key in m.environ:
    print('{:>15s}: {}'.format(key, m.environ[key]) )<jupyter_output>gotmdir_code: /Users/qingli/develop/dev_gotmtool/gotm/code
   gotmdir_data: /Users/qingli/develop/dev_gotmtool/gotm/data
  gotmdir_build: /Users/qingli/develop/dev_gotmtool/gotm/build
    gotmdir_exe: /Users/qingli/develop/dev_gotmtool/gotm/exe
    gotmdir_run: /Users/qingli/develop/dev_gotmtool/gotm/run
 gotmdir_figure: /Users/qingli/develop/dev_gotmtool/gotm/figure
   gotmdir_tool: /Users/qingli/develop/dev_gotmtool/gotmtool<jupyter_text>Build the model<jupyter_code>%%time
m.build()<jupyter_output>[94mWARNING: GOTM is updated. Skipping the build step. Use 'clean=True' to rebuild[0m
CPU times: user 1.05 ms, sys: 6.44 ms, total: 7.49 ms
Wall time: 48.4 ms<jupyter_text>ConfigurationInitialize the GOTM configuration<jupyter_code>cfg = m.init_config()<jupyter_output>Generating default configuration at '/Users/qingli/develop/dev_gotmtool/gotm/run/Idealized_diurnal_cycle/gotm.yaml'...
[92mDone![0m<jupyter_text>Update the configuration<jupyter_code># setup
title = 'Idealized diurnal cycle'
nlev = 80
cfg['title'] = title
cfg['location']['name'] = 'Idealized'
cfg['location']['latitude'] = 45.0
cfg['location']['longitude'] = 0.0
cfg['location']['depth'] = 102.4
cfg['time']['start'] = '2000-01-01 00:00:00'
cfg['time']['stop']  = '2000-01-03 00:00:00'
cfg['time']['dt']    = 600.0
cfg['grid']['nlev']  = nlev

# output
cfg['output']['gotm_out']['title'] = title
cfg['output']['gotm_out']['k1_stop'] = nlev+1
cfg['output']['gotm_out']['k_stop'] = nlev
cfg['output']['gotm_out']['time_unit'] = 'dt'
cfg['output']['gotm_out']['time_step'] = 3

# forcing
datadir = m.environ['gotmdir_data']+'/examples/Idealized_diurnal_cycle'
cfg['temperature']['method'] = 'file'
cfg['temperature']['file'] = datadir+'/t_prof.dat'
cfg['salinity']['method'] = 'file'
cfg['salinity']['file'] = datadir+'/s_prof.dat'
cfg['surface']['fluxes']['heat']['method'] = 'constant'
cfg['surface']['fluxes']['heat']['constant_value'] = -159.15
cfg['surface']['fluxes']['tx']['method'] = 'constant'
cfg['surface']['fluxes']['tx']['constant_value'] = 0.1
cfg['surface']['fluxes']['ty']['method'] = 'constant'
cfg['surface']['fluxes']['ty']['constant_value'] = 0.0
cfg['surface']['swr']['method'] = 'file'
cfg['surface']['swr']['file'] = 'swr.dat'
cfg['surface']['precip']['method'] = 'constant'
cfg['waves']['stokes_drift']['us']['method'] = 'off'
cfg['waves']['stokes_drift']['vs']['method'] = 'off'

# EOS -- use linear
cfg['eq_state']['form'] = 'linear_custom'
cfg['eq_state']['linear']['T0'] = 20.0
cfg['eq_state']['linear']['S0'] = 35.0
cfg['eq_state']['linear']['dtr0'] = -0.2
cfg['eq_state']['linear']['dsr0'] = 0.8
cfg['physical_constants']['rho_0'] = 1000.0

# water type (Jerlov IB)
cfg['light_extinction']['method'] = 'jerlov-ib'
# solar radiation
F0 = 500.
time = pd.date_range(cfg['time']['start'], cfg['time']['stop'], freq='30min')
factor = diurnal_cycle_factor(time)
swr = F0*factor<jupyter_output><empty_output><jupyter_text>We want to run the case using two turbulence closure methods:- The generic length scale (GLS; [Umlauf and Burchard, 2003](https://doi.org/10.1357/002224003322005087)) model in the $k$-$\epsilon$ formulation with the weak-equilibrium stability function by [Canuto et al., 2001](https://doi.org/10.1175/1520-0485(2001)031%3C1413:OTPIOP%3E2.0.CO;2) (C01A).- The K-profile parameterization (KPP; [Large et al., 1994](https://doi.org/10.1029/94RG01872), [Griffies et al., 2015](https://github.com/CVMix/CVMix-description/raw/master/cvmix.pdf)) with [CVMix](http://cvmix.github.io).<jupyter_code># collect the configurations and the labels of the two runs
cfgs = []
labels = []
cfg['turbulence']['turb_method'] = 'second_order'
cfg['turbulence']['tke_method'] = 'tke'
cfg['turbulence']['len_scale_method'] = 'gls'
cfg['turbulence']['scnd']['method'] =  'weak_eq_kb_eq'
cfg['turbulence']['scnd']['scnd_coeff'] =  'canuto-a'
cfg['turbulence']['turb_param']['length_lim'] = 'false'
cfg['turbulence']['turb_param']['compute_c3'] = 'true'
cfg['turbulence']['turb_param']['Ri_st'] = 0.25
cfg['turbulence']['generic']['gen_m'] = 1.5 
cfg['turbulence']['generic']['gen_n'] = -1.0
cfg['turbulence']['generic']['gen_p'] = 3.0 
cfg['turbulence']['generic']['cpsi1'] = 1.44
cfg['turbulence']['generic']['cpsi2'] = 1.92
cfg['turbulence']['generic']['cpsi3minus'] = -0.63
cfg['turbulence']['generic']['cpsi3plus'] = 1.0 
cfg['turbulence']['generic']['sig_kpsi'] = 1.0 
cfg['turbulence']['generic']['sig_psi'] = 1.3
cfgs.append(copy.deepcopy(cfg))
labels.append('GLS-C01A')
cfg['turbulence']['turb_method'] = 'cvmix'
cfg['cvmix']['surface_layer']['kpp']['langmuir_method'] = 'none'
cfgs.append(copy.deepcopy(cfg))
labels.append('KPP-CVMix')
# save solar radiation file
for lb in labels:
    run_dir = m.environ['gotmdir_run']+'/'+m.name+'/'+lb
    os.makedirs(run_dir, exist_ok=True)
    dat_dump_ts(time, [swr], run_dir+'/swr.dat')<jupyter_output><empty_output><jupyter_text>Run the model<jupyter_code>%%time
sims = m.run_batch(configs=cfgs, labels=labels, nproc=2)<jupyter_output>CPU times: user 5.33 ms, sys: 8.81 ms, total: 14.1 ms
Wall time: 327 ms<jupyter_text>ResultsLoad the data into an `xarray.Dataset`.<jupyter_code>fig, axarr = plt.subplots(2, sharex='col')
fig.set_size_inches(8, 8)
data0 = sims[0].load_data()
levels = np.linspace(19.6, 20, 41)
data0.temp.plot(ax=axarr[0], levels=levels)
axarr[0].set_title(labels[0])
axarr[0].set_ylim([-40, 0])
data1 = sims[1].load_data()
data1.temp.plot(ax=axarr[1], levels=levels)
axarr[1].set_title(labels[1])
axarr[1].set_ylim([-40, 0])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.21131909915179878, 'fraction_numerical': 0.04314126937701082, 'mean_word_length': 6.851894374282434, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 34, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '443994', 'license': 'MIT', 'n_tokens_mistral': 2842, 'n_tokens_neox': 2616, 'n_words': 442, 'path': 'examples/Idealized_diurnal_cycle/run_idealized_diurnal_cycle.ipynb', 'repo_name': 'jithuraju1290/gotmtool', 'size': '6838'}