forked from ClimateGlobalChange/DCMIP2016
-
Notifications
You must be signed in to change notification settings - Fork 0
/
DCMIP2016-TestCaseDocument_v1.tex
1592 lines (1388 loc) · 102 KB
/
DCMIP2016-TestCaseDocument_v1.tex
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
\documentclass[times,doublespace]{fldauth}
\usepackage{hyperref}
\usepackage{lineno}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}
\newcommand{\vb}{\mathbf}
\newcommand{\vg}{\boldsymbol}
\newcommand{\diff}[2]{\frac{d #1}{d #2}}
\newcommand{\pdiff}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\avg}[1]{\overline{#1}}
\newcommand{\dblavg}[1]{\overline{\overline{#1}}}
\newcommand{\mat}[1]{\mathbf{\mathsf{#1}}}
\newcommand{\arccot}{\mathrm{arccot}}
\newcommand{\rbd}[1]{\raisebox{-1.5ex}[1.5ex]{#1}}
\newcommand{\rbu}[1]{\raisebox{0.5ex}[-0.5ex]{#1}}
\newcommand{\rb}[1]{\raisebox{2ex}[-2ex]{#1}}
\newcommand\T{\rule{0pt}{2.6ex}}
\newcommand\B{\rule[-1.2ex]{0pt}{0pt}}
\begin{document}
\setcounter{section}{-1}
\title{Dynamical Core Model Intercomparison Project (DCMIP2016) \\
Test Case Document}
\author{Paul Aaron Ullrich \\ Christiane Jablonowski \\ Kevin A. Reed \\ Colin Zarzycki \\ Peter Hjort Lauritzen \\ Ramachandran D. Nair \\ James Kent \\Antonin Verlet-Banide \\ \vspace{3cm} DCMIP Summer School: June 2016}
\maketitle
\begin{center}
Version 0.5 \\
\today
\end{center}
\vspace{2cm}
\begin{center}
Email questions to \textbf{[email protected]}
\end{center}
\clearpage
\section*{Overview of the Test Case Suite}
The set of test cases collected in this document has been developed for the 2016 Dynamical Core Model Intercomparison and (DCMIP2016) in an effort to understand the broad treatment of the equations of motion within a variety of atmospheric General Circulation Models (GCMs). In contrast with the test cases described in the DCMIP2012 and DCMIP2008 test case documents, the suite of test cases that are described here focus on issues related to physics-dynamics coupling, non-hydrostatic scales and modeling systems that support variable resolution. To better isolate these features, the majority of tests described in this document do not include topographic forcing and emphasize features that have a relatively small spatial footprint relative to the global simulation domain. The tests proposed here have been drawn from the recent literature on global model intercomparison testing. To support the integration of these tests into existing models, a collection of initialization routines have been provided that can be used for very quickly setting up the initial conditions. Augmented Kessler physics routines have also been provided to support simple moisture feedbacks, a boundary layer parameterization and surface friction. If models provide both a shallow-atmosphere and deep-atmosphere configuration, we recommend the shallow-atmosphere setup to avoid imbalances in the initial conditions. It is further expected that most models will be run under a non-hydrostatic configuration, which will permit a correct solution to the supercell test. If models can be configured both as a hydrostatic and non-hydrostatic dynamical core, additional hydrostatic simulations might be conducted to evaluate the direct impact of the hydrostatic approximation. Table \ref{tab:TestCases} provides an overview of all test cases described in this document.
\begin{table}[h]
\caption{A list of test cases described in this document.} \label{tab:TestCases}
%\ \\
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}ll}
\hline Test Case \T \B& Description \\
\hline \multicolumn{2}{l}{Tier 1 Test Cases (Required)} \T \B \\
\hline
1-1 \T & Moist baroclinic wave with terminator chemistry \\
1-2 \T & Idealized tropical cyclone \\
1-3 \T & Supercell \\
\hline \multicolumn{2}{l}{Tier 2 Test Cases (Optional)} \T \B \\ \hline
2-1 \T \B & Moist Held-Suarez test \\ \hline
\hline
\end{tabular*}
\end{table}
\section{Practical Considerations}
\subsection{List of Symbols}
Throughout this test case document we will use $\lambda \in [0, 2 \pi)$ to denote longitude, $\varphi \in [-\pi/2, \pi/2]$ to represent latitude, $z$ to represent the height with respect to the mean sea level (assumed to be zero), and $p$ to symbolize the pressure. Table \ref{tab:symbols} lists the symbols used for the initialization of the model.
\begin{table}[h]
\caption{List of symbols for the model initialization} \label{tab:symbols}
\begin{center}
\begin{tabular}{cl}
\hline Symbol & Description \\ \hline
$\lambda$ & Longitude (in radians) \\
$\varphi$ & Latitude (in radians) \\
$z$ & Height with respect to mean sea level (set to zero) \\
$p_s$ & Surface pressure ($p_s$ of moist air if $q>0$) \\
$\Phi_s$ & Surface geopotential \\
$z_s$ & Surface elevation with respect to mean sea level (set to zero) \\
$u$ & Zonal wind \\
$v$ & Meridional wind \\
$w$ & Vertical velocity \\
$\omega$ & Vertical pressure velocity \\
$\delta$ & Divergence\\
$\zeta$ & Relative vorticity\\
$p$ & Pressure (pressure of moist air if $q>0$) \\
$\rho$ & Density (density of moist air if $q>0$)\\
$T$ &Temperature \\
$T_v$ & Virtual temperature \\
$\Theta$ & Potential temperature \\
$\Theta_v$ & Virtual potential temperature \\
$q$ & Specific humidity \\
$P_{ls}$ & Large-scale precipitation rate \\
$q_c$ & Cloud water mixing ratio \\
$q_r$ & Rain water mixing ratio \\
$q_{Cl}$ & Singlet chlorine mixing ratio \\
$q_{Cl2}$ & Chlorine gas mixing ratio \\
\hline
\end{tabular}
\end{center}
\end{table}
\subsection{List of Physical Constants}
A list of physical constants which are used throughout this document is given in Table \ref{tab:PhysicalConstants}. Constants which are specific to each test case are similarly tabulated at the beginning of each section.
\begin{table}[h]
\caption{A list of physical constants used in this document.} \label{tab:PhysicalConstants}
%\ \\
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}lll}
\hline Constant & Description & Value \\
\hline $a_{\tiny \mbox{ref}}$ & Radius of the Earth & $6.37122 \times 10^{6}\ \mbox{m}$ \\
$\Omega_{\tiny \mbox{ref}}$ & Rotational speed of the Earth & $7.292\ \times 10^{-5}\ \mbox{s}^{-1}$ \\
$X$ & Reduced-size Earth reduction factor & variable (default $= 1$) \\
$a$ & Scaled radius of the Earth & $a_{\tiny \mbox{ref}} / X$ \\
$\Omega$ & Scaled rotational speed of the Earth & $\Omega_{\tiny \mbox{ref}} \cdot X$ \\
$g$ & Gravity & $9.80616\ \mbox{m}\ \mbox{s}^{-2}$ \\
$p_0$ & Reference pressure & $1000\ \mbox{hPa}$ \\
$c_p$ & Specific heat capacity of dry air at constant pressure & $1004.5\ \mbox{J}\ \mbox{kg}^{-1}\ \mbox{K}^{-1}$ \\
$c_v$ & Specific heat capacity of dry air at constant volume & $717.5\ \mbox{J}\ \mbox{kg}^{-1}\ \mbox{K}^{-1}$ \\
$R_d$ & Gas constant for dry air & $287.0\ \mbox{J}\ \mbox{kg}^{-1}\ \mbox{K}^{-1}$ \\
$R_\nu$ & Gas constant for water vapor & $461.5$ J kg$^{-1}$ K$^{-1}$ \\
$\kappa$ & Ratio of $R_d$ to $c_p$ & $2/7$ \\
$\varepsilon$ & Ratio of $R_d$ to $R_\nu$ & $0.622$ \\
$M_v$ & Constant for virtual temperature conversion & $0.608$ \\
$\rho_{water}$ & Density of water & 1000 kg m$^{-3}$ \\
\hline
\end{tabular*}
\end{table}
\subsection{Great Circle Distance}
The great circle distance is used throughout the document and is given by
\begin{equation}
R_c(\lambda_1, \varphi_1; \lambda_2, \varphi_2) = a \arccos \left( \sin \varphi_1 \sin \varphi_2 + \cos \varphi_1 \cos \varphi_2 \cos (\lambda_1 - \lambda_2) \right).
\end{equation}
\subsection{Choice of Prognostic Variables and Equation of State}
It is assumed throughout this document that the moist ideal gas law is satisfied,
\begin{equation} \label{eq:idealgaslaw}
p = \rho R_d T_v,
\end{equation} where the virtual temperature is defined in terms of the temperature as
\begin{equation} \label{eq:virtualtemperature}
T_v = T (1 + M_v q).
\end{equation} The virtual potential temperature is similarly related to the potential temperature via the relationship
\begin{equation}
\theta_v = \theta (1 + M_v q).
\end{equation}
\subsection{Small-Planet Experiments}
The test case suite makes use of small-planet experiments that have the potential to expose the differences between hydrostatic and non-hydrostatic modeling approaches at reasonable computational cost. In particular, the small-planet setups allow the evaluation of the model behavior with physical grid spacings down to a few hundred meters. In some instances, we suggest small planets with circumferences of about 40 km and a vertical extent of 30 km which raises questions concerning the validity of the shallow-atmosphere approximation. However, since the experiments are not compared to observations, we still ask for the use of the shallow-atmosphere approach to allow for intercomparisons among the DCMIP models and to avoid imbalances of the initial conditions.
When a non-unity reduction factor $X$ is applied in order to shrink the size of the Earth and thereby the physical grid spacing of the computational grid, a variety of model adjustments become necessary.
Most prominently these include the scaling of the radius, the rotational speed, the model time step and explicit viscosity parameters (if applied). The adjustment steps for small-planet simulations are:
\begin{itemize}
\item Divide the radius of the Earth $a_{\tiny \mbox{ref}}$ by $X$ to obtain the rescaled radius $a = \frac{a_{\tiny \mbox{ref}}}{X}$.
\item Divide the length of the dynamics time step $\Delta t$ by $X$, especially if a CFL condition needs to be obeyed.
\item In case of rotating planets: multiply the Earth's angular velocity $\Omega_{\tiny \mbox{ref}}$ by the factor $X$ to obtain the rescaled angular speed $\Omega = \Omega_{\tiny \mbox{ref}}X$. This guarantees that the characteristics of Rossby waves are comparable in unscaled and scaled model experiments since the Rossby number stays constant.
\item In case of explicit diffusion of type $K_{2k}\nabla^{2k}$ with a prescribed diffusion coefficient $K_{2k}$ (and $k=1,2\ldots$) divide $K_{2k}$ by the factor $X^{2k-1}$. This accounts for a reduction of the e-folding time $\tau$ and the horizontal grid spacing $\Delta x$ according to the relationship $\frac{(\Delta x)^{2k}/X^{2k}}{\tau/X}$. The $K_{2k}$ diffusion coefficient is typically based on such a relationship. Note that some models might provide an automatic scaling of the diffusion coefficients according to the actual dynamics time step and grid spacing. If a model applies Rayleigh friction as a sponge near the model top, the friction coefficient needs to be multiplied by X. Again, this corresponds to a reduction of the e-folding frictional time scale $\frac{1}{\tau/X} = \frac{X}{\tau}$ in small-planet experiments.
\item If physical forcing mechanisms are present on the right hand side of the equations of motion the strengths of the physical forcing must be increased (multiplied) by the factor X. This will not be applicable to the test cases presented here unless Rayleigh friction is considered as such a forcing as outlined above.
\item Any deviations from this prescription should be documented. For example, some modeling groups may have a prescribed diffusion coefficient that depends on grid resolution that can be used in place of the scaled diffusivity described above.
\end{itemize}
\subsection{Notes on the Requested Model Output}
\label{sec:notes_output}
\subsubsection{NetCDF}
A fundamental requirement for the exchange of scientific data is the ability to precisely describe the physical quantities being represented. Therefore, special attention needs to be paid to the representation of the model data in the output files. We require data in the `Network Common Data Form' (netCDF) \cite{netcdf} that adhere to the netCDF Climate and Forecast (CF) metadata convention (if possible to version 1.6 from Dec. 2011 \cite{netcdf-cf}). All netCDF files should have the file name extension `.nc'.
We specify details of the netCDF requirements in Appendix~\ref{sec:netcdf} and will also provide help via NCO operators before the DCMIP event to make the data CF-compliant if necessary. Please communicate your output constraints or concerns (if any) to the DCMIP organizers as soon as possible.
\subsubsection{Naming convention for output file names}
The naming convention for the final output file names is:
\vspace{0.2cm}
\noindent
model.experiment\_id.horizontal\_resolution.levels.grid.equation.description.variable.nc
\vspace{0.2cm}
\noindent
NetCDF-CF compliant files with these naming conventions will be uploaded to the Earth System Grid Federation.
The optional free-text character string at the end of the file name (``description'') might be used to denote a special setting, such as a varied non-default diffusion coefficient or the selection of a variable-resolution configuration.
The standard keywords for the DCMIP models are listed in the left columns of Tables \ref{tab:keywords} and \ref{tab:keywordsgrid}. The keyword for the number of vertical levels, e.g. ``L30'', denotes the number of the full model levels (here 30). In case a model also employs interface levels, this setting indirectly implies that the model has one more interface level. Examples are \\
\vspace{0.0cm}
\noindent
cam-se.161.r200.L60.latlon.hydro.4th\_order\_div\_damping.U.nc \\
tempest.162.r100.L120.interp\_latlon.nonhydro.V.nc \\
fv3.162.r50.L30.cubed.nonhydro.Theta.nc \\
cam-se.162.r50.L30.interp\_latlon.hydro.double\_K4\_2e15.T.nc \\
cam-se.163.r25.L30.interp\_latlon.hydro.variable\_220\_25km.PS.nc \\
The last example expresses the suggested naming convention for a variable-resolution simulation. The keyword `r25' denotes the finest grid spacing in the variable-resolution run which is, in this example, interpolated to a regular latitude-longitude grid with a fine grid spacing of about 28 km. The free-text `description' contains the suggested keyword 'variable' and specifies details about the resolution range in the variable-resolution run. In addition to such an interpolated data set, the original output file of a variable-resolution simulation should also be archived. Its name should be `cam-se.163.r25.L30.cubed.hydro.variable\_220\_25km.nc' for the example above, where the specification `cubed' for the underlying base grid has been chosen.
\begin{table}[h]
\caption{Keywords for the output file naming convention.} \label{tab:keywords}
%\ \\
\begin{tabular*}{\textwidth}{ll}
%\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}ll}
\hline
model \hspace{2cm}\T \B & Organization \\ \hline
cam-se & DOE/SNL \& NCAR \\
chombo & DOE/LBNL \& NCAR \\
csu \B & Colorado State University \\ \hline
dynamico & IPSL \\
endgame & UK Met Office\\
gem & Environment Canada \\
homme-nh & DOE/SNL \& CU Boulder \\
fv3 & NOAA GFDL \\
icon & MPI \& DWD \\
ifs-fv & ECMWF \\
neptune & NRL \& NPS \\
nicam & RIKEN \& University of Tokyo \\
mpas & NCAR \& Los Alamos National Laboratory \\
olam & University of Miami \\
tempest & University of California, Davis \\
\hline
\end{tabular*}
\end{table}
\begin{table}[h]
\caption{Continuation of the list of keywords for the output file naming convention: Grid specifications.} \label{tab:keywordsgrid}
%\ \\
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}ll}
\hline
experiment\_id \T \B & Description of the test case number \\ \hline
161 \T& Moist baroclinic wave \\
161-preciponly \T& Moist baroclinic wave w/ no PBL or surface fluxes \\
162-rjpbl & Idealized tropical cyclone w/ Reed-Jablonowski PBL \\
162-bryanpbl & Idealized tropical cyclone w/ Bryan PBL \\
163 & Supercell \\
\hline horizontal\_resolution \T \B & Approximate grid spacing in degrees or equatorial grid spacing (unscaled)\\ \hline
r200 \T & between 2$^\circ$-3.0$^\circ\quad\;\;\,\,$ or 200-300 km (on unscaled planet)\\
r100 & between 1$^\circ$-1.5$^\circ\quad\;\;\,\,$ or 100-150 km (on unscaled planet)\\
r50 & between 0.5$^\circ$-1.0$^\circ \quad\,$ or 50-100 km (on unscaled planet)\\
r25 \B& between 0.25$^\circ$-0.5$^\circ\,$ or 25-50 km (on unscaled planet)\\ \hline
levels \T \B & Number of vertical levels (full model levels)\\ \hline
L10 \T& 10 levels \\
L15 & 15 levels \\
L20 & 20 levels \\
L30 & 30 levels \\
L60 & 60 levels \\
L120 \B & 120 levels \\
\hline
grid \T \B& Type of native grid or indicator of interpolated grid \\ \hline
cubed \T & cubed sphere \\
hex & hexagonal grid based on an icosahedron, maybe optimized via spring dynamics\\
gauss & Gaussian grid \\
tri & triangular grid based on an icosahedron, maybe optimized via spring dynamics\\
interp\_latlon & interpolated latitude-longitude grid \\
latlon & latitude-longitude \\
oct & octagonal \\
red\_latlon & reduced latitude-longitude grid \\
red\_gauss & reduced Gaussian grid \\
voronoi & spherical centriodal Voronoi tessellation \\
yinyang \B & Yin-Yang \\ \hline
equation \T \B & Indicator for hydrostatic or non-hydrostatic simulation \\ \hline
hydro \T & hydrostatic \\
nonhydro \B & non-hydrostatic \\
\hline
\end{tabular*}
\end{table}
\subsection{Short Note on Data Analysis and Visualization}
We will provide NCAR Command Language (NCL) scripts to help visualize the model results and provide analysis functions. In addition, the DCMIP participants will have access to interactive Graphical User Interfaces (GUIs) to support the visualization and model intercomparison. Among the GUIs are the netCDF viewers Ncview and Panoply which are public domain tools and locally installed on the NCAR mirage server. In addition, we expect to provide some basic online visualization capabilities via NOAA's Live Access Server (LAS) software. A key to the successful visualization is the adherence of the output data sets to the netCDF-CF standard (see Appendix \ref{sec:netcdf}).
\subsubsection{Computational grid}
Most DCMIP models utilize non-orthogonal computational grids like cubed-sphere grids, icosahedral grids, hexagonal grid, Voronoi grids or Yin-Yang grids. Among the mix of DCMIP GCMs are even models that provide provisions for variable-resolution grids. We encourage the use of the variable-resolution configurations as an additional test option whenever possible.
This raises questions concerning the desired representation of the data in the netCDF output files.
If models on non-traditional (non latitude-longitude) grids are used, we ask for two output files that represent the identical model run. The first output file should be written on the native computational grid without any interpolations. In addition, most models will likely provide built-in provisions for interpolated output to a regular (equidistant in degrees) latitude-longitude grid. We therefore also ask for a second output file that represents the data on model levels on the interpolated latitude-longitude grid. We ask for co-located (Arakawa-A type) data on the interpolated grid regardless of the GCM's staggering options. The grid spacing of the interpolated grid should be comparable to the actual resolution of the model run, which might be for example $1^\circ \times 1^\circ$. Using this example, the interpolated grid will have $180 \times 360$ horizontal grid points if the equator and pole points are not part of the interpolated grid (equally spaced in the longitudinal and latitudinal directions). If the equator and pole points are included it yields $181 \times 360$ horizontal grid points. The $180 \times 360$ configuration is preferred for model intercomparison purposes.
Models on regular latitude-longitude or Gaussian grids should only provide a single output file using their native horizontal resolution and model levels. If reduced Gaussian grids are utilized a second file on the full Gaussian grid is requested.
If models are run with variable-resolution grids, we leave the choice of the best suitable interpolation grid to the modeling group. We ask to write all output variables for each experiment to the same file.
\subsection{Short Note on Data Analysis and Visualization}
We will provide NCAR Command Language (NCL) scripts to help visualize the model results and provide analysis functions. In addition, the DCMIP participants will have access to interactive Graphical User Interfaces (GUIs) to support the visualization and model intercomparison. Among the GUIs are the netCDF viewers Ncview and Panoply which are public domain tools and locally installed on the NCAR mirage server. In addition, we expect to provide some basic online visualization capabilities via NOAA's Live Access Server (LAS) software. A key to the successful visualization is the adherence of the output data sets to the netCDF-CF standard (see Appendix \ref{sec:netcdf}).
\subsection{Short Note on the Fortran Templates}
\label{sec:template}
We have provided a set of stand-alone Fortran routines that compute the initial conditions for all test cases. They are named
\begin{itemize}
\item \texttt{baroclinic\_wave\_test.f90}
\item \texttt{tropical\_cyclone\_test.f90}
\item \texttt{supercell\_test.f90}
\end{itemize}
We have also provided a set of stand-alone Fortran routines that compute tendencies associated with the physical parameterizations. They are named
\begin{itemize}
\item \texttt{terminator.f90}
\item \texttt{kessler.f90}
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\section{Moist Baroclinic Wave} \label{sec:baroclinic_wave}
The baroclinic instability test of \cite{ullrich2014proposed} considers a reference state in geostrophic and hydrostatic balance that satisfies the conditions for baroclinic instability. Although a perfect model should be able to maintain this state indefinitely, small truncation errors associated with numerical inaccuracies and grid structure will trigger the development of the wave modes associated with baroclinic development. To control the development of the baroclinic wave, a small perturbation (but one which is large compared with machine truncation) is added to the flow so as to trigger the development of a wave over a period of approximately 10 days. A moist variant of the dry dynamical test of \cite{ullrich2014proposed} is considered here so as to understand the impact of moisture feedbacks on the development of the wave.
This test case is similar in character to the test of \cite{jablonowski2006baroclinic}, but has a number of key differences: (1) this test is an analytical solution of the equations of motion in height ($z$) coordinates, (2) the bottom topography is zero throughout the domain, (3) the new test case does not have a distinct stratosphere (the presence of a stratosphere is largely irrelevant for understanding baroclinic development), and (4) the velocity field goes to zero at the model surface.
\begin{table}[h]
\caption{List of constants used for the Moist Baroclinic Wave test case}
\label{test4:tab}
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}lll}
\hline Constant & Value & Description \\
\hline
$z_{\tiny \mbox{top}}$ & $44000\ \mbox{m}$ & Recommended height position of the model top \\
$p_{\tiny \mbox{top}}$ & $\approx 2.26$ hPa & Recommended pressure at the model top\\
$X$ & $1$ & Reduced-size planet scaling factor, see below\\
$a$ & $a_{\tiny \mbox{ref}}/X$ & Scaled radius of the Earth \\
$\Omega$ & $\Omega_{\tiny \mbox{ref}}X$ & Scaled angular speed of the Earth \\
$p_s$ & $1000\ \mbox{hPa}$ & Surface pressure (constant) \\
$p_0$ & $1000\ \mbox{hPa}$ & Reference pressure (constant) \\
$u_0$ & $35\ \mbox{m\ s}^{-1}$ & Maximum amplitude of the zonal wind \\
$b$ & $2$ & Half-width parameter \\
$K$ & $3$ & Power used for temperature field \\
$T_E$ & $310\ \mbox{K}$ & Horizontal-mean temperature at the surface \\
$T_P$ & $240 \ \mbox{K}$ & Temperature at the polar surface\\
$u_p$ & $1\ \mbox{m\ s}^{-1}$ & Maximum amplitude of the zonal wind perturbation \\
$z_p$ & $15000\ \mbox{m}$ & Maximum height of the zonal wind perturbation \\
$\lambda_p$ & $\pi / 9$ & Longitude of the zonal wind perturbation centerpoint (20$^\circ$ E)\\
$\varphi_p$ & $2 \pi / 9$ & Latitude of the zonal wind perturbation centerpoint (40$^\circ$ N)\\
$R_p$ & $a / 10$ & Radius of the zonal wind perturbation \\
%$T_E$ & $310 \ \mbox{K}$ & Temperature at the equator\\
$\Gamma$ & $0.005\ \mbox{K\ m}^{-1}$ & Temperature lapse rate \\
$\Delta T$ & $4.8 \times 10^{5}\ \mbox{K}$ & Empirical temperature difference \\
$\varphi_w$ & $2 \pi / 9$ & Specific humidity latitudinal width parameter $(40^\circ)$\\
$p_w$ & $340\ \mbox{hPa}$ & Specific humidity vertical pressure width parameter \\
$q_0$ & $0.018$ kg/kg& Maximum specific humidity amplitude \\
$q_t$ & $1.0 \times 10^{-12}$ kg/kg & Specific humidity above artificial tropopause \\
$p_t$ & $10000\ \mbox{hPa}$ & Pressure at artificial tropopause \\
\hline
\end{tabular*}
\end{table}
\subsection{Reference State}
This section describes the analytical form of the reference state for the baroclinic wave. The test case is initialized with a constant surface pressure and with a surface geopotential equal to zero. The meridional wind in the reference state is zero.
In the reference state, the virtual temperature is given by
\begin{equation}
T_v(\varphi, z) = \frac{1}{\tau_1(z)-\tau_2(z) I_T(\varphi)},
\label{virtTemp}
\end{equation} where $I_T(\varphi)$ is defined as
\begin{equation}
I_{T}(\varphi) =(\cos \varphi )^K-\frac{K}{K+2}(\cos \varphi )^{K+2},
\end{equation} and $\tau_1(z)$ and $\tau_2(z)$ are defined as follows:
\begin{align}
\tau_1(z) &= \frac{1}{T_0} \exp\left(\frac{\Gamma z}{T_0}\right) + \left( \frac{T_0-T_P}{T_0T_P} \right)\left[1-2\left(\frac{z g}{b R_d T_0}\right)^2\right] \exp\left[-\left(\frac{z g}{b R_d T_0}\right)^2\right] \\
\tau_2(z) &= \frac{(K + 2)}{2} \left( \frac{T_E-T_P}{T_E T_P} \right) \left[1-2\left(\frac{z g}{b R_d T_0}\right)^2\right] \exp\left[-\left(\frac{z g}{b R_d T_0}\right)^2\right],
\end{align} with $T_0 = \tfrac{1}{2} (T_E + T_P)$. To maintain hydrostatic balance, the pressure is given by:
\begin{equation}
p(\varphi, z) = p_0\exp \left[ -\frac{g}{R_d}(\tau_{\text{int},1}(z) -\tau_{\text{int},2}(z) I_T(\varphi) ) \right]
\end{equation} with $\tau_{\text{int},1}(z)$ and $\tau_{\text{int},2}(z)$ given by
\begin{align}
\tau_{\text{int},1}(z) &=\frac{1}{\Gamma} \left[ \exp\left( \frac{\Gamma z}{T_0} \right)-1 \right] + z \left(\frac{T_0-T_P}{T_0T_P} \right) \exp\left[-\left(\frac{z g}{b R_d T_0}\right)^2\right] \\
\tau_{\text{int},2}(z) &=\frac{(K+2)}{2} \left(\frac{T_E-T_P}{T_E T_P} \right) z \exp\left[-\left(\frac{z g}{b R_d T_0}\right)^2\right].
\end{align} If density is a prognostic variable, it can be obtained from $p$, $T_v$ and the ideal gas law (\ref{eq:idealgaslaw}). Finally, the zonal velocity is
\begin{equation}
u_{\text{ref}}(\varphi, z) = -\Omega_{\text{ref}} a_{\text{ref}} \cos(\varphi)+\sqrt{(\Omega_{\text{ref}} a_{\text{ref}} \cos(\varphi))^2+ a_{\text{ref}} \cos(\varphi)U(z,\varphi))},
\end{equation} where
\begin{equation}
U(z, \varphi) = \frac{g K}{a_{\text{ref}}} \tau_{\text{int}_2}(z) \left[ (\cos \varphi)^{K - 1} - (\cos \varphi)^{K + 1} \right] T_v(\varphi, z).
\end{equation}
\subsection{Perturbations}
To trigger the development of the baroclinic wave, a perturbation is applied to the zonal velocity field that takes the form of a simple exponential bell with a vertical taper:
\begin{equation}
u^\prime(\lambda, \varphi, z) = \left\{ \begin{array}{ll} \displaystyle u_p Z_p(z) \exp \left[ - \left( \frac{R(\lambda, \varphi; \lambda_p, \varphi_p)}{R_p} \right)^2 \right], & \mbox{if $R(\lambda, \varphi; \lambda_p, \varphi_p) < R_p$,} \\ 0, & \mbox{otherwise,} \end{array} \right.
\end{equation} where
\begin{equation}
Z_p(z) = \left\{ \begin{array}{ll} \displaystyle 1 - 3 \left( \frac{z}{z_p} \right)^2 + 2 \left( \frac{z}{z_p} \right)^3, & \mbox{if $z \leq z_p$,} \\ 0, & \mbox{otherwise.} \end{array} \right.
\end{equation} Consequently, the perturbed velocity field takes the form
\begin{equation}
u(\lambda, \varphi, z) = u_{\text{ref}}(\varphi, z) + u^\prime(\lambda, \varphi, z).
\end{equation}
\subsection{Moist initial conditions}
We define the vertical $\eta$ coordinate as
\begin{equation}
\eta(\lambda, \varphi, z) = p(\lambda, \varphi, z) / p_s.
\end{equation} Since the surface pressure of the moist air $p_s$ is constant with $p_s = p_0 = 1000$ hPa the vertical coordinate $\eta$ is represented by $\eta = p/p_0$. Specific humidity is specified in terms of $\eta$ as
\begin{eqnarray}
q(\lambda, \varphi, \eta) &=& \left\{ \begin{array}{ll} \displaystyle q_0 \exp\Bigg[- \Big(\frac{\varphi}{\varphi_{w}}\Big)^4 \Bigg] \exp\Bigg[- \Bigg(\frac{(\eta-1)p_0}{p_{w}}\Bigg)^2 \Bigg], & \mbox{if $\eta > p_t / p_s$,} \\ q_{t}, & \mbox{otherwise.} \end{array} \right.
\end{eqnarray} The functional form of $q$ and its parameters were inspired by observations. This moisture fields leads to maximum relative humidities around 85\% in the lower levels of the midlatitudes.
Note that the moist temperature is colder than the temperature one would obtain with $q = 0$. However, note that in the moist case the virtual temperature and moist pressure determine the strength of the pressure gradient term in the momentum equations. Since these are identical to the temperature and pressure in the dry case, the forcing by the pressure gradient term is the same in both the dry and moist variant of the baroclinic wave. The moist variant of the baroclinic wave without the temperature forcing from large-scale condensation should lead to almost identical results when compared to the dry version. Very small variations are expected since the moisture gets independently transported as a passive tracer in this case and some models utilize the moist variant of the physical constant $c_p$. If possible, the dry $c_p$ should be used. Comparing the evolution of the dry baroclinic wave to its moist variant (without large-scale condensation) can serve as a first sensibility check.
\begin{figure}[tb]
\center\includegraphics[width=\linewidth]{plot_baroclinicwave_init.pdf}
\caption{Initial state for the moist baroclinic wave test.}\label{fig:baroclinicwave_init}
\end{figure}
%
%
%
%
%
\subsection{Terminator `toy'-chemistry}
The terminator `toy'-chemistry is presented in \cite{LCLVT2015GMD} and mimics photolysis-driven processes near the solar terminator. Two passive{\footnote{i.e. tracers do not feed back on the flow}} tracers, Cl and Cl$_2$, that are chemically reactive are transported. The sources and sinks are given by a simple, but non-linear, `toy' chemistry. As a result, strong gradients in the spatial distribution of the species develop near the edge of the terminator. Despite the large spatial variations in Cl and Cl$_2$ the weighted sum Cl$_y$=Cl+2Cl$_2$ should always be preserved in any flow field (if the initial condition for Cl$_y$ is constant). An overview of the `toy' terminator chemistry is given in Appendix \ref{sec:ToyChemistry}. The terminator test demonstrates how well the advection/transport scheme and/or physics-dynamics coupling preserves linear correlations.
The initial conditions for Cl and Cl$_2$ (in terms of {\bf{dry}} mixing ratios $q_{Cl}$ and $q_{Cl2}$, respectively) are the steady-state solutions to the terminator chemistry with no flow \cite{LCLVT2015GMD} (see Figure \ref{fig:terminator-2d-se-day0}):
\begin{figure}[t]
\includegraphics[width=\linewidth]{terminator-2d-se-day0.pdf}
\caption{Contour lines of $<q_{Cl}>/(4.0\times 10^{-6})$, (upper left), $<q_{Cl2}>/(4.0\times 10^{-6})$ (upper right), $<q_{Cly}>$ (lower left) and surface pressure (lower right) at day 0 (initial conditions).}
\label{fig:terminator-2d-se-day0}
\end{figure}
\begin{eqnarray}
q_{Cl}(\lambda,\theta,z,t=0)&=&D-r, \\
q_{Cl2}(\lambda,\theta,z,t=0)&=&\frac{1}{2}\left( q_{Cly} - D + r\right),
\end{eqnarray}
where $q_{Cly}=4.0\times 10^{-6}$kg/kg,
\begin{align}
r &= \frac{k_1}{4 k_2}, \\
D &= \sqrt{r^2 + 2 r q_{Cly}},
\end{align}
and reaction coefficients $k_1$ and $k_2$ are given in appendix \ref{sec:ToyChemistry} (equations \eqref{app:term_k1} and \eqref{app:term_k2}). Please note that the mixing ratios are dry, i.e. the ratio between the density of the species and the density of dry air.
The forcing terms are computed analytically (assuming no flow) over a (physics/chemistry) time-step $\Delta t$:
\begin{equation}
F_{Cl}^n= - L_{\Delta t} \frac{ (q_{Cl}^n - D + r ) ( {q_{Cl}}^n + D + r ) }{ 1+e^{-4 k_2 D \Delta t} + \Delta t L_{\Delta t} (q_{Cl}^n + r ) }.
\end{equation}
where $q_{Cl}^n$ is the value of $q_{Cl}$ at the beginning of the n'th time step and
\begin{equation}
L_{\Delta t} =
\left\{
\begin{array}{cc}
\frac{1-e^{-4 k_2 D \Delta t}}{D \Delta t} &\mbox{if } D > 0 \\[1ex]
4k_2 &\mbox{if } D = 0.
\end{array}
\right. \label{eq:linvsexp}
\end{equation}
and by conservation,
\begin{equation}
F_{\mathrm{Cl}_2 }^n = -\frac{1}{2} F_{\mathrm{Cl} }^n .
\end{equation}
In implementation, $L_{\Delta t} $ needs some care. As $4 k_2 D \Delta t$ approaches machine precision, it is useful to simply use the formula for $D=0$ rather than the expression for $D>0$. The chemistry/physics updated mixing ratios are given by
\begin{eqnarray}
q^n_{Cl}+\Delta t\, F_{Cl}^n,\\
q^n_{Cl2}+\Delta t\, F_{Cl2}^n.
\end{eqnarray}
In terms of Fortran code the analytical forcing is given by:
\begin{verbatim}
! dt is size of physics time step
cly = cl + 2.0*cl2
r = k1 / (4.0*k2)
d = sqrt( r*r + 2.0*r*cly )
e = exp( -4.0*k2*d*dt )
if( abs(d*k2*dt) .gt. 1e-16 )
el = (1.0-e) / (d*dt)
else
el = 4.0*k2
endif
f_cl = -el * (cl-d+r) * (cl+d+r) / (1.0 + e + dt*el*(cl+r))
f_cl2 = -f_cl / 2.0
\end{verbatim}
The reaction rates are defined by
\begin{verbatim}
! k1 and k2 are reaction rates
k1_lat_center = 20.0 ! degrees
k1_lon_center = 300.0 ! degrees
k1 = max(0.d0,sin(lat)*sin(k1_lat_center)
+ cos(lat)*cos(k1_lat_center)*cos(lon-k1_lon_center))
k2 = 1.0
\end{verbatim}
The initial condition is defined by
\begin{verbatim}
cly = 4.0e-6
r = k1 / (4.0*k2)
d = sqrt( r*r + 2.0*cly*r )
cl = d-r
cl2 = cly / 2.0 - (d-r) / 2.0
\end{verbatim}
Fortran subroutines that, given $(\lambda,\theta)$ will return the tendencies for $q_{Cl}$ (mixing ratio) and $q_{Cl2}$ (mixing ratio), are provided as supplemental material in \cite{LCLVT2015GMD}. Similarly for the forcing terms.
The `toy' terminator chemistry test here uses the baroclinic wave initialization (described in section \ref{sec:baroclinic_wave}) based on the moist setup and the flow will transport the two species that interact non-linearly with each other through the toy chemistry. A physics time-step of 30 minutes and 15 minutes, respectively, is used. Note that if the user has the baroclinic wave setup the only additional work is to initialize two tracers and implement the chemistry. Example solutions are shown for HOMME (High-Order Method Modeling Environment) \cite{DetAl2012IJHPCA} spectral elements and HOMME-CSLAM \cite{LetAl2016MWR}. The latter model is based on CSLAM \cite{LNU2010JCP} transport of tracers consistently coupled with spectral element dynamics.
\subsubsection{Diagnostics}
If the initial conditions for $q_{Cl}$ and $2q_{Cl2}$ add up to a constant (as is the case in this setup) then no matter how the individual species evolve the weighted sum $q_{Cly}=q_{Cl}+2q_{Cl2}$ should be constant in space and time. Hence the analytical solution for $q_{Cly}$ is known. The terminator chemistry preserves the linear relationship between $q_{Cl}$ and $q_{Cl2}$ so the only causes for this relationship to break are:
\begin{itemize}
\item the transport operator (usually the limiter/filter) does not exactly preserve linear relations, and/or,
\item physics-dynamics coupling breaks the relationship (see, e.g., Figure \ref{fig:terminator-diags}).
\end{itemize}
The following diagnostics are used in this test case:
\begin{itemize}
\item Average column integrated mixing ratio (two-dimensional variable):
\begin{equation}
<q>(\lambda,\theta)=\frac{\int_{z=0}^{z_{top}} q(\lambda,\theta,z)\, \rho(\lambda,\theta,z)\, dz}{\int_{z=0}^{z_{top}} \rho(\lambda,\theta,z) dz}.
\end{equation}
where $q=q_{Cl},q_{Cl2}$ and $\rho$ is density. The global integrals should be computed consistently with the numerical method (preferably `inline' in the source code on the native grid and not using interpolated data).
\item $\ell_2(t)$, $\ell_\infty(t)$ and relative mass change $\Delta M(t)$ error norms for Cl$_y$:
\begin{align}
\ell_2(t)&=\frac{\sqrt{I\left[ \left(<q_{Cly}>-4.0\times 10^{-6}\right)^2\right] }}{\sqrt{I\left[\left(4.0\times 10^{-6}\right)^2\right]}},\\
\ell_\infty(t)&=\frac{\max_{\text{all }\lambda, \theta} |<q_{Cly}>-4.0\times 10^{-6}|}{4.0\times 10^{-6}},\\
\Delta M(t)&=\frac{\int_{z=0}^{z_{top}} \left\{ I\left[q_{Cly}\rho\right]\right\}\, dz -M_0}{M_0}
\end{align}
respectively, where
\begin{equation}
<q_{Cly}>=<q_{Cl}+2q_{Cl2}>.
\end{equation}
and $M_0$ is the initial mass of Cl$_y$
\begin{equation}
M_0=\int_{z=0}^{z_{top}} \left\{ I\left[4.0\times 10^{-6}\rho \right]\, \right\} dz.
\end{equation}
\end{itemize}
$I\left[ \cdot \right]$ is the global horizontal integral consistent with the numerical method.\\
Note that if the physics-dynamics coupling procedure breaks $Cl_y$ conservation it can be very subtle but shows in the $\Delta M(t)$ diagnostic (see, e.g., Figure \ref{fig:terminator-diags}).
%\begin{equation}
%q_{Cly}=q_{Cl}+2q_{Cl2}.
%\end{equation}
%The exact solution is that $q_{Cly}$ remains constant in time:
%\begin{equation}
%q_{Cly}(t)=q_{Cly}(t=0)=4.0\times 10^{-6}\, kg/kg.
%\end{equation}
\begin{figure}[tb]
\center\includegraphics[width=\linewidth]{terminator_diags-t_phys900.pdf}
\center\includegraphics[width=\linewidth]{terminator_diags-t_phys1800.pdf}
\caption{Global error norms $\ell_2(t)$ (column 1), $\ell_\infty(t)$ (column 2) and relative mass change $\Delta M(t)$ (column 3) for $q_{Cly}$ for HOMME (row 1 and 3) and HOMME-CSLAM (row 2 and 4) with physics/chemistry time-step of 15 minutes (row 1 and 2) and 30 minutes (row 3 and 4), respectively. The non-conservation of Cl$_y$ mass ($\Delta M(t)$) for physics time-step of 30 minutes is due to physics-dynamics coupling in which the tendencies are altered if they result in negative mixing ratios. For HOMME-CSLAM this happens at one point around day 16 and 18, respectively, whereas it happens frequently (in time and space) for HOMME.}\label{fig:terminator-diags}
\end{figure}
%\cite{LCLVT2015GMD}
%\subsection{Grid spacings, simulation time, output and diagnostics}
%
%\begin{itemize}
%\item Dry simulations with the Toy Chemistry module (Appendix \ref{sec:ToyChemistry}) should be performed at 1$^\circ$ resolution with 30 vertical levels for 12 days.
%\item Plots of $q_{Cl}$, $q_{Cl2}$ and $q_T$ at 5000m altitude should be produced after $1$, $5$ and $12$ days.
%\item Experiments could address the coupling frequency between the dynamics and physics.
%\end{itemize}
\subsection{Grid spacings, simulation time, output and diagnostics}
\begin{itemize}
\item Reference simulations (dry and moist) should be performed at 1$^\circ$ resolution with 30 vertical levels for 30 days. Plots should be produced for the moist simulation and the anomaly between moist and dry simulations at day 9, 12 and 15.
\item Plots of minimum surface pressure over the duration of the simulation for both dry and moist configurations.
\item Experiments could address the coupling frequency between the dynamics and physics.
\item A variable resolution simulation should be performed that (a) studies the effect of the baroclinic wave transitioning from coarse resolution to fine resolution and (b) studies the effect of enhanced resolution near the front.
\item Terminator chemistry (use physics time-step of 30 minutes and 15 minutes):
\begin{itemize}
\item Please plot contour lines for $<q_{Cl}>/(4.0\times 10^{-6})$ at day 9. Contour interval must be 0.1 with zero contour. Please offset the zero contour \verb+-1.0E-12+ to avoid contouring round-off undershoots. Contour levels used in Figure \ref{fig:terminator-2d-se} are\\
(\verb+-0.1,-1.0E-12,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2+)
\item Please plot contour lines for $<q_{Cl2}>/(4.0\times 10^{-6})$ at day 9. Contour interval must be 0.04 with round-off offset zero contour. Contour levels used in Figure \ref{fig:terminator-2d-se} are\\
(\verb+-0.04,-1.0E-12,0.04,0.08,0.12,0.16,0.20,0.24,0.28,+\\
\verb+0.32,0.36,0.40,0.44,0.48+)
\item Please plot contour lines for $<q_{Cly}>/(4.0\times 10^{-6})$ at day 9. Contour interval must be 0.04 (if errors are less than 4\% and non-zero please use a contour interval that shows the errors) excluding the 1.0 contour but symmetric about 1.0. Contour levels used in Figure \ref{fig:terminator-2d-se} are\\
(\verb+0.78,0.82,0.86,0.90,0.94,0.98,1.02,1.06,1.1,1.14,1.18,1.22,1.24,+
\verb+1.08,1.12,1.16,1.20+)
\item Please plot global error norms $\ell_2(t)$, $\ell_\infty(t)$ and relative mass change $\Delta M(t)$ for $q_{Cly}$ as a function of time from day 0 to 30 preferably with 3 hourly time spacing. Use vertical axis adjusted to your data. See example on Figure \ref{fig:terminator-diags}.
\end{itemize}
\end{itemize}
\begin{figure}[t]
\includegraphics[width=0.9\linewidth]{terminator-2d-se-day9.pdf}
\includegraphics[width=0.9\linewidth]{terminator-2d-cslam-day9.pdf}
\caption{Upper panel of 4 plots is contour lines of $<q_{Cl}>/(4.0\times 10^{-6})$, (upper left), $<q_{Cl2}>/(4.0\times 10^{-6})$ (upper right), $<q_{Cly}>$ (lower left) and surface pressure (lower right) at day 9 simulated with HOMME at approximately $1^\circ$ resolution (30x30 elements on each cubed-sphere panel with 4x4 quadrature points in each element). Lower panel of 4 plots is the same as the upper panel but simulated with HOMME-CSLAM with 3x3 control volumes within each element. The results are based on the moist baroclinic wave setup (with no moist processes) and a physics-coupling time-step of 15 minutes.}
\label{fig:terminator-2d-se}
\end{figure}
\clearpage
\section{Tropical cyclone} \label{sec:tropical_cyclone}
The simplified tropical cyclone test case on a regular-size Earth is based on the work of \cite{reed2012idealized, reed2011analytic,reed2011impact, reed2011assessing}. In this test an analytic vortex is initialized in a background environment which is tractable to a rapid intensification of tropical cyclones.
\begin{table}[h]
\caption{List of constants used for the Ideaized Tropical Cyclone test}
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}lll}
\hline Constant & Value & Description \\
\hline
$X$ & $1$ & small-planet scaling factor (regular-size Earth)\\
$z_t$ & $15000$ m & Tropopause height \\
$q_0$ & $0.021$ kg/kg & Maximum specific humidity amplitude \\
$q_t$ & $10^{-11}$ kg/kg & Specific humidity in the upper atmosphere \\
$T_0$ & $302.15$ K & Surface temperature of the air \\
$T_s$ & $302.15$ K & Sea surface temperature (SST), 29 C$^\circ$\\
$z_{q1}$ & $3000$ m & Height related to the linear decrease of $q$ with height \\
$z_{q2}$ & $8000$ m & Height related to the quadratic decrease of $q$ with height \\
$\Gamma$ & $0.007$\ K\ m$^{-1}$ & Virtual temperature lapse rate \\
$p_{b}$ & $1015$ hPa & Background surface pressure \\
$\varphi_c$ & $\pi / 18$ & Initial latitude of vortex center (radians) \\
$\lambda_c$ & $\pi$ & Initial longitude of vortex center (radians) \\
$\Delta p$ & $11.15$ hPa & Pressure perturbation at vortex center \\
$r_p$ & $282000$ m & Horizontal half-width of pressure perturbation \\
$z_p$ & $7000$ m & Height related to the vertical decay rate of $p$ perturbation \\
$\epsilon$ & $10^{-25}$ & Small threshold value \\
\hline
\end{tabular*}
\end{table}
\subsection{Initialization}
The background state consists of a prescribed specific humidity profile, virtual temperature and pressure profile. The initial profile is defined to be in approximate gradient wind balance. The vertical sounding is chosen to roughly match an observed tropical sounding documented in \cite{jordan1958mean}. The background specific humidity profile $\overline{q}(z)$ as a function of height $z$ is
\begin{equation}
\begin{split}
\overline{q}(z)&=q_0 \exp\left(- \frac{z}{z_{q1}}\right)\exp\left[-\left(\frac{z}{z_{q2}}\right)^2\right] \text{ ~~for } 0 \leq z \leq z_t \\
\overline{q}(z)&=q_t \text{ ~~for } z_t \leq z
\end{split}
\end{equation}
The background virtual temperature sounding $\overline{T}_v(z)$ is split into two different representations for the lower and upper atmosphere. It is given by
\begin{equation}
\begin{array}{ll} \label{eq:tc_virtualtemperaturebg}
%\phantom{T_{vt} = }\overline{T}_v(z) = T_{v0} - \Gamma z & \mbox{for} \; 0 \le z \le z_t, \\
\overline{T}_v(z) = T_{v0} - \Gamma z & \mbox{for} \; 0 \le z \le z_t, \\
\overline{T}_v(z) = T_{vt} = T_{v0} - \Gamma z_t & \mbox{for} \; z_t < z,
\end{array}
\end{equation} with the virtual temperature at the surface $T_{v0}$ = $T_0 (1+0.608 \, q_0)$ and the virtual temperature at the tropopause level $T_{vt}$ = $T_{v0} - \Gamma z_t$. The background temperature profile can be obtained from (\ref{eq:virtualtemperature}).
The background vertical pressure profile $\overline{p}(z)$ of the moist air is computed using the hydrostatic balance and (\ref{eq:tc_virtualtemperaturebg}). The profile is given by:
\begin{equation}
\begin{array}{ll}\label{eq4}
\displaystyle \overline{p}(z) = p_b \left( \frac{T_{v0} - \Gamma z}{T_{v0}} \right )^{g / R_d \Gamma} & \mbox{for} \; 0 \le z \le z_t, \\
\displaystyle \overline{p}(z) = p_t \, \exp{\left(\frac{g (z_t - z)}{R_d T_{vt}} \right)} & \mbox{for} \; z_t < z.
\end{array}
\end{equation} The pressure at the tropopause level $z_t$ is continuous and given by
\begin{equation}\label{eq4.5}
p_t = p_b \left( \frac{T_{vt}}{T_{v0}} \right )^{\frac{g}{R_d \Gamma}},
\end{equation}
which, for the given set of parameters, is approximately 130.5 hPa.
\subsubsection{Axisymmetric Vortex}
The pressure equation $p(r,z)$ for the moist air is comprised of the background pressure profile (\ref{eq4}) plus a 2D pressure perturbation $p'(r,z)$,
\begin{equation} \label{eq5}
p(r,z) = \overline{p}(z) + p^\prime(r,z),
\end{equation} where $r$ symbolizes the radial distance (or radius) to the center of the prescribed vortex. On the sphere $r$ is defined using the great circle distance
\begin{equation}
r = a \arccos{ \left ( \sin{\varphi_c} \, \sin{\varphi} + \cos{\varphi_c} \, \cos{\varphi} \, \cos{(\lambda - \lambda_c)} \right )}.
\end{equation} The perturbation pressure is defined as
\begin{align} \label{test5:p_pert}
p^\prime(r,z) & = -\Delta p \, \exp\left[{-\left (\frac{r}{r_p} \right ) ^{3/2}} {-\left (\frac{z}{z_p} \right ) ^{2}}\right] \left ( \frac{T_{v0} - \Gamma z}{T_{v0}} \right )^{\frac{g}{R_d \Gamma}} & & \mbox{for $\; 0 \le z \le z_t$}, \nonumber \\
p^\prime(r,z) & = 0 & & \mbox{for} \; z_t < z.
\end{align} The pressure perturbation depends on the pressure difference $\Delta p$ between the background surface pressure $p_b$ and the pressure at the center of the initial vortex, the pressure change in the radial direction $r_p$ and the pressure decay with height within the vortex $z_p$. The moist surface pressure $p_s(r)$ is computed by setting $z = 0$ m in (\ref{eq5}), which gives
\begin{equation}
\label{eq:ps}
p_s(r) = p_b - \Delta p \, \exp\left[{-\left (\frac{r}{r_p} \right ) ^{3/2}}\right].
\end{equation}
The axisymmetric virtual temperature $T_v(r,z)$ is computed using the hydrostatic equation and ideal gas law
\begin{equation}
T_v(r,z) = -\frac{g p(r,z)}{R_d} \left( \frac{\partial p(r,z)}{ \partial z} \right)^{-1}.
\end{equation} Again it can be written as a sum of the background state and a perturbation,
\begin{equation} \label{eq:virt_temp}
T_v(r,z) = \overline{T}_v(z) + T_v^\prime(r,z),
\end{equation} where the virtual temperature perturbation is defined as
\begin{align}
\label{eq:Tv}
T_v^\prime(r,z) &= (T_{v0} - \Gamma z ) \left\{ \left [1+ \frac{2R_d(T_{v0} - \Gamma z)z}{gz_p^2 \left[ 1 - \frac{p_b}{\Delta p}\exp\left({\left (\frac{r}{r_p} \right ) ^{3/2}} + {\left (\frac{z}{z_p} \right ) ^{2}} \right) \right] }\right]^{-1} - 1 \right\} & & \mbox{for} \; 0 \le z \le z_t, \nonumber \\
T_v^\prime(r,z) &= 0 & & \mbox{for} \; z_t < z.
\end{align}
The axisymmetric specific humidity $q(r,z)$ is set to the background profile everywhere
\begin{eqnarray}
\label{eq:q}
q(r,z) = \overline{q}(z).
\end{eqnarray} Consequently, the temperature can be written as
\begin{equation} \label{test5:T_eqn}
T(r,z) = \overline{T}(z) + T^\prime(r,z),
\end{equation} with the temperature perturbation
\begin{align} \label{eq:temperature}
T^\prime(r,z) &= \frac{T_{v0} - \Gamma z}{1+0.608\overline{q}(z)} \left\{ \left [1+ \frac{2R_d(T_{v0} - \Gamma z)z}{gz_p^2 \left[ 1 - \frac{p_b}{\Delta p}\exp\left({\left (\frac{r}{r_p} \right ) ^{3/2}} + {\left (\frac{z}{z_p} \right ) ^{2}} \right) \right] }\right]^{-1} - 1 \right\} & & \mbox{for} \; 0 \le z \le z_t, \nonumber \\
T^\prime(r,z) &= 0 & & \mbox{for} \; z_t < z.
\end{align}
Due to the small specific humidity value in the upper atmosphere (10$^{-11}$ kg/kg for $z > z_t$) the virtual temperature equals the temperature to a very good approximation in this region. The formulation presented here is equivalent to the one presented in \cite{reed2012idealized}.
If the density of the moist air needs to be initialized its formulation is based on the ideal gas law
\begin{equation} \label{eq:density}
\rho(r,z) = \frac{p(r,z)}{R_d T_v(r,z)}
\end{equation}
which utilizes the moist pressure (\ref{eq5}) and virtual temperature (\ref{eq:virt_temp}). The surface elevation $z_s$ and thereby the surface geopotential $\Phi_s=g z_s$ are set to zero.
Finally, the tangential velocity field $v_T(r,z)$ of the axisymmetric vortex is defined by utilizing the gradient-wind balance, which depends on the pressure (\ref{eq5}) and the virtual temperature (\ref{eq:Tv}). The tangential velocity is given by
\begin{eqnarray}
\label{eq:gradient-wind}
v_T(r,z) = -\frac{f_cr}{2}+\sqrt{ \frac{f_c^2r^2}{4}+\frac{R_d \, T_v(r,z) \, r}{p(r,z)} \frac{\partial p(r,z)}{\partial r}},
\end{eqnarray}
where $f_c = 2 \Omega \sin(\varphi_c)$ is the Coriolis parameter at the constant latitude $\varphi_c$. Substituting $T_v(r,z)$ and $p(r,z)$ into ({\ref{eq:gradient-wind}) gives
\begin{align}
\label{eq:gradient-wind-expr}
v_T(r,z) & = -\frac{f_cr}{2}+\sqrt{ \frac{f_c^2r^2}{4}-\frac{\frac{3}{2} \left( \frac{r}{r_p}\right)^{3/2} (T_{v0}-\Gamma z) R_d}{1+\frac{2R_d(T_{v0}-\Gamma z)z}{g z_p^2}-\frac{p_b}{\Delta p}\exp\left({\left (\frac{r}{r_p} \right ) ^{3/2}} + {\left (\frac{z}{z_p} \right ) ^{2}} \right)}} & & \mbox{for} \; 0 \le z \le z_t, \nonumber \\
v_T(r,z) & = 0 & & \mbox{for} \; z_t < z.
\end{align} The last step is to split the tangential velocity (\ref{eq:gradient-wind-expr}) into its zonal and meridional wind components $u(\lambda,\varphi,z)$ and $v(\lambda,\varphi,z)$. Similar to \cite{nair2008moving} these are computed using the following expressions,
\begin{eqnarray}
d_1 &=& \sin\varphi_c \, \cos\varphi - \cos\varphi_c \, \sin\varphi \, \cos(\lambda-\lambda_c) \\
d_2 &=& \cos\varphi_c \, \sin(\lambda-\lambda_c) \\
d &=& \max \big({\epsilon,\sqrt{ {d_1}^2 + {d_2}^2} } \big),
\end{eqnarray}
which are utilized in the projections
\begin{eqnarray}
\label{eqn:u_wind}
u(\lambda,\varphi,z) &=& \frac{v_T(\lambda,\varphi,z) \, d_1}{d}\\ \label{eqn:v_wind}
v(\lambda,\varphi,z) &=& \frac{v_T(\lambda,\varphi,z) \, d_2}{d} \,.
\end{eqnarray}
A small $\epsilon = 10^{-25}$ value avoids divisions by zero. The vertical velocity is set to zero.
\ \\
\noindent \begin{tabular}{|p{\textwidth}|}
\hline \textbf{Note:} We are currently investigating a test case specification that places the idealized tropical cyclone on an $f$-plane. This configuration would remove issues associated with $\beta$ drift of the cyclone and allow for a more direct intercomparison of the simulated storm. \\
\hline
\end{tabular}
\begin{figure}[tb]
\center\includegraphics[width=\linewidth]{plot_tropicalcyclone_init.pdf}
\caption{Initial state for the tropical cyclone test.}\label{fig:tropicalcyclone_init}
\end{figure}
\subsection{Grid spacings, simulation time, output and diagnostics}
\begin{itemize}
\item Moist simulations should be performed at 0.5$^\circ$ resolution with 30 vertical levels for 10 days.
\item Plots of minimum surface pressure over the duration of the simulation.
\item Experiments could address the coupling frequency between the dynamics and physics.
\item A variable resolution simulation should be performed that (a) studies the effect of the tropical cyclone transitioning from fine resolution to coarse resolution and (b) high resolution simulations down to 0.125$^\circ$ over the tropical cyclone.
\end{itemize}
\clearpage
\section{Supercell} \label{sec:supercell}
The supercell test permits the study of a non-hydrostatic moist feature with strong vertical velocities and associated precipitation and is based on the work of \cite{klemp1978simulation}.
\begin{table}[h]
\caption{List of constants used for the Supercell test}
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}lll}
\hline Constant & Value & Description \\
\hline
$X$ & $120$ & Small-planet scaling factor (reduced Earth)\\
$\theta_{tr}$ & $343$ K & Temperature at the tropopause \\
$\theta_0$ & $300$ K & Temperature at the equatorial surface \\
$z_{tr}$ & $12000$ m & Altitude of the tropopause \\
$T_{tr}$ & $213$ K & Temperature at the tropopause \\
$U_s$ & $30$ m/s & Maximum zonal wind velocity \\
$U_c$ & $15$ m/s & Coordinate reference velocity \\
$z_{s}$ & $5000$ m & Lower altitude of maximum velocity \\
$\Delta z_{u}$ & $1000$ m & Transition distance of velocity \\
$\Delta \theta$ & $3$ K & Thermal perturbation magnitude \\
$\lambda_p$ & $0$ & Thermal perturbation longitude \\
$\varphi_p$ & $0$ & Thermal perturbation latitude \\
$r_p$ & $X \times 10000$ m & Perturbation horizontal half-width \\
$z_{c}$ & $1500$ m & Perturbation center altitude \\
$z_{p}$ & $1500$ m & Perturbation vertical half-width \\
\hline
\end{tabular*}
\end{table}
It is assumed that the saturation mixing ratio is given by
\begin{equation}
q_{vs}(p,T) = \left( \frac{380.0}{p} \right) \exp\left(17.27 \times \frac{T-273.0}{T-36.0}\right)
\end {equation} The definition of this test case relies on hydrostatic and gradient wind balance, written in terms of Exner pressure $\pi$ and virtual potential temperature $\theta_v$ as
\begin{equation} \label{eq:BalanceEq}
\pdiff{\pi}{z} = - \frac{g}{c_p \theta_v}, \quad \mbox{and} \quad u^2 \tan \varphi = - c_p \theta_v \pdiff{\pi}{\varphi}.
\end{equation} These equations can be combined to eliminate $\pi$, leading to
\begin{equation} \label{eq:CombinedBalanceEq}
\pdiff{\theta_v}{\varphi} = \frac{\sin (2 \varphi)}{2 g} \left( u^2 \pdiff{\theta_v}{z} - \theta_v \pdiff{u^2}{z} \right).
\end{equation}
The wind velocity is analytically defined throughout the domain. Meridional and vertical wind is initially set to zero. The zonal wind is obtained from
\begin{equation}
\overline{u}(\varphi,z) = \left\{ \begin{array}{ll}
\displaystyle \left(U_s\frac{z}{z_s}-U_c\right)\cos(\varphi), & \mbox{for $z < z_s - \Delta z_u$}, \\[2.0ex]
\displaystyle \left[ \left(-\frac{4}{5}+3\frac{z}{z_s}-\frac{5}{4}\frac{z^2}{z_s^2}\right)U_s-U_c \right] \cos(\varphi), & \mbox{for $\vert z-z_s \vert \leq \Delta z_u$} \\[2.0ex]
\displaystyle \left(U_s-U_c\right)\cos(\varphi), & \mbox{for $z > z_s + \Delta z_u$}
\end{array} \right.
\end{equation}
The equatorial profile is determined through numerical iteration. Potential temperature at the equator is specified via
\begin{equation}
\theta_{\text{eq}}(z) = \left\{ \begin{array}{ll} \displaystyle \theta_0 + (\theta_{tr} - \theta_0)\left(\frac{z}{z_{tr}}\right)^{\frac{5}{4}} & \mbox{for $0 \leq z \leq z_{tr}$,} \\[2.0ex]
\displaystyle \theta_{tr} \exp\left(\frac{g(z-z_{tr})}{c_pT_{tr}}\right) & \mbox{for $z_{tr} \leq z$} \end{array} \right.
\end{equation} And relative humidity is given by
\begin{equation}
\overline{H}(z) = \left\{ \begin{array}{ll} \displaystyle 1 + \frac{3}{4}\left(\frac{z}{z_{tr}}\right)^{5/4} & \mbox{for $0 \leq z \leq z_{tr}$,} \\[2.0ex]
\displaystyle \frac{1}{4} & \mbox{for $z_{tr} \leq z$.} \end{array} \right.
\end{equation} Pressure and temperature at the equator are obtained by iterating on hydrostatic balance with initial state
\begin{equation}
\theta_{v,\text{eq}}^{(0)}(z) = \theta_{\text{eq}}(z),
\end{equation} and iteration procedure
\begin{align}
\pi_{\text{eq}}^{(i)} &= 1 - \int_{0}^{z} \frac{g}{c_p \theta_{v,\text{eq}}^{(i)}} dz \\
p_{\text{eq}}^{(i)} &= p_0 (\pi^{(i)})^{c_p / R_d} \\
T_{\text{eq}}^{(i)} &= \theta_{\text{eq}}(z) \pi_{\text{eq}}^{(i)} \\
q^{(i)}_{\text{eq}} &= H(z) q_{vs}(p_{\text{eq}}^{(i)}, T_{\text{eq}}^{(i)}) \\
\theta_{v,\text{eq}}^{(i+1)} &= \theta_{\text{eq}}(z) (1 + M_v q^{(i)}_{\text{eq}})
\end{align} This iteration procedure appears to converge to machine epsilon after approximately 10 iterations. The equatorial moisture profile is then extended through the entire domain,
\begin{equation}
q(z, \varphi) = q_{\text{eq}}(z).
\end{equation}
Once the equatorial profile has been constructed, the virtual potential temperature through the remainder of the domain can be computed by iterating on (\ref{eq:CombinedBalanceEq}),
\begin{align}
\theta_v^{(i+1)}(z,\varphi) = \theta_{v,\text{eq}}(z) + \int_{0}^{\varphi} \frac{\sin(2 \phi)}{2 g} \left(\overline{u}^2 \pdiff{\theta_{v}^{(i)}}{z} - \theta_v^{(i)} \pdiff{\overline{u}^2}{z} \right) d\varphi.
\end{align} Again, approximately 10 iterations are needed for convergence to machine epsilon. Once virtual potential temperature has been computed throughout the domain, Exner pressure throughout the domain can be obtained from (\ref{eq:BalanceEq}),
\begin{equation}
\pi(z,\varphi) = \pi_{eq}(z) - \int_{0}^{\varphi} \frac{u^2 \tan \varphi}{c_p \theta_v} d\varphi,
\end{equation} and so
\begin{align}
p(z,\varphi) &= p_0 \pi(z,\varphi)^{c_p / R_d}, \\
T_v(z,\varphi) &= \theta_v(z,\varphi) (p / p_0)^{R_d / c_p}.
\end{align}
\subsection{Potential temperature perturbation}
To initiate convection, a thermal perturbation is introduced in the initial potential temperature field:
\begin{equation}
\theta^\prime(\lambda,\phi,z) = \left\{ \begin{array}{ll} \displaystyle \Delta \theta \cos^2\left(\frac{\pi}{2}R_{\theta}(\lambda, \varphi, z)\right) & \mbox{for $R_{\theta}(\lambda, \varphi, z) < 1$,} \\[2.0ex]
0 & \mbox{for $R_{\theta}(\lambda, \varphi, z) \geq 1$,} \end{array} \right.
\end{equation} where
\begin{equation}
R_{\theta}(\lambda, \varphi, z) = \left[ \left( \frac{R_c(\lambda, \varphi; \lambda_p, \varphi_p)}{r_p} \right)^2 + \left( \frac{z - z_c}{z_p} \right)^2 \right]^{1/2}.
\end{equation}
\ \\
\noindent \begin{tabular}{|p{\textwidth}|}
\hline \textbf{Note:} An additional iterative step will be required here to bring the potential temperature perturbation into hydrostatic balance. Without this additional iteration, large vertical velocities will be generated as the model rapidly adjusts to hydrostatic balance. \\
\hline
\end{tabular}
\begin{figure}[tb]
\center\includegraphics[width=\linewidth]{plot_supercell_init.pdf}
\caption{Initial state for the supercell test.}\label{fig:supercell_init_p1}
\end{figure}
\begin{figure}[tb]
\center\includegraphics[width=\linewidth]{plot_supercell_init_p2.pdf}
\caption{Initial state for the supercell test (cont'd).}\label{fig:supercell_init_p2}
\end{figure}
\subsection{Uniform Diffusion} \label{sec:supercell_uniform_diffusion}
The supercell test case requires uniform diffusion to be continuously applied to all prognostic variables. This diffusion takes the form
\begin{align}
\pdiff{\psi}{t} &= \ldots + \nu_s \nabla^2 \psi,
\end{align} or
\begin{align}
\pdiff{}{t} (\rho \psi) &= \ldots + \nu_s \nabla \cdot ( \rho \nabla \psi ),
\end{align} for scalar quantities $\psi$ and $q$, and
\begin{align}
\pdiff{\vb{u}}{t} &= \ldots + \nu_v \nabla^2 \vb{u},
\end{align} for the vector velocity $\vb{u}$. The diffusion coefficients used in this test are $\nu_s = 1500\ \mbox{m}^2\ \mbox{s}^{-1}$ and $\nu_v = 500\ \mbox{m}^2\ \mbox{s}^{-1}$.
\subsection{Physical Parameterizations}
For initial verification, all physical parameterizations should be disabled (except for uniform diffusion, as described in section \ref{sec:supercell_uniform_diffusion}).
\subsection{Grid spacings, simulation time, output and diagnostics}
\begin{itemize}
\item Moist simulations should be performed at 0.5$^\circ$, 1$^\circ$, 2$^\circ$ and 4$^\circ$ resolution with 40 uniformly spaced vertical levels for 120 minutes (7200 seconds).
\item Maximum vertical velocity and instantaneous precipitation rate should be reported at 60 second increments.
\item Plots of vertical velocity and rainwater should be produced at 5 km altitude after 30, 60, 90 and 120 minutes over the domain $[30W, 30E] \times [40S, 40N]$ (script provided).
\item Experiments could address the coupling frequency between the dynamics and physics.
\item A variable resolution simulation should be performed that (a) studies the effect of the supercell transitioning from fine resolution to coarse resolution and (b) high resolution simulations down to 0.125$^\circ$ over the supercell.
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliographystyle{wileyj}
\bibliography{DCMIP2016}
\begin{appendix}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Kessler Physics} \label{sec:KesslerPhysics}
The cloud microphysics update according to the following equation set:
\begin{alignat}{5}
\frac{\Delta \theta}{\Delta t} = & - \frac{L}{c_p \pi} & \Big( \frac{\Delta q_{vs}}{\Delta t} & + E_r \Big) & \\
\frac{\Delta q_v}{\Delta t} = & & \frac{\Delta q_{vs}}{\Delta t} & + E_r \\
\frac{\Delta q_c}{\Delta t} = & & - \frac{\Delta q_{vs}}{\Delta t} & & - A_r & - C_r \\
\frac{\Delta q_r}{\Delta t} = & & & - E_r & + A_r & + C_r & - V_r \pdiff{q_r}{z},
\end{alignat} where $L$ is the latent heat of condensation, $A_r$ is the autoconversion rate of cloud water to rain water, $C_r$ is the collection rate of rain water, $E_r$ is the rain water evaporation rate, and $V_r$ is the rain water terminal velocity.
The pressure follows from the equation of state
\begin{equation}
p=\rho R_dT(1+0.61q_v)
\end{equation} with $p$ the pressure, $\rho$ the density of moist air, $R_d$ the gas constant for dry air, $T$ the temperature and $q_v$ the mixing ratio of water vapor. The equation is rewritten as a nondimensional pressure $\Pi$ equation.
\begin{equation}
\pi = \left(\frac{p}{p_0}\right)^{\frac{R_dT}{cp}}
\end{equation}
To determine the saturation vapor mixing ratio the Teten's formula is used,
\begin{equation}
q_{vs}(p,T) = \left( \frac{380.0}{p} \right) \exp\left(17.27 \times \frac{T-273.0}{T-36.0}\right)
\end {equation}
The autoconvection rate ($A_r$) and collection rate ($C_r$) follow Kessler parametrization and are defined by:
\begin{align}
A_r &= k_1(q_c-a) \\
C_r &= k_2q_cq_r^{0.875}
\end{align} With $k_1=0.001 \text{s}^{-1}$, $a=0.001 \text{g}.\text{g}^{-1}$ and $k_2=2.2 \text{s}^{-1}$
Deriving from \cite{klemp1978simulation} description of cloud water,rain water and water vapor mixing ratios. they are define as followed:
\begin{equation}
q_c^{n+1}=\mbox{max}(q_c^r-\Delta q_r,0)
\end{equation}
\begin{equation}
q_r^{n+1}=\mbox{max}(q_r^r-\Delta q_r+S,0)
\end{equation} where $S$ is the sedimentation term and $\Delta q_r$ is defined as
\begin{equation}
\Delta q_r=q_c^n-\frac{q_c^n-\Delta \text{t}\ \mbox{max}(A_r,0)}{1+\Delta \text{t} C_r}
\end{equation}
The Rain evaporation equation is defined similarly to \cite{ogura1971numerical} description:
\begin{equation}
E_r=\frac{1}{\rho}\frac{\left(1-\frac{q_v}{q_{vs}}\right)C(\rho q_r)^{0.525}}{5.4\times10^5+\frac{2.55\times10^6}{pq_{vs}}}
\end{equation} With ventilation factor C define as
\begin{equation}
C_r=1.6+124.9(\rho q_r)^{0.2046}
\label{venti}
\end{equation}
The liquid water terminal velocity is similar to \cite{soong1973comparison} description with a mean density adjustment as suggested by \cite{kessler1969distribution}:
\begin{equation}
V_r = 36349(\rho q_r)^{0.1346}\left(\frac{\rho}{\rho_0}\right)^{-\frac{1}{2}}
\end{equation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{`Toy' Chemistry} \label{sec:ToyChemistry}
The toy chemistry module represents a simple photolysis-driven chemical reaction that incorporates combination and the dissociation of a chemical species:
\begin{align}
Cl + Cl &= Cl_2 && \mbox{(reaction rate $k_1$)}\\
Cl_2&=Cl+Cl && \mbox{(reaction rate $k_2$)}.
\end{align}
Observe that the total number of molecules of the chemical species are conserved in this reaction,
\begin{equation}
Cl_T=Cl+2 Cl_2.
\end{equation} Representing the mixing ratios of these species as $q_{Cl}$ and $q_{Cl2}$, we can define the total mixing ratio of Chlorine atoms as
\begin{equation}
q_{Cly} = q_{Cl} + 2 q_{Cl2}.
\end{equation}
The differential equations describing the evolution of $Cl$ and $Cl_2$ under this reaction take the form
\begin{align}
\frac{Dq_{Cl}}{Dt} &= 2k_1 q_{Cl2} -2k_2 q_{Cl}^2, \\
\frac{Dq_{Cl2}}{Dt} &= -k_1 q_{Cl2} + k_2 q_{Cl}^2,
\end{align} where $D/Dt$ denotes the Lagrangian derivative. Observe that the total mixing ratio of Chlorine atoms then satisfies
\begin{equation}
\frac{Dq_{Cly}}{Dt} = \frac{Dq_{Cl}}{Dt} + 2 \frac{Dq_{Cl2}}{Dt} = 0,
\end{equation} and so the total mixing ratio of Chlorine is held constant.
The two reaction rate coefficient $k_1$ and $k_2$, representing the the photolytic dissociation and recombination of Chlorine gas are defined as
\begin{align}
k_1(\lambda,\theta)&= \mbox{max} \left[ 0,\sin\theta\sin\theta_c+\cos\theta\cos\theta_c, \label{app:term_k1}
\cos(\lambda-\lambda_c) \right] \\
k_2(\lambda,\theta)&=1,\label{app:term_k2}
\end{align} where $(\lambda_c, \theta_c)=(20^\circ N, 300^\circ E)$ denote the sub-solar point on the Earth's surface.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Surface Fluxes on an Aqua-Planet with Prescribed Sea Surface Temperatures} \label{sec:OceanSurfaceFluxes}
The forcing by surface fluxes from an idealized ocean is described in \cite{reed2012idealized} and is partly reproduced here. We use a model configuration which corresponds to an aqua-planet setup with prescribed sea surface temperatures (SSTs). This forcing by the surface fluxes is applied to the state variables in the lowermost model level using a partially implicit formulation to avoid numerical instabilities. Throughout this section we use the subscript $a$ to denote variables defined on the lowermost model level.
The surface fluxes depend on the \textit{drag coefficient} $C_d$, defined as
\begin{equation} \label{eq:DragCoefficient}
\begin{array}{ll}
C_d = C_{d0} + C_{d1} \vert \vec{v}_a \vert & \mbox{for $\vert \vec{v}_a \vert < 20$ m s$^{-1}$} \\