access/SCM_testcases: make_rcde_nml.py

File make_rcde_nml.py, 5.9 KB (added by Robert Warren, 3 weeks ago)
Line 
1# Created by Rob Warren (rob.warren@monash.edu)
2
3# Creates Fortran 90 namelist for running UM Single Column Model in
4# radiative-convective equilibrium (RCE) or radiative-convective-dynamical
5# equilibrium (RCDE)
6
7# User must specifiy the output file name, a UM vertical levels namelist, the
8# run length and model time step, and four parameters relating to the initial
9# temperature, specific humidity, and vertical velocity profiles
10
11import f90nml as nml
12from f90nml.namelist import Namelist
13import numpy as np
14
15# Output namelist file
16nmlfile = 'namelist.scm'
17
18# Vertical levels namelist
19levsfile = 'vertlevs_L85_50t_35s_85km'
20
21# User-specified run parameters
22nday = 100       # run length (days)
23dt   = 1200.0    # model time step (s)
24
25# User-specified sounding parameters
26t0   = 300.0     # surface temperature (K)
27q0   = 18.65e-3  # surface specific humidity (kg/kg)
28wmax = 0.05      # maximum vertical velocity (m/s)
29hw   = 12500.0   # depth of vertical velocity profile (m)
30
31# Additional sounding parameters
32p0   = 101480.0  # surface pressure (Pa)
33zt   = 15000.0   # height of tropopause (m)
34zq1  = 4000.0    # height scale 1 for specific humidity profile (m)
35zq2  = 7500.0    # height scale 2 for specific humidity profile (m)
36lr   = 0.0067    # lapse rate (K/m)
37qt   = 1.0e-14   # specific humidity above tropopause (kg/kg)
38
39# Other parameters
40lon    = 0.0     # longitude (deg)
41lat    = 0.0     # latitude (deg)
42albedo = 0.07    # surface albedo
43co2    = 348     # co2 concentration (ppmv)
44
45# Physical constants
46rd = 287.04      # gas constant for dry air (J/kg/K)
47= 9.79764     # acceleration due to gravity (m2/s2)
48cp = 1004.64     # specific heat capacity for dry air (J/kg/K)
49
50########################################################################
51
52# Read in vertical levels namelist
53levs = nml.read(levsfile)
54ztop = levs['vertlevs']['z_top_of_model']
55eta_theta = np.array(levs['vertlevs']['eta_theta'])
56eta_rho = np.array(levs['vertlevs']['eta_rho'])
57
58# Convert eta levels to physical heights
59z_theta = ztop*eta_theta
60z_rho = ztop*eta_rho
61
62# Note the number of levels
63nz_theta = len(z_theta)
64nz_rho = len(z_rho)
65nlev = nz_rho
66
67# Get boolean indices of points in troposphere and stratosphere
68itr_theta = (z_theta < zt)
69ist_theta = (z_theta >= zt)
70itr_rho = (z_rho < zt)
71ist_rho = (z_rho >= zt)
72
73# Create specific humidity profile (theta levels)
74q = np.zeros(nz_theta)
75q[itr_theta] = q0*np.exp(-1*z_theta[itr_theta]/zq1)*np.exp(-1*(z_theta[itr_theta]/zq2)**2)
76q[ist_theta] = qt
77
78# Compute virtual temperature at surface and tropopause
79tv0 = t0*(1+0.608*q0)
80tvt = tv0-lr*zt
81
82# Create virtual temperature profile (theta levels)
83tv = np.zeros(nz_theta)
84tv[itr_theta] = tv0-lr*z_theta[itr_theta]
85tv[ist_theta] = tvt
86
87# Create temperature profile (theta levels)
88t = tv/(1+0.608*q)
89
90# Compute pressure at tropopause
91pt = p0*(tvt/tv0)**(g/(rd*lr))
92
93# Create pressure profile on theta levels
94p_theta = np.zeros(nz_theta)
95p_theta[itr_theta] = p0*((tv0-lr*z_theta[itr_theta])/tv0)**(g/(rd*lr))
96p_theta[ist_theta] = pt*np.exp(-1*(g*(z_theta[ist_theta]-zt)/(rd*tvt)))
97
98# Compute potential temperature profile
99th = t*(100000.0/p_theta)**(rd/cp)
100
101# Create pressure profile (rho levels)
102p = np.zeros(nz_rho)
103p[itr_rho] = p0*((tv0-lr*z_rho[itr_rho])/tv0)**(g/(rd*lr))
104p[ist_rho] = pt*np.exp(-1*(g*(z_rho[ist_rho]-zt)/(rd*tvt)))
105
106# Create vertical velocity profile (theta levels)
107w = wmax*np.sin(np.pi*(z_theta)/hw)
108w[z_theta == 0] = 0.0
109w[z_theta >= hw] = 0.0
110
111# Create ozone profile
112ozone = (1e-6)*3.6478*((p/100)**0.83209)*np.exp(-1*(p/100)/11.3515)   # ppmv
113ozone *= (48.00/28.96)   # kg/kg
114
115# Convert CO2 concentration to kg/kg
116co2 *= (1e-6)*(44.01/28.96)
117
118# Create wind profiles
119u = np.zeros(nz_rho)
120v = np.zeros(nz_rho)
121
122# Add extra value at top of pressure profile
123p = np.append(p, 0.0)
124
125# Set number of forcing times and forcing period
126nfor = 2
127obs_pd = (nday*86400.)/(nfor-1)
128
129# Create namelist
130new = Namelist({})
131
132cats = ['cntlscm',
133        'indata',
134        'rundata',
135        'logic',
136        'ingwd',
137        'injules',
138        'inprof',
139        'ingeofor',
140        'physwitch',
141        'radcloud',
142        'inobsfor']
143
144for cat in cats:
145#for ii in range(len(cats)):
146#    cat = cats[ii]
147    new[cat] = Namelist({})
148
149new['cntlscm']['model_levels_nml'] = nz_rho
150new['cntlscm']['land_points']      = 0
151new['cntlscm']['nfor']             = nfor
152
153new['indata']['lat']  = lon
154new['indata']['long'] = lat
155
156new['rundata']['ndayin']   = nday
157new['rundata']['nminin']   = 0
158new['rundata']['nsecin']   = 0
159new['rundata']['timestep'] = dt
160new['rundata']['sice_alb'] = albedo
161new['rundata']['co2start'] = co2
162new['rundata']['co2end']   = co2
163new['rundata']['ozone']    = list(ozone)
164
165new['logic']['obs']           = True
166new['logic']['obs_surf']      = True
167new['logic']['land_sea_mask'] = False
168new['logic']['land_ice_mask'] = False
169new['logic']['soil_mask']     = False
170new['logic']['local_time']    = False
171new['logic']['l_qpos_for']    = True
172
173new['inprof']['tstari'] = t0
174new['inprof']['kill_interp'] = True   # Model will crash without this!
175new['inprof']['z_tom_nml'] = ztop
176new['inprof']['p_in']   = list(p)
177new['inprof']['theta']  = list(th[1:nlev+1])   # exclude surface
178new['inprof']['qi']     = list(q[1:nlev+1])    # exclude surface
179new['inprof']['ui']     = list(u)
180new['inprof']['vi']     = list(v)
181new['inprof']['wi']     = list(w)
182new['inprof']['w_advi'] = list(w)
183
184new['inobsfor']['l_vertadv'] = True
185new['inobsfor']['obs_pd'] = obs_pd
186new['inobsfor']['obs_bot'] = 0.0
187new['inobsfor']['obs_top'] = ztop
188new['inobsfor']['rlx_t'] = 0 
189new['inobsfor']['rlx_q'] = 0
190new['inobsfor']['rlx_u'] = 0
191new['inobsfor']['rlx_v'] = 0
192new['inobsfor']['rlx_w'] = 3
193new['inobsfor']['tstar_forcing'] = list(np.zeros(nfor)+t0)
194new['inobsfor']['t_inc']  = list(np.zeros(nfor*nlev))
195new['inobsfor']['q_star'] = list(np.zeros(nfor*nlev))
196new['inobsfor']['u_inc']  = list(np.zeros(nfor*nlev))
197new['inobsfor']['v_inc']  = list(np.zeros(nfor*nlev))
198new['inobsfor']['w_inc']  = list(np.zeros(nfor*nlev))
199
200# Write namelist to file
201new.write(nmlfile,force=True)