Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_ra_cam_support.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 143.0 KB

Line
1	MODULE module_ra_cam_support
2	use module_cam_support, only: endrun
3	implicit none
4	integer, parameter :: r8 = 8
5	real(r8), parameter:: inf = 1.e20 ! CAM sets this differently in infnan.F90
6	integer, parameter:: bigint = O'17777777777' ! largest possible 32-bit integer
7
8	integer :: ixcldliq
9	integer :: ixcldice
10	! integer :: levsiz ! size of level dimension on dataset
11	integer, parameter :: nbands = 2 ! Number of spectral bands
12	integer, parameter :: naer_all = 12 + 1
13	integer, parameter :: naer = 10 + 1
14	integer, parameter :: bnd_nbr_LW=7
15	integer, parameter :: ndstsz = 4 ! number of dust size bins
16	integer :: idxSUL
17	integer :: idxSSLT
18	integer :: idxDUSTfirst
19	integer :: idxCARBONfirst
20	integer :: idxOCPHO
21	integer :: idxBCPHO
22	integer :: idxOCPHI
23	integer :: idxBCPHI
24	integer :: idxBG
25	integer :: idxVOLC
26
27	integer :: mxaerl ! Maximum level of background aerosol
28
29	! indices to sections of array that represent
30	! groups of aerosols
31
32	integer, parameter :: &
33	numDUST = 4, &
34	numCARBON = 4
35
36	! portion of each species group to use in computation
37	! of relative radiative forcing.
38
39	real(r8) :: sulscl_rf = 0._r8 !
40	real(r8) :: carscl_rf = 0._r8
41	real(r8) :: ssltscl_rf = 0._r8
42	real(r8) :: dustscl_rf = 0._r8
43	real(r8) :: bgscl_rf = 0._r8
44	real(r8) :: volcscl_rf = 0._r8
45
46	! "background" aerosol species mmr.
47	real(r8) :: tauback = 0._r8
48
49	! portion of each species group to use in computation
50	! of aerosol forcing in driving the climate
51	real(r8) :: sulscl = 1._r8
52	real(r8) :: carscl = 1._r8
53	real(r8) :: ssltscl = 1._r8
54	real(r8) :: dustscl = 1._r8
55	real(r8) :: volcscl = 1._r8
56
57	!From volcrad.F90 module
58	integer, parameter :: idx_LW_0500_0650=3
59	integer, parameter :: idx_LW_0650_0800=4
60	integer, parameter :: idx_LW_0800_1000=5
61	integer, parameter :: idx_LW_1000_1200=6
62	integer, parameter :: idx_LW_1200_2000=7
63
64	! First two values represent the overlap of volcanics with the non-window
65	! (0-800, 1200-2200 cm^-1) and window (800-1200 cm^-1) regions.\| Coefficients
66	! were derived using crm_volc_minimize.pro with spectral flux optimization
67	! on first iteration, total heating rate on subsequent iterations (2-9).
68	! Five profiles for HLS, HLW, MLS, MLW, and TRO conditions were given equal
69	! weight. RMS heating rate errors for a visible stratospheric optical
70	! depth of 1.0 are 0.02948 K/day.
71	!
72	real(r8) :: abs_cff_mss_aer(bnd_nbr_LW) = &
73	(/ 70.257384, 285.282943, &
74	1.0273851e+02, 6.3073303e+01, 1.2039569e+02, &
75	3.6343643e+02, 2.7138528e+02 /)
76
77	!From radae.F90 module
78	real(r8), parameter:: min_tp_h2o = 160.0 ! min T_p for pre-calculated abs/emis
79	real(r8), parameter:: max_tp_h2o = 349.999999 ! max T_p for pre-calculated abs/emis
80	real(r8), parameter:: dtp_h2o = 21.111111111111 ! difference in adjacent elements of tp_h2o
81	real(r8), parameter:: min_te_h2o = -120.0 ! min T_e-T_p for pre-calculated abs/emis
82	real(r8), parameter:: max_te_h2o = 79.999999 ! max T_e-T_p for pre-calculated abs/emis
83	real(r8), parameter:: dte_h2o = 10.0 ! difference in adjacent elements of te_h2o
84	real(r8), parameter:: min_rh_h2o = 0.0 ! min RH for pre-calculated abs/emis
85	real(r8), parameter:: max_rh_h2o = 1.19999999 ! max RH for pre-calculated abs/emis
86	real(r8), parameter:: drh_h2o = 0.2 ! difference in adjacent elements of RH
87	real(r8), parameter:: min_lu_h2o = -8.0 ! min log_10(U) for pre-calculated abs/emis
88	real(r8), parameter:: min_u_h2o = 1.0e-8 ! min pressure-weighted path-length
89	real(r8), parameter:: max_lu_h2o = 3.9999999 ! max log_10(U) for pre-calculated abs/emis
90	real(r8), parameter:: dlu_h2o = 0.5 ! difference in adjacent elements of lu_h2o
91	real(r8), parameter:: min_lp_h2o = -3.0 ! min log_10(P) for pre-calculated abs/emis
92	real(r8), parameter:: min_p_h2o = 1.0e-3 ! min log_10(P) for pre-calculated abs/emis
93	real(r8), parameter:: max_lp_h2o = -0.0000001 ! max log_10(P) for pre-calculated abs/emis
94	real(r8), parameter:: dlp_h2o = 0.3333333333333 ! difference in adjacent elements of lp_h2o
95	integer, parameter :: n_u = 25 ! Number of U in abs/emis tables
96	integer, parameter :: n_p = 10 ! Number of P in abs/emis tables
97	integer, parameter :: n_tp = 10 ! Number of T_p in abs/emis tables
98	integer, parameter :: n_te = 21 ! Number of T_e in abs/emis tables
99	integer, parameter :: n_rh = 7 ! Number of RH in abs/emis tables
100	real(r8):: c16,c17,c26,c27,c28,c29,c30,c31
101	real(r8):: fwcoef ! Farwing correction constant
102	real(r8):: fwc1,fwc2 ! Farwing correction constants
103	real(r8):: fc1 ! Farwing correction constant
104	real(r8):: amco2 ! Molecular weight of co2 (g/mol)
105	real(r8):: amd ! Molecular weight of dry air (g/mol)
106	real(r8):: p0 ! Standard pressure (dynes/cm**2)
107
108	! These are now allocatable. JM 20090612
109	real(r8), allocatable, dimension(:,:,:,:,:) :: ah2onw ! (n_p, n_tp, n_u, n_te, n_rh) ! absorptivity (non-window)
110	real(r8), allocatable, dimension(:,:,:,:,:) :: eh2onw ! (n_p, n_tp, n_u, n_te, n_rh) ! emissivity (non-window)
111	real(r8), allocatable, dimension(:,:,:,:,:) :: ah2ow ! (n_p, n_tp, n_u, n_te, n_rh) ! absorptivity (window, for adjacent layers)
112	real(r8), allocatable, dimension(:,:,:,:,:) :: cn_ah2ow ! (n_p, n_tp, n_u, n_te, n_rh) ! continuum transmission for absorptivity (window)
113	real(r8), allocatable, dimension(:,:,:,:,:) :: cn_eh2ow ! (n_p, n_tp, n_u, n_te, n_rh) ! continuum transmission for emissivity (window)
114	real(r8), allocatable, dimension(:,:,:,:,:) :: ln_ah2ow ! (n_p, n_tp, n_u, n_te, n_rh) ! line-only transmission for absorptivity (window)
115	real(r8), allocatable, dimension(:,:,:,:,:) :: ln_eh2ow ! (n_p, n_tp, n_u, n_te, n_rh) ! line-only transmission for emissivity (window)
116
117	!
118	! Constant coefficients for water vapor overlap with trace gases.
119	! Reference: Ramanathan, V. and P.Downey, 1986: A Nonisothermal
120	! Emissivity and Absorptivity Formulation for Water Vapor
121	! Journal of Geophysical Research, vol. 91., D8, pp 8649-8666
122	!
123	real(r8):: coefh(2,4) = reshape( &
124	(/ (/5.46557e+01,-7.30387e-02/), &
125	(/1.09311e+02,-1.46077e-01/), &
126	(/5.11479e+01,-6.82615e-02/), &
127	(/1.02296e+02,-1.36523e-01/) /), (/2,4/) )
128	!
129	real(r8):: coefj(3,2) = reshape( &
130	(/ (/2.82096e-02,2.47836e-04,1.16904e-06/), &
131	(/9.27379e-02,8.04454e-04,6.88844e-06/) /), (/3,2/) )
132	!
133	real(r8):: coefk(3,2) = reshape( &
134	(/ (/2.48852e-01,2.09667e-03,2.60377e-06/) , &
135	(/1.03594e+00,6.58620e-03,4.04456e-06/) /), (/3,2/) )
136
137	integer, parameter :: ntemp = 192 ! Number of temperatures in H2O sat. table for Tp
138	real(r8) :: estblh2o(0:ntemp) ! saturation vapor pressure for H2O for Tp rang
139	integer, parameter :: o_fa = 6 ! Degree+1 of poly of T_e for absorptivity as U->inf.
140	integer, parameter :: o_fe = 6 ! Degree+1 of poly of T_e for emissivity as U->inf.
141
142	!-----------------------------------------------------------------------------
143	! Data for f in C/H/E fit -- value of A and E as U->infinity
144	! New C/LT/E fit (Hitran 2K, CKD 2.4) -- no change
145	! These values are determined by integrals of Planck functions or
146	! derivatives of Planck functions only.
147	!-----------------------------------------------------------------------------
148	!
149	! fa/fe coefficients for 2 bands (0-800 & 1200-2200, 800-1200 cm^-1)
150	!
151	! Coefficients of polynomial for f_a in T_e
152	!
153	real(r8), parameter:: fat(o_fa,nbands) = reshape( (/ &
154	(/-1.06665373E-01, 2.90617375E-02, -2.70642049E-04, & ! 0-800&1200-2200 cm^-1
155	1.07595511E-06, -1.97419681E-09, 1.37763374E-12/), & ! 0-800&1200-2200 cm^-1
156	(/ 1.10666537E+00, -2.90617375E-02, 2.70642049E-04, & ! 800-1200 cm^-1
157	-1.07595511E-06, 1.97419681E-09, -1.37763374E-12/) /) & ! 800-1200 cm^-1
158	, (/o_fa,nbands/) )
159	!
160	! Coefficients of polynomial for f_e in T_e
161	!
162	real(r8), parameter:: fet(o_fe,nbands) = reshape( (/ &
163	(/3.46148163E-01, 1.51240299E-02, -1.21846479E-04, & ! 0-800&1200-2200 cm^-1
164	4.04970123E-07, -6.15368936E-10, 3.52415071E-13/), & ! 0-800&1200-2200 cm^-1
165	(/6.53851837E-01, -1.51240299E-02, 1.21846479E-04, & ! 800-1200 cm^-1
166	-4.04970123E-07, 6.15368936E-10, -3.52415071E-13/) /) & ! 800-1200 cm^-1
167	, (/o_fa,nbands/) )
168
169
170	real(r8) :: gravit ! Acceleration of gravity (cgs)
171	real(r8) :: rga ! 1./gravit
172	real(r8) :: gravmks ! Acceleration of gravity (mks)
173	real(r8) :: cpair ! Specific heat of dry air
174	real(r8) :: epsilo ! Ratio of mol. wght of H2O to dry air
175	real(r8) :: epsqs ! Ratio of mol. wght of H2O to dry air
176	real(r8) :: sslp ! Standard sea-level pressure
177	real(r8) :: stebol ! Stefan-Boltzmann's constant
178	real(r8) :: rgsslp ! 0.5/(gravit*sslp)
179	real(r8) :: dpfo3 ! Voigt correction factor for O3
180	real(r8) :: dpfco2 ! Voigt correction factor for CO2
181	real(r8) :: dayspy ! Number of days per 1 year
182	real(r8) :: pie ! 3.14.....
183	real(r8) :: mwdry ! molecular weight dry air ~ kg/kmole (shr_const_mwdair)
184	real(r8) :: scon ! solar constant (not used in WRF)
185	real(r8) :: co2mmr
186	real(r8) :: mwco2 ! molecular weight of carbon dioxide
187	real(r8) :: mwh2o ! molecular weight water vapor (shr_const_mwwv)
188	real(r8) :: mwch4 ! molecular weight ch4
189	real(r8) :: mwn2o ! molecular weight n2o
190	real(r8) :: mwf11 ! molecular weight cfc11
191	real(r8) :: mwf12 ! molecular weight cfc12
192	real(r8) :: cappa ! R/Cp
193	real(r8) :: rair ! Gas constant for dry air (J/K/kg)
194	real(r8) :: tmelt ! freezing T of fresh water ~ K
195	real(r8) :: r_universal ! Universal gas constant ~ J/K/kmole
196	real(r8) :: latvap ! latent heat of evaporation ~ J/kg
197	real(r8) :: latice ! latent heat of fusion ~ J/kg
198	real(r8) :: zvir ! R_V/R_D - 1.
199	integer plenest ! length of saturation vapor pressure table
200	parameter (plenest=250)
201	!
202	! Table of saturation vapor pressure values es from tmin degrees
203	! to tmax+1 degrees k in one degree increments. ttrice defines the
204	! transition region where es is a combination of ice & water values
205	!
206	real(r8) estbl(plenest) ! table values of saturation vapor pressure
207	real(r8) tmin ! min temperature (K) for table
208	real(r8) tmax ! max temperature (K) for table
209	real(r8) pcf(6) ! polynomial coeffs -> es transition water to ice
210	!real(r8), allocatable :: pin(:) ! ozone pressure level (levsiz)
211	!real(r8), allocatable :: ozmix(:,:,:) ! mixing ratio
212	!real(r8), allocatable, target :: abstot_3d(:,:,:,:) ! Non-adjacent layer absorptivites
213	!real(r8), allocatable, target :: absnxt_3d(:,:,:,:) ! Nearest layer absorptivities
214	!real(r8), allocatable, target :: emstot_3d(:,:,:) ! Total emissivity
215
216	!From aer_optics.F90 module
217	integer, parameter :: idxVIS = 8 ! index to visible band
218	integer, parameter :: nrh = 1000 ! number of relative humidity values for look-up-table
219	integer, parameter :: nspint = 19 ! number of spectral intervals
220
221	! These are now allocatable, JM 20090612
222	real(r8), allocatable, dimension(:,:) :: ksul ! (nrh, nspint) ! sulfate specific extinction ( m^2 g-1 )
223	real(r8), allocatable, dimension(:,:) :: wsul ! (nrh, nspint) ! sulfate single scattering albedo
224	real(r8), allocatable, dimension(:,:) :: gsul ! (nrh, nspint) ! sulfate asymmetry parameter
225	real(r8), allocatable, dimension(:,:) :: ksslt ! (nrh, nspint) ! sea-salt specific extinction ( m^2 g-1 )
226	real(r8), allocatable, dimension(:,:) :: wsslt ! (nrh, nspint) ! sea-salt single scattering albedo
227	real(r8), allocatable, dimension(:,:) :: gsslt ! (nrh, nspint) ! sea-salt asymmetry parameter
228	real(r8), allocatable, dimension(:,:) :: kcphil ! (nrh, nspint) ! hydrophilic carbon specific extinction ( m^2 g-1 )
229	real(r8), allocatable, dimension(:,:) :: wcphil ! (nrh, nspint) ! hydrophilic carbon single scattering albedo
230	real(r8), allocatable, dimension(:,:) :: gcphil ! (nrh, nspint) ! hydrophilic carbon asymmetry parameter
231
232	real(r8) :: kbg(nspint) ! background specific extinction ( m^2 g-1 )
233	real(r8) :: wbg(nspint) ! background single scattering albedo
234	real(r8) :: gbg(nspint) ! background asymmetry parameter
235	real(r8) :: kcphob(nspint) ! hydrophobic carbon specific extinction ( m^2 g-1 )
236	real(r8) :: wcphob(nspint) ! hydrophobic carbon single scattering albedo
237	real(r8) :: gcphob(nspint) ! hydrophobic carbon asymmetry parameter
238	real(r8) :: kcb(nspint) ! black carbon specific extinction ( m^2 g-1 )
239	real(r8) :: wcb(nspint) ! black carbon single scattering albedo
240	real(r8) :: gcb(nspint) ! black carbon asymmetry parameter
241	real(r8) :: kvolc(nspint) ! volcanic specific extinction ( m^2 g-1)
242	real(r8) :: wvolc(nspint) ! volcanic single scattering albedo
243	real(r8) :: gvolc(nspint) ! volcanic asymmetry parameter
244
245	real(r8) :: kdst(ndstsz, nspint) ! dust specific extinction ( m^2 g-1 )
246	real(r8) :: wdst(ndstsz, nspint) ! dust single scattering albedo
247	real(r8) :: gdst(ndstsz, nspint) ! dust asymmetry parameter
248	!
249	!From comozp.F90 module
250	real(r8) cplos ! constant for ozone path length integral
251	real(r8) cplol ! constant for ozone path length integral
252
253	!From ghg_surfvals.F90 module
254	real(r8) :: co2vmr = 3.550e-4 ! co2 volume mixing ratio
255	real(r8) :: n2ovmr = 0.311e-6 ! n2o volume mixing ratio
256	real(r8) :: ch4vmr = 1.714e-6 ! ch4 volume mixing ratio
257	real(r8) :: f11vmr = 0.280e-9 ! cfc11 volume mixing ratio
258	real(r8) :: f12vmr = 0.503e-9 ! cfc12 volume mixing ratio
259
260	integer, parameter :: cyr = 233 ! number of years of co2 data
261
262	integer :: yrdata(cyr) = &
263	(/ 1869, 1870, 1871, 1872, 1873, 1874, 1875, &
264	1876, 1877, 1878, 1879, 1880, 1881, 1882, &
265	1883, 1884, 1885, 1886, 1887, 1888, 1889, &
266	1890, 1891, 1892, 1893, 1894, 1895, 1896, &
267	1897, 1898, 1899, 1900, 1901, 1902, 1903, &
268	1904, 1905, 1906, 1907, 1908, 1909, 1910, &
269	1911, 1912, 1913, 1914, 1915, 1916, 1917, &
270	1918, 1919, 1920, 1921, 1922, 1923, 1924, &
271	1925, 1926, 1927, 1928, 1929, 1930, 1931, &
272	1932, 1933, 1934, 1935, 1936, 1937, 1938, &
273	1939, 1940, 1941, 1942, 1943, 1944, 1945, &
274	1946, 1947, 1948, 1949, 1950, 1951, 1952, &
275	1953, 1954, 1955, 1956, 1957, 1958, 1959, &
276	1960, 1961, 1962, 1963, 1964, 1965, 1966, &
277	1967, 1968, 1969, 1970, 1971, 1972, 1973, &
278	1974, 1975, 1976, 1977, 1978, 1979, 1980, &
279	1981, 1982, 1983, 1984, 1985, 1986, 1987, &
280	1988, 1989, 1990, 1991, 1992, 1993, 1994, &
281	1995, 1996, 1997, 1998, 1999, 2000, 2001, &
282	2002, 2003, 2004, 2005, 2006, 2007, 2008, &
283	2009, 2010, 2011, 2012, 2013, 2014, 2015, &
284	2016, 2017, 2018, 2019, 2020, 2021, 2022, &
285	2023, 2024, 2025, 2026, 2027, 2028, 2029, &
286	2030, 2031, 2032, 2033, 2034, 2035, 2036, &
287	2037, 2038, 2039, 2040, 2041, 2042, 2043, &
288	2044, 2045, 2046, 2047, 2048, 2049, 2050, &
289	2051, 2052, 2053, 2054, 2055, 2056, 2057, &
290	2058, 2059, 2060, 2061, 2062, 2063, 2064, &
291	2065, 2066, 2067, 2068, 2069, 2070, 2071, &
292	2072, 2073, 2074, 2075, 2076, 2077, 2078, &
293	2079, 2080, 2081, 2082, 2083, 2084, 2085, &
294	2086, 2087, 2088, 2089, 2090, 2091, 2092, &
295	2093, 2094, 2095, 2096, 2097, 2098, 2099, &
296	2100, 2101 /)
297
298	! A2 future scenario
299	real(r8) :: co2(cyr) = &
300	(/ 289.263, 289.263, 289.416, 289.577, 289.745, 289.919, 290.102, &
301	290.293, 290.491, 290.696, 290.909, 291.129, 291.355, 291.587, 291.824, &
302	292.066, 292.313, 292.563, 292.815, 293.071, 293.328, 293.586, 293.843, &
303	294.098, 294.35, 294.598, 294.842, 295.082, 295.32, 295.558, 295.797, &
304	296.038, 296.284, 296.535, 296.794, 297.062, 297.338, 297.62, 297.91, &
305	298.204, 298.504, 298.806, 299.111, 299.419, 299.729, 300.04, 300.352, &
306	300.666, 300.98, 301.294, 301.608, 301.923, 302.237, 302.551, 302.863, &
307	303.172, 303.478, 303.779, 304.075, 304.366, 304.651, 304.93, 305.206, &
308	305.478, 305.746, 306.013, 306.28, 306.546, 306.815, 307.087, 307.365, &
309	307.65, 307.943, 308.246, 308.56, 308.887, 309.228, 309.584, 309.956, &
310	310.344, 310.749, 311.172, 311.614, 312.077, 312.561, 313.068, 313.599, &
311	314.154, 314.737, 315.347, 315.984, 316.646, 317.328, 318.026, 318.742, &
312	319.489, 320.282, 321.133, 322.045, 323.021, 324.06, 325.155, 326.299, &
313	327.484, 328.698, 329.933, 331.194, 332.499, 333.854, 335.254, 336.69, &
314	338.15, 339.628, 341.125, 342.65, 344.206, 345.797, 347.397, 348.98, &
315	350.551, 352.1, 354.3637, 355.7772, 357.1601, 358.5306, 359.9046, &
316	361.4157, 363.0445, 364.7761, 366.6064, 368.5322, 370.534, 372.5798, &
317	374.6564, 376.7656, 378.9087, 381.0864, 383.2994, 385.548, 387.8326, &
318	390.1536, 392.523, 394.9625, 397.4806, 400.075, 402.7444, 405.4875, &
319	408.3035, 411.1918, 414.1518, 417.1831, 420.2806, 423.4355, 426.6442, &
320	429.9076, 433.2261, 436.6002, 440.0303, 443.5168, 447.06, 450.6603, &
321	454.3059, 457.9756, 461.6612, 465.3649, 469.0886, 472.8335, 476.6008, &
322	480.3916, 484.2069, 488.0473, 491.9184, 495.8295, 499.7849, 503.7843, &
323	507.8278, 511.9155, 516.0476, 520.2243, 524.4459, 528.7127, 533.0213, &
324	537.3655, 541.7429, 546.1544, 550.6005, 555.0819, 559.5991, 564.1525, &
325	568.7429, 573.3701, 578.0399, 582.7611, 587.5379, 592.3701, 597.2572, &
326	602.1997, 607.1975, 612.2507, 617.3596, 622.524, 627.7528, 633.0616, &
327	638.457, 643.9384, 649.505, 655.1568, 660.8936, 666.7153, 672.6219, &
328	678.6133, 684.6945, 690.8745, 697.1569, 703.5416, 710.0284, 716.6172, &
329	723.308, 730.1008, 736.9958, 743.993, 751.0975, 758.3183, 765.6594, &
330	773.1207, 780.702, 788.4033, 796.2249, 804.1667, 812.2289, 820.4118, &
331	828.6444, 828.6444 /)
332
333	integer :: ntoplw ! top level to solve for longwave cooling (WRF sets this to 1 for model top below 10 mb)
334
335	logical :: masterproc = .true.
336	logical :: ozncyc ! true => cycle ozone dataset
337	! logical :: dosw ! True => shortwave calculation this timestep
338	! logical :: dolw ! True => longwave calculation this timestep
339	logical :: indirect ! True => include indirect radiative effects of sulfate aerosols
340	! logical :: doabsems ! True => abs/emiss calculation this timestep
341	logical :: radforce = .false. ! True => calculate aerosol shortwave forcing
342	logical :: trace_gas=.false. ! set true for chemistry
343	logical :: strat_volcanic = .false. ! True => volcanic aerosol mass available
344
345	real(r8) retab(95)
346	!
347	! Tabulated values of re(T) in the temperature interval
348	! 180 K -- 274 K; hexagonal columns assumed:
349	!
350	data retab / &
351	5.92779, 6.26422, 6.61973, 6.99539, 7.39234, &
352	7.81177, 8.25496, 8.72323, 9.21800, 9.74075, 10.2930, &
353	10.8765, 11.4929, 12.1440, 12.8317, 13.5581, 14.2319, &
354	15.0351, 15.8799, 16.7674, 17.6986, 18.6744, 19.6955, &
355	20.7623, 21.8757, 23.0364, 24.2452, 25.5034, 26.8125, &
356	27.7895, 28.6450, 29.4167, 30.1088, 30.7306, 31.2943, &
357	31.8151, 32.3077, 32.7870, 33.2657, 33.7540, 34.2601, &
358	34.7892, 35.3442, 35.9255, 36.5316, 37.1602, 37.8078, &
359	38.4720, 39.1508, 39.8442, 40.5552, 41.2912, 42.0635, &
360	42.8876, 43.7863, 44.7853, 45.9170, 47.2165, 48.7221, &
361	50.4710, 52.4980, 54.8315, 57.4898, 60.4785, 63.7898, &
362	65.5604, 71.2885, 75.4113, 79.7368, 84.2351, 88.8833, &
363	93.6658, 98.5739, 103.603, 108.752, 114.025, 119.424, &
364	124.954, 130.630, 136.457, 142.446, 148.608, 154.956, &
365	161.503, 168.262, 175.248, 182.473, 189.952, 197.699, &
366	205.728, 214.055, 222.694, 231.661, 240.971, 250.639/
367	!
368	save retab
369	contains
370
371
372
373	subroutine sortarray(n, ain, indxa)
374	!-----------------------------------------------
375	!
376	! Purpose:
377	! Sort an array
378	! Alogrithm:
379	! Based on Shell's sorting method.
380	!
381	! Author: T. Craig
382	!-----------------------------------------------
383	! use shr_kind_mod, only: r8 => shr_kind_r8
384	implicit none
385	!
386	! Arguments
387	!
388	integer , intent(in) :: n ! total number of elements
389	integer , intent(inout) :: indxa(n) ! array of integers
390	real(r8), intent(inout) :: ain(n) ! array to sort
391	!
392	! local variables
393	!
394	integer :: i, j ! Loop indices
395	integer :: ni ! Starting increment
396	integer :: itmp ! Temporary index
397	real(r8):: atmp ! Temporary value to swap
398
399	ni = 1
400	do while(.TRUE.)
401	ni = 3*ni + 1
402	if (ni <= n) cycle
403	exit
404	end do
405
406	do while(.TRUE.)
407	ni = ni/3
408	do i = ni + 1, n
409	atmp = ain(i)
410	itmp = indxa(i)
411	j = i
412	do while(.TRUE.)
413	if (ain(j-ni) <= atmp) exit
414	ain(j) = ain(j-ni)
415	indxa(j) = indxa(j-ni)
416	j = j - ni
417	if (j > ni) cycle
418	exit
419	end do
420	ain(j) = atmp
421	indxa(j) = itmp
422	end do
423	if (ni > 1) cycle
424	exit
425	end do
426	return
427
428	end subroutine sortarray
429	subroutine trcab(lchnk ,ncol ,pcols, pverp, &
430	k1 ,k2 ,ucfc11 ,ucfc12 ,un2o0 , &
431	un2o1 ,uch4 ,uco211 ,uco212 ,uco213 , &
432	uco221 ,uco222 ,uco223 ,bn2o0 ,bn2o1 , &
433	bch4 ,to3co2 ,pnm ,dw ,pnew , &
434	s2c ,uptype ,dplh2o ,abplnk1 ,tco2 , &
435	th2o ,to3 ,abstrc , &
436	aer_trn_ttl)
437	!-----------------------------------------------------------------------
438	!
439	! Purpose:
440	! Calculate absorptivity for non nearest layers for CH4, N2O, CFC11 and
441	! CFC12.
442	!
443	! Method:
444	! See CCM3 description for equations.
445	!
446	! Author: J. Kiehl
447	!
448	!-----------------------------------------------------------------------
449	! use shr_kind_mod, only: r8 => shr_kind_r8
450	! use ppgrid
451	! use volcrad
452
453	implicit none
454
455	!------------------------------Arguments--------------------------------
456	!
457	! Input arguments
458	!
459	integer, intent(in) :: lchnk ! chunk identifier
460	integer, intent(in) :: ncol ! number of atmospheric columns
461	integer, intent(in) :: pcols, pverp
462	integer, intent(in) :: k1,k2 ! level indices
463	!
464	real(r8), intent(in) :: to3co2(pcols) ! pressure weighted temperature
465	real(r8), intent(in) :: pnm(pcols,pverp) ! interface pressures
466	real(r8), intent(in) :: ucfc11(pcols,pverp) ! CFC11 path length
467	real(r8), intent(in) :: ucfc12(pcols,pverp) ! CFC12 path length
468	real(r8), intent(in) :: un2o0(pcols,pverp) ! N2O path length
469	!
470	real(r8), intent(in) :: un2o1(pcols,pverp) ! N2O path length (hot band)
471	real(r8), intent(in) :: uch4(pcols,pverp) ! CH4 path length
472	real(r8), intent(in) :: uco211(pcols,pverp) ! CO2 9.4 micron band path length
473	real(r8), intent(in) :: uco212(pcols,pverp) ! CO2 9.4 micron band path length
474	real(r8), intent(in) :: uco213(pcols,pverp) ! CO2 9.4 micron band path length
475	!
476	real(r8), intent(in) :: uco221(pcols,pverp) ! CO2 10.4 micron band path length
477	real(r8), intent(in) :: uco222(pcols,pverp) ! CO2 10.4 micron band path length
478	real(r8), intent(in) :: uco223(pcols,pverp) ! CO2 10.4 micron band path length
479	real(r8), intent(in) :: bn2o0(pcols,pverp) ! pressure factor for n2o
480	real(r8), intent(in) :: bn2o1(pcols,pverp) ! pressure factor for n2o
481	!
482	real(r8), intent(in) :: bch4(pcols,pverp) ! pressure factor for ch4
483	real(r8), intent(in) :: dw(pcols) ! h2o path length
484	real(r8), intent(in) :: pnew(pcols) ! pressure
485	real(r8), intent(in) :: s2c(pcols,pverp) ! continuum path length
486	real(r8), intent(in) :: uptype(pcols,pverp) ! p-type h2o path length
487	!
488	real(r8), intent(in) :: dplh2o(pcols) ! p squared h2o path length
489	real(r8), intent(in) :: abplnk1(14,pcols,pverp) ! Planck factor
490	real(r8), intent(in) :: tco2(pcols) ! co2 transmission factor
491	real(r8), intent(in) :: th2o(pcols) ! h2o transmission factor
492	real(r8), intent(in) :: to3(pcols) ! o3 transmission factor
493
494	real(r8), intent(in) :: aer_trn_ttl(pcols,pverp,pverp,bnd_nbr_LW) ! aer trn.
495
496	!
497	! Output Arguments
498	!
499	real(r8), intent(out) :: abstrc(pcols) ! total trace gas absorptivity
500	!
501	!--------------------------Local Variables------------------------------
502	!
503	integer i,l ! loop counters
504
505	real(r8) sqti(pcols) ! square root of mean temp
506	real(r8) du1 ! cfc11 path length
507	real(r8) du2 ! cfc12 path length
508	real(r8) acfc1 ! cfc11 absorptivity 798 cm-1
509	real(r8) acfc2 ! cfc11 absorptivity 846 cm-1
510	!
511	real(r8) acfc3 ! cfc11 absorptivity 933 cm-1
512	real(r8) acfc4 ! cfc11 absorptivity 1085 cm-1
513	real(r8) acfc5 ! cfc12 absorptivity 889 cm-1
514	real(r8) acfc6 ! cfc12 absorptivity 923 cm-1
515	real(r8) acfc7 ! cfc12 absorptivity 1102 cm-1
516	!
517	real(r8) acfc8 ! cfc12 absorptivity 1161 cm-1
518	real(r8) du01 ! n2o path length
519	real(r8) dbeta01 ! n2o pressure factor
520	real(r8) dbeta11 ! "
521	real(r8) an2o1 ! absorptivity of 1285 cm-1 n2o band
522	!
523	real(r8) du02 ! n2o path length
524	real(r8) dbeta02 ! n2o pressure factor
525	real(r8) an2o2 ! absorptivity of 589 cm-1 n2o band
526	real(r8) du03 ! n2o path length
527	real(r8) dbeta03 ! n2o pressure factor
528	!
529	real(r8) an2o3 ! absorptivity of 1168 cm-1 n2o band
530	real(r8) duch4 ! ch4 path length
531	real(r8) dbetac ! ch4 pressure factor
532	real(r8) ach4 ! absorptivity of 1306 cm-1 ch4 band
533	real(r8) du11 ! co2 path length
534	!
535	real(r8) du12 ! "
536	real(r8) du13 ! "
537	real(r8) dbetc1 ! co2 pressure factor
538	real(r8) dbetc2 ! co2 pressure factor
539	real(r8) aco21 ! absorptivity of 1064 cm-1 band
540	!
541	real(r8) du21 ! co2 path length
542	real(r8) du22 ! "
543	real(r8) du23 ! "
544	real(r8) aco22 ! absorptivity of 961 cm-1 band
545	real(r8) tt(pcols) ! temp. factor for h2o overlap factor
546	!
547	real(r8) psi1 ! "
548	real(r8) phi1 ! "
549	real(r8) p1 ! h2o overlap factor
550	real(r8) w1 ! "
551	real(r8) ds2c(pcols) ! continuum path length
552	!
553	real(r8) duptyp(pcols) ! p-type path length
554	real(r8) tw(pcols,6) ! h2o transmission factor
555	real(r8) g1(6) ! "
556	real(r8) g2(6) ! "
557	real(r8) g3(6) ! "
558	!
559	real(r8) g4(6) ! "
560	real(r8) ab(6) ! h2o temp. factor
561	real(r8) bb(6) ! "
562	real(r8) abp(6) ! "
563	real(r8) bbp(6) ! "
564	!
565	real(r8) tcfc3 ! transmission for cfc11 band
566	real(r8) tcfc4 ! transmission for cfc11 band
567	real(r8) tcfc6 ! transmission for cfc12 band
568	real(r8) tcfc7 ! transmission for cfc12 band
569	real(r8) tcfc8 ! transmission for cfc12 band
570	!
571	real(r8) tlw ! h2o transmission
572	real(r8) tch4 ! ch4 transmission
573	!
574	!--------------------------Data Statements------------------------------
575	!
576	data g1 /0.0468556,0.0397454,0.0407664,0.0304380,0.0540398,0.0321962/
577	data g2 /14.4832,4.30242,5.23523,3.25342,0.698935,16.5599/
578	data g3 /26.1898,18.4476,15.3633,12.1927,9.14992,8.07092/
579	data g4 /0.0261782,0.0369516,0.0307266,0.0243854,0.0182932,0.0161418/
580	data ab /3.0857e-2,2.3524e-2,1.7310e-2,2.6661e-2,2.8074e-2,2.2915e-2/
581	data bb /-1.3512e-4,-6.8320e-5,-3.2609e-5,-1.0228e-5,-9.5743e-5,-1.0304e-4/
582	data abp/2.9129e-2,2.4101e-2,1.9821e-2,2.6904e-2,2.9458e-2,1.9892e-2/
583	data bbp/-1.3139e-4,-5.5688e-5,-4.6380e-5,-8.0362e-5,-1.0115e-4,-8.8061e-5/
584	!
585	!--------------------------Statement Functions--------------------------
586	!
587	real(r8) func, u, b
588	func(u,b) = u/sqrt(4.0 + u*(1.0 + 1.0 / b))
589	!
590	!------------------------------------------------------------------------
591	!
592	do i = 1,ncol
593	sqti(i) = sqrt(to3co2(i))
594	!
595	! h2o transmission
596	!
597	tt(i) = abs(to3co2(i) - 250.0)
598	ds2c(i) = abs(s2c(i,k1) - s2c(i,k2))
599	duptyp(i) = abs(uptype(i,k1) - uptype(i,k2))
600	end do
601	!
602	do l = 1,6
603	do i = 1,ncol
604	psi1 = exp(abp(l)tt(i) + bbp(l)tt(i)*tt(i))
605	phi1 = exp(ab(l)tt(i) + bb(l)tt(i)*tt(i))
606	p1 = pnew(i)*(psi1/phi1)/sslp
607	w1 = dw(i)*phi1
608	tw(i,l) = exp(-g1(l)p1(sqrt(1.0 + g2(l)*(w1/p1)) - 1.0) - &
609	g3(l)ds2c(i)-g4(l)duptyp(i))
610	end do
611	end do
612	!
613	do i=1,ncol
614	tw(i,1)=tw(i,1)(0.7aer_trn_ttl(i,k1,k2,idx_LW_0650_0800)+&! l=1: 0750--0820 cm-1
615	0.3*aer_trn_ttl(i,k1,k2,idx_LW_0800_1000))
616	tw(i,2)=tw(i,2)*aer_trn_ttl(i,k1,k2,idx_LW_0800_1000) ! l=2: 0820--0880 cm-1
617	tw(i,3)=tw(i,3)*aer_trn_ttl(i,k1,k2,idx_LW_0800_1000) ! l=3: 0880--0900 cm-1
618	tw(i,4)=tw(i,4)*aer_trn_ttl(i,k1,k2,idx_LW_0800_1000) ! l=4: 0900--1000 cm-1
619	tw(i,5)=tw(i,5)*aer_trn_ttl(i,k1,k2,idx_LW_1000_1200) ! l=5: 1000--1120 cm-1
620	tw(i,6)=tw(i,6)*aer_trn_ttl(i,k1,k2,idx_LW_1000_1200) ! l=6: 1120--1170 cm-1
621	end do ! end loop over lon
622	do i = 1,ncol
623	du1 = abs(ucfc11(i,k1) - ucfc11(i,k2))
624	du2 = abs(ucfc12(i,k1) - ucfc12(i,k2))
625	!
626	! cfc transmissions
627	!
628	tcfc3 = exp(-175.005*du1)
629	tcfc4 = exp(-1202.18*du1)
630	tcfc6 = exp(-5786.73*du2)
631	tcfc7 = exp(-2873.51*du2)
632	tcfc8 = exp(-2085.59*du2)
633	!
634	! Absorptivity for CFC11 bands
635	!
636	acfc1 = 50.0(1.0 - exp(-54.09du1))tw(i,1)abplnk1(7,i,k2)
637	acfc2 = 60.0(1.0 - exp(-5130.03du1))tw(i,2)abplnk1(8,i,k2)
638	acfc3 = 60.0(1.0 - tcfc3)tw(i,4)tcfc6abplnk1(9,i,k2)
639	acfc4 = 100.0(1.0 - tcfc4)tw(i,5)*abplnk1(10,i,k2)
640	!
641	! Absorptivity for CFC12 bands
642	!
643	acfc5 = 45.0(1.0 - exp(-1272.35du2))tw(i,3)abplnk1(11,i,k2)
644	acfc6 = 50.0(1.0 - tcfc6) tw(i,4) * abplnk1(12,i,k2)
645	acfc7 = 80.0(1.0 - tcfc7) tw(i,5) * tcfc4*abplnk1(13,i,k2)
646	acfc8 = 70.0(1.0 - tcfc8) tw(i,6) * abplnk1(14,i,k2)
647	!
648	! Emissivity for CH4 band 1306 cm-1
649	!
650	tlw = exp(-1.0*sqrt(dplh2o(i)))
651	tlw=tlw*aer_trn_ttl(i,k1,k2,idx_LW_1200_2000)
652	duch4 = abs(uch4(i,k1) - uch4(i,k2))
653	dbetac = abs(bch4(i,k1) - bch4(i,k2))/duch4
654	ach4 = 6.00444sqti(i)log(1.0 + func(duch4,dbetac))tlwabplnk1(3,i,k2)
655	tch4 = 1.0/(1.0 + 0.02*func(duch4,dbetac))
656	!
657	! Absorptivity for N2O bands
658	!
659	du01 = abs(un2o0(i,k1) - un2o0(i,k2))
660	du11 = abs(un2o1(i,k1) - un2o1(i,k2))
661	dbeta01 = abs(bn2o0(i,k1) - bn2o0(i,k2))/du01
662	dbeta11 = abs(bn2o1(i,k1) - bn2o1(i,k2))/du11
663	!
664	! 1285 cm-1 band
665	!
666	an2o1 = 2.35558sqti(i)log(1.0 + func(du01,dbeta01) &
667	+ func(du11,dbeta11))tlwtch4*abplnk1(4,i,k2)
668	du02 = 0.100090*du01
669	du12 = 0.0992746*du11
670	dbeta02 = 0.964282*dbeta01
671	!
672	! 589 cm-1 band
673	!
674	an2o2 = 2.65581sqti(i)log(1.0 + func(du02,dbeta02) + &
675	func(du12,dbeta02))th2o(i)tco2(i)*abplnk1(5,i,k2)
676	du03 = 0.0333767*du01
677	dbeta03 = 0.982143*dbeta01
678	!
679	! 1168 cm-1 band
680	!
681	an2o3 = 2.54034sqti(i)log(1.0 + func(du03,dbeta03))* &
682	tw(i,6)tcfc8abplnk1(6,i,k2)
683	!
684	! Emissivity for 1064 cm-1 band of CO2
685	!
686	du11 = abs(uco211(i,k1) - uco211(i,k2))
687	du12 = abs(uco212(i,k1) - uco212(i,k2))
688	du13 = abs(uco213(i,k1) - uco213(i,k2))
689	dbetc1 = 2.97558abs(pnm(i,k1) + pnm(i,k2))/(2.0sslp*sqti(i))
690	dbetc2 = 2.0*dbetc1
691	aco21 = 3.7571sqti(i)log(1.0 + func(du11,dbetc1) &
692	+ func(du12,dbetc2) + func(du13,dbetc2)) &
693	to3(i)tw(i,5)tcfc4tcfc7*abplnk1(2,i,k2)
694	!
695	! Emissivity for 961 cm-1 band
696	!
697	du21 = abs(uco221(i,k1) - uco221(i,k2))
698	du22 = abs(uco222(i,k1) - uco222(i,k2))
699	du23 = abs(uco223(i,k1) - uco223(i,k2))
700	aco22 = 3.8443sqti(i)log(1.0 + func(du21,dbetc1) &
701	+ func(du22,dbetc1) + func(du23,dbetc2)) &
702	tw(i,4)tcfc3tcfc6abplnk1(1,i,k2)
703	!
704	! total trace gas absorptivity
705	!
706	abstrc(i) = acfc1 + acfc2 + acfc3 + acfc4 + acfc5 + acfc6 + &
707	acfc7 + acfc8 + an2o1 + an2o2 + an2o3 + ach4 + &
708	aco21 + aco22
709	end do
710	!
711	return
712	!
713	end subroutine trcab
714
715
716
717	subroutine trcabn(lchnk ,ncol ,pcols, pverp, &
718	k2 ,kn ,ucfc11 ,ucfc12 ,un2o0 , &
719	un2o1 ,uch4 ,uco211 ,uco212 ,uco213 , &
720	uco221 ,uco222 ,uco223 ,tbar ,bplnk , &
721	winpl ,pinpl ,tco2 ,th2o ,to3 , &
722	uptype ,dw ,s2c ,up2 ,pnew , &
723	abstrc ,uinpl , &
724	aer_trn_ngh)
725	!-----------------------------------------------------------------------
726	!
727	! Purpose:
728	! Calculate nearest layer absorptivity due to CH4, N2O, CFC11 and CFC12
729	!
730	! Method:
731	! Equations in CCM3 description
732	!
733	! Author: J. Kiehl
734	!
735	!-----------------------------------------------------------------------
736	!
737	! use shr_kind_mod, only: r8 => shr_kind_r8
738	! use ppgrid
739	! use volcrad
740
741	implicit none
742
743	!------------------------------Arguments--------------------------------
744	!
745	! Input arguments
746	!
747	integer, intent(in) :: lchnk ! chunk identifier
748	integer, intent(in) :: ncol ! number of atmospheric columns
749	integer, intent(in) :: pcols, pverp
750	integer, intent(in) :: k2 ! level index
751	integer, intent(in) :: kn ! level index
752	!
753	real(r8), intent(in) :: tbar(pcols,4) ! pressure weighted temperature
754	real(r8), intent(in) :: ucfc11(pcols,pverp) ! CFC11 path length
755	real(r8), intent(in) :: ucfc12(pcols,pverp) ! CFC12 path length
756	real(r8), intent(in) :: un2o0(pcols,pverp) ! N2O path length
757	real(r8), intent(in) :: un2o1(pcols,pverp) ! N2O path length (hot band)
758	!
759	real(r8), intent(in) :: uch4(pcols,pverp) ! CH4 path length
760	real(r8), intent(in) :: uco211(pcols,pverp) ! CO2 9.4 micron band path length
761	real(r8), intent(in) :: uco212(pcols,pverp) ! CO2 9.4 micron band path length
762	real(r8), intent(in) :: uco213(pcols,pverp) ! CO2 9.4 micron band path length
763	real(r8), intent(in) :: uco221(pcols,pverp) ! CO2 10.4 micron band path length
764	!
765	real(r8), intent(in) :: uco222(pcols,pverp) ! CO2 10.4 micron band path length
766	real(r8), intent(in) :: uco223(pcols,pverp) ! CO2 10.4 micron band path length
767	real(r8), intent(in) :: bplnk(14,pcols,4) ! weighted Planck fnc. for absorptivity
768	real(r8), intent(in) :: winpl(pcols,4) ! fractional path length
769	real(r8), intent(in) :: pinpl(pcols,4) ! pressure factor for subdivided layer
770	!
771	real(r8), intent(in) :: tco2(pcols) ! co2 transmission
772	real(r8), intent(in) :: th2o(pcols) ! h2o transmission
773	real(r8), intent(in) :: to3(pcols) ! o3 transmission
774	real(r8), intent(in) :: dw(pcols) ! h2o path length
775	real(r8), intent(in) :: pnew(pcols) ! pressure factor
776	!
777	real(r8), intent(in) :: s2c(pcols,pverp) ! h2o continuum factor
778	real(r8), intent(in) :: uptype(pcols,pverp) ! p-type path length
779	real(r8), intent(in) :: up2(pcols) ! p squared path length
780	real(r8), intent(in) :: uinpl(pcols,4) ! Nearest layer subdivision factor
781	real(r8), intent(in) :: aer_trn_ngh(pcols,bnd_nbr_LW)
782	! [fraction] Total transmission between
783	! nearest neighbor sub-levels
784	!
785	! Output Arguments
786	!
787	real(r8), intent(out) :: abstrc(pcols) ! total trace gas absorptivity
788
789	!
790	!--------------------------Local Variables------------------------------
791	!
792	integer i,l ! loop counters
793	!
794	real(r8) sqti(pcols) ! square root of mean temp
795	real(r8) rsqti(pcols) ! reciprocal of sqti
796	real(r8) du1 ! cfc11 path length
797	real(r8) du2 ! cfc12 path length
798	real(r8) acfc1 ! absorptivity of cfc11 798 cm-1 band
799	!
800	real(r8) acfc2 ! absorptivity of cfc11 846 cm-1 band
801	real(r8) acfc3 ! absorptivity of cfc11 933 cm-1 band
802	real(r8) acfc4 ! absorptivity of cfc11 1085 cm-1 band
803	real(r8) acfc5 ! absorptivity of cfc11 889 cm-1 band
804	real(r8) acfc6 ! absorptivity of cfc11 923 cm-1 band
805	!
806	real(r8) acfc7 ! absorptivity of cfc11 1102 cm-1 band
807	real(r8) acfc8 ! absorptivity of cfc11 1161 cm-1 band
808	real(r8) du01 ! n2o path length
809	real(r8) dbeta01 ! n2o pressure factors
810	real(r8) dbeta11 ! "
811	!
812	real(r8) an2o1 ! absorptivity of the 1285 cm-1 n2o band
813	real(r8) du02 ! n2o path length
814	real(r8) dbeta02 ! n2o pressure factor
815	real(r8) an2o2 ! absorptivity of the 589 cm-1 n2o band
816	real(r8) du03 ! n2o path length
817	!
818	real(r8) dbeta03 ! n2o pressure factor
819	real(r8) an2o3 ! absorptivity of the 1168 cm-1 n2o band
820	real(r8) duch4 ! ch4 path length
821	real(r8) dbetac ! ch4 pressure factor
822	real(r8) ach4 ! absorptivity of the 1306 cm-1 ch4 band
823	!
824	real(r8) du11 ! co2 path length
825	real(r8) du12 ! "
826	real(r8) du13 ! "
827	real(r8) dbetc1 ! co2 pressure factor
828	real(r8) dbetc2 ! co2 pressure factor
829	!
830	real(r8) aco21 ! absorptivity of the 1064 cm-1 co2 band
831	real(r8) du21 ! co2 path length
832	real(r8) du22 ! "
833	real(r8) du23 ! "
834	real(r8) aco22 ! absorptivity of the 961 cm-1 co2 band
835	!
836	real(r8) tt(pcols) ! temp. factor for h2o overlap
837	real(r8) psi1 ! "
838	real(r8) phi1 ! "
839	real(r8) p1 ! factor for h2o overlap
840	real(r8) w1 ! "
841	!
842	real(r8) ds2c(pcols) ! continuum path length
843	real(r8) duptyp(pcols) ! p-type path length
844	real(r8) tw(pcols,6) ! h2o transmission overlap
845	real(r8) g1(6) ! h2o overlap factor
846	real(r8) g2(6) ! "
847	!
848	real(r8) g3(6) ! "
849	real(r8) g4(6) ! "
850	real(r8) ab(6) ! h2o temp. factor
851	real(r8) bb(6) ! "
852	real(r8) abp(6) ! "
853	!
854	real(r8) bbp(6) ! "
855	real(r8) tcfc3 ! transmission of cfc11 band
856	real(r8) tcfc4 ! transmission of cfc11 band
857	real(r8) tcfc6 ! transmission of cfc12 band
858	real(r8) tcfc7 ! "
859	!
860	real(r8) tcfc8 ! "
861	real(r8) tlw ! h2o transmission
862	real(r8) tch4 ! ch4 transmission
863	!
864	!--------------------------Data Statements------------------------------
865	!
866	data g1 /0.0468556,0.0397454,0.0407664,0.0304380,0.0540398,0.0321962/
867	data g2 /14.4832,4.30242,5.23523,3.25342,0.698935,16.5599/
868	data g3 /26.1898,18.4476,15.3633,12.1927,9.14992,8.07092/
869	data g4 /0.0261782,0.0369516,0.0307266,0.0243854,0.0182932,0.0161418/
870	data ab /3.0857e-2,2.3524e-2,1.7310e-2,2.6661e-2,2.8074e-2,2.2915e-2/
871	data bb /-1.3512e-4,-6.8320e-5,-3.2609e-5,-1.0228e-5,-9.5743e-5,-1.0304e-4/
872	data abp/2.9129e-2,2.4101e-2,1.9821e-2,2.6904e-2,2.9458e-2,1.9892e-2/
873	data bbp/-1.3139e-4,-5.5688e-5,-4.6380e-5,-8.0362e-5,-1.0115e-4,-8.8061e-5/
874	!
875	!--------------------------Statement Functions--------------------------
876	!
877	real(r8) func, u, b
878	func(u,b) = u/sqrt(4.0 + u*(1.0 + 1.0 / b))
879	!
880	!------------------------------------------------------------------
881	!
882	do i = 1,ncol
883	sqti(i) = sqrt(tbar(i,kn))
884	rsqti(i) = 1. / sqti(i)
885	!
886	! h2o transmission
887	!
888	tt(i) = abs(tbar(i,kn) - 250.0)
889	ds2c(i) = abs(s2c(i,k2+1) - s2c(i,k2))*uinpl(i,kn)
890	duptyp(i) = abs(uptype(i,k2+1) - uptype(i,k2))*uinpl(i,kn)
891	end do
892	!
893	do l = 1,6
894	do i = 1,ncol
895	psi1 = exp(abp(l)tt(i)+bbp(l)tt(i)*tt(i))
896	phi1 = exp(ab(l)tt(i)+bb(l)tt(i)*tt(i))
897	p1 = pnew(i) * (psi1/phi1) / sslp
898	w1 = dw(i) * winpl(i,kn) * phi1
899	tw(i,l) = exp(- g1(l)p1(sqrt(1.0+g2(l)*(w1/p1))-1.0) &
900	- g3(l)ds2c(i)-g4(l)duptyp(i))
901	end do
902	end do
903	!
904	do i=1,ncol
905	tw(i,1)=tw(i,1)(0.7aer_trn_ngh(i,idx_LW_0650_0800)+&! l=1: 0750--0820 cm-1
906	0.3*aer_trn_ngh(i,idx_LW_0800_1000))
907	tw(i,2)=tw(i,2)*aer_trn_ngh(i,idx_LW_0800_1000) ! l=2: 0820--0880 cm-1
908	tw(i,3)=tw(i,3)*aer_trn_ngh(i,idx_LW_0800_1000) ! l=3: 0880--0900 cm-1
909	tw(i,4)=tw(i,4)*aer_trn_ngh(i,idx_LW_0800_1000) ! l=4: 0900--1000 cm-1
910	tw(i,5)=tw(i,5)*aer_trn_ngh(i,idx_LW_1000_1200) ! l=5: 1000--1120 cm-1
911	tw(i,6)=tw(i,6)*aer_trn_ngh(i,idx_LW_1000_1200) ! l=6: 1120--1170 cm-1
912	end do ! end loop over lon
913
914	do i = 1,ncol
915	!
916	du1 = abs(ucfc11(i,k2+1) - ucfc11(i,k2)) * winpl(i,kn)
917	du2 = abs(ucfc12(i,k2+1) - ucfc12(i,k2)) * winpl(i,kn)
918	!
919	! cfc transmissions
920	!
921	tcfc3 = exp(-175.005*du1)
922	tcfc4 = exp(-1202.18*du1)
923	tcfc6 = exp(-5786.73*du2)
924	tcfc7 = exp(-2873.51*du2)
925	tcfc8 = exp(-2085.59*du2)
926	!
927	! Absorptivity for CFC11 bands
928	!
929	acfc1 = 50.0(1.0 - exp(-54.09du1)) * tw(i,1)*bplnk(7,i,kn)
930	acfc2 = 60.0(1.0 - exp(-5130.03du1))tw(i,2)bplnk(8,i,kn)
931	acfc3 = 60.0(1.0 - tcfc3)tw(i,4)tcfc6 bplnk(9,i,kn)
932	acfc4 = 100.0(1.0 - tcfc4) tw(i,5) * bplnk(10,i,kn)
933	!
934	! Absorptivity for CFC12 bands
935	!
936	acfc5 = 45.0(1.0 - exp(-1272.35du2))tw(i,3)bplnk(11,i,kn)
937	acfc6 = 50.0(1.0 - tcfc6)tw(i,4)*bplnk(12,i,kn)
938	acfc7 = 80.0(1.0 - tcfc7) tw(i,5)tcfc4 bplnk(13,i,kn)
939	acfc8 = 70.0(1.0 - tcfc8)tw(i,6)*bplnk(14,i,kn)
940	!
941	! Absorptivity for CH4 band 1306 cm-1
942	!
943	tlw = exp(-1.0*sqrt(up2(i)))
944	tlw=tlw*aer_trn_ngh(i,idx_LW_1200_2000)
945	duch4 = abs(uch4(i,k2+1) - uch4(i,k2)) * winpl(i,kn)
946	dbetac = 2.94449 * pinpl(i,kn) * rsqti(i) / sslp
947	ach4 = 6.00444sqti(i)log(1.0 + func(duch4,dbetac)) * tlw * bplnk(3,i,kn)
948	tch4 = 1.0/(1.0 + 0.02*func(duch4,dbetac))
949	!
950	! Absorptivity for N2O bands
951	!
952	du01 = abs(un2o0(i,k2+1) - un2o0(i,k2)) * winpl(i,kn)
953	du11 = abs(un2o1(i,k2+1) - un2o1(i,k2)) * winpl(i,kn)
954	dbeta01 = 19.399 * pinpl(i,kn) * rsqti(i) / sslp
955	dbeta11 = dbeta01
956	!
957	! 1285 cm-1 band
958	!
959	an2o1 = 2.35558sqti(i)log(1.0 + func(du01,dbeta01) &
960	+ func(du11,dbeta11)) * tlw * tch4 * bplnk(4,i,kn)
961	du02 = 0.100090*du01
962	du12 = 0.0992746*du11
963	dbeta02 = 0.964282*dbeta01
964	!
965	! 589 cm-1 band
966	!
967	an2o2 = 2.65581sqti(i)log(1.0 + func(du02,dbeta02) &
968	+ func(du12,dbeta02)) * tco2(i) * th2o(i) * bplnk(5,i,kn)
969	du03 = 0.0333767*du01
970	dbeta03 = 0.982143*dbeta01
971	!
972	! 1168 cm-1 band
973	!
974	an2o3 = 2.54034sqti(i)log(1.0 + func(du03,dbeta03)) * &
975	tw(i,6) * tcfc8 * bplnk(6,i,kn)
976	!
977	! Absorptivity for 1064 cm-1 band of CO2
978	!
979	du11 = abs(uco211(i,k2+1) - uco211(i,k2)) * winpl(i,kn)
980	du12 = abs(uco212(i,k2+1) - uco212(i,k2)) * winpl(i,kn)
981	du13 = abs(uco213(i,k2+1) - uco213(i,k2)) * winpl(i,kn)
982	dbetc1 = 2.97558 * pinpl(i,kn) * rsqti(i) / sslp
983	dbetc2 = 2.0 * dbetc1
984	aco21 = 3.7571sqti(i)log(1.0 + func(du11,dbetc1) &
985	+ func(du12,dbetc2) + func(du13,dbetc2)) &
986	* to3(i) * tw(i,5) * tcfc4 * tcfc7 * bplnk(2,i,kn)
987	!
988	! Absorptivity for 961 cm-1 band of co2
989	!
990	du21 = abs(uco221(i,k2+1) - uco221(i,k2)) * winpl(i,kn)
991	du22 = abs(uco222(i,k2+1) - uco222(i,k2)) * winpl(i,kn)
992	du23 = abs(uco223(i,k2+1) - uco223(i,k2)) * winpl(i,kn)
993	aco22 = 3.8443sqti(i)log(1.0 + func(du21,dbetc1) &
994	+ func(du22,dbetc1) + func(du23,dbetc2)) &
995	* tw(i,4) * tcfc3 * tcfc6 * bplnk(1,i,kn)
996	!
997	! total trace gas absorptivity
998	!
999	abstrc(i) = acfc1 + acfc2 + acfc3 + acfc4 + acfc5 + acfc6 + &
1000	acfc7 + acfc8 + an2o1 + an2o2 + an2o3 + ach4 + &
1001	aco21 + aco22
1002	end do
1003	!
1004	return
1005	!
1006	end subroutine trcabn
1007
1008
1009
1010	subroutine trcems(lchnk ,ncol ,pcols, pverp, &
1011	k ,co2t ,pnm ,ucfc11 ,ucfc12 , &
1012	un2o0 ,un2o1 ,bn2o0 ,bn2o1 ,uch4 , &
1013	bch4 ,uco211 ,uco212 ,uco213 ,uco221 , &
1014	uco222 ,uco223 ,uptype ,w ,s2c , &
1015	up2 ,emplnk ,th2o ,tco2 ,to3 , &
1016	emstrc , &
1017	aer_trn_ttl)
1018	!-----------------------------------------------------------------------
1019	!
1020	! Purpose:
1021	! Calculate emissivity for CH4, N2O, CFC11 and CFC12 bands.
1022	!
1023	! Method:
1024	! See CCM3 Description for equations.
1025	!
1026	! Author: J. Kiehl
1027	!
1028	!-----------------------------------------------------------------------
1029	! use shr_kind_mod, only: r8 => shr_kind_r8
1030	! use ppgrid
1031	! use volcrad
1032
1033	implicit none
1034
1035	!
1036	!------------------------------Arguments--------------------------------
1037	!
1038	! Input arguments
1039	!
1040	integer, intent(in) :: lchnk ! chunk identifier
1041	integer, intent(in) :: ncol ! number of atmospheric columns
1042	integer, intent(in) :: pcols, pverp
1043
1044	real(r8), intent(in) :: co2t(pcols,pverp) ! pressure weighted temperature
1045	real(r8), intent(in) :: pnm(pcols,pverp) ! interface pressure
1046	real(r8), intent(in) :: ucfc11(pcols,pverp) ! CFC11 path length
1047	real(r8), intent(in) :: ucfc12(pcols,pverp) ! CFC12 path length
1048	real(r8), intent(in) :: un2o0(pcols,pverp) ! N2O path length
1049	!
1050	real(r8), intent(in) :: un2o1(pcols,pverp) ! N2O path length (hot band)
1051	real(r8), intent(in) :: uch4(pcols,pverp) ! CH4 path length
1052	real(r8), intent(in) :: uco211(pcols,pverp) ! CO2 9.4 micron band path length
1053	real(r8), intent(in) :: uco212(pcols,pverp) ! CO2 9.4 micron band path length
1054	real(r8), intent(in) :: uco213(pcols,pverp) ! CO2 9.4 micron band path length
1055	!
1056	real(r8), intent(in) :: uco221(pcols,pverp) ! CO2 10.4 micron band path length
1057	real(r8), intent(in) :: uco222(pcols,pverp) ! CO2 10.4 micron band path length
1058	real(r8), intent(in) :: uco223(pcols,pverp) ! CO2 10.4 micron band path length
1059	real(r8), intent(in) :: uptype(pcols,pverp) ! continuum path length
1060	real(r8), intent(in) :: bn2o0(pcols,pverp) ! pressure factor for n2o
1061	!
1062	real(r8), intent(in) :: bn2o1(pcols,pverp) ! pressure factor for n2o
1063	real(r8), intent(in) :: bch4(pcols,pverp) ! pressure factor for ch4
1064	real(r8), intent(in) :: emplnk(14,pcols) ! emissivity Planck factor
1065	real(r8), intent(in) :: th2o(pcols) ! water vapor overlap factor
1066	real(r8), intent(in) :: tco2(pcols) ! co2 overlap factor
1067	!
1068	real(r8), intent(in) :: to3(pcols) ! o3 overlap factor
1069	real(r8), intent(in) :: s2c(pcols,pverp) ! h2o continuum path length
1070	real(r8), intent(in) :: w(pcols,pverp) ! h2o path length
1071	real(r8), intent(in) :: up2(pcols) ! pressure squared h2o path length
1072	!
1073	integer, intent(in) :: k ! level index
1074
1075	real(r8), intent(in) :: aer_trn_ttl(pcols,pverp,pverp,bnd_nbr_LW) ! aer trn.
1076
1077	!
1078	! Output Arguments
1079	!
1080	real(r8), intent(out) :: emstrc(pcols,pverp) ! total trace gas emissivity
1081
1082	!
1083	!--------------------------Local Variables------------------------------
1084	!
1085	integer i,l ! loop counters
1086	!
1087	real(r8) sqti(pcols) ! square root of mean temp
1088	real(r8) ecfc1 ! emissivity of cfc11 798 cm-1 band
1089	real(r8) ecfc2 ! " " " 846 cm-1 band
1090	real(r8) ecfc3 ! " " " 933 cm-1 band
1091	real(r8) ecfc4 ! " " " 1085 cm-1 band
1092	!
1093	real(r8) ecfc5 ! " " cfc12 889 cm-1 band
1094	real(r8) ecfc6 ! " " " 923 cm-1 band
1095	real(r8) ecfc7 ! " " " 1102 cm-1 band
1096	real(r8) ecfc8 ! " " " 1161 cm-1 band
1097	real(r8) u01 ! n2o path length
1098	!
1099	real(r8) u11 ! n2o path length
1100	real(r8) beta01 ! n2o pressure factor
1101	real(r8) beta11 ! n2o pressure factor
1102	real(r8) en2o1 ! emissivity of the 1285 cm-1 N2O band
1103	real(r8) u02 ! n2o path length
1104	!
1105	real(r8) u12 ! n2o path length
1106	real(r8) beta02 ! n2o pressure factor
1107	real(r8) en2o2 ! emissivity of the 589 cm-1 N2O band
1108	real(r8) u03 ! n2o path length
1109	real(r8) beta03 ! n2o pressure factor
1110	!
1111	real(r8) en2o3 ! emissivity of the 1168 cm-1 N2O band
1112	real(r8) betac ! ch4 pressure factor
1113	real(r8) ech4 ! emissivity of 1306 cm-1 CH4 band
1114	real(r8) betac1 ! co2 pressure factor
1115	real(r8) betac2 ! co2 pressure factor
1116	!
1117	real(r8) eco21 ! emissivity of 1064 cm-1 CO2 band
1118	real(r8) eco22 ! emissivity of 961 cm-1 CO2 band
1119	real(r8) tt(pcols) ! temp. factor for h2o overlap factor
1120	real(r8) psi1 ! narrow band h2o temp. factor
1121	real(r8) phi1 ! "
1122	!
1123	real(r8) p1 ! h2o line overlap factor
1124	real(r8) w1 ! "
1125	real(r8) tw(pcols,6) ! h2o transmission overlap
1126	real(r8) g1(6) ! h2o overlap factor
1127	real(r8) g2(6) ! "
1128	!
1129	real(r8) g3(6) ! "
1130	real(r8) g4(6) ! "
1131	real(r8) ab(6) ! "
1132	real(r8) bb(6) ! "
1133	real(r8) abp(6) ! "
1134	!
1135	real(r8) bbp(6) ! "
1136	real(r8) tcfc3 ! transmission for cfc11 band
1137	real(r8) tcfc4 ! "
1138	real(r8) tcfc6 ! transmission for cfc12 band
1139	real(r8) tcfc7 ! "
1140	!
1141	real(r8) tcfc8 ! "
1142	real(r8) tlw ! h2o overlap factor
1143	real(r8) tch4 ! ch4 overlap factor
1144	!
1145	!--------------------------Data Statements------------------------------
1146	!
1147	data g1 /0.0468556,0.0397454,0.0407664,0.0304380,0.0540398,0.0321962/
1148	data g2 /14.4832,4.30242,5.23523,3.25342,0.698935,16.5599/
1149	data g3 /26.1898,18.4476,15.3633,12.1927,9.14992,8.07092/
1150	data g4 /0.0261782,0.0369516,0.0307266,0.0243854,0.0182932,0.0161418/
1151	data ab /3.0857e-2,2.3524e-2,1.7310e-2,2.6661e-2,2.8074e-2,2.2915e-2/
1152	data bb /-1.3512e-4,-6.8320e-5,-3.2609e-5,-1.0228e-5,-9.5743e-5,-1.0304e-4/
1153	data abp/2.9129e-2,2.4101e-2,1.9821e-2,2.6904e-2,2.9458e-2,1.9892e-2/
1154	data bbp/-1.3139e-4,-5.5688e-5,-4.6380e-5,-8.0362e-5,-1.0115e-4,-8.8061e-5/
1155	!
1156	!--------------------------Statement Functions--------------------------
1157	!
1158	real(r8) func, u, b
1159	func(u,b) = u/sqrt(4.0 + u*(1.0 + 1.0 / b))
1160	!
1161	!-----------------------------------------------------------------------
1162	!
1163	do i = 1,ncol
1164	sqti(i) = sqrt(co2t(i,k))
1165	!
1166	! Transmission for h2o
1167	!
1168	tt(i) = abs(co2t(i,k) - 250.0)
1169	end do
1170	!
1171	do l = 1,6
1172	do i = 1,ncol
1173	psi1 = exp(abp(l)tt(i)+bbp(l)tt(i)*tt(i))
1174	phi1 = exp(ab(l)tt(i)+bb(l)tt(i)*tt(i))
1175	p1 = pnm(i,k) * (psi1/phi1) / sslp
1176	w1 = w(i,k) * phi1
1177	tw(i,l) = exp(- g1(l)p1(sqrt(1.0+g2(l)*(w1/p1))-1.0) &
1178	- g3(l)s2c(i,k)-g4(l)uptype(i,k))
1179	end do
1180	end do
1181
1182	! Overlap H2O tranmission with STRAER continuum in 6 trace gas
1183	! subbands
1184
1185	do i=1,ncol
1186	tw(i,1)=tw(i,1)(0.7aer_trn_ttl(i,k,1,idx_LW_0650_0800)+&! l=1: 0750--0820 cm-1
1187	0.3*aer_trn_ttl(i,k,1,idx_LW_0800_1000))
1188	tw(i,2)=tw(i,2)*aer_trn_ttl(i,k,1,idx_LW_0800_1000) ! l=2: 0820--0880 cm-1
1189	tw(i,3)=tw(i,3)*aer_trn_ttl(i,k,1,idx_LW_0800_1000) ! l=3: 0880--0900 cm-1
1190	tw(i,4)=tw(i,4)*aer_trn_ttl(i,k,1,idx_LW_0800_1000) ! l=4: 0900--1000 cm-1
1191	tw(i,5)=tw(i,5)*aer_trn_ttl(i,k,1,idx_LW_1000_1200) ! l=5: 1000--1120 cm-1
1192	tw(i,6)=tw(i,6)*aer_trn_ttl(i,k,1,idx_LW_1000_1200) ! l=6: 1120--1170 cm-1
1193	end do ! end loop over lon
1194	!
1195	do i = 1,ncol
1196	!
1197	! transmission due to cfc bands
1198	!
1199	tcfc3 = exp(-175.005*ucfc11(i,k))
1200	tcfc4 = exp(-1202.18*ucfc11(i,k))
1201	tcfc6 = exp(-5786.73*ucfc12(i,k))
1202	tcfc7 = exp(-2873.51*ucfc12(i,k))
1203	tcfc8 = exp(-2085.59*ucfc12(i,k))
1204	!
1205	! Emissivity for CFC11 bands
1206	!
1207	ecfc1 = 50.0(1.0 - exp(-54.09ucfc11(i,k))) * tw(i,1) * emplnk(7,i)
1208	ecfc2 = 60.0(1.0 - exp(-5130.03ucfc11(i,k)))* tw(i,2) * emplnk(8,i)
1209	ecfc3 = 60.0(1.0 - tcfc3)tw(i,4)tcfc6emplnk(9,i)
1210	ecfc4 = 100.0(1.0 - tcfc4)tw(i,5)*emplnk(10,i)
1211	!
1212	! Emissivity for CFC12 bands
1213	!
1214	ecfc5 = 45.0(1.0 - exp(-1272.35ucfc12(i,k)))tw(i,3)emplnk(11,i)
1215	ecfc6 = 50.0(1.0 - tcfc6)tw(i,4)*emplnk(12,i)
1216	ecfc7 = 80.0(1.0 - tcfc7)tw(i,5)* tcfc4 * emplnk(13,i)
1217	ecfc8 = 70.0(1.0 - tcfc8)tw(i,6) * emplnk(14,i)
1218	!
1219	! Emissivity for CH4 band 1306 cm-1
1220	!
1221	tlw = exp(-1.0*sqrt(up2(i)))
1222
1223	! Overlap H2O vibration rotation band with STRAER continuum
1224	! for CH4 1306 cm-1 and N2O 1285 cm-1 bands
1225
1226	tlw=tlw*aer_trn_ttl(i,k,1,idx_LW_1200_2000)
1227	betac = bch4(i,k)/uch4(i,k)
1228	ech4 = 6.00444sqti(i)log(1.0 + func(uch4(i,k),betac)) tlw emplnk(3,i)
1229	tch4 = 1.0/(1.0 + 0.02*func(uch4(i,k),betac))
1230	!
1231	! Emissivity for N2O bands
1232	!
1233	u01 = un2o0(i,k)
1234	u11 = un2o1(i,k)
1235	beta01 = bn2o0(i,k)/un2o0(i,k)
1236	beta11 = bn2o1(i,k)/un2o1(i,k)
1237	!
1238	! 1285 cm-1 band
1239	!
1240	en2o1 = 2.35558sqti(i)log(1.0 + func(u01,beta01) + &
1241	func(u11,beta11))tlwtch4*emplnk(4,i)
1242	u02 = 0.100090*u01
1243	u12 = 0.0992746*u11
1244	beta02 = 0.964282*beta01
1245	!
1246	! 589 cm-1 band
1247	!
1248	en2o2 = 2.65581sqti(i)log(1.0 + func(u02,beta02) + &
1249	func(u12,beta02)) * tco2(i) * th2o(i) * emplnk(5,i)
1250	u03 = 0.0333767*u01
1251	beta03 = 0.982143*beta01
1252	!
1253	! 1168 cm-1 band
1254	!
1255	en2o3 = 2.54034sqti(i)log(1.0 + func(u03,beta03)) * &
1256	tw(i,6) * tcfc8 * emplnk(6,i)
1257	!
1258	! Emissivity for 1064 cm-1 band of CO2
1259	!
1260	betac1 = 2.97558pnm(i,k) / (sslpsqti(i))
1261	betac2 = 2.0 * betac1
1262	eco21 = 3.7571sqti(i)log(1.0 + func(uco211(i,k),betac1) &
1263	+ func(uco212(i,k),betac2) + func(uco213(i,k),betac2)) &
1264	* to3(i) * tw(i,5) * tcfc4 * tcfc7 * emplnk(2,i)
1265	!
1266	! Emissivity for 961 cm-1 band
1267	!
1268	eco22 = 3.8443sqti(i)log(1.0 + func(uco221(i,k),betac1) &
1269	+ func(uco222(i,k),betac1) + func(uco223(i,k),betac2)) &
1270	* tw(i,4) * tcfc3 * tcfc6 * emplnk(1,i)
1271	!
1272	! total trace gas emissivity
1273	!
1274	emstrc(i,k) = ecfc1 + ecfc2 + ecfc3 + ecfc4 + ecfc5 +ecfc6 + &
1275	ecfc7 + ecfc8 + en2o1 + en2o2 + en2o3 + ech4 + &
1276	eco21 + eco22
1277	end do
1278	!
1279	return
1280	!
1281	end subroutine trcems
1282
1283	subroutine trcmix(lchnk ,ncol ,pcols, pver, &
1284	pmid ,clat, n2o ,ch4 , &
1285	cfc11 , cfc12 )
1286	!-----------------------------------------------------------------------
1287	!
1288	! Purpose:
1289	! Specify zonal mean mass mixing ratios of CH4, N2O, CFC11 and
1290	! CFC12
1291	!
1292	! Method:
1293	! Distributions assume constant mixing ratio in the troposphere
1294	! and a decrease of mixing ratio in the stratosphere. Tropopause
1295	! defined by ptrop. The scale height of the particular trace gas
1296	! depends on latitude. This assumption produces a more realistic
1297	! stratospheric distribution of the various trace gases.
1298	!
1299	! Author: J. Kiehl
1300	!
1301	!-----------------------------------------------------------------------
1302	! use shr_kind_mod, only: r8 => shr_kind_r8
1303	! use ppgrid
1304	! use phys_grid, only: get_rlat_all_p
1305	! use physconst, only: mwdry, mwch4, mwn2o, mwf11, mwf12
1306	! use ghg_surfvals, only: ch4vmr, n2ovmr, f11vmr, f12vmr
1307
1308	implicit none
1309
1310	!-----------------------------Arguments---------------------------------
1311	!
1312	! Input
1313	!
1314	integer, intent(in) :: lchnk ! chunk identifier
1315	integer, intent(in) :: ncol ! number of atmospheric columns
1316	integer, intent(in) :: pcols, pver
1317
1318	real(r8), intent(in) :: pmid(pcols,pver) ! model pressures
1319	real(r8), intent(in) :: clat(pcols) ! latitude in radians for columns
1320	!
1321	! Output
1322	!
1323	real(r8), intent(out) :: n2o(pcols,pver) ! nitrous oxide mass mixing ratio
1324	real(r8), intent(out) :: ch4(pcols,pver) ! methane mass mixing ratio
1325	real(r8), intent(out) :: cfc11(pcols,pver) ! cfc11 mass mixing ratio
1326	real(r8), intent(out) :: cfc12(pcols,pver) ! cfc12 mass mixing ratio
1327
1328	!
1329	!--------------------------Local Variables------------------------------
1330
1331	real(r8) :: rmwn2o ! ratio of molecular weight n2o to dry air
1332	real(r8) :: rmwch4 ! ratio of molecular weight ch4 to dry air
1333	real(r8) :: rmwf11 ! ratio of molecular weight cfc11 to dry air
1334	real(r8) :: rmwf12 ! ratio of molecular weight cfc12 to dry air
1335	!
1336	integer i ! longitude loop index
1337	integer k ! level index
1338	!
1339	! real(r8) clat(pcols) ! latitude in radians for columns
1340	real(r8) coslat(pcols) ! cosine of latitude
1341	real(r8) dlat ! latitude in degrees
1342	real(r8) ptrop ! pressure level of tropopause
1343	real(r8) pratio ! pressure divided by ptrop
1344	!
1345	real(r8) xn2o ! pressure scale height for n2o
1346	real(r8) xch4 ! pressure scale height for ch4
1347	real(r8) xcfc11 ! pressure scale height for cfc11
1348	real(r8) xcfc12 ! pressure scale height for cfc12
1349	!
1350	real(r8) ch40 ! tropospheric mass mixing ratio for ch4
1351	real(r8) n2o0 ! tropospheric mass mixing ratio for n2o
1352	real(r8) cfc110 ! tropospheric mass mixing ratio for cfc11
1353	real(r8) cfc120 ! tropospheric mass mixing ratio for cfc12
1354	!
1355	!-----------------------------------------------------------------------
1356	rmwn2o = mwn2o/mwdry ! ratio of molecular weight n2o to dry air
1357	rmwch4 = mwch4/mwdry ! ratio of molecular weight ch4 to dry air
1358	rmwf11 = mwf11/mwdry ! ratio of molecular weight cfc11 to dry air
1359	rmwf12 = mwf12/mwdry ! ratio of molecular weight cfc12 to dry air
1360	!
1361	! get latitudes
1362	!
1363	! call get_rlat_all_p(lchnk, ncol, clat)
1364	do i = 1, ncol
1365	coslat(i) = cos(clat(i))
1366	end do
1367	!
1368	! set tropospheric mass mixing ratios
1369	!
1370	ch40 = rmwch4 * ch4vmr
1371	n2o0 = rmwn2o * n2ovmr
1372	cfc110 = rmwf11 * f11vmr
1373	cfc120 = rmwf12 * f12vmr
1374
1375	do i = 1, ncol
1376	coslat(i) = cos(clat(i))
1377	end do
1378	!
1379	do k = 1,pver
1380	do i = 1,ncol
1381	!
1382	! set stratospheric scale height factor for gases
1383	dlat = abs(57.2958 * clat(i))
1384	if(dlat.le.45.0) then
1385	xn2o = 0.3478 + 0.00116 * dlat
1386	xch4 = 0.2353
1387	xcfc11 = 0.7273 + 0.00606 * dlat
1388	xcfc12 = 0.4000 + 0.00222 * dlat
1389	else
1390	xn2o = 0.4000 + 0.013333 * (dlat - 45)
1391	xch4 = 0.2353 + 0.0225489 * (dlat - 45)
1392	xcfc11 = 1.00 + 0.013333 * (dlat - 45)
1393	xcfc12 = 0.50 + 0.024444 * (dlat - 45)
1394	end if
1395	!
1396	! pressure of tropopause
1397	ptrop = 250.0e2 - 150.0e2coslat(i)*2.0
1398	!
1399	! determine output mass mixing ratios
1400	if (pmid(i,k) >= ptrop) then
1401	ch4(i,k) = ch40
1402	n2o(i,k) = n2o0
1403	cfc11(i,k) = cfc110
1404	cfc12(i,k) = cfc120
1405	else
1406	pratio = pmid(i,k)/ptrop
1407	ch4(i,k) = ch40 * (pratio)**xch4
1408	n2o(i,k) = n2o0 * (pratio)**xn2o
1409	cfc11(i,k) = cfc110 * (pratio)**xcfc11
1410	cfc12(i,k) = cfc120 * (pratio)**xcfc12
1411	end if
1412	end do
1413	end do
1414	!
1415	return
1416	!
1417	end subroutine trcmix
1418
1419	subroutine trcplk(lchnk ,ncol ,pcols, pver, pverp, &
1420	tint ,tlayr ,tplnke ,emplnk ,abplnk1 , &
1421	abplnk2 )
1422	!-----------------------------------------------------------------------
1423	!
1424	! Purpose:
1425	! Calculate Planck factors for absorptivity and emissivity of
1426	! CH4, N2O, CFC11 and CFC12
1427	!
1428	! Method:
1429	! Planck function and derivative evaluated at the band center.
1430	!
1431	! Author: J. Kiehl
1432	!
1433	!-----------------------------------------------------------------------
1434	! use shr_kind_mod, only: r8 => shr_kind_r8
1435	! use ppgrid
1436
1437	implicit none
1438	!------------------------------Arguments--------------------------------
1439	!
1440	! Input arguments
1441	!
1442	integer, intent(in) :: lchnk ! chunk identifier
1443	integer, intent(in) :: ncol ! number of atmospheric columns
1444	integer, intent(in) :: pcols, pver, pverp
1445
1446	real(r8), intent(in) :: tint(pcols,pverp) ! interface temperatures
1447	real(r8), intent(in) :: tlayr(pcols,pverp) ! k-1 level temperatures
1448	real(r8), intent(in) :: tplnke(pcols) ! Top Layer temperature
1449	!
1450	! output arguments
1451	!
1452	real(r8), intent(out) :: emplnk(14,pcols) ! emissivity Planck factor
1453	real(r8), intent(out) :: abplnk1(14,pcols,pverp) ! non-nearest layer Plack factor
1454	real(r8), intent(out) :: abplnk2(14,pcols,pverp) ! nearest layer factor
1455
1456	!
1457	!--------------------------Local Variables------------------------------
1458	!
1459	integer wvl ! wavelength index
1460	integer i,k ! loop counters
1461	!
1462	real(r8) f1(14) ! Planck function factor
1463	real(r8) f2(14) ! "
1464	real(r8) f3(14) ! "
1465	!
1466	!--------------------------Data Statements------------------------------
1467	!
1468	data f1 /5.85713e8,7.94950e8,1.47009e9,1.40031e9,1.34853e8, &
1469	1.05158e9,3.35370e8,3.99601e8,5.35994e8,8.42955e8, &
1470	4.63682e8,5.18944e8,8.83202e8,1.03279e9/
1471	data f2 /2.02493e11,3.04286e11,6.90698e11,6.47333e11, &
1472	2.85744e10,4.41862e11,9.62780e10,1.21618e11, &
1473	1.79905e11,3.29029e11,1.48294e11,1.72315e11, &
1474	3.50140e11,4.31364e11/
1475	data f3 /1383.0,1531.0,1879.0,1849.0,848.0,1681.0, &
1476	1148.0,1217.0,1343.0,1561.0,1279.0,1328.0, &
1477	1586.0,1671.0/
1478	!
1479	!-----------------------------------------------------------------------
1480	!
1481	! Calculate emissivity Planck factor
1482	!
1483	do wvl = 1,14
1484	do i = 1,ncol
1485	emplnk(wvl,i) = f1(wvl)/(tplnke(i)*4.0(exp(f3(wvl)/tplnke(i))-1.0))
1486	end do
1487	end do
1488	!
1489	! Calculate absorptivity Planck factor for tint and tlayr temperatures
1490	!
1491	do wvl = 1,14
1492	do k = ntoplw, pverp
1493	do i = 1, ncol
1494	!
1495	! non-nearlest layer function
1496	!
1497	abplnk1(wvl,i,k) = (f2(wvl)*exp(f3(wvl)/tint(i,k))) &
1498	/(tint(i,k)*5.0(exp(f3(wvl)/tint(i,k))-1.0)**2.0)
1499	!
1500	! nearest layer function
1501	!
1502	abplnk2(wvl,i,k) = (f2(wvl)*exp(f3(wvl)/tlayr(i,k))) &
1503	/(tlayr(i,k)*5.0(exp(f3(wvl)/tlayr(i,k))-1.0)**2.0)
1504	end do
1505	end do
1506	end do
1507	!
1508	return
1509	end subroutine trcplk
1510
1511	subroutine trcpth(lchnk ,ncol ,pcols, pver, pverp, &
1512	tnm ,pnm ,cfc11 ,cfc12 ,n2o , &
1513	ch4 ,qnm ,ucfc11 ,ucfc12 ,un2o0 , &
1514	un2o1 ,uch4 ,uco211 ,uco212 ,uco213 , &
1515	uco221 ,uco222 ,uco223 ,bn2o0 ,bn2o1 , &
1516	bch4 ,uptype )
1517	!-----------------------------------------------------------------------
1518	!
1519	! Purpose:
1520	! Calculate path lengths and pressure factors for CH4, N2O, CFC11
1521	! and CFC12.
1522	!
1523	! Method:
1524	! See CCM3 description for details
1525	!
1526	! Author: J. Kiehl
1527	!
1528	!-----------------------------------------------------------------------
1529	! use shr_kind_mod, only: r8 => shr_kind_r8
1530	! use ppgrid
1531	! use ghg_surfvals, only: co2mmr
1532
1533	implicit none
1534
1535	!------------------------------Arguments--------------------------------
1536	!
1537	! Input arguments
1538	!
1539	integer, intent(in) :: lchnk ! chunk identifier
1540	integer, intent(in) :: ncol ! number of atmospheric columns
1541	integer, intent(in) :: pcols, pver, pverp
1542
1543	real(r8), intent(in) :: tnm(pcols,pver) ! Model level temperatures
1544	real(r8), intent(in) :: pnm(pcols,pverp) ! Pres. at model interfaces (dynes/cm2)
1545	real(r8), intent(in) :: qnm(pcols,pver) ! h2o specific humidity
1546	real(r8), intent(in) :: cfc11(pcols,pver) ! CFC11 mass mixing ratio
1547	!
1548	real(r8), intent(in) :: cfc12(pcols,pver) ! CFC12 mass mixing ratio
1549	real(r8), intent(in) :: n2o(pcols,pver) ! N2O mass mixing ratio
1550	real(r8), intent(in) :: ch4(pcols,pver) ! CH4 mass mixing ratio
1551
1552	!
1553	! Output arguments
1554	!
1555	real(r8), intent(out) :: ucfc11(pcols,pverp) ! CFC11 path length
1556	real(r8), intent(out) :: ucfc12(pcols,pverp) ! CFC12 path length
1557	real(r8), intent(out) :: un2o0(pcols,pverp) ! N2O path length
1558	real(r8), intent(out) :: un2o1(pcols,pverp) ! N2O path length (hot band)
1559	real(r8), intent(out) :: uch4(pcols,pverp) ! CH4 path length
1560	!
1561	real(r8), intent(out) :: uco211(pcols,pverp) ! CO2 9.4 micron band path length
1562	real(r8), intent(out) :: uco212(pcols,pverp) ! CO2 9.4 micron band path length
1563	real(r8), intent(out) :: uco213(pcols,pverp) ! CO2 9.4 micron band path length
1564	real(r8), intent(out) :: uco221(pcols,pverp) ! CO2 10.4 micron band path length
1565	real(r8), intent(out) :: uco222(pcols,pverp) ! CO2 10.4 micron band path length
1566	!
1567	real(r8), intent(out) :: uco223(pcols,pverp) ! CO2 10.4 micron band path length
1568	real(r8), intent(out) :: bn2o0(pcols,pverp) ! pressure factor for n2o
1569	real(r8), intent(out) :: bn2o1(pcols,pverp) ! pressure factor for n2o
1570	real(r8), intent(out) :: bch4(pcols,pverp) ! pressure factor for ch4
1571	real(r8), intent(out) :: uptype(pcols,pverp) ! p-type continuum path length
1572
1573	!
1574	!---------------------------Local variables-----------------------------
1575	!
1576	integer i ! Longitude index
1577	integer k ! Level index
1578	!
1579	real(r8) co2fac(pcols,1) ! co2 factor
1580	real(r8) alpha1(pcols) ! stimulated emission term
1581	real(r8) alpha2(pcols) ! stimulated emission term
1582	real(r8) rt(pcols) ! reciprocal of local temperature
1583	real(r8) rsqrt(pcols) ! reciprocal of sqrt of temp
1584	!
1585	real(r8) pbar(pcols) ! mean pressure
1586	real(r8) dpnm(pcols) ! difference in pressure
1587	real(r8) diff ! diffusivity factor
1588	!
1589	!--------------------------Data Statements------------------------------
1590	!
1591	data diff /1.66/
1592	!
1593	!-----------------------------------------------------------------------
1594	!
1595	! Calculate path lengths for the trace gases at model top
1596	!
1597	do i = 1,ncol
1598	ucfc11(i,ntoplw) = 1.8 * cfc11(i,ntoplw) * pnm(i,ntoplw) * rga
1599	ucfc12(i,ntoplw) = 1.8 * cfc12(i,ntoplw) * pnm(i,ntoplw) * rga
1600	un2o0(i,ntoplw) = diff * 1.02346e5 * n2o(i,ntoplw) * pnm(i,ntoplw) * rga / sqrt(tnm(i,ntoplw))
1601	un2o1(i,ntoplw) = diff * 2.01909 * un2o0(i,ntoplw) * exp(-847.36/tnm(i,ntoplw))
1602	uch4(i,ntoplw) = diff * 8.60957e4 * ch4(i,ntoplw) * pnm(i,ntoplw) * rga / sqrt(tnm(i,ntoplw))
1603	co2fac(i,1) = diff * co2mmr * pnm(i,ntoplw) * rga
1604	alpha1(i) = (1.0 - exp(-1540.0/tnm(i,ntoplw)))**3.0/sqrt(tnm(i,ntoplw))
1605	alpha2(i) = (1.0 - exp(-1360.0/tnm(i,ntoplw)))**3.0/sqrt(tnm(i,ntoplw))
1606	uco211(i,ntoplw) = 3.42217e3 * co2fac(i,1) * alpha1(i) * exp(-1849.7/tnm(i,ntoplw))
1607	uco212(i,ntoplw) = 6.02454e3 * co2fac(i,1) * alpha1(i) * exp(-2782.1/tnm(i,ntoplw))
1608	uco213(i,ntoplw) = 5.53143e3 * co2fac(i,1) * alpha1(i) * exp(-3723.2/tnm(i,ntoplw))
1609	uco221(i,ntoplw) = 3.88984e3 * co2fac(i,1) * alpha2(i) * exp(-1997.6/tnm(i,ntoplw))
1610	uco222(i,ntoplw) = 3.67108e3 * co2fac(i,1) * alpha2(i) * exp(-3843.8/tnm(i,ntoplw))
1611	uco223(i,ntoplw) = 6.50642e3 * co2fac(i,1) * alpha2(i) * exp(-2989.7/tnm(i,ntoplw))
1612	bn2o0(i,ntoplw) = diff * 19.399 * pnm(i,ntoplw)*2.0 n2o(i,ntoplw) * &
1613	1.02346e5 * rga / (sslp*tnm(i,ntoplw))
1614	bn2o1(i,ntoplw) = bn2o0(i,ntoplw) * exp(-847.36/tnm(i,ntoplw)) * 2.06646e5
1615	bch4(i,ntoplw) = diff * 2.94449 * ch4(i,ntoplw) * pnm(i,ntoplw)*2.0 rga * &
1616	8.60957e4 / (sslp*tnm(i,ntoplw))
1617	uptype(i,ntoplw) = diff * qnm(i,ntoplw) * pnm(i,ntoplw)*2.0 &
1618	exp(1800.0(1.0/tnm(i,ntoplw) - 1.0/296.0)) rga / sslp
1619	end do
1620	!
1621	! Calculate trace gas path lengths through model atmosphere
1622	!
1623	do k = ntoplw,pver
1624	do i = 1,ncol
1625	rt(i) = 1./tnm(i,k)
1626	rsqrt(i) = sqrt(rt(i))
1627	pbar(i) = 0.5 * (pnm(i,k+1) + pnm(i,k)) / sslp
1628	dpnm(i) = (pnm(i,k+1) - pnm(i,k)) * rga
1629	alpha1(i) = diff * rsqrt(i) * (1.0 - exp(-1540.0/tnm(i,k)))**3.0
1630	alpha2(i) = diff * rsqrt(i) * (1.0 - exp(-1360.0/tnm(i,k)))**3.0
1631	ucfc11(i,k+1) = ucfc11(i,k) + 1.8 * cfc11(i,k) * dpnm(i)
1632	ucfc12(i,k+1) = ucfc12(i,k) + 1.8 * cfc12(i,k) * dpnm(i)
1633	un2o0(i,k+1) = un2o0(i,k) + diff * 1.02346e5 * n2o(i,k) * rsqrt(i) * dpnm(i)
1634	un2o1(i,k+1) = un2o1(i,k) + diff * 2.06646e5 * n2o(i,k) * &
1635	rsqrt(i) * exp(-847.36/tnm(i,k)) * dpnm(i)
1636	uch4(i,k+1) = uch4(i,k) + diff * 8.60957e4 * ch4(i,k) * rsqrt(i) * dpnm(i)
1637	uco211(i,k+1) = uco211(i,k) + 1.153.42217e3 alpha1(i) * &
1638	co2mmr * exp(-1849.7/tnm(i,k)) * dpnm(i)
1639	uco212(i,k+1) = uco212(i,k) + 1.156.02454e3 alpha1(i) * &
1640	co2mmr * exp(-2782.1/tnm(i,k)) * dpnm(i)
1641	uco213(i,k+1) = uco213(i,k) + 1.155.53143e3 alpha1(i) * &
1642	co2mmr * exp(-3723.2/tnm(i,k)) * dpnm(i)
1643	uco221(i,k+1) = uco221(i,k) + 1.153.88984e3 alpha2(i) * &
1644	co2mmr * exp(-1997.6/tnm(i,k)) * dpnm(i)
1645	uco222(i,k+1) = uco222(i,k) + 1.153.67108e3 alpha2(i) * &
1646	co2mmr * exp(-3843.8/tnm(i,k)) * dpnm(i)
1647	uco223(i,k+1) = uco223(i,k) + 1.156.50642e3 alpha2(i) * &
1648	co2mmr * exp(-2989.7/tnm(i,k)) * dpnm(i)
1649	bn2o0(i,k+1) = bn2o0(i,k) + diff * 19.399 * pbar(i) * rt(i) &
1650	* 1.02346e5 * n2o(i,k) * dpnm(i)
1651	bn2o1(i,k+1) = bn2o1(i,k) + diff * 19.399 * pbar(i) * rt(i) &
1652	* 2.06646e5 * exp(-847.36/tnm(i,k)) * n2o(i,k)*dpnm(i)
1653	bch4(i,k+1) = bch4(i,k) + diff * 2.94449 * rt(i) * pbar(i) &
1654	* 8.60957e4 * ch4(i,k) * dpnm(i)
1655	uptype(i,k+1) = uptype(i,k) + diff qnm(i,k) &
1656	exp(1800.0(1.0/tnm(i,k) - 1.0/296.0)) pbar(i) * dpnm(i)
1657	end do
1658	end do
1659	!
1660	return
1661	end subroutine trcpth
1662
1663
1664
1665	subroutine aqsat(t ,p ,es ,qs ,ii , &
1666	ilen ,kk ,kstart ,kend )
1667	!-----------------------------------------------------------------------
1668	!
1669	! Purpose:
1670	! Utility procedure to look up and return saturation vapor pressure from
1671	! precomputed table, calculate and return saturation specific humidity
1672	! (g/g),for input arrays of temperature and pressure (dimensioned ii,kk)
1673	! This routine is useful for evaluating only a selected region in the
1674	! vertical.
1675	!
1676	! Method:
1677	! <Describe the algorithm(s) used in the routine.>
1678	! <Also include any applicable external references.>
1679	!
1680	! Author: J. Hack
1681	!
1682	!------------------------------Arguments--------------------------------
1683	!
1684	! Input arguments
1685	!
1686	integer, intent(in) :: ii ! I dimension of arrays t, p, es, qs
1687	integer, intent(in) :: kk ! K dimension of arrays t, p, es, qs
1688	integer, intent(in) :: ilen ! Length of vectors in I direction which
1689	integer, intent(in) :: kstart ! Starting location in K direction
1690	integer, intent(in) :: kend ! Ending location in K direction
1691	real(r8), intent(in) :: t(ii,kk) ! Temperature
1692	real(r8), intent(in) :: p(ii,kk) ! Pressure
1693	!
1694	! Output arguments
1695	!
1696	real(r8), intent(out) :: es(ii,kk) ! Saturation vapor pressure
1697	real(r8), intent(out) :: qs(ii,kk) ! Saturation specific humidity
1698	!
1699	!---------------------------Local workspace-----------------------------
1700	!
1701	real(r8) omeps ! 1 - 0.622
1702	integer i, k ! Indices
1703	!
1704	!-----------------------------------------------------------------------
1705	!
1706	omeps = 1.0 - epsqs
1707	do k=kstart,kend
1708	do i=1,ilen
1709	es(i,k) = estblf(t(i,k))
1710	!
1711	! Saturation specific humidity
1712	!
1713	qs(i,k) = epsqses(i,k)/(p(i,k) - omepses(i,k))
1714	!
1715	! The following check is to avoid the generation of negative values
1716	! that can occur in the upper stratosphere and mesosphere
1717	!
1718	qs(i,k) = min(1.0_r8,qs(i,k))
1719	!
1720	if (qs(i,k) < 0.0) then
1721	qs(i,k) = 1.0
1722	es(i,k) = p(i,k)
1723	end if
1724	end do
1725	end do
1726	!
1727	return
1728	end subroutine aqsat
1729	!===============================================================================
1730	subroutine cldefr(lchnk ,ncol ,pcols, pver, pverp, &
1731	landfrac,t ,rel ,rei ,ps ,pmid , landm, icefrac, snowh)
1732	!-----------------------------------------------------------------------
1733	!
1734	! Purpose:
1735	! Compute cloud water and ice particle size
1736	!
1737	! Method:
1738	! use empirical formulas to construct effective radii
1739	!
1740	! Author: J.T. Kiehl, B. A. Boville, P. Rasch
1741	!
1742	!-----------------------------------------------------------------------
1743
1744	implicit none
1745	!------------------------------Arguments--------------------------------
1746	!
1747	! Input arguments
1748	!
1749	integer, intent(in) :: lchnk ! chunk identifier
1750	integer, intent(in) :: ncol ! number of atmospheric columns
1751	integer, intent(in) :: pcols, pver, pverp
1752
1753	real(r8), intent(in) :: landfrac(pcols) ! Land fraction
1754	real(r8), intent(in) :: icefrac(pcols) ! Ice fraction
1755	real(r8), intent(in) :: t(pcols,pver) ! Temperature
1756	real(r8), intent(in) :: ps(pcols) ! Surface pressure
1757	real(r8), intent(in) :: pmid(pcols,pver) ! Midpoint pressures
1758	real(r8), intent(in) :: landm(pcols)
1759	real(r8), intent(in) :: snowh(pcols) ! Snow depth over land, water equivalent (m)
1760	!
1761	! Output arguments
1762	!
1763	real(r8), intent(out) :: rel(pcols,pver) ! Liquid effective drop size (microns)
1764	real(r8), intent(out) :: rei(pcols,pver) ! Ice effective drop size (microns)
1765	!
1766
1767	!++pjr
1768	! following Kiehl
1769	call reltab(ncol, pcols, pver, t, landfrac, landm, icefrac, rel, snowh)
1770
1771	! following Kristjansson and Mitchell
1772	call reitab(ncol, pcols, pver, t, rei)
1773	!--pjr
1774	!
1775	!
1776	return
1777	end subroutine cldefr
1778
1779
1780	subroutine background(lchnk, ncol, pint, pcols, pverr, pverrp, mmr)
1781	!-----------------------------------------------------------------------
1782	!
1783	! Purpose:
1784	! Set global mean tropospheric aerosol background (or tuning) field
1785	!
1786	! Method:
1787	! Specify aerosol mixing ratio.
1788	! Aerosol mass mixing ratio
1789	! is specified so that the column visible aerosol optical depth is a
1790	! specified global number (tauback). This means that the actual mixing
1791	! ratio depends on pressure thickness of the lowest three atmospheric
1792	! layers near the surface.
1793	!
1794	!-----------------------------------------------------------------------
1795	! use shr_kind_mod, only: r8 => shr_kind_r8
1796	! use aer_optics, only: kbg,idxVIS
1797	! use physconst, only: gravit
1798	!-----------------------------------------------------------------------
1799	implicit none
1800	!-----------------------------------------------------------------------
1801	!#include <ptrrgrid.h>
1802	!------------------------------Arguments--------------------------------
1803	!
1804	! Input arguments
1805	!
1806	integer, intent(in) :: lchnk ! chunk identifier
1807	integer, intent(in) :: ncol ! number of atmospheric columns
1808	integer, intent(in) :: pcols,pverr,pverrp
1809
1810	real(r8), intent(in) :: pint(pcols,pverrp) ! Interface pressure (mks)
1811	!
1812	! Output arguments
1813	!
1814	real(r8), intent(out) :: mmr(pcols,pverr) ! "background" aerosol mass mixing ratio
1815	!
1816	!---------------------------Local variables-----------------------------
1817	!
1818	integer i ! Longitude index
1819	integer k ! Level index
1820	!
1821	real(r8) mass2mmr ! Factor to convert mass to mass mixing ratio
1822	real(r8) mass ! Mass of "background" aerosol as specified by tauback
1823	!
1824	!-----------------------------------------------------------------------
1825	!
1826	do i=1,ncol
1827	mass2mmr = gravmks / (pint(i,pverrp)-pint(i,pverrp-mxaerl))
1828	do k=1,pverr
1829	!
1830	! Compute aerosol mass mixing ratio for specified levels (1.e3 factor is
1831	! for units conversion of the extinction coefficiant from m2/g to m2/kg)
1832	!
1833	if ( k >= pverrp-mxaerl ) then
1834	! kaervs is not consistent with the values in aer_optics
1835	! this ?should? be changed.
1836	! rhfac is also implemented differently
1837	mass = tauback / (1.e3 * kbg(idxVIS))
1838	mmr(i,k) = mass2mmr*mass
1839	else
1840	mmr(i,k) = 0._r8
1841	endif
1842	!
1843	enddo
1844	enddo
1845	!
1846	return
1847	end subroutine background
1848
1849	subroutine scale_aerosols(AEROSOLt, pcols, pver, ncol, lchnk, scale)
1850	!-----------------------------------------------------------------
1851	! scale each species as determined by scale factors
1852	!-----------------------------------------------------------------
1853	integer, intent(in) :: ncol, lchnk ! number of columns and chunk index
1854	integer, intent(in) :: pcols, pver
1855	real(r8), intent(in) :: scale(naer_all) ! scale each aerosol by this amount
1856	real(r8), intent(inout) :: AEROSOLt(pcols, pver, naer_all) ! aerosols
1857	integer m
1858
1859	do m = 1, naer_all
1860	AEROSOLt(:ncol, :, m) = scale(m)*AEROSOLt(:ncol, :, m)
1861	end do
1862
1863	return
1864	end subroutine scale_aerosols
1865
1866	subroutine get_int_scales(scales)
1867	real(r8), intent(out)::scales(naer_all) ! scale each aerosol by this amount
1868	integer i ! index through species
1869
1870	!initialize
1871	scales = 1.
1872
1873	scales(idxBG) = 1._r8
1874	scales(idxSUL) = sulscl
1875	scales(idxSSLT) = ssltscl
1876
1877	do i = idxCARBONfirst, idxCARBONfirst+numCARBON-1
1878	scales(i) = carscl
1879	enddo
1880
1881	do i = idxDUSTfirst, idxDUSTfirst+numDUST-1
1882	scales(i) = dustscl
1883	enddo
1884
1885	scales(idxVOLC) = volcscl
1886
1887	return
1888	end subroutine get_int_scales
1889
1890	subroutine vert_interpolate (Match_ps, aerosolc, m_hybi, paerlev, naer_c, pint, n, AEROSOL_mmr, pcols, pver, pverp, ncol, c)
1891	!--------------------------------------------------------------------
1892	! Input: match surface pressure, cam interface pressure,
1893	! month index, number of columns, chunk index
1894	!
1895	! Output: Aerosol mass mixing ratio (AEROSOL_mmr)
1896	!
1897	! Method:
1898	! interpolate column mass (cumulative) from match onto
1899	! cam's vertical grid (pressure coordinate)
1900	! convert back to mass mixing ratio
1901	!
1902	!--------------------------------------------------------------------
1903
1904	! use physconst, only: gravit
1905
1906	integer, intent(in) :: paerlev,naer_c,pcols,pver,pverp
1907	real(r8), intent(out) :: AEROSOL_mmr(pcols,pver,naer) ! aerosol mmr from MATCH
1908	real(r8), intent(in) :: Match_ps(pcols) ! surface pressure at a particular month
1909	real(r8), intent(in) :: pint(pcols,pverp) ! interface pressure from CAM
1910	real(r8), intent(in) :: aerosolc(pcols,paerlev,naer_c)
1911	real(r8), intent(in) :: m_hybi(paerlev)
1912
1913	integer, intent(in) :: ncol,c ! chunk index and number of columns
1914	integer, intent(in) :: n ! prv or nxt month index
1915	!
1916	! Local workspace
1917	!
1918	integer m ! index to aerosol species
1919	integer kupper(pcols) ! last upper bound for interpolation
1920	integer i, k, kk, kkstart, kount ! loop vars for interpolation
1921	integer isv, ksv, msv ! loop indices to save
1922
1923	logical bad ! indicates a bad point found
1924	logical lev_interp_comp ! interpolation completed for a level
1925
1926	real(r8) AEROSOL(pcols,pverp,naer) ! cumulative mass of aerosol in column beneath upper
1927	! interface of level in column at particular month
1928	real(r8) dpl, dpu ! lower and upper intepolation factors
1929	real(r8) v_coord ! vertical coordinate
1930	real(r8) m_to_mmr ! mass to mass mixing ratio conversion factor
1931	real(r8) AER_diff ! temp var for difference between aerosol masses
1932
1933	! call t_startf ('vert_interpolate')
1934	!
1935	! Initialize index array
1936	!
1937	do i=1,ncol
1938	kupper(i) = 1
1939	end do
1940	!
1941	! assign total mass to topmost level
1942	!
1943
1944	do i=1,ncol
1945	do m=1,naer
1946	AEROSOL(i,1,m) = AEROSOLc(i,1,m)
1947	enddo
1948	enddo
1949	!
1950	! At every pressure level, interpolate onto that pressure level
1951	!
1952	do k=2,pver
1953	!
1954	! Top level we need to start looking is the top level for the previous k
1955	! for all longitude points
1956	!
1957	kkstart = paerlev
1958	do i=1,ncol
1959	kkstart = min0(kkstart,kupper(i))
1960	end do
1961	kount = 0
1962	!
1963	! Store level indices for interpolation
1964	!
1965	! for the pressure interpolation should be comparing
1966	! pint(column,lev) with M_hybi(lev)*M_ps_cam_col(month,column,chunk)
1967	!
1968	lev_interp_comp = .false.
1969	do kk=kkstart,paerlev-1
1970	if(.not.lev_interp_comp) then
1971	do i=1,ncol
1972	v_coord = pint(i,k)
1973	if (M_hybi(kk)Match_ps(i) .lt. v_coord .and. v_coord .le. M_hybi(kk+1)Match_ps(i)) then
1974	kupper(i) = kk
1975	kount = kount + 1
1976	end if
1977	end do
1978	!
1979	! If all indices for this level have been found, do the interpolation and
1980	! go to the next level
1981	!
1982	! Interpolate in pressure.
1983	!
1984	if (kount.eq.ncol) then
1985	do i=1,ncol
1986	do m=1,naer
1987	dpu = pint(i,k) - M_hybi(kupper(i))*Match_ps(i)
1988	dpl = M_hybi(kupper(i)+1)*Match_ps(i) - pint(i,k)
1989	AEROSOL(i,k,m) = &
1990	(AEROSOLc(i,kupper(i) ,m)*dpl + &
1991	AEROSOLc(i,kupper(i)+1,m)*dpu)/(dpl + dpu)
1992	enddo
1993	enddo !i
1994	lev_interp_comp = .true.
1995	end if
1996	end if
1997	end do
1998	!
1999	! If we've fallen through the kk=1,levsiz-1 loop, we cannot interpolate and
2000
2001	! must extrapolate from the bottom or top pressure level for at least some
2002	! of the longitude points.
2003	!
2004
2005	if(.not.lev_interp_comp) then
2006	do i=1,ncol
2007	do m=1,naer
2008	if (pint(i,k) .lt. M_hybi(1)*Match_ps(i)) then
2009	AEROSOL(i,k,m) = AEROSOLc(i,1,m)
2010	else if (pint(i,k) .gt. M_hybi(paerlev)*Match_ps(i)) then
2011	AEROSOL(i,k,m) = 0.0
2012	else
2013	dpu = pint(i,k) - M_hybi(kupper(i))*Match_ps(i)
2014	dpl = M_hybi(kupper(i)+1)*Match_ps(i) - pint(i,k)
2015	AEROSOL(i,k,m) = &
2016	(AEROSOLc(i,kupper(i) ,m)*dpl + &
2017	AEROSOLc(i,kupper(i)+1,m)*dpu)/(dpl + dpu)
2018	end if
2019	enddo
2020	end do
2021
2022	if (kount.gt.ncol) then
2023	call endrun ('VERT_INTERPOLATE: Bad data: non-monotonicity suspected in dependent variable')
2024	end if
2025	end if
2026	end do
2027
2028	! call t_startf ('vi_checks')
2029	!
2030	! aerosol mass beneath lowest interface (pverp) must be 0
2031	!
2032	AEROSOL(1:ncol,pverp,:) = 0.
2033	!
2034	! Set mass in layer to zero whenever it is less than
2035	! 1.e-40 kg/m^2 in the layer
2036	!
2037	do m = 1, naer
2038	do k = 1, pver
2039	do i = 1, ncol
2040	if (AEROSOL(i,k,m) < 1.e-40_r8) AEROSOL(i,k,m) = 0.
2041	end do
2042	end do
2043	end do
2044	!
2045	! Set mass in layer to zero whenever it is less than
2046	! 10^-15 relative to column total mass
2047	! convert back to mass mixing ratios.
2048	! exit if mmr is negative
2049	!
2050	do m = 1, naer
2051	do k = 1, pver
2052	do i = 1, ncol
2053	AER_diff = AEROSOL(i,k,m) - AEROSOL(i,k+1,m)
2054	if( abs(AER_diff) < 1e-15*AEROSOL(i,1,m)) then
2055	AER_diff = 0.
2056	end if
2057	m_to_mmr = gravmks / (pint(i,k+1)-pint(i,k))
2058	AEROSOL_mmr(i,k,m)= AER_diff * m_to_mmr
2059	if (AEROSOL_mmr(i,k,m) < 0) then
2060	write(6,*)'vert_interpolate: mmr < 0, m, col, lev, mmr',m, i, k, AEROSOL_mmr(i,k,m)
2061	write(6,*)'vert_interpolate: aerosol(k),(k+1)',AEROSOL(i,k,m),AEROSOL(i,k+1,m)
2062	write(6,*)'vert_interpolate: pint(k+1),(k)',pint(i,k+1),pint(i,k)
2063	write(6,*)'n,c',n,c
2064	call endrun()
2065	end if
2066	end do
2067	end do
2068	end do
2069
2070	! call t_stopf ('vi_checks')
2071	! call t_stopf ('vert_interpolate')
2072
2073	return
2074	end subroutine vert_interpolate
2075
2076
2077	!===============================================================================
2078	subroutine cldems(lchnk ,ncol ,pcols, pver, pverp, clwp ,fice ,rei ,emis )
2079	!-----------------------------------------------------------------------
2080	!
2081	! Purpose:
2082	! Compute cloud emissivity using cloud liquid water path (g/m**2)
2083	!
2084	! Method:
2085	! <Describe the algorithm(s) used in the routine.>
2086	! <Also include any applicable external references.>
2087	!
2088	! Author: J.T. Kiehl
2089	!
2090	!-----------------------------------------------------------------------
2091
2092	implicit none
2093	!------------------------------Parameters-------------------------------
2094	!
2095	real(r8) kabsl ! longwave liquid absorption coeff (m**2/g)
2096	parameter (kabsl = 0.090361)
2097	!
2098	!------------------------------Arguments--------------------------------
2099	!
2100	! Input arguments
2101	!
2102	integer, intent(in) :: lchnk ! chunk identifier
2103	integer, intent(in) :: ncol ! number of atmospheric columns
2104	integer, intent(in) :: pcols, pver, pverp
2105
2106	real(r8), intent(in) :: clwp(pcols,pver) ! cloud liquid water path (g/m**2)
2107	real(r8), intent(in) :: rei(pcols,pver) ! ice effective drop size (microns)
2108	real(r8), intent(in) :: fice(pcols,pver) ! fractional ice content within cloud
2109	!
2110	! Output arguments
2111	!
2112	real(r8), intent(out) :: emis(pcols,pver) ! cloud emissivity (fraction)
2113	!
2114	!---------------------------Local workspace-----------------------------
2115	!
2116	integer i,k ! longitude, level indices
2117	real(r8) kabs ! longwave absorption coeff (m**2/g)
2118	real(r8) kabsi ! ice absorption coefficient
2119	!
2120	!-----------------------------------------------------------------------
2121	!
2122	do k=1,pver
2123	do i=1,ncol
2124	kabsi = 0.005 + 1./rei(i,k)
2125	kabs = kabsl(1.-fice(i,k)) + kabsifice(i,k)
2126	emis(i,k) = 1. - exp(-1.66kabsclwp(i,k))
2127	end do
2128	end do
2129	!
2130	return
2131	end subroutine cldems
2132
2133	!===============================================================================
2134	subroutine cldovrlap(lchnk ,ncol ,pcols, pver, pverp, pint ,cld ,nmxrgn ,pmxrgn )
2135	!-----------------------------------------------------------------------
2136	!
2137	! Purpose:
2138	! Partitions each column into regions with clouds in neighboring layers.
2139	! This information is used to implement maximum overlap in these regions
2140	! with random overlap between them.
2141	! On output,
2142	! nmxrgn contains the number of regions in each column
2143	! pmxrgn contains the interface pressures for the lower boundaries of
2144	! each region!
2145	! Method:
2146
2147	!
2148	! Author: W. Collins
2149	!
2150	!-----------------------------------------------------------------------
2151
2152	implicit none
2153	!
2154	! Input arguments
2155	!
2156	integer, intent(in) :: lchnk ! chunk identifier
2157	integer, intent(in) :: ncol ! number of atmospheric columns
2158	integer, intent(in) :: pcols, pver, pverp
2159
2160	real(r8), intent(in) :: pint(pcols,pverp) ! Interface pressure
2161	real(r8), intent(in) :: cld(pcols,pver) ! Fractional cloud cover
2162	!
2163	! Output arguments
2164	!
2165	real(r8), intent(out) :: pmxrgn(pcols,pverp)! Maximum values of pressure for each
2166	! maximally overlapped region.
2167	! 0->pmxrgn(i,1) is range of pressure for
2168	! 1st region,pmxrgn(i,1)->pmxrgn(i,2) for
2169	! 2nd region, etc
2170	integer nmxrgn(pcols) ! Number of maximally overlapped regions
2171	!
2172	!---------------------------Local variables-----------------------------
2173	!
2174	integer i ! Longitude index
2175	integer k ! Level index
2176	integer n ! Max-overlap region counter
2177
2178	real(r8) pnm(pcols,pverp) ! Interface pressure
2179
2180	logical cld_found ! Flag for detection of cloud
2181	logical cld_layer(pver) ! Flag for cloud in layer
2182	!
2183	!------------------------------------------------------------------------
2184	!
2185
2186	do i = 1, ncol
2187	cld_found = .false.
2188	cld_layer(:) = cld(i,:) > 0.0_r8
2189	pmxrgn(i,:) = 0.0
2190	pnm(i,:)=pint(i,:)*10.
2191	n = 1
2192	do k = 1, pver
2193	if (cld_layer(k) .and. .not. cld_found) then
2194	cld_found = .true.
2195	else if ( .not. cld_layer(k) .and. cld_found) then
2196	cld_found = .false.
2197	if (count(cld_layer(k:pver)) == 0) then
2198	exit
2199	endif
2200	pmxrgn(i,n) = pnm(i,k)
2201	n = n + 1
2202	endif
2203	end do
2204	pmxrgn(i,n) = pnm(i,pverp)
2205	nmxrgn(i) = n
2206	end do
2207
2208	return
2209	end subroutine cldovrlap
2210
2211	!===============================================================================
2212	subroutine cldclw(lchnk ,ncol ,pcols, pver, pverp, zi ,clwp ,tpw ,hl )
2213	!-----------------------------------------------------------------------
2214	!
2215	! Purpose:
2216	! Evaluate cloud liquid water path clwp (g/m**2)
2217	!
2218	! Method:
2219	! <Describe the algorithm(s) used in the routine.>
2220	! <Also include any applicable external references.>
2221	!
2222	! Author: J.T. Kiehl
2223	!
2224	!-----------------------------------------------------------------------
2225
2226	implicit none
2227
2228	!
2229	! Input arguments
2230	!
2231	integer, intent(in) :: lchnk ! chunk identifier
2232	integer, intent(in) :: ncol ! number of atmospheric columns
2233	integer, intent(in) :: pcols, pver, pverp
2234
2235	real(r8), intent(in) :: zi(pcols,pverp) ! height at layer interfaces(m)
2236	real(r8), intent(in) :: tpw(pcols) ! total precipitable water (mm)
2237	!
2238	! Output arguments
2239	!
2240	real(r8) clwp(pcols,pver) ! cloud liquid water path (g/m**2)
2241	real(r8) hl(pcols) ! liquid water scale height
2242	real(r8) rhl(pcols) ! 1/hl
2243
2244	!
2245	!---------------------------Local workspace-----------------------------
2246	!
2247	integer i,k ! longitude, level indices
2248	real(r8) clwc0 ! reference liquid water concentration (g/m**3)
2249	real(r8) emziohl(pcols,pverp) ! exp(-zi/hl)
2250	!
2251	!-----------------------------------------------------------------------
2252	!
2253	! Set reference liquid water concentration
2254	!
2255	clwc0 = 0.21
2256	!
2257	! Diagnose liquid water scale height from precipitable water
2258	!
2259	do i=1,ncol
2260	hl(i) = 700.0*log(max(tpw(i)+1.0_r8,1.0_r8))
2261	rhl(i) = 1.0/hl(i)
2262	end do
2263	!
2264	! Evaluate cloud liquid water path (vertical integral of exponential fn)
2265	!
2266	do k=1,pverp
2267	do i=1,ncol
2268	emziohl(i,k) = exp(-zi(i,k)*rhl(i))
2269	end do
2270	end do
2271	do k=1,pver
2272	do i=1,ncol
2273	clwp(i,k) = clwc0hl(i)(emziohl(i,k+1) - emziohl(i,k))
2274	end do
2275	end do
2276	!
2277	return
2278	end subroutine cldclw
2279
2280
2281	!===============================================================================
2282	subroutine reltab(ncol, pcols, pver, t, landfrac, landm, icefrac, rel, snowh)
2283	!-----------------------------------------------------------------------
2284	!
2285	! Purpose:
2286	! Compute cloud water size
2287	!
2288	! Method:
2289	! analytic formula following the formulation originally developed by J. T. Kiehl
2290	!
2291	! Author: Phil Rasch
2292	!
2293	!-----------------------------------------------------------------------
2294	! use physconst, only: tmelt
2295	implicit none
2296	!------------------------------Arguments--------------------------------
2297	!
2298	! Input arguments
2299	!
2300	integer, intent(in) :: ncol
2301	integer, intent(in) :: pcols, pver
2302	real(r8), intent(in) :: landfrac(pcols) ! Land fraction
2303	real(r8), intent(in) :: icefrac(pcols) ! Ice fraction
2304	real(r8), intent(in) :: snowh(pcols) ! Snow depth over land, water equivalent (m)
2305	real(r8), intent(in) :: landm(pcols) ! Land fraction ramping to zero over ocean
2306	real(r8), intent(in) :: t(pcols,pver) ! Temperature
2307
2308	!
2309	! Output arguments
2310	!
2311	real(r8), intent(out) :: rel(pcols,pver) ! Liquid effective drop size (microns)
2312	!
2313	!---------------------------Local workspace-----------------------------
2314	!
2315	integer i,k ! Lon, lev indices
2316	real(r8) rliqland ! liquid drop size if over land
2317	real(r8) rliqocean ! liquid drop size if over ocean
2318	real(r8) rliqice ! liquid drop size if over sea ice
2319	!
2320	!-----------------------------------------------------------------------
2321	!
2322	rliqocean = 14.0_r8
2323	rliqice = 14.0_r8
2324	rliqland = 8.0_r8
2325	do k=1,pver
2326	do i=1,ncol
2327	! jrm Reworked effective radius algorithm
2328	! Start with temperature-dependent value appropriate for continental air
2329	! Note: findmcnew has a pressure dependence here
2330	rel(i,k) = rliqland + (rliqocean-rliqland) * min(1.0_r8,max(0.0_r8,(tmelt-t(i,k))*0.05))
2331	! Modify for snow depth over land
2332	rel(i,k) = rel(i,k) + (rliqocean-rel(i,k)) * min(1.0_r8,max(0.0_r8,snowh(i)*10.))
2333	! Ramp between polluted value over land to clean value over ocean.
2334	rel(i,k) = rel(i,k) + (rliqocean-rel(i,k)) * min(1.0_r8,max(0.0_r8,1.0-landm(i)))
2335	! Ramp between the resultant value and a sea ice value in the presence of ice.
2336	rel(i,k) = rel(i,k) + (rliqice-rel(i,k)) * min(1.0_r8,max(0.0_r8,icefrac(i)))
2337	! end jrm
2338	end do
2339	end do
2340	end subroutine reltab
2341	!===============================================================================
2342	subroutine reitab(ncol, pcols, pver, t, re)
2343	!
2344
2345	integer, intent(in) :: ncol, pcols, pver
2346	real(r8), intent(out) :: re(pcols,pver)
2347	real(r8), intent(in) :: t(pcols,pver)
2348	real(r8) corr
2349	integer i
2350	integer k
2351	integer index
2352	!
2353	do k=1,pver
2354	do i=1,ncol
2355	index = int(t(i,k)-179.)
2356	index = min(max(index,1),94)
2357	corr = t(i,k) - int(t(i,k))
2358	re(i,k) = retab(index)*(1.-corr) &
2359	+retab(index+1)*corr
2360	! re(i,k) = amax1(amin1(re(i,k),30.),10.)
2361	end do
2362	end do
2363	!
2364	return
2365	end subroutine reitab
2366
2367	function exp_interpol(x, f, y) result(g)
2368
2369	! Purpose:
2370	! interpolates f(x) to point y
2371	! assuming f(x) = f(x0) exp a(x - x0)
2372	! where a = ( ln f(x1) - ln f(x0) ) / (x1 - x0)
2373	! x0 <= x <= x1
2374	! assumes x is monotonically increasing
2375
2376	! Author: D. Fillmore
2377
2378	! use shr_kind_mod, only: r8 => shr_kind_r8
2379
2380	implicit none
2381
2382	real(r8), intent(in), dimension(:) :: x ! grid points
2383	real(r8), intent(in), dimension(:) :: f ! grid function values
2384	real(r8), intent(in) :: y ! interpolation point
2385	real(r8) :: g ! interpolated function value
2386
2387	integer :: k ! interpolation point index
2388	integer :: n ! length of x
2389	real(r8) :: a
2390
2391	n = size(x)
2392
2393	! find k such that x(k) < y =< x(k+1)
2394	! set k = 1 if y <= x(1) and k = n-1 if y > x(n)
2395
2396	if (y <= x(1)) then
2397	k = 1
2398	else if (y >= x(n)) then
2399	k = n - 1
2400	else
2401	k = 1
2402	do while (y > x(k+1) .and. k < n)
2403	k = k + 1
2404	end do
2405	end if
2406
2407	! interpolate
2408	a = ( log( f(k+1) / f(k) ) ) / ( x(k+1) - x(k) )
2409	g = f(k) * exp( a * (y - x(k)) )
2410
2411	end function exp_interpol
2412
2413	function lin_interpol(x, f, y) result(g)
2414
2415	! Purpose:
2416	! interpolates f(x) to point y
2417	! assuming f(x) = f(x0) + a * (x - x0)
2418	! where a = ( f(x1) - f(x0) ) / (x1 - x0)
2419	! x0 <= x <= x1
2420	! assumes x is monotonically increasing
2421
2422	! Author: D. Fillmore
2423
2424	! use shr_kind_mod, only: r8 => shr_kind_r8
2425
2426	implicit none
2427
2428	real(r8), intent(in), dimension(:) :: x ! grid points
2429	real(r8), intent(in), dimension(:) :: f ! grid function values
2430	real(r8), intent(in) :: y ! interpolation point
2431	real(r8) :: g ! interpolated function value
2432
2433	integer :: k ! interpolation point index
2434	integer :: n ! length of x
2435	real(r8) :: a
2436
2437	n = size(x)
2438
2439	! find k such that x(k) < y =< x(k+1)
2440	! set k = 1 if y <= x(1) and k = n-1 if y > x(n)
2441
2442	if (y <= x(1)) then
2443	k = 1
2444	else if (y >= x(n)) then
2445	k = n - 1
2446	else
2447	k = 1
2448	do while (y > x(k+1) .and. k < n)
2449	k = k + 1
2450	end do
2451	end if
2452
2453	! interpolate
2454	a = ( f(k+1) - f(k) ) / ( x(k+1) - x(k) )
2455	g = f(k) + a * (y - x(k))
2456
2457	end function lin_interpol
2458
2459	function lin_interpol2(x, f, y) result(g)
2460
2461	! Purpose:
2462	! interpolates f(x) to point y
2463	! assuming f(x) = f(x0) + a * (x - x0)
2464	! where a = ( f(x1) - f(x0) ) / (x1 - x0)
2465	! x0 <= x <= x1
2466	! assumes x is monotonically increasing
2467
2468	! Author: D. Fillmore :: J. Done changed from r8 to r4
2469
2470	implicit none
2471
2472	real, intent(in), dimension(:) :: x ! grid points
2473	real, intent(in), dimension(:) :: f ! grid function values
2474	real, intent(in) :: y ! interpolation point
2475	real :: g ! interpolated function value
2476
2477	integer :: k ! interpolation point index
2478	integer :: n ! length of x
2479	real :: a
2480
2481	n = size(x)
2482
2483	! find k such that x(k) < y =< x(k+1)
2484	! set k = 1 if y <= x(1) and k = n-1 if y > x(n)
2485
2486	if (y <= x(1)) then
2487	k = 1
2488	else if (y >= x(n)) then
2489	k = n - 1
2490	else
2491	k = 1
2492	do while (y > x(k+1) .and. k < n)
2493	k = k + 1
2494	end do
2495	end if
2496
2497	! interpolate
2498	a = ( f(k+1) - f(k) ) / ( x(k+1) - x(k) )
2499	g = f(k) + a * (y - x(k))
2500
2501	end function lin_interpol2
2502
2503
2504	subroutine getfactors (cycflag, np1, cdayminus, cdayplus, cday, &
2505	fact1, fact2)
2506	!---------------------------------------------------------------------------
2507	!
2508	! Purpose: Determine time interpolation factors (normally for a boundary dataset)
2509	! for linear interpolation.
2510	!
2511	! Method: Assume 365 days per year. Output variable fact1 will be the weight to
2512	! apply to data at calendar time "cdayminus", and fact2 the weight to apply
2513	! to data at time "cdayplus". Combining these values will produce a result
2514	! valid at time "cday". Output arguments fact1 and fact2 will be between
2515	! 0 and 1, and fact1 + fact2 = 1 to roundoff.
2516	!
2517	! Author: Jim Rosinski
2518	!
2519	!---------------------------------------------------------------------------
2520	implicit none
2521	!
2522	! Arguments
2523	!
2524	logical, intent(in) :: cycflag ! flag indicates whether dataset is being cycled yearly
2525
2526	integer, intent(in) :: np1 ! index points to forward time slice matching cdayplus
2527
2528	real(r8), intent(in) :: cdayminus ! calendar day of rearward time slice
2529	real(r8), intent(in) :: cdayplus ! calendar day of forward time slice
2530	real(r8), intent(in) :: cday ! calenar day to be interpolated to
2531	real(r8), intent(out) :: fact1 ! time interpolation factor to apply to rearward time slice
2532	real(r8), intent(out) :: fact2 ! time interpolation factor to apply to forward time slice
2533
2534	! character(len=*), intent(in) :: str ! string to be added to print in case of error (normally the callers name)
2535	!
2536	! Local workspace
2537	!
2538	real(r8) :: deltat ! time difference (days) between cdayminus and cdayplus
2539	real(r8), parameter :: daysperyear = 365. ! number of days in a year
2540	!
2541	! Initial sanity checks
2542	!
2543	! if (np1 == 1 .and. .not. cycflag) then
2544	! call endrun ('GETFACTORS:'//str//' cycflag false and forward month index = Jan. not allowed')
2545	! end if
2546
2547	! if (np1 < 1) then
2548	! call endrun ('GETFACTORS:'//str//' input arg np1 must be > 0')
2549	! end if
2550
2551	if (cycflag) then
2552	if ((cday < 1.) .or. (cday > (daysperyear+1.))) then
2553	write(6,*) 'GETFACTORS:', ' bad cday=',cday
2554	call endrun ()
2555	end if
2556	else
2557	if (cday < 1.) then
2558	write(6,*) 'GETFACTORS:', ' bad cday=',cday
2559	call endrun ()
2560	end if
2561	end if
2562	!
2563	! Determine time interpolation factors. Account for December-January
2564	! interpolation if dataset is being cycled yearly.
2565	!
2566	if (cycflag .and. np1 == 1) then ! Dec-Jan interpolation
2567	deltat = cdayplus + daysperyear - cdayminus
2568	if (cday > cdayplus) then ! We are in December
2569	fact1 = (cdayplus + daysperyear - cday)/deltat
2570	fact2 = (cday - cdayminus)/deltat
2571	else ! We are in January
2572	fact1 = (cdayplus - cday)/deltat
2573	fact2 = (cday + daysperyear - cdayminus)/deltat
2574	end if
2575	else
2576	deltat = cdayplus - cdayminus
2577	fact1 = (cdayplus - cday)/deltat
2578	fact2 = (cday - cdayminus)/deltat
2579	end if
2580
2581	if (.not. validfactors (fact1, fact2)) then
2582	write(6,*) 'GETFACTORS: ', ' bad fact1 and/or fact2=', fact1, fact2
2583	call endrun ()
2584	end if
2585
2586	return
2587	end subroutine getfactors
2588
2589	logical function validfactors (fact1, fact2)
2590	!---------------------------------------------------------------------------
2591	!
2592	! Purpose: check sanity of time interpolation factors to within 32-bit roundoff
2593	!
2594	!---------------------------------------------------------------------------
2595	implicit none
2596
2597	real(r8), intent(in) :: fact1, fact2 ! time interpolation factors
2598
2599	validfactors = .true.
2600	if (abs(fact1+fact2-1.) > 1.e-6 .or. &
2601	fact1 > 1.000001 .or. fact1 < -1.e-6 .or. &
2602	fact2 > 1.000001 .or. fact2 < -1.e-6) then
2603
2604	validfactors = .false.
2605	end if
2606
2607	return
2608	end function validfactors
2609
2610	subroutine get_rf_scales(scales)
2611
2612	real(r8), intent(out)::scales(naer_all) ! scale aerosols by this amount
2613
2614	integer i ! loop index
2615
2616	scales(idxBG) = bgscl_rf
2617	scales(idxSUL) = sulscl_rf
2618	scales(idxSSLT) = ssltscl_rf
2619
2620	do i = idxCARBONfirst, idxCARBONfirst+numCARBON-1
2621	scales(i) = carscl_rf
2622	enddo
2623
2624	do i = idxDUSTfirst, idxDUSTfirst+numDUST-1
2625	scales(i) = dustscl_rf
2626	enddo
2627
2628	scales(idxVOLC) = volcscl_rf
2629
2630	end subroutine get_rf_scales
2631
2632	function psi(tpx,iband)
2633	!
2634	! History: First version for Hitran 1996 (C/H/E)
2635	! Current version for Hitran 2000 (C/LT/E)
2636	! Short function for Hulst-Curtis-Godson temperature factors for
2637	! computing effective H2O path
2638	! Line data for H2O: Hitran 2000, plus H2O patches v11.0 for 1341 missing
2639	! lines between 500 and 2820 cm^-1.
2640	! See cfa-www.harvard.edu/HITRAN
2641	! Isotopes of H2O: all
2642	! Line widths: air-broadened only (self set to 0)
2643	! Code for line strengths and widths: GENLN3
2644	! Reference: Edwards, D.P., 1992: GENLN2, A General Line-by-Line Atmospheric
2645	! Transmittance and Radiance Model, Version 3.0 Description
2646	! and Users Guide, NCAR/TN-367+STR, 147 pp.
2647	!
2648	! Note: functions have been normalized by dividing by their values at
2649	! a path temperature of 160K
2650	!
2651	! spectral intervals:
2652	! 1 = 0-800 cm^-1 and 1200-2200 cm^-1
2653	! 2 = 800-1200 cm^-1
2654	!
2655	! Formulae: Goody and Yung, Atmospheric Radiation: Theoretical Basis,
2656	! 2nd edition, Oxford University Press, 1989.
2657	! Psi: function for pressure along path
2658	! eq. 6.30, p. 228
2659	!
2660	real(r8),intent(in):: tpx ! path temperature
2661	integer, intent(in):: iband ! band to process
2662	real(r8) psi ! psi for given band
2663	real(r8),parameter :: psi_r0(nbands) = (/ 5.65308452E-01, -7.30087891E+01/)
2664	real(r8),parameter :: psi_r1(nbands) = (/ 4.07519005E-03, 1.22199547E+00/)
2665	real(r8),parameter :: psi_r2(nbands) = (/-1.04347237E-05, -7.12256227E-03/)
2666	real(r8),parameter :: psi_r3(nbands) = (/ 1.23765354E-08, 1.47852825E-05/)
2667
2668	psi = (((psi_r3(iband) * tpx) + psi_r2(iband)) * tpx + psi_r1(iband)) * tpx + psi_r0(iband)
2669	end function psi
2670
2671	function phi(tpx,iband)
2672	!
2673	! History: First version for Hitran 1996 (C/H/E)
2674	! Current version for Hitran 2000 (C/LT/E)
2675	! Short function for Hulst-Curtis-Godson temperature factors for
2676	! computing effective H2O path
2677	! Line data for H2O: Hitran 2000, plus H2O patches v11.0 for 1341 missing
2678	! lines between 500 and 2820 cm^-1.
2679	! See cfa-www.harvard.edu/HITRAN
2680	! Isotopes of H2O: all
2681	! Line widths: air-broadened only (self set to 0)
2682	! Code for line strengths and widths: GENLN3
2683	! Reference: Edwards, D.P., 1992: GENLN2, A General Line-by-Line Atmospheric
2684	! Transmittance and Radiance Model, Version 3.0 Description
2685	! and Users Guide, NCAR/TN-367+STR, 147 pp.
2686	!
2687	! Note: functions have been normalized by dividing by their values at
2688	! a path temperature of 160K
2689	!
2690	! spectral intervals:
2691	! 1 = 0-800 cm^-1 and 1200-2200 cm^-1
2692	! 2 = 800-1200 cm^-1
2693	!
2694	! Formulae: Goody and Yung, Atmospheric Radiation: Theoretical Basis,
2695	! 2nd edition, Oxford University Press, 1989.
2696	! Phi: function for H2O path
2697	! eq. 6.25, p. 228
2698	!
2699	real(r8),intent(in):: tpx ! path temperature
2700	integer, intent(in):: iband ! band to process
2701	real(r8) phi ! phi for given band
2702	real(r8),parameter :: phi_r0(nbands) = (/ 9.60917711E-01, -2.21031342E+01/)
2703	real(r8),parameter :: phi_r1(nbands) = (/ 4.86076751E-04, 4.24062610E-01/)
2704	real(r8),parameter :: phi_r2(nbands) = (/-1.84806265E-06, -2.95543415E-03/)
2705	real(r8),parameter :: phi_r3(nbands) = (/ 2.11239959E-09, 7.52470896E-06/)
2706
2707	phi = (((phi_r3(iband) * tpx) + phi_r2(iband)) * tpx + phi_r1(iband)) &
2708	* tpx + phi_r0(iband)
2709	end function phi
2710
2711	function fh2oself( temp )
2712	!
2713	! Short function for H2O self-continuum temperature factor in
2714	! calculation of effective H2O self-continuum path length
2715	!
2716	! H2O Continuum: CKD 2.4
2717	! Code for continuum: GENLN3
2718	! Reference: Edwards, D.P., 1992: GENLN2, A General Line-by-Line Atmospheric
2719	! Transmittance and Radiance Model, Version 3.0 Description
2720	! and Users Guide, NCAR/TN-367+STR, 147 pp.
2721	!
2722	! In GENLN, the temperature scaling of the self-continuum is handled
2723	! by exponential interpolation/extrapolation from observations at
2724	! 260K and 296K by:
2725	!
2726	! TFAC = (T(IPATH) - 296.0)/(260.0 - 296.0)
2727	! CSFFT = CSFF296(CSFF260/CSFF296)*TFAC
2728	!
2729	! For 800-1200 cm^-1, (CSFF260/CSFF296) ranges from ~2.1 to ~1.9
2730	! with increasing wavenumber. The ratio <CSFF260>/<CSFF296>,
2731	! where <> indicates average over wavenumber, is ~2.07
2732	!
2733	! fh2oself is (<CSFF260>/<CSFF296>)**TFAC
2734	!
2735	real(r8),intent(in) :: temp ! path temperature
2736	real(r8) fh2oself ! mean ratio of self-continuum at temp and 296K
2737
2738	fh2oself = 2.0727484**((296.0 - temp) / 36.0)
2739	end function fh2oself
2740
2741	! from wv_saturation.F90
2742
2743	subroutine esinti(epslon ,latvap ,latice ,rh2o ,cpair ,tmelt )
2744	!-----------------------------------------------------------------------
2745	!
2746	! Purpose:
2747	! Initialize es lookup tables
2748	!
2749	! Method:
2750	! <Describe the algorithm(s) used in the routine.>
2751	! <Also include any applicable external references.>
2752	!
2753	! Author: J. Hack
2754	!
2755	!-----------------------------------------------------------------------
2756	! use shr_kind_mod, only: r8 => shr_kind_r8
2757	! use wv_saturation, only: gestbl
2758	implicit none
2759	!------------------------------Arguments--------------------------------
2760	!
2761	! Input arguments
2762	!
2763	real(r8), intent(in) :: epslon ! Ratio of h2o to dry air molecular weights
2764	real(r8), intent(in) :: latvap ! Latent heat of vaporization
2765	real(r8), intent(in) :: latice ! Latent heat of fusion
2766	real(r8), intent(in) :: rh2o ! Gas constant for water vapor
2767	real(r8), intent(in) :: cpair ! Specific heat of dry air
2768	real(r8), intent(in) :: tmelt ! Melting point of water (K)
2769	!
2770	!---------------------------Local workspace-----------------------------
2771	!
2772	real(r8) tmn ! Minimum temperature entry in table
2773	real(r8) tmx ! Maximum temperature entry in table
2774	real(r8) trice ! Trans range from es over h2o to es over ice
2775	logical ip ! Ice phase (true or false)
2776	!
2777	!-----------------------------------------------------------------------
2778	!
2779	! Specify control parameters first
2780	!
2781	tmn = 173.16
2782	tmx = 375.16
2783	trice = 20.00
2784	ip = .true.
2785	!
2786	! Call gestbl to build saturation vapor pressure table.
2787	!
2788	call gestbl(tmn ,tmx ,trice ,ip ,epslon , &
2789	latvap ,latice ,rh2o ,cpair ,tmelt )
2790	!
2791	return
2792	end subroutine esinti
2793
2794	subroutine gestbl(tmn ,tmx ,trice ,ip ,epsil , &
2795	latvap ,latice ,rh2o ,cpair ,tmeltx )
2796	!-----------------------------------------------------------------------
2797	!
2798	! Purpose:
2799	! Builds saturation vapor pressure table for later lookup procedure.
2800	!
2801	! Method:
2802	! Uses Goff & Gratch (1946) relationships to generate the table
2803	! according to a set of free parameters defined below. Auxiliary
2804	! routines are also included for making rapid estimates (well with 1%)
2805	! of both es and d(es)/dt for the particular table configuration.
2806	!
2807	! Author: J. Hack
2808	!
2809	!-----------------------------------------------------------------------
2810	! use pmgrid, only: masterproc
2811	implicit none
2812	!------------------------------Arguments--------------------------------
2813	!
2814	! Input arguments
2815	!
2816	real(r8), intent(in) :: tmn ! Minimum temperature entry in es lookup table
2817	real(r8), intent(in) :: tmx ! Maximum temperature entry in es lookup table
2818	real(r8), intent(in) :: epsil ! Ratio of h2o to dry air molecular weights
2819	real(r8), intent(in) :: trice ! Transition range from es over range to es over ice
2820	real(r8), intent(in) :: latvap ! Latent heat of vaporization
2821	real(r8), intent(in) :: latice ! Latent heat of fusion
2822	real(r8), intent(in) :: rh2o ! Gas constant for water vapor
2823	real(r8), intent(in) :: cpair ! Specific heat of dry air
2824	real(r8), intent(in) :: tmeltx ! Melting point of water (K)
2825	!
2826	!---------------------------Local variables-----------------------------
2827	!
2828	real(r8) t ! Temperature
2829	real(r8) rgasv
2830	real(r8) cp
2831	real(r8) hlatf
2832	real(r8) ttrice
2833	real(r8) hlatv
2834	integer n ! Increment counter
2835	integer lentbl ! Calculated length of lookup table
2836	integer itype ! Ice phase: 0 -> no ice phase
2837	! 1 -> ice phase, no transition
2838	! -x -> ice phase, x degree transition
2839	logical ip ! Ice phase logical flag
2840	logical icephs
2841	!
2842	!-----------------------------------------------------------------------
2843	!
2844	! Set es table parameters
2845	!
2846	tmin = tmn ! Minimum temperature entry in table
2847	tmax = tmx ! Maximum temperature entry in table
2848	ttrice = trice ! Trans. range from es over h2o to es over ice
2849	icephs = ip ! Ice phase (true or false)
2850	!
2851	! Set physical constants required for es calculation
2852	!
2853	epsqs = epsil
2854	hlatv = latvap
2855	hlatf = latice
2856	rgasv = rh2o
2857	cp = cpair
2858	tmelt = tmeltx
2859	!
2860	lentbl = INT(tmax-tmin+2.000001)
2861	if (lentbl .gt. plenest) then
2862	write(6,9000) tmax, tmin, plenest
2863	call endrun ('GESTBL') ! Abnormal termination
2864	end if
2865	!
2866	! Begin building es table.
2867	! Check whether ice phase requested.
2868	! If so, set appropriate transition range for temperature
2869	!
2870	if (icephs) then
2871	if (ttrice /= 0.0) then
2872	itype = -ttrice
2873	else
2874	itype = 1
2875	end if
2876	else
2877	itype = 0
2878	end if
2879	!
2880	t = tmin - 1.0
2881	do n=1,lentbl
2882	t = t + 1.0
2883	call gffgch(t,estbl(n),itype)
2884	end do
2885	!
2886	do n=lentbl+1,plenest
2887	estbl(n) = -99999.0
2888	end do
2889	!
2890	! Table complete -- Set coefficients for polynomial approximation of
2891	! difference between saturation vapor press over water and saturation
2892	! pressure over ice for -ttrice < t < 0 (degrees C). NOTE: polynomial
2893	! is valid in the range -40 < t < 0 (degrees C).
2894	!
2895	! --- Degree 5 approximation ---
2896	!
2897	pcf(1) = 5.04469588506e-01
2898	pcf(2) = -5.47288442819e+00
2899	pcf(3) = -3.67471858735e-01
2900	pcf(4) = -8.95963532403e-03
2901	pcf(5) = -7.78053686625e-05
2902	!
2903	! --- Degree 6 approximation ---
2904	!
2905	!-----pcf(1) = 7.63285250063e-02
2906	!-----pcf(2) = -5.86048427932e+00
2907	!-----pcf(3) = -4.38660831780e-01
2908	!-----pcf(4) = -1.37898276415e-02
2909	!-----pcf(5) = -2.14444472424e-04
2910	!-----pcf(6) = -1.36639103771e-06
2911	!
2912	if (masterproc) then
2913	write(6,)' SATURATION VAPOR PRESSURE TABLE COMPLETED *'
2914	end if
2915
2916	return
2917	!
2918	9000 format('GESTBL: FATAL ERROR *********************************',/, &
2919	' TMAX AND TMIN REQUIRE A LARGER DIMENSION ON THE LENGTH', &
2920	' OF THE SATURATION VAPOR PRESSURE TABLE ESTBL(PLENEST)',/, &
2921	' TMAX, TMIN, AND PLENEST => ', 2f7.2, i3)
2922	!
2923	end subroutine gestbl
2924
2925	subroutine gffgch(t ,es ,itype )
2926	!-----------------------------------------------------------------------
2927	!
2928	! Purpose:
2929	! Computes saturation vapor pressure over water and/or over ice using
2930	! Goff & Gratch (1946) relationships.
2931	! <Say what the routine does>
2932	!
2933	! Method:
2934	! T (temperature), and itype are input parameters, while es (saturation
2935	! vapor pressure) is an output parameter. The input parameter itype
2936	! serves two purposes: a value of zero indicates that saturation vapor
2937	! pressures over water are to be returned (regardless of temperature),
2938	! while a value of one indicates that saturation vapor pressures over
2939	! ice should be returned when t is less than freezing degrees. If itype
2940	! is negative, its absolute value is interpreted to define a temperature
2941	! transition region below freezing in which the returned
2942	! saturation vapor pressure is a weighted average of the respective ice
2943	! and water value. That is, in the temperature range 0 => -itype
2944	! degrees c, the saturation vapor pressures are assumed to be a weighted
2945	! average of the vapor pressure over supercooled water and ice (all
2946	! water at 0 c; all ice at -itype c). Maximum transition range => 40 c
2947	!
2948	! Author: J. Hack
2949	!
2950	!-----------------------------------------------------------------------
2951	! use shr_kind_mod, only: r8 => shr_kind_r8
2952	! use physconst, only: tmelt
2953	! use abortutils, only: endrun
2954
2955	implicit none
2956	!------------------------------Arguments--------------------------------
2957	!
2958	! Input arguments
2959	!
2960	real(r8), intent(in) :: t ! Temperature
2961	!
2962	! Output arguments
2963	!
2964	integer, intent(inout) :: itype ! Flag for ice phase and associated transition
2965
2966	real(r8), intent(out) :: es ! Saturation vapor pressure
2967	!
2968	!---------------------------Local variables-----------------------------
2969	!
2970	real(r8) e1 ! Intermediate scratch variable for es over water
2971	real(r8) e2 ! Intermediate scratch variable for es over water
2972	real(r8) eswtr ! Saturation vapor pressure over water
2973	real(r8) f ! Intermediate scratch variable for es over water
2974	real(r8) f1 ! Intermediate scratch variable for es over water
2975	real(r8) f2 ! Intermediate scratch variable for es over water
2976	real(r8) f3 ! Intermediate scratch variable for es over water
2977	real(r8) f4 ! Intermediate scratch variable for es over water
2978	real(r8) f5 ! Intermediate scratch variable for es over water
2979	real(r8) ps ! Reference pressure (mb)
2980	real(r8) t0 ! Reference temperature (freezing point of water)
2981	real(r8) term1 ! Intermediate scratch variable for es over ice
2982	real(r8) term2 ! Intermediate scratch variable for es over ice
2983	real(r8) term3 ! Intermediate scratch variable for es over ice
2984	real(r8) tr ! Transition range for es over water to es over ice
2985	real(r8) ts ! Reference temperature (boiling point of water)
2986	real(r8) weight ! Intermediate scratch variable for es transition
2987	integer itypo ! Intermediate scratch variable for holding itype
2988	!
2989	!-----------------------------------------------------------------------
2990	!
2991	! Check on whether there is to be a transition region for es
2992	!
2993	if (itype < 0) then
2994	tr = abs(float(itype))
2995	itypo = itype
2996	itype = 1
2997	else
2998	tr = 0.0
2999	itypo = itype
3000	end if
3001	if (tr > 40.0) then
3002	write(6,900) tr
3003	call endrun ('GFFGCH') ! Abnormal termination
3004	end if
3005	!
3006	if(t < (tmelt - tr) .and. itype == 1) go to 10
3007	!
3008	! Water
3009	!
3010	ps = 1013.246
3011	ts = 373.16
3012	e1 = 11.344*(1.0 - t/ts)
3013	e2 = -3.49149*(ts/t - 1.0)
3014	f1 = -7.90298*(ts/t - 1.0)
3015	f2 = 5.02808*log10(ts/t)
3016	f3 = -1.3816(10.0*e1 - 1.0)/10000000.0
3017	f4 = 8.1328(10.0*e2 - 1.0)/1000.0
3018	f5 = log10(ps)
3019	f = f1 + f2 + f3 + f4 + f5
3020	es = (10.0*f)100.0
3021	eswtr = es
3022	!
3023	if(t >= tmelt .or. itype == 0) go to 20
3024	!
3025	! Ice
3026	!
3027	10 continue
3028	t0 = tmelt
3029	term1 = 2.01889049/(t0/t)
3030	term2 = 3.56654*log(t0/t)
3031	term3 = 20.947031*(t0/t)
3032	es = 575.185606e10*exp(-(term1 + term2 + term3))
3033	!
3034	if (t < (tmelt - tr)) go to 20
3035	!
3036	! Weighted transition between water and ice
3037	!
3038	weight = min((tmelt - t)/tr,1.0_r8)
3039	es = weightes + (1.0 - weight)eswtr
3040	!
3041	20 continue
3042	itype = itypo
3043	return
3044	!
3045	900 format('GFFGCH: FATAL ERROR ******************************',/, &
3046	'TRANSITION RANGE FOR WATER TO ICE SATURATION VAPOR', &
3047	' PRESSURE, TR, EXCEEDS MAXIMUM ALLOWABLE VALUE OF', &
3048	' 40.0 DEGREES C',/, ' TR = ',f7.2)
3049	!
3050	end subroutine gffgch
3051
3052	real(r8) function estblf( td )
3053	!
3054	! Saturation vapor pressure table lookup
3055	!
3056	real(r8), intent(in) :: td ! Temperature for saturation lookup
3057	!
3058	real(r8) :: e ! intermediate variable for es look-up
3059	real(r8) :: ai
3060	integer :: i
3061	!
3062	e = max(min(td,tmax),tmin) ! partial pressure
3063	i = int(e-tmin)+1
3064	ai = aint(e-tmin)
3065	estblf = (tmin+ai-e+1.)* &
3066	estbl(i)-(tmin+ai-e)* &
3067	estbl(i+1)
3068	end function estblf
3069
3070
3071	function findvalue(ix,n,ain,indxa)
3072	!-----------------------------------------------------------------------
3073	!
3074	! Purpose:
3075	! Subroutine for finding ix-th smallest value in the array
3076	! The elements are rearranged so that the ix-th smallest
3077	! element is in the ix place and all smaller elements are
3078	! moved to the elements up to ix (with random order).
3079	!
3080	! Algorithm: Based on the quicksort algorithm.
3081	!
3082	! Author: T. Craig
3083	!
3084	!-----------------------------------------------------------------------
3085	! use shr_kind_mod, only: r8 => shr_kind_r8
3086	implicit none
3087	!
3088	! arguments
3089	!
3090	integer, intent(in) :: ix ! element to search for
3091	integer, intent(in) :: n ! total number of elements
3092	integer, intent(inout):: indxa(n) ! array of integers
3093	real(r8), intent(in) :: ain(n) ! array to search
3094	!
3095	integer findvalue ! return value
3096	!
3097	! local variables
3098	!
3099	integer i,j
3100	integer il,im,ir
3101
3102	integer ia
3103	integer itmp
3104	!
3105	!---------------------------Routine-----------------------------
3106	!
3107	il=1
3108	ir=n
3109	do
3110	if (ir-il <= 1) then
3111	if (ir-il == 1) then
3112	if (ain(indxa(ir)) < ain(indxa(il))) then
3113	itmp=indxa(il)
3114	indxa(il)=indxa(ir)
3115	indxa(ir)=itmp
3116	endif
3117	endif
3118	findvalue=indxa(ix)
3119	return
3120	else
3121	im=(il+ir)/2
3122	itmp=indxa(im)
3123	indxa(im)=indxa(il+1)
3124	indxa(il+1)=itmp
3125	if (ain(indxa(il+1)) > ain(indxa(ir))) then
3126	itmp=indxa(il+1)
3127	indxa(il+1)=indxa(ir)
3128	indxa(ir)=itmp
3129	endif
3130	if (ain(indxa(il)) > ain(indxa(ir))) then
3131	itmp=indxa(il)
3132	indxa(il)=indxa(ir)
3133	indxa(ir)=itmp
3134	endif
3135	if (ain(indxa(il+1)) > ain(indxa(il))) then
3136	itmp=indxa(il+1)
3137	indxa(il+1)=indxa(il)
3138	indxa(il)=itmp
3139	endif
3140	i=il+1
3141	j=ir
3142	ia=indxa(il)
3143	do
3144	do
3145	i=i+1
3146	if (ain(indxa(i)) >= ain(ia)) exit
3147	end do
3148	do
3149	j=j-1
3150	if (ain(indxa(j)) <= ain(ia)) exit
3151	end do
3152	if (j < i) exit
3153	itmp=indxa(i)
3154	indxa(i)=indxa(j)
3155	indxa(j)=itmp
3156	end do
3157	indxa(il)=indxa(j)
3158	indxa(j)=ia
3159	if (j >= ix)ir=j-1
3160	if (j <= ix)il=i
3161	endif
3162	end do
3163	end function findvalue
3164
3165
3166	subroutine radini(gravx ,cpairx ,epsilox ,stebolx, pstdx )
3167	!-----------------------------------------------------------------------
3168	!
3169	! Purpose:
3170	! Initialize various constants for radiation scheme; note that
3171	! the radiation scheme uses cgs units.
3172	!
3173	! Method:
3174	! <Describe the algorithm(s) used in the routine.>
3175	! <Also include any applicable external references.>
3176	!
3177	! Author: W. Collins (H2O parameterization) and J. Kiehl
3178	!
3179	!-----------------------------------------------------------------------
3180	! use shr_kind_mod, only: r8 => shr_kind_r8
3181	! use ppgrid, only: pver, pverp
3182	! use comozp, only: cplos, cplol
3183	! use pmgrid, only: masterproc, plev, plevp
3184	! use radae, only: radaeini
3185	! use physconst, only: mwdry, mwco2
3186	#if ( defined SPMD )
3187	! use mpishorthand
3188	#endif
3189	implicit none
3190
3191	!------------------------------Arguments--------------------------------
3192	!
3193	! Input arguments
3194	!
3195	real, intent(in) :: gravx ! Acceleration of gravity (MKS)
3196	real, intent(in) :: cpairx ! Specific heat of dry air (MKS)
3197	real, intent(in) :: epsilox ! Ratio of mol. wght of H2O to dry air
3198	real, intent(in) :: stebolx ! Stefan-Boltzmann's constant (MKS)
3199	real(r8), intent(in) :: pstdx ! Standard pressure (Pascals)
3200	!
3201	!---------------------------Local variables-----------------------------
3202	!
3203	integer k ! Loop variable
3204
3205	real(r8) v0 ! Volume of a gas at stp (m**3/kmol)
3206	real(r8) p0 ! Standard pressure (pascals)
3207	real(r8) amd ! Effective molecular weight of dry air (kg/kmol)
3208	real(r8) goz ! Acceleration of gravity (m/s**2)
3209	!
3210	!-----------------------------------------------------------------------
3211	!
3212	! Set general radiation consts; convert to cgs units where appropriate:
3213	!
3214	gravit = 100.*gravx
3215	rga = 1./gravit
3216	gravmks = gravx
3217	cpair = 1.e4*cpairx
3218	epsilo = epsilox
3219	sslp = 1.013250e6
3220	stebol = 1.e3*stebolx
3221	rgsslp = 0.5/(gravit*sslp)
3222	dpfo3 = 2.5e-3
3223	dpfco2 = 5.0e-3
3224	dayspy = 365.
3225	pie = 4.*atan(1.)
3226	!
3227	! Initialize ozone data.
3228	!
3229	v0 = 22.4136 ! Volume of a gas at stp (m**3/kmol)
3230	p0 = 0.1*sslp ! Standard pressure (pascals)
3231	amd = 28.9644 ! Molecular weight of dry air (kg/kmol)
3232	goz = gravx ! Acceleration of gravity (m/s**2)
3233	!
3234	! Constants for ozone path integrals (multiplication by 100 for unit
3235	! conversion to cgs from mks):
3236	!
3237	cplos = v0/(amdgoz) 100.0
3238	cplol = v0/(amdgozp0)0.5100.0
3239	!
3240	! Derived constants
3241	! If the top model level is above ~90 km (0.1 Pa), set the top level to compute
3242	! longwave cooling to about 80 km (1 Pa)
3243	! WRF: assume top level > 0.1 mb
3244	! if (hypm(1) .lt. 0.1) then
3245	! do k = 1, pver
3246	! if (hypm(k) .lt. 1.) ntoplw = k
3247	! end do
3248	! else
3249	ntoplw = 1
3250	! end if
3251	! if (masterproc) then
3252	! write (6,*) 'RADINI: ntoplw =',ntoplw, ' pressure:',hypm(ntoplw)
3253	! endif
3254
3255	call radaeini( pstdx, mwdry, mwco2 )
3256	return
3257	end subroutine radini
3258
3259	subroutine oznini(ozmixm,pin,levsiz,num_months,XLAT, &
3260	ids, ide, jds, jde, kds, kde, &
3261	ims, ime, jms, jme, kms, kme, &
3262	its, ite, jts, jte, kts, kte)
3263	!
3264	! This subroutine assumes uniform distribution of ozone concentration.
3265	! It should be replaced by monthly climatology that varies latitudinally and vertically
3266	!
3267
3268	IMPLICIT NONE
3269
3270	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
3271	ims,ime, jms,jme, kms,kme, &
3272	its,ite, jts,jte, kts,kte
3273
3274	INTEGER, INTENT(IN ) :: levsiz, num_months
3275
3276	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: XLAT
3277
3278	REAL, DIMENSION( ims:ime, levsiz, jms:jme, num_months ), &
3279	INTENT(OUT ) :: OZMIXM
3280
3281	REAL, DIMENSION(levsiz), INTENT(OUT ) :: PIN
3282
3283	! Local
3284	INTEGER, PARAMETER :: latsiz = 64
3285	INTEGER, PARAMETER :: lonsiz = 1
3286	INTEGER :: i, j, k, itf, jtf, ktf, m, pin_unit, lat_unit, oz_unit
3287	REAL :: interp_pt
3288	CHARACTER*256 :: message
3289
3290	REAL, DIMENSION( lonsiz, levsiz, latsiz, num_months ) :: &
3291	OZMIXIN
3292
3293	REAL, DIMENSION(latsiz) :: lat_ozone
3294
3295	jtf=min0(jte,jde-1)
3296	ktf=min0(kte,kde-1)
3297	itf=min0(ite,ide-1)
3298
3299
3300	!-- read in ozone pressure data
3301
3302	WRITE(message,*)'num_months = ',num_months
3303	CALL wrf_debug(50,message)
3304
3305	pin_unit = 27
3306	OPEN(pin_unit, FILE='ozone_plev.formatted',FORM='FORMATTED',STATUS='OLD')
3307	do k = 1,levsiz
3308	READ (pin_unit,*)pin(k)
3309	end do
3310	close(27)
3311
3312	do k=1,levsiz
3313	pin(k) = pin(k)*100.
3314	end do
3315
3316	!-- read in ozone lat data
3317
3318	lat_unit = 28
3319	OPEN(lat_unit, FILE='ozone_lat.formatted',FORM='FORMATTED',STATUS='OLD')
3320	do j = 1,latsiz
3321	READ (lat_unit,*)lat_ozone(j)
3322	end do
3323	close(28)
3324
3325
3326	!-- read in ozone data
3327
3328	oz_unit = 29
3329	OPEN(oz_unit, FILE='ozone.formatted',FORM='FORMATTED',STATUS='OLD')
3330
3331	do m=2,num_months
3332	do j=1,latsiz ! latsiz=64
3333	do k=1,levsiz ! levsiz=59
3334	do i=1,lonsiz ! lonsiz=1
3335	READ (oz_unit,*)ozmixin(i,k,j,m)
3336	enddo
3337	enddo
3338	enddo
3339	enddo
3340	close(29)
3341
3342
3343	!-- latitudinally interpolate ozone data (and extend longitudinally)
3344	!-- using function lin_interpol2(x, f, y) result(g)
3345	! Purpose:
3346	! interpolates f(x) to point y
3347	! assuming f(x) = f(x0) + a * (x - x0)
3348	! where a = ( f(x1) - f(x0) ) / (x1 - x0)
3349	! x0 <= x <= x1
3350	! assumes x is monotonically increasing
3351	! real, intent(in), dimension(:) :: x ! grid points
3352	! real, intent(in), dimension(:) :: f ! grid function values
3353	! real, intent(in) :: y ! interpolation point
3354	! real :: g ! interpolated function value
3355	!---------------------------------------------------------------------------
3356
3357	do m=2,num_months
3358	do j=jts,jtf
3359	do k=1,levsiz
3360	do i=its,itf
3361	interp_pt=XLAT(i,j)
3362	ozmixm(i,k,j,m)=lin_interpol2(lat_ozone(:),ozmixin(1,k,:,m),interp_pt)
3363	enddo
3364	enddo
3365	enddo
3366	enddo
3367
3368	! Old code for fixed ozone
3369
3370	! pin(1)=70.
3371	! DO k=2,levsiz
3372	! pin(k)=pin(k-1)+16.
3373	! ENDDO
3374
3375	! DO k=1,levsiz
3376	! pin(k) = pin(k)*100.
3377	! end do
3378
3379	! DO m=1,num_months
3380	! DO j=jts,jtf
3381	! DO i=its,itf
3382	! DO k=1,2
3383	! ozmixm(i,k,j,m)=1.e-6
3384	! ENDDO
3385	! DO k=3,levsiz
3386	! ozmixm(i,k,j,m)=1.e-7
3387	! ENDDO
3388	! ENDDO
3389	! ENDDO
3390	! ENDDO
3391
3392	END SUBROUTINE oznini
3393
3394
3395	subroutine aerosol_init(m_psp,m_psn,m_hybi,aerosolcp,aerosolcn,paerlev,naer_c,shalf,pptop, &
3396	ids, ide, jds, jde, kds, kde, &
3397	ims, ime, jms, jme, kms, kme, &
3398	its, ite, jts, jte, kts, kte)
3399	!
3400	! This subroutine assumes a uniform aerosol distribution in both time and space.
3401	! It should be modified if aerosol data are available from WRF-CHEM or other sources
3402	!
3403	IMPLICIT NONE
3404
3405	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
3406	ims,ime, jms,jme, kms,kme, &
3407	its,ite, jts,jte, kts,kte
3408
3409	INTEGER, INTENT(IN ) :: paerlev,naer_c
3410
3411	REAL, intent(in) :: pptop
3412	REAL, DIMENSION( kms:kme ), intent(in) :: shalf
3413
3414	REAL, DIMENSION( ims:ime, paerlev, jms:jme, naer_c ), &
3415	INTENT(INOUT ) :: aerosolcn , aerosolcp
3416
3417	REAL, DIMENSION(paerlev), INTENT(OUT ) :: m_hybi
3418	REAL, DIMENSION( ims:ime, jms:jme), INTENT(OUT ) :: m_psp,m_psn
3419
3420	REAL :: psurf
3421	real, dimension(29) :: hybi
3422	integer k ! index through vertical levels
3423
3424	INTEGER :: i, j, itf, jtf, ktf,m
3425
3426	data hybi/0, 0.0065700002014637, 0.0138600002974272, 0.023089999333024, &
3427	0.0346900001168251, 0.0491999983787537, 0.0672300010919571, &
3428	0.0894500017166138, 0.116539999842644, 0.149159997701645, &
3429	0.187830001115799, 0.232859998941422, 0.284209996461868, &
3430	0.341369986534119, 0.403340011835098, 0.468600004911423, &
3431	0.535290002822876, 0.601350009441376, 0.66482001543045, &
3432	0.724009990692139, 0.777729988098145, 0.825269997119904, &
3433	0.866419970989227, 0.901350021362305, 0.930540025234222, &
3434	0.954590022563934, 0.974179983139038, 0.990000009536743, 1/
3435
3436	jtf=min0(jte,jde-1)
3437	ktf=min0(kte,kde-1)
3438	itf=min0(ite,ide-1)
3439
3440	do k=1,paerlev
3441	m_hybi(k)=hybi(k)
3442	enddo
3443
3444	!
3445	! mxaerl = max number of levels (from bottom) for background aerosol
3446	! Limit background aerosol height to regions below 900 mb
3447	!
3448
3449	psurf = 1.e05
3450	mxaerl = 0
3451	! do k=pver,1,-1
3452	do k=kms,kme-1
3453	! if (hypm(k) >= 9.e4) mxaerl = mxaerl + 1
3454	if (shalf(k)*psurf+pptop >= 9.e4) mxaerl = mxaerl + 1
3455	end do
3456	mxaerl = max(mxaerl,1)
3457	! if (masterproc) then
3458	write(6,*)'AEROSOLS: Background aerosol will be limited to ', &
3459	'bottom ',mxaerl,' model interfaces.'
3460	! 'bottom ',mxaerl,' model interfaces. Top interface is ', &
3461	! hypi(pverp-mxaerl),' pascals'
3462	! end if
3463
3464	DO j=jts,jtf
3465	DO i=its,itf
3466	m_psp(i,j)=psurf
3467	m_psn(i,j)=psurf
3468	ENDDO
3469	ENDDO
3470
3471	DO j=jts,jtf
3472	DO i=its,itf
3473	DO k=1,paerlev
3474	! aerosolc arrays are upward cumulative (kg/m2) at each level
3475	! Here we assume uniform vertical distribution (aerosolc linear with hybi)
3476	aerosolcp(i,k,j,idxSUL)=1.e-7*(1.-hybi(k))
3477	aerosolcn(i,k,j,idxSUL)=1.e-7*(1.-hybi(k))
3478	aerosolcp(i,k,j,idxSSLT)=1.e-22*(1.-hybi(k))
3479	aerosolcn(i,k,j,idxSSLT)=1.e-22*(1.-hybi(k))
3480	aerosolcp(i,k,j,idxDUSTfirst)=1.e-7*(1.-hybi(k))
3481	aerosolcn(i,k,j,idxDUSTfirst)=1.e-7*(1.-hybi(k))
3482	aerosolcp(i,k,j,idxDUSTfirst+1)=1.e-7*(1.-hybi(k))
3483	aerosolcn(i,k,j,idxDUSTfirst+1)=1.e-7*(1.-hybi(k))
3484	aerosolcp(i,k,j,idxDUSTfirst+2)=1.e-7*(1.-hybi(k))
3485	aerosolcn(i,k,j,idxDUSTfirst+2)=1.e-7*(1.-hybi(k))
3486	aerosolcp(i,k,j,idxDUSTfirst+3)=1.e-7*(1.-hybi(k))
3487	aerosolcn(i,k,j,idxDUSTfirst+3)=1.e-7*(1.-hybi(k))
3488	aerosolcp(i,k,j,idxOCPHO)=1.e-7*(1.-hybi(k))
3489	aerosolcn(i,k,j,idxOCPHO)=1.e-7*(1.-hybi(k))
3490	aerosolcp(i,k,j,idxBCPHO)=1.e-9*(1.-hybi(k))
3491	aerosolcn(i,k,j,idxBCPHO)=1.e-9*(1.-hybi(k))
3492	aerosolcp(i,k,j,idxOCPHI)=1.e-7*(1.-hybi(k))
3493	aerosolcn(i,k,j,idxOCPHI)=1.e-7*(1.-hybi(k))
3494	aerosolcp(i,k,j,idxBCPHI)=1.e-8*(1.-hybi(k))
3495	aerosolcn(i,k,j,idxBCPHI)=1.e-8*(1.-hybi(k))
3496	ENDDO
3497	ENDDO
3498	ENDDO
3499
3500	call aer_optics_initialize
3501
3502
3503	END subroutine aerosol_init
3504
3505	subroutine aer_optics_initialize
3506
3507	USE module_wrf_error
3508
3509	! use shr_kind_mod, only: r8 => shr_kind_r8
3510	! use pmgrid ! masterproc is here
3511	! use ioFileMod, only: getfil
3512
3513	!#if ( defined SPMD )
3514	! use mpishorthand
3515	!#endif
3516	implicit none
3517
3518	! include 'netcdf.inc'
3519
3520
3521	integer :: nrh_opac ! number of relative humidity values for OPAC data
3522	integer :: nbnd ! number of spectral bands, should be identical to nspint
3523	real(r8), parameter :: wgt_sscm = 6.0 / 7.0
3524	integer :: krh_opac ! rh index for OPAC rh grid
3525	integer :: krh ! another rh index
3526	integer :: ksz ! dust size bin index
3527	integer :: kbnd ! band index
3528
3529	real(r8) :: rh ! local relative humidity variable
3530
3531	integer, parameter :: irh=8
3532	real(r8) :: rh_opac(irh) ! OPAC relative humidity grid
3533	real(r8) :: ksul_opac(irh,nspint) ! sulfate extinction
3534	real(r8) :: wsul_opac(irh,nspint) ! single scattering albedo
3535	real(r8) :: gsul_opac(irh,nspint) ! asymmetry parameter
3536	real(r8) :: ksslt_opac(irh,nspint) ! sea-salt
3537	real(r8) :: wsslt_opac(irh,nspint)
3538	real(r8) :: gsslt_opac(irh,nspint)
3539	real(r8) :: kssam_opac(irh,nspint) ! sea-salt accumulation mode
3540	real(r8) :: wssam_opac(irh,nspint)
3541	real(r8) :: gssam_opac(irh,nspint)
3542	real(r8) :: ksscm_opac(irh,nspint) ! sea-salt coarse mode
3543	real(r8) :: wsscm_opac(irh,nspint)
3544	real(r8) :: gsscm_opac(irh,nspint)
3545	real(r8) :: kcphil_opac(irh,nspint) ! hydrophilic organic carbon
3546	real(r8) :: wcphil_opac(irh,nspint)
3547	real(r8) :: gcphil_opac(irh,nspint)
3548	real(r8) :: dummy(nspint)
3549
3550	LOGICAL :: opened
3551	LOGICAL , EXTERNAL :: wrf_dm_on_monitor
3552
3553	CHARACTER*80 errmess
3554	INTEGER cam_aer_unit
3555	integer :: i
3556
3557	! read aerosol optics data
3558
3559	IF ( wrf_dm_on_monitor() ) THEN
3560	DO i = 10,99
3561	INQUIRE ( i , OPENED = opened )
3562	IF ( .NOT. opened ) THEN
3563	cam_aer_unit = i
3564	GOTO 2010
3565	ENDIF
3566	ENDDO
3567	cam_aer_unit = -1
3568	2010 CONTINUE
3569	ENDIF
3570	CALL wrf_dm_bcast_bytes ( cam_aer_unit , IWORDSIZE )
3571	IF ( cam_aer_unit < 0 ) THEN
3572	CALL wrf_error_fatal ( 'module_ra_cam: aer_optics_initialize: Can not find unused fortran unit to read in lookup table.' )
3573	ENDIF
3574
3575	IF ( wrf_dm_on_monitor() ) THEN
3576	OPEN(cam_aer_unit,FILE='CAM_AEROPT_DATA', &
3577	FORM='UNFORMATTED',STATUS='OLD',ERR=9010)
3578	call wrf_debug(50,'reading CAM_AEROPT_DATA')
3579	ENDIF
3580
3581	#define DM_BCAST_MACRO(A) CALL wrf_dm_bcast_bytes ( A , size ( A ) * r8 )
3582
3583	IF ( wrf_dm_on_monitor() ) then
3584	READ (cam_aer_unit,ERR=9010) dummy
3585	READ (cam_aer_unit,ERR=9010) rh_opac
3586	READ (cam_aer_unit,ERR=9010) ksul_opac
3587	READ (cam_aer_unit,ERR=9010) wsul_opac
3588	READ (cam_aer_unit,ERR=9010) gsul_opac
3589	READ (cam_aer_unit,ERR=9010) kssam_opac
3590	READ (cam_aer_unit,ERR=9010) wssam_opac
3591	READ (cam_aer_unit,ERR=9010) gssam_opac
3592	READ (cam_aer_unit,ERR=9010) ksscm_opac
3593	READ (cam_aer_unit,ERR=9010) wsscm_opac
3594	READ (cam_aer_unit,ERR=9010) gsscm_opac
3595	READ (cam_aer_unit,ERR=9010) kcphil_opac
3596	READ (cam_aer_unit,ERR=9010) wcphil_opac
3597	READ (cam_aer_unit,ERR=9010) gcphil_opac
3598	READ (cam_aer_unit,ERR=9010) kcb
3599	READ (cam_aer_unit,ERR=9010) wcb
3600	READ (cam_aer_unit,ERR=9010) gcb
3601	READ (cam_aer_unit,ERR=9010) kdst
3602	READ (cam_aer_unit,ERR=9010) wdst
3603	READ (cam_aer_unit,ERR=9010) gdst
3604	READ (cam_aer_unit,ERR=9010) kbg
3605	READ (cam_aer_unit,ERR=9010) wbg
3606	READ (cam_aer_unit,ERR=9010) gbg
3607	READ (cam_aer_unit,ERR=9010) kvolc
3608	READ (cam_aer_unit,ERR=9010) wvolc
3609	READ (cam_aer_unit,ERR=9010) gvolc
3610	endif
3611
3612	DM_BCAST_MACRO(rh_opac)
3613	DM_BCAST_MACRO(ksul_opac)
3614	DM_BCAST_MACRO(wsul_opac)
3615	DM_BCAST_MACRO(gsul_opac)
3616	DM_BCAST_MACRO(kssam_opac)
3617	DM_BCAST_MACRO(wssam_opac)
3618	DM_BCAST_MACRO(gssam_opac)
3619	DM_BCAST_MACRO(ksscm_opac)
3620	DM_BCAST_MACRO(wsscm_opac)
3621	DM_BCAST_MACRO(gsscm_opac)
3622	DM_BCAST_MACRO(kcphil_opac)
3623	DM_BCAST_MACRO(wcphil_opac)
3624	DM_BCAST_MACRO(gcphil_opac)
3625	DM_BCAST_MACRO(kcb)
3626	DM_BCAST_MACRO(wcb)
3627	DM_BCAST_MACRO(gcb)
3628	DM_BCAST_MACRO(kvolc)
3629	DM_BCAST_MACRO(wvolc)
3630	DM_BCAST_MACRO(kdst)
3631	DM_BCAST_MACRO(wdst)
3632	DM_BCAST_MACRO(gdst)
3633	DM_BCAST_MACRO(kbg)
3634	DM_BCAST_MACRO(wbg)
3635	DM_BCAST_MACRO(gbg)
3636
3637	IF ( wrf_dm_on_monitor() ) CLOSE (cam_aer_unit)
3638
3639	! map OPAC aerosol species onto CAM aerosol species
3640	! CAM name OPAC name
3641	! sul or SO4 = suso sulfate soluble
3642	! sslt or SSLT = 1/7 ssam + 6/7 sscm sea-salt accumulation/coagulation mode
3643	! cphil or CPHI = waso water soluble (carbon)
3644	! cphob or CPHO = waso @ rh = 0
3645	! cb or BCPHI/BCPHO = soot
3646
3647	ksslt_opac(:,:) = (1.0 - wgt_sscm) * kssam_opac(:,:) + wgt_sscm * ksscm_opac(:,:)
3648
3649	wsslt_opac(:,:) = ( (1.0 - wgt_sscm) * kssam_opac(:,:) * wssam_opac(:,:) &
3650	+ wgt_sscm * ksscm_opac(:,:) * wsscm_opac(:,:) ) &
3651	/ ksslt_opac(:,:)
3652
3653	gsslt_opac(:,:) = ( (1.0 - wgt_sscm) * kssam_opac(:,:) * wssam_opac(:,:) * gssam_opac(:,:) &
3654	+ wgt_sscm * ksscm_opac(:,:) * wsscm_opac(:,:) * gsscm_opac(:,:) ) &
3655	/ ( ksslt_opac(:,:) * wsslt_opac(:,:) )
3656
3657	do i=1,nspint
3658	kcphob(i) = kcphil_opac(1,i)
3659	wcphob(i) = wcphil_opac(1,i)
3660	gcphob(i) = gcphil_opac(1,i)
3661	end do
3662
3663	! interpolate optical properties of hygrospopic aerosol species
3664	! onto a uniform relative humidity grid
3665
3666	nbnd = nspint
3667
3668	do krh = 1, nrh
3669	rh = 1.0_r8 / nrh * (krh - 1)
3670	do kbnd = 1, nbnd
3671	ksul(krh, kbnd) = exp_interpol( rh_opac, &
3672	ksul_opac(:, kbnd) / ksul_opac(1, kbnd), rh ) * ksul_opac(1, kbnd)
3673	wsul(krh, kbnd) = lin_interpol( rh_opac, &
3674	wsul_opac(:, kbnd) / wsul_opac(1, kbnd), rh ) * wsul_opac(1, kbnd)
3675	gsul(krh, kbnd) = lin_interpol( rh_opac, &
3676	gsul_opac(:, kbnd) / gsul_opac(1, kbnd), rh ) * gsul_opac(1, kbnd)
3677	ksslt(krh, kbnd) = exp_interpol( rh_opac, &
3678	ksslt_opac(:, kbnd) / ksslt_opac(1, kbnd), rh ) * ksslt_opac(1, kbnd)
3679	wsslt(krh, kbnd) = lin_interpol( rh_opac, &
3680	wsslt_opac(:, kbnd) / wsslt_opac(1, kbnd), rh ) * wsslt_opac(1, kbnd)
3681	gsslt(krh, kbnd) = lin_interpol( rh_opac, &
3682	gsslt_opac(:, kbnd) / gsslt_opac(1, kbnd), rh ) * gsslt_opac(1, kbnd)
3683	kcphil(krh, kbnd) = exp_interpol( rh_opac, &
3684	kcphil_opac(:, kbnd) / kcphil_opac(1, kbnd), rh ) * kcphil_opac(1, kbnd)
3685	wcphil(krh, kbnd) = lin_interpol( rh_opac, &
3686	wcphil_opac(:, kbnd) / wcphil_opac(1, kbnd), rh ) * wcphil_opac(1, kbnd)
3687	gcphil(krh, kbnd) = lin_interpol( rh_opac, &
3688	gcphil_opac(:, kbnd) / gcphil_opac(1, kbnd), rh ) * gcphil_opac(1, kbnd)
3689	end do
3690	end do
3691
3692	RETURN
3693	9010 CONTINUE
3694	WRITE( errmess , '(A35,I4)' ) 'module_ra_cam: error reading unit ',cam_aer_unit
3695	CALL wrf_error_fatal(errmess)
3696
3697	END subroutine aer_optics_initialize
3698
3699
3700	subroutine radaeini( pstdx, mwdryx, mwco2x )
3701
3702	USE module_wrf_error
3703
3704	!
3705	! Initialize radae module data
3706	!
3707	!
3708	! Input variables
3709	!
3710	real(r8), intent(in) :: pstdx ! Standard pressure (dynes/cm^2)
3711	real(r8), intent(in) :: mwdryx ! Molecular weight of dry air
3712	real(r8), intent(in) :: mwco2x ! Molecular weight of carbon dioxide
3713	!
3714	! Variables for loading absorptivity/emissivity
3715	!
3716	integer ncid_ae ! NetCDF file id for abs/ems file
3717
3718	integer pdimid ! pressure dimension id
3719	integer psize ! pressure dimension size
3720
3721	integer tpdimid ! path temperature dimension id
3722	integer tpsize ! path temperature size
3723
3724	integer tedimid ! emission temperature dimension id
3725	integer tesize ! emission temperature size
3726
3727	integer udimid ! u (H2O path) dimension id
3728	integer usize ! u (H2O path) dimension size
3729
3730	integer rhdimid ! relative humidity dimension id
3731	integer rhsize ! relative humidity dimension size
3732
3733	integer ah2onwid ! var. id for non-wndw abs.
3734	integer eh2onwid ! var. id for non-wndw ems.
3735	integer ah2owid ! var. id for wndw abs. (adjacent layers)
3736	integer cn_ah2owid ! var. id for continuum trans. for wndw abs.
3737	integer cn_eh2owid ! var. id for continuum trans. for wndw ems.
3738	integer ln_ah2owid ! var. id for line trans. for wndw abs.
3739	integer ln_eh2owid ! var. id for line trans. for wndw ems.
3740
3741	! character*(NF_MAX_NAME) tmpname! dummy variable for var/dim names
3742	character(len=256) locfn ! local filename
3743	integer tmptype ! dummy variable for variable type
3744	integer ndims ! number of dimensions
3745	! integer dims(NF_MAX_VAR_DIMS) ! vector of dimension ids
3746	integer natt ! number of attributes
3747	!
3748	! Variables for setting up H2O table
3749	!
3750	integer t ! path temperature
3751	integer tmin ! mininum path temperature
3752	integer tmax ! maximum path temperature
3753	integer itype ! type of sat. pressure (=0 -> H2O only)
3754	integer i
3755	real(r8) tdbl
3756
3757	LOGICAL :: opened
3758	LOGICAL , EXTERNAL :: wrf_dm_on_monitor
3759
3760	CHARACTER*80 errmess
3761	INTEGER cam_abs_unit
3762
3763	!
3764	! Constants to set
3765	!
3766	p0 = pstdx
3767	amd = mwdryx
3768	amco2 = mwco2x
3769	!
3770	! Coefficients for h2o emissivity and absorptivity for overlap of H2O
3771	! and trace gases.
3772	!
3773	c16 = coefj(3,1)/coefj(2,1)
3774	c17 = coefk(3,1)/coefk(2,1)
3775	c26 = coefj(3,2)/coefj(2,2)
3776	c27 = coefk(3,2)/coefk(2,2)
3777	c28 = .5
3778	c29 = .002053
3779	c30 = .1
3780	c31 = 3.0e-5
3781	!
3782	! Initialize further longwave constants referring to far wing
3783	! correction for overlap of H2O and trace gases; R&D refers to:
3784	!
3785	! Ramanathan, V. and P.Downey, 1986: A Nonisothermal
3786	! Emissivity and Absorptivity Formulation for Water Vapor
3787	! Journal of Geophysical Research, vol. 91., D8, pp 8649-8666
3788	!
3789	fwcoef = .1 ! See eq(33) R&D
3790	fwc1 = .30 ! See eq(33) R&D
3791	fwc2 = 4.5 ! See eq(33) and eq(34) in R&D
3792	fc1 = 2.6 ! See eq(34) R&D
3793
3794	IF ( wrf_dm_on_monitor() ) THEN
3795	DO i = 10,99
3796	INQUIRE ( i , OPENED = opened )
3797	IF ( .NOT. opened ) THEN
3798	cam_abs_unit = i
3799	GOTO 2010
3800	ENDIF
3801	ENDDO
3802	cam_abs_unit = -1
3803	2010 CONTINUE
3804	ENDIF
3805	CALL wrf_dm_bcast_bytes ( cam_abs_unit , IWORDSIZE )
3806	IF ( cam_abs_unit < 0 ) THEN
3807	CALL wrf_error_fatal ( 'module_ra_cam: radaeinit: Can not find unused fortran unit to read in lookup table.' )
3808	ENDIF
3809
3810	IF ( wrf_dm_on_monitor() ) THEN
3811	OPEN(cam_abs_unit,FILE='CAM_ABS_DATA', &
3812	FORM='UNFORMATTED',STATUS='OLD',ERR=9010)
3813	call wrf_debug(50,'reading CAM_ABS_DATA')
3814	ENDIF
3815
3816	#define DM_BCAST_MACRO(A) CALL wrf_dm_bcast_bytes ( A , size ( A ) * r8 )
3817
3818	IF ( wrf_dm_on_monitor() ) then
3819	READ (cam_abs_unit,ERR=9010) ah2onw
3820	READ (cam_abs_unit,ERR=9010) eh2onw
3821	READ (cam_abs_unit,ERR=9010) ah2ow
3822	READ (cam_abs_unit,ERR=9010) cn_ah2ow
3823	READ (cam_abs_unit,ERR=9010) cn_eh2ow
3824	READ (cam_abs_unit,ERR=9010) ln_ah2ow
3825	READ (cam_abs_unit,ERR=9010) ln_eh2ow
3826
3827	endif
3828
3829	DM_BCAST_MACRO(ah2onw)
3830	DM_BCAST_MACRO(eh2onw)
3831	DM_BCAST_MACRO(ah2ow)
3832	DM_BCAST_MACRO(cn_ah2ow)
3833	DM_BCAST_MACRO(cn_eh2ow)
3834	DM_BCAST_MACRO(ln_ah2ow)
3835	DM_BCAST_MACRO(ln_eh2ow)
3836
3837	IF ( wrf_dm_on_monitor() ) CLOSE (cam_abs_unit)
3838
3839	! Set up table of H2O saturation vapor pressures for use in calculation
3840	! effective path RH. Need separate table from table in wv_saturation
3841	! because:
3842	! (1. Path temperatures can fall below minimum of that table; and
3843	! (2. Abs/Emissivity tables are derived with RH for water only.
3844	!
3845	tmin = nint(min_tp_h2o)
3846	tmax = nint(max_tp_h2o)+1
3847	itype = 0
3848	do t = tmin, tmax
3849	! call gffgch(dble(t),estblh2o(t-tmin),itype)
3850	tdbl = t
3851	call gffgch(tdbl,estblh2o(t-tmin),itype)
3852	end do
3853
3854	RETURN
3855	9010 CONTINUE
3856	WRITE( errmess , '(A35,I4)' ) 'module_ra_cam: error reading unit ',cam_abs_unit
3857	CALL wrf_error_fatal(errmess)
3858	end subroutine radaeini
3859
3860
3861	end MODULE module_ra_cam_support

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format