forked from pressel/pycles
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathTracers.pyx
211 lines (171 loc) · 8.56 KB
/
Tracers.pyx
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
#!python
#cython: boundscheck=False
#cython: wraparound=False
#cython: initializedcheck=False
#cython: cdivision=True
cimport Grid
cimport ReferenceState
cimport PrognosticVariables
cimport DiagnosticVariables
cimport ParallelMPI
from NetCDFIO cimport NetCDFIO_Stats
import cython
from libc.math cimport fmax
cimport numpy as np
import numpy as np
include "parameters.pxi"
import cython
cdef extern from "thermodynamic_functions.h":
double pv_c(double p0, double qt, double qv) nogil
def TracersFactory(namelist):
try:
use_tracers = namelist['tracers']['use_tracers']
except:
use_tracers = False
if use_tracers:
return UpdraftTracers(namelist)
else:
return TracersNone()
cdef class TracersNone:
def __init__(self):
return
cpdef initialize(self, Grid.Grid Gr, PrognosticVariables.PrognosticVariables PV, NetCDFIO_Stats NS, ParallelMPI.ParallelMPI Pa):
return
cpdef update(self, Grid.Grid Gr, ReferenceState.ReferenceState Ref, PrognosticVariables.PrognosticVariables PV,
DiagnosticVariables.DiagnosticVariables DV,ParallelMPI.ParallelMPI Pa):
return
cpdef update_cleanup(self, Grid.Grid Gr, ReferenceState.ReferenceState Ref, PrognosticVariables.PrognosticVariables PV,
DiagnosticVariables.DiagnosticVariables DV,ParallelMPI.ParallelMPI Pa):
return
cpdef stats_io(self, Grid.Grid Gr, NetCDFIO_Stats NS, ParallelMPI.ParallelMPI Pa):
return
cdef class UpdraftTracers:
def __init__(self, namelist):
if namelist['microphysics']['scheme'] == 'None_SA' or namelist['microphysics']['scheme'] == 'SB_Liquid':
self.lcl_tracers = True
else:
self.lcl_tracers = False
self.index_lcl = 0
return
cpdef initialize(self, Grid.Grid Gr, PrognosticVariables.PrognosticVariables PV, NetCDFIO_Stats NS, ParallelMPI.ParallelMPI Pa):
# Assemble a dictionary with the tracer information
# Can be expanded for different init heights or timescales
self.tracer_dict = {}
self.tracer_dict['surface'] = {}
self.tracer_dict['surface']['c_srf_15'] = {}
self.tracer_dict['surface']['c_srf_15']['timescale'] = 15.0 * 60.0
if self.lcl_tracers:
self.tracer_dict['lcl'] = {}
self.tracer_dict['lcl']['c_lcl_15'] = {}
self.tracer_dict['lcl']['c_lcl_15']['timescale'] = 15.0 * 60.0
for var in self.tracer_dict['surface'].keys():
PV.add_variable(var, '-', var, 'tracer diagnostics' , "sym", "scalar", Pa)
if self.lcl_tracers:
for var in self.tracer_dict['lcl'].keys():
PV.add_variable(var, '-', var, 'tracer diagnostics', "sym", "scalar", Pa)
NS.add_ts('grid_lcl', Gr, Pa )
return
cpdef update(self, Grid.Grid Gr, ReferenceState.ReferenceState Ref, PrognosticVariables.PrognosticVariables PV,
DiagnosticVariables.DiagnosticVariables DV,ParallelMPI.ParallelMPI Pa):
cdef:
Py_ssize_t istride = Gr.dims.nlg[1] * Gr.dims.nlg[2]
Py_ssize_t jstride = Gr.dims.nlg[2]
Py_ssize_t i,j,k,ishift,jshift,ijk
Py_ssize_t var_shift
double tau
# Set the source term
for level in self.tracer_dict.keys():
for var in self.tracer_dict[level].keys():
var_shift = PV.get_varshift(Gr, var)
tau = self.tracer_dict[level][var]['timescale']
with nogil:
for i in xrange(Gr.dims.npg):
PV.tendencies[var_shift + i] += -fmax(PV.values[var_shift + i],0.0)/tau
#Below we assume domain uses only x-y decomposition, for the general case we should use pencils
#but we don't address that complication for now
# Set the value of the surface based tracers
for var in self.tracer_dict['surface'].keys():
var_shift = PV.get_varshift(Gr, var)
with nogil:
for i in xrange(Gr.dims.nlg[0]):
for j in xrange(Gr.dims.nlg[1]):
ijk = i * istride + j * jstride + Gr.dims.gw
PV.values[var_shift + ijk] = 1.0
# Find surface thermodynamic properties to compute LCL (lifting condensation level)
# Ref: M. G. Lawrence, "The relationship between relative humidity and the dew point
# temperature in moist air: A simple conversion and applications", Bull. Am. Meteorol. Soc., 86, 225-233, 2005
cdef:
double [:] T_dew = np.zeros(Gr.dims.nl[0]*Gr.dims.nl[1],dtype=np.double,order='c')
double [:] z_lcl = np.zeros(Gr.dims.nl[0]*Gr.dims.nl[1],dtype=np.double,order='c')
Py_ssize_t t_shift, qv_shift, qt_shift
Py_ssize_t count = 0
double B1 = 243.04 # C
double A1 = 17.625 # no units
double C1 = 610.94 # Pa
double CtoK = 273.15 # additive celsius to kelvin conversion factor
double vapor_pressure
double nxny_i =1.0/( Gr.dims.n[0] * Gr.dims.n[1])
double lcl_, lcl
if self.lcl_tracers:
t_shift = DV.get_varshift(Gr,'temperature')
qt_shift = PV.get_varshift(Gr, 'qt')
qv_shift = DV.get_varshift(Gr,'qv')
with nogil:
for i in xrange(Gr.dims.gw, Gr.dims.nlg[0]-Gr.dims.gw):
for j in xrange(Gr.dims.gw, Gr.dims.nlg[1]-Gr.dims.gw):
ijk = i * istride + j * jstride + Gr.dims.gw
vapor_pressure = pv_c(Ref.p0_half[Gr.dims.gw], PV.values[qt_shift + ijk], DV.values[qv_shift + ijk])
T_dew[count] = B1 * log(vapor_pressure/C1)/(A1-log(vapor_pressure/C1)) + CtoK
z_lcl[count] = 125.0 * (DV.values[t_shift+ijk] - T_dew[count])
count += 1
lcl_ = np.sum(z_lcl) * nxny_i
lcl = Pa.domain_scalar_sum(lcl_)
# Pa.root_print('LCL is ' + str(lcl))
for k in xrange(Gr.dims.nlg[2]-Gr.dims.gw, Gr.dims.gw-1, -1):
if Gr.zl_half[k] <= lcl:
self.index_lcl = k
break
# Pa.root_print('Grid LCL ' + str(Gr.zl_half[self.index_lcl]))
for var in self.tracer_dict['lcl'].keys():
var_shift = PV.get_varshift(Gr, var)
with nogil:
for i in xrange( Gr.dims.nlg[0]):
for j in xrange( Gr.dims.nlg[1]):
ijk = i * istride + j * jstride + self.index_lcl
PV.values[var_shift + ijk] = 1.0
return
cpdef update_cleanup(self, Grid.Grid Gr, ReferenceState.ReferenceState Ref, PrognosticVariables.PrognosticVariables PV,
DiagnosticVariables.DiagnosticVariables DV,ParallelMPI.ParallelMPI Pa):
cdef:
Py_ssize_t istride = Gr.dims.nlg[1] * Gr.dims.nlg[2]
Py_ssize_t jstride = Gr.dims.nlg[2]
Py_ssize_t i,j,k,ishift,jshift,ijk
Py_ssize_t var_shift
#Below we assume domain uses only x-y decomposition, for the general case we should use pencils
#but we don't address that complication for now
# Set the value of the surface based tracers
for var in self.tracer_dict['surface'].keys():
var_shift = PV.get_varshift(Gr, var)
with nogil:
for i in xrange(Gr.dims.nlg[0]):
for j in xrange(Gr.dims.nlg[1]):
ijk = i * istride + j * jstride + Gr.dims.gw
PV.tendencies[var_shift + ijk] = 0.0
for k in xrange( Gr.dims.nlg[2]):
ijk = i * istride + j * jstride + k
PV.values[var_shift + ijk] = fmax(PV.values[var_shift + ijk],0.0)
for var in self.tracer_dict['lcl'].keys():
var_shift = PV.get_varshift(Gr, var)
with nogil:
for i in xrange( Gr.dims.nlg[0]):
for j in xrange( Gr.dims.nlg[1]):
ijk = i * istride + j * jstride + self.index_lcl
PV.tendencies[var_shift + ijk] = 0.0
for k in xrange( Gr.dims.nlg[2]):
ijk = i * istride + j * jstride + k
PV.values[var_shift + ijk] = fmax(PV.values[var_shift + ijk],0.0)
return
cpdef stats_io(self, Grid.Grid Gr, NetCDFIO_Stats NS, ParallelMPI.ParallelMPI Pa):
if self.lcl_tracers:
NS.write_ts('grid_lcl',Gr.zl_half[self.index_lcl], Pa)
return