-
Notifications
You must be signed in to change notification settings - Fork 7
PreparingInputFiles
Prev: Example two-stream flux calculations in cloudy sky
This tutorial supplements the information given on producing input files in chapter 5 of the Socrates User Guide (pdf can be downloaded from the latest release).
The radiation code can accept input in either netCDF format (using the l_run_cdf driver) or the ASCII equivalent CDL format (using the l_run_cdl) driver. There are utilities available in the code to convert between these formats with options given in the man pages:
man Ccdf2cdl
man Ccdl2cdf
The ''cdfdump'' utility will also dump the CDL version of a netCDF file to the screen.
A number of python utilities are available to create netCDF files in the required format from arrays of data. These can be found in $RAD_DIR/python/. Further utilities written in Fortran are available to manipulate profiles. These are called from the scripts found in sbin/ and generally work with CDL format files.
An example utility called ncprofiles can be found in the python/ directory that creates many of the files needed by the radiation code. This requires an existing .t profile and will then create other profiles using the same level set. As an example we can copy the tro.t McClatchey profile to our example directory (and convert to netCDF). The McClatchey profiles are standard reference atmospheric profiles taken from the book: Optical Properties of the Atmosphere, 1972 by R. A. McClatchey et al.
mkdir $RAD_DIR/examples/example3
cd $RAD_DIR/examples/example3
cp $RAD_DATA/mcc_profiles/one_km/tro.t cdl_tro.t
Ccdl2cdf -o tro.t cdl_tro.t
Then run the python ncprofiles.py script:
cp tro.t tro_python.t
python $RAD_DIR/python/ncprofiles.py tro_python
The ncprofiles script calls various utilities from the nctools package to create the individual profiles. The full list of utilities is given in the user guide (section 5.1). These provide an easy way to create all the required profiles and allows data manipulation to be done in python.
For UM (or other GCM) users, there is currently no standard way to convert atmospheric profiles from UM "STASH" output to Socrates input files. The best way is to import the required arrays into python and then create the input files using these tools. Note that the Socrates code (and the python nctools) can handle multicolumn data, although you may encounter memory limitations if you attempt to solve for very large arrays. GCMs (including the UM) often do not use pressure as a vertical coordinate and require multicolumn arrays of pressure on layers/levels. These can be supplied to the Socrates code as extra input files with the extensions .p and .pl (equivalent to the .t and .tl temperature files). If these are present the vertical coordinate given in the files becomes a dummy field and is not used.
This Fortran routine takes an input file of ASCII data in columns and produces the required files representing each column in CDL format (which can then be converted to netCDF with Ccdl2cdf). In the input file, columns of data are given between *PROFILE and *END directives. At the head of each column the variable name is given followed by the units in brackets. The variable names and units recognized by the code are given in the appendix to the Socrates User Guide. This code will perform a unit conversion from the units given in the header to the SI units required by the radiative transfer code.
Examples of using the raw_input routine are given in examples/raw_input/run_me and examples/prsc/. The first example converts trace gas amounts in parts-per-million by volume to the units of kg/kg (or "none") required by the code. To view the input file:
cd $RAD_DIR/examples/raw_input
less atm_profile
Then run raw_input
raw_input
answering the prompts:
atm_profile
mls
n
The output files are named mls_1.* The _1 extension is added to distinguish this as the first *PROFILE directive given in the input file. The input file can contain as many *PROFILE directives as required.
There are a number of Fortran utilities called from scripts in sbin/ that may be used to create and manipulate the input files. Some are documented in the user guide and man pages, although most are not:
man Cinterp
To find out more about these utilities currently requires looking directly at the scripts in sbin/ and then delving further into the Fortran routines themselves.
An example script in examples/setup_atmos/run_me demonstrates the use of some of these utilities. We'll run through this here. (Note: a slightly different version of this example is given in the user guide section 7.1.)
mkdir $RAD_DIR/examples/example4
cd $RAD_DIR/examples/example4
First we copy in temperature fields from the standard mid-latitude summer McClatchey profiles:
ATM=mls
BASE=wk_$ATM
cp $RAD_DATA/mcc_profiles/one_km/$ATM.t $BASE.tl
cp $RAD_DATA/mcc_profiles/one_km/$ATM.tstar $BASE.tstar
The raw McClatchey profile of temperature is used to define the edges of the atmospheric layers (suffix .tl). Other atmospheric quantities are defined at the midpoints of layers, so we make the appropriate mid-points:
Cmid_point -o $BASE.null $BASE.tl
We have now defined the grid. We can interpolate the temperatures at the mid-points of layers (suffix .t) from the edge temperatures. The best option for interpolation appears to be linear interpolation of the temperature with the logarithm of the pressure (option -lgn).
Cinterp -g $BASE.null -o $BASE.t -n "t" -u "K" -L "Central Temperatures" -lgn $BASE.tl
Interpolate the specific humidity to these levels. For ozone and water vapour interpolation of the log of the specific humidity in the log of the pressure (option -lgg) seems to perform best.
Cinterp -g $BASE.null -o $BASE.q -n "q" -u "None" -L "Specific humidity" \
-lgg $RAD_DATA/mcc_profiles/one_km/$ATM.q
Cinterp -g $BASE.null -o $BASE.o3 -n "o3" -u "None" -L "Ozone mmr" \
-lgg $RAD_DATA/mcc_profiles/one_km/$ATM.o3
For other gases we can add a constant mixing ratio by incrementing the null field. For CO2 we add a value of 5.241e-4 kg/kg for all pressures between 0.0 and 1.2e5 Pa:
Cinc_field -R 0.0,1.2e5:5.241e-4 -o $BASE.co2 \
-n "co2" -u "None" -L "CO2 mmr" $BASE.null
Cinc_field -R 0.0,1.2e5:0.2314 -o $BASE.o2 \
-n "o2" -u "None" -L "O2 mmr" $BASE.null
Surface albedos are required. A value of 0.06 is appropriate for the SW for an ocean surface (though is used for both SW and LW in this example). The -N and -T options just specify the arbitrary longitude and latitude respectively:
Cgen_surf_cdl -o $BASE.surf -n alb \
-L "Surface Albedo" -u "None" \
-b 0.06 -N 0.0 -T 0.0
A file of solar zenith angles is required, along with the incident solar flux at the top-of-the atmosphere. Fields which have no vertical structure can be generated using gen_horiz_cdl:
Cgen_horiz_cdl -o $BASE.szen -n szen \
-L "Solar zenith angle" -u "Degrees" \
-F 10.0 -N 0.0 -T 0.0
Cgen_horiz_cdl -o $BASE.stoa -n stoa \
-L "Solar Irradiance" -u "W.m-2" \
-F 1365.0 -N 0.0 -T 0.0
Some simple ice cloud fields are now generated. We set a constant ice mass fraction and effective dimension to be used wherever there is cloud:
IWC=0.000025
Cinc_field -R 0.0,1.2e5:$IWC -o $BASE.iwm -n "iwm" -u "kg/kg" \
-L "Ice Water Content" $BASE.null
IRE=0.000010
Cinc_field -R 0.0,1.2e5:$IRE -o $BASE.ire -n "dl" -u "m" \
-L "Mean maximum dimension" $BASE.null
The cloud fraction is set to one between given cloud base and top pressures:
P_CLTOP=5.0e4
P_CLBASE=6.0e4
Cinc_field -R $P_CLTOP,$P_CLBASE:1.0 -o $BASE.clfr \
-n "clfr" -u "None" -L "Cloud Fraction" $BASE.null
We now have the full set of input fields needed for this example. The examples/setup_atmos/run_me file goes on to run a number of different SW and LW calculations using these input profiles. One of these compares the fluxes calculated using the two-stream and spherical harmonics methods:
echo "Calculating SW fluxes using spherical harmonics"
Cl_run_cdl -B $BASE -s $RAD_DATA/spectra/ga7/sp_sw_ga7 -R 1 6 \
-S -r -g 2 -c -C 3 -K 1 -i 8 +R \
+S 3 7 0 0 -G 5 0 -H -T +F
fmove wk_$ATM sw_sph_cloudy_$ATM
echo "Calculating SW fluxes using two-stream approximation"
Cl_run_cdl -B $BASE -s $RAD_DATA/spectra/ga7/sp_sw_ga7 -R 1 6 \
-S -r -g 2 -c -C 2 -K 1 -i 8 +R
fmove wk_$ATM sw_ts_cloudy_$ATM
The heating rates from each method can then be compared:
ncplot sw_sph_cloudy_$ATM.hrts &
ncplot sw_ts_cloudy_$ATM.hrts &
Ccdl2cdf sw_sph_cloudy_$ATM.hrts
Ccdl2cdf sw_ts_cloudy_$ATM.hrts
ncdiff nc_sw_sph_cloudy_$ATM.hrts nc_sw_ts_cloudy_$ATM.hrts nc_sw_diff_cloudy_$ATM.hrts
ncplot nc_sw_diff_cloudy_$ATM.hrts &