Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_ra_cam_support.F @ 3431

Last change on this file since 3431 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 138.2 KB

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

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

Download in other formats:

Original Format