-
Notifications
You must be signed in to change notification settings - Fork 12
/
Copy pathconstants.py
214 lines (194 loc) · 7.81 KB
/
constants.py
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
import math
import os
from enum import Enum
from ndsl.logging import ndsl_log
# The FV3GFS model ships with two sets of constants, one used in the GFS physics
# package and the other used for the Dycore. Their difference are small but significant
# In addition the GSFC's GEOS model as its own variables
class ConstantVersions(Enum):
GFDL = "GFDL" # NOAA's FV3 dynamical core constants (original port)
GFS = "GFS" # Constant as defined in NOAA GFS
GEOS = "GEOS" # Constant as defined in GEOS v13
SHiELD = "SHiELD" # Constant as defined in NOAA SHiELD
CONST_VERSION_AS_STR = os.environ.get("PACE_CONSTANTS", "GFS")
try:
CONST_VERSION = ConstantVersions[CONST_VERSION_AS_STR]
ndsl_log.info(f"Constant selected: {CONST_VERSION}")
except KeyError as e:
raise RuntimeError(f"Constants {CONST_VERSION_AS_STR} is not implemented, abort.")
#####################
# Common constants
#####################
ROOT_RANK = 0
X_DIM = "x"
X_INTERFACE_DIM = "x_interface"
Y_DIM = "y"
Y_INTERFACE_DIM = "y_interface"
Z_DIM = "z"
Z_INTERFACE_DIM = "z_interface"
Z_SOIL_DIM = "z_soil"
TILE_DIM = "tile"
X_DIMS = (X_DIM, X_INTERFACE_DIM)
Y_DIMS = (Y_DIM, Y_INTERFACE_DIM)
Z_DIMS = (Z_DIM, Z_INTERFACE_DIM)
HORIZONTAL_DIMS = X_DIMS + Y_DIMS
INTERFACE_DIMS = (X_INTERFACE_DIM, Y_INTERFACE_DIM, Z_INTERFACE_DIM)
SPATIAL_DIMS = X_DIMS + Y_DIMS + Z_DIMS
WEST = 0
EAST = 1
NORTH = 2
SOUTH = 3
NORTHWEST = 4
NORTHEAST = 5
SOUTHWEST = 6
SOUTHEAST = 7
INTERIOR = 8
EDGE_BOUNDARY_TYPES = (NORTH, SOUTH, WEST, EAST)
CORNER_BOUNDARY_TYPES = (NORTHWEST, NORTHEAST, SOUTHWEST, SOUTHEAST)
BOUNDARY_TYPES = EDGE_BOUNDARY_TYPES + CORNER_BOUNDARY_TYPES
N_HALO_DEFAULT = 3
#######################
# Tracers configuration
#######################
# nq is actually given by ncnst - pnats, where those are given in atmosphere.F90 by:
# ncnst = Atm(mytile)%ncnst
# pnats = Atm(mytile)%flagstruct%pnats
# here we hard-coded it because 8 is the only supported value, refactor this later!
if CONST_VERSION == ConstantVersions.GEOS:
# 'qlcd' is exchanged in GEOS
NQ = 9
elif (
CONST_VERSION == ConstantVersions.GFS
or CONST_VERSION == ConstantVersions.SHiELD
or CONST_VERSION == ConstantVersions.GFDL
):
NQ = 8
else:
raise RuntimeError("Constant selector failed, bad code.")
#####################
# Physical constants
#####################
if CONST_VERSION == ConstantVersions.GEOS:
RADIUS = 6.371e6
PI = 3.14159265358979323846
OMEGA = 2.0 * PI / 86164.0
GRAV = 9.80665
RGRAV = 1.0 / GRAV
RDGAS = 8314.47 / 28.965
RVGAS = 8314.47 / 18.015
HLV = 2.4665e6
HLF = 3.3370e5
KAPPA = RDGAS / (3.5 * RDGAS)
CP_AIR = RDGAS / KAPPA
TFREEZE = 273.15
SAT_ADJUST_THRESHOLD = 1.0e-6
elif (
CONST_VERSION == ConstantVersions.GFS
or CONST_VERSION == ConstantVersions.SHiELD
):
RADIUS = 6.3712e6 # Radius of the Earth [m]
PI = 3.1415926535897931
OMEGA = 7.2921e-5 # Rotation of the earth
GRAV = 9.80665 # Acceleration due to gravity [m/s^2].04
RGRAV = 1.0 / GRAV # Inverse of gravitational acceleration
RDGAS = 287.05 # Gas constant for dry air [J/kg/deg] # 287.04
RVGAS = 461.50 # Gas constant for water vapor [J/kg/deg]
HLV = 2.5e6 # Latent heat of evaporation [J/kg]
HLF = 3.3358e5 # Latent heat of fusion [J/kg] # 3.34e5
CP_AIR = 1004.6
KAPPA = RDGAS / CP_AIR # Specific heat capacity of dry air at
TFREEZE = 273.15
SAT_ADJUST_THRESHOLD = 1.0e-8
elif CONST_VERSION == ConstantVersions.GFDL:
RADIUS = 6371.0e3 # Radius of the Earth [m] #6371.0e3
PI = 3.14159265358979323846 # 3.14159265358979323846
OMEGA = 7.292e-5 # Rotation of the earth # 7.292e-5
GRAV = 9.80 # Acceleration due to gravity [m/s^2].04
RGRAV = 1.0 / GRAV # Inverse of gravitational acceleration
RDGAS = 287.04 # Gas constant for dry air [J/kg/deg] # 287.04
RVGAS = 461.50 # Gas constant for water vapor [J/kg/deg]
HLV = 2.500e6 # Latent heat of evaporation [J/kg]
HLF = 3.34e5 # Latent heat of fusion [J/kg] # 3.34e5
KAPPA = 2.0 / 7.0
CP_AIR = RDGAS / KAPPA # Specific heat capacity of dry air at
TFREEZE = 273.16 # Freezing temperature of fresh water [K]
SAT_ADJUST_THRESHOLD = 1.0e-8
else:
raise RuntimeError("Constant selector failed, bad code.")
DZ_MIN = 2.0
CV_AIR = CP_AIR - RDGAS # Heat capacity of dry air at constant volume
RDG = -RDGAS / GRAV
CNST_0P20 = 0.2
K1K = RDGAS / CV_AIR
CNST_0P20 = 0.2
CV_VAP = 3.0 * RVGAS # Heat capacity of water vapor at constant volume
ZVIR = RVGAS / RDGAS - 1 # con_fvirt in Fortran physics
TICE = 273.16 # Freezing temperature
TICE0 = TICE - 0.01
CP_VAP = 4.0 * RVGAS # Heat capacity of water vapor at constant pressure
if CONST_VERSION == ConstantVersions.SHiELD:
C_ICE = 2.106e3 # Heat capacity of ice at 0 degrees Celsius
C_LIQ = 4.218e3 # Heat capacity of water at 0 degrees Celsius
DC_ICE = C_LIQ - C_ICE # Isobaric heating / cooling (J/kg/K)
DC_VAP = CP_VAP - C_LIQ # Isobaric heating / cooling (J/kg/K)
D2ICE = DC_VAP + DC_ICE # Isobaric heating / cooling (J/kg/K)
LV0 = (
HLV - DC_VAP * TICE0
) # 3148711.3338762247, evaporation latent heat coefficient at 0 degrees Kelvin
LI00 = (
HLF - DC_ICE * TICE0
) # -242413.92000000004, fusion latent heat coefficient at 0 degrees Kelvin
LI2 = (
LV0 + LI00
) # 2906297.413876225, sublimation latent heat coefficient at 0 degrees Kelvin
LI0 = HLF - DC_ICE * TICE0
else:
C_ICE = 1972.0 # Heat capacity of ice at -15 degrees Celsius
C_LIQ = 4.1855e3 # Heat capacity of water at 15 degrees Celsius
DC_ICE = C_LIQ - C_ICE # Isobaric heating / cooling (J/kg/K)
DC_VAP = CP_VAP - C_LIQ # Isobaric heating / cooling (J/kg/K)
D2ICE = DC_VAP + DC_ICE # Isobaric heating / cooling (J/kg/K)
LV0 = (
HLV - DC_VAP * TICE
) # 3.13905782e6, evaporation latent heat coefficient at 0 degrees Kelvin
LI00 = (
HLF - DC_ICE * TICE
) # -2.7105966e5, fusion latent heat coefficient at 0 degrees Kelvin
LI2 = (
LV0 + LI00
) # 2.86799816e6, sublimation latent heat coefficient at 0 degrees Kelvin
LI0 = HLF - DC_ICE * TICE
EPS = RDGAS / RVGAS
E00 = 611.21 # Saturation vapor pressure at 0 degrees Celsius
T_WFR = TICE - 40.0 # homogeneous freezing temperature
TICE0 = TICE - 0.01
T_MIN = 178.0 # Minimum temperature to freeze-dry all water vapor
T_SAT_MIN = TICE - 160.0
LAT2 = (HLV + HLF) ** 2 # used in bigg mechanism
RHO_0 = 1.0 # reference air density (kg/m^3), ref: IFS
RHO_W = 1.0e3 # density of cloud water (kg/m^3)
RHO_I = 9.17e2 # density of cloud ice (kg/m^3)
RHO_R = 1.0e3 # density of rain (Lin et al. 1983) (kg/m^3)
RHO_S = 1.0e2 # density of snow (Lin et al. 1983) (kg/m^3)
RHO_G = 4.0e2 # density of graupel (Rutledge and Hobbs 1984) (kg/m^3)
RHO_H = 9.17e2 # density of hail (Lin et al. 1983) (kg/m^3)
VISD = 1.717e-5 # dynamics viscosity of air at 0 deg C and 1000 hPa
# (Mason, 1971) (kg/m/s)
VISK = 1.35e-5 # kinematic viscosity of air at 0 deg C and 1000 hPa
# (Mason, 1971) (m^2/s)
VDIFU = 2.25e-5 # diffusivity of water vapor in air at 0 deg C and 1000 hPa
# (Mason, 1971) (m^2/s)
TCOND = 2.40e-2 # thermal conductivity of air at 0 C and 1000 hPa
# (Mason, 1971) (J/m/s/K)
SCM3 = math.exp(1.0 / 3 * math.log(VISK / VDIFU)) # Schmidt number, Sc ** (1 / 3)
# Lin et al. (1983)
QCMIN = 1.0e-15 # min value for cloud condensates (kg/kg)
QFMIN = 1.0e-8 # min value for sedimentation (kg/kg)
DT_FR = 8.0 # t_wfr - dt_fr: minimum temperature water can exist
# (Moore and Molinero 2011)
CDG = 3.15121 # drag coefficient of graupel (Locatelli and Hobbs, 1974)
CDH = 0.5 # drag coefficient of hail (Heymsfield and Wright, 2014)
DZ_MIN_FLIP = 1.0e-2 # used for correcting flipped height (m)
# Terminal Velocity Parameters, Lin et al. (1983)
GCON = (4.0 * GRAV * RHO_G / (3.0 * CDG * RHO_0)) ** 0.5
HCON = (4.0 * GRAV * RHO_H / (3.0 * CDH * RHO_0)) ** 0.5